The study was set in the Pugnetto hypogean complex, in the nearby to the hamlet of Pugnetto, municipality of Mezzenile, Lanzo Valley, Graian Alps, Piedmont (NW-Italy). The site is protected under the European Habitat Directive 43/92 (S.C.I. IT 1110048) and hosts five natural caves classified as “Caves not open to the public” (H 8310), namely the Borna di Pugnetto (cadastre number Pi 1501, entrance at N45°16'19", E5°02'26"; altitude 820 m a.s.l.), the Tana del Lupo (Pi 1502, N45°16'19", E5°02'22"; 813 m a.s.l.), the Creusa d‘le Tampe (Pi 1503, N45°16'12", E5°02'33"; 870 m a.s.l.), the Tana della Volpe (Pi 1504, N45°16'13", E5°02'34", 885 m a.s.l.) and the Cavernetta (N45°16'17", E5°02'36", 895 m a.s.l.) which still lacks a cadastre number but was mentioned in the original description of the hypogean site by Muratore (1946).

Our study was conducted in the Borna di Pugnetto (hereinafter Borna) and in the Creusa d‘le Tampe (hereinafter Creusa) caves (Figs 1a, b, d; 2), as well as in the surrounding MSS. The Borna has a planimetric development of 765 m. It consists of a main gallery of 300 m directed north-south, which splits at the end into two main branches – “Ramo della Madonna” and “Ramo della Fontana”. The Creusa has a planimetric development of 47 m. It consists of a single tunnel, which leads to the final chamber after a narrow passage. The area surrounding the caves is characterized by a well-developed colluvial MSS made of fragmented rocks (calcschist, marble and phyllite) which form a network of habitat spaces with sizes between 0.1 and 10–20 cm, covered by 30–40 cm of leaf litter and soil (Fig. 1c, e). Surface conditions are characterized by a close and shaded Luzulo-Fagetum beech forest (Sindaco et al. 2009).

Figure 1. 

a Main entrance of the Borna di Pugnetto (photo credit: Alberto Chiarle and Mauro Paschetta, 2014) b Main entrance of the Creusa d’le Tampe (photo credit: Elena Piano, 2013) c exposed soil/MSS profile in a fresh-cut along a slope in the vicinity of the Borna di Pugnetto (photo credit: Jacopo Orlandini, 2014) d the typical cave geo-morphology within the Borna di Pugnetto (photo credit: Alberto Chiarle and Mauro Paschetta, 2014) e detail of the MSS geo-morphological structure (photo credit: Jacopo Orlandini, 2014).

Figure 2. 

Map of the study area. The shape and the topographic position of the four caves (Borna Maggiore di Pugnetto, Tana del Lupo, Creusa d’le Tampe, Tana della Volpe) was obtained from the original planimetric drawings of Muratore (1946). The position of the sampling plots in caves (“cave triplets”, C1–C8), in the MSS (“MSS triplets”, M1–M8) and in the leaf litter (“epigean”, L1–L6) are represented by coloured dots. The different sectors of the cave are coloured with different shades of grey representing the subjacency – i.e., vertical distance from the surface – according to Motta and Motta (2015).

We used twenty-four pitfall traps (diameter 9 cm, volume 40 ml) to collect invertebrates in caves. The pitfall traps were arranged in groups of three (hereinafter cave-triplets), at a distance of ~5 m to one another. Six cave-triplets (code: C1–C6) were placed in the Borna at 4, 30, 90, 150, 230 (“Ramo della Madonna”) and 350 (“Ramo della Fontana”) meters from the main entrance. Two cave-triplets (code: C7, C8) were placed in the Creusa at 5 and 25 meters from the main entrance (Fig. 2). We derived the subjacency (DEPTH; i.e., the direct vertical distance from the surface) of each cave triplet from the geological survey of Motta and Motta (2015). Specifically, we created three categories of subjacency (Fig. 2): 0–20 m (triplets C1, C5, C7, C8), 20–40 m (C2, C6), and 60–80 m (C3, C4) [see Mammola et al. (2015) for details].

For collecting invertebrates in the MSS, we utilized twenty-four Subterranean Sampling Devices (SSD; after López and Oromí 2010). Since the description of the MSS by Juberthie et al. (1980a, 1980b, 1981) and Uéno (1980), there has been considerable effort to construct and improve effective sampling techniques in this habitat (reviewed in López and Oromí 2010, Domingo-Quero and Alonso-Zarangara 2010, Mammola et al. 2016, Růžička and Dolanský 2016). In this study, we constructed SSDs relying on the prototypes of Nitzu et al. (2010, 2014) and Deltshev et al. (2011). SSDs consisted of PVC pipe (diameter 12.5 cm) with a silicone cap closing the top. A straight row of 25 holes (diameter: ~1.5 cm) was pierced at a distance of 8 cm from the lower-end of each pipe (Fig. 3a). At the same height, in correspondence of the level of holes row, we placed a standard pitfall trap. The holes drilled along the surface of the PVC pipe allowed the fauna to access to the inside and eventually fall into the pitfall trap; at the same time, our SSD allowed to recover, empty, refill and replace the pitfall trap trough the closing cap, ensuring that the substrate was not mechanically excavated at each sampling session (Fig. 3f, g). In addition, we inserted a styrofoam cylinder within the PVC tube, filling completely the pipe and insulating the trap from the external meteorological fluctuations. A screw was placed just above the sampling holes, to prevent the styrofoam cylinder from sliding down (Fig. 3b).

Figure 3. 

a sampling holes (details) b Blocking screw c installation of an MSS-triplet d SSD of three different length e MSS-triplet buried in the ground f, g renewing the pitfall trap inside the SSD. Photo credits: Elena Piano.

Three adjacent SSDs of three different length (40, 60 and 80 cm) were buried vertically in the ground (hereinafter MSS-triplet; Fig. 3c, d, e). The 40 cm SSD intercepted the substrate at the interface soil/MSS, while the 60 and 80 cm SSDs intercepted the MSS. Six MSS-triplets (18 SSDs; code: M1–M6) were set in the nearby of the Borna, and two MSS-triplets (6 SSDs; code: M7, M8) were placed in the nearby of the Creusa (Fig. 2). We installed the traps in the MSS following the advices reported in literature (Domingo-Quero and Alonso-Zarangara 2010, Giachino and Vailati 2010, López and Oromí 2010, Mammola et al. 2016; see Fig. 3). Pitfall traps in caves and in the MSS were baited with chicken meat and filled with brine (supersaturated preserving solution of water and NaCl; Giachino and Vailati 2010).

Six pitfall traps (with brine, not baited) were also placed in leaf litter habitat (epigean; code: L1–L6), as a reference to discriminate correctly between the epigean fauna and the specialized fauna – see paragraph “Specimens sorting and ecological classification” for details.

Figure 2 shows the position of the traps in the area. We replaced the traps approximately once a month from June 2012 to June 2013. However, the access to the Borna is forbidden from 1st November to 31th March in order to protect the roosting bats inhabiting the cave. In order to reduce disturbance, in this period we accessed the Borna only twice (December and March), resulting in ten total sampling sessions over the year. In order to evaluate the effect of the bait on the sampling probability, we replaced it every two sampling sessions, thus resulting in 5 sampling sessions with fresh and 5 with exhausted baits.

To characterize the subterranean microclimate, we placed one Hygrochron temperature and humidity datalogger in correspondence of each pitfall trap in the cave and at the lower-end of each SSD. Hygrochron were programmed to sample temperature (T) and relative humidity (RH) every three hours for the whole sampling period (accuracy of ± 0.5°C and ± 1%, respectively). We also derived the mean daily outside temperature to the same periods from the nearest thermo-hygro-pluviometric weather station (Fua, Lanzo Torinese, Cod. 111; N45°17'23", E7°29'38"; 550 m a.s.l.). The temperature values recorded by the weather station were corrected with the standard environmental lapse rate – the change of temperature with altitude for the stationary atmosphere. In all analyses relating the abundance and species richness of external and adapted elements with the explanatory parameters (see later sections), we used the pseudo-replicates of each trap as basic sample units.

Trapped individual were sorted, identified and classified either as epigean (category: “external”) or subterranean elements (category: “adapted”). According to our expertise (MI and SM: Araneae; PMG: Coleoptera and some other orders of Insecta) and the availability of specialists for additional taxa (see Acknowledgments), identification of the species levels was possible for the orders Araneae, Opiliones, Pseudoscorpiones, Chilopoda, Isopoda (one species), and for most orders of insects (especially Coleoptera; see Appendix 1).

In subterranean biology, species are often classified into ecological categories (e.g. trogloxenes, troglophiles, troglobionts) according to their preferred habitat and general association with the subterranean domain (Sket 2008, Trajano and Carvalho 2017). However, most ecological classifications, such as the Schiner-Racovitza system, are deemed to oversimplify real cases given that boundaries between the categories are often vague and difficult to attribute (Martìn et al. 2001, Giachino and Vailati 2016). Alternatively, in this work, to categorise our samples we followed the general idea based on a source-sink population model recently proposed by Trajano (2012). By examining the prevalence of individuals either in the leaf litter traps or in the traps installed in the subterranean habitats, we were able to classify species as “external” (species likely having source populations in the epigean environment) or “adapted” (species likely having source populations in the hypogean environment). Following this criterion, we calculated the abundance and species richness of external and subterranean elements (Next, Rext, Nad, Rad, respectively) for each trap in each sampling survey. Species represented by less than three individuals were not classified, and thus excluded from the analysis, unless if we possessed a solid literature background testifying the species ecological requirements – e.g., Lana (2001) for a few subterranean species in Piedmont and Mammola et al. (2017) for most spider species.

We performed regression-type analysis following the general advices of Zuur and Ieno (2016). All analysis were conducted in R (R development core team 2015). In order to test statistically our working hypothesis, we computed three distinct sets of regression models. Firstly, we compared the cave and the MSS habitats in term of species richness and abundance of external and adapted organisms (aim 1). Subsequently, we analysed the spatial and temporal gradients in the MSS and cave environments separately, aiming at investigate the relevant abiotic factors driving richness and abundance in these two hypogean compartments (aims 2 and 3).

The first set of models was computed through mixed-design analysis of variance with Poisson distributed data (i.e. generalized linear mixed models, GLMMs; equation 1), whereas for the second (equation 2) and third (equation 3) sets, we primarily relied on Poisson generalized additive mixed models (GAMMs). The mixed part of both GLMMs and GAMMs was introduced in order to account for multiple observations from the same triplet over time, by specifying the triplet as random factor. The latter variable was included as random factor in order to account for the variation it introduced in our samples – and thus to correctly estimate the regression coefficients, – rather than to test for its direct effect on the dependent variables.

Prior to model fitting, we explored the three datasets following the standard protocol for data exploration proposed by Zuur et al. (2010). According to Zuur et al. (2009, 2010), the inclusion of outliers and highly correlated predictors in the regression analysis leads to misleading results – type I and II statistical errors. We thus used Cleveland‘s dotplots to assess the presence of outliers in dependent and independent variables and we investigated multi-collinearity among covariates. For each basic sample unit and for each set of models, we chose the abundance of external elements (Next), the species richness of external elements (Rext) and the abundance of adapted elements (Nad) as dependent variables. Such variables were selected aiming at investigate the gradient of specialization trough the hypogean environment (aims 2 and 3). Instead, considering the generally low species richness of adapted elements per trap (Rad), we were forced to exclude this latter variable from our analysis. For both abundance (Next, Nad) and richness (Rext) count data we assumed a Poisson distribution, but we tested for over-dispersion after model fitting and we switched to a negative binomial distribution when the over-dispersion parameter was higher than 1.2.

Comparing MSS and cave communities (aim 1)

Species richness of external elements and abundance of adapted and external elements were analysed in relation to the explanatory two-levels factor HABITAT (levels: “Cave” and “MSS”) using the glmer R command in the package ‘lme4’ (Bates et al. 2013).

The fixed structure of the model was:

(1) y ~ HABITAT

Where y is one of Next, Rext and Nad. The random part of the model allowed us to deal with repeated observations and measurements of the same triplet (temporal dependence) and the clumped distribution of the traps within the triplet (spatial dependence).

Spatial and temporal gradients in MSS (aim 2) and cave (3) communities

Species richness of external elements and abundance of adapted and external elements were analysed in relation to the explanatory variables using the gamm R command in the package ‘mgcv’ (Wood 2015). Generalized additive mixed model were used in order to account for possible non-linear trends of the sampling period (Serie_i; continuous variable). The optimum amount of smoothing was estimated trough generalized cross-validation (GCV). However, whenever the effect was linear and/or not significant, we dropped the smoothed term and fitted a new model including only parametric terms, thus adopting a linear approach (Poisson or negative binomial GLMM).

For both the cave and the MSS compartments, we also included the two-levels factor BAIT (level= “Fresh” and “Exhausted”) as covariates and the two-levels factor SITE (levels= “Borna” and “Creusa”). The first variable was included to evaluate the effect of the ageing of the bait on our dependent variables. The second variable was introduced to take into account for possible local effects, since we pooled together records from two different caves and associated MSS.

For the analysis of the MSS compartment, in addition to the above mentioned variables we also included in the models the three-level categorical variable sampling depth (DEPTH; levels: “0.4m”, “0.6m” and “0.8m”). For the cave habitat, in addition to the aforementioned variables, we also included the distance from the cave entrance (Dst_i; continuous variable) and the three-levels categorical variable subjacency (SUBJ; levels: “0–20m”, “20–40m” and “60–80m”). We excluded from the regression analysis the microclimatic variables, given that due to malfunctioning, several dataloggers did not recorded reliable measurements of temperature and relative humidity during the sampling period (see results, further details in Mammola et al. 2015).

The (fixed) structures of the initial models were (aims 2 and 3, respectively):

(2) y ~ DEPTH + s(Serie_i) + BAIT + SITE

(3) y ~ Dst_i + SUBJ + s(Serie_i) + BAIT + SITE

Where y is one of Next, Rext and Nad, and s(Serie_i) indicate the smoothing term. The random part of the model is equal to the previous model (1). For the 2^nd and 3^rd models, we adopted a statistical hypothesis testing framework, whereby model reduction was carried out on the full model by sequentially deleting non-significant terms and potential interactions according to AIC values (Zuur et al. 2009), until a minimum adequate model of significant fixed effects remained (i.e. best model supported by observations). For all models, p-values for parametric term were based on z-tests; t-test in the case of GAMMs and GLMM with negative binomial distribution. P-values for smoother terms represent approximate significance based on F statistic. Model validation was carried out on the final models, following the approach of Zuur et al. (2009).

Approximately 15,700 Arthropoda (20 Orders), 117 Mollusca (order Pulmonata), 868 Crustacea (order Isopoda), and 14 Anellida (Lumbricus sp.) were collected. We report the complete list of the taxa and number of specimen collected in each habitat in Appendix 1. Diptera were the most represented, with 7,469 individuals collected (adults and larvae), followed by Coleoptera (3,543 individuals; adults and larvae), Collembola (1,755 individuals), Hymenoptera (1,046 individuals, mostly Formicidae) and Acari (631 individuals). MSS showed higher values of order diversity compared to caves, in which, abundance of specimens was higher.

Due to malfunction, several dataloggers did not sample either temperature or relative humidity during the sampling period. Specifically, nine out of eighteen in the Borna, two out of six in the Creusa and ten out of twenty-four in the MSS. Due to this significant loss of data, it was not possible to include climatic data in the regression analyses.

Relative humidity in the caves proved to be almost constantly close to saturation, with mean monthly values ranging from 85 to 100 %. However, in the vicinity of the entrance the relative humidity dropped down to 70–75 % in winter. In the MSS the mean monthly relative humidity was always above 90 %. With regard to temperature, changes and max–min ranges were attenuated with increasing distance from the entrance and delayed in respect to the outside values. The mean annual temperature values deep inside the two caves was comparable (Tmean ± SD: Borna = 9.0 ± 0.4 °C; Creusa = 8.9 ± 0.8 °C) and presented little variations over the year. In the outermost sections, temperature was nearly stable in summer, spring and autumn, with a max–min range around 4.5 °C in the vicinity of the entrance for both caves. However, the microclimate at the entrance zone drastically changed during winter, when we observed a drop in the mean temperature values (mean values always below 6 °C). The coldest temperature values were recorded in December and January (Tmin: Borna = –2.0 °C and Creusa = –0.9 °C).

The range of temperature variation, both daily and monthly, was lower in the MSS in respect to surface (Fig. 4). In the MSS the thermal variability at the surface seemed to moderately affect the microclimatic conditions of MSS, resulting in seasonal variations in temperature at all sampled depths. The maximum temperature in the MSS did not exceed 19 °C even though air temperatures reached as high as 32.3 °C (June). Minimum temperature values rarely fell below 3 °C, with occasional drops below 0 °C (lower Tmin recorded in the MSS in December: –2.0 °C).

Figure 4. 

Annual trends of temperatures in the Pugnetto hypogean complex. The shade of blues indicate the relative position of the dataloggers at each cave-triplet, from the outermost (lighter blues) to the innermost sections (darker blues). Records from only one MSS-triplet are shown.

Initial data exploration revealed the presence of a few outlying values (mostly due to higher prevalence of Diptera and Hymenoptera in certain traps), which were removed from the dataset. Both abundance and species richness of external elements were lower in caves in respect to MSS (Next: t = –4.37, p < 0.001; Rext: z = –4.93, p < 0.001). Conversely, the abundance of subterranean elements was higher in cave (Nad: t = 5.39, p < 0.001; Table 1).

The best model structures resulting from model selection concerning the analysis of MSS and caves (equations 2 and 3) are reported in Table 1. The analysis of the MSS data revealed contrasting patterns of species richness and abundance of external and subterranean elements according to sampling depth. The abundance and species richness of external elements was lower at –0.80 m (test relative to “0.40 m”; Next: t = –3.217, p = 0.001; Rext: t = –3.283, p = 0.001), whereas no significant effect was detected in respect to –0.60 m (Fig. 5a, b). The abundance of specialized elements was higher at –0.60 m (test relative to “0.40m”; Nad: z = 2.481, p = 0.013) and –0.80 m (test relative to “0.40m”; Nad: z = 4.333, p < 0.001; Fig. 5c). The ageing of the bait negatively influenced the abundance and richness of external elements, with higher values with fresh bait (test relative to “Exhausted”; Next: t = 4.489, p < 0.001; Rext: t = 3.005, p = 0.002). No effect was detected in respect to the abundance of adapted elements (Nad: z = 4.333, p = 0.507).

Figure 5. 

Boxplots showing the results of the regression analysis of the MSS (a–c) and the cave (d–f) data. Outlying values are not shown. Significance codes: < 0.001 ***; < 0.005 **; < 0.05 *.

Abundance and richness of external elements showed a significant non-linear U-shaped trend in respect to the sampling series (Next: F = 9.02, p < 0.001; Rext: F = 10.68, p < 0.001), with higher values in summer and early autumn, followed by a drastic decline in winter and an uprising in spring (Fig. 6 a, b). There was no significant relationship between the sampling series and the abundance of adapted elements.

Figure 6. 

Predicted values (black line) and 95% confidence intervals (grey surface) of the effect of the sampling series (Serie_i) on the abundance of external elements in the MSS (a), on the species richness of external elements in the MSS (b) and on the abundance of external elements in the cave at subjacency of 0–20m (c) derived from GAMM analyses. Only fixed effects are shown.

Concerning the cave habitat, we detected an higher abundance of external and adapted elements in the Creusa cave (test relative to “Borna”; Next: t = 3.693, p < 0.001; Nad: t = 2.965, p = 0.003), whereas the richness of external elements was not significantly influenced by this parameter (Fig. 5d, f). We detected a significant interaction between the sampling series and subjacency, with non-linear effects of the series in respect the abundance of external elements at 0–20m subjacency (Next at “0–20m”: F = 5.544, p < 0.001). In the proximity of the surface, higher abundances were predicted during summer, autumn and spring, while a drastic decline was observed in winter (Fig. 6c). There was no significant interaction between the sampling serie and the other subjacency levels.

We observed a decrease in the richness of external elements at increasing distance from the cave entrance (Rext; z = –2.575, p = 0.010). Moreover, the richness of external elements was lower at 60–80m subjacency (test relative to “0–20 m”; Rext: z = -2.298, p = 0.021). No significant effect was detected in respect to the 20–40m subjacency (Fig. 5e). There was also a significant increase in the abundance of subterranean elements at higher subjacency in respect to 0–20m (Nad at “20–40m”: t = 4.525, p < 0.001; Nad at “60–80m”: t = 2.245, p = 0.024; Fig. 5f). There was no significant relationship between the distance from the entrance and the abundance of adapted elements.