# Analysis of complex trophic networks reveals the signature of land-use intensification on soil communities in agroecosystems

Juliette Bloor, Sara Si-Moussi, Pierre Taberlet, Pascal Carrère, Mickael Hedde

► To cite this version:

Juliette Bloor, Sara Si-Moussi, Pierre Taberlet, Pascal Carrère, Mickael Hedde. Analysis of complex trophic networks reveals the signature of land-use intensification on soil communities in agroecosystems. Scientific Reports, 2021, 11 (1), 10.1038/s41598-021-97300-9. hal-03345139

HAL Id: hal-03345139

https://hal.inrae.fr/hal-03345139

Submitted on 26 Apr 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

OPEN

# Analysis of complex trophic networks reveals the signature of land-use intensification on soil communities in agroecosystems

Juliette M. G. Bloor<fnref n="1" />,*, Sara Si-Moussi²,³,⁵, Pierre Taberlet³,⁴, Pascal Carrère<fnref n="1" /> & Mickaël Hedde²

<figure><img src="image_3.png" /></figure>

Increasing evidence suggests that agricultural intensification is a threat to many groups of soil biota, but how the impacts of land-use intensity on soil organisms translate into changes in comprehensive soil interaction networks remains unclear. Here for the first time, we use environmental DNA to examine total soil multi-trophic diversity and food web structure for temperate agroecosystems along a gradient of land-use intensity. We tested for response patterns in key properties of the soil food webs in sixteen fields ranging from arable crops to grazed permanent grasslands as part of a long-term management experiment. We found that agricultural intensification drives reductions in trophic group diversity, although taxa richness remained unchanged. Intensification generally reduced the complexity and connectance of soil interaction networks and induced consistent changes in energy pathways, but the magnitude of management-induced changes depended on the variable considered. Average path length (an indicator of food web redundancy and resilience) did not respond to our management intensity gradient. Moreover, turnover of network structure showed little response to increasing management intensity. Our data demonstrates the importance of considering different facets of trophic networks for a clearer understanding of agriculture-biodiversity relationships, with implications for nature-based solutions and sustainable agriculture.

Soil is a key reservoir of biological diversity on Earth, supporting complex food webs and biotic interactions that underlie biogeochemical cycling and the provision of fundamental ecosystem services such as plant production, carbon storage and the biological regulation of pest species<fnref n="1" />,². Indeed, biodiverse soils are considered to play a pivotal role in sustainable agricultural systems and food security<fnref n="1" />,³. In recent years, widespread concerns over global biodiversity loss in aboveground organisms have heightened the interest in soil biodiversity patterns in response to global change⁴⁻⁶. Nevertheless, holistic understanding of the impacts of global change pressures on the diversity and interactions within complex networks of soil organisms remains limited⁷.

Soil systems are particularly sensitive to changes in land-use and management practices which modify soil physico-chemical properties⁶,⁸. Agricultural intensification can affect soil biodiversity both directly (through the application of pesticides, fertilizers and tillage) and indirectly, via changes in plant diversity and production which affect inputs of organic matter to the soil⁵,⁹,¹⁰. A large body of research suggests that intensive agricultural management practices and mineral fertilizer inputs can have negative effects on the species richness of many soil functional groups including earthworms, springtails, nematodes and oribatid mites⁹,¹¹⁻¹³. Increasing land-use intensity has also been associated with decreases in microbial biomass¹⁰,¹⁴ and/or shifts in the relative abundance of different microbial groups (bacterial, archaeal and fungal taxa)¹⁵,¹⁶. For example, increased inorganic nitrogen inputs and organic matter with a low C:N ratio are commonly reported to reduce the fungal:bacteria ratio in intensively-managed agroecosystems¹⁰,¹⁶. Such responses to agricultural intensification at the species- or functional group-level can have implications for soil food web structure and trophic networks by modifying trophic resource availability, and hence the potential for trophic interactions and the prevalence of highly-connected taxa. Indeed, recent evidence suggests that intensive land-use management decreases the complexity of soil bacterial

<footnote n="1">¹Université Clermont Auvergne, INRAE, VetAgro-Sup, UREP, Clermont-Ferrand, France. ²Eco&Sols, Université Montpellier, CIRAD, INRAE, Institut Agro, IRD, Montpellier, France. ³Laboratoire d'Ecologie Alpine (LECA), CNRS, Université Grenoble Alpes, Grenoble, France. ⁴UiT – The Arctic University of Norway, Tromsø Museum, Tromsø, Norway. ⁵Laboratoire TIMC-IMAG, CNRS, Grenoble INP, Université Grenoble Alpes, Grenoble, France. ⁶email: juliette.bloor@inrae.fr</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Inference of trophic network structure for the present study using Next-Generation Sequencing technologies and a database of ecological interactions.</figcaption></figure>

networks¹⁷, as well as that of soil networks based on a relatively-limited number of micro-/mesofaunal groups⁸,⁹. However, agricultural intensification is not consistently harmful to soil fauna diversity¹⁶,¹⁸,¹⁹, and variation or non-linearity in species responses may therefore lead to idiosyncratic effects on complete food web structure. Given that loss of interactions among species may also supersede the effects of species loss within trophic groups on ecosystem function²⁰, there is a pressing need for data on the responses of complex soil interaction networks to land-use intensification.

To date, the construction of comprehensive food webs and entire interaction networks for soil communities has encountered a number of technical difficulties. As soil organisms show a 10⁶-fold difference in both body size and abundance², traditional taxonomic sampling methods show considerable variation in spatio-temporal resolution and are difficult to standardize across groups. In addition, simple co-occurrence analysis of soil organisms may not accurately reflect or represent trophic interactions²¹. Consequently the construction of interaction networks requires detailed information on the nature of interactions and predator–prey relationships between different taxonomic groups which may be difficult to acquire for complex communities²². Over the last decade, environmental DNA and metabarcoding technology (DNA taxonomy combined with high-throughput DNA sequencing) have emerged as a promising tool for soil biodiversity monitoring²³⁻²⁵. Environmental DNA (eDNA) refers to the DNA present in environmental samples (such as soil, water or air), and includes a mixture of genomic DNA from both living cells/organisms (e.g. from skin, mucus or secretions) as well as extracellular DNA linked to cell death and degradation. Compared to the traditional description of taxonomic-based trophic networks, eDNA allows a faster and efficient construction of more comprehensive networks based on presence data²³,²⁶. Moreover, the identification and assignment of trophic groups in complex networks now benefits from the increasing availability of open knowledge bases on biotic interactions²⁷ and functional guilds²⁸. It is therefore possible to combine inherently-qualitative eDNA data with information on trophic interactions, and overcome some of the traditional limitations of soil food web analysis. However, although eDNA metabarcoding has clear potential to provide novel insights into the drivers of changes in soil biotic communities, this methodology has rarely been tested in the context of agricultural intensification²⁹.

