Impact of forest disturbance on microarthropod communities depends on underlying ecological gradients and species traits

Sampling design and site selection

Sampling was carried out within the Eastern Italian Alps, between Trento and Vicenza provinces. Our study area was severely hit by Vaia windstorm in October 2018 causing large windfall areas (Chirici et al., 2019). Windfalls mainly occurred on spruce forests and with a patchy distribution. Within our study area, windfalls were firstly mapped with high-resolution satellite images. A Digital Elevation Model (DEM, 25 m resolution) was retrieved from https://land.copernicus.eu/ and was used to compute topography-related variables, i.e., elevation, slope (inclination degree), aspect (radiant distance from South). Data from 61 local climate stations were used to obtain precipitation time-series of the previous 10 years (retrieved from https://www.meteotrentino.it and https://www.arpa.veneto.it/). Mean annual precipitation values were used to compute a continuous map using kriging interpolation methods in SAGA version 7.8.2 (Conrad et al., 2015). These data were used in the sampling design and as explanatory variables in the models. Field inspections were carried out in candidate sites to assess operator accessibility, salvage logging operations and the presence of undisturbed spruce forest in their proximity with similar topographic and stand conditions. Only sites with a clear predominance of Norway spruce (Picea abies (L.) H. Karst.) forest and already logged windfalls (i.e., only stumps) were considered eligible. In total, we selected 29 sites, from four geographic zones, according to independent gradients of elevation (from 1,100 to 1,950 m a.s.l.), precipitation (from 1,077 to 1,677 mm), slope (from 3 to 34 degree), and aspect (Fig. 1). The maximum distance between sites is about 50 km. Sites were arranged in four different geographic areas consisting of four different mountainous massifs. Vulcanic bedrocks are present in the North-West, carbonate bedrock is the predominant substrate in the South and in the East.

Then, we applied a nested pair-based sampling design. Once in the field, for each selected site, we chose two sampling points including a point in the windfall area and a point in the undisturbed forest nearby. A minimum distance of about 30 m from boundaries and ecotones was considered to avoid edge effects. The average distance between paired sampling points (within pair) was 114 m (min: 50 m; max: 265 m). Geographic positions of each sampling point were recorded with a GPS receiver. Given the precise location of the two sampling points (within each site), point-specific elevation, slope, and aspect were computed to check if there were difference between the two sampling points. In each site, the two sampling points (undisturbed forest and windfall) did not differ for elevation, slope, and aspect.

Within each sampling point, we selected three pseudo-replicated small plots of 20 × 20 cm. Pseudoreplicates were chosen few meters from the center of the sampling point. In windfall habitat, bare soil due to uprooting and machinery trails was avoided, and samples were taken from undamaged soil. For each pseudo-replicate, surface moss vegetation and soil until 5 cm depth were collected using a shovel and soil corer (Hishi, Kawakami & Katayama, 2022). Percentage of moss cover was recorded by visual inspection estimates. For each sampling point, the three pseudo-replicates were mixed and pooled in a unique sample before extraction of invertebrates. Tools used for soil collection were cleaned before changing sampling point using a bleach solution. We performed two rounds of sampling during 2020, in June and in September, which is almost 2 years after the storm. For each sampling point, a sample of sieved soil (c. 10 g) without roots and moss, was preserved from the last round for organic matter, hereafter OM, analysis. OM content was analyzed by dry combustion using a CN analyzer (Primacs SNC100, Skalar, Breda, Netherlands) on dried, sieved and ground (<2 mm) samples. In detail, organic carbon released up to 400 °C (TOC 400) and residual oxidizable carbon (ROC) released up to 600 °C were analyzed in sequence using a temperature gradient from 150 °C to 600 °C after dry combustion and in presence of oxygen. OM content was then calculated as the sum of these two forms. Detailed information on sampling design is given in Supplemental Materials, Appendix A. For those locations belonging to protected areas, we obtained a sampling permission from the authorities.

Sample preparation and DNA sequencing

