When no color pattern is available: Application of double observer methods to estimate population size of the Alpine salamander

ABSTRACT

Monitoring wild populations is an essential tool to assess the conservation status and the ecological requirements of a species. Capture–mark–recapture (CMR), based on individual recognition, is the most commonly used and most effective technique. However, in cases of species with no individual color pattern, tracing the encounter history of individuals without invasive marking methods is impossible. In this study we aimed to (1) estimate population abundance and density using a less effort-intensive and nonstressful technique, (2) test a long-term monitoring protocol, and (3) assess the fine-scale ecological requirements of a black-colored amphibian, Salamandra atra, in the Italian Alps. For three populations we applied an N-mixture model on data collected using a dependent double-observer approach. To understand ecological requirements, we assessed the relative importance of a set of environmental and topographical variables. The double-observer approach was a cost-effective technique that provided reliable demographic estimates of population density. Our results suggest that the most important fine-scale ecological variables positively associated with salamander abundance were canopy cover and terrain ruggedness, which are strictly related to shelter availability and soil moisture.

KEYWORDS:

• Abundance estimation
• Alps
• amphibians
• double observer
• salamanders

Introduction

There are several reasons for monitoring wild animal populations. Among them, to (1) determine the conservation status of species, (2) understand their ecological needs in order to plan appropriate conservation strategies, (3) assess the effectiveness of management and biodiversity actions, and (4) know the effects of anthropogenic pressures in a given territory (Yoccoz, Nichols, and Boulinier 2001; Witmer 2005; Reynolds, Thompson, and Russell 2011; Shrestha and Lapeyre 2018). However, obtaining accurate estimates of population abundance and density is time- and resource-consuming, mainly because of the imperfect detection of wild animals in their habitat (Seber 1982; Williams, Nichols, and Conroy 2002). Imperfect detection originates from our inability to capture or sight all the individuals of a wild population and from the fact that the available ones for sampling are just a fraction of the entire population (Nichols and Conroy 1996; Borchers et al. 2002; Schmidt 2003).

A rich array of analytical methods, known as capture–mark–recapture (CMR) protocols, have been developed to analyze data from marked populations, and these are considered the most informative techniques to estimate population abundance and other demographic parameters (Seber 1982; Williams, Nichols, and Conroy 2002; Barker et al. 2017). CMR models are based on the encounter history of captured individuals in successive recapture occasions, allowing the estimation of the population size (Seber 1982; Schmidt 2003). Traditionally, individual recognition is accomplished by capturing animals and placing unique marks on them. Usually these methods are not invasive or minimally invasive, but in some cases unexpected and negative effects on marked individuals have been observed (e.g., Saraux et al. 2011). Accordingly, these possible animal welfare consequences (McMahon, van den Hoff, and Burton 2005) and the effort associated to capture and physically mark individuals have pressed researchers to develop alternative methods (Zemanova 2021). Therefore, in recent years, digital photography and image processing software have gained popularity in providing photographic methods (i.e., photographic mark–recapture) to recognize individuals by their unique color marks or patterns (Bolger et al. 2012; Zemanova 2021). The fundamental prerequisite for the photographic marking application is that individuals must exhibit color patterns on some region of their coat, skin, or eyes sufficiently distinctive to discriminate among them and that this pattern must be stable over the duration of the study period (Bolger et al. 2012). However, several methods of photo identification often require the capture and handling of individuals, in particular when the focal species is small, cryptic, and possesses its characteristic pattern in hidden parts such as the gular, ventral, or abdominal regions.