Here we use eDNA metabarcoding, data mining on trophic interactions and machine-learning in order to infer comprehensive, interaction-based soil food webs (Fig. 1), and to examine the trophic structure of soil communities for sixteen agroecosystems along a gradient of agricultural pressure (from arable crops to grazed permanent grasslands) in a long-term management trial. Analyses of food web properties and trophic networks were used to assess the potential for different network metrics to discriminate between agricultural practices, and to determine the sensitivity of trophic networks to land-use intensity. Our overarching hypothesis was that land-use intensification would decrease the trophic complexity of soil communities. Specifically, we hypothesized that increasing land-use intensity would reduce trophic network size, modifying soil food web connectance (number of realized interactions) and average centrality (prevalence of highly-connected species). In addition, we examined impacts of increasing land-use intensity on energy-flow pathways within trophic networks, predicting an increase in the dominance of bacteria over fungi but a decrease in the importance of parasitic interactions (the latter response due to a decrease in trophic resources). Finally, we explored the linkages between multi-trophic β-diversity (i.e. turnover in trophic group networks between fields) and land-use intensification, in order to test for possible divergence between the responses of local assemblages (α-diversity) and compositional turnover
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Relationships between land-use intensity and average trophic network properties across experimental fields. Regression lines and 95% CI are shown where significant. Network metrics are given in Table 2. The colour code emphasizes the land-use intensity gradient from least intensive (purple) to most intensive (yellow).</figcaption></figure>

to agricultural pressure. Previous work on a limited number of soil trophic groups has suggested that this phenomenon occurs in temperate grassland communities along a land-use intensification gradient¹⁹. Such empirical evidence could demonstrate the role of biotic interactions in determining the outcome of ecosystem responses to land-use intensification, and highlight which network metrics may be useful indicators of soil biodiversity in managed ecosystems.

## Results

**Taxonomic composition and trophic groups.** Metabarcoding data yielded a final database (and metabweb) comprised of 21,667 pairwise interactions (6% parasitic, 0.5% symbiotic), covering 99.7% of the taxa recorded (total of 1049 families/species/genera). 981 taxa were retained after filtering marine organisms and taxa from non-temperate terrestrial ecosystems (Supplementary Data table). Assessment of the taxonomic diversity among all fields showed that 14%, 48% and 38% of OTUs were assigned to family, genus and species level respectively. Across all fields, five out of 44 phyla represented about 50% of all taxa. The Proteobacteria were the most represented phylum (17.1% of the taxonomically-identified sequences), followed by Arthropoda, Actinobacteria, Nematoda and Ciliophora comprising 9.1%, 7.7%, 7.7% and 6.5% of all taxa respectively. On average we detected 369 taxa in each experimental field. Overall, 79.1% of phyla and 36% of families were found to be common across all experimental fields. Trophic inference generated a total of 16 trophic groups across experimental fields (range 12–16 per field, Supplementary Fig. S1).

**Food web responses to land-use intensification.** Assessment of food web properties indicated considerable variation in diversity and structure across experimental fields (Fig. 2). Eight out of the 12 network
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Principal Coordinates Analysis based on food web dissimilarities for interaction networks constructed with Operational Taxonomic Units (OTUs) in experimental fields. Colours illustrate the gradient of land-use intensity (LUI) across experimental fields. Black diamonds represent the barycenter of each study sites.</figcaption></figure>

indices representing α-diversity and network topology showed significant linear responses to agricultural intensification gradient (Fig. 2, Supplementary Table S1). Total number of trophic groups and trophic group equitability (entropy) declined along the gradient of intensification, as did average node degree, omnivory level, link density and prevalence of parasitic links. The magnitude of decrease in response to increasing land-use intensity was greatest for average degree, trophic group entropy, network link density and parasitic links (standardized effect sizes, SES, ranging from -0.15 to -0.22); decreases in trophic group richness and omnivory index were more limited (SES -0.03 and -0.10 respectively). In contrast, both the bacteria-to-fungi path ratio and the detritivory-to-root herbivory path ratio increased along the gradient of intensification, with a two-fold increase from values in 'low-intensity' fields to those at the 'high-intensity' end of the gradient (SES of 0.16 and 0.09 for the bacteria-to-fungi path ratio and the detritivory-to-root herbivory path ratio respectively). Total number of taxa (node richness), mean trophic level, mean and maximum path length showed no response to the agricultural intensification gradient (Fig. 2).

Analysis of network dissimilarity across experimental fields indicated a significant turnover in node composition between fields (DF = 15, F = 10.08, *p* < 0.01) and sites (DF = 3, F = 3.33, *p* < 0.01) for taxa-based food webs (Fig. 3), but no significant turnover in trophic group networks. Permanent grassland fields at the low end of the land-use intensity gradient were characterized by greater β-diversity compared to the more-intensively managed fields at the Mons and Lusignan sites (Fig. 3), but network dissimilarity showed no clear relationship with land use intensity (LUI). Site rankings along PCoA1 were negatively correlated to soil pH across experimental fields (*r* = -0.96, *p* < 0.05), but were unrelated to any other pedoclimatic variables.

## Discussion

The construction of data-based food webs including multiple interaction types is essential for the accurate appraisal of ecosystem vulnerability to global change<sup>23,30,31</sup>. Recent work on complex aboveground trophic networks has shown that agricultural intensification drives changes in network structure and robustness<sup>32</sup>, but studies of complex soil networks are lacking and current understanding of soil network responses to land-use is based on a relatively-limited number of taxonomic groups and network metrics<sup>8,9</sup>. Our study of complex soil food webs generated by eDNA metabarcoding in 16 agricultural fields along a broad LUI gradient resulted in three major findings: (1) trophic interaction networks show a range of consistent qualitative responses to agricultural intensification; (2) soil network complexity is more sensitive to agricultural intensification than network size; (3) local assemblages (metrics of α-diversity) show stronger statistical relationships with agricultural intensification than does composition turnover (β-diversity). Together our results bring new insights into the importance of different facets of food web structure for biodiversity preservation, and the contribution of biological interactions to agriculture-biodiversity relationships.

