Distribution, habitat suitability, conservation state and natural history of endangered salamander Bolitoglossa pandi

Ethics statement

Sex was not determined on living salamanders due to the high risk of injury to the animal. Fieldwork was done under the scientific research non-commercial purpose permit of collection of wild specimens of biological diversity issued by the National University of Colombia (Research Project 38615), and the Colombian National Environmental Licensing Authority (ANLA) by resolution No. 0255 of 14 March 2014. This study was conducted following the Colombian animal welfare law and the collection of wild specimens of the biological diversity acts (Ley 1774, 2016; Decreto 1376, 2013), as well as considering the Universal Declaration on Animal Welfare (UDAW) endorsed by Colombia in 2007.

Study area

We searched for salamanders at nine localities in four municipalities of Cundinamarca located on the western slope of the Cordillera Oriental of Colombia: Guaduas, Supatá, Venecia and Villeta (Fig. 1). Searches at each locality were carried out within an elevational gradient ranging from 1,638 to 2,315 m a.s.l (Table 1). The sampled area includes sub-Andean and Andean forests, as well as areas transformed by urban growth, agriculture, and cattle ranching. The sampling area is characterized by a bimodal climate [high dry season (from mid-December to mid-March); high rainy season (from mid-March through June); low dry season (from July to mid-September), low rainy season (from mid-September to mid-December)]. We monitored the environmental temperature (ET) and relative humidity (RH) at sampling sites using Ebro^® thermo-hygrometers (model EBI 20-TH1).

Sampling and data collection

We conducted surveys for salamanders using two approaches. First, during the rainy season in April–May 2013, we performed visual encounter surveys (VES; Crump & Scott, 1994) in three localities associated with cloud forests throughout an altitudinal gradient (1,648–2,002 m a.s.l) in the municipalities of Guaduas and Villeta (Table 1). Two researchers surveyed day and night for five days, for a total of 100 h of sampling effort.

Second, we employed the basic sampling protocol described by Mackenzie et al. (2003) for single-species occupancy modeling in six localities, three in the municipality of Supatá and three in Venecia. We randomly selected a total of 296 plots (5 m ×  5 m), which were located throughout an altitudinal gradient at each sampling locality (1,637–2,315 m a.s.l), grouping the following vegetation covers: Andean forest fragments, restored Andean riparian forest, pastures and roadsides (Table 1). During three sampling occasions (September–October 2017, March 2018, and July 2018), each plot was surveyed day and night for five consecutive days by ten researchers, resulting in a total of 2100 h of sampling effort. During each survey, the detection/non-detection of B. pandi specimens were recorded. When a salamander was present, we measured its perch height using a measuring tape (±0.1 cm). Once salamanders were caught, we recorded their weight with a Pesola^® dynamometer of 50 g (±0.1 g) and took photographs to measure their body size (SVL, snout-vent length (mm); TL, Tail length (mm)) using the software Image–J v. 1.52 (Bourne, 2010). All specimens of B. pandi were subsequently released near the plot where they were sighted.

Based on 18 specimens collected, euthanized using 2% lidocaine, and fixed in 10% formalin (Chen & Combs, 2001), we described the morphological variability of B. pandi. We made a small incision in the groin region to identify their sex and sexual maturity through macroscopic observation of the gonads. All the morphological terminology employed follows several contributions (Brame & Wake, 1962; Brame & Wake, 1963; Acosta-Galvis & Restrepo, 2001; Lynch, 2001; Acosta-Galvis & Hoyos, 2006; Acosta-Galvis & Gutiérrez-Lamus, 2012; Brcko, Hoogmoed & Neckel-Oliveira, 2013; Bingham et al., 2018). All specimens were deposited in the amphibian collection at Instituto de Investigación de Recursos Biológicos Alexander von Humboldt, (IAvH-Am), as well as in the amphibian collection at the Instituto de Ciencias Naturales de la Universidad Nacional de Colombia (ICN).

Habitat structure data collection

We used the Point Intercept Method described by Elzinga, Salzer & Willoughby (1998) to estimate the percent vegetation cover. We grouped plants into eight growth forms: graminoids, forbs, palm trees, mosses, lichens, vines, shrubs, and trees. We divided vegetation into layers: 0–0.1 m, 0.1–1 m, 1–1.5 m, 1.5–3 m, 3–5 m, and 5–12 m. We also considered ground characteristics such as leaf litter, bare soil and bare rocks. We estimated the vegetation cover of each of the plots where salamanders were surveyed, employing a set of 15 intercept points distributed in three parallel lines of five points separated by one meter of distance. At each point, we used a sampling bar of 1.5 m to register the contact of the growth forms of each vegetation layer below 1.5 m. This provided us with a 6.67% cover resolution by layer. We assessed the percentage cover of the vegetation layers above 1.5 m (mostly trees) using five intercept points: the corners of the quadrant and the central point of the third line. In this way, we reached a 20% cover resolution for upper layers.