In situations where invasive marking methods (e.g., visible implant elastomers, passive integrated transponders, tattoos, etc.) are not desirable and individuals do not have discriminating color patterns—for instance, in completely monochromatic black species—alternative methods such as distance sampling (DS; e.g., Buckland et al. 2004) and an adaptation of the CMR framework known as “multiple observers” (MO) method (Southwell 1996, table 16) can be applied. In DS animals are counted along a line transect or a point, and the distance from the center of the transect is recorded for each observation. Population density from line transect surveys may be therefore estimated by modeling the way detectability decreases with distance from the observer. Unfortunately, DS has been judged only moderately effective for amphibians and reptiles, because these animals may be undetected even on the transect line (i.e., one of the fundamental assumptions of DS), leading to assumption violation and imprecise density estimates (see Dodd [2010], [2016] and references therein; see also Marques 2009). In the MO approach, animals are sighted and recorded by two or more observers in succession, with a protocol that allows partitioning of the sightings into those scored by each observer. MO has been mainly applied to animals readily observable from long distances, such as large terrestrial or marine mammals (Cook and Jacobson 1979; Langtimm et al. 2011; Bröker et al. 2019) and waterfowl (Koneff et al. 2008; Vrtiska and Powell 2011), and also to animal traces and signs (amphibian egg masses; Grant et al. 2005) or animals calls (frogs and toads: Shirose et al. 1997; birds: Forcey et al. 2006). Recently, multinomial N-mixture models, which are a generalization of the binomial N-mixture model (Royle 2004), have been developed to accommodate sampling frameworks yielding multinomial outcomes, such as the MO protocol under a metapopulation design (Chandler, Royle, and King 2011; Kéry 2018; Costa, Romano, and Salvidio 2020). These multinomial N-mixture models have been successfully applied mainly on large mammals with high detection probability (e.g., manatees: Langtimm et al. 2011; narwhals: Bröker et al. 2019) but also on terrestrial amphibians (Costa, Romano, and Salvidio 2020). Amphibians and reptiles include several mimetic, nocturnal, small-sized, and cryptic species that for these reasons are typically characterized by a low detection probability (e.g., McDiarmid et al. 2012; Griffiths et al. 2015).

In this study, we used the multinomial N-mixture model on data collected using a dependent double-observer approach to estimate population abundance and density of a completely black-colored and fully terrestrial salamander (Figure 1). To date, population abundance estimates of the alpine salamander (Salamandra atra Laurenti, 1768) have been obtained using invasive or highly stressful methods, such as toe clipping (Luiselli et al. 2001) and transponders after anesthetization (Helfer, Broquet, and Fumagalli 2012), and a moderately stressful approach; that is, removal methods (Romano et al. 2018a). In addition, we evaluated whether this approach can be used as an expeditious method in a long-term monitoring program to detect possible demographic trends of alpine salamander populations, which is among the aims of the monitoring program of the Natura 2000 Network.

Figure 1. The alpine salamander, Salamandra atra, a fully terrestrial and completely black amphibian, endemic of the Central and Eastern Alps and Dinaric Alps.

Materials and methods

Study species

The alpine salamander, Salamandra atra, is a fully terrestrial and viviparous salamander endemic to the Central and Eastern Alps and the Dinaric Alps where some isolated populations have been recorded (Sillero et al. 2014). The alpine salamander is a polytypic species: two Italian endemics possess a very narrow distribution range and yellowish spots on the body (S. atra aurorae and S. atra pasubiensis). Two other subspecies, S. a. atra and S. a. prenjensis, both with wider distributions, are characterized by completely black coloration. This study was performed on the most diffused subspecies, the nominal one, for which the typical habitats are mixed coniferous and deciduous forests, alpine meadows, and rocky tundra-like areas, mainly on limestone substrates. Annual activity pattern is concentrated in the warmest months (April–October), whereas during the rest of the year the salamanders are inactive (Klewen 1988). Moreover, these salamanders are active during rainy conditions and in particular at night.

Study area, study framework, and sampling design