Five out of nine network metrics describing node, community and network properties showed significant negative relationships with land-use intensity in our study, providing support for the hypothesis that increasing management intensity decreases trophic diversity and the complexity of soil interaction networks. In general, we found that the magnitude of negative effects induced by our LUI gradient was greater for soil network complexity
than network size; agricultural intensification was associated with markedly-simpler networks (lower interactions and distribution of interactions between nodes) compared to relatively-limited changes observed in trophic group richness. Such a loss of network complexity and cohesion has previously only been demonstrated for soil microbial networks in relation to cropping intensity<sup>17,33</sup>, although our observed negative relationship between link density and LUI also matches results from comparisons between soil networks in grasslands and arable crop systems<sup>8</sup>. Reduced diversity of interactions may stem from less structurally-complex habitats, and shifts in environmental and/or biotic filtering under intensive management practices<sup>34,35</sup>. Decreases in trophic network complexity are considered to have repercussions for community stability and soil functioning in a changing environment<sup>31</sup>, and recent work with combined above- and belowground trophic groups has also demonstrated linkages between network complexity and the provision of ecosystem services<sup>36</sup>.

We found greater effects of LUI on trophic group richness/equitability compared with taxa richness, in line with the idea that species number is less critical for soil ecosystem function than functional diversity due to functional redundancy in soils<sup>37</sup>. Where newly-occurring or “lost” species exhibit unique functional attributes, changes in functional richness can markedly exceed the change in taxonomic diversity. In systems with relatively low functional redundancy, increases in LUI may therefore cause strong functional shifts. Limited responses of taxa richness to LUI observed in the present work could partly reflect the dominance of smaller-bodied organisms in our food webs, since previous studies have suggested that larger soil organisms and top predators show greater sensitivity to agricultural intensification compared with smaller-sized soil biota<sup>7</sup>. Low rates of occurrence of large top predators in our dataset may also contribute to the limited responses of average path length and mean trophic level to agricultural intensification in our soil networks; additional studies with biodiversity ‘soups’ and complementary sampling methods could help test the generality of our findings<sup>38</sup>.

Energy-flow pathways within trophic networks reflect the quantity and quality of resources, stoichiometric constraints and 'top-down' regulatory processes<sup>39,40</sup>. In the present study, increasing management intensity increased the relative importance of both the bacterial energy channel compared to the fungal channel, and the detrital energy channel compared to the root herbivory channel. These results are consistent with changes in plant biomass allocation and rates of decomposition, root exudation and nitrogen mineralization under agricultural intensification, which promote the detrital energy channel<sup>39,41</sup>, but see<sup>10</sup>. Moreover, the 'slow cycle' fungal energy channel tends to be more dominant in extensively-managed systems where the ratio of carbon to nitrogen (C:N) in organic matter is high. Although some interconnections between energy channels exist<sup>42</sup>, dominance of fungal, bacterial or root energy channels in soil trophic networks has significant implications for subsequent nutrient turnover rates and ecosystem stability<sup>39</sup>. Our results highlighted a decrease in the prevalence of parasitic interactions with increasing agricultural pressure. This response pattern likely reflects a decrease in trophic resources (vector abundance) due to the increased use of pesticides and fertilizers along our LUI gradient. Data on soil parasites and agricultural intensification is scarce in the literature, but decreases in the prevalence of parasitic nematodes and mites have previously been observed following the conversion of grassland to arable land<sup>12</sup>.

In contrast to our metrics of network size and structure or energy channels, turnover in trophic networks (β-diversity of nodes) showed no clear pattern with increasing LUI across our study systems. Moreover, we found no significant turnover in trophic groups, supporting the idea of a common “backbone” of interactions across food webs<sup>43</sup>. These results suggest that agricultural intensification has divergent effects on local interaction networks and the network turnover of soil food webs across different sites, with greater effects of LUI on α- rather than β-diversity. Indeed, the marked responses of local assemblages (α-diversity) to agricultural intensification in the present study could have served as a stabilizing mechanism for broader-scale patterns of diversity<sup>44</sup>. Limited β-diversity responses and absence of biotic homogenization to LUI may also be linked to confounding effects of pedoclimatic conditions<sup>15,45</sup>, or to functional redundancy and phylogenetic niche conservatism<sup>46</sup>. Irrespective of the underlying mechanisms, insensitivity of beta diversity to increasing agricultural pressures may promote the buffering capacity of soil food webs against land-use change, with cascading effects on the connections between detritus-based and primary production-based food webs, and the regulation of aboveground and belowground dynamics<sup>40</sup>.

With regards to our methodological approach, eDNA metabarcoding can meet challenges faced by traditional taxonomic monitoring, use high taxonomic resolution to confront current paradigms of ecosystem dynamics, and provide a valuable tool to support evidence-based decision making<sup>47</sup>. In the present study we show that eDNA metabarcoding is an effective technique for visualizing complex soil communities, and for analysing the structural and functional relationships within soil food webs. To date, very few studies on soil food web structure integrating microbes, micro- and mesofauna have documented large numbers of taxa (> 500), and those rare examples involved spatially-distinct and labour-intensive diversity protocols<sup>8</sup>. Not only were we able to identify a very high number of taxa from microbes and key micro- and meso-faunal groups using a simple field protocol, but the analysis of standard soil samples also provided good assurance of species/ taxa co-occurrence in space.

Unlike bulk-community DNA sampling which targets whole organisms, eDNA samples detect biological signals in the environment from traces of intracellular and extracellular DNA. Given that DNA fragments may persist for periods ranging from several days to several months in the environment, this increases the temporal uncertainty associated with the diversity assessment<sup>47</sup>. It is therefore possible that taxa detected in soil samples in the present study were not all present at exactly the same time, and hence our food webs represent the potential complexity of trophic interactions in our experimental fields, rather than a snapshot of trophic interactions at a single date. Further work using additional sampling dates, and combined bulk- and eDNA sampling, would provide a clearer picture of the degree of temporal uncertainty associated with eDNA metabarcoding results at our study sites.

