Environmental suitability models predict population density, performance and body condition for microendemic salamanders

Species can show strong variation of local abundance across their ranges. Recent analyses suggested that variation in abundance can be related to environmental suitability, as the highest abundances are often observed in populations living in the most suitable areas. However, there is limited information on the mechanisms through which variation in environmental suitability determines abundance. We analysed populations of the microendemic salamander Hydromantes flavus, and tested several hypotheses on potential relationships linking environmental suitability to population parameters. For multiple populations across the whole species range, we assessed suitability using species distribution models, and measured density, activity level, food intake and body condition index. In high-suitability sites, the density of salamanders was up to 30-times higher than in the least suitable ones. Variation in activity levels and population performance can explain such variation of abundance. In high-suitability sites, salamanders were active close to the surface, and showed a low frequency of empty stomachs. Furthermore, when taking into account seasonal variation, body condition was better in the most suitable sites. Our results show that the strong relationship between environmental suitability and population abundance can be mediated by the variation of parameters strongly linked to individual performance and fitness.

. Potential relationships between climatic suitability, population parameters and species abundance. In blue the proposed pathways, in orange pathways for which no data are available in this study. Filled arrows represent hypotheses tested in the present study.

Results
Environmental suitability. Overall, we obtained Hydromantes flavus records in 25 grid cells, which represent all the known localities of the species (Fig. 2). An ensemble SDM showed an excellent performance, with an overall true skill statistic = 0.879 (sensitivity: 100%; specificity: 87.7%; standard deviation of the true skill statistic across the models: 0.159). The SDM suggested that environmental suitability (ES) increased in the highest elevation areas of the Monte Albo (North Eastern Sardinia), where mean temperature was low and annual precipitation was high (Fig. 2, Supplementary Figs S1, S2).
Model calibration was good. The coefficient of determination between the proportion of positive cases and the identity line was high (R 2 = 0.56) (Fig. S3). This suggested that the model is useful both for discriminating and for ranking, even though the discrimination capacity of the model was higher than in a perfectly calibrated model 32 . The point-biserial correlation calculated on the basis of the re-calibrated POC-plot was high (COR = 0.52), indicating very good calibration 33 .

Distribution of salamanders.
We measured the distance from the surface (depth) of 178 salamanders from ten populations ( Table 1). The average depth was highly variable among populations (Table 1), ranging from 8.5 to ~150 m. The average depth was significantly higher in sites with low ES; results were identical if we analyse data from all caves (r = −0.74, N = 10, P = 0.014; Fig. 3a), or if we only considered caves for which at least 5 depth measurements were available (r = −0.73, N = 8, P = 0.038).
Stomach contents. We collected stomach contents from 212 salamanders in three populations (Table 1).
One population inhabited an area with limited suitability (ES = 0.431), while the other two were found in areas highly suitable (ES > 0.53). Thirty-seven stomach contents showed only unidentifiable material and were therefore discarded from the analyses (13 for Site 1, 21 for Site 5 and 3 for Site 7). Twenty-three salamanders had empty stomach, while in 152 we recognized at least one prey item. The frequency of salamanders with empty stomach strongly differed among populations (χ 2 = 14.52, df = 2, P < 0.001); empty stomachs were more frequent in salamanders living in less suitable sites (Fig. 3b). The frequency of empty stomachs was unrelated to the survey period (χ 2 = 1.39, df = 1, P = 0.239). Bayesian credible intervals (CIs) confirmed that the salamanders of the two sites with the highest ES had a similar low frequency of empty stomachs, as their 95% CIs showed wide overlap, while the pattern was clearly different in the least suitable site (Fig. 3b).
Abundance and density of populations. Population abundance was estimated using a N-mixture model with Poisson error, as it showed lower AIC than zero-inflated models (Poisson model: AIC = 555.74, zero-inflated Poisson: AIC = 557.74). N-mixture models indicated a high detection probability of individuals (mean ± SE: 0.504 ± 0.029). The estimated number of salamanders was highly variable among sites, ranging from five to 103  (Table 1); the average of density estimates across all sites was 0.06 ± 0.03 individuals/m 2 . Salamanders density strongly differed among sites, being significantly higher in areas showing the highest environmental suitability (r = 0.893, N = 7, P = 0.007; Table 1 and Fig. 3c).

