Print
A Salamander tale: Relative abundance, morphometrics and microhabitat of the critically endangered Mexican salamander Pseudoeurycea robertsi (Taylor, 1939)
expand article infoArmando Sunny, Hublester Domínguez-Vega§, Carmen Caballero-Viñas, Fabiola Ramírez-Corona|, Marco Suárez-Atilano, Andrea González-Fernández#
‡ Universidad Autónoma del Estado de México, Toluca, Mexico
§ Universidad Intercultural del Estado de México, San Felipe del Progreso, Mexico
| Universidad Nacional Autónoma de México, Ciudad de México, Mexico
¶ University of California, Berkeley, United States of America
# Universidad Autónoma Metropolitana-Unidad, Lerma, Mexico
Open Access

Abstract

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.

Key Words

Abies religiosa forest, amphibians, conservation, deforestation, endangered species, Nevado de Toluca Volcano, old-growth forest, Pinus sp. forest, Plethodontidae

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 old-growth 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.

Materials and methods

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 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 MgCl2, 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: MK357639MK357709.

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 × log10 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 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).

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

Figure 1. 

Locations of the 14 study sites in the Nevado de Toluca Volcano, located in the Trans-Mexican Volcanic Belt (dark-grey on the inset map), Mexico. The blue number are the P. robertsi individuals found in each sampling site.

Figure 2. 

Number of individuals per species in each sampling site.

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

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.

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

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.

Table 1.

Mean linear measurements (in mm) of juvenile and adult P. robertsi. Measurements are: 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.