The study area is located in the Italian Dolomites, in the Paneveggio–Pale di San Martino Natural Park (province of Trento, Northern Italy), and is characterized by open habitats (pastures, other grasslands, and rocky areas) mixed with coniferous woodlands dominated by European larch (Larix decidua) and Norway spruce (Picea abies). Three sites were chosen, in similar habitat conditions, to start population monitoring, based on a series of inspections performed in 2019, which also took into account the accessibility and safety of operators at night. Each study site was a transect 300 m long. The first, OTA, is at about 1,900 m.a.s.l., near the locality Malga Venegiota (cadastral municipality of Tonadico; 46°18′48″ N, 11°48′53″ E); the second, GIA, is at about 1,850 m.a.s.l., near the locality Malga Venegia (cadastral municipality of Tonadico; 46°19′21″ N, 11°47′46″ E), and the last one, FOS, is at about 1,950 m.a.s.l., near the locality Malga Fosse (cadastral municipality of Siror, 46°17′20″ N, 11°47′57″ E). The linear distance between transects ranged from 2.3 km (OTA-GIA) to 3.9 km (GIA-FOS). The transects were on mountain slopes with different exposures: OTA and FOS had a southwest aspect and GIA had a southeast aspect. Each transect was divided in thirty subtransects of 10 m each (measured by means of a range finder), which were individually flagged a few days before salamander sampling and constituted sampling units of the dependent double-observer surveys. Several environmental variables describing salamander microhabitat were also measured for each subtransect. Three variables were detected on the field: average width of the subtransect (width, measured with a metric wheel tape), canopy cover (canopy) divided in three categories (high = cover >70%, medium = cover >20% and <70%; low = cover <20%), and presence/absence of shelters (shelter; e.g., stones, deadwood). Two additional variables were derived from a digital elevation model of the study area (1-m mesh size; Lidar 2009–PAT) and calculated within SAGA-GIS software v.7.9 (Conrad et al. 2015): the ruggedness of the subtransect (ruggedness; i.e., the relative measure of elevation difference between adjacent elevation cells), calculated as the Terrain Ruggedness Index (Riley, DeGloria, and Elliot 1999), and the Topographic Position Index (TPI), which expresses the topographic position of each cell within the landscape (Guisan, Weiss, and Weiss 1999). Although Lauber (2004) found no main effect of air humidity on salamander counts, our pilot studies (Roner et al. in press) showed significant differences in the number of active salamanders in different meteorological conditions and times of the day, confirming Günther and Grossenbacher (1996) and Geiger (2006). Consequently, salamanders were counted when they were more active; that is, during a rainy night at the end of August 2020, between 9:00 p.m. and 02:00 a.m. The dependent double-observer technique was used, in which observer 1 points to and counts out all salamanders in each subtransect to observer 2 (Cook and Jacobson 1979; Nichols et al. 2000). Observer 2 records what observer 1 reports but also records in a separate column any additional salamanders detected by him and missed by observer 1, avoiding any possible influence. At each subtransect the two observers exchanged roles. Each transect was covered in a time ranging from 30 to 45 minutes.

Data analysis

We analyzed dependent double-observer data using multinomial N-mixture models, employing each subtransect as a site (n = 90). For the double-dependent observer protocol, multinomial N-mixture models include two parameters: one for the latent variable (i.e., mean abundance at sampling sites; lambda, λ) and another one for the detection process (i.e., detection probability, p; Royle and Dorazio 2006; Chandler, Royle, and King 2011; Costa, Romano, and Salvidio 2020). Both abundance and detection parameters can be conveniently modeled as a function of environmental or sampling covariates through a log or logit link, respectively (Kéry and Royle 2016). Prior to model building, all continuous predictors were scaled by means of function scale() within the R environment (mean = 0; SD = 1). We built a global model; that is, the more complex model on which other models are nested, with a Poisson error distribution for λ and including all sub-transect-specific covariates on the abundance side of the formula. We modeled p with a fixed observer-dependent effect, which we retained in all models, and considered p to be also affected by width. Model building occurred within the R package “unmarked” using function multinomPois(), set for data type “depDouble” (Fiske and Chandler 2011). We conducted model selection within the R package “AICcmodavg” (Mazerolle 2020) by simplifying the global model using stepwise backwards elimination, dropping one term at a time, beginning from the one with the highest p value and removing terms until Akaike’s informative criterion corrected for small sample size (AICc; using number of sites as sample size for a more conservative approach; Burnham and Anderson 2002; MacKenzie et al. 2017) no longer decreased (Kéry and Royle 2020). From the best model, we derived abundance estimates (N-hat)for each population from the posterior distribution of the latent abundance using function ranef() in the R package “unmarked.” From the posterior distribution of the latent abundance, we also derived density estimates (D-hat; individuals/hectares) for each population by calculating the total surveyed area for each population using the subtransect specific width. Confidence intervals for N-hat and D-hat were obtained by parametric bootstrap, using function parboot() in the R package “unmarked” (Fiske and Chandler 2011; Kéry and Royle 2016).

