Cave features, seasonality and subterranean distribution of non-obligate cave dwellers

Surveys

We surveyed 16 caves in the Northern Tuscan Apennines (Central Italy, between 43°52′29″N, 11°13′04″E and 43°59′51″N, 10°13′48″E) monthly from January to December of 2013. We completely explored 13 sites, while the remaining three caves were explored until the point at which speleological equipment was indispensable and exploration became too difficult (minimum explored length = 6 m, maximum = 60 m, average ± SE = 23.44 ± 3.94 m). All caves were divided into sectors of 3 m in length, which resulted in 125 cave sectors, considering all the caves (Table S1). We decided to divide the inner environment into three meters cave sectors, as this subdivision allows adequate characterization of microclimatic variation within caves (Lunghi, Manenti & Ficetola, 2015). During each survey, we recorded morphological and environmental features for each sector. On a few occasions, some sectors were not accessible because of temporary inaccessibility (e.g., flooding) and missing data represents 2.7% of the observations. At the end of each sector, we measured maximum height and width using a metric wheel. Furthermore, we estimated average wall heterogeneity (i.e., richness of clefts; Camp & Jensen, 2007) by placing a string of 1 m in the middle of the sector and measuring the distance between the two string extremities using a measuring tape. We performed the measures for both right and left walls and calculated the average (Lunghi, Manenti & Ficetola, 2014). Air temperature and humidity were recorded using a Lafayette TDP92 thermo-hygrometer (accuracy: 0.1 °C and 0.1%) adopting precautions to avoid influence on cave microclimate (Lopes Ferreira et al., 2015), while the average incident light of sectors was estimated using a Velleman DVM1300 light meter (minimum recordable light: 0.1 lux).

Detection of species

We performed 12 monthly surveys in each cave and the interval between two consecutive surveys in the same cave varied between nine and 45 days, to have surveys for every month of the year. All 3 m-sectors were re-sampled by the same person, who dedicated 7.5 min to each sector to assess the presence of species using visual encounter surveys (Crump & Scott, 1994).

During the surveys, we assessed the presence and distribution of eight trogloxenic and troglophilic taxa within each cave sector: one orthopteran (Dolichopoda laetitiae), three spiders (Meta menardi, Metellina merianae, and Tegenaria sp.), two gastropods (Chilostoma planospira and Limax sp.) and two anuran amphibians (Rana italica and Bufo bufo). The species were selected due to their common presence in caves and their easy identification without invasive methods, giving us the opportunity to collect data limiting disturbance and changes in individual behaviour. Furthermore these species show significant variation in life history and ecological traits. The cave cricket Dolichopoda laetitiae is generally abundant in subterranean environments of Central Italy (Allegrucci et al., 2014). These crickets are scavengers of subterranean environments but can also feed outside caves, thus increasing the supply of allochthonous organic matter into the caves (Lavoie, Helf & Poulson, 2007; Weckerly, 2012). The three arachnid species exploit subterranean environments for a significant part of their lives and are among the major cave predators (Manenti, Lunghi & Ficetola, 2015; Novak et al., 2010). We grouped all detected individuals as Tegenaria sp. because identifying spiders at the species level in the genus Tegenaria is difficult without using invasive methods (Bolzern, Burckhardt & Hänggi, 2013). Chilostoma planospira is a snail that feeds on vegetal organic matter and even lichens (Albano et al., 2014; Baur et al., 2007), while the genus Limax contains large slugs (Limax were only identified to genus) that feed on residual organic matter (Horsák, Zelený & Hájek, 2014). The frog Rana italica and the toad Bufo bufo frequently exploit caves but mostly use these habitats as shelter (Bonini, Razzetti & Barbieri, 1999; Lunghi, Manenti & Ficetola, 2014).

Statistical analyses

The detection of species is imperfect, and detection probability <1 can bias the results of standard regression analyses (MacKenzie et al., 2006). Approaches have been developed to take this into account in regression models (MacKenzie et al., 2006), but they generally require multiple surveys being performed within periods without changes due to birth, death, and movement of individuals to and from the population (closed populations; MacKenzie et al., 2006). Considering the time interval between most of surveys, the assumption of closed populations was likely not met (MacKenzie et al., 2006), making the application of these methods difficult. We therefore used two complementary approaches to control for the potential effects of imperfect detection. First, caves were surveyed using a standardized study design to limit variation in detection rate among sectors, and we dedicated the same sampling effort to all sectors during all seasons. Under standardized monitoring efforts, analyses taking and not taking into account imperfect detection generally produce consistent results (Banks-Leite et al., 2014). We thus analyzed data with generalized linear mixed models (GLMMs). This approach allows for heterogeneity of species distribution among caves and sectors, even though it does not integrate imperfect detection. Subsequently, we repeated analyses using approaches allowing to integrate imperfect detection, to assess the robustness of GLMMs analyses.