Juvenile (N = 48) Adult (N = 137)
HL 5.48 10.80
HW 4.15 7.29
LED 1.08 2.22
WBE 1.80 2.71
MBW 4.03 7.24
PFL 1.75 3.94
SVL 18.54 45.41
TLO 14.76 43.74
TW 1.98 3.81
TL 33.30 89.15

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 mm2 (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.

Table 2.

Dorsal pattern and color of P. robertsi in each sampling site (sorted by elevation). Coloration pattern letters (A–C, 1–7) refer to Fig. 6.

Site Altitude (m a.s.l.) Number of segments on the dorsal stripe Mean area of the segments on the dorsal stripe (mm2) Coloration on the dorsal stripe pattern Predominant color of the dorsal stripe
Agua Bendita 2 800–2 890 1–2 2.9 4, 5, 6 B, C
Rancho Viejo 2 880–2 900 1–2 1.9 4, 5, 6, 7 C, No dorsal stripe
Amanalco B 2 890–2 950 1–3 2.1 1, 4, 5 A, C
Amanalco C 2 930 1–2 2.5 2, 4, 7 C
Meson Viejo 2 960–2 990 1–2 3.1 2, 4, 5, 7 B, No dorsal stripe
Amanalco A 2 970–3 000 1–3 1.9 2, 3, 4, 5, 6 B, C
Las Lagrimas 2 982 1–2 1.4 1, 4, 5 B
El Contadero 3 088
Huacal Viejo-Agua Bendita 3 100–3 200
Santa Cruz 3 200–3 220 1–2 3 1, 4, 5 B
Palo Seco 3 230–3 240 1–2 1.9 4, 5, 6 B, C
Raices 3 235–3 287 1–2 1.9 2, 6 C
Carretera 3 300–3 325 1–3 2.3 1, 2, 3 A
San Juan de las Huertas 3 440–3 460 1–2 1.9 1, 2, 6 C

Micro-habitat characterisation

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 (χ2 = 0.22, DF = 1, P = 0.64). However, we found a highly significant positive association between fallen/cut logs and the presence of salamanders (χ2 = 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 m3; the largest volume of logs were in Las Lagrimas (19.63 m3), Palo Seco (12.96 m3), Meson Viejo (11.68 m3), Raices (11.59 m3) and Rancho Viejo (11.07 m3; 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 (R2 = -0.47, P = 0.09).

Figure 8. 

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

Discussion

We estimate that P. robertsi has a wider distribution than previously reported (8 km2, 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 old-grown Abies forests of the NTV, which occupy approximately 130 km2 (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) (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.

Acknowledgements

We are grateful to Yurii Kornilev and two anonymous reviewers for their comments. We appreciate the help received by Giovanny González-Desales in the editing of some photos and the help provided in determining the colours of the salamanders. This paper was completed while M.S.A was on his post-doctoral stay at the Museum of Vertebrate Zoology at UC Berkeley (UC MEXUS-CONACYT: 10327268) and A.G.F was on his post-doctoral stay at UAM Lerma (PODEP: 511-6/2019-15932). A.S received financial support from the Secretary of Research and Advanced Studies (SYEA) of the Autonomous University of the State of Mexico (Grant: 4732/2019CIB).

References

  • Abràmoff MD, Magalhães PJ, Ram SJ (2004) Image processing with ImageJ. Biophotonics International 11(7): 36–42.
  • Adams DC, Rohlf FJ (2000) Ecological character displacement in Plethodon: biomechanical differences found from a geometric morphometric study. Proceedings of the National Academy of Sciences 97(8): 4106–4111. https://doi.org/10.1073/pnas.97.8.4106
  • AmphibiaWeb (2020) AmphibiaWeb. University of California, Berkeley. https://amphibiaweb.org [Accessed 10 Nov 2020]
  • Arriaga L, Espinoza JM, Aguilar C, Martínez E, Gómez L, Loa E (2000) Regiones Terrestres Prioritarias de México. México. Comisión Nacional para el Conocimiento y Uso de la Biodiversidad (CONABIO), 609 pp.
  • Best ML, Welsh Jr HH (2014) The trophic role of a forest salamander: impacts on invertebrates, leaf litter retention, and the humification process. Ecosphere 5(2): 1–19. https://doi.org/10.1890/ES13-00302.1
  • Bille T (2009) Field observations on the salamanders (Caudata: Ambystomatidae, Plethodontidae) of Nevado de Toluca, Mexico. Salamandra 45(3): 155–164.
  • Bookstein FL (1991) Morphometrics Tools for Landmarks Data: Geometry and Biology. Cambridge University Press, Cambridge, 456 pp.
  • Broadbent EN, Asner GP, Keller M, Knapp DE, Oliveira PJ, Silva JN (2008) Forest fragmentation and edge effects from deforestation and selective logging in the Brazilian Amazon. Biological Conservation 141(7): 1745–1757. https://doi.org/10.1016/j.biocon.2008.04.024
  • Cope ED (1865) Third contribution to the herpetology of tropical America. Proceedings of the Academy of Natural Sciences of Philadelphia 17: 185–198.
  • Cope ED (1869) A review of the species of the Plethodontidae and Desmognathidae. Proceedings of the Academy of Natural Sciences of Philadelphia 21: 93–118. https://www.jstor.org/stable/4060103
  • Crockatt ME, Bebber DP (2015) Edge effects on moisture reduce wood decomposition rate in a temperate forest. Global Change Biology 21(2): 698–707. https://doi.org/10.1111/gcb.12676
  • Dehnhard N, Achurch H, Clarke J, Michel LN, Southwell C, Sumner MD, Eens M, Emmerson L (2020) High inter- and intraspecific niche overlap among three sympatrically breeding, closely related seabird species: Generalist foraging as an adaptation to a highly variable environment?. Journal of Animal Ecology 89(1): 104–119. https://doi.org/10.1111/1365-2656.13078
  • Depraz S, Sanial E, Catalán AKR, Rojas AS (2017) Less protection for better conservation? A politicised relationship between a city and its protected area in the vicinity of Nevado de Toluca (Mexico). Journal of Urban Research (16). https://doi.org/10.4000/articulo.3261
  • Diallo D, Sall AA, Buenemann M, Chen R, Faye O, Diagne CT, Faye O, Ba Y, Dia I, Watts D, Weaver SC (2012) Landscape ecology of sylvatic chikungunya virus and mosquito vectors in southeastern Senegal. PLOS Neglected Tropical Diseases 6(6): e1649. https://doi.org/10.1371/journal.pntd.0001649
  • DOF [Diario Oficial de la Federación Mexicana] (2013) Decreto que reforma, deroga y adiciona diversas disposiciones del diverso publicado el 25 de enero de 1936, por el que se declaró Parque Nacional la montaña denominada “Nevado de Toluca” que fue modificado por el diverso publicado el 19 de febrero de 1937. http://dof.gob.mx/nota_detalle_popup.php?codigo=5315889 [Accessed 5 Nov 2020]
  • Flores-Villela O (1993) Annotated list of the species of amphibians and reptiles of Mexico, recent taxonomic changes, and new species. Carnegie Museum of Natural History, Pittsburgh, 73 pp. https://doi.org/10.2307/1447111
  • Franco-Maass S, Regil-García HH, González-Esquivel C, Nava-Bernal G (2006) Cambio de uso del suelo y vegetación en el Parque Nacional Nevado de Toluca, México, en el periodo 1972–2000. Investigaciones geográficas 61: 38–57. https://doi.org/10.14350/rig.29996
  • Frías-Alvarez P, Zúniga-Vega JJ, Flores-Villela O (2010) A general assessment of the conservation status and decline trends of Mexican amphibians. Biodiversisty and Conservation 19(13): 3699–3742. https://doi.org/10.1007/s10531-010-9923-9
  • García E (2004) Modificaciones al sistema de clasificación climática de Köppen. Universidad Nacional Autónoma de México, 246 pp.
  • Garcia TS, Sih A (2003) Color change and color-dependent behavior in response to predation risk in the salamander sister species Ambystoma barbouri and Ambystoma texanum. Oecologia 137(1): 131–139. https://doi.org/10.1007/s00442-003-1314-4
  • García-Bañuelos P, Rovito SM, Pineda E (2019) Representation of threatened biodiversity in protected areas and identification of complementary areas for their conservation: Plethodontid Salamanders in Mexico. Tropical Conservation Science 12: e1940082919834156. https://doi.org/10.1177/1940082919834156
  • González-Fernández A, Arroyo-Rodríguez V, Ramírez-Corona F, Manjarrez J, Aguilera-Hernández A, Sunny A (2019) Local and landscape drivers of the number of individuals and genetic diversity of a microendemic and critically endangered salamander. Landscape Ecology 34(8): 1989–2000. https://doi.org/10.1007/s10980-019-00871-2
  • González-Fernández A, Segarra J, Sunny A, Couturier S (in prep.) Forest cover loss in the Nevado de Toluca volcano (México) after the change to a less restrictive category.
  • Gray JE (1850) Catalogue of the Specimens of Amphibia in the Collection of the British Museum. Part II. Batrachia Gradientia, etc. Printed by Order of the Trustees. Spottiswoodes and Shaw, London.
  • Hall TA (1999) BioEdit: a user-friendly biological sequence alignment editor and analysis program for Windows 95/98/NT. Nucleic Acids Symposium Series (Vol. 41, No. 41, pp. 95–98). [London]: Information Retrieval Ltd., c1979–c2000.
  • Hyndman RJ (1995) The problem with Sturges Rule for constructing histograms. Monash University, Melbourne.
  • IUCN (2020) The IUCN Red List of Threatened Species. Version 2020–1. [Accessed 14 Nov 2020]
  • Kozak KH, Wiens JJ (2016) What explains patterns of species richness? The relative importance of climatic-niche evolution, morphological evolution, and ecological limits in salamanders. Ecology and Evolution 6(16): 5940–5949. https://doi.org/10.1002/ece3.2301
  • Legendre P, Legendre L (2003) Numerical Ecology, Amsterdam, Netherlands, Elsevier.
  • Lindström J, Reeve R, Salvidio S (2010) Bayesian salamanders: analyzing the demography of an underground population of the European Plethodontid Speleomantes strinatii with state-space modelling. BMC Ecology 10: 1–4. https://doi.org/10.1186/1472-6785-10-4
  • Lynch JF, Wake DB, Yang SY (1983) Genic and morphological differentiation in Mexican Pseudoeurycea (Caudata: Plethodontidae). Copeia 1983: 884–894. https://doi.org/10.2307/1445090
  • Magalhaes MF, Batalha DC, Collares-Pereira MJ (2002) Gradients in stream fish assemblages across a Mediterranean landscape: contributions of environmental factors and spatial structure. Freshwater Biology 47(5): 1015–1031. https://doi.org/10.1046/j.1365-2427.2002.00830.x
  • Mastretta-Yanes A, Cao R, Nicasio-Arzeta S, Quadri P, Escalante-Espinosa T, Arredondo L, Piñero D (2014) ¿Será exitosa la estrategia del cambio de categoría para mantener la biodiversidad del Nevado de Toluca? Oikos 12: 7–17.
  • Mastretta-Yanes A, Moreno-Letelier A, Piñero D, Jorgensen TH, Emerson BC (2015) Biodiversity in the Mexican highlands and the interaction of geology, geography and climate within the Trans-Mexican Volcanic Belt. Journal of Biogeography 42(9): 1586–1600. https://doi.org/10.1111/jbi.12546
  • Moritz C, Schneider CJ, Wake DB (1992) Evolutionary relationships within the Ensatina eschscholtzii complex confirm the ring species interpretation. Systematic Biology 41(3): 273–291. https://doi.org/10.1093/sysbio/41.3.273
  • Ollivier L, Welsh Jr HH (2003) Determining sex and life stage of Del Norte Salamanders from external cues. Northwestern Naturalist 129–134. https://doi.org/10.2307/3536538
  • Ordoñez-Ifarraguerri A, Siliceo-Cantero HH, Suazo-Ortuño I, Alvarado-Díaz J (2017) Does a frog change its diet along a successional forest gradient? The case of the Shovel-Nosed Treefrog (Diaglena spatulata) in a tropical dry forest in Western Mexico. Journal of Herpetology 51(3): 411–416. https://doi.org/10.1670/16-024
  • Orozco Hernández ME (2007) Entre la competitividad local y la competitividad global: floricultura comercial en el Estado de México, Convergencia. Revista de Ciencias Sociales 14(45): 111–160.
  • Parra-Olea G (2002) Molecular phylogenetic relationships of neotropical salamanders of the genus Pseudoeurycea. Molecular Phylogenetics and Evolution 22(2): 234–246. https://doi.org/10.1006/mpev.2001.1048
  • Parra-Olea G, Flores-Villela O, Mendoza-Almeralla C (2014) Biodiversidad de anfibios en México. Revista Mexicana de Biodiversidad 85: 460–466. https://doi.org/10.7550/rmb.32027
  • Pough FH, Smith EM, Rhodes DH, Collazo A (1987) The abundance of salamanders in forest stands with different histories of disturbance. Forest Ecology and Management 20(1–2): 1–9. https://doi.org/10.1016/0378-1127(87)90146-0
  • R Development Core Team (2017) R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna. http://www.r-project.org [on May 18, 2019]
  • Ramírez-Bautista A, Hernández-Salinas U, García-Vázquez UO, Leyte-Manrique A, Canseco-Márquez L (2009) Herpetofauna del Valle de México: diversidad y conservación. Universidad Autónoma del Estado de Hidalgo, 213 pp.
  • Rohlf FJ (2010) TpsDig 2-thin plate spline digitizer. Ecology and Evolution, State University at Stony Brook, New York.
  • Rohlf FJ, Slice D (1990) Extensions of the Procrustes method for the optimal superimposition of landmarks. Systematic Biology 39(1): 40–59. https://doi.org/10.2307/2992207
  • Romano A, Amat F, Rivera X, Sotgiu G, Carranza S (2010) Evidence of tail autotomy in the European plethodontid Hydromantes (Atylodes) genei (Temmick and Schlegel, 1838) (Amphibia: Urodela: Plethodontidae). Acta Herpetologica 5(2): 199–206. https://digital.csic.es/handle/10261/43602
  • Rzedowski J (2006) Vegetación de México. 1ra. Edición digital, Comisión Nacional para el Conocimiento y Uso de la Biodiversidad, México, 504 pp.
  • Saunders SP, Ries L, Oberhauser KS, Thogmartin WE, Zipkin EF (2018) Local and cross-seasonal associations of climate and land use with abundance of monarch butterflies Danaus plexippus. Ecography 41(2): 278–290. https://doi.org/10.1111/ecog.02719
  • Schindelin J, Arganda-Carreras I, Frise E, Kaynig V, Longair M, Pietzsch T, Preibisch S, Rueden C, Saalfeld S, Schmid B, Tinevez JY (2012) Fiji: an open-source platform for biological-image analysis. Nature Methods 9(7): 676–682. https://doi.org/10.1038/nmeth.2019
  • Schmidt T, Arens P, Smulders MJ, Billeter R, Liira J, Augenstein I, Durka W (2009) Effects of landscape structure on genetic diversity of Geum urbanum L. populations in agricultural landscapes. Flora-Morphology, Distribution, Functional Ecology of Plants 204(7): 549–559. https://doi.org/10.1016/j.flora.2008.07.005
  • SEMARNAT (2010) Norma Oficial Mexicana NOM-059-SEMARNAT-2010, Protección ambiental-Especies nativas de México de flora y fauna silvestres-Categorías de riesgo y especificaciones para su inclusión, exclusión o cambio. Lista de especies en riesgo. Diario Oficial de la Federación, 10 diciembre 2010, México.
  • Sheets HD (2003) MakeFan Version 6, a tool for drawing alignment “fans” at equal angular spacing.
  • Signorell A, Aho K, Alfons A, Anderegg N, Aragon T, Arppe A (2016) DescTools: Tools for descriptive statistics. R package version 0.99, 18.
  • Soto VH, Yoshikawa K, Schörghofer N (2020) Climatic variation in the high mountains of central Mexico: Temperature and precipitation indices at Nevado de Toluca volcano. Atmósfera 33(4): 301–318. https://doi.org/10.20937/ATM.52768
  • Sunny A, Duarte-deJesus L, Aguilera-Hernández A, Ramírez-Corona F, Suárez-Atilano M, Percino-Daniel R, Manjarrez J, Monroy-Vilchis O, González-Fernández A (2019) Genetic diversity and demography of the critically endangered Roberts’ False Brook Salamander (Pseudoeurycea robertsi) in Central Mexico. Genetica 147(2): 149–164. https://doi.org/10.1007/s10709-019-00058-2
  • Sunny A, Monroy-Vilchis O, Reyna-Valencia C, Zarco-González MM (2014) Microhabitat types promote the genetic structure of a micro-endemic and critically endangered mole salamander (Ambystoma leorae) of Central Mexico. PLoS ONE 9: e103595. https://doi.org/10.1371/journal.pone.0103595
  • Vilaseca C, Méndez MA, Pinto CF, Benítez HA (2020) Assessment of Shape Variation Patterns in Triatoma infestans (Klug 1834) (Hemiptera: Reduviidae: Triatominae): A First Report in Populations from Bolivia. Insects 11: e274. https://doi.org/10.3390/insects11050274
  • Wake DB (1987) Adaptive radiation of salamanders in Middle American cloud forests. Annals of the Missouri Botanical Garden 242–264. https://doi.org/10.2307/2399397
  • Walls SC, Belanger SS, Blaustein AR (1993) Morphological variation in a larval salamander: dietary induction of plasticity in head shape. Oecologia 96(2): 162–168. https://doi.org/10.1007/BF00317728
  • Wickham H, Chang W, Henry L, Pedersen TL, Takahashi K, Wilke C, Woo K (2019) “ggplot2: Create Elegant Data Visualizations Using the Grammar of Graphics,” R Package Version 3.1.1. https://CRAN.R-project.org/package=ggplot2