The study area is part of the UNESCO Middle Elbe Biosphere Reserve in northern Germany, and covers ca. 100 km of the Elbe River. The area is characterized by a Central European temperate climate with a mean annual temperature of 9.3 °C and mean annual precipitation of 615 mm (measured at the Lenzen weather station, 53.08°N, 11.48°E). The Elbe floodplain is an anthropogenically altered landscape with a history of dyking, deforestation and agricultural use, with the studied hardwood floodplain forests representing small remnants of the historically extended floodplain forest ecosystem. The arboral vegetation in the study sites is typical for central European hardwood floodplain forests and is characterized by oak (Quercus robur) and elm (Ulmus laevis); typical understory vegetation is hawthorn (Crataegus monogyna). For further information on the study area see Shupe et al. (2021).

In the research programme in which the current study took place, two gradients of ecological habitats are being investigated: a gradient of hydrological regimes and a forest development gradient. For the hydrological regimes, nine forest sites in the study area were investigated (Fig. 1), of which three were located in the active floodplain (i.e. the Elbe floodplain area in front of a dyke flooded during sporadic inundation events, usually after snowmelt in early spring; although none took place during the study period), three in the seepage water zone (i.e. the zone behind the Elbe dykes currently affected by occasional high groundwater levels), and three in the floodplains of tributaries independent from the active Elbe flooding regime. For the forest development gradient, a total of 16 sites of the re-forestation programme were investigated, all of which are located in the active floodplain (Fig. 1). Of these sites four represented dense older forest stands (mean tree age 80–200 years with multi-layered canopy cover and a rich shrub layer), four were in sparse older forest stands (mean tree age 80–200 years with a less developed canopy and shrub layer), and four established in young forest plantations (mean tree age 17–26 years with a strong grass-dominated herbaceous layer). Additionally, four sites were located in nearby agricultural grasslands and served as a control, since these represented the starting point of reforestation. Old dense, old sparse, young forests and grassland are referred to as forest development stages in the following.

Environmental data on vegetation (e.g., cover and species number of vegetation layers, litter cover and thickness, leaf area index) were collected by botanical project partners (see Suppl. material 1. Mean tree age was estimated by measuring annual tree rings of four dominant oak trees per forest site (Shupe et al. 2021). Data on soil properties (e.g., pH, texture, C and N content) of individual forest sites, as well as precipitation data were obtained from soil-science project partners (Vásconez Navas et al. 2023, see Suppl. material 1). The average number of flooding days per year as a 25-year-mean (1990–2016) was estimated for the individual active floodplain forest sites using a digital model developed for the Elbe river by Weber and Hatz (2020) and Weber and Rosenzweig (2020). Density and species number of soil fauna as used as explanatory variables were obtained from 10 MacFadyen soil cores per forest site (Scheunemann et al., in prep.).

Soil temperature and moisture were assessed at the starting day of each Bait Lamina test at 10 measurement points per forest site. Soil temperature at 10 cm depth was measured using a digital soil thermometer, while soil moisture was calculated from soil cores up to 5 cm depth taken adjacent to the soil thermometer. These soil cores had been used for extraction of soil fauna by heat (data not presented, Scheunemann et al., in prep.). The fresh soil cores were weighed, then dried during animal extraction over 10 days at maximally 55 °C and weighed again. The gravimetric water content was calculated as follows:

%_H2O = ((fresh weight – dry weight) / fresh weight) * 100)

Calibration with additional drying at 105 °C for two days revealed essentially no difference to the method above. Mean values of soil temperature and moisture measurements per forest site and season were then used as explanatory variables in the statistical analyses.

Feeding activity of soil mesofauna was assessed by the Bait Lamina test according to the ISO standardisation 18311:2016 (ISO 18311 2018). Test strips (12 × 0.5 × 0.1 cm) are made of PVC and contain a row of 16 drilled holes of 1 mm diameter at a distance of 5 mm from each other. The holes were filled with a mixture of powdered cellulose (70%), wheat bran (27%; < 500 µm particle size) and activated charcoal (3%) as standard bait for soil invertebrates. Filled Bait Lamina strips were obtained from terra protecta GmbH, Berlin, Germany.