Analysis with mixed models

We identified the relationship between species and cave features using GLMMs assuming binomial error. The presence of each species was considered as a dependent variable while cave features (microclimatic and morphological) were considered as independent variables. As each cave was visited 12 times, cave and sector identity were used as random factors, while month of survey was included as a categorical variable. To test if species tend to frequent different microclimatic zones during a given time, we also considered the interaction between month (sampling period) and microclimatic features as independent variables. All possible combinations of independent variables were built and ranked using the Akaike’s Information Criterion corrected for small sample size (AICc) (Fang, 2011). For each species, we first built GLMMs using all possible combinations of independent variables and the model with the lowest AICc values was considered the best AICc model. Complex models showing AICc values greater than simpler, nested models were not considered in the set of candidate models (Richards, Whittingham & Stephens, 2011). Variables were transformed using logarithms or square root to better fit the normal distribution when necessary. Likelihood ratio tests were used to assess the significance in terms of the best AICc models.

Analyses taking into account imperfect detection

Species tend to move across sectors during different periods of the year because of the change of microhabitats (see results). However, analyses integrating detection probability generally assume that multiple surveys are performed during short periods in which occupancy is constant (MacKenzie et al., 2006). We therefore used occupancy modelling to calculate the detection probability for each species, on the basis of pairs of surveys that were performed within 14 days, as these periods may better ensure meeting the assumption of closed populations (MacKenzie et al., 2006). To calculate detection probability of species we used 35 pairs of surveys, performed in all caves and across seasons within an interval which did not exceed 14 days. Subsequently, given that we generally performed only one survey per month, we analyzed data using standard general linear models, while weighing absences on the basis of the detection probability of species (Gómez-Rodríguez et al., 2012). Such an approach allows us to integrate imperfect detection even when only one survey per period is performed (Gómez-Rodríguez et al., 2012).

Integrating random factors into this analysis is not possible, thus, cave identity was included as a fixed factor. For each species, we built GLMs with all the possible combinations of independent variables, and ranked them following AICc. We tested the significance of variables included in the best-AICc model using a likelihood ratio test (Bolker et al., 2008).

Species richness

To analyze the effect of seasonality on species richness, we used GLMMs with Poisson error. In addition to the eight focal taxa, we also integrated data on the distribution of two species that were frequently observed in the caves, but were not considered in single-species analyses: the salamander Hydromantes italicus and the spider Amaurobius ferox. The salamander Hydromantes italicus was the focus of a previously published study performed in the same area (Lunghi, Manenti & Ficetola, 2015) and represents one of the top predators of the studied taxa. The spider Amaurobius ferox was recorded during surveys, the number of detections (mean: 3.67 detections per month) was too small for building robust distribution models. Species richness calculated for each sector represented the dependent variable, while microclimatic features (temperature, humidity and light) and their relative interaction with the sampling period (month) were the independent variables. The individual caves and their sectors were used as random factors. All statistical analyses were performed with program R using packages lme4, MuMIn, MASS and unmarked (Bartoń, 2016; Bates et al., 2015; Fiske & Chandler, 2011; R Core Team, 2016; Venables & Ripley, 2002).

(a) Species distribution

The distribution of species varied among cave sectors over time. For all taxa, we detected significant relationships between occurrences, survey period, and the environmental features recorded (Table 1), but these relationships were highly variable amongst species. For all taxa, results of both mixed models and GLMs considering imperfect detection were generally consistent, although some relationships were detected by just one of these analyses. The occurrence of the cricket D. laetitiae was negatively related to sector height and light, and positively related to wall heterogeneity and temperature. The best AICc models also included the interactions between month and humidity (Table 2A; Table S3A). From winter to early spring this cricket occupied sectors characterized by high humidity, while in late summer and autumn it was associated with less humid sectors (Fig. 1A). Furthermore, the best mixed model also included the interaction between month and light (Table 2A), while the best model considering imperfect detection also included a positive relationship with width of sector and the interaction between month and temperature (Table S3A).