Results

During the surveys we counted a total of 109 salamanders (see Supplementary Materials). After backwards model selection, the resulting best model accounted for the effect of width and observer on p and of canopy and ruggedness on λ (all models, also including the global model, are presented in Table 1). Detection probability was negatively affected by subtransect width (Table 2, Figure 2). Estimates of detection probability were high and not significantly different between observers (observer 1; i.e., the observer pointing to detected salamanders, p = .89, 95 percent confidence interval [CI] [0.76–0.96]; observer 2; i.e., the observer recording salamanders pointed by observer 1 and also those detected by himself but missed by observer 1, p = .85, 95 percent CI [0.24–0.99]), and the combined detection probability—that is, the probability of an individual being detected by at least one observer (see Nichols et al. 2000)—was p = .98 (95 percent CI [0.82–0.99]). Mean expected abundance (λ = 1.04, 95 percent CI [0.78–1.39]) was positively affected by medium and high canopy cover, negatively affected by low canopy cover (Table 2; Figure 2), and positively affected by terrain ruggedness (Table 2, Figure 2). Population-level abundance estimates, derived from the posterior distribution of the latent abundance, greatly differed among populations. Abundance and density estimates in GIA (N-hat = 26, 95 percent CI [24–36]; D-hat = 231, 95 percent CI [134–376]) and FOS (N-hat = 22, 95 percent CI [22–24]; D-hat = 299, 95 percent CI [190–480]) were similar, whereas a higher abundance was estimated in OTA (N-hat = 66, 95 percent CI [63–87]; D-hat = 570, 95 percent CI [435–883]).

Table 1. List of candidate models, ranked by AICc

Model	Parameters	AICc	ΔAICc
p(Observer+Width) λ(Canopy+Ruggedness)	7	314.41	—
p(Observer+Width) λ(Canopy+Ruggedness+TPI)	8	314.67	0.26
p(Observer+Width) λ(Canopy+Ruggedness+TPI+Shelter)	9	316.90	2.49
p(Observer+Width) λ(Canopy+Ruggedness+TPI+Width)	10	319.40	4.99

Table 2. Parameter estimates and 95 percent confidence intervals from the best model. Detection parameters are on the logit scale. Abundance parameters are on the log scale

Parameter		Estimate	95% Confidence interval
Detection (p)	Intercept (observer 1)	2.13	1.16–3.1
	β Observer 2	−0.38	−2.63 to 1.86
	β Width	−1.00	−1.91 to −0.09
Abundance (λ)	Intercept (canopy: high)	0.04	−0.24 to 0.33
	β Canopy: low	−0.61	−1.35 to 0.12
	β Canopy: medium	0.40	0.02–0.80
	β Ruggedness	0.52	0.36–0.67

Figure 2. Plots representing the relationship between expected values of abundance and detection and covariates included in the best model. Transect width is expressed in meters. Three plots for abundance represent different levels of canopy. Blue and red lines represent the mean for abundance and detection, respectively. Grey lines represent 95 percent confidence intervals.

