Planting food forests can increase soil biodiversity in agricultural landscapes of Northwest Europe

Soil organisms drive ecosystem processes that underlie key benefits to humans such as nutrient cycling, plant production and soil carbon storage^1,2. Soils and their biota are, however, subject to a variety of threats, of which intensive agricultural practices are considered among the most important^3,4. Therefore, there is a strong interest in shifting to food production systems that preserve or even increase soil biodiversity, either by using less intensive management practices⁵, or by implementing radically different approaches to food production. In the temperate zone, one of these alternative approaches that is currently gaining traction is the use of complex multilayered polycultures called food forests⁶.

Temperate food forests are agroforestry systems structurally designed in analogy to semi-open forests or forest edges that are characterised by a high diversity of perennial plants occupying different vegetation layers and a high habitat heterogeneity^7,8. They are typically managed without synthetic fertilisers or pesticides and with limited soil disturbance⁷. Many pioneering food forests in Northwestern Europe have passed an initial establishment phase and are now older than five years. This provides a timely opportunity to assess if installing food forests could help increase soil biodiversity in agricultural landscapes. To evaluate this, we need insight into how the biodiversity of a range of taxonomic groups of soil organisms compares between food forests and the dominant land use types in these landscapes^9,10: the intensively managed permanent grasslands and croplands they replace and the forests to which they bear structural resemblance⁷.

In temperate regions, many taxonomic groups of soil organisms, including fungi and soil fauna, typically reach a higher total biomass or abundance in systems with lower levels of soil disturbance, higher levels of litter input, and a more humid microclimate^11,12. Therefore, food forests can be expected to harbour a higher total biomass or abundance than grass- and croplands for these groups. For the same reason, the total biomass or abundance of these groups may still be lower in food forests than in forests, where the canopy is usually more closed and higher litter inputs could provide more habitat and resources^13,14.

Food forests may not only differ from grasslands, croplands and forests in terms of the total biomass or abundance of soil organism groups, but also their community composition, as species differ in resource and microhabitat requirements and tolerances to factors like low pH, desiccation and soil disturbance^15,16,17. The presence of shrubs and trees and the relatively low soil disturbance in food forests could provide suitable habitat for typical forest floor inhabitants such as oribatid mites¹⁸ and filamentous fungi¹⁹. However, food forests are typically maintained as semi-open systems with herbaceous patches^7,8, and their soils might still hold the biotic and abiotic legacies of the agricultural land they were planted on²⁰. They therefore potentially harbour species more associated with grass- and croplands as well (e.g. smaller, fast-growing species²¹ and more drought-tolerant species¹⁶). As a result, soil communities in food forests may consist of a mix of grass- and cropland- as well as forest-related species.

By offering a high diversity of resources and microhabitats, food forests are expected to also support a high diversity of soil organisms²². At a small spatial scale, the high plant diversity and complex vegetation structure in food forests could allow co-existence of a large number of species (i.e. high alpha-diversity)^23,24. At the scale of the entire food forest, these systems are typically also highly heterogeneous, consisting of patches with diverging vegetation structure and composition⁷. This may lead to a higher difference in soil community composition among locations within a food forest (i.e. a higher beta-diversity) than within forests, grasslands and croplands which may be more homogeneous. However, how these aspects of soil biodiversity differ between temperate food forests and other land use types is in fact almost fully unknown (but see Kreitzman et al. (2022)²⁵).

We examined the possible contribution of food forests to soil biodiversity in intensively used agricultural landscapes by comparing the biodiversity of a wide range of taxonomic groups of soil organisms between 15 food forests and nearby forests, grasslands and croplands in Belgium and the Netherlands. We evaluated the total biomass or number of individuals (depending on the group studied), alpha-diversity, community composition and beta-diversity of 12 taxonomic groups that cover different phylogenetic domains, size classes and trophic levels in the soil mineral and litter layer: bacteria, fungi, protists, nematodes, microarthropods (mites and springtails), macroarthropods (terrestrial isopods, millipedes, centipedes, harvestmen and carabid beetles) and earthworms. To study this comprehensive set of organism groups, we combined amplicon sequencing techniques for the identification of microorganisms and micro- and mesofaunal taxa with classical expert-based species identification of macrofaunal groups. We hypothesised that for most groups (1) food forests harbour a lower total biomass or number of individuals than forests, but higher than grass- and croplands; (2) alpha-diversity is higher in food forests than in forests, grasslands and croplands; (3) the community composition in food forests is intermediate between the composition in forests on the one hand and grass- and croplands on the other and (4) beta-diversity is higher in food forests than in forests, grasslands and croplands. Because soil physicochemical properties play an important role in the influence of land use on soil communities²⁶, we also compared key physicochemical properties between food forests and the reference land use types and explored how variation in these properties relates to variation in soil community composition.

Soil physicochemical properties

