Woody encroachment in grassland elicits complex changes in the functional structure of above- and belowground biota

. Woody plant encroachment affects dry grasslands globally. To predict changes in biodiversity and ecosystem processes, it is important to understand how this process affects the functional composition of grassland organism groups. In this context, seminatural wooded meadows represent a form of experimental manipulation — where open grassland and woody patches co-occur in homogeneous environmental conditions due to human management decisions — which provides an opportunity to address the effect of woody plant encroachment on vegetation and soil biota. We used environmental DNA metabarcoding to address variation in plant, soil fungal, and soil animal communities in parallel. We also addressed functional groups of fungi — animal and plant pathogens, saprotrophs, decomposers, arbuscular mycorrhizal, ectomycorrhizal, endophytic, and other symbiotrophic fungi — and of soil animals — fungivores, bacteri-vores, litter feeders, root feeders, macro plant feeders, algal/lichen feeders, predators, and parasites. Co-variation between communities was detected from aboveground vegetation plots and metabarcoding of soil DNA, in terms of estimated richness and compositional patterns. Differences between open and wooded patches were most pronounced among plants and symbiotic fungi, whereas soil animals exhibited less marked differences. For most organisms, mean richness, as well as total richness per habitat type, was higher in open than wooded patches, but ectomycorrhizal fungi exhibited the opposite pattern. The functional structure of the soil biotic community, as characterized by the proportion of DNA sequences attributed to different functional groups, differed signi ﬁ cantly between open and wooded grassland patches, with symbiotic fungi (arbuscular mycorrhizal, ectomycorrhizal, and other symbiotrophic [mostly orchid mycorrhizal] fungi) contributing most to the difference. This study supports the notion that a soil DNA-based metabarcoding approach can provide insights into the diversity and composition of multiple taxonomic groups in natural ecosystems. It also provides a ﬁ rst demonstration of the complex changes to the functional structure of the belowground community that accompany woody plant encroachment in grasslands.

and ecosystem processes, it is important to understand how this process affects the functional composition of grassland organism groups. In this context, seminatural wooded meadows represent a form of experimental manipulation-where open grassland and woody patches co-occur in homogeneous environmental conditions due to human management decisions-which provides an opportunity to address the effect of woody plant encroachment on vegetation and soil biota. We used environmental DNA metabarcoding to address variation in plant, soil fungal, and soil animal communities in parallel. We also addressed functional groups of fungi-animal and plant pathogens, saprotrophs, decomposers, arbuscular mycorrhizal, ectomycorrhizal, endophytic, and other symbiotrophic fungi-and of soil animals-fungivores, bacterivores, litter feeders, root feeders, macro plant feeders, algal/lichen feeders, predators, and parasites. Covariation between communities was detected from aboveground vegetation plots and metabarcoding of soil DNA, in terms of estimated richness and compositional patterns. Differences between open and wooded patches were most pronounced among plants and symbiotic fungi, whereas soil animals exhibited less marked differences. For most organisms, mean richness, as well as total richness per habitat type, was higher in open than wooded patches, but ectomycorrhizal fungi exhibited the opposite pattern. The functional structure of the soil biotic community, as characterized by the proportion of DNA sequences attributed to different functional groups, differed significantly between open and wooded grassland patches, with symbiotic fungi (arbuscular mycorrhizal, ectomycorrhizal, and other symbiotrophic [mostly orchid mycorrhizal] fungi) contributing most to the difference. This study supports the notion that a soil DNAbased metabarcoding approach can provide insights into the diversity and composition of multiple taxonomic groups in natural ecosystems. It also provides a first demonstration of the complex changes to the functional structure of the belowground community that accompany woody plant encroachment in grasslands.