Within the arachnofauna, occurrence of M. menardi was positively related with sector width, and negatively with height and light. Interactions between month and light and between month and temperature were also included in the best AICc model (Table 2B; Table S3B). From late summer to winter this spider occupied sectors with warmer temperatures, while sectors with low light were preferred during the whole year (Figs. 1B and 1C). Furthermore, the best mixed model also included a positive relationship with sector humidity (Table 2B). The frequency of M. menardi was higher in autumn. The occurrence of M. merianae was positively related with sector temperature; the highest frequencies of the species were detected between spring and autumn (Table 2C; Table S3C). The best mixed model also suggested a negative relationship with sector height and light, and detected an interaction between month and light (Table 2C), while the best model considering imperfect detection also included a negative relationship with humidity (Table S3C). The occurrence of Tegenaria was negatively related to light. Furthermore, there was a significant interaction between month and light (Table 2D; Table S3D), as the relationships between these spiders and low light was particularly evident in winter (Fig. 1D), a period in which the species was particularly frequent in caves. Furthermore, the best mixed model also included a positive relationship with sector temperature and a negative one with width (Table 2D), while the best GLM also included a negative relationship with sector humidity and height, and the interaction between humidity and month (Table S3D).

The interaction between month and microclimatic features was never included in the best-AICc model of gastropods, while the month of survey was always included (Tables 2E and 2F). The occurrence of C. planospira was positively related with sector height (Table 2E; Table S3E) and occurrences were higher during spring and summer. The best mixed models also included a negative relationship with sector width and light (Table 2E). For Limax slugs we did not detect any relationship with the examined cave features, however, its occurrence was higher during summer and early autumn (Table 2F; Table S3F).

The occurrence of Rana italica was negatively related to light and wall heterogeneity, and positively related to sector width and height (Table 2G ; Table S3G). The best mixed model also included the interaction between month and light (Table 2G), while the best model considering imperfect detection included a negative relationship with sector humidity (Table S3G). Finally, the occurrence of B. bufo was positively related to sector height (Table 2H; Table S3H); the best mixed model included also a negative relationship with temperature (Table 2H).

In this study, we did not include cave depth as an additional variable in models assessing species distribution, because distance from the entrance was strongly related to multiple features of caves, such as light and humidity, and including strongly correlated variables into regression models can bias the regression outcomes (Dormann et al., 2013), making model selection particularly hard (Cade, 2015). For instance, the negative correlation between sector and incident light was very strong (r =  − 0.94). However, distance from cave entrance often affects the distribution of non-trogblobitic species (Lavoie, Helf & Poulson, 2007), thus we repeated analyses adding distance from the entrance to the best AICc-models. Adding distance reduced the AICc values of the best models of all species except Bufo bufo, and the relationship between depth and species occurrence was negative for most of taxa (Table S5). This suggests that the majority of study taxa are actually found nearby the cave entrance. Nevertheless, adding cave depth to our models did not modify the coefficients of relationships between species and microhabitat features (Table S5), therefore our conclusions are not affected by the decision to not include depth in regression models.

(b) Species richness

Species richness was mainly related to microclimatic features (Table S4). The richest communities were positively related to sector humidity (${\chi }_1^2=9.7$, P = 0.002) and temperature (${\chi }_1^2=98.3$, P < 0.001), while negatively related with light (${\chi }_1^2=62.8$, P < 0.001). These relationships determined a shift of the sectors with highest richness through the year. During winter, species richness was more evenly distributed among all cave sectors (Fig. 2A), while in spring and summer the highest richness tended to be more concentrated toward the cave entrance (Figs. 2B–2C).

Species distributions

All the species tend to occupy areas in which their physiological requirements are satisfied (Kearney & Porter, 2009), thus they are associated with areas where the environmental conditions (e.g., microclimate) are suitable. Survey month was a major factor determining the presence of species in caves, as it was present in nearly all best models (Table 2 and Table S3). This fact underlines the high seasonality of species distributions, even for populations exploiting subterranean habitats. In the extreme case of Limax slugs, month was the only variable included in the best AICc model (Table 2F). These slugs did not show any significant relationship with the considered microhabitat features, but the occurrence of these gastropods was strongly affected by period of the year. For nearly all the species, the highest occurrence was detected during spring and summer months, periods in which ectothermal animals are generally more active, and food availability is high (Bale & Hayward, 2010).

Physical and morphological cave features (e.g., width, height and wall heterogeneity of inner environment) were very important, especially for invertebrate predators. Meta menardi showed an association with sectors characterized by a low cave roof. This predator hunts through both active search and sit-and-wait strategies, so constructing its web in sectors where the ceiling is not far from the ground might increase rates of prey capture in webs. However, cave morphology did not only influence predators. Non-predator species could also be influenced by sector structure, such as Dolichopoda laetitiae and Chilostoma planospira. The cricket D. laetitiae was associated with sectors characterized by high irregularity, and many clefts which may have given the opportunity to better avoid predators. Chilostoma planospira was related to sectors in which food availability might be higher, as this species frequented areas characterized by high passage which in turn have more wall surface, as this feature is positively related with the abundance of lichens (Manenti, 2014).