In the laboratory we extracted arthropods from about 4 L of soil for each sampling site using a Berlese apparatus maintained active for 3 weeks. Animals were collected in propylene glycol and then were stored in absolute ethanol. For each sample, using decontaminated forceps and small brush, we retained only those invertebrates that were morphologically recognized as belonging to soil mites (orders Sarcoptiformes, Trombidiformes, and Mesostigmata) or springtails (class Collembola). We selected only these groups because they show similar average sizes and our sampling design is appropriate for retrieving representative communities of them. Other groups of larger arthropods, such as Carabidae, Chilopoda or Araneae, were removed since they would have been occurred randomly within the samples. Abundances of mites and springtails were counted under a dissecting microscope. All animals from each site, collected during the first and second round of sampling, were pooled together and stored in absolute ethanol until DNA extraction. Before DNA extraction, samples were dried in a Eppendorf Vacufuge centrifuge with a Qiagen vacuum pump (KNF Laboport Vacuum pump 1.0 bar) to remove ethanol. During counting and dehydration steps, some microcentrifuge tubes with 1 mL of ethanol were left open to check for potential contamination (negative controls). DNA extraction was performed using DNeasy^® PowerMax^® Soil extraction kit (Qiagen, Hilden, Germany) to ensure an optimum DNA extraction even in the presence of soil particles. DNA extractions of negative controls were routinely processed to assess contaminants from the lab procedure. For metabarcoding library preparation, we amplified a 313-base pair (bp) fragment of the mitochondrial DNA Cytochrome Oxidase c subunit I (COI) gene. PCR amplification was done using the mlCOIintF forward primer 5′-GGWACWGGWTGAACWGTWTAYCCYCC-3′ and jgHCO2198 reverse primer 5′-TAIACYTCIGGRTGICCRAARAAYCA-3′ (Leray et al., 2013). Prior to PCR amplification, a PCR condition setting up was performed in order to optimize amplicon yield. 5 µL of DNA was added to 45 µL of PCR mix, including 23.75 µL of water, 10 µL of FlexiBuffer, 6 µL of MgCl solution (25 mM), 2 µL of each primer (10 ρmol/μL), 1 µL of dNTPs (10 mM), and 0.25 µL of Promega GoTaq Hot Start G2 (Promega, Madison, WI, USA). PCRs were run under the following conditions: one pre-incubation step of 95 °C for 2 min, followed by 45 cycles of 95 °C for 30 s, 47 °C for 60 s and 72 °C for 60 s, with a final extension at 72 °C for 5 min. All DNA amplifications were performed in a Veriti 96-Well Fast Thermal Cycler (AB Applied Biosystems).

For each sample three PCR replicates were performed, and the resulting amplicons were pooled before the purification step. Purification was carried out using QIAquick PCR purification kit (Qiagen, Hilden, Germany). Amplifications of negative controls belonging to the same PCR run were pooled. Library indexing and amplicon sequencing (300 bp) were performed at IGA (IGA Technology Services, Jacopo Linussio 51, Udine, Italy), using NexteraXT Index Kit (FC-131-1001/FC-131-1002) and an Illumina NOVASEQ6000 (300-bp paired-end mode). Please refer to the Supplemental Materials, Appendix A for the scheme of the overall sampling and sample processing protocol.

Bioinformatic pipeline

Demultiplexed paired-end raw reads were used to retrieve ASVs (Amplicon Sequence Variants) using a custom pipeline. First, primers were cut using cutadapt (Martin, 2011) and reads without primers were discarded. Second, we used a denoising approach to retrieve Amplicon Sequence Variants (ASVs). Reads were filtered, denoised, merged, and chimeras were removed using dada2 version 1.16 package (Callahan et al., 2016).