Statistical analysis

We evaluated the association between habitat structure and the natural history traits of B. pandi by multiple correlation analysis, with P < 0.05 as the significance level. The following variables were considered: SVL (mm), weight (g), perch height (mm), leaf litter depth (mm), vegetation layers, vegetation life form, percent vegetation cover and elevation (meters above sea level). Using the habitat variables, we performed a principal component analysis (PCA) to explore which of these variables presented greater variability between plots with regards the detection/non-detection of B. pandi and, therefore, which of these could explain the observed differences between plots. The variable suitability for PCA was tested performing a Kaiser–Meyer–Olkin test (KMO > 0.5, P < 0.05). Afterward, a quadratic discriminant analysis was performed to determine which of the habitat variables had the greatest discrimination capacity between plots where B. pandi was detected/non-detected.

We assessed the variability in salamander abundance observed in the Supatá population through multiple regression analysis. First, we considered the salamander abundance as a dependent variable and the habitat structure variables recorded at each sampling plot as independent variables: leaf litter depth (mm), vegetation layers, vegetation growth forms, percent vegetation cover, elevation (meters above sea level), temperature (°C), and environment relative humidity. All variables were Ln–transformed prior to perform the statistical analysis.

Second, we evaluated assumptions of normality, autocorrelation, and homoscedasticity using Kolmogorov–Smirnov test, Durbin–Watson test and Breusch–Pagan test, respectively. Given that the p–value of the Durbin–Watson test can easily be less than 0.05 when sample size is very large, we used the Durbin–Watson statistic test (DW) as an autocorrelation criterion. According to Durbin & Watson (1950), a DW of less than 1 indicates a strong positive autocorrelation, a DW greater than 4 indicates a strong negative autocorrelation, values between 1 and 3 suggest a moderate autocorrelation, and a value close to 2 means that there is no autocorrelation.

Third, we tested for multicollinearity between the variables using the variance inflation factor (VIF) with a threshold of 10. Fourth, we selected the “best” regression model employing the Akaike Information Criterion (AIC; Akaike, 1973), considering that models with ΔAIC values of less than two are equally plausible (White & Burnham, 1999). Finally, we used the hierarchical partitioning method to evaluate the contribution of all the independent variables of the regression model (Chevan & Sutherland, 1991).

Activity pattern and population stage-structure

We only assessed the activity pattern and the population stage-structure of the B. pandi population at the Cuzcungos locality given the high abundance observed (Table 1). We estimated the activity pattern through kernel density curves for each sampling occasion, as well as combining all sampling occasions. We explored the salamander abundance variability during their activity period by performing a nested ANOVA. Hence, the activity period of B. pandi was divided into nine time-intervals (from H1 = 18:30–19:30, H2 = 19:31–20:30, H3 = …until H9 = 02:31–03:30), and each salamander sighting was allocated into its respective interval. The sampling occasion was used as the primary factor, and the time intervals as the secondary factor nested in the primary factor. We evaluated assumptions of normality and homogeneity of variances using a Shapiro–Wilk test and Levene test, respectively. Additionally, we analyzed the variability in observed body size over the nine time-intervals through a non-parametric ANOVA using a Kruskal-Wallis test (KW) as a measure of the central tendency of the samples (Sokal & Rohlf, 1981).

We used SVL as a descriptive variable of the population stage-structure of B. pandi. We compared the variability in population stage-structure among sampling occasions through a Wilcoxon test, with the null hypothesis being that population median stage-structure was the same across all sampling occasions. According to the categories proposed by Acosta-Galvis & Gutiérrez-Lamus (2012) and Del Río-García, Serrano-Cardozo & Ramírez-Pinilla (2014), as well as the development stage of the collected salamanders, the population was divided into stage classes as follows: neonates (≤23 mm), juveniles (24–30 mm), and adults (≥30 mm).

All statistical analyses were performed using the software Rwizard 4.3 (Guisande et al., 2014) and the following R packages: car (Fox & Weisberg, 2019), hier.part (Nally & Walsh, 2004), lawstat (Hui, Gel & Gastwirth, 2008), nortest (Gross & Ligges, 2015), overlap (Ridout & Linkie, 2009) stat (Bolar, 2019) and usdm (Naimi et al., 2014).