Finally, quantifying land-use intensity is a non-trivial issue in studies of agricultural intensification, since multiple processes relating to inputs, outputs and mechanical disturbance need to be considered<sup>48,49</sup>. As yet, there is no widely-accepted, purely quantitative measure for ranking a set of agricultural fields along a gradient
<table><caption>**Table 1.** Soil type, management and land-use intensity (LUI) ranking for the sixteen experimental fields in France. Values for LUI are based on a multicriteria evaluation (Supplementary Fig. S2); more positive values indicate more intensive management practices. Mean annual temperature and rainfall of sites as follows: Mons: 11.7 °C, 614 mm; Lusignan: 11.7 °C, 800 mm; Theix: 8.7 °C, 780 mm; Laqueuille: 8 °C, 1000 mm.</caption><thead><tr><th>Site</th><th>Soil type</th><th>Agroecosystem</th><th>Management</th><th>LUI</th></tr></thead><tbody><tr><td rowspan="4">Mons</td><td rowspan="4">Silty loam (Luvisol)</td><td rowspan="4">Arable crops (6-years rotation: spring pea; winter wheat; rapeseed; spring barley; maize; winter wheat)</td><td>Conventional tillage, reference N inputs</td><td>1.70</td></tr><tr><td>Conventional tillage, reduced N inputs</td><td>1.46</td></tr><tr><td>Reduced tillage</td><td>-0.11</td></tr><tr><td>Reduced tillage, residue removal</td><td>0.26</td></tr><tr><td rowspan="4">Lusignan</td><td rowspan="4">Loamy clay (Cambisol)</td><td>Arable crops (3-years rotation: maize; wheat; barley)</td><td>Conventional tillage, reference N inputs</td><td>2.14</td></tr><tr><td>Crop-grass rotation (6-years rotation: maize; wheat; barley; 3-years sown grassland)</td><td>Tillage for crops, mown grassland, reference N inputs</td><td>1.52</td></tr><tr><td rowspan="2">Grass-crop rotation (9-years rotation: 6-years sown grassland; maize; wheat; barley)</td><td>Tillage for crops, mown grassland, reference N inputs (2 fields, one with additional pesticide treatments)</td><td>0.98</td></tr><tr><td>Tillage for crops, mown grassland, reduced N inputs</td><td>0.87</td></tr><tr><td rowspan="2">Theix</td><td rowspan="2">Sandy clay loam (Cambisol)</td><td rowspan="2">Mown permanent grassland</td><td>Tillage for crops, Grazed grassland</td><td>0.57</td></tr><tr><td>Reference NPK inputs</td><td>-1.05</td></tr><tr><td rowspan="4">Laqueuille</td><td rowspan="4">Clay loam (Cambisol)</td><td rowspan="4">Mown permanent grassland</td><td>No fertilizer inputs</td><td>-1.91</td></tr><tr><td>Reference NPK inputs</td><td>-1.13</td></tr><tr><td>No fertilizer inputs</td><td>-1.96</td></tr><tr><td>High stocking rate, Reference N inputs</td><td>-1.86</td></tr><tr><td rowspan="2">Silt loam (Andosol)</td><td rowspan="2">Grazed permanent grassland</td><td>High stocking rate, Reference N inputs</td><td>-1.86</td></tr><tr><td>Low stocking rate, no N inputs</td><td>-2.08</td></tr></tbody></table>

of management intensity, although indices of LUI do exist for functionally-similar agroecosystems such as grasslands¹⁹,⁴⁸ or arable cropping systems³³. In the present study, we used an unweighted, multicriteria approach⁴⁸ and calculated a continuous LUI index based on management practices over the preceding 6-year period. Our approach and choice of criteria allowed the positioning of very different agroecosystems (permanent grasslands, annual cropping systems) along a standardised gradient, and the site ranking obtained was subsequently validated by expert opinion (ANAEE-F ACBB site managers). Our calculation of the intensification gradient also provided an objective method to quantify the differences in agricultural intensification between sites, which is not possible to achieve based on expert knowledge alone. Of course, the absolute values of such an index will change depending on the exact criteria used and the cohort of sites examined, making direct comparisons between studies difficult. The development and validation of a standardised, additive index for use across the entire range of cropping systems and managed grasslands would facilitate the identification of broad-scale response patterns and assessment of the relative importance of drivers of agricultural intensification in future studies.

## Conclusion

Together, the results of this study extend previous work on multiple taxa and land-use, and highlight the need to consider the response patterns of complex trophic networks for a clearer understanding of linkages between biodiversity and agricultural practices. Our results indicate that agricultural intensification has profound impacts on trophic interactions and soil trophic network structure, driving decreases in trophic group diversity (richness, equitability) and complexity of interactions (average degree, link density). Observed changes in energy flow pathways were consistent with shifts in resource availability (labile C, living roots) and increases in pesticide use due to agricultural intensification. In contrast, eDNA metabarcoding did not highlight sensitivity to agricultural intensification for either taxa richness, maximum path length or mean trophic level. Soil trophic networks in our studied agroecosystems showed divergent responses of α- and β-diversity to intensification, and we also found evidence for a common backbone of trophic group interactions across food webs, with implications for both ecological assembly rules and diversity conservation strategies. Although more broad-scale studies which include the temporal dynamics of complex soil food webs are needed, the present results reveal that agricultural intensification simplifies entire soil trophic networks, with likely consequences for biodiversity-ecosystem function relationships.

## Methods

**Field sites and sampling.** Sixteen experimental fields were selected from within four sites of the long-term ANAEE-F ACBB agroecosystems management trial in France (http://www.soere-acbb.com/): Mons (49° 52' N, 3° 1' E, 85 m a.s.l.), Lusignan (46° 25' N, 0° 7' E, 151 m a.s.l.), Theix (45° 43' N, 3° 1' E, 880 m a.s.l.) and Laqueuille (45° 38' N, 2° 44' E, 1100 m a.s.l.). These four temperate sites cover a range of soil types and climate conditions, and represent a gradient of management intensity from permanent grasslands to arable cropping systems (Table 1). At each site, we selected fields in order to generate a broad gradient of management intensity across sites (Table 1); at the time of the study, fields had been subjected to experimental treatments for over 10 years. The gradient of management intensity for the 16 fields was quantified using a multi-criteria evaluation based on soil tillage (depth, frequency), N inputs (mineral N, N fixation), C exports (biomass harvests,
residue management) and frequency of crop protection treatments for the period 2010–2016; fields were ranked according to their scores from the first axis of a principal components analysis using the management criteria (Supplementary Fig. S2).