Our comparison of five key soil properties between food forests and the reference land use types (Table S1, Fig. S2.1 and Table S2.1) and a principal component analysis (PCA) based on these properties (Fig. 1) show that food forest soils were roughly intermediate between those in forests on the one hand and grass- and, especially, croplands on the other. Compared to forests, soil bulk density, pH and total phosphorus (P) content were higher in food forests, while total carbon (C) content was lower. Conversely, food forests had a lower bulk density and higher total C content than grass- and croplands, as well as a lower pH and total P content than croplands. Soil pH was somewhat higher in food forests than in grasslands, while the difference in total P content was limited and not statistically significant (p ≥ 0.05). Overall, total nitrogen (N) content was on the higher side in food forests, but only the difference with croplands was statistically significant. Summary statistics of additional chemical soil properties are given in Table S2.2.

Biomass or the number of individuals

We found the following statistically significant (p < 0.05) differences between food forests and the reference land use types (Table 1). In comparison to forests, the biomass of non-arbuscular mycorrhizal (non-AM) fungi was lower in food forests (Fig. 2b). In contrast, the biomass of bacteria and AM fungi and especially the number of millipedes and earthworms were higher in food forests than in forests (Fig. 2a, c, h, l). Compared to grasslands, we found a higher biomass, or higher number of individuals, in food forests for non-AM fungi and mites and particularly for isopods, millipedes, centipedes and harvestmen (Fig. 2b, f–j). However, we found a lower biomass of AM fungi and fewer nematodes and earthworms in food forests than in grasslands (Fig. 2c, d, l). Lastly, compared to croplands, we found a higher biomass of bacteria, non-AM fungi and AM fungi, and higher numbers of nematodes, as well as much higher numbers of springtails, mites, isopods and millipedes in food forests (Fig. 2a–h). In contrast, fewer carabids were captured in food forests than in croplands (Fig. 2k).

Alpha-diversity: the effective number of species/taxa (¹D)

For the effective number of species/taxa, the following differences were statistically significant (Table 2). Compared to forests, the effective number of species/taxa of bacteria, non-AM fungi and earthworms was higher in food forests (Fig. 3a, b, l), but for millipedes and harvestmen it was lower (Fig. 3h, j). We found a higher effective number of species/taxa in food forests than grasslands for bacteria, non-AM fungi, isopods, millipedes and harvestmen (Fig. 3a, b, g, h, j). In comparison to croplands, we found a higher effective number of species/taxa in food forests for non-AM fungi, isopods, millipedes and harvestmen (Fig. 3b, g, h, j). Overall, patterns in species/taxon richness ⁰D were similar to those in the effective number of species/taxa ¹D (for details, see S3).

Community composition

Variation in community composition was related to differences in both land use type and the geographical location of the study sites, but the proportion of variation attributable to differences in land use type (R² = 1-15%) was always lower than the variation attributable to differences in location (R² = 21-65%) (Table S4). For earthworms, the PERMANOVA analysis did not reveal differences between food forests and forests (p ≥ 0.05), but did indicate differences in community composition between food forests and grass- and croplands, while the ordination plot showed overlap of the land use types, with only a few distinct food forest and forest samples (Fig. 4g). For all other groups, there was a statistically significant difference in community composition between food forests and each of the three reference land use types. In the case of non-AM fungi, macroarthropods and to a lesser extent protists and microarthropods, food forest communities were generally intermediate between those of the forests on the one hand, and those of the grass- and croplands on the other hand (Fig. 4b, d, e, f). Bacteria and nematode communities in food forests generally showed a greater similarity to those in grass- and croplands than to those in forests, which were often markedly different (Fig. 4a, c). These patterns are also visible in differences in the mean proportion of particular taxa (classes, families or species) between food forest sites and sites of the reference land use types (S6). Notable differences are the higher proportion of non-AM fungi from the Ascomycota classes and lower proportion of Agaricomycetes in the food forests than the forests and, for macroarthropods, the lower proportion of Carabidae and higher proportion of Isopoda families in food forests than grass- and croplands (S6).

Our exploration of the relationships between bulk density, pH, total C, N and P content and the patterns in community composition summarised in the ordination (Fig. 4) is shown in Figs. S7.1 and S7.2. In some cases, variation in the soil property aligned with the main direction of differences in community composition between land use types. In particular differences in pH, C and P seemed associated to the difference between food forests and forests for bacteria, while for nematodes this was the case for bulk density and total P. For non-AM fungi, especially variation in bulk density and total P appeared related to the overall differences between land use types, and for protists and microarthropods this was most notable for bulk density.

Beta-diversity: average pairwise Bray-Curtis dissimilarity

When comparing the beta-diversity between food forests and forests, we only found a statistically significant difference for the pooled macroarthropod group, for which the beta-diversity was higher in food forests than in forests (mean 0.61 vs. 0.45, Fig. S8f). The difference between food forests and grasslands was only statistically significant for the nematodes, for which we also found a somewhat higher beta-diversity in food forests (mean 0.64 vs. 0.53, Fig. S8d). The statistical comparison between food forests and croplands could only be conducted for the bacteria, non-AM fungi and macroarthropod groups. Here, only the beta-diversity of non-AM fungi was statistically significantly higher in food forests than croplands (mean 0.62 vs. 0.46, Fig. S8b).

