The Nevado de Toluca Volcano (NTV; 19°09'N, 99°45'W) is located in the central-south part of the Trans-Mexican Volcanic Belt (TMVB), about 23 km southwest from Toluca, the capital of the State of Mexico. The southern flank of the Volcano is relatively flat and with multiple river valleys; the north, east and west are made up of structures of recent volcanic activity. The climate in NTV forests is a temperate semi-cold climate, with precipitation during the summer and an average annual temperature of 10 °C and an annual rainfall of 1186 mm (García 2004; Soto et al. 2020). The dominant land cover types at the Volcano are old-growth and secondary forest patches of Sacred Fir (A. religiosa) that grow between 2800–3500 m a.s.l. forming sky-islands in the TMVB (Rzedowski 2006; Mastretta-Yanes et al. 2015) and pines (Pinus hartwegii, Pinus montezumae and P. pseudostrobus). In these forests, naturally-fallen logs are abundant, especially in A. religiosa forests, which are in steep areas subjected to wind and storms. There are also some broad-leaved forests (Quercus and Alnus) at lower elevations in the eastern part of the Volcano and alpine grasslands at the highest elevations. There is an important extension of agricultural lands at the north-eastern part of the Volcano, including cattle pastures (Franco-Maass et al. 2006) and human settlements (Toscana-Aparicio and Granados-Ramírez 2015). The NTV is a priority terrestrial region for biodiversity conservation (Arriaga et al. 2000) and was declared a National Park in 1936. However, in 2013, the Mexican government changed the protection status of the NTV from National Park (a highly restrictive category in terms of resource use) to a less restrictive protection category (Flora and Fauna Protection Area, DOF 2013) – a controversial decision that could lead to further degradation of the last well-preserved Abies forests (Mastretta-Yanes et al. 2014). The NTV is impacted by several human activities of varying intensity: legal and illegal logging, where drug cartels are also involved, commercial export-orientated floriculture, high chemical input due to potato production and mineral extraction, leading to habitat loss and fragmentation (Orozco-Hernández 2007; Toscana-Aparicio and Granados-Ramírez 2015; Depraz et al. 2017).

Pseudoeurycea robertsi (Taylor, 1939) is a plethodontid salamander closely related to P. altamontana (Taylor, 1939) and P. longicauda (Lynch, Wake & Yang, 1983) (Parra-Olea 2002). Taylor (1938) reported snout-vent lengths of 35–51 mm for females and 49 mm for a single male. The tail is laterally compressed and almost equal to snout-vent length or shorter. The head is broad, rather flattened and with truncate snout. The limbs are well developed and without interdigital membrane as this is a terrestrial salamander. The first digit is very short. Bille (2009) reported three dorsal colouration patterns: a well-defined dorsal stripe, dorsal mottling and almost without pattern with few scattered spots.

The study species reproduces by direct development, lives in Pinus-Abies forests at elevations between 2900–3600 m a.s.l. and can be found under rocks, logs and loose bark of fallen logs and stumps. It has not been found in severely-disturbed habitats and can tolerate only very slight selective logging (Taylor 1938; Bille 2009; IUCN SSC Amphibian Specialist Group 2016). It was a relatively common species, but due to ongoing decline in the extent and quality of its habitat, the population is suspected to be decreasing (Taylor 1938; Bille 2009; IUCN SSC Amphibian Specialist Group 2016).

In literature, it is reported that P. robertsi is present in Pinus sp.-A. religiosa forests (IUCN SSC Amphibian Specialist Group 2016), thus we looked for the species in forests dominated by both Pinus sp., A. religiosa and mixed forests of Pinus sp. and A. religiosa. To identify the species in the field, we used the phenotypic characteristics described by Bille (2009) and field guides (Flores-Villela 1993; Ramírez-Bautista et al. 2009). We also obtained genetic data to corroborate the correct identification of the species, since Bille (2009) reported that P. robertsi was morphologically similar to P. leprosa (Cope 1869), also present in the NTV and has similar dorsal patterns, mainly in low-altitude areas (Bille 2009).

