We sampled two populations of each species (Fig. 1B); the first population of A. altamirani was in Organillos (19°31'38.17"N, 99°28'39.92"W (the datum in all coordinates described are: WGS–84), with an altitude of 3,335 MASL), and the second population was in Sehuayán (19°31'31"N, 99°26'09.52"W with an altitude of 3,185 MASL), both of which are in Sierra de las Cruces. The first population of A. rivulare was in Corral de Piedra (19°13'6.60"N, 99°57'54.77"W, with an altitude of 2,836 MASL), and the second population was in Raíces (19°9'37.26"N, 99°49'32.11"W, with an altitude of 3,225 MASL), both of which are in the Nevado de Toluca Volcano (NTV). We sampled the individuals with a fishing net, and we sampled 2 mm² of tail clips of adult mole salamanders. This methodology is a low-impact method that does not affect the survival or growth of the mole salamanders (Arntzen et al. 1999; Polich et al. 2013). Tissue was preserved in 90% ethanol and then frozen at -20 °C until processed. Finally, we released the mole salamanders after tail clipping. Our study received the approval of the ethics committee from Universidad Autónoma del Estado de México (3047-2011E, 4732/2019CIB and 9855714) and the collection permits of SEMARNAT (SEMARNAT: SGPA/DGVS/001777/18).

We extracted DNA following the manufacturer’s instructions for the GF-1 nucleic acid extraction kit (Vivantis Technologies, Subang Jaya, Malaysia), and we used it as a template for amplification of nine microsatellite loci following published protocols (Parra-Olea et al. 2007). PCR microsatellite products were multiplexed and run on an ABI Prism3730xl (Applied Biosystems, Foster City, CA, USA) with Rox-500 as an internal size standard. We obtained allele sizes with PEAKSCANNER 1.0 software (Applied Biosystems), and the fragment lengths were obtained with TANDEM 1.08 software (Matschiner and Salzburger 2009). In all runs we included negative controls in at least two runs and we repeat 2 samples per plate to account for genotyping error and to guarantee reproducibility.

We tested the presence of null alleles and large allele dropout in the MICROCHECKER 2.2.3 software (Van Oosterhout et al. 2004). In addition, in POPPR 2.4.1 (Kamvar et al. 2014) for R software (version 3.4.0; R Development Core Team 2017), we made an analysis to create a genotype accumulation curve we used for determining the minimum number of loci necessary to discriminate between individuals in each population of the species studied (Kamvar et al. 2014). This function randomly samples loci without replacement and counts the number of observed multilocus genotypes (Kamvar et al. 2014).

All analyses were done for each study location (it was determined that the sampling sites were independent populations by the STRUCTURE results, see below) and species. We calculated the observed (H_o) and expected (H_e) heterozygosity, the number of alleles (N_a), effective number of alleles (N_e), number of genotypes and the number of heterozygotes and homozygote genotypes in STRATA G 2.0.2 (Archer et al. 2017) and GENALEX. We calculated departures from Hardy-Weinberg equilibrium (HWE) and Linkage Disequilibrium (LD) between pairs of microsatellite loci in PEGAS package (Paradis 2010) implemented in R. These calculations were evaluated for each sampled locality and locus with a Markov chain approximation considering 10,000 dememorizations, 1,000 batches and 10,000 iterations per batch. In order to correct the P values, we used a False Discovery Rate (FDR) approach according to Benjamini and Hochberg (1995) implemented in the package FDRTOOL 1.2.15 (Strimmer 2008; Klaus and Strimmer 2013) for R.

We searched for a genetic structure pattern using several algorithms for each of the species and sampling sites. First, we used a Bayesian algorithm implemented in the STRUCTURE 2.3.4 software (Pritchard et al. 2000; Falush et al. 2003; Hubisz et al. 2009). The chosen parameters were correlated allele frequencies with 1,000,000 burn-in periods and 1,000,000 MCMC iterations (Falush et al. 2003). For the degree of admixture, a Dirichlet parameter was applied with correlated allele frequencies. Therefore, the most credible number of populations were estimated using the maximum value of ΔK (Evanno et al. 2005), applied in the STRUCTURE HARVESTER 0.6.92 software (Earl and vonHoldt 2012). The second method was the analysis of molecular variance (AMOVA) based on F_ST and R_ST as implemented by GENALEX 6 (Peakall and Smouse 2006). A Wilcoxon test with 30,000 permutations was applied to search significance, using the degree of similarity of the populations based on the populations’ genotypes in GENALEX 6. In order to detect the similarity degree of the populations of each species, we applied a Discriminant Principal Components Analysis (DPCA) based on the genotypes; we used the populations as priors and we ran the analysis in ADEGENET 2.0.1 (Jombart et al. 2016) and ADE4 1.7-6 (Dray and Siberchicot 2017) packages for the R software, and we calculated F_ST based on Weir and Cockerham (1984) in GENALEX 6. Finally, we tested the existence of the population structure by computing Minimum Spanning Networks (MSN) with Bruvo’s (Bruvo et al. 2004) and Nei´s distance algorithm (Nei 1972) with 1,000 bootstraps in POPPR 2.4.1 and MAGRITTR 1.5 (Bache and Wickham 2016) for the R package. This analysis visualizes the relationships among individuals and it can be a more adequate visualization tool than trees (Bache and Wickham 2016).

The historical signal of demographic fluctuations was explored for each population by applying a Bayesian algorithm implemented in MSVAR 0.4.1 software (Beaumont 1999). We estimated the rate of change (r) of the effective population size, defined as N_crnt/N_stbl (where N_crnt was the current inbreeding effective population size and N_stbl was the ancestral stable inbreeding effective population size). The r ratio was expressed in log10. Therefore, the population declined if we had a negative r value, was stable if r is equal to zero, and the population was expanded if the r value was positive (Gasca-Pineda et al. 2013; Sunny et al. 2015). In order to test for a genetic signature of recent bottlenecks, we used the BOTTLENECK 5.1.26 software (Cournet and Luikart 1996; Piry et al. 1999). We estimated the observed and expected heterozygosity under the two-phase model (TPM) because the TPM is an intermediate model of evolution which is considered more appropriate for microsatellites. The settings applied were for a 90% step-wise mutation model and 10% variance and were run with 10,000 replicates. Excess heterozygosity was tested using a Wilcoxon test. In order to explore the actual effective population size (N_e), we used the LD method implemented in the NEESTIMATOR 2 software (Do et al. 2014). We calculated the F_IS inbreeding values in GENALEX 6. Also as an inbreeding measure, we used the relatedness estimator (rqg) of Queller and Goodnight (1989), which was calculated by the GENALEX software. To test for significant differences among mean population relatedness, we calculated the upper and lower 95% confidence intervals for the expected range of rqg using 9999 permutations. These intervals corresponded to the range of rqg that would be expected if reproduction was random across populations. Additionally, we calculated confidence intervals for estimates of mean relatedness within a population to 95% by bootstrap resampling (9999 permutations). Population rqg values that fall above the 95% expected values indicate that processes such as inbreeding or genetic drift are increasing relatedness. Finally, relatedness among individuals was evaluated using the ML-RELATE software (Kalinowski et al. 2006).

One hundred and fifty-two individuals were sampled from two Ambystoma species (A. altamirani and A. rivulare), two locations were sampled and 38 tissues were collected from each locality.

We did not find evidence of null alleles or large allele dropout in the populations of each species. The genotype accumulation curve found that the minimum number of loci necessary to discriminate between individuals was eight (Suppl. material 1: Fig. S1). Therefore, we concluded that our study had enough loci (N = 9).

Across the nine loci in the A. altamirani populations in Organillos we found 3–7 alleles per locus and a total of 30 alleles; Sehuayán had 2–7 alleles per locus with a total of 26 alleles (Table 1, Suppl. material 1: Fig. S2). In Corral de Piedra the first A. rivulare population we found 2–5 alleles per locus with a total of 20 alleles, and in Raíces the second A. rivulare population we found 2–6 alleles per locus with a total of 33 alleles (Table 1, Suppl. material 1: Fig. S3). In the A. altamirani Organillos population we found 29 homozygote genotypes and 35 heterozygote genotypes, In the second A. altamirani population Sehuayán we found 16 homozygote genotypes and 22 heterozygote genotypes. In the A. rivulare Corral de Piedra population we found 26 homozygote genotypes and 19 heterozygote genotypes, and Raíces had 21 Homozygote genotypes and 20 heterozygote genotypes (Suppl. material 1: Table S1). In relation to the observed and expected heterozygosity values, the A. altamirani Organillos population showed lower observed heterozygosity values (H_o = 0.719±0.033) compared to Sehuayán (H_o = 0.857±0.029); the A. rivulare population Corral de Piedra showed lower observed heterozygosity values (H_o = 0.576±0.034) as compared to Raíces (H_o = 0.754±0.059; Table 1). False discovery rate correction tests found departures from HWE due to heterozygote deficiency in one locus in the populations of Organillos (A. altamirani), Corral de Piedra (A. rivulare) and Raíces (A. rivulare) (Suppl. material 1: Table S2). We did not find LD between any loci of either population for each species.

Bayesian assignment analyses corroborated high population divergences among populations (Fig. 2). The highest log likelihood given by STRUCTURE and ∆K method was K = 2 (LnPr = -1793.2) for A. altamirani and K = 2 (LnPr = -1457.6) for A. rivulare. The populations of A. altamirani present admixia with medium genetic differentiation among them (F_ST = 0.053; Table 2), but the populations of A. rivulare do not present admixia with high genetic differentiation among them (F_ST = 0.211; Table 2). In relation to the AMOVA results, for A. altamirani the results revealed the majority of genetic variation resided within populations (92%) followed by among populations (8%), with an F_ST fixation index of 0.084 and a p-value = 0.001, for A. rivulare the AMOVA results showed significant levels of genetic variation within populations (67%) and among populations (33%), with an F_ST fixation index of 0.332 and a p-value = 0.001(Suppl. material 1: Tables S3, S4). The DPCA and the MSN found the same patterns of population structure where the populations of A. altamirani are more similar among them, whereas the populations of A. rivulare are very different among them (Figs 3, 4).

Figure 2. 

A) Population genetic structure of Ambystoma altamirani analyzed with STRUCTURE. B) Evanno et al. (2005) plots for detecting the number of K groups that best fit the data. C) Population genetic structure of Ambystoma rivulare analyzed with STRUCTURE. D) Evanno et al. (2005) plots for detecting the number of K groups that best fit the data.