We investigated how planting food forests could contribute to increasing soil biodiversity in agricultural landscapes by comparing total biomass or number of individuals, alpha-diversity, community composition and beta-diversity of 12 groups of soil organisms between food forests and nearby forests, grasslands and croplands. There was a high variability within land use types that can be related to differences in environmental conditions, but also to variation in age, size, land use history, vegetation and management. Nevertheless, we detected clear differences between food forests and the other land use types in almost all soil properties and in many soil organism groups’ total biomass or number of individuals and community composition, as well as some differences in alpha- and beta-diversity.

We expected to mostly find a lower total biomass or number of individuals in food forests than forests. However, we instead found similar or higher values in food forests than forests for all groups except non-AM fungi. This suggests that although food forests typically have a more open vegetation structure and canopy than forests, this does not (per se) lead to an overall reduced availability of resources or habitat for soil organisms. Possibly, a lower quantity of tree litter in food forests than forests is compensated for by a higher quality (e.g. lower C/N ratio) of these inputs, resulting from the high proportion of plants with N-fixing symbionts and high-quality litter^7,16,27. A higher litter quality may not only provide easier access to nutrition for detritivores. It may also provide more suitable soil conditions for certain groups, such as the observed higher pH that can help explain the higher bacterial biomass and higher number of earthworm individuals in the food forests²⁸. The higher biomass of AM fungi may additionally be explained by a higher abundance of host species in food forests (e.g. Rosaceae spp.)^29,30. We did find a lower biomass of non-AM fungi in food forests than forests. This could result from inherent land use differences (i.e. in vegetation and management). However, it could also be due to the food forests overall being younger than the forests (Table S10), so that their soil communities have had less time to develop from more bacterial-dominated to fungal-dominated states, a trend that is often observed during forest development on agricultural land^21,31,32.

In line with our hypothesis, the total biomass or number of individuals in food forests was higher than in grasslands and often much higher than in croplands for more than half of the studied groups. These groups, i.e. non-AM fungi, mites, isopods, millipedes, centipedes and harvestmen, are typically found in higher quantities in less disturbed environments with higher litter availability and humidity^4,12,33, which is what food forests may offer. We did find higher AM fungal biomass and more earthworms in grasslands than in food forests. Presumably, the high abundance of fine roots in grasslands can support a greater biomass of AM fungi³⁴, whereas earthworms may benefit from a still higher input of easily degradable organic matter (e.g. manure) in grasslands than food forests³⁵. The higher numbers of carabids in croplands than in food forests may be explained by higher prey densities or a higher mobility/activity and, hence, higher capture rates in open landscapes³⁶. In an earlier study on young food forests (half were 1-5 years old), the number of individuals of litter-dwelling arthropods differed little between food forests and croplands²⁵. We now show that, in terms of both the number of individuals of soil fauna and microorganismal biomass, food forests older than 5 years do appear to differ from grass– and croplands and rather resemble forests.

We found little support for our hypothesis that the high plant and structural diversity in food forests translates into a greater alpha-diversity of soil organisms in food forests than in all reference land use types. We did find a higher effective number of taxa (¹D) of non-AM fungi in food forests than in all three reference land use types. This could be explained by the high diversity of symbiotic partners and substrates provided by food forests’ diverse plant communities³⁷ and the sensitivity of many fungal species to physical disturbance^12,38. Additionally, we found a higher effective number of species of isopods, millipedes, and harvestmen in food forests than in grasslands and croplands, which accords with the findings of Kreitzman et al. (2022)²⁵. While this could partly be influenced by a greater number of individuals captured in food forests (S9), food forests may offer habitat to a wider range of species thanks to their more diverse habitat, more humid conditions and higher diversity of litter³⁹.

Overall, the limited differences in the effective number of species/taxa between food forests and the other land use types contrast with our expectations but align with earlier studies reporting inconsistent relationships between belowground alpha-diversity and land use intensity²¹ or plant diversity⁴⁰. It could also be that, so far, the influx of species into the food forests has been limited due to a small species pool in the surrounding landscapes or that species have a low ability to reach the food forests due to habitat fragmentation⁴¹. In contrast to the effective number of taxa (¹D), the taxon richness (⁰D) of nematodes and protists was higher in food forests than other land use types. This could indicate that taxa have colonised the food forests which were not (yet) present in sufficiently high numbers to also noticeably contribute to the effective number of taxa.

