The study area was located in the forest Monte el Viejo in Palencia (Northern Spain; UTM 30N coordinates: X: 370265, Y: 4646195), which covers around 1,500 ha in a moor at 865 m.a.s.l. The area is considered as a regressive phase of a Mediterranean oaks forest (i.e. Quercus ilex L. and Quercus faginea Lam.) (Figure 1). The area has been progressively transformed from forest to pastures (sheep grazing) and cereal croplands, from >3,100 ha in the mid-19^th century to 1,435 ha in 1968. More recently, some high-sloped areas unproductive for croplands because of soil erosion were afforested with Pinus spp. and Cupressus spp. to serve as intensely managed oak forests for firewood production specially until late 1960s (Alario et al., 1981; Zorraquino, 2015). Currently, this area mainly hosts recreative activities (hiking, cycling or birdwatching) since it is located ~6 km from the city of Palencia (76,738 inhabitants in 2024). Consequently, the area has intense silviculture interventions aimed to reduce the incidence of wildfires, as well as to keep an irregular state of stands (landscape silviculture). This forest is mainly formed by stands of both oak species (Q. ilex and Q. faginea) with patches of Mediterranean forest understory. Local soils are classified as Leptosol-lithic (LPli) and Regosol-calcaric (RGca) according to FAO´s World Reference Base (WRB) (IUSS Working Group WRB, 2015) and ITACYL (https://suelos.itacyl.es/mapas).

Monte el Viejo landscape is in a moor and in general lacks a predominant orientation, the study area is located on an average slope of 21% and southeast and northwest orientations (Source: Instituto Geográfico Nacional -IGN-; https://centrodedescargas.cnig.es/CentroDescargas/index.jsp#). Average altitude is around 822 m.a.s.l. in our research area. The local climate is characterized by intense summer drought (73 mm average) and cold winters with little rainfall (134 mm as average), the mean annual temperature is 11.1 ºC and the mean annual precipitation is 471 mm approximately (http://agroclimap.aemet.es/#; automatic station Palencia-Autilla del Pino, call sign: 2400E).

The sampling was carried out in April 2023. This period is characterized by mild temperatures (~15 ºC) at the beginning of the rainy season. The methodology for selecting each soil sample followed a systematic sampling (point) approach so that each selected plot (size 40x20 m each plot with 3 points in a straight line separated 10 m from each other; Figure 1) met the following criteria proposed by Menta et al. (2018a): (i) to take at least 3 replicates of each sampling plot, (ii) to keep >10 m from any anthropic constructions, and (iii), to keep ≥10 m between soil sampling points to record the heterogeneity of soil. Soil samples for microarthropods identification (1 disturbed sample/point; 3/plot: 9 samples in total) and for soil characterization (3 undisturbed and 1 disturbed/point; 36 samples in total) were collected. The sampled plots were selected in three biotopes: (I) secondary grassland in abandoned cropland that over time have been spontaneously colonized by herbaceous community; (II) shrub community (mainly Thymus spp., Salvia spp., and Cistus laurifolius L.) over a hillside that lies midway between herbaceous layer at the bottom of the slope and a more settled holm oak (Q. ilex) woodland at higher altitudes; and (III) Mediterranean young oak forest dominated by Q. faginea as canopy species, with presence of regrowth and scarce understory layer.

Sampling plots (plot codes are shown in Table S1 [suppl]) were initially characterized assessing the following geographical and topographical variables using QGIS (https://qgis.org): point coordinates (ETRS89 UTM 30N); A: altitude using a Digital Terrain Model (Source: Instituto Geográfico Nacional -IGN-; https://centrodedescargas.cnig.es/CentroDescargas/index.jsp#); S: slope; and O: orientation in respect of north of the sampling point. In addition, five stand variables were also measured in each sampling plot: the fraction of tree canopy cover (Fcc) estimated with the most recent Light Detection and Ranging (LiDAR) coverage layer (2^nd coverage in 2019). More specifically, the plugins FUSION (http://forsys.cfr.washington.edu/FUSION/fusion_overview.html) and LASTOOLS (http://rapidlasso.com/LAStools) were used to calculate a grid surface (DTM) of 10x10 m cells over the study area and on this grid the forest canopy model was created using classes of the DTM to calculate the percentage of land covered by trees; D: the average of diameter at chest height of the 5 nearest trees to each sampling point measured with a manual forceps; H: the average height of the 5 nearest trees measured with Suunto PM-5/1520 PC Clinometer (Suunto, Finland); Hs: the average height of the 5 nearest shrubs, and percentage of alive herbaceous coverage (Hc) measured in a 1 m² units randomly located surrounding the sampling points (see abbreviations in Table S1 [suppl]).

The study of the soil microarthropods was performed according to Parisi et al. (2005) and Menta et al. (2018a) with minor modifications. Briefly, the soil samples (1 per sample point, 3 per plot, 10 cm³) were taken with a small shovel after removing the litter layer and separately introduced in labelled hermetic sealing plastic bags. Samples were immediately taken to the laboratory at room temperature and transferred to Berlese-Tüllgren extraction devices (i.e. one soil sample in each device) where they were subjected to natural light photoperiod and constant heating through 60 W incandescent infrared lamps (GiganTerra, France) placed 30 cm above each sample for one week. Each Berlese-Tüllgren devices included 2 mm light rectangular metal mesh and a plastic tray. Microarthropods naturally present in the soil migrate deeper falling in a preservative solution (2/3 ethanol 96% v/v and 1/3 glycerol 99% v/v) until morphological identification.

Taxonomic identification of the soil microarthropods [that is, arthropods with body size between 0.1 and 2 mm, (Nielsen, 2019)] after 7 days in Berlese-Tüllgren devices (Barrientos et al., 2004) was performed. The soil microarthropods preserved were transferred to 4.5 cm diameter Petri dishes to facilitate their observation under a stereomicroscope (MOTIC SMZ-168, 7.5-50x). The sampled individuals were photographed with a Moticam 580 5.0 MP camera connected to the stereomicroscope before being transferred to 2 ml tubes containing ethanol 70% v/v using a Pasteur pipette. Finally, the liquid was poured into zooplankton counting chambers (Aquatic BioTechnology S.L., Spain) where the specimens were classified up to Order or Class level.

The QBS-ar soil quality index was determined based on the presence of specific taxonomic groups, with each group assigned an ecomorphological index (EMI) score reflecting its adaptation to the soil environment from 0 to 20. Euedaphic taxa, fully adapted to belowground life, received the highest score of 20, while hemi-edaphic and epi-edaphic organisms received progressively lower scores (Mantoni et al., 2020). If different adaptation levels were recorded within the same taxonomic group, the highest EMI value was used (Parisi et al., 2005). The QBS-ar was the sum of all EMI values for each sample. For Collembola, the scores provided by Menta et al. (2018a) were adapted to the taxa found in the study area (Rodríguez-Pajares, 2023).

Absolute abundance was estimated after microarthropods identification (Table 1). The number of taxonomic groups (that is, taxonomic richness; Sobs) per sample, Shannon´s diversity index (Hsw), and Pielou´s evenness index (Ep) were calculated using the package “vegan” (https://CRAN.R-project.org/package=vegan) in R programming environment (https://www.R-project.org/); all these indexes are shown in Table 2. Dominance of each taxonomic group in each plot was determined using the Camargo´s index (Camargo, 1993), indicated in bold in Table 1. The rank-abundance curves were computed by sampled plot in R to identify those taxonomic groups that provided higher absolute abundances within biotopes. In parallel, the IndVal index (Dufrêne & Legendre, 1997) was computed to detect whether any taxonomic group played a decisive role as indicator of the biotope sampled. This index was calculated using the package “indicspecies” (https://cran.r-project.org/web/packages/indicspecies/index.html) in R with 999 permutations.

The variations of community composition by biotopes were evaluated by computing a Non-metric Multidimensional Scaling (NMDS) based on Bray-Curtis dissimilarity using the package “vegan”. Compositional differences were visualized in the ordination space using the package “ggplot2” (https://cran.r-project.org/web/packages/ggplot2/index.html) in R. In addition, Bray-Curtis dissimilarity by biotope and soil texture was further investigated by computing PERMANOVAs using the package “vegan” with 999 permutations.

Physical features of each sampled soil were evaluated. First, the surface parameters in situ [i.e. surface stoniness (SuS); leaf litter weight (Lw) in 0.50x0.50 m area, and leaf litter cover (Lc) in the same surface] were taken by placing 4 steel rods of 1 m in length in a representative area forming a square of 1 m².

Two types of soil samples, i.e. undisturbed and disturbed, were taken to the laboratory for further analysis. The former were extracted using Blake's cylinder (5 cm length and 4.75 cm diameter) from the first 10 cm of the soil, selected at random around the same sampling point used for the collection of soil arthropods (3 per sampling point, i.e. 9 in total per plot). This method is used to calculate the bulk density (Bd) according to Blake (1986), as well as the water content of soil (Hu), following the method proposed by MAPA (1994).

Otherwise, the disturbed soil samples (3/plot) were taken into the first 10 cm soil using a small shovel obtaining quite disturbed soil together with the coarse elements, avoiding leaf litter. Fine soil was obtained following the protocol of MAPA (1994) to obtain soil stoniness (SoS), as well as texture (T) according to feel procedure (USDA, 1999). Real density (Rd) was also calculated using pycnometers as described by Blake (1986). In addition, porosity (P) was measured as soil volume occupied by the pores by difference between Bd and Rd. The organic matter content (OM) was determined according to Schulte (1996).

All the categorical variables from plot characterization were transformed into binary to develop a principal component analysis (PCA). The PCA was performed using “FactoMineR” (https://cran.r-project.org/web/packages/FactoMineR/index.html) and “factoextra” (https://cran.r-project.org/web/packages/factoextra/index.html) packages in R to detect tendencies that cluster sampling plots by environmental gradients. The percentage of explained variability was evaluated by each dimension, as well as the relative contribution of each variable in explained variability of the corresponding dimensions.

The possible bias caused by sample weight was evaluated computing a generalized linear model (GLM; Gaussian distribution of errors) in R considering the indicators Sobs, Hsw, Ep, and QBS-ar as response variables in the corresponding GLMs, being the sample weight the explicative variable. The possible variations in ecological indicators by biotope (plot) and environmental factors previously identified by PCA (i.e. Hs, Fcc, and H) were evaluated by computing additional GLMs. Environmental factors were included in the models as explicative variables alone, in combination, as well as interaction whenever it retained biological meaning. Eighty-four models were computed (21 per indicator) being then compared among them and against a null model (a simple model where the ecological indicator was explained only by the plot) through Akaike's Information Criteria (AIC) using the “AICcmodavg” package (https://cran.r-project.org/web/packages/AICcmodavg/index.html). Final selected model of each indicator was the most parsimonious one (lowest AIC value). The fitness quality of selected model was evaluated using the package “DHARMa” (https://cran.r-project.org/web/packages/DHARMa/index.html) and predicted values were visualized with “effects” package (https://cran.r-project.org/web/packages/effects/index.html). When an explicative factor resulted in a significate effect, the LSD Fisher test was performed using “agricolae” package (https://cran.r-project.org/web/packages/agricolae/index.html) as post-hoc analysis to detect variations by variable levels (see abbreviations in Table S1 [suppl]).

The sampled biotopes (i.e. secondary grassland, shrublands, and young oak forest) showed noticeable differences in terms of soil and stand features. Specifically, the woody vegetation strata are more developed as one moves from grassland to shrubland to forest, in that order. However, the opposite is true for grasses: the higher the cover of shrub and tree strata, the lower their density. The understory layer was more developed in the shrubland than in forest (Table S2 [suppl]), and in turn the trees in the shrubland were younger, as evidenced by tree diameter and height.

In terms of edaphic description, Lc, Lw, Hu, and OM, were higher in the forest than in the other biotopes (Table 3 and Table S3 [suppl]), due to the greater presence of trees that provide OM coming mainly from their leaves. The texture, although variable, shows in all cases to be rather loamy, accompanied by some stoniness in all cases. Porosity is high in all cases, but more noticeable in the forest, which allows for more aerated and spongy soils.

Soil microarthropods collected from Berlese-Tüllgren extraction devices were identified as belonging to 14 taxonomic groups (Table 1). Acari were the most abundant (44%), followed by Hymenoptera (34%, primarily epigeic ants from the Formicidae family) and Collembola (15%). The remaining taxa (Araneae, Pseudoescorpionida, Chilopoda, Diplopoda, Pauropoda, Protura, Psocoptera, Hemiptera, Lepidoptera, Diptera and Coleoptera) were underrepresented (< 2%). Nematodes and annelids (Oligochaeta) were also found in the samples, but they were not included in the QBS-ar method, and thus not considered in the study.

The analysis of community composition revealed mites (Acari: 243.33 ± 58.52 individuals/m², mean value and standard error), springtails (Collembola: 80.00 ± 26.03 individuals/m²) and beetles (Coleoptera: 11.11 ± 4.23 individuals/m²) as ubiquitous in the sampled biotopes. However, Collembola with a higher degree of adaptation to edaphic life were especially abundant in the forest (116.67 ± 41.77 individuals/m²). On the other hand, Hymenoptera, also quite abundant (i.e. 169.03 ± 104.73 individuals/m²) were almost exclusive to more open areas (shrubland and grassland). Otherwise, proturans (Protura) and the one caterpillar (Lepidoptera) specimen only occurred in the oak forest. Specimens of Psocoptera were found in the grassland and in the forest. Pseudoscorpionida, Chilopoda and Diptera were found in the forest and shrubland meanwhile Hemiptera was the only group recorded both, in grassland and shrubland. Finally, the underrepresented groups were Araneae, Pauropoda, Diplopoda and Lepidoptera (i.e. mean abundance per group ≤ 1.11 individuals/m². Table 1 and Figure S1 [suppl]).

The IndVal index highlighted the ants as indicators in secondary grassland and shrubland (p-value=0.01). This indicates high specificity, as most Formicidae ants are predominantly found in grasslands compared to other habitats, and high fidelity, as these ants are present in most sites within the grassland, indicating a strong association with this habitat type.

Regarding ecological indicators, QBS-ar, Sobs and Hsw showed almost homogeneous values (Table 2). More specifically, the shrubland yielded the highest QBS-ar value from 35 to 101 with a mean of 68.67. On the contrary, the lowest one between 28 and 77 was found in the grassland (Figure 2), between 28 and 77 with 52.33 average points. The forest plot was the one that varied the least.

Sobs and Hsw were close among the different biotopes sampled (5.44 ± 0.60 and 1.08 ± 0.08 points, respectively). Sobs reached its maximum (8) and minimum (3) value, both repeated in the grassland and in the shrubland (Table 2). The difference between these values was 5 points, with a mean of 5.67 points in the grassland and somewhat lower (5.33) in the shrubland while in the forest the values have been much more homogeneous (5 and 6 points, with a mean of 5.33). On the other hand, the highest Hsw was in the shrubland with a value of 1.59, which decreased in this same plot to 0.79, giving a mean value of 1.20. The range of values of the plot with the lowest Hsw, grassland (1.19-0.79, with a mean of 1.03), was higher than forest, which hardly changed (1.04-0.93 points, with a mean of 1.00) (Table 2).

Pielou´s evenness (Ep) values were similar across sampling points and habitats but varied the most in grasslands. The shrubland had the highest average Ep value, which was 0.75 (from 0.72 to 0.76), about 0.11 points higher than the grassland and 0.16 points higher than the forest, and, which had the lowest average value.

The analysis of PCA revealed that H represents 98.2% of the variability (first dimension of PCA; Figure 3) being the rest of variables less explicative (< 1.8% of explained variability each one). The variables Hs and Fcc represented 73.9% and 8.3% of explained variability in dimension 2, respectively. Accordingly, these stand parameters were selected as plausible explicative variables in data modelling. Otherwise, the NMDS showed a noticeable clustering of forest samples in terms of Bray-Curtis dissimilarity while shrubland and secondary grasslands were almost similar in terms of community composition (Table S4 [suppl] and Figure S2 [suppl]). The analysis of dissimilarity in sampled communities by soil texture and biotope using PERMANOVA did not show significant variations (p-value > 0.06 in both cases).

According to the fitted models (Table 4), the indicator Ep significantly varied according to shrub height (Hs; model M76; p-value < 0.01). Model M76 resulted the most explicative showing noticeable lower AIC value than M64 and the null model (> 2 AIC in both cases; Table 4). According to the model M76, the value of Ep was expected to be higher as Hs increases (Figure 4).

The original hypothesis of this study considered that ecological indicators could by driven by the biotope, nevertheless, the inclusion of this variable (plot) in M76 did not increase the quality of the model (see description of M66, Table 4). To evaluate this possible source of variability, the variation of Ep by this factor was evaluated through a null model (M64) showing a significant effect according to model description (p-value < 0.01) and post-hoc analysis despite M64 resulted in a higher AIC value than selected model M76 (ΔAIC = 4.30, Table 4). The weight of the soil sample did not reveal any significant variation for the four ecological indicators (p-value ≥ 0.57 in all cases) according to computed GLMs, therefore this variable was excluded from models’ computation.