Responses of soil microarthropod taxon (Hexapoda: Protura) to natural disturbances and management practices in forest-dominated subalpine lake catchment areas

Disturbances are intrinsic drivers of structure and function in ecosystems, hence predicting their effects in forest ecosystems is essential for forest conservation and/or management practices. Yet, knowledge regarding belowground impacts of disturbance events still remains little understood and can greatly vary by taxonomic and functional identity, disturbance type and local environmental conditions. To address this gap in knowledge, we conducted a survey of soil-dwelling Protura, across forests subjected to different disturbance regimes (i.e. windstorms, insect pest outbreaks and clear-cut logging). We expected that the soil proturan assemblages would differ among disturbance regimes. We also hypothesized that these differences would be driven primarily by variation in soil physicochemical properties thus the impacts of forest disturbances would be indirect and related to changes in food resources. To verify that sampling included two geographically distant subalpine glacial lake catchments that differed in underlying geology, each having four different types of forest disturbance, i.e. control, bark beetle outbreak (BB), windthrow + BB (wind + BB) and clear-cut. As expected, forest disturbance had negative effects on proturan diversity and abundance, with multiple disturbances having the greatest impacts. However, differences in edaphic factors constituted a stronger driver of variability in distribution and abundance of proturans assemblages. These results imply that soil biogeochemistry and resource availability can have much stronger effects on proturan assemblages than forest disturbances.

The CL catchment covers an area of 89 ha (including the lake area of 10.7 ha) and the PL catchment covers 67 ha (with the lake area of 7.6 ha) 33 . The bedrock in the catchments is predominantly composed of mica-schist (muscovite gneiss), with quartzite intrusions in the CL catchment and granites in the PL catchment. Cambisols and Haplic Podzols predominate at both catchments 34,35 . The soil cover in the CL catchment is dominated (ca. 58%) by Cambisols, (0.49 ± 0.20 m deep), whereas in the PL catchment undeveloped, thin, organic-rich soils (O and A horizons, 0.20 ± 0.13 m deep) that are found across ca. 38% of the area are most frequent 34,35 .
Proturan sampling and laboratory analyses. Soil samples were collected twice a year for three consecutive years, i.e., in July and October 2012, June and October 2013, and June and September 2014. During each sampling session, five replicate soil cores were taken randomly with a steel corer. Each sample represented a soil core 3.6 cm in diameter (10 cm 2 in area each), to a maximum depth of 7-12 cm (depending on the soil). The soil microarthropods were subsequently extracted in a modified high-gradient apparatus in the laboratory for seven days. All specimens of proturans were separated from extracted material under a dissecting microscope, fixed in ethanol and mounted in permanent microscope slides. Species identification was based on recent taxonomic keys [36][37][38][39] . All specimens were deposited in the collection of Institute of Systematics and Evolution of Animals PAS (Kraków, Poland). Abundance of proturan assemblages per each site was expressed as the density (D) of total number of specimens per m 2 . Subsequently, a set of assemblage parameters was calculated per each sampling occasion and site including, species richness (S) and Shannon's diversity index (H').