As there is no published information about abundance and/or seasonality for P. robertsi, we carried out a pilot survey for plethodontids in 2015, carrying out monthly 2–3-day surveys per sampling site from the beginning of the rainy season (late April) until the beginning of the dry season (late October). We found most of the salamanders in fallen logs and stumps in A. religiosa forest and mixed forests (under the bark of fallen trees, inside decaying fallen logs, under logs and under the bark of stumps), thus we focused on these microhabitats for the formal abundance survey. Moreover, it is easier to focus on fallen logs than to survey the entire forest floor. We did not find any individuals in Pinus sp. forests. In April and October, we did not detect any plethodontid salamanders. During May and September, we found an increasing and decreasing abundance of plethodontid salamanders, respectively (for example, at a site where we found around 20 individuals in each visit in June, July and August, we hardly found 3–4 in May and September and none in April and October). Thus, to minimise seasonality effects, the main study was based on a survey carried out in mid-June to early August, which are the wettest months of the year when we found the highest numbers of individuals in the pilot study.

From mid-June to early August 2016, we carried out diurnal surveys for terrestrial salamanders, focusing on P. robertsi individuals in 14 sampling sites of 10 ha (an extension large enough to find enough fallen logs) distributed across the volcano, between 2850 and 3450 m a.s.l. (Fig. 1). Site selection was arbitrary, but trying to cover the entire gradient of A. religiosa-Pinus sp. forest. We did not carry out night-time sampling due to safety concerns. Linear transects are difficult to perform due to the density of the forests and, therefore, we standardised our sampling effort based on time. Each site was visited once and sampled by two people for two effective hours of sampling (we excluded the time we spent handling individuals and measuring the fallen logs and stumps); thus, fieldwork was from 09:00 to 16:00 h. For each individual, we collected morphometric data, took a tissue sample and characterised the microhabitat where it was found (see subsections below). After measurements were completed, individuals were released back to their corresponding capture sites.

We removed 2 mm of the tail tip of adult salamanders for DNA extraction. Given that Pseudoeurycea robertsi exhibits caudal autotomy as a defence mechanism, the removal of only a small portion of the distal part of the tail reduces the disadvantages associated with caudal autotomy (Romano et al. 2010). Moreover, to avoid spreading pathogens between individuals (e.g. Batrachochytrium dendrobatidis), we used new latex gloves when handling each individual. Tissues were preserved in 90% ethanol in the field and then frozen in the laboratory within the same day at -20 °C until processed. We extracted DNA with the GF-1 nucleic acid extraction kit (Vivantis) following the manufacturer’s instructions and used it as the template for amplification of the mitochondrial cytochrome b gene (cyt b) with the following primers: MVZ15 (5´G A A C T A A T G G C C C A C A C Q Q T A C G N A A -3´) and MVZ16 (5´A A A T A G G A A R T A T C A Y T C T G G T T T R A T -3) from Moritz et al. (1992). PCR amplifications were performed in a total volume of 15 μl. We used 0.15 μl of 5U Taq polymerase (QIAGEN), 0.75 μl of 10X PCR buffer with 2.5 mM MgCl₂, 0.3 mM of deoxynucleotide triphosphates (dNTPs) and 0.3 mM of forward and reverse primers. The PCR conditions were as follows: initial 4 min denaturation at 94 °C, followed by 35 cycles, each cycle consisting of 94 °C denaturing for 1 min, 50 °C annealing temperature for 1 min and extension at 72 °C for 2 min and finally at 72 °C for 3 min. Purification and pair-end sequencing of DNA was performed by the UW High Throughput Genomics Center of the CINVESTAV Irapuato, Mexico. The forward and reverse sequences for each individual were aligned and edited manually using BIOEDIT 7.1.3 (Hall 1999). The retrieved sequences of cyt b were blasted in GenBank (Morgulis et al. 2008) and optimised for the 100 most similar sequences (Megablast) to confirm they matched to P. robertsi, considering the E value and the percentage of identity. Likewise, we constructed rooted network using the Neighbour-Net algorithm with SPLITSTREE 4.6 (Huson and Bryant 2006), based on the patristic distance corrected by the GTR+I+G model of evolution to observe if our samples were grouped with P. robertsi or P. leprosa individuals; as well, a sample of P. altamontana was used to root the haplotype network, as this species is the most phylogenetically related to P. robertsi (Parra-Olea 2002). The sequences were deposited in the GenBank database with access numbers: MK357639–MK357709.

