Vegetation, topography, and soil depth drive microbial community structure in two Swedish grasslands

Abstract Soil microbial diversity and community composition are shaped by various factors linked to land management, topographic position, and vegetation. To study the effects of these drivers, we characterized fungal and bacterial communities from bulk soil at four soil depths ranging from the surface to below the rooting zone of two Swedish grasslands with differing land-use histories, each including both an upper and a lower catenary position. We hypothesized that differences in plant species richness and plant functional group composition between the four study sites would drive the variation in soil microbial community composition and correlate with microbial diversity, and that microbial biomass and diversity would decrease with soil depth following a decline in resource availability. While vegetation was identified as the main driver of microbial community composition, the explained variation was significantly higher for bacteria than for fungi, and the communities differed more between grasslands than between catenary positions. Microbial biomass derived from DNA abundance decreased with depth, but diversity remained relatively stable, indicating diverse microbial communities even below the rooting zone. Finally, plant-microbial diversity correlations were significant only for specific plant and fungal functional groups, emphasizing the importance of functional interactions over general species richness.


Introduction
Soil micr obial comm unities play an important r ole in gr assland ecosystem functioning, from nutrient cycling to carbon (C) sequestr ation (Pr ommer et al. 2020 ). By breaking down and conv erting or ganic matter into plant-av ailable nutrients and by establishing symbioses with plant roots, they can enhance soil fertility and contribute to ecosystem health (Morris et al. 2020 ). Ho w e v er, these pr ocesses can be affected by biodiversity loss and by shifts in microbial community composition driven by land use change (Sala et al. 2000, Cao et al. 2017. While human land use affects ca. 70% of Earth's surface and exerts pr essur e on global ecosystems, its effects on grasslands are particularly strong. Gr asslands cov er ca. 40% of the ice-fr ee terr estrial surface (Sala et al. 2017 ), and grazing and mowing are among the most spatiall y extensiv e land use pr actices (Erb et al. 2017 ). Conv ersions fr om gr assland to cr opland, fertilization or intensiv e gr azing can have long-lasting consequences for microbial communities and on the soil processes they mediate (Buckley and Schmidt 2001, Steenwerth et al. 2002, Habekost et al. 2008 ). Sweden has a long tradition of grassland management (Eriksson and Cousins 2014 ). Historicall y, gr asslands hav e been used to support liv estoc k for dairy and meat products, and the manure from livestock to fertilize crop fields. Ho w ever, over the past century, the abandonment of traditional grassland management practices, coupled with the conv ersion of v ast ar eas to cr oplands and for ests, has significantl y r educed the extent of pr esent-day gr asslands in Sweden to a fraction of their original size (Eriksson and Cousins 2014 ).
Additionall y, former cr op fields are sometimes r epur posed and mana ged as gr asslands, r esulting in the emer gence of a ne w type of gr assland. These cr opland-to-gr assland conv ersions could lead to incr eased or ganic matter inputs and plant-soil inter actions (Yang et al. 2021 ). While grassland plant communities have been extensiv el y studied, less is known about the status of microbial communities in Swedish grasslands on former arable fields and how they relate to existing plant communities.
Interactions between soil microorganisms and plant communities thr ough positiv e or negativ e feedbac k depend on plant species composition and diversity (Ettema and Wardle 2002, van der Heijden et al. 2008, Be v er et al. 2013. While gener al corr elations between abov egr ound and belowgr ound species div ersity may not be evident (Hooper et al. 2000 ), certain functional groups of fungi and bacteria, such as pathogens and root-associated organisms like arbuscular mycorrhizal fungi, have shown a positive association with grassland plant biomass and species richness (Wardle et al. 2004, Van Der Heijden et al. 2006, Be v er et al. 2013, Hiiesalu et al. 2014. This suggests that the nature and strength of the correlation between plant species richness and microbial diversity might differ between microbial functional groups . Furthermore , bacterial and fungal communities in soils can be expected to exhibit high variation across scales (Nunan et al. 2020 ), across both natural and anthropogenic (land use) gradients (Ettema 2002). For example, topogr a phicall y driv en v ariability in soil conditions explains differences in microbial community composition expressed in r elativ e abundance of specific fungal functional gr oups (Sc hlat-ter et al. 2018 ), and poorl y dr ained lo w er catenary positions can have higher microbial biomass compared to summits and slopes (Watts et al. 2010 ). Micr obial composition differ ences may then be explained by higher soil moisture and increased availability of C and nutrients due to both transport processes and local retention mechanisms at lo w er catenary positions. Micr obial comm unities also vary vertically in the soil profile. Microbial abundance is gener all y higher in the upper topsoil layers where root densities are highest (Berg et al. 1998 ) and in the rhizosphere (Lv et al. 2023 ) compared to bulk soil and decreases with soil depth as root biomass declines (Fierer et al. 2003 ). The r esulting v ertical str atification of micr obial comm unities is determined by differences in substrate availability, soil horizon properties, and soil a ggr egate sizes (Fox et al. 2018 ). This variability is lar gel y dependent on plant species identity and r oot distribution and activity, since plants provide soil microbial communities with substrates through aboveground and belowground litter input and root exudation (Wardle et al. 2004, Zhao et al. 2021. As different microbial functional groups vary in their response to substr ate av ailability, the decline in r oot biomass and in available C and nitrogen (N) with depth leads to shifts in microbial community composition, and may, for example, result in decreasing fungal div ersity (Sc hlatter et al. 2018 ). Ho w e v er, onl y a fe w studies have investigated how soil fungal and bacterial communities change with soil depth below the rooting zone in grasslands, although communities at different depths may play distinct roles for C and nutrient cycling.
Land management, plant community composition, topography, and soil depth may all be important in shaping microbial diversity and function. Her e, we c har acterized the micr obial comm unity composition of two grasslands with different land-use histories, whic h r eflected differ ences in the plant species ric hness, plant and functional group composition, and proximity to the forest border. By selecting two colocated grasslands with the same climate and similar in soil c har acteristics, suc h as pH, soil texture, and nutrient availability, w e w ere able to control for potentially confounding effects. In each grassland, the fungal and bacterial comm unities wer e obtained fr om bulk soil at an upper and a lo w er catenary position, i.e. sites with similar soil parent material but different slope positions. Soils were sampled at four depths ranging from 0 to 50 cm, with the deepest soil layer containing < 5% of the total root biomass. We investigated fungal and bacterial community composition in relation to edaphic properties and vegetation to identify potential drivers of fungal and bacterial diversity in grasslands and to test hypotheses that: (i) vegetation, catenary position, and soil depth affect the variation in soil fungal and bacterial community composition; (ii) with reference to soil depth in particular, microbial biomass abundance and diversity decrease with soil depth following a decline in resource availability; and (iii) micr obial comm unity div ersity corr elates with plant div ersity and plant functional group composition.

Site description and sampling
The sampling plots were established in 2019 in two grasslands in proximity to Tovetorp Research Station (south-central Sweden, Table S1). The two grasslands (Tov etor p and Ämtvik) are located < 4 km apart with similar soil parent material (post-glacial silty clay) and a humid continental climate with a mean annual precipitation of 619 mm and a mean annual temper atur e of 7.3 • C ( http://weather.zoologi.su.se ). Both gr asslands ar e former plo w ed arable fields, but with differences in time and management since cessation of a gricultur e . T he To v etor p location has been a conventionally plo w ed agricultural field until 1993, and later occasionally mo w ed and grazed b y fallo w deer. The Ämtvik location used to be a crop field belonging to a small farm, but was converted into a grassland and extensiv el y gr azed by cattle or mo w ed for at least 40 years prior to the study (information r etrie v ed fr om local farmers). While the abov egr ound biomass does not differ significantly between the two grasslands , To vetorp has a higher mean root biomass (Table S6). The vegetation in both grasslands is dominated by per ennial gr asses, but Ämtvik hosts consider abl y more forbs and Leguminosae (26% of the total aboveground plant biomass) than Tov etor p (9%), and has a higher mean plant species richness (Table S6). The sampling sites in Ämtvik are also closer to the edge of a mixed deciduous and coniferous forest (6-20 m, compared to 50-90 m in Tovetorp).
In eac h gr assland, 12 plots (2 × 2 m) were placed 0.5-2 m apart at a higher catenary position, and 12 similar plots were placed at a lo w er catenary position (hereafter called high/low ele v ation sites; see Fig. S1). While these paired high/low elevation sites were only ∼50 m apart and elevation differed only by ∼6 m, the lower elevation sites had on average ca. 10% higher volumetric soil moistur e measur ed at 60 cm depth. Vegetation inventories were carried out throughout the growing season of 2019 (April to October) to ca ptur e the maxim um cov er a ge of eac h species . T he presence of all plant species was recorded in each plot, and their cov er a ge was estimated to the nearest 1%. In addition to plant species richness and diversity calculated as Shannon-Wiener div ersity index ( H '), abov egr ound biomass (r eported as dry weight) was measured in mid-July 2019 by harvesting all the vegetation from one quarter (1 m 2 ) of e v ery plot by cutting at ground level, including moss and dead plant biomass. In half of all plots, the biomass from a 50 × 50 cm square was sorted into grasses , forbs , and legumes to obtain their r elativ e abundance in dry biomass (moss , woody plants , and dead plant biomass were measured, but excluded from the analyses). Distance from the closest trees (m) was measured to the center of each plot. Root biomass in the upper 30 cm was obtained from one soil core (8 cm diameter) sampled from each plot in September 2019, and the dry weight was recorded after rinsing the roots from soil by rinsing them with water on a 0.5 mm mesh sie v e, follo w ed b y drying at 60 • C for at least 48 h.
Soil moisture was measured every 3 weeks throughout the gr owing season fr om one access tube (1 m long) permanently installed in each plot, using a PR2 profile probe (Delta-T Devices Ltd, Cambridge , UK). T he values used in the analyses are growing season av er a ges of volumetric soil water content (%) in each plot, measur ed e v ery 10 cm in the first 50 cm. One soil sample for elemental and DNA-based analyses was collected in each plot in July and August 2019 with a Pürckhauer soil corer (2.5 cm diameter; Eijkelkamp, The Netherlands) from four depths (0-10 cm, 10-20 cm, 20-30 cm, and 40-50 cm) and immediatel y stor ed fr ozen. A portion of eac h sample was used for pH measurements (Mantech Automax 73, Guelph, ON, Canada), while the rest was fr eeze-dried, gr ound, and sie v ed. A subsample was sent to the Stable Isotope Facility at UC Davis (California) for estimation of total C and N contents ( μg/mg dry weight soil). Soil organic matter (SOM; %) was estimated from total organic matter content after loss on ignition (550 • C) from three additional soil samples collected from each site in September. The samples for the DN A analyses w ere obtained b y grouping the 12 plots of each grassland and elevation into three blocks and pooling the soil samples of the four adjacent plots of each block for each depth, thus resulting in a total of 48 pooled samples (three replicate samples per gr assland, ele v ation, and depth).

DN A extr action and quantifica tion of fungi and bacteria
DN A w as extr acted fr om a subsample of ∼100 mg of soil for soil depths 0-10 cm, 10-20 cm, and 20-30 cm and 400 mg for the deepest soil samples (40-50 cm) using the Mac her ey-Na gel Nu-cleoSpin Soil Kit (Dür en, German y), with the following specifications: For all samples from 0 to 30 cm depth, 700 μl lysis buffer SL1 and 150 μl enhancer SX were used, and for the deep soil samples, r ea gent quantities wer e adjusted for the lar ger soil weight to 1000 μl of lysis buffer SL1 and 214 μl of enhancer SX. First, the soil samples went through mechanical lysis with ceramic beads in a bead beater (Pr ecell ys, Bertin Tec hnologies SAS, FR), and the extracted DN A w as quantified using a NanoDr op spectr ophotometer (T hermo Scientific , Wilmington, DE, USA) and diluted to 1 ng/ul for fungal and bacterial PCRs.
The abundances of fungal and bacterial communities were estimated by quantifying the ITS2 region and the 16S rRNA gene, r espectiv el y. Quantitativ e PCR anal yses wer e run in duplicates in 20 μl reactions using 2 ng DNA template, IQ SYBR green supermix, and bovine serum albumin (BSA; 0.1%), as described in Castaño et al. ( 2022 ). Primers used for fungal qPCR w ere forw ar d primer gITS7 (GTGAR TCATCGAR TCTTTG; Ihrmark et al. 2012 ) and r e v erse primers ITS4 (75%; 5 -TCCTCCGCTT A TTGA T A TGC-3 , White et al. 1990 ) and ITS4a (25%; 5 -C AC ACGCTGTCCTCGCCTT A TTGA T A TGC-3 , Sterkenburg et al. 2015 ), and primers for bacterial qPCR were forward primer Pro341F (CCTACGGGNBGC ASC AG; Takahashi et al. 2014 ) and reverse primer 534R (A TT ACCGCGGCTGCTGG; Muyzer et al 1993 ). Before the analyses, PCR inhibition tests were run by adding a known quantity of standard plasmid (pGEM-T plasmid), which sho w ed no significant inhibition by sample extracts . T he reactions were run on a quantitative PCR cycler (BioRad iQ5, Life Technologies, Carlsbad, CA, USA), with cycling conditions of 5 min at 95 • C, 39 cycles of [15 s at 95 • C, 30 s at 56 • C, 40 s at 72 • C, and 5 s at 78 • C]. Standard curves for quantifications were obtained by serial dilutions of linearized plasmids containing either the ITS2 or 16S marker. The resulting copy numbers were corrected for the amount and the concentration of the DNA extracted from the original soil samples and calculated per g of dry soil. The fungal: bacterial ratio was calculated based on ITS and 16S rRNA gene copy numbers obtained through qPCR.

Sequencing of fungal and bacterial communities
Fungal and bacterial libraries were prepared by amplifying the ITS region and 16S rRNA gene, respectively, by using the same primers as for the qPCR for fungal ITS2 amplicons (see above), and forw ar d primer Pro341F (CCTACGGGNBGC ASC AG) and rev erse primer Pr o805R (GA CTA CNV GGGTATCTAATCC; Takahashi et al. 2014 ) for the bacterial V3-V4 region amplicons . T he ITS2 primers contained unique identification tags for each sample (Clemmensen et al. 2016 ), while the 16S rRNA-specific primers contained Nexter a ada ptor sequences (Illumina, San Diego, CA, USA). PCR for fungal DN A amplicons w as performed in duplicates of 50 μl reactions with 25 ng of DNA template, 0.2 mM of each nucleotide, 2.75 mM MgCl 2 , 0.5 μM of gITS7 primer, 0.3 μM of ITS4 primer , 0.1 μM of ITS4a primer , and 0.025 U μl −1 DreamTaq polymerase in 1 X DreamTaq buffer (Thermo Scientific, MA, USA). PCR conditions were 5 min at 95 • C, 25-35 cycles of [30 s at 95 • C, 30 s at 56 • C, and 30 s at 72 • C] and 7 min at 72 • C. The PCR products were purified using the AMPure kit (Beckman Coulter, CA, USA). The bacterial libr aries wer e pr epar ed using a two-step PCR protocol. The first PCR used a 15 μl reaction in duplicates with 4 ng of DNA template, BSA (20 μg μl −1 ), and 0.25 μM of each primer in 1 X Phusion High Fidelity PCR Master Mix (Thermo Scientific, MA, USA), follo w ed b y a second (indexing) PCR using 3 μl of purified product from the first PCR as template , BSA, polymerase , and 0.2 μM of each forw ar d and r e v erse m ultiplexing primers (containing Nextera barcoding regions for dual labeling). PCR conditions for the first r eaction wer e 3 min at 98 • C, 25 cycles of [15 s at 98 • C, 30 s at 55 • C, and 40 s at 72 • C] and 10 min at 72 • C, for the second (in which sample tags are added) 3 min at 98 • C, 8 cycles of [30 s at 98 • C, 30 s at 55 • C, and 45 s at 72 • C] and 5 min at 72 • C. All PCR runs were run on a SimpliAmp™ Thermal Cycler (Thermo Scientific, MA, USA). Bacterial PCR products were purified using Ser a-Ma g™ ma gnetic beads (GE Healthcar e, IL, USA).
The final PCR products for fungi and bacteria were quantified using a Qubit ® 2.0 fluor ometer (Invitr ogen, Carlsbad, CA, USA), pooled in equal DNA quantities, and the resulting pools cleaned using the E.Z.N.A. ® Cycle Pure Kit (Omega Bio-Tek, Norcross, GA, USA). For the fungal pool, adaptor ligation and sequencing were performed by NGI-Uppsala/SciLifeLab (National Genomics Infrastructure, Uppsala, Sweden) on a PacBio Sequel instrument (Pacific Biosciences, Menlo Park, CA, USA) using one Sequel SMRT cell (v3). The bacterial pool was sequenced on the MiSeq platform using the 2 × 250 bp paired-end chemistry (Illumina, San Diego, CA, USA).

Sequence processing and taxonomic classification
Fungal raw sequence counts were analyzed using the bioinformatics pipeline SCA T A ( https://scata.myk opat.slu.se ), wher e sequences were quality filtered, screened for primers and identification ta gs, and cluster ed accor ding to Ky asc henk o et al. ( 2017 ). Sequences with mean quality of < 20 or containing single bases with quality of three or less were removed, as well as all sequences containing mismatched identification tags. Sequences were clustered into study-le v el species hypotheses (SHs; Kõljalg et al. 2013 ) based on single-linkage clustering with a 1.5% threshold dissimilarity for sequences to enter an operational taxonomic unit (OTU). Fungal and nonfungal sequences were identified using BLASTn using the PlutoF platform (Abar enk ov et al. 2010 ) in the UNITE database (fungal database, v ersion 8.4), whic h bases identifications on SHs, and by constructing neighbor-joining trees in the MEGAN community edition (version 6.20.2019; Huson et al. 2011, Bálint et al. 2014 ). Re presentati ve sequences for the fungal OTUs were assigned to a ppr opriate taxonomic le v els depending on the similarity of the r epr esentativ e sequence ( > 97% match for species-level identification, 90% for genus, 85% for family, 80% for order, 75% for class, and 70% for division/phylum).
For fungi, sequencing of the dataset generated a total of 252 114 sequences, of which 163 904 passed quality control (QC) and were clustered into a total of 1619 O TUs . A total of fiv e hundr ed fiftytwo nonfungal OTUs (22 264 sequences, 24% of counts that passed QC) wer e r emov ed, as well as all OTUs in the global dataset with three counts or less, and 11 OTUs pr esent onl y in the negative controls with small count numbers . T he final fungal dataset consisted of 780 OTUs across all samples (69 292 sequence counts), for whic h r elativ e abundances wer e calculated per sample. To account for differences in fungal abundance and thus in the amount of DNA extracted from each soil sample in the analyses, and as a better r epr esentation of fungal biomass, the r elativ e species abundances were also multiplied by the total number of fungal copies per sample resulting from the qPCR (after eliminating the nonfungal proportions), as described in Parker et al. ( 2022 ). Identification of fungal functional groups was done using the FungalTraits database (Põlme et al. 2020 ) follo w ed b y careful manual curation, wher e unmatc hed sequences wer e classified into functional guilds based on taxonomic identity or similarity to database sequences recorded from well-defined substrates , if possible . About half of the OTUs (395 OTU) corresponding to almost 95% of the total sequences (65 769 sequences) were assigned to 12 functional groups based on primary lifestyle (animal parasites, mycoparasites , plant pathogens , endophytes , lic henized and lic hen par asites, wood sa pr otr ophs, litter sa pr otr ophs, dung sa pr otr ophs, soil sa pr otr ophs, undefined sa pr otr ophs, arbuscular mycorrhizae, and ectomycorrhizae), whereas the remaining 385 were labeled as "unkno wn" regar ding their ecology.
The 16S rRNA r aw r eads wer e first trimmed using the FASTX-toolkit ( http://hannonlab.cshl.edu/fastx _ toolkit ) and then merged using PEAR (Zhang et al. 2014 ). Removal of the c himer as was performed with UCHIME (Edgar et al. 2011 ), and dereplication and clustering of the sequences into OTU at 98% similarity threshold was performed with VSEARCH (Rognes et al. 2016 ). OTU clusters with < 2 reads across all samples were discarded. Re presentati ve sequences for each OTU were aligned and classified with the SILVA Incremental Aligner (SINA) using the SILVA 138 database as a r efer ence (Yilmaz et al. 2014 ). OTUs identified as mitochondria and chloroplasts as well as eukaryotes w ere discar ded resulting in a final dataset consisted of 6 975 146 sequence counts (excluding mitochondria and c hlor oplasts) gr ouped cluster ed into 7115 O TUs . Species accum ulation curv es were plotted for all samples (Fig. S2).
Fungal and bacterial sequences are available in the NCBI Sequence Read Arc hiv e ( www.ncbi.nlm.nih.gov/sr a ) under the accession number PRJNA938302.

Sta tistical anal yses
Differences in soil properties and vegetation and differences in microbial DN A cop y numbers obtained thr ough qPCR (her eafter called "abundance") between depths , grasslands , and elevations were tested with three-way ANOVAs (the microbial abundance after being r ank-tr ansformed). The dissimilarities between microbial fungal and bacterial community composition between grasslands, ele v ations, and soil depths were tested in R (version 3.3.3; R core Team 2017 ) by Analysis of Similarity (ANOSIM; R package: vegan, Oksanen et al. 2020 ). To further test which fungal divisions, bacterial phyla, and fungal functional groups differed in abundance between gr asslands, ele v ations, and soil depths, we used the mvabund function (R package: mvabund, Wang et al. 2021 ). This test is used to fit generalized linear models to each response variable (each species or group) and is suited for microbial species-abundance datasets.
The total fungal and bacterial community datasets were further analyzed using ordination methods in CANOCO 5.02 (Microcomputer Po w er, Ithaca, NY, USA). Graphical representation of the variation in community composition between grasslands, ele v ations, and soil depths was obtained with a principal correspondence analysis (PCA) summarizing the variation in species composition among samples. For the fungal communities, the PCA was based on the r elativ e abundance of 780 identified fungal O TUs , and data were not rarefied. For the bacterial communities, the PCA was based on the r elativ e abundances of 7115 identified O TUs , and data were rarefied by scaling with ranked subsampling according to Beule and Karlovsky (2020) to 59 701 sequences per sample. Correlation of fungal or bacterial communities with edaphic and plant variables was evaluated for significance using r edundancy anal ysis (RDA) with forw ar d selection and Monte Carlo permutations (without permutations within plots to account for dependency between soil depths from the same plot), in order to test potential drivers of the microbial communities. RDA anal yses testing onl y for effects of the categories grassland, ele v ation, and soil depth (sample categories), preceded testing the edaphic and plant variables (continuous variables). Vectors in the PCA species plots (unconstr ained anal ysis with supplementary variables function in CANOCO) indicate direction and degree of correlation between the PCA axes and the significant edaphic and plant v ariables, whic h wer e selected based on significance in the RDA analyses. One sample was lost during fungal sequence pr ocessing, r esulting in two replicates instead of three for the tr eatment Tov etor p low ele v ation (at 10-20 cm depth) and 47 fungal communities in total. The missing sample was replaced in the fungal RDA analysis by the mean number of counts and corr esponding mean v alues for plant and soil v ariables based on the two replicates within the same treatment, in order to ac hie v e a balanced split-plot design to account for soil depth dependency. The fungal and bacterial communities at the four soil depths were further tested separately with individual RDA analyses to e v aluate whether potential drivers of community composition depended on soil depth. For all analyses, sequencing output per sample was accounted for by including sequence counts as a co variate . Both fungal and bacterial communities were also analyzed using nonmetric multidimensional scaling (NMDS) summarizing the similarities in species composition between samples for comparison, since this method is commonly used for bacteria.
Fungal and bacterial OTU numbers (hereafter "species richness") and Shannon-Wiener diversity index, which reflects both richness and species relative abundances ( H ', hereafter "diversity"), were calculated for each sample for fungi (rarefied to 592 counts per sample) and bacteria (r ar efied to 59 701 counts per sample), and for the major fungal functional groups (saprotrophs, pathogens, and mycorrhizae) and the most abundant fungal and bacterial phyla (Ascomycota and Basidiomycota for fungi, Actinobacteriota, Pr oteobacteriota, Verrucomicr obiota, Acidobacteriota, and Chloroflexi for bacteria). To test for differences in species richness and diversity between grasslands , elevations , and soil depths, we used t -tests or ANOVAs followed by Tuk e y's HSD test. Corr elations between micr obial abundance and div ersity with v egetation and eda phic v ariables wer e tested with a corr elation matrix (Pearson correlation coefficient).
All pr edictiv e v ariables used in the statistical anal yses wer e selected after c hec king for m ulticollinearity with a corr elogr am and by calculating the variable inflation factor, which excluded the following variables from further analyses: total plant cover, ric hness of gr asses and forbs, pH, and soil C: N r atio. All analyses except for the ordination analyses (see above) were performed using R, and residuals were checked for normality.

Site characteristics
The grasslands at lo w er elevation had higher soil C ( P < .001) and root biomass ( P < .001) compared to the grasslands at higher ele v ation (Fig. 1 ). Differ ences between Ämtvik and Tov etor p wer e mostl y e vident in the plant community composition (more forbs and legumes, P = .001 and P = .02, r espectiv el y) and str onger influ- ence of nearby trees in Ämtvik. Soil C and root biomass decreased consistently with depth in the soil profile, whereas soil moisture and pH increased ( P < .001 for all four variables, Fig. 2 ).

Dissimilarity in microbial communities across sites, ele v a tion, and soil depths
The abundance of fungi and bacteria based on quantifications of marker genes decreased with soil depth (fungi: P < .001, bacteria: P < .001; Fig. 3 A, B), but fungal abundance did not differ between grasslands and elevations. Bacterial abundance was higher in the Tov etor p gr assland compar ed to Ämtvik ( P < .01) and in the lo w er ele v ations ( P = .04). Fungal and bacterial abundances wer e well correlated, so that the fungal: bacterial ratio did not differ between depths and ele v ations, e v en if it was ov er all higher in Ämtvik than in Tov etor p ( t = 2.3, P = .028) (Fig. 3 C).
The gener al tr end of decr easing micr obial abundance with depth but no significant difference in abundance between grasslands and ele v ations was true for most fungal functional groups and fungal and bacterial phyla, with the exception of mycorrhizal fungi (mvabund analyses, relative abundances presented in Fig. 4 ). Arbuscular mycorrhiza and ectomycorrhiza did not differ significantly in abundance through the depths consider ed, and onl y ectomycorrhiza differ ed significantl y between gr asslands, being pr edominantl y pr esent in Ämtvik (details in Table S4).
T he To v etor p and Ämtvik gr asslands and the two catenary positions were colonized by distinct soil fungal (Figs 5 A and S3a) and bacterial communities (Figs 5 B and S3b). For fungi, the variation across depths and replicates was somewhat higher in Ämtvik, and ele v ation a ppear ed to be a str onger driving force in Ämtvik compared to Tovetorp (Fig. 5 A). For bacteria, community variation was higher at Tov etor p (Fig. 5 B), and the high ele v ation emer ged instead as clearly distinct from the other three communities. Grass-land site explained 25.8% of the variation in fungal composition, while ele v ation explained 12.4% (RDA; categorical tr eatments, not shown). For bacterial communities, 43.6% and 20.4% of the variation were explained by grassland site and elevation, respectively (RDA; categorical treatments, not shown). Microbial communities also shifted with soil depth explaining an additional 3.9% of the variation for fungi (Fig. S3c) and an additional 10.7% of the variation for bacteria (Fig. S3d) (RDAs; gr assland, ele v ation, and depth, not shown). Significant differences in community composition among gr asslands, ele v ations, and depths wer e also identified in the NMDS (Fig. S4) and ANOSIM (Table S3) analyses, although the ANOSIM sho w ed higher correlation coefficient values for depth compared to grassland and elevation. RDA analyses of bacterial communities generally explained twice as much of the v ariation in comm unity composition than for fungi. When absolute comm unity v ariation (i.e. r elativ e abundances m ultiplied with total fungal copy numbers in samples derived by qPCR data, Figs S5 and S6) was tested in the same way, the resulting fungal compositional patterns stayed very similar, with the main difference being a clearer separation between Tovetorp high and low ele v ations. Variation in fungal and bacterial communities in the differ ent gr asslands and ele v ations was maintained at eac h individual soil depth (RDAs; Figs S7 and S8).

Dri v ers of microbial community composition
Corr elations between micr obial comm unity composition and edaphic and plant variables partly differed between fungi and bacteria. In the RDA analyses, the adjusted explained variation for plant and soil variables together was 26.5% for fungi and 49.4% for bacteria (Table S2). The drivers of fungal communities with the highest explanatory po w er in the RDA analyses with forw ar d selection w er e distance to the closest tr ees (ex-  plaining 12.5% of the variation) and number of plant species (8.0%) (Fig. 5 , Table S2). Fungal comm unity v ariation also corr elated significantly with plant root biomass, proportion of forb, and grass biomass , abo veground plant biomass , SOM, and soil moistur e (gr owing season av er a ge). pH was significantl y corr elated with SOM, and did not impr ov e the model when replacing SOM in the RDA. For bacterial communities, distance to closest trees (20.9%), number of plant species (15.0%), soil moisture (10.4%), and total soil N (8.5%) were the str ongest driv ers (Fig. 5 , Table S2), and plant root biomass, proportion of forbs, and grasses biomass and soil C were also significant. The other driv ers eac h explained 4.2% or less of the variation in microbial community compositions.
When each soil depth was tested independentl y, m ultiv ariate anal yses r e v ealed that distance to the closest tree and the number of plant species were significant drivers of both fungal and bacterial communities for all soil depths (Figs S7 and S8). For the bacterial communities, soil N was a significant driver  Table S4. for the two upper la yers , and soil C for the two deeper soil layers.

Dri v ers of microbial a bundance, ric hness, and di v ersity
Bacterial and fungal abundances were positively correlated with root biomass and SOM content (and both C and N contents), but negativ el y corr elated with soil moistur e (Fig. 6 ), r eflecting the v ertical tr ends noted in Fig. 3 . Bacterial communities had ov er all higher species richness than fungi (total mean sp. richness bacteria = 3145; total mean sp. richness fungi = 114) and higher diversity ( H ' = 3.7 in fungi, H ' = 6.3 in bacteria; P < .001, Fig. 7 A). Bacterial abundance and div ersity decr eased with depth ( P = .01 and P = .037, r espectiv el y; the differ ence was significant only between the 0-10 cm and 40-50 cm soil depths, data not shown) and did not differ between grasslands. Further, bacterial diversity was negativ el y corr elated with soil moistur e (Fig. 6 , Table S5a). Fungal diversity was similar across sites and soil depths (Fig. 7 A) and was not significantly correlated with any variable, although fungal species richness declined somewhat with depth ( P < .001, data not shown). Bacterial and fungal species richness, ho w ever, correlated to the same variables as bacterial and fungal abundances (positive correlations with SOM, total C, r oot biomass, and negativ e corr elation with soil moistur e, Table  S5).
Div ersity v alues differ ed also between the major fungal functional groups, in particular for mycorrhizal fungi compared to pathogens and sa pr otr ophs. Mycorrhizal fungal div ersity had the lo w est values and differed significantly between the two grasslands ( P < .001, mean H ' mycorrhizae in Ämtvik = 1.4 and Tovetorp = 0.3), while the diversity of pathogens and saprotrophs did not differ between the grasslands and ele v ations (Fig. 7 B). None of the fungal groups differed significantly in diversity between depths. Further, the diversity of mycorrhizal fungi correlated positiv el y with the proportion of forb and legume biomass, abov egr ound plant biomass, and distance from trees. In contrast, fungal sa pr otr ophs corr elated with forb biomass and were the onl y gr oup that sho w ed a positiv e corr elation with plant div ersity. P athogen fungal div ersity onl y corr elated with soil moistur e (Table S5b).
For fungi, diversity differed between the most abundant fungal divisions (mean H ' Ascomycota = 3.5, mean H ' Basidiomycota = 2.2; P < .001) and decreases with depth in Ascomycota ( P < .001). None of the groups, ho w ever, differed betw een grasslands, and only Basidiomycota had higher diversity in the high-elevation sites ( P < .001).

Gr assland site, ca tenary position, and soil depth dri v e microbial community structure
The fungal and bacterial communities were sampled from two gr asslands r elativ el y similar in climate and soil properties but that differed in land-use history, particularly with regard to time and management since cessation of intensive agriculture . T his legacy was expressed in the vegetation characteristics , i.e . root biomass , plant species richness, composition of plant functional groups, and proximity to the forest border. Differences were found in both fungal and bacterial community compositions across grasslands, ele v ations, and soil depths, with gr asslands-likel y deriv ed fr om land use-emerging as the strongest control in the RDA. This cor-  Table S5), with larger dots r epr esenting str onger significance. Fungal and bacterial div ersity indexes are calculated on the rarefied datasets. Correlations with the diversity of broad taxonomic groups and fungal functional groups are reported in the Supplements. r obor ates the first hypothesis that these three factors are important drivers of microbial community structure. Topography can affect micr obial comm unities b y influencing w ater or nutrient av ailability dir ectl y thr ough slope runoff or indir ectl y thr ough plant growth (Zhang et al. 2013 ). While soil moisture was negativ el y corr elated with micr obial abundance in our sites, catenary position had a significant effect on micr obial comm unity composition, but not abundance . T his suggests that certain microbial taxa might have been particularly sensitive even to slight variations in topogr a phy, and highlights the importance of integrating measur es of landsca pe heter ogeneity with land use for a comprehensive understanding of nutrient availability and limitations on soil microbial communities . Moreo ver, while microbial communities tend to be more homogeneous in deeper soil layers (Eilers et al. 2012 ), it is important to recognize that this pattern is influenced by both land use (van Leeuwen et al. 2017 ) and topogr a phy. Ther efore, it becomes essential to also consider soil depth as an additional dimension in the environmental gradient that contributes to the structure and function of microbial communities (Eilers et al. 2012 ). The integration of land use, topogr a phy, and soil depth will allow for a better understanding of the scales that determine microbial dynamics and spatial variation in soils. Ho w ever, these three factors explained twice the variation for bacteria than for fungi. This could partly be explained by the fact that fungal hyphae can form networks that cover extensive areas (Dahlberg andStenlid 1990 , Genney et al. 2006 ), grow across the depths and elev ations consider ed, and hav e longer life spans compared to bacteria. This means that despite the pooling of the soil core samples to attain r epr esentativ e samples at the plot scale, the lo w er explained variation of fungal community composition could be due to fungal communities being more spatially uniform than bacterial communities because of their specific biotic characteristics. In addition, the higher number of bacterial O TUs , resulting also from the different sequencing methods used, improves prediction potential for bacteria.

Partial evidence for plant-soil microbial comm unity feedbac k
In further support of our first hypothesis, the variation in micr obial comm unity composition was explained by plant variables (root biomass, species richness, proportion of grasses and forbs, and distance to the forest edge) for both fungi and bacteria, but also by some of the eda phic v ariables. Despite pr e vious studies (Dassen et al. 2017, Liu et al. 2020 ) suggesting a stronger link between vegetation and fungal community composition compared to bacteria, the total explained variation for fungi remained lo w er than for bacteria. The distance to the closest tree and the number of plant species were the drivers explaining most of the variation for both fungal and bacterial comm unities, and wer e also the only two drivers significant across all depths . T his suggests an important role of plant species identity and plant species richness in shaping both fungal and bacterial communities (Dassen et al. 2017 ) e v en in deeper soil la yers . In our study, plant species richness was determined by observations aboveground. Previous work has shown that plant identification through metabarcoding from root samples can reveal higher species richness compared to the values obtained by conventional plant inventories aboveground (Hiiesalu et al. 2014 ), which could improve the explanatory po w er of analyses of linkages between plant and micr obial comm unities. Further, soil pH is commonly a strong driver of microbial communities (Pellissier et al. 2014, Cao et al. 2017, Hao et al. 2021. Here, since pH varied with soil depth, but not considerably across grasslands or catenary positions, and was significantl y corr elated with SOM, it did not emerge as a significant predictor for bacteria and did not contribute to explaining the fungal community composition in this study. The higher abundance of ectomycorrhizal fungi in Ämtvik was lik ely dri ven by the proximity to the forest edge (6-20 m on aver a ge in Ämtvik, compar ed to 50-90 m in Tov etor p), as they ar e mostly associated to trees (Dickie and Reich 2005 ). While plantmycorrhizal associations may dir ectl y affect plant species distribution b y fav oring taxa that are connected to the hyphal network (Ettema and Wardle 2002 ), in this case, the proximity of trees was negatively correlated with grassland plant species richness, possibly due to shading and competition with tr ee r oots . T his may partly explain why arbuscular mycorrhizal fungi (typically associated with grassland species) were predominantly found in the more open locations, as was also observed by Moora et al. ( 2014 ). In addition, the covariation between arbuscular mycorrhizal fungi communities and plant communities decreases along a succession gradient from open grassland to forest (Neuenkamp et al. 2018 ). Ther efor e, it is not surprising that the distance from the forest edge emerges as an important driver of microbial community composition at our sites, suggesting that future studies on grassland biodiversity should integrate the patch size of grasslands and the proximity to forests in complex landscapes using spatial analyses also when observing soil microbial communities.

Depth effects on microbial abundance and di v ersity
Fungal and bacterial abundances decreased with depth, indicating a decline in microbial biomass and biological activity below the rooting zone, and thus confirming part of our second hypothesis. Suc h decr ease has also been sho wn b y other authors (Fierer et al. 2003, Hao et al. 2021, and has been attributed to a reduction of the availability of C sources with depth (Fierer et al. 2003 ). This interpretation is supported by 3) at four soil depths and for (B) major fungal functional groups (mean H ' mycorrhizae = 0.8; pathogens = 1.9; and sa pr otr ophs = 2.9) in two grasslands in south-eastern Sweden (A = Ämtvik and T = Tov etor p) and catenary positions (H = high and L = lo w). Bars sho w mean and interquartile range (IQR); whiskers extend to 1.5 x IQR; dots in the gr a ph ar e outliers. our r esults, whic h sho w ed a similar reduction in root biomass, soil C, and SOM with soil depth. Fungal abundance decreased with depth in all functional groups except for arbuscular mycorrhiza and ectomycorrhiza, which were unrelated to depth. This finding aligns with Schlatter et al. ( 2018 ) and suggests that the driv ers of comm unity composition of mycorrhizal fungi are distinct from those of other functional groups and might instead be more linked to specific plant host abundance and rooting depth.
Despite a significant decline in species ric hness, onl y bacterial Shannon diversity decreased slightly with depth, as in Eilers et al. ( 2012 ), thus partl y r ejecting our second hypothesis. It is possible that the r elativ el y high div ersity in deeper soil layers is a result of more scarce r esources pr omoting coexistence of microbial forag ing strateg ies, in contrast to richer surface conditions where fewer more competitive species can dominate.

Plant functional group composition is a better predictor of microbial di v ersity than plant di v ersity
Bacterial div ersity corr elated negativ el y onl y to soil moistur e, while fungal diversity sho w ed no significant correlations to edaphic or vegetation factors . T his partly rejects our third hypothesis of a general correlation between microbial diversity and plant diversity. It also indicates that the drivers of fungal and bacterial div ersity ar e separ ate fr om the driv ers of micr obial abundance. Pr e vious studies, ho w e v er, hav e shown that microbial diversity in grasslands may require larger-scale sampling to be accur atel y mapped (Pellissier et al. 2014 ), and that the interactions within microbial and between microbial and plant communities can be affected by the presence and ecosystem functions of soil macrofauna (Wardle et al. 2004 ), which was not included in this study. In addition, the low mycorrhizal abundance, particularly of arbuscular m ycorrhizae, was lik ely due to methodological limitations. First, the microbial DN A w as extr acted fr om bulk soil and not fr om r oots . Second, the primer gITS7 does not co ver all arbuscular mycorrhizal fungi (Cheeke et al. 2021 ), and more specific primers would be needed to ca ptur e that community.
Finally, our study highlights the significance of considering differ ent micr obial functional gr oups and their corr elations with envir onmental v ariables when examining ecosystem functioning and community response to environmental change, as previously suggested b y F rac et al. ( 2018 ). Consistent with earlier studies (Wardle et al. 2004, Hiiesalu et al. 2014, we found a strong association between m ycorrhizal di versity and vegetation properties, likely influenced by forest edge effects. Additionally, sapro-tr oph div ersity was corr elated with plant div ersity and v arious measures of plant community composition, as also observed by Str ec ker et al. ( 2015 ), indicating a potential impact of plant functional groups or species composition on litter quality (e.g. lignin and N contents). These findings suggest that correlations between vegetation and microbial diversity differ among microbial functional groups, and imply that plant functional groups might serve as better predictors of microbial diversity than plant diversity, as pr e viousl y demonstr ated for sa pr otr ophic fungi (Fr ancioli et al. 2021 ).

Conclusion
Fungal and bacterial community composition differed between gr asslands (potentiall y linked to differ ences in v egetation and land-use history), ele v ations, and soil depths . T he tested edaphic and plant variables generally explained twice as much of the comm unity v ariation for bacteria than for fungi. This differ ence could be due to the scale at which the drivers and the communities ar e deriv ed, or because other driv ers, suc h as biotic inter actions among fungi and between fungi and plant species, are relatively more important for fungi than for bacteria. By characterizing micr obial comm unities acr oss ele v ation positions under the same climate and soil, and along the soil depth profile, we identified r esource av ailability (r oot biomass and soil C and N contents), soil moisture, plant functional group composition, and species richness as the main drivers of microbial abundance and diversity. Microbial abundance and bacterial diversity decreased with depth, but fungal diversity remained stable. While there was no gener al corr elation between micr obial div ersity and plant div ersity, diversity within specific microbial groups (mycorrhizal and sa pr otr ophic fungi) was correlated with the occurrence of some plant functional groups, indicating interactions between ecosystem functions abov egr ound and belowgr ound. Finall y, integr ation of land use, topogr a phy, and soil depth will allow for a better understanding of the scales that determine microbial dynamics and spatial variation in grassland soils.

Funding
This r esearc h w as funded b y the Bolin Centre for Climate Resear ch at Stockholm University and by the Carl Mannerfelt Foundation. S.M. r eceiv ed funding fr om the Eur opean Researc h Council (ERC) under the European Union's Horizon 2020 r esearc h and innov ation pr ogr amme (gr ant a gr eement no 101001608). P .F . was supported by the Swedish Research Council FORMAS (2016-01107) and K.E.C. by FORMAS (2020-01110).