Discussion

The MO method is a cost-effective approach to estimate alpine salamander abundance. Bailey et al. (2004) argued that, in the short term, reliable methods (i.e., robust design and removal methods) usually estimate surface population of salamanders, because only in long-term studies the temporary emigration may be taken into account.Nonexposure to sampling can have significant consequences for abundance estimation. Consequently, the results we obtained can be comparable to other short-term methods (e.g., Romano et al. 2018a), rather than to long-term studies. We estimated the abundance of salamanders active on the surface during optimal weather conditions; therefore, we are confident that the population active on the ground surface is a proxy of the true population. In the weather and night conditions in which we performed the surveys, salamanders exhibited a very high detection probability, which led to an estimate close to the observed number of individuals. In the site OTA, salamander density was also estimated using the removal method three years before this study (472 ind./ha; 95 percent CI [310–633], estimated in September 2017; see Romano et al. 2018a) and did not significantly differ from the results we obtained in the present study (D-hat = 570, 95 percent CI [435–883]). Considering that site OTA has not undergone any considerable changes during these three years, the concordance of the two methods in the density estimate highlights their mutual robustness and reliability. Population density was remarkably higher in site OTA than in the other two sites (GIA and FOS).

Available literature shows that population density of the nominal subspecies of the alpine salamander is highly variable among sites, spanning from more than 3,000 ind./ha (Helfer, Broquet, and Fumagalli 2012) to 120 ind./ha (Klewen 1986), also including intermediate values (Klewen 1986, 1988; Helfer, Broquet, and Fumagalli 2012). However, a higher value was obtained extrapolating on a hectare from a small area (625 m²), whereas lower values were extrapolated from larger areas (3,000 m²). Considering that, as for many other small vertebrate species, the spatial distribution of Salamandra atra may be significantly aggregated also within small areas (about 1,300 m²; S. a. aurorae; Bonato and Fracasso 2003), extrapolation of population density from very small areas to hectares should be regarded with caution because it could greatly under- or overestimate the actual densities at larger scale. To clarify this question—that is, whether the densities along the transects are comparable to the surrounding areas (as suggested by our preliminary studies in progress)—further research is needed.

The MO method applied here is confirmed as a useful tool to perform abundance and trend monitoring of this elusive amphibian, for which marking techniques are not a viable option. Furthermore, this technique provides an important reduction of sampling effort compared to CMR and removal methods. Finally, the MO method eliminates the need for capture, marking, and handling of individuals, preventing any possible stress or negative effects on animals.

The included microhabitat variables in our modeling framework allowed us to conveniently evaluate the fine-scale ecological requirements of the alpine salamander. Local abundance was positively affected both by closed canopy cover and by terrain ruggedness. Higher canopy cover could favor salamander occurrence and abundance at two different levels. First, the positive effect of the canopy cover can be directly related to phenomena associated with protection from solar radiation and desiccation. In addition, considering that alpine salamanders are usually active at night and during rainy weather, when they retreat under shelters or in shallow subterranean habitats, they are exposed to climatic variations, linked to the surface environment (Culver and Pipan 2014). Therefore, a higher canopy cover may have a buffering effect on climatic variationsin subterranean shelters. Second, canopy cover is intuitively correlated to the number of trees and their size, which may in turn be directly linked to shelter availability because the basal part of the tree trunks, with their cavities and buttresses, is used by salamanders as shelter (Basile et al., 2017; Piraccini et al. 2017). However, our results show that salamanders were more abundant along subtransects with intermediate canopy cover (Table 2). A preference for areas with intermediate canopy cover can be explained by the effect of canopy on soil moisture, resulting from both the outcome of gross precipitation and canopy interception, and also from the shading effect of the canopy on soil drying. Canopy interception loss ranges between 9 and 48 percent of gross precipitation and tends to be higher for coniferous than for broad-leaf forests (Gerrits and Savenije 2011). Intuitively, soil drying is slower if the soil is shaded and if it is protected from the wind by trees. Subtransects characterized by high conifer canopy cover require extremely abundant rainfall to humidify the soil, whereas subtransects without or with low canopy cover are more exposed to direct solar radiation and wind action, resulting in faster soil drying. Subtransects with intermediate canopy cover are probably the best compromise for a greater maintenance of soil moisture, which is essential for salamanders.