Variation in community composition was more related to differences in geographic location than differences in land use type (S4), pointing to the importance of local factors such as soil type and the regional species pool in shaping soil communities^2,42. Nevertheless, in line with our third hypothesis, we found that community composition of non-AM fungi and macroarthropods, and to a lesser extent protists and microarthropods, in food forests was generally intermediate between the composition in forests on the one hand and grass- and croplands on the other. For non-AM fungi, this likely had to do with, among others, the proportion of Agaricomycetes being intermediate in food forests (Fig. S6.2). This is a fungal class that in addition to plant parasitic taxa contains the majority of ectomycorrhizal taxa as well as fungi that decay woody substrates termed white-rot and brown-rot fungi, which points to the potential role of the vegetation in shaping food forests’ fungal community composition^43,44. Our findings suggest that the intermediate composition in food forests could further be related to levels of P being intermediate in food forests (Fig. S7.1), which concurs with earlier studies finding an important influence of P concentrations on fungal community composition^45,46. For macroarthropods, the intermediate community composition in food forests could largely be due to their mix of more desiccation-resistant and -sensitive taxa. Notably, food forest samples had a lower proportion of carabids and a higher proportion of isopods compared to grass- and croplands³⁶, but still harboured a higher proportion of the desiccation-resistant Armadillidium vulgare than forests⁴⁷ (Fig. S6.6). For earthworms, the difference in community composition between food forests and grass- and croplands was most notable in the higher proportion of litter-dwelling earthworm species, especially Lumbricus castaneus (Savigny, 1826), in food forests⁴⁸.

Overall, we see three main reasons for which food forest communities may have been intermediate between those in grass- and croplands and those in forests. First, it could be due to food forests’ semi-open vegetation leading to a mix of microhabitats (e.g. in terms of microclimate)¹⁶. Second, litter inputs from food forest vegetation, with its high proportion of plants with N-fixing symbionts and resource-acquisitive species, may present an intermediate between the mineral inputs or easily degradable organic inputs in crop- and grasslands, and the typically more complex, recalcitrant substrates in forests⁴⁹. This can affect soil community composition directly as well as indirectly, through its impact on soil properties like N content and pH (S2)^19,49. Third, it may reflect how the food forests were overall planted more recently than the forests. Therefore, their communities may still be more influenced by the biotic and abiotic legacies of the previous agricultural land use (e.g. compaction and fertilisation, S2)^4,12,20. These legacies may also help explain why the community composition of bacteria and nematodes in food forests more closely resembled those in grass- and croplands than those in forests. For example, like in grass- and croplands, bacterial communities in food forests contained higher proportions of Bacilli (Firmicutes) than most forest sites (Fig. S6.1). It seems that differences in pH, C and P – properties which can exert important influence on bacterial community composition^50,51– could play a role in this difference between food forest and forest communities (Fig. S7.1).

We expect that food forest soil communities will differentiate further from grass- and cropland communities as they develop^52,53. Due to the paucity of food forests, notably older sites, in the region, we were not able to disentangle effects of age and other factors such as size, soil type and land use history (see Table S10 for the variation herein in our dataset) on food forest communities. Long-term monitoring in observational and experimental studies or chronosequence approaches are needed to gain insight into successional trajectories of food forest soil communities. Such studies can shed light on the influence of environmental conditions, the landscape context, land use history and management on soil community development in food forests, and on how soil community development compares between food forests and forests.

Lastly, we did not find strong support for our hypothesis that the high heterogeneity in plant community composition and structure in food forests leads to a higher beta-diversity of soil communities in food forest sites than in forest, grassland and cropland sites. Again, it may be that the rather short time since the establishment of the food forests and the potentially limited influx of new species did not allow for communities at different locations in the food forests to differentiate strongly yet⁵². In addition, differences between food forests and reference land use types may be smaller than expected because even in habitats with a homogeneous vegetation such as conventional croplands, the soil environment can be very heterogeneous at small scales²⁶.

Taken together, our findings indicate that planting food forests can modify and partly enrich soil communities in agricultural soils. Despite their relatively young age, the food forests in our study already appeared to harbour higher total biomass or number of individuals and sometimes also a higher diversity of several groups of soil organisms than grasslands and, especially, croplands. Alpha-diversity was, moreover, never found to be significantly lower in food forests than in grass- or croplands. Converting grass- or croplands to food forests may thus locally enhance soil biodiversity, particularly by promoting typically disturbance-sensitive and litter-associated groups such as non-AM fungi, harvestmen and isopods. Food forests further appear to support soil communities with a different composition than those in grass- and croplands. Consequently, planting food forests in landscapes dominated by grass- and/or croplands may also help increase the total soil biodiversity at the landscape scale. Moreover, by providing suitable conditions for drought- and disturbance-sensitive organisms, establishing food forests at strategic locations in the landscape could enhance habitat connectivity for these organisms, and thereby further support the viability of their populations that are under threat in intensively managed agricultural landscapes^4,54,55. On that account, food forestry could be considered as an approach to increase soil biodiversity in intensive agricultural landscapes without taking agricultural land out of production.