Discussion
Analysing multiple population parameters is essential to unravel the complex processes linking environmental suitability and individual fitness, and to understand the factors determining spatial variation of species abundance ( Fig. 1). In this study, intensive sampling on multiple populations provided data on several key features of populations and allowed to test multiple hypotheses on processes that can affect the relationship between the abundance or density of a species, and environmental suitability. Even though we did not consider all the potential factors and population features (e.g., no data on breeding success or survival are available, Fig. 1), our study provides a first insight on how the interplay of multiple processes can determine the variability of abundance that can be observed across a species' range.
Modelling the ES of Hydromantes flavus allowed us to identify strong heterogeneity of suitability within the Monte Albo, with suitability variation at least in part related to altitude (Figs 2 and S1, S2). At the mountain base, the microclimate is drier and warmer compared to the top, and the wet conditions at the top of the mountain probably increase suitability for Hydromantes (Fig. 2). Hydromantes are lungless salamanders that require high environmental humidity for breathing, have a narrow thermal niche, and their activity at the surface is only possible during wet, cold periods 26 . In higher ES areas, local climate shows a prolonged suitability for Hydromantes, a condition that likely reduces their inactivity period 27 . Indeed, when external climate is cool and moist, salamanders can exit from their refuges to prey in environments with high prey abundance 34 . As a consequence, SDM identified a suitability gradient from the lowest to the highest altitudes of the Monte Albo.
We expected a limited impact of interspecific interactions on Hydromantes flavus abundance 24,26 . Nevertheless, the relationship between suitability and population density was not strictly linear, as low densities were also observed in some populations living in areas with high suitability (Fig. 3c). Triangular relationships between ES and species abundance occur when biological interactions hamper populations to reach their maximum achievable abundance 5,11 . For instance, some areas of the Monte Albo are strongly affected by urbanization and mining 35 , thus human activities can determine variation of parameters not measured in this study such as vegetation features or invertebrate distribution, which can in turn affect salamander populations.
The activity hypothesis predicts that a species is more active when environmental conditions are suitable 13 . Hydromantes are epigean terrestrial salamanders that can spend long periods in underground environments to escape unsuitable environmental conditions, particularly during dry and hot seasons 26 . However, underground environments are not homogeneous, and are strongly influenced by external environmental conditions 36 . The activity of salamanders, and their exploitation of microhabitats, is strongly affected by a trade-off between their physiological constraints and the need of food. On the one hand, microclimatic conditions of shallow cave sectors are similar to the ones found outdoor, and these sectors can be highly unsuitable when external conditions are harsh (dry and hot) 27 . On the other hand, food availability is generally the highest close to the surface, therefore the exploitation of these sectors is important, particularly when the need of resources for growth or reproduction is high 29 . Salamander populations inhabiting the areas with most suitable climate (i.e., wet and cool climate) were more frequently found close to the surface than the ones inhabiting the least suitable sites, probably because the outdoor climatic conditions have a direct impact on the microclimate of the sites where salamanders live. Therefore, high climatic suitability favours a prolonged activity and the exploitation of cave sectors where more resources are available 29 . The foraging hypothesis predicts better foraging performance in highly suitable sites, and the high frequency of empty stomachs that we recorded in the less suitable sites confirms that foraging events occurring in these populations may be less frequent compared to high ES sites. Sampling was performed during spring and autumn, which likely are the periods of highest activity 26 . During these seasons, salamanders probably intensify their foraging activity before aestivation and the winter diapause. Several processes can determine differences in foraging success between populations that inhabit different sites. In high ES areas salamanders can devote more time to foraging, while in low ES areas prolonged unsuitable conditions force individuals to remain in shelters 27,29 . Furthermore, it is also possible that resource availability shows spatial heterogeneity, and more resources can be available in the best sites, thus favouring preying activity 37 . Distinguishing between these non-exclusive processes is not easy; moreover, we did not measure variation of prey abundance. However, Hydromantes salamanders are generalist predators consuming a wide range of prey, and are able to feed on most of the underground and outdoor invertebrates 24,25 ; thus we expect that prey is generally present throughout the whole species range. Nevertheless, detailed analyses are needed to assess whether the spatial variation of prey availability determines the differences in feeding performance across populations.
The body condition hypothesis predicts that longer activity period and/or better foraging allow an improvement in the conditions of individuals. Body condition of salamanders showed strong variation among populations, being significantly higher in high ES areas (Fig. 4a). The higher activity and better foraging are likely involved, determining differences in body conditions among salamanders from different areas. However, ES was not the only factor related to BCI variation. In our analyses, body condition was measured using the Residual Index. This index has an excellent performance in limiting the effect of body size on BCI 38 , thus it is not surprising that BCI differences between adults and juveniles were small. Within adults, body condition of females was significantly higher than in males (Fig. 4c), and this probably happens because females accumulate reserves for the breeding activities, which can last several months 39 . Body condition also showed a strong seasonal variation, with better body condition in June (Fig. 4b). Starting from June, outdoor conditions become extremely unsuitable for cave salamanders who move in deep underground shelters where food availability is low 27 . Therefore, Hydromantes likely intensively forages during the previous months in order to store more energy for aestivation 26 , and this explains the good body condition observed in June.
Positive relationships between suitability and activity pattern, feeding and/or body condition, are expected to improve fitness, with potential effects on abundance. Given that the previous hypotheses were confirmed, we expected a positive relationship between the spatial variation and local density across the whole species range. The density of salamanders showed a 30-fold variation among sites and, following our prediction, populations living in high ES areas showed the highest densities (Fig. 3c). We hypothesize that such high density is related to the co-action of multiple processes, such as prolonged activity and better feeding, which in turn improve body condition (Fig. 1). Individual body condition is a key factor which increases individual fitness and, indirectly, can affect species local abundance through different pathways 31 . First, a better body condition makes individuals able to better withstand stressing events (e.g., starvation, adverse environmental conditions), thus improving survival 40,41 . Furthermore, individuals showing better body conditions may devote more resources to breeding and parental care, improving the number and survival of the offspring 42 . These two paths are not mutually exclusive and can work synergistically. Both paths can promote population growth, especially if biological interactions do not represent a major limitation 1 . In the study species, measuring survival and breeding success is challenging, as these amphibians show an elusive breeding behaviour, and monitoring reproduction in nature is rarely possible. For instance, in 40 years of studies on H. flavus, only one egg clutch was observed during a speleological exploration 26,43 . Nevertheless, given that the detection probability of these salamanders in caves is very high, in future studies capture-mark-recapture might provide better information on the differences of survival and individual growth rate across populations 44 .
However, it should be remarked that density-dependent processes occurring at the population level are an additional issue that can complicate the relationships between suitability, local density and fitness. On the one hand, in suboptimal environments densities can be too low (e.g., hampering mating encounters), thus further limiting fitness and reducing abundance below the expected values. On the other hand, in sites with high density competition for resources increases, and this might limit food availability, with potential impacts on body conditions and fitness 45 .
Identifying the processes determining correlations between ES and local species abundance is not easy, as many factors may affect the final outcome 9 . In fact, some studies did not find positive relationships between ES and species abundance 5,7 . Integrative analyses, disentangling the multiple processes that affect population performance, are needed to unravel the complex dynamics acting at a local scale, thus helping to link local-scale population processes to those acting at broader scales, such as range-wide variation of suitability. Measuring multiple population parameters represents a key tool to understand the actual effects of environmental variation on populations. Such approach will allow to move beyond the mere measure of population abundance, to improve our understanding of the variation of fitness and population dynamics across species ranges, and also provide key data to inform explicit and mechanistic modelling of populations 3 . Study species and area. Hydromantes flavus is one of the eight European plethodontid salamanders 46 . It has a small distribution range (<90 km 2 ) and is endemic to the Monte Albo in north-eastern Sardinia 26 , Italy (Fig. 2). Plethodontid salamanders are lungless and breathe mostly through the skin 26 , thus they have a narrow ecophysiological niche, requiring high moisture and relatively cold temperature 26,47 . Such microclimatic requirements are generally found in underground environments (e.g., caves), where the species can be observed throughout the year 27 , especially when outdoor conditions become unsuitable 48,49 . However, H. flavus is not an obligate cave dweller: during suitable seasons it is active outdoor, preying on invertebrates 34,50 . When underground, these salamanders usually occupy sectors not far from the surface, to be closer to food resources 29 . Underground shallow areas are strongly influenced by external climatic conditions, thus external climate influence salamanders even when they are underground 27 . Suitability modelling. We used correlative species distribution models to assess relationships between salamander distribution and major bioclimatic variables, and to obtain measures of broad-scale ES. We considered four bioclimatic variables: annual mean temperature, temperature seasonality, annual precipitation and precipitation seasonality (period: 1979-2013); variables were extracted from the Chelsa-climate dataset at the 30 arc-seconds resolution (approx. 920 × 700 m within the study area), which provides improved climatic estimates in landscapes with complex topography 51 . These variables represent average conditions and their variability across the year, and are major determinants of vertebrate distribution 52,53 . Furthermore, these variables are enough to explain most of the climatic variation 54 , and other important variables (e.g., winter and summer temperatures) are strongly related to linear combinations of the four variables considered. We did not include variables representing the biotic habitat 9 , because these salamanders are mostly related to fine-scale microhabitat variables that are not captured by remote sensing or broad-scale habitat maps 29 . The correlation between variables was weak to moderate (for all pairwise correlations, |r| <0.77). To calibrate models, we used all published presence records of the species (reviewed in 26,55 ), updated with records from our own surveys (Fig. 2). To limit the effect of spatial heterogeneity in sampling efforts, we considered only one presence record per grid cell. SDM were built using five modelling approaches: Generalized Additive Models, Boosting Regression Trees, Classification Trees, Multiple Adaptive Regression Splines and Random Forests. Models were calibrated using a 67% random sample of the presence data and evaluated against the remaining 33% data using the True Skill Statistic (TSS) 56 . To improve reliability of models, we restricted the geographical background to a small region potentially accessible to the species [57][58][59] . Background data were limited to within 17 km from the known presence records (Fig. 2), which is the maximum distance between H. flavus populations belonging to the same genetic cluster 60 . This analysis was repeated five times, thus providing a fivefold cross-validation; models were run using the package biomod2 (version 3.3-7) 61 in R 3.3.3 62 . Given that alternative SDM can provide variable outputs, ensemble forecasting of the different SDM 63 was then used to obtain an overall suitability map. Suitability was calculated as the sum of occurrence-probability projections made by the five modelling techniques run over the three subsamples, weighed by their TSS 64 . It should be remarked that TSS is prevalence independent when making binary predictions, but not when making continuous predictions 65 . We used two approaches to assess the calibration of the model. First, we built calibration plots and calculated the coefficient of determination between the proportion of positive cases and the identity line 32 . Second, as our SDMs were based on presence-only data, we also used POC-plots to build presence-only calibration plots 33 . Calibration was calculated using point-biserial correlation, on the basis of the re-calibrated POC-plot 33 . Surveys. From June 2013 to May 2017, we performed multiple detailed surveys in a total of 14 underground sites (caves) where Hydromantes flavus was known to be present, covering the whole range of the species 26 (Table 1, Fig. 2). Surveys were performed from April to October, the period in which Hydromantes salamanders are more active and show the highest abundance and detectability in underground environments 27,66 . During the first survey, we divided each cave in 3-m longitudinal sectors, and measured the maximum height of the ceiling and the maximum width of each sector 49 . These measures were then used to estimate the explored cave surface. Subsequently, we measured multiple features of the populations: distribution of salamanders within these caves, stomach content, body condition and density. Some of these population features required sampling during specific months and thus, due to sampling constraints (i.e., restricted sampling periods, need of speleological equipment), it was impossible to measure all the population features for all the caves, nevertheless for all the features we gathered data covering the whole range of the species (Table 1).