Geographic distribution

We found a total of 34 B. pandi individuals at three new localities, extending the geographical range of the species by 96.5 Km (airline) northwest from the type-locality, and 33.6 km (airline) west from the northernmost locality in the municipality of Supatá (Fig. 1). These new localities belong to the municipalities of Guaduas and Villeta in the department of Cundinamarca. All the salamanders in the new localities were found at night, within the understory of oak groves dominated by ferns. The salamanders were found in different vertical strata ranging from leaf litter, where they remained hidden, to shrubby substrates up to 2.5 m. Additionally, two salamanders were found in ecotonal areas associated with sugar cane crops and rangeland areas for livestock.

Associations of habitat structure and natural history traits

We observed a significant association between habitat structure and morphological traits. Snout-vent length, tail length, and mass were significantly associated with all the habitat structure variables assessed, but less so with the vegetation layers. SVL, TL, and mass showed a negative correlation with leaf litter and elevation (Table 2). In contrast, SVL and mass showed a positive correlation with perch height and percent vegetation cover (Table 2, Fig. 2).

Similarly, the habitat structure variables were significantly associated with the detection or non-detection of B. pandi throughout the sampling plots. The first two components of the PCA accounted for 62.4% of the variability observed. The habitat structure variability was clustered in two groups associated with the detection or non-detection of B. pandi (Fig. 3A). These groups were moderately overlapping in multivariate space, but they were differentiated by elevation, leaf litter depth, and percent vegetation cover. The presence of B. pandi was positively correlated with highly structured plots and with deep leaf litter (Table 3). The quadratic discriminant analysis confirms that the detection or non-detection of B. pandi depends on habitat variables. The cross-validation percentage was 92.6%, indicating that the quadrants in which B. pandi was detected can be clearly distinguished by habitat variables such as vegetation layers, leaf litter depth and mean environmental temperature (Fig. 3B).

Results of the multiple regression analysis showed that leaf litter depth, mean environmental temperature, percent vegetation cover, and elevation were the variables that contributed the most to explaining the observed variability in the abundance of B. pandi (r² = 60.3%, P < 0.001). These habitat variables also composed the best-fitted regression model (Table 4). However, this model showed a moderate autocorrelation (DW = 1.10), which means that the variance explained by the habitat variables may be close to 60.3%.

Activity pattern

Bolitoglossa pandi is completely nocturnal; its activity period extends throughout the night, from 18:30 h until 05:00 h. The environmental temperature recorded during this activity period ranged from 12.6–25.6 °C ($\bar{X}=15.3$), and relative humidity ranged from 73.8–98.8 ($\bar{X}=80.3$). The activity peaks were all subject to environmental influences showing a significant association with the local weather conditions. The observed salamander abundance was strongly and positively correlated with the environmental temperature (R_ET = 0.251, P < 0.001), whereas it was moderately and negatively correlated with relative humidity (R_RH = − 0.174; P = 0.017). Across sampling occasions, we observed two consistent activity peaks, the first from 20:30 h to 21:30 h, and the second from 23:30 h to 00:30 h (Fig. 4A). Results of the nested ANOVA indicated that salamander abundance between these peaks was significantly different (F₉₋₂₅₉ = 4.57, P < 0.001; Fig. 5A, Table 5), being higher in the second activity peak during the second and third sampling occasions, with the opposite pattern during the first sampling occasion (Figs. 4B–4D).

Likewise, the observed body length of the salamanders showed significant differences between the two activity peaks (K_8df = 99.70, P < 0.001, Fig. 5B), suggesting niche time-partitioning between body size classes. Most of the salamanders observed during the first activity peak had SVL >30 mm (adults), whereas all the salamanders observed at the second activity peak had SVL <27 mm (juveniles and neonates). However, in contrast to the differences in observed abundance between activity peaks, the variation in salamander SVL was between peaks was consistent across the three sampling occasions.

Population stage-structure

A total of 1391 individuals of B. pandi were observed throughout the study in the Supatá population, exhibiting a population density ranging from 0.04–1.44 individuals/m². The population stage-structure showed significant differences between sampling occasions, suggesting seasonal variability between development stages (Table 6). Regardless of the variability observed in the population stage-structure of B. pandi, the population is mostly dominated by juveniles and neonates which represent between 34–64% of individuals (Fig. 6). Neonates were observed across all sampling occasions, indicating constant recruitment. During the second sampling occasion, a bias in the median SVL towards smaller individuals was observed, suggesting that the major recruitment peak occurs during March when the high rainy season start (Fig. 6B). However, no samples were available between October to mid-December when the low rainy season occurs, thus, a second recruitment peak could be possible observed during this period. On the contrary, adults were conspicuous throughout all sampling occasions, especially so in the first sampling occasion.