All tentative P. robertsi individuals were photographed with a 20.2 megapixel digital camera Canon Powershot ELPH 360 HS. We used a tripod to take numerous pictures from a height of 20 cm from the millimetre paper on which individuals were placed ventrally. We obtained morphometric measurements for each individual, analysing more than one picture if necessary, by using the measure tool in software FIJI 1.53c (Schindelin et al. 2012). For each picture, we defined the pixel size, based on the millimetre paper. We recorded the following morphometric measurements: head length (HL), from the tip of the snout to the neck; head width (HW), across the widest point of the head; left eye diameter (LED); width between eyes (WBE); median body width (MBW), across the trunk midway between the front and hind limb insertions; posterior femur length (PFL); snout vent length (SVL); tail length (TLO); tail width (TW); total length (TL), from the posterior left limb insertion to the tip of the longest outstretched toe. We did not sex individuals as only external measurements do not discriminate sufficiently between males and females (Ollivier and Welsh 2003) and we decided against invasive/lethal techniques. We identified age class based on TL; individuals were classified as juvenile or adult, according to the Sturges´ Rule (Lindström et al. 2010), which is amongst the most common strategies to determine the number of classes in a frequency histogram (Legendre and Legendre 2003), by approximating the resulting histogram to a normal distribution (Hyndman 1995). The number of classes of a histogram, using Sturges´ Rule, is determined by the expression k = 1 + 3.3 × log¹⁰ m; where: k = the number of classes and m = the number of observations in the dataset (Legendre and Legendre 2003). This Rule has proven a satisfactory strategy to tackle the problem of “over smooth” or “discontinuous” histograms resulting from an inadequate determination of the number of classes, when “n” is less than 200 (Hyndman 1995). Several studies have used this approach to analyse ecological problems (Magalhaes et al. 2002; Schmidth et al. 2009; Diallo et al. 2012; Sunny et al. 2014).

For geometric morphometrics of the head, the software MakeFan 6 (Sheets 2003) was used to set the support points and fans and the tpsDig 2 (Rohlf 2010) was used to digitise landmark coordinates. All morphometric measurements were made in software MorphoJ 1.06 (Klingenberg 2011). Twenty-three landmarks were digitised (Suppl. material 1: Fig. S1): 11 to contour of the head, three for each eye, two as support points for the layout of fans that covered the entire head and four landmarks in the central area of the head (at the posterior tips of the nasal and the frontal bones and at the anterior and posterior of the interparietal bone).

The shape information was extracted using a Procrustes superposition analysis, which eliminates the effect of size, position and orientation to standardise each specimen according to centroid size (Rohlf and Slice 1990). To characterise head shape variation amongst the sampling localities, a principal component analysis (PCA) was carried out, based on the covariance matrix of shape of the datasets of dorsal view (Vilaseca et al. 2020). In addition, we visualised the variations in the configuration through vectors and deformation grids of each sampled site, which were built by applying the thin plate spline interpolation function (Bookstein 1991).

Based on the photographs, converted into 8-bit images in the RGB colour space (a total of 256 colours), we calculated the colour frequency of the dorsal stripe of each individual and at each sampling site using the 3D Colour Inspector and Colour Histogram plug-ins (Barthel 2007) of FIJI. We calculated the number of segments of the dorsal stripe and their area. We carried out a two-way ANOVA and the post hoc Least Significant Difference test (LSD) to see if there were significant differences between number and area of the segments of the back stripe and between locations, using the package DescTools (Signorell et al. 2016) for R (version 3.4.0; R Development Core Team 2017).