We measured soil fauna feeding activity along the forest development gradient in 2019 and in the hydrological regimes in 2020 (see below for detailed sampling dates). The Bait Lamina tests were carried out at the forest sites in three plots at each site during the average periods of highest soil-faunal activity (and avoiding summer dry periods) in late spring and early autumn to assess seasonal differences. Distance between plots within a forest site was a minimum of 10 m to avoid autocorrelation (Irmler 1998), distance between the two closest sites was minimally 300 m (between forest and comparative grassland sites) and maximally 10 km (between active floodplain and tributary sites; Fig. 1). Test strips were placed vertically in the soil such that the bait holes were located at soil depths between 1 and 8 cm to allow evaluation of bait-material disappearance in the respective soil depths. In each plot of every site, 16 test strips were placed in the soil in a 4 × 4 grid with 10 cm distance between single strips. Test strips were exposed to the soil for four weeks. In the forest-development gradient, the experiment started in spring between April 29 and May 3 2019 and finished on May 27/28 2019, and again started in early autumn between September 23 and 26 2019 and finished on October 21/22 2019. In the hydrological regimes, the experiment started in spring on May 11/12 2020 and finished on June 8/9 2020, while in early autumn it started on September 22/23 2020 and finished on October 20/21 2020.

After removal from the field, test strips of each plot were collected together, placed in a separate plastic bag and frozen at -20 °C to stop microbial decomposition of the bait material until further processing. In the lab, attached soil was carefully removed from defrosted strips and the strips visually inspected for empty bait holes by holding against a diapositive slide viewer as a light source. As stipulated by ISO 18311:2016, data was collected as “actively fed on” if light could be observed to shine through a bait hole, and as “not fed on” with bait holes through which no light emitted (even if small feeding traces were found at the sides). Every individual hole in each test strip was inspected in this way, resulting in a data set consisting of binary data representing the feeding activity at each Bait Lamina position (= soil depth). Strips that had been disturbed (chewed on or removed from soil by wild animals) during the four-week field exposure were excluded from further analyses. From all 16 bait laminae of one plot the percentage of bait holes that had been fed on was calculated for every soil depth. Further, the results of the three plots per forest site were averaged to avoid pseudoreplication.

A number of numerical environmental parameters (i.e. soil properties, vegetation and soil fauna parameters; see Suppl. material 1) were available for all forest sites. In a first model-building step, these were analysed for collinearity, which was performed separately for both datasets (forest development gradient 2019 and hydrological regime 2020), since these partially consisted of different forest sites. Collinearity of numerical parameters was investigated by principal component analysis using the prcomp() function in the “vegan” package (Oksanen et al. 2020) after standardization to zero mean and unit variance using the decostand() function. To achieve dimension reduction, we used the resulting site scores of the first three principle components (PC; all that explained more than 10% of variance) instead of the original 42 site parameters as independent variables in the subsequent analyses. Inspection of the scores of individual environmental factors showed that, in the 2020 data set (hydrological regime), PC1 represented the hydrological regime (i.e., number of flooding days per year), PC2 represented soil parameters and PC3 represented vegetation parameters as well as soil pH (Suppl. material 2). In contrast, for the 2019 data set (forest development), the first PC mainly reflected vegetation parameters (e.g. species numbers of shrubs and trees, leaf area index, litter cover), PC2 reflected soil parameters (e.g., total C and N content, pH, soil moisture) and PC3 reflected soil fauna density and further vegetation parameters (Suppl. material 3). Further, we built a correlation matrix using the corrplot() function in the “corrplot” package in R studio (Wei and Simko 2021, not presented). The correlation matrices were inspected for correlation indices, and factors with correlation indices > 0.5 or <-0.5 were considered to be autocorrelated. We calculated correlation indices using Pearson’s as well as Spearman’s correlation coefficient, which both showed similar results. As expected after the PCA, we found correlations within vegetation parameters, soil properties and soil fauna parameters, supporting the reduction of the measured environmental data to the respective PCs.