Morphological variability and comparisons to other species

Bolitoglossa pandi is a small species, with SVL = 13.4−51 mm (N = 1034), 34.4−48 mm adult males (n = 7), and SVL 36.6 −51 mm adult females (N = 8). Extensively webbed hands and feet with third toes and triangular digits; ventral surfaces of digit tips without terminal flattened tubercles; snout short and rounded in the lateral profile; head length 4.9–10.1 mm; head width 5.1–7.9 (N = 17); snout rounded in dorsal view, irregular white spots, cream-colored nasolabial grooves, and edges of the lips irregularly dark brown with irregular light spotting (Figs. 7A–7B); protruding eyes on dorsal view, brown iris with black reticules (Fig. 7B); well-defined post-cephalic constriction; ventral surfaces (preserved) brown or dark grey with numerous tiny cream guanophores (Figs. 7C–7D); inverted bracket shaped scapular spots (Fig. 7J); males have white testes.

Bolitoglossa pandi can be distinguished from its Colombian congeners by the extensive webbing of its hands and feet (versus webbing of hands and feet reduced in B. adspersa, B. hiemalis, B. hypacra, B. palmata, B. ramosi, B. savagei, B. tamaense, B. tatamae, B. walkeri and B. vallecula). Furthermore, it can be distinguished from species with extensive webbing (such as Bolitoglossa lozanoi and B. nicefori) by having more protruding eyes, and a longer and triangular third digit. Bolitoglossa pandi also differs from B. biseriata and B. silverstonei by having a dark brown or dark grey ventral surface with irregular white dots (versus cream ventral surface with brown suffusions and dots in B. biseriata and B. silverstonei), brown-reddish iris with black reticles (versus golden with brown dots in B. biseriata and B. silverstonei), and white testes in adult males (versus black mesorchium in testes of adult males of B. biseriata). It also differs from B. medemi by the absence of digital depressions in the tips on the digits and toes (present in B. medemi). B. pandi also presents an upper lip with irregular light spotting (versus uniform in B. guaneae). B. pandi can be differentiated from B. altamazonica and B. leandrae by having extensive interdigital webbing with a longer and triangular third digit (complete webbing and tips rounded in B. altamazonica and B. leandrae). B. pandi is morphologically very similar to B. phalarosoma, but it differs by having a dark brown or dark grey ventral surface with diffuse white dots and some white blotches (versus usually light brown ventral surfaces, with some irregular cream spots in B. phalarosoma), plantar and palmar regions in ventral dark brown or dark grey view (versus cream in B. phalarosoma). B. pandi also differs from B. capitana by having a smaller adult size with a longer and triangular third digit (versus rounded third digit in B. capitana) and a shorter head.

Color variability

Bolitoglossa pandi exhibits a wide variation in color, ranging from a uniform dorsal color pattern of different shades of light to dark brown (Fig. 7E); one reddish-brown pattern that can have diffuse grey or scarcely distinguishable dark blotches that extend to the dorsolateral region (Fig. 7F), or an ochre pattern, with some diffuse irregular dark brown, yellow or cream spots (Figs. 7G–7I). Some specimens exhibit a very diffuse band in the paravertebral region covering almost the entire dorsal surface or including a dark brown inverted triangle shape in the interorbital region (Figs. 7J–7K). The caudal region is highly variable, ranging from uniform reddish-brown or completely ochre segments, with very small scattered white dots (Figs. 7D–7F) to irregular cream, yellow or orange patches (Fig.  7I), and very small or irregular longitudinal black spots (Figs. 7I–7K). The distal end of the tail becomes uniform light brown and has cream blotches towards the proximal region in some individuals (Fig. 7H). The cephalic region is dark brown with some white spots with irregular dots up to the supralabial region in lateral view (Fig. 7B); a loreal region with ochre patches and light brown iris with black reticles; nasolabial projections are cream; the ventrolateral surface is dark brown. Ventral surfaces are dark brown with some cream and white blotches with scattered dots; the mental and gular surfaces are uniformly dark brown or dark grey with some irregular white circular spots bordering the maxillary region. In adult males the mental gland is light brown; the palmar and plantar surfaces are always dark brown or dark grey.