We sampled fallen logs and stumps with a diameter > 5 cm and length > 30 cm. For each log/stump, we measured the relative humidity under the bark of the logs of the A. religiosa forest and under the bark of the logs deposited in grasslands during the pilot study with a thermohygrometer (STEREN TER-150). We recorded the geographic coordinates and altitude for the sampling sites and each salamander observation with a hand-held GPS (Garmin etrex 20, with precision of 3 m, altitude being only measured once for each sampling site). For each observed salamander, we categorised the following microhabitat types: tree species, if the log or stump had fallen naturally or was cut down, its length and width and the number of salamanders’ present. All log measurements were taken with a 10-m flexible measuring tape. We summed the number of salamanders and estimated overall logs and stumps volume in each sample site. We applied a Pearson correlation to test if elevation was significantly correlated to the number of salamanders. We applied t-tests to assess if there were significant differences between the length and diameter of logs with and without salamanders. Chi-squared (χ²) was performed to test for an association between tree species and presence of salamanders and between the cause of tree fall (human or natural) and presence of salamanders. We performed an ANOVA to test significant differences in the length and diameter of logs per site. Statistical significance for all tests was set at p-value < 0.05. All statistical tests were performed using R and the graphs were made using ggplot2 R package (Wickham et al. 2019).

The species identification based on cyt b confirmed that, considering all sampling sites, 89% of all salamanders found were P. robertsi and 11% were P. leprosa with an E value of 0.0 and an identity percentage of 98% considering a sequence of 617 base pairs. The Haplotype network analysis showed a clear distance between the clades belonging to P. robertsi, P. altamontana (the most closely related to P. robertsi species, Parra-Olea 2002) and P. leprosa (Fig. 3).

Figure 3. 

Haplotype network obtained with the Neighbor-Net algorithm. P. robertsi clusters are shown in green, P. leprosa are shown in blue and P. altamontana are shown in orange. Scale bar are substitutions per site.

We found P. robertsi in A. religiosa forests and in mixed forests. Only four salamanders were found in stumps; all other individuals were found under the bark of fallen logs. We found nine individuals of Aquiloeurycea cephalica (Cope 1865) (all in Agua Bendita), eight individuals of P. leprosa (five in Palo Seco, two in Carretera and one in Agua Bendita) and four individuals of Isthmura bellii (Gray 1850) (two in Amanalco A, one in Mesón Viejo and one in Rancho Viejo). Pseudoeurycea robertsi was present at all sites where we found other salamander species; sites with no P. robertsi revealed no other species. Agua Bendita was the only sampling site where we found three species of salamanders (Fig. 2; Suppl. material 1: Table S2). On 7 Aug 2016, in Palo Seco, we found a P. robertsi female (TL = 91.1 mm) attending an egg mass of 32 white eggs, under the bark of an A. religiosa log.

We analysed 135 P. robertsi adult individuals from 11 localities: five from Agua Bendita, 25 from Amanalco A, 18 from Amanalco B, eight from Amanalco C, 10 from Carretera, 13 from Las Lagrimas, 16 from Mesón Viejo, 11 from Palo Seco, 18 from Rancho Viejo, three from Raíces and eight from Santa Cruz. Juvenile P. robertsi had a TL = 33.30 mm and adults – TL = 89.15 mm (Table 1). The PCA of the head measurements showed that the three PCs accounted for 42.5% (PC1: 20.9%, PC2: 11.9%, PC3: 9.8%) of the variation in head shape. All individuals of P. robertsi from all sampling sites were grouped into a single group by the PCA (Fig. 4). The deformation grids of the landmarks showed a greater deformation in the nose tip and in the central part of the head, which were displaced posteriorly (Fig. 5). However, according to the PCA analysis, the shape of the dorsal head view amongst sampling sites of P. robertsi was homogeneous (Fig. 4).

Figure 4. 

PCA analysis of head dimensions of adult P. robertsi individuals.

Figure 5. 

Deformation grids of the head landmarks of adult P. robertsi individuals.