Soil property analyses.
To relate the distribution of proturan species to geochemical characteristics soils samples were collected from each of the sampling sites. These samples were collected once during the first sampling events (in July 2012) from the upper soil horizon (O; Ol+Of; up to 10 cm depth) in six pits (15 × 15 cm, 225 cm 2 ). Prior to the analysis samples were stored at 4 °C in the dark for 3-5 days. Map of the sampling sites in the Bohemian Forest National Park and Protected Landscape Area. CL -Čertovo Lake catchment, PL -Plešné Lake catchment: (a) location of the area studied in Czechia, (b) location of the examined catchments area at the Bohemian Forest, (c) study sites in Čertovo Lake (CL) catchment area, (d) study sites in Plešné Lake (PL) catchment area. Numbers refer to the study sites given in Table 1.
In the lab each soil sample was passed through a 5-mm stainless-steel sieve to remove coarse particles, and was then divided into two parts. One was analyzed for concentrations of water-extractable dissolved organic carbon (DOC), dissolved nitrogen (DN), and dissolved phosphorus (TP H2O ) (water extract 1:10 by weight; field moist soil; 1 hour shaking on a horizontal shaker; filtration through glass fibre filters). Concentrations of DOC and DN were determined with a Formacs TOC/TN analyser (Skalar, Netherlands), and TP H2O was determined by perchloric acid digestion and the molybdate method 40 . The pH was measured in distilled water (pH H2O ) from naturally moist soils (1:10).
The other part of the soil sample was air-dried between two sheets of filter paper for 14-21 days at laboratory temperature, sieved through a stainless steel 2-mm sieve, then pooled to form a composite sample for each disturbance type per catchments and analyzed for exchangeable cations (CEC) and exchangeable acidity (Al 3+ ex and H + ex ). Exchangeable base cations (BC ex = sum of Ca 2+ , Mg 2+ , Na + , K + ) and exchangeable acidity (the sum of Al 3+ ex and H + ex ) were determined at natural soil pH by extracting 2.5 g of the soil with 50 ml of 1 M NH 4 Cl and 1 M KCl, respectively. Base cation concentrations were measured via atomic absorption spectrometry (Varian, Australia), and Al 3+ ex and H + ex were determined by titration (phenolphthalein, 0.1 M NaOH) according to Thomas 41 . Effective cation exchange capacity (CEC) was calculated as the sum of BC ex , Al 3+ ex and H + ex . Base saturation (BS) was calculated as the percentage of BC ex in CEC. The soil contents of DOC, DN and TP H2O were based on three replicates, while CEC, BS, and Al 3+ ex and H + ex were based on the analysis of one mixed sample. Soil moisture was measured during each sampling session weighting the soil cores collected for proturan analysis before and after extraction of animals. Gravimetrically determined soil moisture was expressed as the percentage of actual water content per soil dry weight. Soil temperature was measured with dataloggers (TidbiT v2, UTBI, Onset Computer Corporation, USA) throughout the entire research period (2012-2014). In order to calculate the mean soil temperature (T soil °C) for individual study sites, the mean daily temperatures obtained for all sampling days were used.
To evaluate the relationship between variability of proturan assemblages and availability of their food resources the following approaches were undertaken: (1) dissolved organic carbon (DOC) concentration represented the main source of energy rich C substrates in soil for microbial decomposition 42 and was considered as a proxy of microbial biomass and activity 42 ; (2) dissolved nitrogen (DN) and total dissolved phosphorous (TP H2O ) were considered as microbial nutrient limitation and used as a proxy of root colonization by mycorrhizal fungi 43 . Data analysis. The Kruskal-Wallis one-way ANOVA was used to compare differences in soil properties among sampling sites, representing different disturbance regimes. Furthermore, differences in soil properties between catchments were examined using the two-tailed Mann-Whitney U-test. A principal component analysis (PCA) was then performed to test whether variations in centered and standardized (z-scores) of the soil characteristics were able to discriminate the forest stands. In order to determine relevance of differences between the bedrock type, the PCA scores of weights of soil characteristics for the first and second axis were used as a synthetic dependent variable in the Mann-Whitney U-test.
Linear Mixed Model (LMM, Gaussian error, restricted maximum likelihood) was performed to analyze effect of disturbances regimes on value of Shannon diversity index (H'). Normal distribution of residuals was assumed. Forest disturbance regimes were used as a fixed factor and the nested time was used as a random factor (1 | year/ season). The LMM model was run using the lmer function in the nlme R-package 44 and Imer Test R-package 45 . Count data, such as Protura species richness (S) and density (D; expressed as thousands of individuals per m 2 ) were analyzed using Generalized Linear Mixed Models. GLMM with Poisson error and maximum likelihood was applied to analyze the effects of disturbance regime on Protura species richness (S). Density of Protura (D) was analyzed using GLMM model with negative binomial error and maximum likelihood estimator. Model with

Catchment
Plešné Lake (PL) Čertovo Lake (CL) www.nature.com/scientificreports www.nature.com/scientificreports/ negative binomial error was chosen because proturan density data were overdispersed. GLMMs were performed using glmer function of the Ime4 package 45 . LMM and GLMM models were fit with 3 groups of disturbances contrasted to undisturbed (control) sites. Manual model selection on the three possible models (using as random factor: season, year or the nested year/season) was carried out by the bias-corrected Akaike information criterion (AICc).
The effect of lake catchment, time and forest site disturbances on the variation in Protura assemblages composition was assessed with Canonical Correspondence Analysis (CCA). CCA was appropriate because the response data were compositional and had gradient length longer than 4.0 SD as determined by Detrended Correspondence Analysis (DCA) 46 . CCA was constrained by factors such as bedrock type (granite, gneiss), time (year, season), and forest stand disturbance regimes (control, bark beetle outbreak, windthrow and clear-cut). We also used partial CCA (pCCA) with time as covariate to corroborate differences in the composition of proturan assemblages between catchments and forest disturbance regimes. Decomposition of variation based on partial CCA ordination 47,48 was performed to explain the unique effect (conditional) of the lake catchment (C), time (T) and forest disturbance regimes (FS) and the joint overlap of the investigated factors upon proturan assemblages. The adjusted variation using the number of degrees of freedom was applied. The significance of the models was estimated by the unconstrained Monte Carlo permutation test. We performed the multivariate analysis using the data matrix species × forest stands with the row data (number of individuals summed from five soil cores) collected at each of the forest stands during the same sampling session and then standardized to individuals per square meter. In the LMM, GLMM and multivariate analysis, time was treated as repeated measures for adjusting p-values as a function of correlation in the data. Before the multivariate analyses, proturan data were log(x + 1)-transformed to satisfy assumptions of normality and uniformity of variance. The level of significance in all analyses was at α = 0.05. Calculations were made using Statistica 10.0, Canoco 5.0 47 software packages and R v. 3.4.4 49 .

Results
Edaphic factors. Descriptive statistics of the recorded edaphic properties for forest stands of the Plešné (PL) and Čertovo (CL) lake catchments are presented in Table 2 The results of PCA ordination corroborated the separation of soils on granite bedrock (PL sites) from those on mica-schist gneiss (CL sites), with significantly higher loadings for the first PC axis (two-tailed Mann-Whitney U-test for sum rang of soil properties: p = 0.0002; n = 8 in all cases), (Fig. 2). The PCA Axis 1 accounted for 61.8% of the total variance in soil properties and represented the negative correlations with exchangeable acidity (Al 3+ and H + ) and positive with base saturation (BS). The PCA Axis 2 appeared to be related positively to cation exchange capacity (CEC) and soil moisture, and accounted for 18.7% of the total variance (Supplementary  Table S1). In the PCA plot (Fig. 2) the granite sites (PL) are on the right and the sites on gneiss (CL) on the left. Nonetheless, the PL-3 and CL-3 sites with bark beetle outbreak and windthrow (wind + BB) were located at a short distance on the diagram. These sites are characterized by a high content of acidic cations (Al 3+ , H + ) and lower values of base saturation (BS).   Table S2). The majority of the recorded proturans belonged to the genus Eosentomon Berlese, 1909 (10 species), while the genus Acerentomon Silvestri, 2007 was represented by a single species, Acerentomon tuxeni Nosek, 1961. The most widespread species was Eosentomon gramineum Szeptycki, 1986.
The generalized linear mixed model (GLMM) which contrasted disturbance regimes with undisturbed sites (control) demonstrated that combined disturbance regimes (wind + BB), besides disturbances by bark beetle outbreak (BB) and post-disturbance forest management practices (clear-cut), had significant and the strongest negative effect on Protura density (D). Furthermore, the evaluation of disturbance regimes on Protura species richness (S) using the generalized linear mixed model (GLMM), and on Shannon's diversity index (H') using the linear mixed model (LMM), also demonstrated negative impact of the combined disturbance regimes, such as bark beetle + windthrow. There was no statistically significant relationship between the value of Protura species richness (S) and disturbances by bark beetle outbreak (BB) or post-disturbance forest management practices (clear-cut) alone. However, in the case of H' index, the lack of relationship was only found in response to bark beetle outbreak (BB) ( Table 3). The variation in mean values of the basics parameters of Protura assemblages, such as density (D), species richness (S) and Shannon's index of diversity (H') among disturbance regimes in both catchments can be found as Supplementary Table S2.
The total CCA analysis demonstrated that all of the investigated factors (i.e., lake catchment, time and forest stands) accounted for 43.1% of the variation (adjusted explained variation 12.8%) in the species dataset (λ-trace = 2.63, F-ratio = 1.42, p-value = 0.01). Variation decomposition (Table 4) demonstrated that lake catchment Figure 2. The PCA ordination biplot (PC 1 and PC 2) of the forest stands exposed to different disturbances, with subalpine lake catchments on various bedrock type) as the supplementary variable. The model calculated with row data of soil characteristics, interspecies correlation scaling, species score divided by SD and centring by species, not standardized by sample. Site abbreviations and numbering are given in Supplementary Table S1   www.nature.com/scientificreports www.nature.com/scientificreports/ had the greatest unique effects within the model and accounted for 7.2% of the adjusted variation, whereas disturbance regimes and time were not significant and explained only a relatively small part of the variation (1.2% and 3.1%, respectively). Moreover, the joint overlap effect of all three factors was negative. The results of partial CCA (pCCA) with time as a covariate further corroborated differences in the composition of proturan assemblages between catchments and forest stands with different disturbance regimes showing significant species-environment relationships (λ-trace = 1.39, F-ratio = 1.50, p-value = 0.03). For instance, Eosentomon silesiacum Szeptycki, 1985 and E. occidentale Szeptycki, 1985 were associated with catchment sites underlying gneiss and E. semiarmatum Denis, 1927 and Acerentomon tuxeni Nosek, 1961 with sites underlying granite, whereas E. germanicum Prell, 1912 was associated with sites exposed to the effects of multiple disturbances due to both bark beetle outbreak and windthrow (wind + BB) (Fig. 3).

