2.4. Data analysis
We performed all statistical analyses in R (R Core Team, 2021).
In a first step, we analysed the chemical and physical parameters of collected sediments by means of a Principal Component Analysis (PCA) with the function ‘prcomp’. To detect possible shifts in physical and chemical parameters of collected sediments among transects (HP – MP – LP – C), sections (I – II – III – IV) and caves (Bossea, Caudano, Pertosa-Auletta and Vento), we applied a Permutational Multivariate Analysis of Variance (PERMANOVA, Anderson, 2001) based on Euclidean distances, specifying transects, sections and caves as factors, with the function “adonis” from the vegan package (Oksanen et al., 2019). Statistical significance was tested via 9,999 random permutations.
Before proceeding with statistical models, we converted the two factors of our two-factorial sampling design, namely sections and transects, into ordinal variables in order to effectively represent a gradient of progressive distance from the cave entrance (Sections: I = 1, II = 2, III =3) and from the tourist path (Transects: HP = 1, MP = 2, LP = 3). Although air circulation in caves can be extremely complex (Badino, 2018), we here considered the distance from the cave entrance (Dist_entrance) as a proxy of the dispersal of propagules from the external environment vehiculated by air, i.e. dispersal selective mechanisms. The distance from the tourist path (Dist_path) was intended as a proxy of human disturbance determined by visitors, i.e. disturbance selective mechanisms. Physical and chemical parameters of collected sediments, namely granulometric composition (% silt, sand and clay) and concentrations of organic carbon and nitrogen, were considered as proxy of the environmental selection, i.e. habitat filtering. More in detail, we obtained the dissimilarity based on a Bray-Curtis distance of sediment composition of each sample collected in the tourist part of each cave from the sediment composition of the sample collected in the control area. Therefore, the obtained variable (Substrate) summarises the differences in the substrate conditions in touristic sections compared to the control area.
The effect of these three variables and their interactions was tested on the: (i) ASV-ratios of Bacteria, Archaea and Fungi; ii) dissimilarity measures (i.e. β-replacement and β-richness) of Bacteria, Archaea and Fungi, by means of Generalized Linear Mixed Models (GLMMs). In a preliminary step, we graphically checked for potential outliers in our dependent variables (Zuur et al., 2010). Statistical models were performed with the function ‘glmmTMB’ from the glmmTMB package (version 1.1.2.3, Brooks et al., 2017). We assumed a Gamma error distribution for the set of models performed on the ASV-ratios and a Beta error distribution for the set of models performed on the dissimilarity measures. To account for the spatial dependency of samples within the same show cave, a cave identifier (CaveID) was incorporated as a random factor. The interaction factors were retained only if they significantly contributed to the fit of the model based on the outcomes of the ‘anova’ function. Model validation was performed with the function ‘check_model’ from the performance package (Lüdecke et al., 2021) as well as by visually inspecting the distribution of the residuals (Zuur et al., 2016).