A Salamander tale: Relative abundance, morphometrics and microhabitat of the critically endangered Mexican salamander Pseudoeurycea robertsi (Taylor, 1939)

Roberts’ False Brook Salamander (Pseudoeurycea robertsi) is a critically endangered plethodontid salamander, endemic to the Nevado de Toluca Volcano (NTV), Mexico. Little is known about the biology and ecology of this species, including its microhabitats. Thus, this study aimed to collect basic information about P. robertsi. We sampled fourteen forested sites in the NTV; to corroborate the correct identification of the species we used genetic data, we assessed the variation in head morphometric measurements and dorsal colouration patterns amongst localities and the microhabitat features associated with P. robertsi presence. Of the four potential salamander species, P. robertsi was the most abundant (89.80%) and widely distributed (approximately within 130 km2) salamander in the NTV. We did not find significant variations in morphometry; however, we found significant differences in dorsal patterns between populations (in the number and size of segments of the dorsal stripe). The average total length for 185 adults was 89.15 mm (38.7–117.9 mm); we found seven patterns of dorsal stripe. We found 98% of P. robertsi individuals under the bark of fallen logs in Abies religiosa and A. religiosa-Pinus sp. forests, with a higher number of detected salamanders in naturally-fallen logs than in cut logs (34% vs. 10%). Thus, keeping well-preserved A. religiosa forests and retaining fallen logs is essential to P. robertsi conservation.


Introduction
Fifty-one percent of Mexican amphibian species are at risk of extinction due to anthropogenic factors (NOM-059, SEMARNAT 2019) and plethodontids are especially vulnerable because many species occur only in small ranges, with isolated and fragmented populations and require very specific microhabitat characteristics (Wake and Lynch 1976;Wake 1987;IUCN 2020). According to the IUCN (2020), 97% of Mexican plethodontid species are in a risk category, with 44 species Endangered and 49 species Critically Endangered, placing Mexico as the country with the highest number of threatened plethodontid species (García-Bañuelos 2019;IUCN 2020). Moreover, 81% of the species of this group are endemic to Mexico (Frost 2019;AmphibiaWeb 2020). Salamanders are important in forest ecosystems as they represent an important proportion of the vertebrate biomass in oldgrowth forests (Davic and Welsh 2004); they are the most abundant vertebrates in North American forest ecosystems (Pough et al. 1987;Best and Welsh 2014) and are key to ecosystem functions as they feed on many invertebrate species and have an important indirect regulatory role in the processing of detritus-litter by ingestion of detritivore prey. Likewise, they may be important prey of other vertebrates (Davic and Welsh 2004).
Roberts' False Brook Salamander (Pseudoeurycea robertsi, Taylor, 1939) is a critically-endangered plethodontid salamander, endemic to forests of the Nevado de Toluca Volcano (NTV) and its surroundings (IUCN SSC Amphibian Specialist Group 2016). It is one of the most threatened amphibians in the country (Wilson et al. 2013); however, almost nothing is known about the ecological and morphological aspects of this species. Until our team began to study the species, only Bille (2009) had provided some important information about this salamander; he reported the presence of P. robertsi in four locations of the NTV and described three different dorsal patterns. In addition, a lot of controversy has arisen regarding the NTV, as its protection status has been recently changed from National Park to a less restrictive category that allows forest harvesting in most of Sacred Fir (Abies religiosa) forests of the natural protected area (Depraz et al. 2017;González-Fernández in prep.). Here, we aim to provide basic information (including relative abundance, morphometric measurements, variations in colouration and dorsal patterns and microhabitat use) and develop conservation strategies for this critically-endangered salamander which inhabits a natural protected area and its surroundings, subjected to controversial management decisions.

Study site
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).

Study species
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).

Pilot study
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.

Main 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.

Genetic sampling
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  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 2 , 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 SPLIT-STREE 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.

Morphometric measurements
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 10 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).

Variation in dorsal pattern
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 ANO-VA 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

Micro-habitat characterisation
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 (χ 2 ) 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 gg-plot2 R package (Wickham et al. 2019).