Differences between food forests and the other land use types did depend on the taxonomic group and aspect of biodiversity studied. Although the fact that we detected more clear differences for some groups than others could partially be attributed to differences in sample size between groups, we believe it is also driven by taxon-specific responses to land use⁵⁶. This highlights the importance of considering a broad range of taxonomic groups when aiming to understand impacts of land use on soil biodiversity. Future studies should take this a step further and look at the identity (and thus, for example, rarity) and functioning of the species finding habitat in food forests. Differences in soil communities between food forests and other land use types likely lead to differences in key soil functions¹. For example, a higher fungal diversity has been proposed to improve storage of stabilised C in the soil⁵⁷. Additionally, the higher numbers of litter-dwelling earthworms, isopods and millipedes in food forests seem promising for food production systems that rely more on nutrient cycling than mineral fertiliser input for plant nutrient provisioning^3,58. Lower numbers of carabids on the other hand may have implications for pest control⁵⁹. The current study uncovers important differences in soil biodiversity between food forests and dominant land use types in agricultural areas. Future studies can explicitly address how such differences impact soil functioning and, thereby, other aspects of agricultural sustainability.

Study region and study sites

This study was conducted in Belgium and the Netherlands, a region under temperate oceanic climate conditions with a mean annual rainfall of 800–900 mm and mean annual temperature of 10.5–11.5 °C in 1991–2020 (www.kmi.be, www.knmi.nl). We selected 15 food forests at 14 different locations (two food forests were situated at the same location) spread across the study region and occurring on different soil types, that matched the following criteria: a) consisting of a food producing canopy or subcanopy layer integrated with at least two other food producing vegetation layers; b) first planting at least 5 years before sampling; and c) planted on former grass- or cropland (Fig. 5a). When the food forests were visited for sample collection in 2021 and 2022, their ages ranged from 6 to 29 years and sizes from 0.18 to 3 ha (Table S10). We aimed to select three references per food forest: cropland, permanent grassland and mixed forest. In total, 38 reference sites were sampled instead of 45, mostly because cropland references could not be found for each food forest (Table S10). Reference sites were situated within 2 km from the food forest, on soils with the same texture and a similar drainage class as the food forest (see S2 for physicochemical soil properties of the food forests and reference land use types). Within each study site (food forest or reference), six plots of 5 × 5 m were selected using a stratified random design (Fig. 5b). Additionally, two sampling points for pitfall trap sampling were selected within 10 m from the centre of each half of the study site (Fig. 5b).

Soil physicochemical properties

Soil properties were measured in each of the six plots per study site between April and August in 2021 and between May and July in 2022. For bulk density (g cm^-3), an undisturbed soil sample of 100 cm³ was collected from the 0–15 cm mineral soil layer using a cylindrical core (Kopecky ring, Royal Eijkelkamp, the Netherlands) after which it was oven-dried at 105 °C for at least 24 hours and weighed. For the chemical properties, a composite soil sample of the 0–25 cm layer was obtained per plot by mixing 25 soil cores. These samples were dried at 40 °C for 48 h, ground and sieved over a 2 mm-mesh sieve. At the Netherland Institute of Ecology (NIOO-KNAW), pH-KCl and available NO₃^– and NH₄⁺ concentrations of the samples were determined by shaking 10 g of soil in a 1:5 ratio soil/KCl (1 M) mixture for 5 min at 300 rpm. pH was then measured with a pH meter Orion Star A211 with pH electrode model Ross sure-flow 8172 BNWP, Thermo Scientific Orion, USA. Afterwards, sample solutions were poured over a filter (Whatman Grade 1, 85 mm) and available NO₃^– and NH₄⁺ were measured colorimetrically with an AutoAnalyzer (SFA-system, SKALAR, the Netherlands). All other chemical analyses were performed at the Forest & Nature Lab of Ghent University (ForNaLab). Bioavailable phosphorous (P_Olsen) was extracted in NaHCO₃ (ISO 11263:1994(E)) and P_oxalate was extracted in ammonium oxalate-oxalic acid (POx, AlOx and FeOx; according to NEN 5776:2006). Total P was measured after complete destruction of the soil samples with HClO₄ (65%), HNO₃ (70%) and H₂SO₄ (98%) in teflon bombs for 4 h at 150 °C. All P-concentrations were measured through colorimetric measurement according to the malachite green procedure, using a spectrophotometer (Varian Cary 50 UV-Vis) at a wavelength of 700 nm. Exchangeable potassium (K⁺), calcium (Ca²⁺) and magnesium (Mg²⁺) were extracted in BaCl₂ (0.1 M; ISO 11260), after which concentrations were measured through inductively coupled plasma optical emission spectroscopy (Thermo Scientific™ iCAP™ 7400 ICP-OES). To determine the mass percentage of C and N in the soil, the samples were combusted at 1150 °C and the combustion gases were measured by a thermal conductivity detector in a CNS elemental analyzer (vario Macro Cube, Elementar, Germany).

Collection and processing of soil microorganism samples

Soil sampling for bacteria (Bacteria), fungi (Fungi) and protists (Protista) was conducted simultaneously with the sampling for the soil physicochemical properties. In each of the six plots per study site, we collected one composite soil sample by mixing 25 soil cores of the topsoil (0–25 cm). For protists, only eight of the 14 locations were sampled. Samples were cooled during transport, stored at 4 °C and homogenised (plant material and debris was removed) and subsampled within 5 days after collection. Prior to the analyses, subsamples were freeze-dried and stored at -20 °C for microbial biomass measurements and -80 °C for DNA isolation at NIOO-KNAW.