As categorical variables we used hydrological regime (active floodplain, seepage water zone, tributary) or forest development (old dense, old sparse, young forest, grassland), as well as soil texture estimated a priori by visual inspection in the field (the categories sandy and loamy corresponded well with sand content measured in the lab). In case we found high correlation of a categorical variable with the respective PC (GVIF value > 5 after using the vif() function of the “car” package; Fox and Weisberg 2019), we excluded the respective PC from further analyses and only retained the categorical variable. In this way, in the hydrological regime data set PC1 was replaced by the categorical “hydrology” variable, and in the forest development data set the categorical variable “forest development” replaced PC1.

We then applied linear mixed effects models (lmer() in the “lme4” package (Bates et al. 2015)) to the bait-lamina data sets of the mean percentage of “fed upon” baits at a respective soil depth as dependent variable and site identity and season as random terms. The null models only included soil depth and the two random terms as independent variables, while hydrology (or forest development), soil texture, the respective PCs as well as all possible two-fold interactions of these parameters were added consecutively to the model (forward selection) by hand. Based on Akaike’s Information Criterion (AIC) and anova() functions in the “stats” package (implemented in R), after every step models were compared and only the better model was kept. In case of the forest development data, this procedure resulted in forest development being not included in the model as a single term. Stats::anova() comparing the model with and without forest development as single term (next to soil depth and the two random terms) resulted in no significant difference (p = 0.46) between these models. Forest development did not have a significant effect on feeding activity either (car::Anova showed forest development p = 0.54 for the model containing forest development), therefore we kept the simpler model. In a later step of model selection, however, interactions of forest development with other factors were included in the model as they significantly improved the model. The significance level of variables and interactions was estimated via Type III sums-of-squares Wald Chi-square test using the Anova() function in the “car” package. The final models were as follows:

for hydrological regime:

“Feeding activity” ~ Soil_depth + Hydrology + Texture + PC2 + PC3 + Hydrology : Soil_depth + Hydrology : Texture + Hydrology : Season + Soil_depth : Season + Texture : Season + (1|SiteID) + (1|Season);

and for forest development:

“Feeding activity” ~ Soil_depth + Soil_depth : Forest_development + Soil_depth : PC2 + Soil_depth : Season + Soil_depth : Texture + Forest_development : Texture + (1|SiteID) + (1|Season).

Pairwise comparisons of factors were calculated using the emmeans() function in the “emmeans” package (Lenth 2022). Since soil depth was treated as a numerical instead of factorial variable, pairwise comparisons including soil depth were performed using the lstrends() function in the “emmeans” package.

Furthermore, the average depth of faunal activity was evaluated using the Depth index (DI) (Gongalsky et al. 2004), calculated as:

DI = ((∑ n_i d_i) / N),

with DI being the mean depth of feeding activity, n_i the number of all pierced bait holes at depth i, d_i the respective soil depth [in mm] and N being the total number of pierced bait holes over all soil depths and test strips of one plot. Separate linear mixed effects models (lmer()) for the hydrological regime and forest development data set were built with DI as the response variable using the methods described above.

The final models for effects on the Depth Index were as follows:

for hydrological regime:

Depth Index ~ Hydrology + Hydrology : Season + (1|SiteID) + (1|Season);

and for forest development:

Depth Index ~ Forest development + Texture + Forest_development : Texture + Forest_development : Season + (1|SiteID) + (1|Season).