Discussion
One of the main goals of our study was to assess responses of proturan assemblages from montane forests to natural catastrophic disturbances and to post-disturbance forest management practices. To achieve that goal, our sampling design covered forests with four different type of disturbance, distributed across two subalpine glacial lake catchments representing variability in local geology. As expected, forest disturbance had negative effects on proturan diversity and abundance. However, results of the LMM and GLMM models (Table 3) demonstrated that multiple forest disturbance, i.e, bark beetle outbreak (BB) + windthrow significantly impacted proturan abundance, and value of Shannon's diversity index and species richness. On the other hand differences in diversity index and richness were between undisturbed forests and sites affected only by bark beetle outbreak (BB) relatively smaller or statistically insignificant.
The apparent differences in proturan abundance, richness and Shannon's diversity index related to a combination of natural disturbances (i.e., wind+BB) suggest that multiple disturbances can play an important role in shaping the response of proturan assemblages. This finding conforms to the concept that the effects of multiple disturbances can drive non-random and directional changes in the structure and diversity of belowground communities by provoking the decline of species, which is often explained by filtering effects of environmental changes on vulnerable traits 50 . Unfortunately it has been impossible to include in our sampling design sites with only windthrow disturbance, to discard the possibility that the combined BB + windthrow negative effects on proturans were due only or mainly to the latter factor. The results of the pCCA analysis (with time as a covariate) demonstrated that Eosentomon germanicum was the only species related primarily to sites on gneiss with a combined effects of natural disturbances (i.e., wind+BB). It is interesting to note, that Eosentomon germanicum is characterised by a wide ecological plasticity, and is broadly distributed in different habitats throughout Europe and Northern Africa 36 . This finding is in agreement with Coyle et al. 17 , who indicated that soil faunal responses to environmental disturbances are often species-specific and can be affected by complex environmental interactions. Furthermore, the withstand of forest disturbances by a widespread-generalist species of Protura shows that combined effect of natural forests disturbances can shift the specialist -generalist balance 51 . This finding also correspond well with Platt and Connell 52 concept suggesting that directional replacement of species following natural disturbances favors "functional homogenization" of biodiversity 53 .
An important question pertaining to natural disturbances such as insect outbreaks and windthrow regards management intervention practices and their importance in mediating the impacts and restoration of forest ecosystems. Although the discussion relates primarily to the restoration of tree growth, understanding the effects of different forestry practices on belowground components, which mediate the supply of ecosystem services for aboveground components 54 , is also crucial in order to devise informed strategies of ecosystem management and conservation 17 . According to some studies 50 , forest management practices have no significant impact on soil-decomposer communities. Our study, however, highlights significant impact of post-disturbance clear-cut practices on diversity measures, richness and abundance of proturan assemblages. This is consistent with findings of Fischer et al. 1 , who demonstrated that management practices often have stronger impacts on forest ecosystems than the disturbances themselves. Post-disturbance clear-cuts have pronounced impacts on forest floor heterogeneity and considerably decrease natural buffering capacity of the soil by lowering organic matter contents 55 . Protura occur mostly in the tree rhizosphere; hence, the removal of trees from windthrow stands and subsequent environmental changes induced by soil erosion, loss of soil nutrients, modification of the soil moisture regime, increased sun exposure and changes in the soil microbial communities 56 can all explain the negative implications for abundance and diversity of proturans in the soil.
Nonetheless, it is important to indicate that results of our study demonstrated that local differences in soil characteristics of subalpine glacial lake catchments, related to underlying bedrock, constituted more important driver of variability in distribution and abundance of proturan assemblages than forest disturbances. For instance, the variance partitioning showed that the catchment underlying various bedrock had the greatest unique effect on the variation in proturan assemblage composition and explained a more than twice larger part of the variation than the type of disturbance (Table 4). These high values support our hypothesis that differences in the physicochemical properties of soils, derived from different bedrock in the subalpine glacial lake catchments, will affect the species composition of proturan assemblages. This finding is consistent with Hahm et al. 29 who pointed that bedrock composition can act as major bottom-up regulators of the montane ecosystems. But forest disturbances and tree dieback can also significantly influence the soil chemical properties and organic matter content across the study sites, and thus diversity and distribution patterns of the investigated proturan assembladges. The significant joint effect of both catchments and forest stands and negative shared variation implies that these variables overlap and that they have a additive effect on variations within the investigated proturan communities.
Taking into account enhanced values of total phosphorous and lowered of nitrogen content and moisture in soils developed on granite and their negative correlations with exchangeable acidity and positive with base saturation (BS) (i.e., Plešné Lake catchment, Table 2, Fig. 2, Supplementary Table S1) it is reasonable to expect that these variables provide the most important direct and/or indirect predictors explaining variation in proturan species distributions. Though, similarly to other soil microarthropods, differences in proturan assemblages are most likely associated also with differences in soil moisture and pH 22 , changes in the content of nutrients between these unsaturated, shallow and highly permeable catchment soils are linked via the soil microbial communities 57,58 . The same set of variables were usually linked to distribution, abundance and diversity in numerous studies of other soil microfauna e.g. 59 indicating that the relationships found in our study between proturans and soil geochemistry could be generalized.
Protura, similarly to other soil microarthropods, participate in transformation and release of nutrients from organic matter, thus by affecting decomposition processes both directly and indirectly through comminution of litter, feeding on microorganisms and dispersal of microbial propagules 60 . Feeding on soil microbiota proturans comprise an important middle link in the soil food webs. According to Šantrůčková et al. 30 , the higher ability of granite to release phosphorous into the soil supports higher microbiologically mediated P-fluxes. This essential nutrient can directly mediate root colonization by mycorrhizal fungi 43 , which provide a major food resource for proturan 23 . Thus the relationship between proturans and differences in soil geochemistry between the catchments can be related to the indirect influence of bedrock on fungal community composition. However, it is important to highlight that geochemical factors can significantly influence distribution of soil invertebrates by affecting the availability of their food resources 56 . Although the lack of appropriate data on proturan feeding ecology and trophic relationships prevents us from informed connections between species distributions and their potential food resources, the preference of at least some species for ectotrophic mycorrhizal fungi was shown by some authors e.g. 23 . This indicates that proturans can be highly sensitive to soil geochemistry and environmental changes that impact mycorrhizae. For instance, an extraordinary high abundance of Protura was reported from a windfall area of a spruce forest three years after the storm 61 . These high densities were related to the amount and/ or vitality of ectotrophic mycorrhiza associated with the fine roots of spruce seedlings in unsalvaged gaps. On the other hand, significant decline of proturans was recorded in a tree girdling experiment that constrained development of ectotrophic mycorrhizae fungi 62 . Considering the apparent importance of food resources, the high differences in proturan species abundances and their local distributions imply the existence of species-specific feeding preferences. In fact the dietary specialization can be more important driver of proturan distributions than the soil geochemistry and/or forest disturbances, which may act mostly as indirect drivers through their effects on ectomycorrhizal fungi and other decomposer rhizomorph producers 63 . Thus better knowledge of proturan feeding preferences might be the key to our understanding the importance of drivers underlying changes in distribution patterns of their assemblages. Focusing on this little known soil microarthropod taxon, our work provides a crucial baseline for further investigations on proturan diversity and understanding the impacts of environmental changes on distribution of their assemblages.

conclusions
To improve understanding of the impacts of forest disturbances on the diversity and distribution of soil arthropods, we investigated proturan assemblages in the glacial lake catchments of the Šumava Mountains in spruce forests subjected to different types of disturbances, including natural catastrophic events such as bark beetle outbreak and/or windthrow, as well as following forest management (clear-cut). As expected, forest disturbance had detrimental effects on proturan diversity and abundance, with multiple disturbances and clear-cuts having the greatest impacts. This result indicates that current management practices tend to have greater impact on the ecosystem than the natural catastrophic events themselves. Our findings also suggest that the composition of proturan assemblages was primarily driven by edaphic factors directly related to the type of bedrock. It is likely, however, that the observed associations between species distribution and soil geochemistry are a consequence of indirect relationships as geochemical factors can significantly influence distribution of secondary consumers by affecting the availability of their food resources.