Phospholipid fatty-acid (PLFA) and neutral lipid fatty-acid (NLFA) analyses were used to determine bacterial and fungal biomass⁶⁰. Protist biomass was not determined, due to current absence of reliable estimates⁶¹. Phospholipids were extracted twice from 2 g of soil using first 10 and then 5 ml Bligh & Dyer solution (Chloroform:MeOH:citrate buffer 1:2:0.8 v/v/v). Phases were split by dissolving the supernatant in 4 ml chloroform and 4 ml citrate buffer⁶² after which lipids were fractionated based on their polarity (neutral and polar lipids) with SPE Bond Elut 1cc LRC-SI columns (Agilent, USA). Non-polar phospholipids were collected by adding chloroform and polar phospholipids were collected by adding methanol. Both fractions were dried using a RapidVap Vacuum Evaporator System (Labconco, USA) and derivatised into fatty acid methyl esters (FAMEs): neutral-lipid fatty acids (NLFAs) from the chloroform fraction and polar-lipid fatty acids (PLFAs) from the methanol fraction. After addition of an internal standard (FAME standard mix C12:0 sigma 61689 and C19:0 sigma N5377), PLFAs and NLFAs were identified and quantified using a gas-chromatograph with flame ionization detection (Agilent Technologies GC Sampler 80, USA) and their concentrations in µg per gram of freeze-dried soil were calculated. An overview of the fatty acids used as biomarkers for the different microorganismal groups can be found in Table S11.1.

To obtain information on the relative abundance of bacterial, non-arbuscular mycorrhizal fungal (non-AM fungal) and protist taxa in the soil samples, we isolated DNA from 250 mg of soil using the DNeasy PowerSoil Pro kit (Qiagen, USA) following the manufacturer’s standard protocol. AM-fungi were not sequenced for this study. Details of the primers used can be found in Table S11.2. Library preparation, PCR amplification and sequencing were performed by GenomeQuebec (Canada) using the Illumina NextSeq 2000 platform (PE300).

Collection and processing of soil micro- and mesofaunal samples

Nematodes (Nematoda), mites (Sarcoptiformes, Mesostigmata and Trombidiformes) and springtails (Collembola) were sampled in 2021, in eight of the 14 locations. For each of the six plots per study site, nematodes were extracted from a subsample of 100 ml fresh soil ( ~ 80 g) using an Oostenbrink elutriator at NIOO-KNAW and stored in 10 ml 97% ethanol⁶³. The total number of nematodes was counted in a 1 ml subsample using an inverted-light microscope (40x magnification). Counts were expressed as the number of individuals per gram of dry soil. To obtain the relative abundance of nematode taxa in the soil samples, DNA was extracted from a 4 ml subsample using an adapted protocol of the DNeasy PowerSoil kit (Qiagen, USA). DNA was isolated from nematodes in the 97% ethanol solution. After two rounds of centrifuging (1667 rpm and 2307 rpm respectively), the supernatant (ethanol) was discarded. The pellet, including nematodes, was redissolved in C1 buffer of the DNeasy PowerSoil kit after which the manufacturer’s standard protocol was followed. The final solution of DNA material was concentrated to 40 µl with the Concentrator 5301 (Eppendorf, Germany) at 45 °C. See Table S11.2 for an overview of the primers used. Library preparation, PCR amplification and sequencing were performed by GenomeQuebec (Canada) using the Illumina NextSeq 2000 platform (PE300).

To sample mites and springtails, an undisturbed, cylindrical soil sample of 6 cm diameter and 9 cm height was taken in three of the six plots per study site in eight out of the 14 locations. Samples were kept cool and extracted within three days using a Berlese-Tullgren funnel⁶⁴ (Free University of Amsterdam). In closed-off climate chambers, the top of the cores was heated up to 30 °C while the bottom was cooled to 5 °C for three weeks, resulting in soil fauna migrating downwards through the soil core passing a funnel into a 70% ethanol solution, which was afterwards stored at 4 °C. The total number of mites and springtails per sample was counted using a stereomicroscope. Counts were expressed as number of individuals per gram of dry soil. To obtain the relative abundance of springtail and mite taxa in the samples, DNA was extracted from the animals stored in ethanol using an adapted protocol of the DNeasy PowerSoil Pro kit. After counting, samples were transferred to 50 ml tubes and oven-dried (60 °C) to evaporate all ethanol and water to optimise the DNA isolation steps. As a pre-treatment, before using the Powersoil Pro Kit, a hard lysis was performed by adding 3 g of beads (Zirconia 1 mm, Biospec products, USA) and 3 ml of CD1 buffer and vortexing at maximum speed for 10 min. Afterwards, samples were equally divided in four 2.5 ml tubes, whereafter 18.75 μl SDS (20%) and 15 μl proteinase K (20 mg/ml) were added to each of the four subsamples. Samples were then incubated at 900 rpm and 56 °C for 2 h. Afterwards, the four subsamples were again combined (taking 150 µl of each) and the steps of the Powersoil Pro Kit were continued. Library preparation, PCR amplification and sequencing were performed by GenomeQuebec (Canada) using the Illumina NextSeq 2000 platform (PE300).