Distribution of salamanders in caves.
To measure salamander distribution we performed surveys in May-June, as in this period the exploitation of cave environments is the highest 27 . Since several caves received multiple surveys (11 caves surveyed; average: 4.82 surveys per cave; Table 1), for each cave we considered the survey with the highest number of observations. We measured the depth (distance from the cave entrance) of all salamanders detected inside the caves, using a 30-m fibreglass tape meter and a laser-meter (accuracy 2 mm). We used the correlation between the average depth of salamanders (log-transformed) and environmental suitability to test the activity hypothesis. Estimates of average depth might be less robust in caves for which we captured only a few individuals, therefore we repeated such analysis considering only caves with ≥5 individuals detected in at least one survey; results were identical in both cases. Dataset is shown in Supplementary Table S1.
Stomach content. In three populations (Table 1), we performed stomach flushing, which is an unharming technique widely used to check stomach contents of amphibians 67 . For each population, we performed two capture sessions in spring and two in autumn (period: autumn 2015-spring 2017; total: four samplings per site). Stomach flushing was performed using a 5 ml syringe connected to a 1 mm ∅ plastic pipe; the end of the pipe was introduced in the oral cavity of the salamander and 5 ml of water was gently injected in the stomach. Reflux was canalized by a funnel into a plastic jar. The obtained stomach contents were preserved in 75% ethanol and then identified using an optical microscope. Stomach contents were separated in two groups: empty (no items detected) and full (at least one food item was observed) 50 .
To assess if the frequency of empty stomachs differs among caves (i.e., to test the foraging hypothesis), we run the binomial Generalized Linear Models (GLM) using the stomach condition (empty/full) as dependent variable, while season (autumn/spring) and cave identity were used as independent variables. We also used the Bayesian equal-tailed Jeffreys intervals (package MKmisc; 68 ) to estimate 95% CI of the frequency of empty stomachs across populations, as Jeffreys intervals are a robust approach for the estimation of binomial CI 69 . The dataset is published in 50 .
Body condition index. Captured salamanders were weighted (using an electronic scale; precision: 0.01 g) and measured (total length; using a plastic ruler). To include the largest number of sites in this analysis, we also considered a cave where only one salamander was measured (Site 6, Table 1); results did not change if this cave was excluded. For each individual we calculated the Residual Index, which is the difference between the observed and the predicted body mass of animals and is considered among the most reliable body condition indexes 38,70 . To calculate the residual index, we regressed weight against total length of salamanders, and for each individual we extracted the residuals of the regression 38,70 . Weight and length were log-transformed to improve linearity. We considered the total length because these animals often store fat also in the tail 71 . We identified age classes and sex on the basis of secondary sexual characters and body size. Salamanders with male sexual characters (mental glands and premaxillary teeth) were considered adult males; individuals without male characters but ≥80 mm were considered adult females (80 mm is the size of the smallest observed adult males); individuals <80 mm were considered juveniles 26 .
Populations were surveyed in different periods of the year (from October 2015 to April 2016), and salamander body features were measured in 20 surveys performed in 11 populations (average: 1.8 surveys per site). To avoid pseudoreplication, for each population we only considered one survey per month, selecting the one with the highest number of measured salamanders, as individual identification was impossible in most of cases 72 . Body weight can show seasonal variation. To test the body condition hypothesis, we ran a linear mixed model (package lme4; 73 ) considering BCI of salamanders as dependent variable, ES, month of survey and group (m/f/j) as independent variables, and population identity as random factor. Sample size was not homogeneous among populations, thus residual degrees of freedom were approximated following Satterthwaite 74 . Given that we detected significant BCI differences between the tree groups, we used orthogonal contrast to test whether there are differences between juvenile and adults and, within adults, between males and females. Dataset is shown in Supplementary Table S2.
Population abundance and density. In seven caves we performed repeated visual encounter surveys in a short period, to estimate salamander abundance. During each survey, the same person dedicated 7.5 min. of observation to each 3-m long longitudinal sectors counting the active salamanders 66 . Each cave was surveyed five times in 2016, from May 4 th to June 27 th . In this period, cave occupancy and detection probability of salamanders is the highest, immigration/emigration to/from the cave is minimum, and no hatches are known to occur 27,43 , thus allowing to meet the closed population assumptions of N-mixture models 75 .
We used N-mixture models to estimate population size on the basis of repeated counts 76 . This approach provides accurate estimates of actual population size, particularly in species with high detection probability such as cave salamanders 77 . We used Akaike's Information Criterion (AIC) to select the most appropriate error SCIENTIFIC REPORtS | (2018) 8:7527 | DOI:10.1038/s41598-018-25704-1 distribution (Poisson or zero-inflated Poisson); negative binomial models were not considered as they can produce infinite abundance estimates 78 . We used Empirical Bayes methods to estimate the posterior distribution of the abundance (mean and 95% Bayesian credible intervals, CI) 79 . Surveyed surface was strongly variable among caves (33-1229.61 m 2 ), therefore we calculated salamander density to allow comparison among sites. We calculated population density on the basis of abundance estimates and the surveyed surface of caves, and then assessed the correlation between population density (square-root transformed) and ES. N-mixture models were run using the unmarked package in R 80 . Dataset is shown in Supplementary Table S3.
Spatial autocorrelation can influence the outcome of statistical tests, therefore we used Moran's I to assess the occurrence of residual autocorrelation in all regression and correlation analyses 81 . For all regression models, residual autocorrelation was weak and not significant (in all tests, Moran's I ≤ 0.2, all P > 0.05), suggesting that spatial autocorrelation was not an issue to our results. For SDMs, spatial autocorrelation of residuals was significant but weak (Moran's I = 0.12, P < 0.01).