The variation in dorsal patterns was high between individuals and we found significant differences in dorsal patterns between some populations. Although the number of segments on the dorsal stripe varied overall from one to three (mean ± SD = 1.5 ± 0.59), it was consistent within some sites; for example, individuals from Amanalco A, Amanalco B and Carretera only had three segments. We found significant differences between the sampling sites (F = 243.429, DF = 10, P = 0) and after applying the LSD test, we found significant differences amongst the sites of Las Lagrimas-Agua Bendita, Meson Viejo-Amanalco A, Santa Cruz-Amanalco A, Meson Viejo-Amanalco B, Las Lagrimas-Amanalco C, Las Lagrimas-Carretera, Meson Viejo-Las Lagrimas, Santa Cruz-Las Lagrimas, Rancho Viejo-Meson Viejo and Santa Cruz-Rancho Viejo. The area of the dorsal stripe segments varied from 1.4 to 3.1 mm² (mean ± SD = 3.0 ± 1.63). We found significant differences amongst the sampling sites (F = 4.407, DF = 10, P = 0) and after applying the LSD test, we found significant differences between Amanalco C-Amanalco A and Meson Viejo-Amanalco A.

The colour of the dorsal stripes ranged from red to orange (Fig. 6A, B) and from brown to black (Fig. 6C). We found seven patterns of dorsal stripes (Fig. 6), two of them (Fig. 6.1, 6.2) described by Bille (2009); we could not find the third pattern (all black with few scattered spots). The patterns that we found were: 1) well-defined, brick red dorsal stripe (Fig. 6.1), 2) dense dorsal mottling (Fig. 6.2), 3) dorsal red line only present in the tail (Fig. 6.3), 4) semi- well-defined, yellow dorsal stripe (Fig. 6.4), 5) dorsal yellow line only present in the tail (Fig. 6.5), 6) few scattered spots in the dorsal line and in the tail (Fig. 6.6) and 7) all black without spots or dorsal line (Fig. 6.7). The proportions of the patterns per site are provided in Fig. 6D and Table 2.

Figure 6. 

Dorsal stripe coloration types in P. robertsi (A–C) and proportion of colors per site, ranging from red to orange (A, B) and from brown to black (C), different dorsal patterns found: 1) well-defined, brick red dorsal stripe; 2) dense dorsal mottling; 3) dorsal red line only present on the tail; 4) semi well-defined, yellow dorsal stripe; 5) dorsal yellow line only present on the tail; 6) few scattered spots in the dorsal line and in the tail; 7) all black without spots or dorsal line; D) color patterns of the dorsal stripe in the different sampling sites.

A total of 873 fallen logs and stumps (667 A. religiosa logs and 206 Pinus sp. logs, Fig. 7A; Suppl. material 1: Table S3) were measured. We found salamanders in 41 (20%) Pinus sp. logs and 121 (18%) A. religiosa logs. We did not find an association between log species and presence of salamanders (χ² = 0.22, DF = 1, P = 0.64). However, we found a highly significant positive association between fallen/cut logs and the presence of salamanders (χ² = 13.47, DF = 1, P = 0.00); the proportion of salamanders was higher in naturally-fallen logs than in cut logs (34% vs. 10%). The sites with the highest percentage of naturally-fallen trees were Amanalco A (100%) and Meson Viejo (94%), while the sites with the most felled logs were Huacal Viejo-Agua Bendita (100%), El Contadero (100%), Santa Cruz (100%), Palo Seco (100%) and Carretera (98%; Fig. 7B; Suppl. material 1: Table S4). The logs’ volume ranged from 1.3 to 19.6 m³; the largest volume of logs were in Las Lagrimas (19.63 m³), Palo Seco (12.96 m³), Meson Viejo (11.68 m³), Raices (11.59 m³) and Rancho Viejo (11.07 m³; Fig. 7C; Suppl. material 1: Table S5). We found significant differences in the average length of fallen logs per sampling site (F (13, 881) = 8.05, P = 0.00; Fig. 8A; Suppl. material 1: Table S6) and diameter of fallen logs per sampling site (F (13, 857) = 3.56, P = 0.00; Fig. 8B; Suppl. material 1: Table S7). Salamanders were significantly more abundant in longer logs (mean log length with salamanders = 4.2 m vs. 2.7 m without: t = 5.4, DF = 789, P = 0; raw data in Suppl. material 1: Table S8). We found no significant differences between the diameter of logs with and without salamanders (t = 1.02, df = 789, P = 0.31; raw data in Suppl. material 1: Table S9). The average length of the logs ranged from 1 m to 4 m, the average diameter ranged from 0.30 m to 0.42 m (Suppl. material 1: Table S10). The sampling sites that presented the highest percentage of A. religiosa logs were Raices (96%), Palo Seco (96%) and Amanalco A (91%), while the highest percentage of Pinus sp. was in Meson Viejo (69%; Suppl. material 1: Table S11).