Soil sampling was carried out at a standardized plant phenological stage across sites i.e. the peak of the plant vegetative growing period in 2016 (mid-May to mid-June depending on the plant phenology at each site). Four plots of 0.25 m² were selected for sampling in each experimental field; plots were randomly-positioned in each quadrant of the field, at least 10 m apart and at least 5 m away from the edge of the field. One soil core (10 cm depth, 8 cm diameter) was taken from each plot and coarsely-sieved to remove stones/large debris (heat-sterilized sieve, 4 mm mesh size), providing a total of 64 soil cores. A 15 g sample of freshly-sieved soil from each soil core was placed in a sterile plastic container and preserved with silica gel within 30 min of field collection to prevent microbial growth prior to DNA extraction⁵⁰.

**Molecular analysis and sequence curation.** Soil biodiversity was assessed with DNA metabarcoding using eight DNA markers. Three markers targeted regions of the ssu rRNA gene in the Bacteria and Eukaryota domains, whereas five additional markers focused on earthworms, collembola, oligochaetes, insects and vascular plants (primer details for Bact02, Coll01, Euk02, Euk03, Inse01, Lumb01, Oligo01, Sper01 given in Supplementary Table S2). Extracellular DNA was extracted from 15 g of soil following methods described previously⁵¹. Briefly, soil was extracted with a saturated phosphate buffer (Na₂HPO₄; 0.12 M; pH 8); sub-samples of the slurry were centrifuged, and the resulting supernatant was used as starting material for the NucleoSpin Soil kit (Macherey-Nagel, Düren, Germany), following manufacturer's instructions but skipping the lysis step. The resulting DNA extracts were recovered in 100 µL and diluted ten times. Each diluted extract was then PCR-amplified, with all the above primer pairs. All PCR were carried out in a final volume of 20 µL containing 2 µL of DNA extract. The amplification mixture consisted of 10µL of AmpliTaq Gold 360 master mix (Applied Biosystems), 0.5 µM of each primer and 0.16 µL (20 mg/mL) of bovine serum albumin (BSA, Roche Diagnostic). Polymerase activation was performed at 95 °C for 10 min, followed by 32–45 cycles at 95 °C for 30 s (denaturation), 45–55 °C for 30 s (primer annealing) and 72 °C for 60–90 s (extension), followed by a final elongation for 7 min at 72 °C.