Collection and processing of soil macrofaunal samples

In 13 locations (14 of the 15 food forest sites and their respective references), litter-dwelling macroarthropods were collected using two pitfall traps per study site that were set up in four sampling sessions. Using pitfall trap sampling, the number of captured individuals is a composite measure of their density and activity in the area⁶⁵. These pitfall traps were located within 10 m from the centre of each half of the study site, at a location that was considered most representative for this area of the study site in terms of vegetation. In 21 of the 324 observations, the pitfall traps were placed closer to the edge of the field to avoid disturbance by grazers or inconvenience for the managers. Pitfall traps were 8.5 cm in diameter and 11 cm deep and contained a 50% ethylene glycol/water solution and a drop of detergent to reduce water surface tension. The traps were covered by aluminum roofs, leaving a gap of about 3 cm to allow arthropods to enter while sheltering the traps from precipitation, litter and direct sunlight¹⁶. Each sampling session consisted of 14 days. The pitfall traps were present in the study sites in early April, early May, early July and in late September. The collected arthropods retained on a sieve with mesh width of 0.5 mm were stored in 70% ethanol. Isopods (Isopoda), millipedes (Diplopoda), centipedes (Chilopoda), harvestmen (Opiliones) and carabid beetles (Carabidae) were sorted at the ForNaLab and sent to taxonomic experts for identification up to species level. Some pitfall traps were unusable in some of the sampling sessions (in 28 of the 352 cases, or 8%, especially in the fourth sampling session). Therefore, when combining the data of the different sampling sessions (four sessions in the case of harvestmen, only the three first sessions for the other groups), we averaged the number of captured individuals of a given species across the sampling sessions rather than summing them.

In all 14 locations (15 food forests and their respective references), earthworms (Lumbricidae) were collected in three of the six plots. Earthworms were collected in October and November of 2022 when soil moisture was relatively high and temperatures were above freezing point. Within 0.1 m² squares, litter was removed and hand-sorted for earthworms. Subsequently, the soil within this 0.1 m² was excavated to a depth of 20 cm and hand-sorted for earthworms for 50 min (30 min by a first observer, 20 min by another observer). Parallel to the hand-sorting, up to 9 l of a mustard solution (6 g mustard powder in 1 l water) was poured into the pit to bring deep burrowing earthworms to the surface, which were collected during the 50 min hand-sorting⁶⁶. All earthworms were fixed on 5% formaldehyde and further preserved in a 70% ethanol solution. Adult earthworms were identified to species level and counted at ForNaLab. Juvenile earthworms and juvenile harvestmen that could not be identified to species level were included in the total number of individuals but excluded from the analyses of diversity and community composition.

Data processing

All data processing and analysis was performed in R V4.4.2⁶⁷. All packages used for data processing, analysis and visualisation are listed and referenced in S12.

For the sequencing data, demultiplexing and removing of adapters was performed by GenomeQuebec. Further processing of raw reads of 16S, 18S and ITS primer pairs was done with dada2. Filter and trim settings for each primer pair can be found in Table S13.1. Primers were identified and removed using cutadapt. After dereplicating unique sequences, sample inference was done using adjusted code to estimate the sequencing error rate of the datasets, with altered loess arguments (weights and span) and enforcing monotonicity because the data was obtained using the Nextseq 2000 platform. Thereafter, paired ends were merged and chimeras were removed using the consensus method. For the 18S primer pair to identify the protist community only the forward primer was used, because of insufficient overlap of paired sequence ends. The SILVA database V138.1 was used to assign taxonomy for the 16S reads up to genus level⁶⁸. The UNITE V9.0 database was used for the ITS reads assigning up to species level⁶⁹. The PR2 database V5.0.0 was used for the protist and assigned up to species level⁷⁰ and the NemaTaxa V1.0 database was used for the nematode 18S reads and assigned ASVs up to genus level⁷¹. Following Nardi et al. (2023)⁷², for the microarthropods, the ASVs from the dada2 pipelines were first clustered into OTUs with a 97% similarity threshold prior to assigning taxonomy using VSEARCH. Assigning taxonomy was done using taxalogue to combine the BOLD V4⁷³, GBOL⁷⁴ and MIDORI GB260⁷⁵ databases. An overview of the total number of reads throughout the sequence data processing can be found in Table S13.1. We used phyloseq for data exploration and filtering of non-target taxa and to create ASV/OTU abundance matrices used for further analyses.

Alpha- and beta-diversity measures