The occurrence of seven out of eight study taxa was significantly influenced by the month of survey; the only exception was Bufo bufo. This widespread toad has a wide physiological tolerance, thus, the limited microclimatic excursion of cave environments may have a limited effect on the presence of this species. Furthermore, in three studied taxa we detected clear interactions between month and microclimatic features (Table 2, Fig. 1). These interactions may occur for several reasons. In some cases the species may select the same microhabitat across the seasons and mitigate within the caves, selecting areas with suitable microclimate over the year. For instance, a previous study on the same group of caves showed that cave salamanders are consistently associated with sectors having a temperature of 10−15 °C (Lunghi, Manenti & Ficetola, 2015); in summer these conditions represent the coolest sector, while in winter these are the conditions of warmest sectors. As a consequence, the relationship between salamander occurrence and temperature is positive during winter and negative during summer (Lunghi, Manenti & Ficetola, 2015). However, the significant interaction between environmental variables and month suggests taxa were associated with different microhabitats during different periods. Species were associated with the darkest or most humid sectors in winter, while in summer they were found in superficial sectors with more light and less humidity (Fig. 1). Warmer sectors were more often inhabited during winter (Fig. 1). For most of the species, activity peaks are reached in summer, when individuals are generally closer to the surface. Conversely, during winter periods, species become less active and tend to occupy deeper sectors of caves, where prolonged suitable conditions may let them reduce costs linked to acclimatization (Angilletta Jr, Niewiarowski & Navas, 2002) (Fig. 1).

Our study shows that at least for the populations considered, cave exploitation is not random and occurs during the whole year. Moreover, at least for the relatively superficial cave sectors, the distribution of species is strongly related to the different cave microclimatic features, and to their seasonal variation during the year. Seasonal effects are extremely important, and are a major driver of species distribution (Sheldon & Tewksbury, 2014). In caves, the seasonality occurring in shallow sectors is generally higher than in the deepest parts, where temperature shows very limited variation (Ravbar & Kosutnik, 2014). Nevertheless, even at the cave entrance microclimatic conditions are significantly more stable than in outdoor environments (Ray, Beever & Rodhouse, 2016).

Species richness

During the year, the highest species richness was found in sectors with specific microclimatic conditions (i.e., high humidity, warm temperature and low light). These conditions are usually realized in sectors that are not far from the entrance, but are deep enough to prevent excessive light exposure (Culver & White, 2005). However, the cave area in which such a combination is realized is not static, thus its location is not fixed but changes depending on surface environmental influence. For example, the stronger the external influence, the deeper the suitable area. In such zones, individuals are protected from extreme external thermic fluctuations but, when external condition are suitable, they are able to reach external foraging sites (Fang, 2011). Distance from the entrance is certainly important for the distribution of most of study species (Table S5), but it is strongly related to the microclimatic features, which actually explain most of variation in species distribution (Table S5). Therefore, we did not include distance as predictor in our models, given that focusing on microclimatic features allows to better characterize the most likely determinant of species

The distribution of sectors with the highest species richness seems to be bimodal, with peaks of richness at about 15 and 48 m from the entrance (Fig. 2).The farther peak of richness might be due to the presence of cryptic access to the main cave environment. In cave studies the main entrances (where humans can get in) are generally considered to be the sole accesses to inner environments. However, even if that is true for humans, small animals can use small fractures and passages to enter cave environments. In fact, in the two longest examined caves, respectively at 48 and 57 m of depth, the presence of interstitial openings likely permits species to enter and exit from cave environments without passing through the main entrance.

Study limitations

In our study we standardized surveys in order to minimize the impact of variation of detection probability on our results (Banks-Leite et al., 2014). Nevertheless the difficulty of applying occupancy modeling to our data can complicate data interpretation. Indeed, in the results it is not easy to distinguish between true variation in occupancy and variation in detection probability.

To confirm the appropriateness of results, we repeated analyses by integrating the detection probability of each species. The outcome confirmed that mixed model were generally robust, as 70% of relationships detected by mixed models were confirmed by analyses considering imperfect detection. When sites are surveyed only once per period, weighted GLMs can allow us to take into account imperfect detection (Gómez-Rodríguez et al., 2012). However, this approach also has some limitations. Given that standard implementations do not allow including random factors, it hampers our ability to fully consider the effects of spatial heterogeneity, and the potential impacts of cave and sector identity on species distribution. Under these situations, a comparison between the two methods can be the best approach to identify the most consistent relationships, thus obtaining more robust results.