We performed a first taxonomic assignment of ASVs using Bayesian classifier in qiime2 (https://qiime2.org) and blastn on a custom database. Unassigned, low confident (<0.97 for Bayesian confidence or <0.97 identity for blastn) or ambiguous ASVs were mapped against BOLD database (www.boldsystems.org/) to confirm their identity (January 2022). ASVs that were successfully assigned to our target groups (i.e., class Collembola and orders Sarcoptiformes, Mesostigmata, and Trombidiformes) were retained, whereas those not belonging to these groups were discarded. A further curation step was performed on sequence alignment to remove sequences with gaps or stop codons, since the amplicon region was expected to be entirely coding. AliView software (Larsson, 2014) was used for aligning and edit the sequences separately for each group.

Approaches combining denoising and clustering generally produce reliable species-like entities, which are relatively similar to those obtained with non-molecular approaches (Antich et al., 2021). However, different taxa might show different barcoding gap size, by varying intra- and inter-specific distances. Commonly used clustering algorithms (such as vsearch and swarm) depend on a priori settings and they were not optimized for a relaxed taxonomic clustering. Thus, we used an approach from DNA taxonomy, namely ASAP (Puillandre, Brouillet & Achaz, 2021), developed to identify the most appropriate distance threshold for species-level clustering from alignments of DNA sequences. ASAP clustering was performed separately for the ASVs of each group, soil mites and springtails. In order to increase the reliability of the approach, we downloaded all available species-level COI sequences of mites and springtails from BOLD database, aligned and trimmed them to our ASV datasets, and we used the large alignments (i.e., reference sequences together with ASVs) in ASAP (https://bioinfo.mnhn.fr/abi/public/asap/). Following ASAP, ASVs were merged into Operative Taxonomic Units (OTUs) based on the barcode gap distance for each group. This aggregation step decreases intra-specific metabarcoding entities, and increases comparability with data retrieved by traditional identification. Each OTU was identified at the best taxonomic resolution according to ASV taxonomic assignment, and eventually updated according to ASAP clustering results. The final OTU table was also reduced by removing OTUs with read numbers accounting for less than 0.1% in each sample, in order to further diminish the risk of false positives in our DNA metabarcoding pipeline. The resulting incidence OTU table (only presence/absence) was used as the base for the following ecological analyses (see available data and Supplemental Materials, Appendix B).

Statistics

First, we tested the effect of Vaia windstorm on the overall abundance of microarthropods obtained from visual counts. We used linear mixed-effect models (LMMs) with habitat type (two levels: windfall and undisturbed forest) as a predictor, in addition to elevation, precipitation, slope, aspect, OM, and their interaction with habitat type as fixed effects, with the identity of each pair of samples nested in geographic zone (four levels) as random effect. Non-significant interaction terms were removed, and models were run again. LMMs were fitted using lme4 package version 1.27.1 (Bates et al., 2015) in R version 4.1.3 (R Core Team, 2022). The response variable was log-transformed to meet model assumptions. We used DHARMa version 0.4.5 (Harting, 2021) and car version 3.0 (Fox & Weisberg, 2019) packages for checking model assumptions and collinearity among predictors. Analysis-of-variance tables were extracted with function Anova() in car. Spatial autocorrelation of residuals was checked using Montecarlo simulation (999 simulations) of Moran I metric. We did not find issues affecting the models (spatial autocorrelation was not significant: soil mites P-value = 0.40; springtails P-value = 0.36).

Second, we investigated the effect of Vaia windstorm on species (OTU) richness. The models we used had the same structure and rationale of the LMMs used for abundance data, separately for each of the two taxonomic groups. Model assumptions, presence of collinearity and spatial autocorrelation were checked as well. We did not find issues affecting the models (spatial autocorrelation was not significant: soil mites P-value = 0.25; springtails P-value = 0.50).

Third, we investigated the effect of wind disturbance on functional guilds, assigning each OTU to trophic niches based on literature. Although most of the OTUs were not identified to species level, feeding preferences of species are mostly the same within each family (Krantz & Walter, 2009; Potapov et al., 2016). Soil mites were divided in predators (including omnivorous) and no predators (including primary and secondary decomposers) (Schneider et al., 2004; Krantz & Walter, 2009; Maraun et al., 2011; Fischer, Meyer & Maraun, 2014; Maaß et al., 2015; Schaefer & Caruso, 2019; Nae et al., 2021). We preferred using such broader categorization since the trophic guilds of many species are still unknown. Springtails were divided following Potapov et al. (2016) in four guilds: euedaphic microorganism consumers (hereafter euedaphic), hemiedaphic microorganism consumer (hereafter hemiedaphic), epigeic animal and microorganism consumers (hereafter animal consumers), epigeic plant and microorganism consumers (hereafter plant consumers). To test the effect of disturbance on trophic guild we used LMMs with the number of OTUs as response variable; trophic guild, habitat, and interaction between trophic guild and habitat type as predictors; site ID nested in pair ID nested in geographic area as random factor. The response variable was log-transformed for soil mite dataset. We used DHARMa and car packages as for the other models; in addition, pairwise comparisons were extracted with emmeans package (Lenth, 2020).

Fourth, we investigate differences in species composition (i.e., beta diversity) based on OTUs occurrence data. To test the effect of ecological predictors on species composition we used Adonis analysis with marginal effects using adonis2() in vegan version 2.7 (Oksanen et al., 2020). As a response variable, we used occurrence-based Sorensen beta dissimilarity index among communities, separately for springtails and soil mites. As predictors, we used habitat type, OM, and the geographic zone and pair ID as blocks. We used 9,999 permutations for P value computation.

Effects of wind-induced disturbance and interactions with underlying ecological gradients

Our results showed an overall and strong decline of microarthropod abundance and richness in windfalls, compared to near undisturbed forest, suggesting that soil perturbations could strongly modify microarthropod communities and, possibly, their ecological functioning (Figs. 2 and 3). Our results are consistent with previous studies based on morphologically identified microarthropod communities showing negative effects of windstorm and salvage logging on springtails (Čuchta, Miklisová & Kováč, 2019; Sławski & Sławska, 2019), soil mites (Lóšková et al., 2013; Wehner et al., 2021b), and proturans (Sterzyńska et al., 2020). Literature evidence showed that post-disturbance communities are mainly constituted by few high-adaptable species, such as Tectocepheus sp. (mites) and some species belonging to Entomobrydae (springtails) (Maraun et al., 2003), or by species with disturbance-adapted strategies, such as parthenogenetic reproduction (Wehner et al., 2021a). On the contrary, other studies found that richness of springtails increased in disturbed plot (Urbanovičová, Kováč & Miklisová, 2010; Urbanovičová, Miklisová & Kováč, 2014), and others found divergent responses in abundance and species richness (Čuchta, Miklisová & Kováč, 2012). On the one hand, these differences might depend on the time span of the studies since the response may change depending on time from disturbance (Čuchta, Miklisová & Kováč, 2019). For instance, after intense drought mites were slower than springtails in recovery processes (Lindberg & Bengtsson, 2005). On the other hand, different responses might also depend on interactions with site-specific conditions. Indeed, here we showed how underlying ecological gradients (i.e., precipitation) beyond the main effect of disturbance, can change the response of the communities, especially in complex environments such as alpine mountain forests.

Indeed, we showed that interaction between canopy openness (i.e., wind-induced disturbance) and precipitation influenced the response of soil communities to disturbance. We found higher effect-size of disturbance on mite abundance at low mean annual precipitation (Fig. 2). In other words, drier conditions (i.e., low annual precipitations) showed higher loss of mite individuals after forest canopy removal. Precipitation gradient can strongly influence soil biota, especially for forest ecosystems where aridity might have detrimental effects on soil-living communities (Blankinship, Niklaus & Hungate, 2011). Indeed, moisture is known to be an important factor for microarthropods (Hopkin, 1997; Lindberg & Bengtsson, 2005), and habitat alterations resulting in drier conditions might exacerbate negative effects on microarthropod communities. Despite analyzing soil moisture might be more informative for soil arthropods than large-scale rainfall, moisture is expected to change with forest disturbance and cannot be separated from the disturbance gradient. On the contrary, two close sampling points are expected to have the same average rainfall, and rainfall effect can be tested independently. We suggest that disturbance-driven soil impacts might exacerbate under climate change in the future, especially for forest ecosystems. Despite soil disturbance after windstorms could be intensified by topography conditions, such as steep slopes, due to soil instability (Mitchell, 2013), we did not find interactions with other local factors. Probably our large-scale approach cannot effectively catch all the effects at the local scale.

Simple and interacting effects of wind-induced disturbance on trophic guilds

Microarthropods could exhibit a large differentiation in trophic niche and feeding preferences, even within taxonomic groups, such as oribatid mites (Schneider et al., 2004) and springtails (Potapov et al., 2016). However, after disturbance, ecosystems could exhibit trophic web declining with collapsing of certain trophic niches. For instance, oribatid communities in plantations showed missing niches in food webs compared to rainforest habitat (Krause et al., 2021), and functional guilds of microarthropods are differently affected by soil warming (Thakur et al., 2023). Here, we assessed the effect of soil disturbance on trophic guilds to test whether forest disturbances can influence functional traits in the disturbed soils. We found that functional guilds of soil mites and springtails are differently affected by windstorm. In particular, we observed that the richness of soil mites strongly declined regardless of trophic guilds. Previous studies found that the proportion of detritivore oribatid mites might depend on forest management (Farská, Prejzková & Rusek, 2014). Many reasons could be advocate for explaining the lack of this effect in our data. First, trophic guild assignment was based on taxonomic identifications since we could not directly measure it. Since taxa might be assigned with different taxonomic resolution, those with low taxonomic resolution cannot be specifically addressed to precise trophic niche, but we used general trophic niche (mostly at family level). Second, we used a coarse trophic partitioning, in order to include taxa with wide range of trophic preferences (generalists or poorly studied taxa). Moreover, previous studies showed that the same species can slightly change its feeding preferences under certain environmental conditions (Melguizo-Ruiz et al., 2017; Maraun et al., 2020). Specific approach, such as isotopic analyses might provide further insights on the effect of soil food webs after wind disturbances in temperate forests (Maraun et al., 2011).

On the contrary, we found significant interactions between habitat type and trophic guild in springtails. Epigeic plant consumers were the most affected guild in windfalls. This guild-specific response might be due to the dramatic decrease of moss cover in windfalls (Supplemental Materials, Appendix C), as moss represents the main feeding substrate. Indeed, higher insolation due to canopy removal changes micro-conditions and salvage logging might largely contribute by exacerbating the desiccation process (Waldron et al., 2014). Our results agreed with previous studies emphasizing that habitat specialist species of springtails, such as bryophilous species, are the most sensitive ones to environmental changes after forest clearings (Urbanovičová, Miklisová & Kováč, 2013). Actually, the combined effect of canopy openings and salvage logging led to huge changes in litter and superficial soil layers (i.e., organic layer), vegetation, and moss cover (Rumbaitis del Rio, 2006), largely affecting those species related to litter microhabitats, in particular specialists. Despite our studies focused on short-term effects, similar long-term studies revealed that springtail communities can remain affected by the disturbance even for decades (Sławski & Sławska, 2019). Hence, we suggest that our findings on springtail perturbations in windfall communities might have long-term outcomes. Here, we showed that mite and springtail communities might have different response to soil disturbance: an overall decrease for soil mites and guild-related losses for springtails.

Effects of wind-induced disturbance and distance on species composition

In our analyses soil mites and springtails responded to different spatial scales. These results agreed with literature, suggesting that dispersal might play a role in post-disturbance communities with group-related responses (Rousseau et al., 2019). Composition of soil mite communities exhibited a high dissimilarity among samples, suggesting that small-scale conditions, such as the presence of microhabitats, are important drivers for species assemblages in forest communities (Fig. 5). Our results showed that habitat type and OM were the most important factors contributing to differences in species composition for mites. Windfall and undisturbed forest hosted completely different communities and differences in mite communities might occur at very small scale, even within the same habitat type. On the contrary, the distance at large spatial scale (i.e., geographic area) did not play an important role in the model. Moreover, our model revealed a relatively high amount of unexplained variance, probably due to the high diversity and dependence of microscale conditions, which are not always available for such large-scale studies. Our findings agreed with the low mobility and strict environmental niche preferences of this group (Lehmitz et al., 2012). Thus, even short distances may actually reflect huge differences in terms of environmental variables, as similar studies found that soil mite communities can differ in species composition within few metres (Dong et al., 2017).

Much more similar communities are expected in springtails since they have higher mobility than mites, especially those species with a well-developed furca (Potapov et al., 2020). Indeed, in our model the main factor shaping differences between communities was geographic zone, while habitat type, aspect, and OM were less important (Fig. 5). Our findings suggested that geographic distance had a higher explanatory power than soil disturbance on differences in species assemblages for springtails. Our result might depend on two main reasons. First, besides the relevance of microenvironmental parameters, within continuous habitat types springtail communities might be shaped by geographic distance due to their limiting dispersal (Arribas et al., 2021). Second, the higher explanatory power of geographic distance than habitat type might still reflect a short-term response to disturbance. Here, we sampled only one forest type (i.e., pure spruce forests) just 2 years after the disturbance, thus it is not surprising that geographic zone (i.e., related to historical patterns) is still the most important factor. On the other hand, small-scale variables, such as habitat type, still play a great importance as predictors for species composition, in agreement with previous studies (Salmon & Ponge, 2012; Sterzyńska & Skłodowski, 2018; Arribas et al., 2021).

Limitations

Our method allowed us to derive comparable results among samples and thus drawing general conclusions on the ecological patterns after massive wind disturbance in forest. Unfortunately, we could only assess the compound effect of windthrow and salvage logging, since no-logged areas were not accessible. Besides this study limitation, salvage logging is extensively carried out in most of the forests in Europe, being the main management strategy after severe windthrows. Moreover, this dataset should not be considered a complete overview of the forest-living communities, because metabarcoding has some limitations. First, different markers might provide slightly different communities. Second, our samples are mainly constituted by specimens, but they included also soil particles, making them between pure tissue samples (i.e. bulk sample) and a soil sample (i.e. environmental sample). Finally, the taxonomic resolution depends on the reference database, which is often affected by the lack of completeness for the studied groups.