Neoteny is a developmental strategy wherein an organism reaches sexual maturity without associated adult characteristics. In salamanders, neoteny takes the form of individuals retaining aquatic larval characteristics such as external gills upon maturation. Mole salamanders (Ambystoma) occupy a wide range of habitats and areas across the North American continent and display examples of non-neotenic, facultatively neotenic, and obligate neotenic species, providing high variation for investigating the factors influencing the evolution of neoteny. Here, we use phylogenetic comparative methods to test existing hypotheses that neoteny is associated with elevational and latitudinal distribution, cave-associated isolation, and hybridisation-related polyploidy. We also test if neoteny influences the diversity of habitats a species can occupy since the restriction to aquatic life should constrain the availability of different niches. We find that neoteny tends to occur in a narrow latitudinal band between 20-30° North, with particularly narrow latitudinal ranges for obligate compared to facultative neotenic species (16-52° North). We also find that facultatively neotenic species occur at elevations more than twice as high as other species on average and that species with a higher frequency of neoteny typically have lower habitat diversity. Our results suggest that evolutionary transitions between non-neotenic and facultative neoteny states occur at relatively high and approximately equal rates. Moreover, we estimate that obligate neoteny cannot evolve directly from non-neotenic species (and vice versa), such that facultative neoteny acts as an evolutionary ‘stepping stone’ to and from obligate neoteny. However, our transition rate estimates suggest that obligate neoteny is lost > 4 times faster than it evolves, partly explaining the rarity of obligate species. These results support the hypothesis that low latitudes favour the evolution of neoteny, presumably linked to more stable (aquatic) environments due to reduced seasonality, but once evolved it may constrain the diversity of habitats.

Maximum elevation data was taken from the IUCN Red List database (IUCN 2021) for all species except the following; A. barbouri (Micheletti and Storfer 2015), A.bishopi, A. mabeei, A. maculatum, A. talpoideum, A. texanum, A. talpoideum, A.texanum (observations from AmphibiaWeb 2021, checked for elevation on www.freemaptools.com/elevation-finder.htm), A. californiense (Lanoo 2005), A. cingulatum, (animaldiversity.org 2021), A. jeffersonianum (Thompson and Gates 1982), and A. opacum (Klemens 1993). Latitudinal data (northern and southern distributional limits) were collected using the species distributions from the IUCN Red List database (IUCN 2021). The number of habitats occupied by each species was estimated by counting the second-level habitat types (for example 1.4 Forest -Temperate) listed in the IUCN Red List (IUCN 2021) according to their habitat classification scheme. Inclusion in the Ambystoma jeffersonianum polyploid complex was based on information from Robertson et al. (2006), Lowcock et al. (1991), and Bogart and Lichts (1986). The presence of the species in a cave environment was based on information from Briggler (2007), Jones and Thompson (1983), Sunny et al. (2014), Gorički et al. (2019), and Ryk (2019). We did not include reports of the species from undercut riverbanks as cave-living, as although semi-enclosed, these represent a microhabitat within a river quite distinct from caves as relevant to our hypothesis.

Data on neoteny were collected from amphibiaweb.org and Everson et al. (2021) and coded as a categorical variable with three levels: no neoteny (metamorphosing into terrestrial adults), facultative neoteny, and obligate neoteny. This initial three-state variable is herein referred to as “frequency of neoteny”. For our analyses, we considered three further alternative variables by recoding our initial three-state variable. The recoded variables considered were ‘neoteny’ (pooling facultative or obligate neoteny), ‘facultative neoteny’ (binary variable distinguishing facultative neoteny from both no and obligate neoteny), and obligate neoteny (binary variable distinguishing obligate neoteny from both no and facultative neoteny). We chose to use these alternative variables because it is possible that some of the proposed hypotheses might apply to (e.g.) obligate neoteny only, facultative neoteny, or any degree of neoteny (either as if they were the same thing or in an ordered way with obligate neoteny more extreme than facultative). Hence, using our multiple coding scheme strategy enables a more in-depth evaluation of the hypotheses and under what conditions they might apply.