The second significant environmental variable (i.e., the terrain ruggedness) has to be regarded from this perspective. Terrain heterogeneity is indeed an important variable for several reasons and for many species: it provides different information considering different spatial scales (e.g., Nellemann and Reynolds 1997; Riley, DeGloria, and Elliot 1999; Sappington, Longshore, and Thompson 2007; Einoder et al. 2018). We used the Terrain Ruggedness Index at a spatial scale (i.e., 1-m meshes) suitable for describing the relevant elements for salamanders considering their size, behavior, ecology, and vagility. At this scale, a higher index corresponds to a greater abundance of rocky geomorphological elements on the ground that can provide shelters for alpine salamanders. Terrestrial salamanders are indeed highly constrained to a narrow range of environmental conditions due to their physiological requirements (Feder 1983; Peterman and Semlitsch 2014), and they are able to behaviorally buffer nonoptimal environmental conditions, in particular temperature and moisture (e.g., Farallo, Wier, and Miles 2018). One strategy is to be nocturnal, as the alpine salamander is, but this does not seem to be sufficient, and terrestrial salamanders regulate these environmental parameters through microhabitat selection (Farallo, Wier, and Miles 2018). However, considering their limited spatial movements (Bonato and Fracasso 2003; see table 8.1 in Vitt and Caldwell 2013) salamanders are strictly dependent to shelter availability in the surrounding (Popescu and Hunter 2011; Peterman and Semlitsch 2013). In the subspecies S. a. aurorae, the uneven availability of suitable rocky shelters was considered the main factor influencing the heterogeneous distribution of the population studied by Bonato and Fracasso (2003), and the availability of moist shelters was identified as a driving variable in determining species occurrence in forest patches (Romano et al. 2018b). Shelter availability, associated to canopy cover and high terrain ruggedness, seems to play a key role in determining differences in the alpine salamander abundance at a small spatial scale.

Conclusions

In this study we tested a method typically used for large species or signs of presence, and therefore with high detection probability and easily recognizable, on a small, elusive, and not easily markable amphibian with a uniformly black coloration. The dependent double-observer method was reliable and in agreement with the abundance estimate obtained from the removal method in a previous study. Furthermore, the analysis of microhabitat variables provided important information about habitat selection and ecological requirements; these factors are of paramount importance for proper species management and conservation. Ongoing studies on other amphibians and reptiles will allow us to test the validity of the double-observer approach in providing a robust and noninvasive method for those species, such as the alpine salamander, that are difficult to mark or for species whose handling is not recommended. Finally, annual monitoring using the double-observer protocol will allow us to detect demographic trends, leading to a better assessment of the conservation of salamanders in alpine habitats.

Acknowledgments

The Paneveggio–Pale di San Martino Natural Park provided logistic facilities close to the sampling site. In particular, we are grateful to the Technical Coordinator Piergiovanni Partel, the Director Vittorio Ducoli, and the park guard Gilberto Volcan. Thanks to Rachele Gobbi for the contributions to field sampling.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This study is a part of the Natura 2000 Network Habitat Directive monitoring and was partially supported by Servizio Aree Protette e Sviluppo Sostenibile PAT (Piano di monitoraggio fauna vertebrata Direttiva Habitat; Progetto LIFE11/ NAT/IT/000187 “TEN”–Trentino Ecological Network) and by the Paneveggio–Pale di San Martino Natural Park (Fondi Piano di Sviluppo Rurale PAT).