Results
We detected and measured 185 P. robertsi individuals (137 adults and 48 juveniles; Suppl. material 1: Table S1), resulting in the largest sample reported for this species. Individuals were found in 12 of the 14 sites, eight of which were new localities. In two sampling sites (Huacal Viejo-Agua Bendita and El Contadero), we found no salamanders. The number of P. robertsi within each site ranged from 0 to 46 (mean ± SD = 13.2 ± 12.7 individuals; Figs 1, 2; Suppl. material 1: Table S2). We sampled 873 logs and stumps (mean ± SD = 62.36 ± 24.08 logs, ranging from 16 logs in El Contadero to 103 in Palo Seco, Suppl. material 1: Table S3). The altitudinal range of our sampling was from 2800 to 3460 m a.s.l.

Sampling and species identification
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).
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

Morphometric measurements
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).

Variation in dorsal pattern
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 2 (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.
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 2 = -0.47, P = 0.09).

Discussion
We estimate that P. robertsi has a wider distribution than previously reported (8 km 2 , IUCN SSC Amphibian Specialist Group 2016) as we found the species in 12 of the 14 sampling sites, which were distributed throughout the A. religiosa forests of the NTV. Our data suggest that the species is present in most of the well preserved oldgrown Abies forests of the NTV, which occupy approximately 130 km 2 (González-Fernández et al., in prep), which still represents a very restricted distribution. Almost all individuals were found under the loose bark of  A. religiosa and Pinus sp. logs, a microhabitat formed when the bark separates from the rest of the log during the wood degradation process, creating a hollow space that provides shelter and invertebrate prey (Bille 2009). Although some individuals can be found in stumps, inside decaying logs, under the logs or in other microhabitat types (for example during the pilot study we found one individual under a rock), the preferred microhabitat by the species during the wet season is clearly under the bark of fallen logs. The sites Amanalco A and Meson Viejo had the highest percentages of naturally-fallen trees (100% and 94% from a total of 63 and 52 logs, respectively, Suppl. material 1: Table S3). We found more individuals in naturally-fallen logs than in cut logs (although cut logs were more abundant), which may be because naturally-fallen logs provide more suitable microhabitat due to their often-larger size and/or more extensive surface area at the complex break point that facilitates faster decomposition. In sites with moderate amount of logging, without clearings, there were also relatively high salamander abundances. However, we found no salamanders in forests with many clearings like El Contadero and Huacal Viejo-Agua Bendita. We found no higher use of A. religiosa logs over Pinus sp. logs as a microhabitat; however, we could not find salamanders in Pinus sp. only forests, probably because they are more open, with higher solar radiation and lower humidity levels. We also found records of the species at a lower altitude (2800 m a.s.l.) than previously reported (2900-3600 m a.s.l.; Bille 2009; IUCN SSC Amphibian Specialist Group 2016). At this lower altitude site (Agua Bendita), we also found nine individuals of A. cephalica. Our results show that, although four plethodontids inhabit the NTV, P. robertsi is the most widely-distributed species with the highest abundance. Pseudoeurycea robertsi and P. leprosa are the most morphologically-similar salamanders in the NTV (Bille 2009), thus interspecific competition may be higher between these species. Generally, niche overlap between similar species can increase interspecific competition (Navarro et al. 2013;Dehnhard et al. 2020). However, in communities of plethodontid salamanders in North America, morphological and ecological differences between closely-related sympatric species are associated with differences in resource use (Adams and Rohlf 2000;Adams 2010).
Geometric morphometrics of the head showed no significant differences between P. robertsi individuals and sampling sites which may suggest that there have not been differences in food resource use in the long term, unlike other studies that have found variations in head shape related to food resource use in other salamander species (Walls et al. 1993;Maerz et al. 2006;Silva-González et al. 2011).
Due to the great variety in dorsal patterns found so far, considering the differences in number, size, colour and shape of the segments of the dorsal stripe, we speculate that some other patterns may exist that cannot be classified in any of the seven categories described here (Fig. 6). Differences in dorsal colouration patterns and skin colour in amphibians have been reported as a response to different internal or external factors, such as substrate colour (colouration patterns similar to the colouration of the substrate as an anti-predatory response), sexual dimorphism, reproductive state or environmental temperature (Garcia and Sih 2003;Venesky and Anthony 2007;Rudh and Qvarnström 2013). Some of these causes may explain the differences in colouration patterns we observed; however, more studies should be undertaken regarding this topic.

Ecological and conservation implications
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) . 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 , 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.