All statistical analyses were performed using R statistical software, version 4.1.2 (2021-11-01, “Bird Hippie”, R Core Team (2021), with R studio version 1.4.1717. Graphical visualizations were created using the “ggplot2” package (Wickham 2016).

The hydrological regime was related to a number of abiotic site parameters (Table 1; Suppl. material 1: table S1A). Briefly, average number of flooding days per year was 5.2 in the active floodplain, but 0 in the seepage water zone and tributary floodplains. Leaf area index, i.e. density of the forest, increased with distance to the main river, while pH and average sand content of soil decreased. Within the active floodplain, the four forest development stages were rather similar in their environmental factors. However, they differed in leaf area index with old dense and young forests having higher canopy cover (high leaf area index) than old sparse forests (Table 2; Suppl. material 1: table S1B). Comprehensive information on ecological parameters of single sampling sites is given in Suppl. material 1.

The hydrological regime did not directly affect overall feeding activities (p = 0.11), and soil texture, soil parameters (represented by PC2) and vegetation parameters (represented by PC3) alone also did not significantly affect feeding activity. In contrast, soil fauna feeding activity decreased with soil depth in all hydrological regimes (p < 0.001, Table 3), and general feeding activity rates were higher in autumn than in spring (only included as random term in the model). Feeding activities increased from spring to autumn to a higher extent at shallow than at larger soil depths, resulting in a stronger decrease in feeding activity with soil depth in autumn (soil depth x season: p < 0.0001, Fig. 2). However, the decrease in feeding activity with soil depth was also affected by hydrology, as reflected in a stronger depth decrease in tributary sites than in the active floodplain or seepage water zone (interaction soil depth x hydrology: p < 0.0001; pairwise comparisons p = 0.0001 for tributary vs. seepage water zone, p = 0.0004 for tributary vs. active floodplain, but p = 0.88 for active vs. seepage water zone). Activity rates varied according to the season in sites of different soil textures, with sites on sandy soils showing higher activities than those on loamy soils in spring, but not in autumn (interaction soil texture x season: p < 0.0001). The interaction of soil texture and hydrology was not included in the model because all tributary sites were of loamy soil, unbalancing the experimental design in this respect. For comparison of soil fauna feeding activity of single sites see Suppl. material 4.

The Depth Index, i.e. average soil depth of soil fauna feeding activity, was significantly affected by hydrological regime (p = 0.001) with maximum average activity being at larger depths in the seepage water zone compared to tributary sites. In active floodplain sites, the Depth Index was similar to that of the seepage water zone, but the difference to tributary sites was not significant (Fig. 3, Table 3A). Neither season nor soil texture nor other environmental variables (PCs) modified this effect.

Within the active floodplain, forest developmental stage, soil texture, PC2 and PC3 were not included in the final model as single factors (only as interaction terms) due to higher AIC values in models containing these variables. Therefore, the effect of these variables individually was considered to be non-significant. As in the hydrological gradient, soil fauna feeding activity decreased with soil depth in all forest development stages (p < 0.0001, Table 3B). Further, the decrease in feeding activity with soil depth was more pronounced during spring than autumn (soil depth x season interaction p < 0.0001), and in sandy compared to loamy soils (p < 0.0001 for spring vs. autumn and p = 0.0004 for loamy vs. sandy, Fig. 4). General feeding activity rates were higher in autumn than in spring (only included as random term in the model), but significant differences between forest developmental stages were only observed in spring (forest developments x season interaction p < 0.0001) with lower feeding activity rates in grasslands than in old dense and young forests in spring.

Further, in grasslands on sandy soil, feeding activity was lower than in dense old and young forests on sandy soil (p = 0.01), while on loamy soil the differences between forest developmental stages were not significant. In addition, in grasslands on sandy soils feeding activity rates decreased stronger with soil depth than in young forests (interaction soil depth x forest development: p = 0.01; comparison young vs. grassland across depth with p = 0.029; all other comparisons p > 0.05). For soil fauna feeding activity in the individual sampling sites, see Suppl. material 5.

Forest development significantly affected Depth Index, i.e. maximum average soil-fauna activity occurred at a larger average soil depth in young forests than in grasslands (p = 0.04), while average soil depth of the activities in old dense and old sparse forests were between these values (Table 3B, Fig. 5). Soil texture affected Depth Index, resulting in maximum average soil-fauna feeding activity at greater soil depths in loamy compared to sandy soils. The effect of forest development on the Depth Index was modified by season and soil type. While in spring the mean Depth Index was lower in grasslands than old dense and young forests (but not old sparse forests), no difference was found in autumn. While no significant difference was found between forest development stages on loamy soils, the mean Depth Index was significantly lower in grasslands compared to young, old dense and old sparse forests on sandy soils (p = 0.03).

The feeding-activities observed here for the most part conform to other studies, confirming the suitability of the Bait-Lamina test in the current study, although adjustment of exposure time would have been necessary during autumn 2019 (see below). Soil-fauna feeding activity significantly decreased with soil depth and increased from spring to autumn in both gradients, similar to findings of earlier studies (e.g., Geissen et al. 2007; Rozen et al. 2010). This was expected, since decreasing soil fauna density with soil depth results in lower measured activity in deeper soil (Birkhofer et al. 2011). Soil fauna feeding activity has also been observed to increase with temperature (Gongalsky et al. 2004; Rozen et al. 2010; Simpson et al. 2012; Klimek et al. 2015), suggesting potential higher activity in early autumn and at sun-exposed soil surfaces.

Nonetheless, in the forest development gradient, feeding activities in autumn were close to 100% in all forest and grassland sites (see Fig. 4), indicating that the ISO-recommended exposition time of 4 weeks during autumn was too long to reveal differences in feeding activity. In the hydrological gradient, on the other hand, four weeks exposition was appropriate in early autumn, since the measured activities remained below 90%, suggesting that the slightly lower temperatures in this sampling year reduced overall activity (Suppl. material 1 for comparison of soil temperatures). Interestingly, soil fauna density and species richness did not significantly influence observed feeding activities, confirming previous reports that such community parameters cannot serve as proxies for feeding activity (see Introduction).

Besides frequency of inundation, the hydrological regime was strongly associated with further environmental factors, separating (1) forests of low leaf-area index, i.e. sparse canopy cover, on soils of neutral pH in the active floodplain from (2) very dense forests on strongly acidic soils in the tributary floodplains and (3) forests in the seepage water zone with intermediate canopy cover and soil pH. The factor “hydrological regime” therefore represented a number of environmental factors and was expected to significantly affect soil fauna feeding activity as well. However, a general overall effect was not confirmed and no significant differences in average feeding activity were found between hydrological regimes. Soil fauna feeding activity was characterized by high variability between different sites of the same hydrological regime, indicating that further factors besides general inundation type (and associated environmental conditions) affected soil fauna activities in our study. In addition, since no severe flooding event in the Elbe floodplain had occurred for 6 years prior to our study, we suggest that no effect of direct inundation persisted in the active floodplain sites, which would be in line with findings by Russell (2008) showing changes in the soil fauna community structure persisting only for several months after an inundation event.

Although overall feeding activity of soil fauna did not differ between hydrological regimes, the decrease in activity with soil depth was significantly stronger in tributary sites than sites of the seepage water zone or active floodplain, in particular in autumn. This indicates that soil-fauna activity was prevented deeper in the soils of the tributary sites, but not in the other hydrological regimes (cf. Rozen et al. 2010). High ground-water levels and resulting anoxia in deeper soil layers are a plausible explanation for restricted faunal vertical movement – and thus their low activity – in these sites. It is remarkable that anoxia by high ground water level had occurred mainly in the floodplain of tributaries, but not in the active floodplain, which contradicts our expectations. We therefore assume that the stagnant water in sites behind the dyke is much more negatively affecting soil fauna feeding activities than the running (flood) water in the active floodplain. In addition, anoxic conditions reduce mobility of mobile nutrients and increase toxicity of aluminium (Geissen and Brümmer 1999; Geissen et al. 2007), thereby potentially reducing soil fauna density and activity. In our study, hydromorphic soil characteristics indicated that anoxic conditions – due to a high groundwater table – occurred at rather shallow soil depths in sites of tributaries as well as seepage water zone (at 35 and 40 cm, respectively, compared to > 100 cm in the active floodplain soils; Vásconez Navas et al. 2023). Surprisingly, a reduction of soil fauna activity in deeper soil only occurred in tributary floodplain sites but not the seepage water zone, indicating at least temporal differences in soil water saturation despite occasional anoxic conditions in both hydrological regimes. All investigated tributaries are regulated for agricultural purposes (pers. comm. A. Gröngröft). We therefore assumed a complete independence of water levels between the Elbe River and its tributaries with anoxic conditions by high groundwater table probably persisting throughout the summer in the tributary floodplain soils. In sites of the active floodplain and seepage water zone, on the contrary, low water levels of the Elbe River in summer resulted in lower groundwater tables, allowing soil fauna to migrate deeper. Consequently, our results indicate a higher impact of recent hydrological effects (i.e. high groundwater tables) than of long-term inundation effects.

Soil texture alone apparently did not affect soil fauna feeding activity in our study, but the absence of sandy tributary sites hampered a comparison of soil-texture effects between all hydrological regimes. In the seepage water zone higher soil fauna activities were observed in sandy compared to loamy soils. This indicates better abiotic conditions for the soil fauna in sandy soils in the seepage water zone with probably better aeration and warming of the soil surface, in spite of the lower water holding capacity and therefore lower mean soil moisture of sandy soils (Bardgett 2005). However, this effect was restricted to the seepage water zone, since activity rates were similar in sandy and loamy soils in sites of the active floodplain, suggesting that in the active floodplain different parameters might affect soil fauna activity. Furthermore, when comparing activity rates in sites of the seepage water zone to those of the active floodplain, the positive effect of sandy soil becomes only visible during spring. During early autumn, on the contrary, activity rates were similar in the active floodplain and sandy seepage water zone sites, but there was a negative effect of loamy soil in sites of the seepage water zone. Our sampling campaigns took place in years of extreme drought in Europe, leading to very dry soils in the investigated sites. In particular, in the seepage water zone sites with loamy soil, very low soil water content was observed, and in addition loamy soils became very dry and hard during summer. Therefore, the positive impact of increasing temperature on soil fauna during early autumn probably was limited by a decrease in soil moisture, as soil fauna is known to avoid dry conditions by retreating to deeper soil horizons or inactive stages (Gongalsky et al. 2008; Thakur et al. 2018). In conclusion, abiotic soil conditions seemed to affect soil fauna activity in particular in the seepage water zone, indicating a limitation of soil fauna feeding activity by soil temperature in spring, but by soil moisture in autumn.

Since in the present study all sites of the forest development gradient were located in the active floodplain, they represented a habitat diversification in this specific hydrological regime. However, environmental and vegetation parameters varied mainly between grassland sites on the one hand and forest sites on the other hand. This was expected to be reflected in the observed soil fauna activity. Earlier studies did not find differences in soil fauna feeding activity between central European coniferous and deciduous forests of different ages (Geissen et al. 2007; Klimek et al. 2015), but between grasslands and forest sites in general (Hamel et al. 2007). In our study, forest development also did not directly affect soil fauna feeding activity and we found neither a significant difference in overall feeding activity between grasslands and forests, nor between forests of different age or density. Instead, we found that forest development modulated the effect of soil depth, i.e., a steeper decrease of activity with soil depth was observed in grasslands than in young forests. Differences in soil fauna feeding activity between grasslands and young forests only occurred at > 2.5 cm soil depth, and DI values suggested that in young forests mean soil fauna activity occurred at a greater soil depth as compared to grasslands. We therefore conclude that in grasslands soil fauna activity was restricted to the soil surface by environmental factors, probably by rapidly declining soil moisture with soil depth due to intense water uptake of grass roots (pers. obs.). On the other hand, in forests (in particular young forest plantations) soil moisture was likely appropriate for soil fauna feeding activity at least to 9 cm depth.

Soil texture also did not affect soil fauna feeding activity directly in this gradient and, due to very high general feeding activities in autumn, significant differences between soil fauna activities in different forests or grassland sites were only observed in the spring sampling. In old dense and young forests, feeding activity was higher in sandy compared to loamy sites. The most probable explanation for this result is the lower exposition of the soil to sunlight in denser forests. Old dense and young forests exhibited lower soil temperatures (mean difference of ~0.5 °C) and higher moisture (mean difference of ~5%) in loamy sites compared to sandy sites, while in open habitats (sparse forests and grasslands) these differences were less pronounced. This indicates that in shady habitats, sandy soils provided better habitat conditions for soil fauna than loamy soils due to faster warming (before development of canopy cover during spring) and reduced potential anoxia, while the close canopy cover prevented desiccation of the soil surface (Thakur et al. 2018). On the contrary, in grasslands and old sparse forests the effect of soil texture was restricted to deeper soil with a marked decrease in soil fauna feeding activity in sandy but not in loamy soils. This was surprising since the higher porosity of sandy soil, resulting in higher oxygen content in deeper soil, was expected to increase soil fauna activity in deeper soil (cf. Hassink et al. 1993). Soil moisture and temperature were similar in loamy and sandy soils in open habitats, indicating that the causes of similar feeding activity rates at the soil surface, but differences in deeper soil, might not lie in exposure to sunlight. Instead, differences most probably were related to site-specific effects e.g. pH or groundwater level, that resulted in strong differences in feeding activity rates between sites of individual forest developmental stages (Suppl. material 5).