INTRODUCTION
Increasing woody plant encroachment is occurring globally as a result of changing climate and land use (Eldridge et al. 2011, Archer et al. 2017. In central and northern Europe, woody encroachment is common in seminatural grasslands (Eriksson et al. 2002). These grasslands represent ecosystems that have developed and persisted within a natural forest biome due to historical and continuing human activities, including animal husbandry, mowing, and collection of firewood (Poschlod and WallisDeVries 2002). Cessation of management results in woody plant encroachment, usually accompanied by a decrease in floristic diversity (Dengler et al. 2014).
Wooded meadows represent a specific type of seminatural grassland where groups of trees and shrubs have historically been retained, usually to provide a supply of firewood (Kukk and Kull 1997). Such meadows provide an anthropogenic experimental setting for studying the potential effects of woody plant encroachment on local biota, because the presence and location of woody plants have been determined by human management decisions in an otherwise largely uniform habitat.
Biodiversity studies typically consider the responses of single trophic groups to environmental change, neglecting the fact that the response of a given trophic group may depend on the abundance and diversity of other trophic groups (Soliveres et al. 2016). Although the diversities of different taxonomic and functional groups tend to be positively correlated (Wolters et al. 2006), some groups may still exhibit different patterns (Gossner et al. 2016, Banerjee et al. 2018, Delgado-Baquerizo et al. 2019. The use of environmental DNA (eDNA) metabarcoding allows us to study large fractions of the entire biota (Calder on-Sanou et al. 2020) and to address co-variation among different taxonomic and functional groups (Prober et al. 2015, Janssen et al. 2018, Neuenkamp et al. 2018, Wubs et al. 2019. Such a multitaxon approach makes it possible to identify groups that exhibit high turnover rates along the spatial and temporal gradients of interest, thus providing indirect information about changes to ecosystem function. For instance, Ohlmann et al. (2018) found that fungi make a major contribution to the structuring of beta diversity across trophic groups, while saprophytic fungi in particular have been found to be strongly related to the turnover of the provision of multiple ecosystem functions along an elevational gradient (Martinez-Almoyna et al. 2019). Analogous studies are needed in order to understand which organisms are potentially the largest drivers of community structure change along different ecological gradients.
We were interested in the effect of woody vegetation on co-variation between different taxonomic and functional groups of soil biota in a wooded meadow. We used eDNA metabarcoding to address several taxonomic and functional groups in parallel, considering eukaryotes (in particular soil micro-and mesofauna), fungi, and more specifically arbuscular mycorrhizal (AM) fungi (Kress et al. 2015, Creer et al. 2017). In addition, we described plant communities using conventional vegetation plots in parallel with soil metabarcoding of the plant community to assess the extent to which the aboveground vegetation can be described by soil metabarcoding. We hypothesized that woody plants act as ecosystem engineers (sensu Jones et al. 1994Jones et al. , 1996 and drive the diversity and composition of the remainder of the community, in a meadow where initial climatic, topographic, and soil conditions are largely uniform. We expected parallel variation in diversity among most major taxonomic and functional groups, with correlated variation in community composition. Fungi have a significant role in structuring the biotic community in general (Ohlmann et al. 2018); for instance, ectomycorrhizal fungi (Becklin et al. 2012) and plants (Tedersoo et al. 2020) specifically exhibit competitive superiority over AM fungi and plants. Therefore, we expected that fungal communities respond to the presence or absence of (ectomycorrhizal) woody vegetation, with AM fungi exhibiting the strongest response because of the higher abundance of ectomycorrhizal competitors in wooded habitat patches. On the other hand, we only expected small to moderate differences in soil faunal composition, because earthworm, enchytraeid, macroarthropod, and microarthropod studies in similar climatic conditions have indicated fairly high similarity in abundance, biomass, and taxonomic composition between grasslands and deciduous woodlands (Rosswall et al. 1977, Axelsson et al. 1984, Andren et al. 1990, Taylor et al. 2019).

Study site
Laelatu wooded meadow is located in Western Estonia (58.59°N, 23.72°E). Mean annual precipitation in the region is 628 mm, mean annual temperature is 6.6°C, and mean monthly temperatures range from 17.8°C in July to À3.9°C in February (latest published climatic time series 1981-2010; Estonian Environment Agency 2020). The meadow of 153 ha has been a hayfield since at least the beginning of the 18th century, but probably for centuries earlier (Kukk andKull 1997, Kukk 2004). The mowed area began to decrease in the 1940s, and currently, only about 15 ha is used as a hayfield and mown annually. According to information from local landowners, currently wooded patches developed around large tree individuals during the 1940s and further expanded during the 1960s and 1970s.
The physiography and soil conditions at the Laelatu meadow area were thoroughly described by Sepp and Rooma (1970), who provided a soil map. At sampling plots, the soil type was rendzic leptosol, with a humus layer of~20-25 cm depth and pH 6.5-7.0 covering the underlying limestone-rich parent material. The vegetation has been classified as a "Sesleria caerulea-Filipendula vulgaris" community (Krall and Pork 1970). Small-scale plant richness is very high in treeless patches (>40 vascular plant species/m 2 ; Wilson et al. 2012; but lower under the woody vegetation; Aavik et al. 2008). The standing crop is between 200 and 300 g/m 2 , and 98% of the plant species are perennials (Zobel and Liira 1997).

Sampling
Sampling was conducted in four plots of 30 9 30 m each in the central part of the meadow. Two plots were located in open parts of the meadow without any trees or shrubs (except single 1-yr-old saplings) within a radius of 30 m. The remaining two plots were located under small woody groves, mainly of Quercus robur and Corylus avellana, standing within a grassland matrix. It was not possible to establish more plots of this size because of the intermingled structure of wooded and open patches in the rest of the meadow, which resulted from a temporary cessation of mowing during the 1960s and 1970s (Aavik et al. 2008).
In both of the two plots in each of the habitat types (open and wooded), a 5-g soil sample was collected from nine points on a regularly spaced 30 9 30 m sampling grid in early July 2018, resulting in a total of 36 samples (2 habitat types 9 2 plots 9 9 samples). Up to 40% of the DNA extracted from soil samples has been noted to originate from extracellular material, that is, remnants of dead organisms (Carini et al. 2016), thus providing not only a strict temporal snapshot, but rather a more stable view of the structure of the belowground community. Each soil sample consisted of approximately 5 g of soil collected from a single spade-sized sample from the 2-10 cm topsoil. Soil samples were dried with silica gel and stored airtight at room temperature as in Saks et al. (2014) before further analysis. In addition, we collected topsoil (2-10 cm) samples for chemical analysis from each sampling point and pooled samples within each plot (altogether about 500 g per plot), which were subsequently air-dried at room temperature before analyses. For each sampling point, we also described a 1 9 1 m vegetation subplot where all vascular plant species in the ground-layer community were recorded and their abundance was estimated as percentage cover (Peet and Roberts 2013). Percentage cover estimates of shrubs and trees were based on the 10 9 10 m area around each subplot.

Soil chemical properties
We determined soil pH, organic N and C, and content of plant available P, K, Mg, and Ca in soils. All soil samples were air-dried after sampling at room temperature and sieved thereafter, using a sieve with 2-mm openings (Retch, Haan, Germany). Soil pH was measured in 1 mol/L KCl solution following ISO 10390:2005 using the Mettler Toledo pH meter SevenEasy with electrode Mettler Toledo InLab Expert Pro. The content of organic N in soil was determined by Kjeldahl method using the digestion block DK-20 and distillation unit UDK-126 produced by Velp Scientifica Srl (Usmate Velate, Italy). For the determination of organic carbon content in the soil, Tyurin's method was used, in which oxidation was provided by boiling the soil samples in v www.esajournals.org K 2 Cr 2 O 7 solution acidified with sulfuric acid (Vorobyova 1998). For determination of soil plant available P, K, Mg, and Ca, the Mehlich III extraction method was used (Soil and Plant Analysis Council 1992). The content of elements in the Mehlich III extract was determined with a microwave plasma atomic emission spectrometer MP-4200 (Agilent, Santa Clara, California, USA). Soil chemical analyses were performed at the Institute of Agricultural and Environmental Sciences, Estonian University of Life Sciences, Tartu.

Molecular methods
DNA was extracted from each individual 5-g sample of dried soil with the MoBio PowerMax Soil DNA Isolation Kit (MoBio Laboratories, Carlsbad, California, USA), with the following modifications described by Gazol et al. (2016): (1) Bead solution tubes were shaken at 60°C for 10 min at 100 rpm in the shaking incubator, and (2) the samples were dried for 10 min at room temperature under a fume hood before adding the final elution buffer.
AM fungal sequences were amplified from soil DNA extracts using the AM fungal-specific primers for the small subunit (SSU) ribosomal RNA gene V4 region: WANDA (Dumbrell et al. 2011) and AML2 (Lee et al. 2008). First, PCR was carried out with amplicon-specific primers linked to Illumina Nextera XT (Illumina, San Diego, California, USA) sequencing adapters (Illumina forward primer adaptor: 5 0 -TCGTCGGCAGCGTCA-GATGTGTATAAGAGACAG-3 0 ; Illumina reverse primer adaptor: 5 0 -GTCTCGTGGGCTCGGA-GATGTGTATAAGAGACAG-3 0 ). The reaction mix contained 12.5 µL 2x KAPA HiFi HotStart PCR mix (Kapa Biosystems, Wilmington, Massachusetts, USA); 1 µL of each 10 µmol/L primer; 6 µL of template DNA; and Milli-Q water to reach a total reaction volume of 25 µL. The PCR was performed under the following cycling conditions: 95°C for 3 min; 35 cycles of 95°C for 30 s, 55°C for 30 s, and 72°C for 30 s; followed by 72°C for 5 min. PCR products were purified using Agencourt AMPure XP beads (Beckman Coulter, Brea, California, USA). A second PCR was performed with Nextera XT index adapters for multiplexing (index PCR). The reaction mix contained 15 µL 2x KAPA HiFi HotStart PCR mix; 5 µL of Nextera XT index 1 Primer (N7xxx) and 5 µL of Nextera XT index 2 Primer (E5xxx); and 5 µL of purified PCR amplicons. PCR was performed under the following cycling conditions: 95°C for 3 min; 7 cycles of 95°C for 30 s, 55°C for 30 s, and 72°C for 30 s; followed by 72°C for 5 min. Indexed amplicons were purified with Agencourt AMPure XP beads. Indexed and purified amplicons were quantified with Qubit 2.0 Fluorometer (Invitrogen, Grand Island, New York, USA) and pooled in equal amounts (each 3 µL) to a library. Libraries were submitted to quality control in the 2200 TapeStation (Agilent) and subsequently diluted to 4 nmol/L and quantified with KAPA qPCR. Final libraries were denaturated with NaOH (0.1 N) and sequenced on the Illumina MiSeq platform with 2 9 300 bp paired-end sequencing chemistry.
The total fungal community was identified by sequencing the ITS2 region with degenerate primer pair fITS7 and fITS7o (forward) and ITS4 (reverse primer; Ihrmark et al. 2012, Kohout et al. 2014. Degenerate primers were used to allow for maximum detection of various fungal taxa. Plants were identified by sequencing the chloroplast trnL region with primers trnL(UAA)g and trnL(UAA) h (Taberlet et al. 2007). The total eukaryotic community was identified by sequencing the 18S SSU rRNA gene V4 region with primers F574 and R952 (Hadziavdic et al. 2014). PCR conditions for total fungal and total eukaryotic communities were the same as above. The PCR cycling conditions for DNA extracted with the plant-specific primer were as follows: an activation step of 10 min at 95°C; 35 cycles of 30 s at 95°C, 30 s at 55°C, and 30 s at 72°C; followed by final extension step of 72°C for 10 min. The PCR cycling conditions for total eukaryotes were as follows: 95°C for 15 min; followed by 35 cycles consisting of 95°C for 45 s, 60°C for 45 s, and 72°C for 1 min; and a final extension step of 72°C for 10 min. Library preparation and sequencing were the same as above. Sequencing was conducted on the Illumina MiSeq platform with 2x250 bp chemistry. All sequencing was conducted at Asper Biogene LLC (Tartu, Estonia).

Bioinformatic methods
Illumina paired-end raw reads were demultiplexed and cleaned using a series of bioinformatics steps (Vasar et al. 2017). In general, reads were demultiplexed into samples based on double barcode sequences, while allowing one v www.esajournals.org mismatch for both reads. Additionally, primer sequences were checked, allowing one mismatch. Barcode and primer sequences were removed, and only paired reads where both reads had an average quality score of >30 were retained. Paired-end reads were combined with FLASh (v1.2.10, Mago c and Salzberg 2011), using the default thresholds: overlap between 10 bp and 300 bp and overlap identity at least 75%, orphaned reads were removed. Putative chimeric sequences were identified and removed using vsearch in de novo mode, when using GenBank (Clark et al. 2016) for taxonomy assignment and reference database mode, when using other databases. Sequences were then submitted to a BLAST+ (v2.5.0, Camacho et al. 2009) search with the blastn algorithm against databases described below.
In the cleaning and quality filtering stage, additional trimming was performed for fungal ITS and chloroplast trnL reads: ITS reads were trimmed at 220 bases for forward and 180 bases for reverse reads, as suggested by Lekberg et al. (2018); trnL reads contained Illumina-specific overhang adapter sequences which occur when the insert size is shorter than the read length; the reads were therefore trimmed up to the adapter sequences. Before the chimera removal step, sequences from the fungal ITS, chloroplast trnL, and general eukaryotic amplicons were preclustered with 97% identity using vsearch (v2.14.1, Rognes et al. 2016) into operational taxonomic units (OTU) to reduce computational complexity and time. Cluster information was stored to allow mapping clusters back to individual reads.
AM fungal sequences (SSU) were matched against virtual taxa (VT; phylogenetically defined taxonomic units) in the MaarjAM database (status June 2019, € Opik et al. 2010), with 97% identity and 95% alignment thresholds. General fungal sequences (ITS2) were assigned to taxa in the UNITE database (status November 2018, Nilsson et al. 2018) with 97% identity and 95% alignment thresholds. To find additional hits for ITS2, sequences that did not get a sufficiently good match against UNITE were submitted to a GenBank search using a similar 97% identity and 95% alignment length threshold, and using the 500 best hits to determine a consensus taxonomy for each query (i.e., the lowest taxonomic rank common to >50% of the hits was accepted).
Matches from the two complementary searches were then concatenated to a single fungal ITS2 dataset. Plant trnL and general eukaryotic SSU sequences were identified using the following consensus taxonomy approach: (1) Sequences were submitted to a search against GenBank with 95% identity, 80% alignment, and 1e-10 BLAST score thresholds; (2) for each individual query sequence, the 500 best hits were retained; (3) among these hits, the lowest taxonomic rank common to >50% of the hits was recorded.
The approach of using a common taxonomic rank for each sequence cluster produced more reliable estimates of the taxonomic identity of the clusters, but has the caveat of producing a nested taxonomy that contains some very crude taxa and others with very fine resolution. The exact lists of plant and soil animal taxa should be interpreted in this light and taken only as a rough indication of which groups of taxa were present. From the perspective of community analyses, however, this approach does not pose a problem, because the sequence clusters (operational taxonomic units; OTUs) are generated in the same way as for the standard BLAST-based approach (see the approach used for AM fungi), with the difference only in how taxonomy is assigned to the OTUs.
We addressed three taxonomic groups: plants, soil micro-and mesofauna, and fungi. The two last groups were further divided based on their distinctive roles in the functioning of ecosystems. Animal taxa were grouped based on dietary traits according to expert opinion, namely fungivores, bacterivores, litter feeders, root feeders, macro plant feeders, algal/lichen feeders, predators, and parasitic animals. Fungal taxa were further classified into eight functional groups, based on the FUNGuild database (Nguyen et al. 2016), namely animal pathogens, plant/fungal pathogens (including parasites), saprotrophs (i.e., fungi whose main autecological niche is saprotrophy), fungal decomposers (i.e., fungi that fill the decomposition niche in a community-including saprotrophs, but also certain mycorrhizal fungi; Lindahl and Tunlid 2015), AM fungi, EcM fungi, fungal endophytes (but excluding fungi with recorded pathotrophic mode), and other symbiotrophs (mostly ecto-and orchid mycorrhizal and lichenized fungi). Only FUNGuild assignments with confidence levels of probable and highly probable were retained, whereas remaining fungal taxa were considered as undefined fungi. For functional group analysis, we aggregated the sequence data to the group level. In cases where only a phylogenetically high-level consensus classification was possible for soil animals, several dietary traits were attributed to a single OTU (based on the range of dietary traits which could possibly occur in the phylogenetic group as a whole). Even if an animal OTU could be identified at a high taxonomic resolution, these could frequently encompass several dietary trait groups. Those OTUs comprising multiple dietary groups were thus counted as members of each group. Sequences for fungal groups were extracted from the ITS2 dataset, with the exception of AM fungi, for which we used the SSU dataset; soil animal sequences were extracted from the general eukaryotic SSU dataset. For community analyses, we also created composite groups of all animals, that is, all OTUs from the general eukaryotic dataset that were identified to Metazoa; and all fungi (ITS), that is, all OTUs from the ITS2 dataset that were identified to the fungal kingdom. Overview tables of the OTUs from different major organism groups detected in the study, as well as the functional and dietary groupings of the fungal and animal OTUs, respectively, are available in Appendix S2.
Amplicon sequencing of 18 soil samples from wooded patches and 18 soil samples from open patches of the wooded meadow yielded, respectively, 1237 and 1345 fungal OTUs (ITS region), 73 and 80 AM fungal VT (AMF specific SSU region), 200 and 249 soil animal OTUs, and 38 and 53 plant OTUs. Vegetation plots yielded 33 and 77 plant species from wooded and open patches, respectively. More detailed information about cumulative OTU numbers and sequence counts can be found in Appendix S1: Table S3.
Across the entire study, in total, 94 plant species were detected in vegetation plots, while sequencing of plant DNA from soil yielded 59 roughly species-level OTUs that could be classified among representatives of land plants (Streptophyta, with the exclusion of bryophytes which were not assessed in the vegetation plots). Of the 59 OTUs from soil DNA, three could be identified to a species-level taxon that was also present aboveground and 26 could be identified to genera represented aboveground, while 25 plant families detected from soil were also present aboveground, representing 86% of sequences (see Appendix S1: Supporting Results S1 for detailed comparison of plants detected aboveground and from soil).
Representative sequences for each habitat type and taxon combination are accessible from the European Nucleotide Archive under Study ID PRJEB34711.

Statistical methods
All statistical analyses were performed in R 4.0.0 (R Core Team 2020). To compare richness of soil biota, plant community richness, and compositional differences between the two habitat types, we fitted linear mixed models (LMM) and used permutational multivariate ANOVAs (PER-MANOVA). LMMs (function lmer() from R package lme4; Bates et al. 2015) used taxon richness as the dependent variable, habitat as an independent variable, and plot as a random factor. P values were obtained with the package multcomp (Hothorn et al. 2008), and R package MuMIn (Barto n 2020) was used to calculate conditional pseudo-R 2 values for the mixed-effects models. We confirmed the results of richness models by running parallel analyses with Chao extrapolated Shannon diversity (Chao et al. 2014) as the dependent variable. Since the patterns remained the same, we present the results obtained using observed richness. To visualize compositional differences, we used NMDS ordination (function metaMDS() from R package vegan; Oksanen et al. 2019) of sample-wise Bray-Curtis distances derived from Hellinger-transformed (Legendre and Gallagher 2001) sequence counts or, in the case of plant survey data, Bray-Curtis distances of percentage cover estimates. When the communities of certain organism groups were not able to be ordinate sufficiently reliably (i.e., NMDS stress values > 0.2), we iteratively subset the data matrix by excluding an increasingly larger number of the least abundant taxa, until sufficient stress values were obtained. When stress values remained >0.2 after reaching a threshold of retaining at least 50% of the most abundant taxa, or when stress was below the default S min = 10 À4 threshold of the metaMDS() function (Appendix S1: Supporting Results S2), the dataset was excluded from ordination and subsequent co-variation analyses. We tested the effect of habitat on the composition of soil eukaryotic communities using PERMANOVA with 999 permutations (function adonis() from R package vegan), assuming a nested data structure of plots within habitat types. For data aggregated to the functional group level, we used two parallel approaches for testing compositional differences: (1) using the number of sequences belonging to each functional group in a sample and (2) using the number of taxa belonging to a functional group in a sample. For a coarse-scale overview of the biotic changes associated with woody plant encroachment, we performed a v 2 test (function chisq.test() in R) on the cumulative sequence count table of the groups of soil biota (fungal functional groups, animal dietary groups, and plants) in either habitat type. The standardized residuals (residual divided by its estimated standard error under H 0 of the v 2 test) of individual cells (abundances of a certain group in a particular habitat type) were then used to infer the relative effect of each group on the difference between the two habitat types.
Pairwise correlation (i.e., a scaled measure of covariance) among groups of soil organisms and the plant community was assessed using Procrustean randomization tests (Peres-Neto and Jackson 2001), using the functions procrustes() and protest() from the vegan package. To ensure equal dimensionality of community data, the matrices were ordinated to the first two axes with NMDS.

Above-and belowground plant community patterns
Vegetation plots indicated differences between wooded and open patches in mean plant species richness (open 28.7 AE 1 SE and wooded 10.2 AE 0.5 SE species; LMM P < 0.001) and in plant community composition (PERMANOVA R 2 = 0.5, P < 0.001). Using plant DNA sequences from soil, a similar pattern was detected, with wooded and open patches differing in mean richness (open 19.9 AE 1.4 SE and wooded 13.7 AE 0.8 SE; LMM P = 0.022) and community composition (PERMANOVA R 2 = 0.44, P = 0.002). Strong co-variation between plant communities detected by vegetation survey and from plant DNA in the soil was also revealed by Procrustes analysis (r = 0.879; P < 0.001; Fig. 1).

Habitat patterns among all studied groups
Thirteen out of eighteen studied organism groups exhibited different mean richness between the two studied habitat types (Table 1). In all cases, mean richness was higher in open patches, except for EcM fungi, which exhibited significantly higher richness in wooded patches.
In addition to the differences observed among plant communities (presented above), the two habitats hosted different communities of general soil fungi (ITS2; PERMANOVA R 2 = 0.2, P = 0.002) and soil animals (SSU V4; PERMA-NOVA R 2 = 0.13, P = 0.002), with the soil fungal communities of wooded habitat being distinct between the two sampling plots (Fig. 2). With the exceptions of litter and macro plant-feeding soil animals, all soil organism groups included differed highly significantly in their community compositions between the two habitat types v www.esajournals.org (Table 1), with habitat accounting for 9% (bacterivorous animals) to 44% (plants) of compositional variation.
Out of 171 tested pairwise compositional covariations between communities of the groups included in the study (Procrustes Protest analysis), 99 indicated significant correlation at the sample level among different organism groups (Appendix S1: Supporting Results S2). Notable significant results included correlations between community compositions of plants and all fungi (ITS2 general fungal amplicon; r = 0.9, P = 0.003; Appendix S1: Fig. S1), plants and all soil animals (r = 0.7, P = 0.003), plants and AM fungi (AM fungal-specific amplicon; r = 0.6, P = 0.003), and all fungi and all soil animals (r = 0.8, P = 0.003; Appendix S1: Fig. S1). All studied groups were included in the Procrustes analyses, with the exception of other symbiotrophic fungi, which occurred so sparsely that the community matrix could not be reliably ordinated.

Soil organism functional groups
In terms of habitat-wise differences at the organism group level, there was a marked effect of habitat type on the relative abundance (based on sequence count) of different soil organism groups (contingency table v 2 = 2.39 9 10 4 , df = 17, P < 0.001; Fig. 3), with AM, EcM, other symbiotrophic fungi, and fungal decomposers contributing most to the v 2 score. At group level, soil animal and fungal communities both differed in the incidence-based (number of taxa belonging to a group present in a sample) community composition (PERMANOVA R 2 = 0.42, P = 0.02 for animals; and R 2 = 0.53, P = 0.001 for fungi). Group abundance-based (number of sequences belonging to a group in a sample) community composition was different between habitats for fungi (PERMANOVA R 2 = 0.15, P = 0.02). Cumulative species richness was significantly different between the habitat types for plants (based on survey data), EcM fungi, and undefined fungi (i.e., fungi which could not be assigned to a functional group; Appendix S1: Table S3).

DISCUSSION
We showed that the presence of wooded patches in an otherwise homogeneous open grassland ecosystem impacts the diversity and composition of the wider eukaryotic community. There was a general tendency that mean richness was higher in open than wooded patches. EcM fungi exhibited the opposite pattern, indicating that some of the soil community might be linked to host presence. The functional structure of the belowground community, as characterized by the proportion of DNA sequences attributed to different functional groups, differed significantly between open and wooded grassland patches. To this difference, symbiotic fungi (AM, EcM, and other symbiotrophic fungi) contributed the most. Green arrows and triangles represent open habitats, and blue arrows and circles represent wooded habitats. Each point (triangle/circle) represents a single soil sample or the corresponding vegetation plot. Arrow length indicates the extent of mismatch between communities detected with the two approaches. Arrows originate from the community detected by sequencing and terminate at the communities detected from vegetation plots. The correlation-like statistic derived from the symmetric Procrustes sum of squares (r) and the significance of the Procrustes statistic (p) are reported.
v www.esajournals.org A majority of the soil eukaryotic groups exhibited significant co-variation in community structure, either with the plant community, or with other groups. Overall, our study represents the first demonstration of complex changes to the functional structure of soil biota accompanying woody plant encroachment in grasslands.
A decrease in plant diversity due to woody plant encroachment has been repeatedly reported from European seminatural grasslands WallisDeVries 2002, Dengler et al. 2014) including Laelatu wooded meadow (Aavik et al. 2008). Our results confirm these established patterns of plant diversity and demonstrate analogous patterns among other soil biota-higher diversity of fungi and soil animals in open than wooded grassland patches. Current evidence concerning correlated patterns of richness among different taxonomic groups is inconsistent (Wolters et al. 2006, Gossner et al. 2016, Banerjee et al. 2018, Delgado-Baquerizo et al. 2019, Noreika et al. 2019). In our system, the reverse diversity pattern-higher diversity in wooded than in open grassland patches-emerged only among EcM fungi, which reflects the presence of their host plant species (Schwob et al. 2017) in wooded patches.
Woody plant encroachment influences grassland plant community composition (Kahmen and Poschlod 2008, Eldridge et al. 2011, Ratajczak et al. 2012, Jakobsson et al. 2019. There is also evidence that woody plant encroachment causes changes in the composition of fungal (Yannarell et al. 2014, Schwob et al. 2017, Neuenkamp et al. 2018, Grau et al. 2019, nematode (Biederman and Boutton 2009), and collembolan (Kwok and Eldridge 2016) communities. Such widespread impacts were evidenced by simultaneous responses in multiple groups. Among others, we recorded significant pairwise correlations between community matrices of plants and fungi, and plants and soil animals. All fungal functional groups significantly co-varied with the plant community. Where most animal dietary groups likewise co-varied significantly with plant community, some of them-such as fungi and bacterivorous, parasitic, and litter-feeding animals-were less tightly interlinked with the vegetation (see Appendix S1: Supporting Results S2). Such dampening of cascading knock-on effects between trophic layers may perhaps be explained by patterns of generality or specific trophic associations. Considering, for example, fungal-feeding animals, we note that some or many may be more dependent on saprotrophic than on mycorrhizal fungi for food (Rem en et al. 2010)-with the former fungi showing less response to woody encroachment than the latter (Fig. 3). At the community level, the main difference in response to woody encroachment was still observed between kingdoms; where fungi exhibited clear responses, animals revealed less distinct patterns, with plants in between (Fig. 3). Inconsistent and/or nondirectional responses to successional change have been previously detected in chronosequences of soil fauna on abandoned arable land (e.g., Kardol et al. 2005). Again, such differential patterns among animals and fungi may be attributable to the trophic wiring of the system (see the previous paragraph), with plants and fungi being more intimately dependent on each other, but soil animals being characterized by omnivorous or generalistic feeding associations (see Digel et al. 2014).
As methodological caveats, we note that any approach to species detection is subject to limitations. For soil animals, eDNA methods can be impressively accurate in detecting a large part of species detected by conventional methods (e.g., Bienert et al. 2012) and thus provide a promising avenue for detecting successional impacts on soil fauna (Fernandes et al. 2018, Rota et al. 2020). Yet, the specific resolution achieved varies with a wealth of methodological choices (including storage and extraction methods and primer choice (Horton et al. 2017, Fernandes et al. 2018, Marquina et al. 2019 Fig. 3. The relative contribution of each studied functional group to the difference between wooded and open habitats. Contribution of each functional group and habitat combination to the overall v 2 score in a contingency table of cumulative organism group frequency (sequence counts of the organism groups among all samples from each habitat) is presented. Red colors indicate positive association with a habitat type; blue colors indicate negative association. Color intensity indicates the size of the standardized residual between expected and observed frequency (sequence counts). Plants, fungi, and soil micro-and mesofauna were addressed. reviews, see, e.g., Garlapati et al. 2019, Beng andCorlett 2020). Such methodological impacts preclude general statements about the overall fraction of species reliably detected by eDNA techniques, or about what constitutes a sufficient sample for detection of some specific fraction of species. Beyond technical constraints on sampling and laboratory protocols, species identification is also highly dependent on the reference data available for accurate species identification. In the current study, fungal sequences were identified against systematically curated databases, whereas animal and plant sequences were compared to the motley data available in GenBank. This very difference could perhaps partly explain the more distinct patterns detected among fungi than among animals or plants (see Fig. 3), as a lack of taxonomic resolution may act to smudge ecological patterns.
To gain a comprehensive insight into community differences, we addressed the relative abundances of fifteen functional groups, representing plants, seven fungal guilds and eight guilds of soil fauna. Wooded and open grassland patches varied in functional group structure, with mainly symbiotic fungi contributing most to the variation. Our results differ from those of Martinez-Almoyna et al. (2019) who addressed altitudinal changes in alpine grassland ecosystems and found that saprotrophic fungal communities exhibited the highest turnover. This discrepancy highlights the dependence of patterns of soil fungal composition on the particular ecological setting. It must be considered, however, that the results may be dependent on the decisions made in classifying fungi based on databases-for example, many symbiotrophic fungi in the FUN-Guild database are dual classified as saprotrophic-symbiotrophic, highlighting the importance of clear descriptions of the methods for interpreting database results. In the studied wooded meadow, the structure of symbiotic fungal communities exhibited the strongest contrast between wooded and open habitats. Indeed, mycorrhizal fungal communities differ between grasslands and forests (Tedersoo et al. 2014, Davison et al. 2015, P€ artel et al. 2017. Notably, the encroachment of woody plants into grassland will simultaneously change soil conditions (Montan e et al. 2010, Grau et al. 2019) and the availability of particular host plants (Schwob et al. 2017), which in turn drives the changes in fungal communities. We therefore suggest that the strongest contrasts of community composition among the multi-kingdom compositional shift in communities may be indicative of the most important ecological interactions or processes underlying community changes. In other words, by simultaneously resolving community structure across multiple, potentially interacting taxa at the same time, we are able to provide additional higher order information about the relative impact of environmental change on particular organism groups within the context of an overall shift.
We found co-variation between vegetation plots and metabarcoding of soil DNA in terms of estimated richness and compositional patterns among plant communities. This result corroborates earlier findings (Yoccoz et al. 2012, Edwards et al. 2018, Janssen et al. 2018, Tr€ ager et al. 2019, Calder on-Sanou et al. 2020. However, the molecular approach could not assign plants to species level, because the short trnL P6 loop marker region used in the analysis does not distinguish all species in certain groups including Asteraceae and Poaceae (Hiiesalu et al. 2014, Tr€ ager et al. 2019) and due to the incompleteness of reference databases. However, where the aim is to identify plant diversity dynamics rather than exact species lists, an approach based on plant eDNA sequencing appears as informative as traditional plant surveys.

CONCLUSIONS
Woody plant encroachment is a common phenomenon in many grasslands. Here, we find that woody plants have a major effect on multiple kingdoms of soil eukaryotic biota. In particular, we show largely parallel changes in the diversity of different taxonomic and functional groups, with the exception of EcM fungi. Moreover, we show that symbiotic fungal communities exhibit the strongest contrast between open grassland and woody patches, indicating that they, along with their host plants, are key elements of the biotic change in this system. Our results suggest that a multitrophic, cross-kingdom perspective adds much power to distinguishing the major community changes along ecological gradients of interest. Given the challenges of species detection and identification across diverse branches in v www.esajournals.org the tree of life, such a perspective can only be achieved though DNA-based methods.

ACKNOWLEDGMENTS
We are grateful to Maria Viketoft and Tryggve Persson for their expert classification of soil animals, and to the latter for added insightful comments on the faunal patterns detected. This research was supported by the University of Tartu (PLTOM20903) and the European Regional Development Fund (Centre of Excellence EcolChange).