Figure 7. 

Characteristics of the logs per site: A) percentage of tree species; B) percentage of trees naturally fallen, and logs felled; C) volume of logs per site.

The relative humidity in the microhabitats of A. religiosa forest varied from 84 to 99% even during the warmest time. It is important to mention that the logs we found at the forest edges were completely dry and tough and those inside the forests were wet and decomposing. We only found salamanders in wet logs with intermediate levels of decomposition. Finally, altitude and number of P. robertsi individuals did not correlate significantly (R² = -0.47, P = 0.09).

Figure 8. 

Average length (A) and average diameter (B) of logs per site.

Our results suggest that A. religiosa forest and mixed forests with naturally-fallen logs hold high salamander abundance in the NTV. Alterations in these forests like logging can have negative impacts on this endemic and critically-endangered salamander, decreasing its abundance and even leading the species to extinction. Even selective logging, which is usually considered preferable to deforestation to maintain natural ecosystem functions, services or biodiversity can also be detrimental for this forest specialist species, as this kind of logging extends more deeply into the interior core of remaining forest areas (Broadbent et al. 2008) and can alter the microclimatic conditions of the forest. Drier forest conditions reduce wood decomposition rate (Crockatt and Bebber 2015), which may affect salamander refuges. We have found that the main drivers for the number of P. robertsi individuals were logs volume (positive effect) and edge density (negative effect) (González-Fernández et al. 2019). In 2013 (DOF 2013), the Mexican government changed the protection status of the NTV from National Park to a less restrictive protection category (Flora and Fauna Protection Area) allowing forest harvesting practices with commercial purposes in more than a half of the A. religiosa forest (Mastretta-Yanes et al. 2014, González-Fernández et al., in prep). This change in protection status may negatively impact the habitat of the P. robertsi and the other forest-specialist salamanders that are also present in the NTV. Moreover, A. religiosa forests, including those in the NTV, receive the migratory colonies of the Monarch Butterfly (Danaus plexippus; Saunders 2018). Therefore, A. religiosa forest preservation must be a priority in the NTV, especially in the areas of Amanalco A, B and C, Meson Viejo and Rancho Viejo, since they have higher levels of genetic diversity (Sunny et al. 2019), higher abundance of P. robertsi and contain the best preserved old-growth A. religiosa and mixed forests within the NTV. Therefore, we suggest key steps necessary for the future preservation of these forests: 1) Banning commercial logging, but allowing a subsistence and low-intensity forest use by local communities, 2) keeping fallen logs in the forests and including this measure within the management plans of the area 3) avoiding reforestation practices with nursery trees, which have limited genetic diversity and usually come from other localities. Local seeds can be used instead when reforestation practices are necessary; however, in most of the natural protected area, reforestation is carried out very close to the forest where regeneration can easily occur without human intervention, especially because A. religiosa forests have high regeneration capacity (Franco-Maass 2006) and 4) finally, implementing environmental education highlighting the value of the natural forests, not only to local communities, but also to the authorities to avoid inadequate forest management practices and practices that are only focused on obtaining economic benefits.