As measures of the alpha-diversity (at the species or ASV/OTU (hereafter referred to as ‘taxon’) level) of the sample communities, we computed the Hill numbers ⁰D (species/taxon richness) and ¹D (the effective number of species/taxa given by the exponent of Shannon entropy ((H)) with (H=-mathop{sum }limits_{i=1}^{R}{p}_{i}mathrm{ln}left({p}_{i}right),) where R is the total number of species and p_i is the proportion of individuals belonging to the i-th species)⁷⁶. For the sequencing data, we used iNext’s estimateD to compute the mean ⁰D and ¹D, using individual-based rarefaction to control for uneven library sizes across samples⁷⁷. For this, we removed all samples with a read depth below a certain threshold (see Table S13.2 and Fig. S13). For each of the remaining samples, we used 50 iterations of random subsampling to this read depth, calculating the diversity metric for each of those 50 random subsamples. We used the average of these 50 values as our estimate of the diversity metric for the respective sample to mitigate the risk of excluding or oversampling rare taxa⁷⁷. For non-sequencing data, we used hillR’s hill_taxa to compute the alpha-diversity measures based on the raw count data. Beta-diversity was calculated at the level of a study site as the average pairwise Bray-Curtis dissimilarity between plots within this study site⁷⁸, using vegan’s avgdist for the sequencing data and vegdist for the non-sequencing data.

Data analysis

We compared key soil physicochemical properties (bulk density, pH, total C, total N and total P) and soil communities’ biomass or number of individuals, alpha-diversity (⁰D and ¹D), community composition and beta-diversity between food forests and each of the reference land use types. For the biomass or number of individuals and alpha-diversity measures, this comparison was conducted for each of the 12 groups separately. For the analysis of the community composition and beta-diversity, we pooled the two microarthropod groups (springtails and mites), as well as the five macroarthropod groups (isopods, millipedes, centipedes, harvestmen and carabids), because of the low number of reads or counts for the separate groups. We only included communities that contained at least five individuals for these analyses⁷⁹. For each of the comparisons, we tested the null hypotheses that there was no difference between the food forests and forests, between the food forests and grasslands and between the food forests and croplands.

To test the null hypotheses for the soil physicochemical properties, biomass or number of individuals and alpha-diversity, we formulated generalised linear mixed models (GLMMs) with land use type as a fixed effect and the location of the study sites as a random effect in glmmTMB. For the number of springtail individuals and the effective number of species/taxa of protists and mites the data did not allow for a reliable estimation of the location random effect variance, so we used generalised linear models (GLMs) with land use type as the sole predictor. For the beta-diversity, we also used GLMs with land use as the sole predictor and excluded croplands from the analysis for protists, nematodes, microarthropods and earthworms, because they had less than six observations. We modelled the soil physicochemical properties, biomass, number of individuals and effective number of species with a Tweedie distribution when the data contained zeros and a gamma distribution with log link function when it did not. The species/taxon richness was modelled with a gamma distribution for the sequencing data and a Poisson distribution for the non-sequencing data. We modelled the beta-diversity with a beta distribution after transformation of the values to y’ = (y(n-1) + 0.5)/n, with n being the total number of observations, to fit the (0,1) interval⁸⁰. Model formulae and additional information on the fitted models is given in S1. There were no strong deviations from model assumptions for most models, as evaluated by the diagnostics on the simulated residuals using the DHARMa package. For the models of the effective number of species/taxa of the protists, isopods, earthworms and especially centipedes, there was some underdispersion (dispersion as determined by DHARMa’s testDispersion 0.50 to 0.63), leading to a conservative bias. To estimate the variance of the model estimated means to obtain 95% confidence intervals (CIs) and p-values, we used parametric bootstrapping with 5000 simulations using lme4’s bootMer. Estimated means were back-transformed from the log-scale with a bias adjustment for the random effect variance to report means per land use type and the ratios thereof. Asymptotic Wald confidence intervals were used for these ratios. We used the Nakagawa-R²⁸¹ or Nagelkerke-R² as a coefficient of determination for the GLMMs and GLMs respectively, computed using performance. Variation in key soil properties among study sites was visually summarised using a principal component analysis.

To test the null hypotheses for the community composition, we used PERMANOVA analysis with adonis2 in vegan. Prior to the PERMANOVA analysis, we conducted a PERMDISP with betadisper in vegan to evaluate whether we could reliably use the PERMANOVA analysis (S5)⁸². For the comparison of the non-AM fungi, protists, nematodes and earthworms between food forests and croplands, the dispersion was statistically significantly (p < 0.05) lower in croplands, which also had a lower number of observations. This leads to some conservative bias, but we considered the difference in dispersion (summarised as the average distance from the centroid) sufficiently limited to proceed with the PERMANOVA (S4). To visualise how the community composition of each group compared between the food forests and reference land use types, we used Nonmetric Multi-Dimensional Scaling, reducing the dimensionality to two dimensions, using the metaMDS function in vegan. For the taxa captured in the same sample plots where soil properties were measured (bacteria, non-AM fungi, protists, nematodes, microarthropods and earthworms), we additionally explored how variation in key soil properties (bulk density, pH, total C, total N and total P) related to the variation in community composition represented by the ordination. We did so by modelling these relationships with generalised additive mixed models using mgcv and visualising the obtained relationships in contour plots (S7).

For all tests, p-values were corrected for multiple testing using a Holm-adjustment. We considered differences between land use types statistically significant when p < 0.05.