Figure 3. 

Scatterplot showing the relationships among the populations of each species based on a Discriminant Principal Components Analysis of 9 microsatellite genotypes.

Figure 4. 

A) and B) Minimum Spanning Networks with the Bruvo’s distance algorithm, representing the relationships among individuals and populations of each species. C) and D) Tree constructed by the NJ method using the estimated standardized genetic distances using the Nei´s distance algorithm (Nei 1972) with 1,000 bootstraps.

MSVAR results suggested that there has been a significant population size reduction in all the studied Ambystoma populations: Organillos, r = -0.972; Sehuayán, r = -1.403; Corral de Piedra, r = -1.031 and Raíces, r = -1.399. The bottleneck analysis detected genetic signs of recent demographic changes typical of bottleneck events, associated with a heterozygote excess in all populations: Organillos, P = 0.008; Sehuayán, P = 0.002; Corral de Piedra, P = 0.002; Raíces, P = 0.002 and Corral de Piedra and Sehuayán had a shifted distribution. The effective population size (N_e) estimated from LD was N_e = 34.7 (20.9–21.7, 95% CI) for Organillos, N_e = 44.1 (21.3–36.0, 95% CI) for Sehuayán, N_e = 57.6 (21.3–37.6, 95% CI) for Corral de Piedra and N_e = 41.5 (20.1–23.6, 95% CI) for Raíces. The F_IS statistic as an indicator of inbreeding for the A. altamirani populations showed negative and low inbreeding values (F_IS = -0.128; Table 2) and showed positive, low inbreeding values for the A. rivulare populations (F_IS = -0.095; Table 2). We found that mean pairwise relatedness (r) within populations (Fig. 5) was generally in accordance with that observed in other Ambystoma populations (Parra-Olea et al. 2012; Sunny et al. 2014a; Percino-Daniel et al. 2016). Organillos had low values of inbreeding [mean rqg = 0.247, confidence interval (CI) = 0.024 -(-0.025)], Sehuayán had medium to high levels of inbreeding [mean rqg = 0.458, confidence interval (CI) = 0.029 -(-0.025)], Corral de Piedra had the highest inbreeding values [mean rqg = 0.505, confidence interval (CI) = 0.028 -(-0.026)] and Raíces had medium to high inbreeding values [mean rqg=0.471, confidence interval (CI) = 0.032 -(-0.026). In all the populations, the rqg values fell above the 95% expected values from permutations which indicates that inbreeding or drift are increasing the relatedness; they also fell outside the range expected under panmixia. Furthermore, there were low values of inbreeding in the relatedness analysis, the proportion of relatedness of individuals within each population was similar and most of the individuals were unrelated, followed by full siblings, half-siblings and parents/offspring (Suppl. material 1: Table S5).

Figure 5. 

Mean within-lake pairwise relatedness coefficient rqg across the four Ambystoma populations studied. The green bars are 95% upper and lower expected values for a null distribution generated from 9999 permutations of data from all populations, and enclose the values expected if breeding were panmictic across all populations; relatedness in all sampled populations fell outside the range expected under panmixia. Black bar represents the observed mean relatedness in each population, the black bars are the upper and lower bootstrap value for each population.

The observed heterozygosity values were medium to high, and most of the genotypes were heterozygous with the exception of the two populations of A. rivulare (Corral de Piedra and Raíces) (Table 1, Suppl. material 1: Table S1). Ambystoma species have high levels of genetic diversity (Goprenko et al. 2007; Greenwald et al. 2009; Sunny et al. 2014a; Percino-Daniel et al. 2016) despite having a fragmented and limited distribution. However, the human population growth in the towns around the two populations of A. rivulare (Corral de Piedra and Raíces) has been increasing in recent decades, causing an influx of sewage, waste from local agriculture and pollutants from trout farms into the rivers. Also, in Corral de Piedra we did not find larvae, and each mole salamander had tail bites. This could be due to trout predation (Zamora et al. 2018) because of trout escape from farms, which favors native species reduction (Gamradt et al. 1997; Kiesecker and Blaustein 1997; Kiesecker et al. 2001). Introduced fish species have been linked to reductions in amphibians’ population sizes (Pearson and Goater 2009; Zambrano et al. 2010; Alcaraz et al. 2015; Zamora et al. 2018), sometimes to the point of extinction from direct consumption (Watson et al. 1991; Tyler et al. 1998) increasing the loss of genetic diversity. We found a significant deviation from the HWE proportions due to a heterozygote deficiency in Organillos, Corral de Piedra and Raíces. This is a common result when microsatellites are not specific for the species and for threatened species with fragmented populations (Degne et al. 2007; Spear and Storfer 2010; Vázquez-Domínguez et al. 2012; Sunny et al. 2014a).

Structure analysis suggests two populations for each species. The populations of A. rivulare showed no signs of admixia, although the populations of A. altamirani were more admixed (Figs 2–4). The populations of A. altamirani and A. rivulare we studied are ~67 km apart in linear distance, and the two populations of A. rivulare are ~16 km apart in linear distance, so the genetic exchange is extremely reduced. The two populations of A. altamirani are closer, ~2.5 km apart in linear distance so gene flow could still take place. In other studies, mole salamander migrations occur between temporal ponds and lakes as individuals look for food-rich habitats (Percino-Daniel et al. 2016). The known maximum dispersal distance in mole salamanders is less than 2 km (Smith and Green 2005; Percino-Daniel et al. 2016). Migration among other populations could be limited as a result of physical barriers like discontinuity of rivers, forests, roads, and towns, as well as the strong philopatric tendencies of the mole salamanders for breeding sites favoring high genetic structuring, even at small scales (Spear et al. 2005; Zamudio and Wieczorek 2007; Richardson 2012; Richardson and Urban 2013). These populations in the TMVB occupy some of the most disturbed areas of the country, mainly due to habitat fragmentation stream pollution and fish introduction (Recuero et al. 2010; Zambrano et al. 2010; Rodríguez-Amador et al. 2013; Lemos-Espinal et al. 2016). Since there are different factors affecting Ambystoma populations, we recommend some measures to try to minimize the impact of the anthropogenic activities on these amphibian populations.

The studied populations are isolated from other populations of mole salamanders. This phenomenon could explain the low N_e values (N_e = 34.7–57.6) found in all populations and the asymmetry in the proportions of males and females and differences in the reproductive success between individuals favour low N_e values (Wang 2009; Savage et al. 2010; Parra-Olea et al. 2012; Percino-Daniel et al. 2016) due to high asymmetry in reproductive success among members of a population (Savage et al. 2010). Some years, only a few individuals successfully breed, and the variance in mating success may contribute strongly to lower overall effective population sizes (Savage et al. 2010). Another explanation of the low N_e values found could be the bottleneck effect caused by the introduction of trout; however, we do not know when this introduction happened, as trout predate the early stages of development of Ambystomas and also eat the tail of adult Ambystomas, a phenomenon which has already been reported by other studies (Pilliod and Peterson 2001; Welsh et al. 2006; Zambrano et al. 2010; Martín-Torrijos et al. 2016). It has also been reported that trout can transmit pathogens and certain emerging infectious diseases (Johnson and Speare 2005; Fernández-Benéitez et al. 2008; Van den Berg et al. 2013; Sandoval-Sierra et al. 2014). Likewise, trout can compete for food with Ambystomas (Werner and Anholt 1996; Tyler et al. 1998). All these features can lead the population to a process of genetic bottleneck or genetic drift, which in turn reduces genetic diversity and makes the population lose fitness and the ability to adapt to changes in the environment (Frankham et al. 2005). In all populations, the rqg values were above the 95% expected values from permutations (Fig. 5), indicating that inbreeding or genetic drift are increasing the relatedness, and they fell outside the expected range under panmixia. Furthermore, there were low values of inbreeding in the relatedness analysis; the proportion of relatedness of individuals within each population was similar. Despite the lack of strong signs of inbreeding and relatedness, inbreeding and genetic drift are acting in these populations.

In order to conserve this species and all the species that live in the coniferous forests of TMVB, it is necessary to avoid excessive legal and illegal logging and give support to the local communities with incentives such as payments for ecosystem services. Also, we consider the implementation of an environmental education program to be fundamental to avoid excess logging; maintaining and increasing forest core areas in order to minimize the forest edges, also, preventing the loss of the largest forest patches in order to avoid deviations from circularity in patch shapes to increase the area of core habitat (Ewers and Didham 2006; González-Fernández et al. 2019). Therefore, we must conserve the endangered Abies-Pinus forest. In the TMVB, there are only 1346.9 km² of Abies forest (1.1%) and 6507.7 km² of Pinus forest (5.4%). It is therefore extremely important to conserve these forests, considering that Mexico had the seventh-highest net annual forest loss of any country in the world, and the clearing of primary forests averaged more than 1 percent per year (FAO 2011; Blackman et al. 2015), placing Mexico in the first levels of global deforestation (FAO 2006; Ellis and Porter-Bolland 2008). Likewise, the amphibians and reptiles have a very limited distribution, sometimes in very small areas with special characteristics that generate certain microhabitats suitable for their survival and reproduction (Sunny et al. 2014b; Lemos-Espinal et al. 2016, 2017; Zamora et al. 2018). So, it is necessary to start considering the creation of natural reserves with landscape corridors that include these microhabitats. Also, it is important to implement better reforestation and assisted regeneration practices, with germplasm of the area. Currently, it is still a common practice that authorities reforest the areas of Abies forests with other conifer species, such as Pinus ayacahuite, P. pseudostrobus, P. patula or Cupressus lusitanica, which in turn can change the microclimate and the environmental and habitat conditions that the Abies forest generates, as high humidity and low temperatures, characteristics that the high-mountain Ambystomas are adapted to (Sunny et al. 2014b; Monroy-Vilchis et al. 2015; Lemos-Espinal et al. 2016). This is important to highlight since, with the recent change of protection level of the NTV forest, harvesting practices are being allowed in almost all Abies forest extensions (Mastretta-Yanes et al. 2014; González-Fernández et al. 2019). Likewise, it is necessary to stop the influx of sewage and decrease waste from local agriculture and pollutants into the rivers, avoid the excess of free livestock, grassland fires, restrict ecotourism activities and prohibit the use of Ambystoma salamanders as traditional medicine for the lungs or traditional food (Monroy-Vilchis et al. 2015). It is also necessary to have more control over trout breeding because trout escape from farms to the streams, primarily due to low-quality fences and meshes dividing the farms from streams. We recommend divisions to be built higher and periodically revised to avoid holes from which trout can escape. Trout farms are affecting the populations of mole salamanders; several studies have demonstrated that the predation by introduced fish has been linked to reductions in the population’s size, survival, growth, egg predation, reproduction and organic waste pollution of the rivers (Valiente 2006; Mercado-Silva et al. 2008; Drake et al. 2014; Kenison et al. 2016) since mole salamanders do not recognize trout as predators and their survival decreases proportionately (Petranka et al. 1987; Petren et al. 1993; Relyea 2001; Pilliod et al. 2013). Trout farms have been encouraged by Mexican governmental agencies and recommended as a potential conservation tool for native forests (López-García et al. 2014; Zamora et al. 2018), so the best strategies may include isolating trout farms from streams containing native amphibians, increasing efforts to prevent trout escaping and attempting to eradicate populations of escaped trout from streams, thereby balancing the economic and conservation value of trout farms with their potential negative effects on native amphibians (Zamora et al. 2018). Also, it is necessary that each trout farm has an area of wastewater treatment since these waters when released to the stream can change the physical, chemical and bacteriological parameters and contaminate the stream for up to 3–12 km (Boaventura et al. 1997; Johnson and Speare 2005). Certain sediments turn completely anoxic and, with hydrogen sulfide by-products of microbial sulfate reduction, can increase to toxic levels (Valdemarsen et al. 2009, 2010, 2015). Also, it is very important that management agencies restore a fishless state in all lentic habitats that have the potential to host mole salamanders (Pearson 2004). Finally, additional research is needed, as well as effective communication between scientists and managers and sensible management actions to carry out a coordinated strategy between scientists, managers, local people and fishermen to ensure the long-term conservation of the mole salamanders. If these conservation strategies are not carried out, the expected heterozygosity values found will continue to decline, inbreeding will be higher, the effective population size will be lower and the mole salamanders will not have enough genetic diversity to be able to adapt to changes in the environment (Lande 1988; DeYoung and Honeycutt 2005; Frankham et al. 2005; Hedrick 2011).