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).