We carried-out four technical PCR replicates for each sample and for each primer pair. In order to minimize possible biases in the experimental workflow, each technical replicate for a primer pair included three extraction controls, 12 blanks (no primer, no template), nine PCR negative controls (ultrapure water), and eight positive controls (positive controls only for Sper01 and Euka02 primer pairs). PCR products were purified using the MinElute PCR purification kit (Qiagen GmbH). Library preparation was performed using the MetaFast protocol by Fasteris (https://www.fasteris.com/dna/?q=content/metastar-protocol-amplicon-metagenomic-analysis), which significantly limits the tag-jump problem⁵⁰. For Bact02, sequencing was performed by 2×250-bp paired-end sequencing on the Illumina MiSeq platform, while for all other primer pairs sequencing was performed by 2×125-bp paired-end sequencing on the Illumina HiSeq 2500 platform using default settings at Fasteris. Sequence data were processed using OBITools software⁵² to (1) assemble and dereplicate reads, (2) match sequences to the original samples, (3) remove noise from the data by removing singletons, low-quality sequences, putative PCR/sequencing artefacts (criteria used to remove low-quality reads⁵³), and (4) taxonomically-assign the remaining sequences. At the end of the data curation process we obtained a total of 5880 robust OTUs (Operational Taxonomic Units).

**Construction of trophic interaction networks and inference of trophic groups.** Heuristic food webs were constructed from DNA-generated taxa lists⁵⁴; we identified all possible pair-wise trophic interactions using existing structured trophic knowledge, information mined from the literature and a rule-based approach (Fig. 1). We first determined the complete taxonomic tree for the OTUs using DNA reference databases and the Name Parser API of the Global Biodiversity Information Facility (GBIF)⁵⁵. We next created a database of pair-wise biotic interactions (predation, parasitism, symbiosis) for our complete taxonomic tree by compiling information from available databases, including NemaGuild, FunGuild²⁸, BETSI (http://betsi.cesab.org/), GloBI²⁷ and systematic literature searches. In a limited number of cases where data was lacking, we generalized the information on interactions within families. For example, in the absence of available data on species-level interactions, members of a fungivorous genus of nematodes were assumed to feed on all fungal taxa in the network. The database was trimmed based on species co-occurrence and habitat data; pairwise 'resource-consumer' interactions were only included for organisms known to co-occur in a soil layer (soil surface, 0–10 cm, > 10 cm). The database of plant and animal taxa was also trimmed based on geographical range, targeting taxa known to occur in temperate Europe (taxa of microorganisms were considered to have a broad geographical distribution). Local food webs for each plot were extracted from the metaweb based on the list of taxa present in the four sets of DNA sequences for each plot.

Given that many taxa share the same sets of resources and consumers, partitioning food webs into trophic groups can be a useful tool to better identify network structure and function⁵⁶,⁵⁷. Aggregating taxa into trophic groups allows the simplification of complex food webs whilst preserving the information content of the initial network structure. With this approach, information on intra-specific plasticity is lost, but we assume that inter-specific plasticity is greater than intraspecific plasticity and drives network structure responses at the trophic group level. We assigned taxa to trophic groups using an automatic inference technique and machine-learning approach to avoid limitations of subjective, user-defined groups⁵⁸; stochastic block models (SBM)⁵⁹ were fitted to the heuristic metaweb using the "blockmodels" package in R (Supplementary Fig. S3). This approach generates trophic groups based on a common resource/food source or predator irrespective of taxonomic affiliation (for example, nematodes and protists which feed on the same bacterial groups will be clustered together). The optimal number of groups was determined using the Integrated Complete-data Likelihood criterion (ICL)⁶⁰.
<table><caption>**Table 2.** Network metrics applied at experimental fields to assess soil biodiversity, trophic complexity and energy pathways of the soil food web<sup>63,64</sup>. The interaction network is described by nodes (N), which are either OTUs or trophic groups (G), and edges (L, links from a resource node to a consumer node). Edges are described by their type (predator, parasitic, symbiotic), and weighted by the probability of occurrence of the resource and consumer pair.</caption><thead><tr><td>Metric</td><td>Scope</td><td>Definition</td><td>Biological relevance</td></tr></thead><tbody><tr><td>Node richness</td><td>Community</td><td>Number of nodes S</td><td>Food web size. Indicator of potential relationships</td></tr><tr><td>Trophic group richness</td><td>Community</td><td>Number of distinct trophic groups present in a community|G|</td><td>Functional diversity</td></tr><tr><td>Trophic group entropy</td><td>Community</td><td>Distribution of read counts among groups in an observed community using Hill numbers framework</td><td>Indicator of equitability</td></tr><tr><td>Average degree</td><td>Node</td><td>Number of interactions in which the OTU or group is involved</td><td>Centrality index: indicator of complexity</td></tr><tr><td>Mean trophic level</td><td>Node</td><td>Ordinal classification based on the relative position in the network, defined as 1 + the weighted average of the trophic levels of its resources/hosts</td><td>Centrality index: Indicator of the variety of trophic levels/ resources exploited as a food source</td></tr><tr><td>Omnivory level</td><td>Node</td><td>Variety in the trophic levels of a consumer's food</td><td>Centrality index: Indicator of the degree of specialization</td></tr><tr><td>Link density</td><td>Network</td><td>Average number of edges per node L/S</td><td>Connectivity index: indicator of quantity of relationships</td></tr><tr><td>Average path length</td><td>Network</td><td>Mean distance (d) between any pair of nodes in the network Σ<sub>ij</sub> d<sub>ij</sub>/S<sup>2</sup></td><td>Connectivity index: indicator of trophic redundancy and food web resilience</td></tr><tr><td>Maximum path length</td><td>Network</td><td>Longest path between any pair of nodes in the network max<sub>i≠j</sub> (d<sub>ij</sub> ∈ 1..S) ⌊d<sub>ij</sub>⌋</td><td>Connectivity index: indicator of network length</td></tr><tr><td>% parasitic links</td><td>Energy pathways</td><td>Proportion of total links that are parasitic</td><td>Diversity of interaction types</td></tr><tr><td>Bacteria:Fungi path ratio</td><td>Energy pathways</td><td>Ratio of decomposer-bacteria chains to saprotroph-fungi total pathway lengths</td><td>Indicator of bacterial dominance and N cycling efficiency</td></tr><tr><td>Detritivore:Herbivore path ratio</td><td>Energy pathways</td><td>Ratio between detritivore and root herbivory total pathway lengths</td><td>Indicator of importance of the detrital food path (dead vs. living organisms as a resource)</td></tr><tr><td>Node turnover (hill number: q = 0, 1,2)</td><td>Network stability</td><td>Turnover in the composition of the network (OTU or group)</td><td>β-diversity</td></tr></tbody></table>

In addition to partitioning taxa into groups, the stochastic block model 'learns' the probability of interactions between trophic groups that generate the trophic group metaweb. Heuristic food webs based on DNA metabarcoding have previously been validated as a powerful tool in freshwater systems<sup>61</sup>.

**Statistical and network analyses.** Standard metrics in connection with α-diversity and network robustness were used to describe the topology and properties of the 64 food webs inferred in the 16 experimental fields (Table 2, see Supplementary Fig. S4 for correlations between indices); in particular we focused on trophic network size, diversity of taxa or trophic groups (α-diversity), connectance, centrality and the nature of energy flow pathways (an indicator of network functioning) using the R package *Netindices*. In addition, we assessed multi-trophic β-diversity using the Simpson dissimilarity index and calculations of pairwise network dissimilarity between the 16 experimental fields. We computed the turnover in the composition (nodes) of both the trophic group networks and the taxa-based food webs using Hill numbers<sup>62</sup>) and the R package *econetwork*. We then applied generalized linear mixed effect models (GLMM) to test whether individually, these characteristics differed along our experimental land-use intensity gradient. Land-use intensity score and site identifier were used as fixed and random effects respectively:

$$m_i \sim \text{LUI}_i + 1|\text{site}_i$$

To visualize the overall effects of land-use intensity on β-diversity, we computed the first two components of the principal coordinate analysis (pcoa function from package *ape*) on the matrix of plot-by-plot network dissimilarities (Hill number = 0).

## Data availability

The complete list of taxa and associated trophic groups is provided as a data file. The datasets generated during and/or analysed during the current study are available on request. All soil samples were collected with permission from ANAEE-F ACBB site managers, and in accordance with ANAEE-F ACBB guidelines. No plant parts were used in the present study.

Received: 16 March 2021; Accepted: 20 August 2021
Published online: 14 September 2021

## References

1. Brussaard, L., Ruiter, P. C. & Brown, G. G. Soil biodiversity for agricultural sustainability. *Agric. Ecosyst. Environ.* **121**, 233–244 (2007).

2. Nielsen, U. N., Wall, D. H. & Six, J. Soil biodiversity and the environment. *Annu. Rev. Environ. Resour.* **40**, 63–90 (2015).
3. El Mujtar, V., Muñoz, N., McCormick, B. P., Pulleman, M. & Tittonell, P. Role and management of soil biodiversity for food security and nutrition; where do we stand?. *Glob. Food Secur.* **20**, 132–144 (2019).

4. Bardgett, R. D. & Putten, W. H. Belowground biodiversity and ecosystem functioning. *Nature* **515**, 505–511 (2014).

5. Kardol, P. & De Long, J. R. How anthropogenic shifts in plant community composition alter soil food webs. *F1000Res* **7**, 4 (2018).

6. Smith, P. et al. Global change pressures on soils from land-use and management. *Glob. Change Biol.* **22**, 1008–1028 (2016).

7. Geisen, S. et al. A methodological framework to embrace soil biodiversity. *Soil Biol. Biochem.* **136**, 107536 (2019).

8. Creamer, R. E. et al. Ecological network analysis reveals the inter-connection between soil biodiversity and ecosystem function as affected by land use across Europe. *Appl. Soil. Ecol.* **97**, 112–124 (2016).

9. Tsiafouli, M. A. et al. Intensive agriculture reduces soil biodiversity across Europe. *Glob. Change Biol.* **21**, 973–985 (2015).

10. de Vries, F. T. et al. Soil food web properties explain ecosystem services across European land use systems. *PNAS* **110**, 14296–14301 (2013).

11. Ponge, J. F. et al. Collembolan communities as bioindicators of land-use intensification. *Soil Biol. Biochem.* **35**, 813–826 (2003).

12. Postma-Blaauw, M. B., de Goede, R. G. M., Bloem, J., Faber, J. H. & Brussaard, L. Soil biota community structure and abundance under agricultural intensification and extensification. *Ecology* **91**, 460–473 (2010).

13. Decaëns, T. & Jiménez, J. J. Earthworm communities under an agricultural intensification gradient in Colombia. *Plant Soil* **240**, 133–143 (2002).

14. Dequiedt, S. et al. Biogeographical patterns of soil molecular microbial biomass as influenced by soil characteristics and management. *Globa. Ecol. Biogeogr.* **20**, 641–652 (2011).

15. Thomson, B. C. et al. Soil conditions and land-use intensification effects on soil microbial communities across a range of European field sites. *Soil Biol. Biochem.* **88**, 403–413 (2015).

16. de Graaff, M. A., Hornslein, N., Throop, H., Kardol, P. & van Diepen, L. T. A. Effects of agricultural intensification on soil biodiversity and implications for ecosystem functioning: A meta-analysis. *Adv. Agron.* **155**, 1–44 (2019).

17. Karimi, B. et al. Biogeography of soil bacterial networks along a gradient of cropping intensity. *Sci. Rep.* **9**, 3812 (2019).

18. Wardle, D. A., Nicholson, K. S., Bonner, K. I. & Yeates, G. W. Effects of agricultural intensification on soil-associated arthropod population dynamics, community structure, diversity and temporal variability over a seven-year period. *Soil Biol. Biochem.* **31**, 1691–1706 (1999).

19. Gossner, M. M. et al. Land-use intensification causes multitrophic homogenization of grassland communities. *Nature* **540**, 266–269 (2016).

20. Valiente-Banuet, A. et al. Beyond species loss: The extinction of ecological interactions in a changing world. *Funct. Ecol.* **29**, 299–307 (2015).

21. Freilich, M. A., Wieters, E., Broitman, B. R., Marquet, P. A. & Navarrete, S. A. Species co-occurrence networks: Can they reveal trophic and non-trophic interactions in ecological communities?. *Ecology* **99**, 690–699 (2018).

22. Gray, C. et al. FORUM: Ecological networks: The missing links in biomonitoring science. *J. Appl. Ecol.* **51**, 1444–1449 (2014).

23. Evans, D., Kitson, J., Lunt, D., Straw, N. & Pocock, M. Merging DNA metabarcoding and ecological network analysis to understand and build resilient terrestrial ecosystems. *Funct. Ecol.* **30**, 1904–1916 (2016).

24. Ruppert, K. M., Kline, R. J. & Rahman, M. S. Past, present, and future perspectives of environmental DNA (eDNA) metabarcoding: A systematic review in methods, monitoring, and applications of global eDNA. *Glob. Ecol. Conserv.* **17**, 00547 (2019).

25. Taberlet, P., Coissac, E., Hajibabaei, M. & Rieseberg, L. H. Environmental DNA. *Mol. Ecol.* **21**, 1789–1793 (2012).

26. Vacher, C. et al. Chapter one—Learning ecological networks from next-generation sequencing data. *Adv. Ecol. Res.* **54**, 1–39 (2016).

27. Poelen, J. H., Simons, J. D. & Mungall, C. J. Global biotic interactions: An open infrastructure to share and analyze species-interaction datasets. *Ecol. Inform.* **24**, 148–159 (2014).

28. Nguyen, N. H. et al. FUNGuild: An open annotation tool for parsing fungal community datasets by ecological guild. *Fungal Ecol.* **20**, 241–248 (2016).

29. Dopheide, A. et al. Rarity is a more reliable indicator of land-use impacts on soil invertebrate communities than other diversity metrics. *Elife* **9**, e52787 (2020).

30. García-Callejas, D., Molowny-Horas, R. & Araújo, M. B. Multiple interactions networks: Towards more realistic descriptions of the web of life. *Oikos* **127**, 5–22 (2018).

31. Delmas, E. et al. Analysing ecological networks of species interactions. *Biol. Rev.* **94**, 16–36 (2019).

32. Morrison, B. M. L., Brosi, B. J. & Dirzo, R. Agricultural intensification drives changes in hybrid network robustness by modifying network structure. *Ecol. Lett.* **23**, 359–369 (2020).

33. Banerjee, S. et al. Agricultural intensification reduces microbial network complexity and the abundance of keystone taxa in roots. *ISME J.* **13**, 1722–1736 (2019).

34. Thakur, M. P. & Wright, A. J. Environmental filtering, niche construction, and trait variability: The missing discussion. *Trends Ecol. Evol.* **32**, 884–886 (2017).

35. Xue, L. et al. Long term effects of management practice intensification on soil microbial community structure and co-occurrence network in a non-timber plantation. *For. Ecol. Manag.* **459**, 117805 (2020).

36. Felipe-Lucia, M. R. et al. Land-use intensity alters networks between biodiversity, ecosystem functions, and services. *PNAS* https://doi.org/10.1073/pnas.2016210117 (2020).

37. Heemsbergen, D. A. et al. Biodiversity effects on soil processes explained by interspecific functional dissimilarity. *Science* **306**, 1019–1020 (2004).

38. Erdozain, M. et al. Metabarcoding of storage ethanol vs. conventional morphometric identification in relation to the use of stream macroinvertebrates as ecological indicators in forest management. *Ecol. Indic.* **101**, 173–184 (2019).

39. Moore, J. C., McCann, K., Setälä, H. & De Ruiter, P. C. Top-down is bottom-up: Does predation in the rhizosphere regulate above-ground dynamics?. *Ecology* **84**, 846–857 (2003).

40. Wollrab, S., Diehl, S. & De Roos, A. M. Simple rules describe bottom-up and top-down control in food webs with alternative energy pathways. *Ecol. Lett.* **15**, 935–946 (2012).

41. de Vries, F. T. & Wallenstein, M. D. Below-ground connections underlying above-ground food production: A framework for optimising ecological connections in the rhizosphere. *J. Ecol.* **105**, 913–920 (2017).

42. de Vries, F. T. & Caruso, T. Eating from the same plate? Revisiting the role of labile carbon inputs in the soil food web. *Soil Biol. Biochem.* **102**, 4–9 (2016).

43. Bramon Mora, B., Gravel, D., Gilarranz, L. J., Poisot, T. & Stouffer, D. B. Identifying a common backbone of interactions underlying food webs from different ecosystems. *Nat. Commun.* **9**, 2603 (2018).

44. Socolar, J. B., Gilroy, J. J., Kunin, W. E. & Edwards, D. P. How should beta-diversity inform biodiversity conservation?. *Trends Ecol. Evol.* **31**, 67–80 (2016).

45. Malik, A. A. et al. Land-use driven change in soil pH affects microbial carbon cycling processes. *Nat. Commun.* **9**, 3591 (2018).

46. Jia, Y. & Whalen, J. K. Functional redundancy and phylogenetic niche conservatism in the soil microbial community. *Pedosphere* **30**, 18–24 (2020).

47. Bush, A. et al. DNA metabarcoding reveals metacommunity dynamics in a threatened boreal wetland wilderness. *Proc. Natl. Acad. Sci.* **117**, 8539–8545 (2020).

48. Blüthgen, N. et al. A quantitative index of land-use intensity in grasslands: Integrating mowing, grazing and fertilization. *Basic Appl. Ecol.* **13**, 207–220 (2012).
49. Ruiz-Martinez, I., Marraccini, E., Debolini, M. & Bonari, E. Indicators of agricultural intensity and intensification: A review of the literature. *Ital. J. Agron.* **10**, 74–84 (2015).

50. Taberlet, P., Bonin, A., Zinger, L. & Coissac, E. *Environmental DNA: For Biodiversity Research and Monitoring* (OUP Oxford, Oxford, 2018).

51. Taberlet, P. et al. Soil sampling and isolation of extracellular DNA from large amount of starting material suitable for metabarcoding studies. *Mol. Ecol.* **21**, 1816–1820 (2012).

52. Boyer, F. et al. obitools: A unix-inspired software package for DNA metabarcoding. *Mol. Ecol. Resour.* **16**, 176–182 (2016).

53. Zinger, L. et al. Extracellular DNA extraction is a fast, cheap and reliable alternative for multi-taxa surveys based on soil DNA. *Soil Biol. Biochem.* **96**, 16–19 (2016).

54. Compson, Z. G. et al. Chapter two—Linking DNA metabarcoding and text mining to create network-based biomonitoring tools: A case study on boreal wetland macroinvertebrate communities. *Adv. Ecol. Res.* **59**, 33–74 (2018).

55. G.B.I.F. GBIF backbone taxonomy. (2017).

56. Allesina, S. & Pascual, M. Food web models: A plea for groups. *Ecol. Lett.* **12**, 652–662 (2009).

57. Gauzens, B., Thébault, E., Lacroix, G. & Legendre, S. Trophic groups and modules: Two levels of group detection in food webs. *J. R. Soc. Interface* **12**, 20141176 (2015).

58. Makiola, A. et al. Key questions for next-generation biomonitoring. *Front. Environ. Sci.* **7**, 197 (2020).

59. Nowicki, K. & Snijders, T. A. B. Estimation and prediction for stochastic block structures. *J. Am. Stat. Assoc.* **96**, 1077–1087 (2001).

60. Biernacki, C., Celeux, G. & Govaert, G. Assessing a mixture model for clustering with the integrated completed likelihood. *IEEE Trans. Pattern Anal. Mach. Intell.* **22**, 719–725 (2000).

61. Compson, Z. G. et al. Network-based biomonitoring: Exploring freshwater food webs with stable isotope analysis and DNA metabarcoding. *Front. Ecol. Evol.* **7**, 395 (2019).

62. Ohlmann, M. et al. Diversity indices for ecological networks: A unifying framework using Hill numbers. *Ecol. Lett.* **22**, 737–747 (2019).

63. Gauzens, B., Legendre, S., Lazzaro, X. & Lacroix, G. Intermediate predation pressure leads to maximal complexity in food webs. *Oikos* **125**, 595–603 (2016).

64. Lau, M. K., Borrett, S. R., Baiser, B., Gotelli, N. J. & Ellison, A. M. Ecological network metrics: Opportunities for synthesis. *Ecosphere* **8**, 01900 (2017).

## Acknowledgements

This research was financed by the INRA ECOSERV initiative and has benefited from the ANAEE-F ACBB infrastructures of Mons, Lusignan and Theix-Laqueuille (supported by ANR-11-INBS-0001) as well as CAP 20-25 support at Theix-Laqueuille (16-IDEX-0001). We thank K. Klumpp, F. Louault, G. Vitte and F. Gastal, principal investigators of these sites for their assistance in this work. Thanks also to A. Salcedo, J.-N. Gaillot, M. Mattei and A. Pinot for help with fieldwork, and to J. Amossé (CEREMA, France) and E. Lara (Real Jardín Botánico, CSIC, Spain) for taxonomic advice.

## Author contributions

J.M.G.B., P.C., P.T. and M.H. designed the study. J.M.G.B., P.C. and P.T. contributed to the fieldwork and/or to soil analysis. P.T. produced the metabarcoding data. M.H. and S.S. built the trophic groups and led the statistical analysis. J.M.G.B. and M.H. wrote the first draft of the manuscript with help from S.S., and all authors reviewed the manuscript.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1038/s41598-021-97300-9.

**Correspondence** and requests for materials should be addressed to J.M.G.B.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_7.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2021