Site-Dependent Relationships Between Fungal Community Composition, Plant Genotypic Diversity and Environmental Drivers in a Salix Biomass System

Soil fungi are strongly affected by plant species or genotypes since plants modify their surrounding environment, but the effects of plant genotype diversity on fungal diversity and function have not been extensively studied. The interactive responses of fungal community composition to plant genotypic diversity and environmental drivers were investigated in Salix biomass systems, posing questions about: (1) How fungal diversity varies as a function of plant genotype diversity; (2) If plant genotype identity is a strong driver of fungal community composition also in plant mixtures; (3) How the fungal communities change through time (seasonally and interannually)?; and (4) Will the proportion of ECM fungi increase over the rotation? Soil samples were collected over 4 years, starting preplanting from two Salix field trials, including four genotypes with contrasting phenology and functional traits, and genotypes were grown in all possible combinations (four genotypes in Uppsala, Sweden, two in Rostock, Germany). Fungal communities were identified, using Pacific Biosciences sequencing of fungal ITS2 amplicons. We found some site-dependent relationships between fungal community composition and genotype or diversity level, and site accounted for the largest part of the variation in fungal community composition. Rostock had a more homogenous community structure, with significant effects of genotype, diversity level, and the presence of one genotype (“Loden”) on fungal community composition. Soil properties and plant and litter traits contributed to explaining the variation in fungal species composition. The within-season variation in composition was of a similar magnitude to the year-to-year variation. The proportion of ECM fungi increased over time irrespective of plant genotype diversity, and, in Uppsala, the 4-mixture showed a weaker response than other combinations. Species richness was generally higher in Uppsala compared with that in Rostock and increased over time, but did not increase with plant genotype diversity. This significant site-specificity underlines the need for consideration of diverse sites to draw general conclusions of temporal variations and functioning of fungal communities. A significant increase in ECM colonization of soil under the pioneer tree Salix on agricultural soils was evident and points to changed litter decomposition and soil carbon dynamics during Salix growth.

Soil fungi are strongly affected by plant species or genotypes since plants modify their surrounding environment, but the effects of plant genotype diversity on fungal diversity and function have not been extensively studied. The interactive responses of fungal community composition to plant genotypic diversity and environmental drivers were investigated in Salix biomass systems, posing questions about: (1) How fungal diversity varies as a function of plant genotype diversity; (2) If plant genotype identity is a strong driver of fungal community composition also in plant mixtures; (3) How the fungal communities change through time (seasonally and interannually)?; and (4) Will the proportion of ECM fungi increase over the rotation? Soil samples were collected over 4 years, starting preplanting from two Salix field trials, including four genotypes with contrasting phenology and functional traits, and genotypes were grown in all possible combinations (four genotypes in Uppsala, Sweden, two in Rostock, Germany). Fungal communities were identified, using Pacific Biosciences sequencing of fungal ITS2 amplicons. We found some site-dependent relationships between fungal community composition and genotype or diversity level, and site accounted for the largest part of the variation in fungal community composition. Rostock had a more homogenous community structure, with significant effects of genotype, diversity level, and the presence of one genotype ("Loden") on fungal community composition. Soil properties and plant and litter traits contributed to explaining the variation in fungal species composition. The within-season variation in composition was of a similar magnitude to the year-to-year variation. The proportion of ECM fungi increased over time irrespective of plant genotype diversity, and, in Uppsala, the 4-mixture showed a weaker response than other combinations. Species richness was generally higher in Uppsala compared with that in Rostock and increased over time, but did not increase with plant genotype diversity. This significant site-specificity underlines the need for consideration of diverse sites to draw general conclusions of temporal variations and functioning of fungal

INTRODUCTION
Fungi are principal drivers of biogeochemical cycling, linking above-and belowground ecosystem components through their roles as symbiotic mediators of plant nutrient uptake (Smith and Read, 2008), and as decomposers of litter (Schneider et al., 2012). Soil fungal communities are strongly influenced by plant genotypes or species (Prescott and Grayston, 2013;Gallart et al., 2018b), since plants modify their surrounding environment (De Deyn et al., 2008), for example by providing living and dead organic matter, and by feeding soil microbial communities via rhizodeposition (Steinauer et al., 2016). All these plant-derived substrates are likely to be affected by plant genotypes, species, and diversity (Weih and Nordh, 2002;Hoeber et al., 2020). The relationships between plant diversity and ecosystem functions are complex (Gamfeldt et al., 2013;van der Plas et al., 2016;Baeten et al., 2019), but because of the interactions between fungi and plants, these relationships are likely mediated and possibly explained by the composition and function of soil fungi. In this contribution, we explore the links between fungal communities and plant genotype diversity.
Most studies of plant genotypes impact on fungal communities have focused on differences between individual genotypes, and very few on genotype mixture or diversity effects (e.g., Baum et al., 2018). Host genotypes may be crucial to structuring ectomycorrhizal (ECM) fungal communities, both compositional (Püttsepp et al., 2004;Korkama et al., 2006;Hrynkiewicz et al., 2010;Velmala et al., 2013;Gehring et al., 2017) and structural (Lang et al., 2013), and may affect the ability to form mycorrhiza (Tagu et al., 2005;Gherghel et al., 2014). However, no support for host genotype structuring ECM assemblages has been reported (Bubner et al., 2013;Lang et al., 2013). A plant genotype can also affect other non-symbiotic root-associated fungal communities Gallart et al., 2018b;Bonito et al., 2019). Many previous studies were conducted on seedlings or younger trees; thus, the extent to which variation in tree genetics influences the structure and function of fungal communities within populations of adult trees remains an open question.
In addition to the fungal compositional changes, a plant genotype also affects the functional outcomes of the observed community shifts. Utilizing natural genetic variation in Populus hybrid populations, Schweitzer et al. (2011) tested whether stand gene diversity structures soil microbial communities, finding correlations between tree gene diversity, plant secondary chemistry and microbial community composition, impacting soil nitrogen availability. These results demonstrated that the effects of plant genetic diversity on other organisms may be mediated by the variation in plant functional traits (Schweitzer et al., 2011). In another study, rhizosphere and root-associated fungal communities in P. radiata differed in their response to both tree genotype and organic or inorganic nitrogen additions, suggesting that microbial communities more closely associated with roots were more sensitive to genotype-specific responses (Gallart et al., 2018a). Furthermore, fine root growth and soil and rootassociated fungal abundance and activity changed under host genotype mixtures caused by changed competitive conditions for individual plant genotypes (Elferjani et al., 2014;Baum et al., 2018). Since soil fungal communities are central to plant nutrient supply, plant biomass production, and litter decomposition, they can significantly affect the carbon (C) storage in plant biomass and soil organic matter. Therefore, analyses of the impact of plant genotype diversity on soil fungal diversity can contribute to develop a rationale for the definition of genotypespecific optimization of tree diversity to combine high biomass productivity with high total C storage.
Tree species or genotype identity can sometimes have stronger impacts on ecosystem processes and organisms than diversity per se (Scherer-Lorenzen et al., 2005). The presence of a given genotype with effects on fungal communities can thus determine soil fungal community responses to genotype mixtures at different diversity levels. However, there are also studies that have shown no or little effect of host identity on fungal communities associated with willows (Erlandson et al., 2016(Erlandson et al., , 2018Arraiano-Castilho et al., 2020). Apart from genotype, fungal communities are strongly driven by environmental factors, including those changing during the growing season; especially soil pH and nutrient status (Lauber et al., 2008;Tedersoo et al., 2014). Several studies have shown that, for a given plant species, soil properties have a stronger effect on the soil and root-associated fungal community composition than the host genotype (Karlinski et al., 2013;Cregger et al., 2018;Bonito et al., 2019). Additionally, temporal trends in fungal community composition in relation to plant host genotype occur, both between and within seasons. Although less is known about seasonal variation, seasonal differences in the composition of, e.g., saprotrophic microfungi associated with roots of two Salix clones were observed previously by Baum and Hrynkiewicz (2006). Due to the dependence of soil fungal communities on climatic and edaphic factors, as well as plant community composition and/or the tree host, assessment of the interactive responses to biotic and abiotic factors has been put forward as a major challenge in fungal ecology (van der Heijden et al., 2015).
We addressed this challenge, using fast-growing Salix genotypes as a case study. Salix genotypes are cultivated in a short rotation coppice to produce biomass for energy purposes (Weih, 2004), and the cultivation of Salix is considered a sustainable biomass source with a positive greenhouse gas balance (Cunniff et al., 2015). Instead of planting single Salix genotypes, increased plant diversity has been proposed to be advantageous in terms of decreased pest, disease, and abiotic stress, and increased plant productivity and nutrient status (McCracken et al., 2011). Due to the fast growth and great variability among genotypes or closely related species with different plant traits, Salix is a suitable woody perennial model system to test for effects of both genotype and environmental factors on different ecosystem processes (Weih et al., 2019). Salix is a dual mycorrhizal plant genus, forming ECM and arbuscular mycorrhizal (AM) associations, with the two types of fungi colonizing roots with an antagonistic behavior (Lodge and Wentworth, 1990;Baum et al., 2018). In northern Europe, species of Salix have been shown to predominantly form ECM associations (Püttsepp et al., 2004). We used a short-rotation coppice system, where four genotypes of Salix were grown in all possible combinations, offering an opportunity to study the effect of plant genotypic diversity on soil fungal community composition. Previous studies from those trials showed genotype identity to have a more important influence than genotype diversity on shoot biomass productivity , and leaf litter quality coupled to genotype drive litter decomposition more than stand diversity or climate (Hoeber et al., 2020). The trials were established on former arable land  where we expect ECM fungal colonization to be sparse at the establishment and increase over time (see Kalucka and Jagodzinski, 2016). The aim of the present study was to describe the soil fungal communities, posing the following research questions with a focus on the potential drivers of soil fungal community structure and diversity: (1) How does fungal diversity vary as a function of plant diversity; (2) Is plant genotype identity a strong driver of fungal community composition even in tree mixtures?; (3) How does the fungal community change through time (seasonally and interannually)?; and (4) Will the proportion of ECM fungi within the total fungal community increase over the rotation, from the establishment on arable land and up until the first biomass harvest? Finally, the significance of plant genotype diversity on the soil fungal community composition was tested in relation to other drivers (plant biomass, litter chemistry, litter decomposition, soil properties, fungal biomass).

Study Sites
Two field trials were established on arable land in May, 2014, in Uppsala, Sweden (59 • 49' N 17 • 39' E) and Rostock, Northern Germany (54 • 02' N 12 • 05' E) within the ECOLINK-Salix trial , which is part of the TreeDivNet global tree diversity network . Four genotypes of Salix, partly differing at species or intraspecific levels, were used as the stand components: "Björn" (Salix schwerinii E. Wolf. × S. viminalis L.), "Jorr" (S. viminalis), "Loden" (S. dasycladosWimm.), and "Tora" (S. schwerinii × S. viminalis). "Björn" and "Tora" are taxonomically closely related as they derive from the same parent material and are siblings but differ in some morphological and functional characteristics. "Jorr" is more closely related to the siblings "Björn" and "Tora, " whereas "Loden" belongs to a different species and thus is taxonomically most distant from the other three varieties. The genotypes were planted in pure cultures and various mixtures (2-, 3-, and 4-mixtures, giving four different plant diversity levels) and arranged in a randomized block design with a total of 45 plots over three blocks in Uppsala (15 plots per block). In Rostock, only the two genotypes "Loden" and "Tora" were grown, resulting in a total of nine plots. Plots measured 9.6 × 9.6 m and contained 12 rows with 12 plants in each row with offset every second row, resulting in a hexagonal planting pattern with equal distances of 0.8 m between individuals. This spacing corresponds to approximately 15,600 plants ha −1 . Further details about the establishment of the sites, soil, and climatic conditions can be found in Hoeber et al. (2018). Soil temperature and precipitation for the experimental period at both sites are reported in Supplementary Figure 1.

Soil Sampling
For molecular analysis, soil cores (3 × 10 cm) were sampled in April 2014, 2015, and 2017, and in April, August, and November 2016, covering one short rotation coppice cycle from planting to harvesting (2014-2016), one sampling point after (2017), and seasonal variations for 1 year (2016). In each plot, a total of nine soil cores were collected and pooled, generating a total of 270 samples from Uppsala (45 samples per sampling occasion), and 54 samples from Rostock (nine samples per sampling occasion). Immediately after sampling, soil cores were either stored in a freezer and later freeze dried (Uppsala), or stored in the freezer and oven-dried at 45 • C for 48 h (Rostock); after which, each dried sample was homogenized, using a pestle and mortar.

Quantification of Fungal and Root Biomass
For soil samples in 2017, total fungal biomass was quantified from 0.3 g of a soil sample, using the fungal biomarker ergosterol. Ergosterol is the most common sterol of Ascomycota and Basidiomycota. Using established methods (Nylund and Wallander, 1992), esterified ergosterol was extracted with 10% KOH in MeOH, filtered through a 0.45-µm Teflon filter, and 50 µl of each sample was analyzed, using high performance liquid chromatography, with a C18 reverse-phase column (Nova-Pak; 3.9 × 150 mm; Waters, Milford, CT, USA), preceded by a C18 reverse-phase guard column (Nova-Pak; 3.9 × 20 mm; Waters) as in Cheeke et al. (2017). The ergosterol peak was detected at 282 nm, using a UV detector. Fungal biomass was determined from ergosterol concentrations, using a conversion factor of a 3µg ergosterol mg −1 dry sample (Salmanowicz and Nylund, 1988) and a correction factor (1/0.62) to compensate for unextracted ergosterol (Montgomery et al., 2000).

Molecular Analysis
DNA was extracted from a subsample of approximately 400 µl by volume from each homogenized sample, using Macherey-Nagel NucleoSpin Soil kit (Düren, Germany), with the following modification; lysis buffer SL2 700 µl with 150 µl enhancer and quantified using a Thermo NanoDrop spectrophotometer (Wilmington, DE, USA). The internal transcribed spacer (ITS) region is the universal barcode for fungi, and ITS2 amplicons were produced, using the forward primer gITS7 (GTGARTCATCGARTCTTTG; (Ihrmark et al., 2012) and the two mixed reverse primers ITS4 (75%; 5 ′ -TCCTCCGCTTAT TGATATGC-3 ′ ; (White et al., 1990) and ITS4arch (25%; 5 ′ -CACACGCTGTCCTCGCCTTATTGATATGC-3 ′ ), elongated with unique identification tags (Clemmensen et al., 2016). PCR was performed in 50 ml reactions with 25 µl DNA template (diluted to 0.5 ng/µl), 0.2 mM of each nucleotide, 0.75 mM MgCl 2 , forward primer at 0.5 µM, reverse primer at 0.3 µM, and 0.5 U polymerase (DreamTaq, Thermo Scientific, MA, USA) in PCR buffer on a 2720 Thermal Cycler (Life Technologies, California, USA). PCR conditions were 5 min at 94 • C, 25-35 cycles of (30 s at 95 • C, 30 s at 56 • C, 30 s at 72 • C) and 7 min at 72 • C, and cycle numbers were adapted for each sample to give similar band strength to avoid oversaturation of the PCR pool. Two pooled PCR replicates from each sample were purified, using the AMPure kit (Beckman Coulter, California, USA) according to the instructions of the manufacturer, and quantified, using a Qubit R 2.0 fluorometer (Invitrogen, Carlsbad, California, USA). The products were mixed in equal amounts into four pools and cleaned, using the E.Z.N.A. R Cycle Pure Kit (Omega biotek, Norcross, Georgia, USA). Adaptor ligation and sequencing on PacBio Sequel instrument (Pacific Biosciences, Menlo park, California, USA) were performed by NGI-Uppsala/SciLifeLab (National Genomics Infrastructure, Uppsala, Sweden), using six Sequel SMRT cells (v2; pool 1) and two Sequel SMRT cells (v3; pools 2-4).

Sequence Analyses
Raw sequence reads were analyzed, using the bioinformatics pipeline SCATA (https://scata.mykopat.slu.se) (Ihrmark et al., 2012). Sequences were quality filtered and screened for primers and identification tags as described in Kyaschenko et al. (2017) with some adjustments. After removal of sequences with mean quality of 20 bases and lower or containing bases with a quality of 3, sequences (complementary reversed, if needed) were searched for primers and identification tags. Only sequences containing matching tags at both ends were retained. All sequences were clustered into study-level species hypotheses (SHs; Koljalg et al., 2013), using a 1.2% threshold distance for sequences to enter an SH. A reference database from UNITE (UNITE community general FASTA release, version 7.2, release date, 10.10.2017) was included in the clustering. Sequences are available in the NCBI Sequence Read Archive (www.ncbi.nlm.nih.gov/sra) under the accession number PRJNA703824. To remove non-fungal sequences, the representative sequences from each SH were compared against GenBank nucleotide database using BLASTn, after which we used MEGAN (Huson et al., 2011) for the BLAST results and fasta file to assign the lowest common ancestor and to identify sequences that were not fungal (Balint et al., 2014). All singletons and the small clusters with five or less counts were removed. The representative sequences from each study-level SH were then compared to all global SHs, using the massBLASTer through the PlutoF platform in UNITE (Abarenkov et al., 2010) and assigned to appropriate a taxonomic level (at least 97% similarity was required for species-level identification, 90% for genus, 85% for family, 80% for order, 75% for class, and 70% for division/phylum), in order of decreasing global relative abundance until 80% of the sequences were covered. To ensure that the difference in sequencing chemistry (pool 1 vs. other pools) did not affect the outcome of community compositional analyses due to different sequencing depth, both the original study-level SHs table for all clusters and counts data, and the identified study-level SHs, were rarefied (100 random subsampling) to 320 and 168 sequences per sample, respectively. Samples with less sequences were removed (eight samples for all study-level SHs, 10 for the identified). Additionally, we applied centered log-ratio (clr) transformations (Gloor et al., 2017;Quinn et al., 2018) to the original study-level SHs table to obtain values that were scale-invariant and to account for differences in count numbers. Study-level SHs were further assigned to ecological functions, using FUNGuild software (Nguyen et al., 2016). In the functional group bar charts, all available functional groups (guilds) from FUNGuild were included with some modifications: SHs classified as fungal parasite-saprotroph and fungal parasiteundefined saprotroph were called fungal parasite/saprotroph; fungal parasite-protistan parasite was called putative fungal parasite; plant pathogen-plant saprotroph/undefined saprotroph/wood saprotroph combinations were called plant pathogen/saprotroph; endophyte-litter saprotroph/soil saprotroph/undefined saprotroph/wood saprotroph were called putative endophyte/saprotroph; combinations of different saprotrophic classification were called saprotrophs, and any other SHs with a functional classification with several possibilities were assigned as unknown. See Supplementary Table 1 for functional groups and exploration type assignments (Agerer, 2001(Agerer, , 2006Tedersoo et al., 2006;Lilleskov et al., 2011;Katanic et al., 2014), and Supplementary Table 2 for relative abundances of the 446 most common identified SHs (rarefied and not rarefied data). In multivariate plots, showing the most common SHs (see below), functional groups were simplified, using ECM, saprotroph, plant pathogen, animal pathogen, and unknown.

Statistical Analyses
Ordination analyses were performed using CANOCO version 5.02 (Microcomputer Power, Ithaca, NY, USA). Variation in soil fungal community composition (total fungal community; 446 SHs, and ECM fungal community; 18 SHs, see Supplementary Tables 2, 3) was visualized, using principal components analysis (PCA) and detrended correspondence analysis (DCA). Separate multivariate analyses were done for Rostock and Uppsala, as well as Rostock and Uppsala together, based on the individual samples, using both the rarefied (reported in the results) and the original data (supplemental data in Supplementary Material). In order to further control for differences in count numbers between samples, we also did multivariate analyses using clr-transformed data (not reported). The PCAs and DCAs reported were based on identified communities resolved at the SHs level, but additional analyses (not reported) were done with total fungal communities, including all study-level SHs to confirm that they produced the same pattern as the subset of identified SHs. Correlation between Soil sampling for fungal biomass was done in April 2017, at the end of the first short rotation coppice cycle. a Genotypes correspond to A = "Björn", B = "Jorr", C = "Loden", and D = "Tora". b Fungal biomass was estimated from esterified ergosterol concentrations.
plant, fungal, and soil parameters with site, block, genotypes, and diversity levels, and their interaction with the identified soil fungal community (rarefied data) were tested, using redundancy analysis (RDA). We used combinations of previously published data from the same field trials on plant biomass ( Table 4) and fungal biomass from 2017 (Table 1). These potential drivers were tested against the fungal community data for the corresponding sampling time.
Simple term effects are reported after Bonferroni corrections. Diversity measurements (species richness; i.e., total counts of SHs, Shannon-Wiener index and Simpson's index of diversity 1-D) based on the rarefied data for all study-level SHs are reported both as averages (± SD) per treatment (Supplementary Table 5; site × year × harvest time × diversity level × genotype × block) and for all individual samples (Supplementary Table 6). Species data were arcsine transformed before all multivariate analyses, with the exception of the identified SHs from Uppsala (original data, not rarefied), where species data were log transformed. Rare SHs were down-weighted in the DCA of the total fungal community analysis. We also used the multiresponse permutation procedure (MRPP), a non-parametric procedure in PC-ORD version 5.33 software (McCune, 2006) for testing the hypothesis of no difference between two or more a priori assigned groups (McCune and Grace, 2002). This was done to test for the effects of genotype combination, presence of one genotype ("Björn, " "Jorr, " "Loden" or "Tora"), diversity level (monoculture, 2-, 3-, and 4-mixture), year (2014, 2015, 2016, and 2017), harvest time (April, August, and November), year × harvest time, and block on compositional data for Uppsala and Rostock (total and identified community, rarefied data) separately. The effects of site, diversity level, year, genotype combination, block, and harvest time on fungal biomass and diversity were tested, using R (version 3.6.1, R Core Team, 2019). Multiple comparison analyses were performed, using two-way ANOVA, followed by Tukey's HSD test. To summarize general temporal trends in community composition, functional groups were aggregated into two broad root-associated and saprotroph categories by summing up abundances as follows: root-associated included arbuscular mycorrhizal (two clusters), ectomycorrhizal, and endophyte (one cluster) groups; saprotrophs included endophyte-undefined saprotroph-wood saprotroph, all plant saprotroph, undefined saprotroph, and wood saprotroph groups. The abundances in these two groups were then normalized by the total abundance to calculate proportions of root-associated fungi and saprotrophs. Temporal trends were assessed for the rootassociated and saprotroph proportions and the root-associated: saprotroph ratio, using a linear mixed effect model (fitlme function in Matlab R2018b, The Mathworks) initially including year, month, and number of Salix genotypes as fixed effects and the site as random effect. Since the number of genotypes had no predictive power, it was removed from the regression.

Fungal Biomass
Fungal biomass for the genotypes "Loden" and "Tora" and their mixture present at both sites were significantly higher in Uppsala (1. ±0.02 mg g −1 soil) compared to Rostock (0.8 ±0.02 mg g −1 soil) ( Table 1 and Supplementary Table 7). There were no differences in fungal biomass between different plant genotypes or diversity levels for either site.

Sequencing Output
Sequencing generated a total of 2,294,094 reads, of which 1,340,101 reads passed quality control (QC). After removal of non-fungal sequences (427,706 reads, 32% of reads passed QC) and singletons (23,483 reads), the remaining 888,912 reads clustered into 4,839 study-level SHs. Clusters with five or fewer reads (2,400 clusters and 6,568 reads) were removed, resulting in a total of 2,421 study-level SHs; out of which, 446 SHs (80% of the fungal sequences) were subjected to taxonomic and functional identification (Supplementary Tables 1, 2). Each sample had an average of 287 reads (3,023 maximum, 200 minimum). After rarefying the dataset, 1,943 and 446 study-level SHs remained for all and identified clusters, respectively. Although the number of sequences per sample varied greatly (Supplementary Figure 2), it was clear from the multivariate analyses that the rarefied and the original dataset (Figure 1, Supplementary Figures 3, 4), as FIGURE 1 | Variation in soil fungal community composition in Salix genotype trials with the four genotypes "Björn," "Jorr," "Loden," and "Tora" planted in monoculture and various genotype mixtures (2-, 3-, and 4-mixtures) at two different sites (Rostock, Germany, and Uppsala, Sweden). Community (Continued) well as the clr-transformed dataset (not reported), showed very similar community compositional patterns.

Fungal Functional Groups and Diversity
Uppsala and Rostock were colonized by distinct soil fungal communities (  Table 2). For the functional groups, the proportion of ECM fungi increased significantly over time (from preplanting until the first short rotation coppice cycle ended) irrespective of the plant diversity level (monoculture, 2-, 3-, or 4-mixture), and was followed by a significant decrease in the proportion of saprotrophic fungi (Figure 2, Table 2). ECM fungi increased about 6-fold in Uppsala and 70-fold in Rostock 3 years after planting at both sites ( Figure 3A), but this large increase in Uppsala was only observed in monocultures 2-and 3-mixtures, while the 4-mixture showed a weaker response ( Figure 3C). The 2-mixture of genotypes "Loden" and "Tora" in Rostock showed an increase of ECM fungi by ca 55% compared to the monocultures ( Figure 3B). Fungal diversity (Supplementary Tables 5, 6) did not differ depending on the plant diversity level or genotype mixture (Supplementary Table 8). Species richness was significantly higher in Uppsala compared to Rostock (Figure 1B,  Supplementary Tables 5, 6, 8), with an average of 94 ± 1 and 69 ± 2 SHs, respectively, taken over all treatments and times. Although fungal species richness was variable between years within individual genotype treatments, there was an increase in species diversity over time for some of the genotypes for both sites (Figure 1, Supplementary Tables 5, 8).

Potential Drivers of Community Responses
The overall community composition was significantly related to site, soil properties of pre-planted fields, fungal biomass in soil, plant shoot biomass, some of the leaf chemistry variables, and some of the litter decomposition variables (see Figure 4 for analyses of 2014 and 2017 data, 2015 and 2016 not shown). The potential drivers together explained 20-40% of the variation in fungal community composition. Site was always the explanatory variable, accounting for most of the variation in fungal community composition (p ≤ 0.026 up to 13.5%; Figure 4). Soil properties of pre-planted fields significantly explained part of the variation (p = 0.018, ca 5% of the variation for each variable) for the fungal community in 2014 ( Figure 4A) and the following years, although explaining less of the variation over time (not shown). Leaf N and P concentrations after one growing season were the only leaf chemical characteristics that weakly explained community composition in 2015 (p = 0.032 TABLE 2 | Summary of results from linear mixed effect models to predict root-associated and saprotrophic contributions to the total sequence counts, and the ratio of root-associated over saprotrophic counts.

Year-to-Year and Seasonal Variation in Community Composition and Diversity
For both sites, the fungal community changed over time; in Uppsala ( Figure 5 and Supplementary Figure 5), 1 year after planting the composition was similar to pre-planting and dominated by SHs in Ascomycota (Helotiales, Torrubiella sp., Ciliolarinia sp., Crocicreas furvum, and Penicillium sp.); after which, other ascomycetes became more common alongside with Mortierella sp. and basidiomycetes (Inocybe curviceps and Hebeloma alpinum). In Rostock, the community in pre-planted soil was separated from planted and developed over time (Figure 6 and Supplementary Figure 6) but with a more similar community structure compared to Uppsala. Pre-planting, SHs beloning predominantly to Ascomycota (Helotiales, Cladoaporium sp., and Fusarium sp.; Figure 6B) were most common. After planting, ECM fungi (Hebeloma alpinum, Hebeloma mesophaeum, and Laccariatortilis) and other basidiomycetes (Lentaria afflata and Clitopilus hobsonii) became more common (Figure 6B). The within-season variation was of a similar magnitude as the year-to-year variation in Uppsala ( Figure 5A) and was also apparent in Rostock, with August showing the largest spread in composition ( Figure 6A). Species richness (Supplementary Table 8) and the proportion of ECM fungi (Supplementary Figure 7) were highest in November in Uppsala. In Rostock, the proportion of ECM fungi varied strongly between years, and after the third growing season ECM fungi contributed almost 25% to the total fungal community ( Figure 3A). There were significant effects of year (MRPP Uppsala: p < 10 −9 , A = 0.11; Rostock: p < 10 −9 , A = 0.10), and block (Uppsala: p < 10 −9 , A = 0.0074) on community composition of a priori assigned groups. Block 3 diverged somewhat in community composition from other FIGURE 3 | Distribution of fungal functional groups in soil from Salix genotype trials in Rostock, Germany, and Uppsala, Sweden: (A) over time (preplanting in 2014 to 2017) (B) in monocultures of "Loden" or "Tora" and in 2-mixtures of "Loden" and "Tora", and (C) for Uppsala only, in monoculture, 2-, 3-, and 4-mixtures, as estimated by PacBio sequencing of amplified ITS2 markers. Abundances are given as percent of the identified amplicon sequences, based on rarefied data (accounting for 80% of total sequences). There was a large increase in the proportion of ECM fungi three years after planting at both sites, especially in monoculture, 2-and 3-mixtures in Uppsala but with a smaller increase in 4-mixtures.
blocks (Supplementary Figure 8). For the ECM community in Uppsala, development over time was similar to the fungal community, but less pronounced, and ECM species richness also increased over time (Supplementary Figure 9). Specifically, species richness was significantly higher in 2017 compared with both 2015 and 2016 (Supplementary Table 8).

DISCUSSION
In this study, we investigated the interactive responses of soil fungal community composition to plant genotypic diversity and environmental drivers in a Salix biomass system. Although there was no overall effect of genotype, we found some site-dependent relationships between fungal community composition and plant genotypic diversity, and site accounted for the largest part of the variation in fungal community composition. In Rostock, where the fungal community structure was more homogenous compared to in Uppsala, there were significant effects of genotype, diversity level, and the presence of one genotype ("Loden") on fungal community composition. The proportion of ECM fungi increased over time, and, in Uppsala, the 4-mixture showed a weaker response than other genotype combinations.
Overall, fungal communities were thus site-dependent in relation to plant community diversity. However, the ECOLINK-Salix trials would need to be studied longer and also include more genotypes, since potential genotype effects on plant, microbial, and soil parameters may become apparent in later short rotation coppice cycles.

Few Effects of Plant Genotype Diversity on Soil Fungal Diversity
For answering our first research question about how fungal diversity varies as a function of plant diversity, we found some support for plant genotype effects in the site-specific analyses where fungal diversity increased in the 2-mixture of "Loden" and "Tora" compared to the monocultures in Rostock (MRPP analysis), suggesting that the plant genotype diversity effect is site or context dependent Bonito et al., 2019). A common garden experiment, using a Populus model system, found a host genotype to explain as much as 70% of the variation in microbial community composition based on PLFAs (Schweitzer et al., 2008); however, this large difference in plant genotype effects on microbial communities was most likely explained by minimization of site bias and the use of the method to characterize the microbial community. In the present study, both sites were characterized by distinct soil fungal communities with a high number of SHs, comparable to what previous community studies have reported in terms of composition, dominating taxa and mean number of taxa (Perez-Izquierdo et al., 2017;Pinus pinaster), although Rostock had significantly lower species richness compared with Uppsala.
Since soil microbial diversity in general is linked to the clay content (Xue et al., 2018), which is much higher in Uppsala compared with Rostock (Hoeber et al., 2020), we would still expect higher fungal diversity in Uppsala (Essene et al., 2017). A low ECM and endophyte diversity in roots in all genotype treatments in the Rostock trial were also shown previously , but the mycorrhizal root tips did not show a strong increase in ECM colonization in the 2-mixture compared to monocultures as the increase in the proportion ECM fungi in soil in the Rostock 2-mixture in the present study. In Uppsala, site conditions were more heterogeneous (block effects). The large increase in the proportion of ECM fungi over time was confined to monocultures, 2-and 3-mixtures, and much lower in the 4-mixtures, suggesting that there could be an "optimal" genotype mixture or a plant diversity level for sustaining ECM fungal communities. Similarly, the 4-mixture tended to exhibit slower litter decomposition compared to monocultures, 2-and 3-mixtures, in Uppsala and Freiburg (the third ECOLINK site, not included in the present study; Hoeber et al., 2020). Tree productivity results for Freiburg also showed that shoot biomass production was more variable and lower in the 4-mixture compared with the less diverse stands . These findings also corroborate previous studies, showing that higher tree or genotype diversity can have diverse impacts on fungal diversity and composition (e.g., Baum et al., 2018;Bonito et al., 2019).

Plant Genotype Identity and Other Drivers of Fungal Communities
Plant genotype identity was not a strong driver of fungal community composition in line with previous studies (Erlandson et al., 2016(Erlandson et al., , 2018Arraiano-Castilho et al., 2020) -the only evidence was the significant effect of "Loden" on fungal community composition in Rostock. The different genotypes in the ECOLINK-Salix trial partly belong to different Salix species and have contrasting phenology, functional traits, and ecophysiology (Weih and Nordh, 2002;Weih et al., 2021). Although increasing genotype diversity so far did not significantly affect aboveground plant productivity in young stands, admixing of two genotypes ("Jorr" and "Loden") was predicted to enhance the total shoot biomass production in a mixed stand, while two other genotypes ("Tora" and "Björn") were more likely to reduce the total shoot biomass production in mixtures . Hence, we expected to find similar effects on soil fungal communities, because of the expected links between productivity and belowground communities. The effect on fungal communities of having "Loden" present (in Rostock) corroborates the potential beneficial effect of "Loden" on biomass production in mixed stands, since a larger tree biomass production in stands in which "Loden" (as opposed to "Tora") is admixed may possibly be linked to more belowground C, supporting soil fungal communities. This possible mechanistic link would, however, require verification by means of further experimental work.
Since the relative importance of interactions between stand characteristics, soil properties, and climatic conditions in relation to tree genotype diversity or identity is hardly known so far, we also tested the significance of plant genotypic diversity on the soil fungal community composition in relation to other drivers (soil properties, plant biomass, litter chemistry, litter decomposition, and fungal biomass). Although several of the plant derived properties and soil properties were significantly correlated with soil fungal communities, we found no significant overall effects of plant genotype mixtures. These results support the idea that the different drivers acted in concert dependent on site -but not always via plant genotype differences -to shift fungal community composition and diversity in the present study. This contrasts with the idea that genotypic variation in aboveground plant traits is important to soil microbial dynamics and nutrient availability (Schweitzer et al., 2004). ECM fungal communities have, for example, been linked to plant biomass (growth and size of host plants; Korkama et al., 2006;Velmala et al., 2013) and leaf chemistry (ECM community covary with senescent leaf chemistry; Lamit et al., 2016). One study also found no plant genotype effect, reporting on diverse and evenly distributed soil fungal communities in different Populus genotypes (Danielsen et al., 2012). Soil chemistry also can be a strong driver of fungal communities (Gallart et al., 2018b;Bonito et al., 2019), especially nutrient availability and soil pH.

Year-to-Year Variation in Fungal Community Composition as Large as Seasonal Variation
We expected the variation in fungal community composition within the season to be larger than the year-to-year variation, since within-season variation in environmental conditions and plant inputs is often higher than year-to-year changes. Instead, the within-season variation was of a similar magnitude as the year-to-year variation, although community composition in Uppsala varied more than in Rostock, possibly due to a wider seasonal temperature range and greater spatial heterogeneity. Large within-season variation in fungal community composition has been frequently reported (Voriskova et al., 2014;Haas et al., 2018). For example, Baum and Hrynkiewicz (2006) found a seasonal shift in saprotrophic fungi and enzymatic activities associated with two Salix genotypes, and, although community composition shifted over the growing season, diversity remained larger for one of the genotypes. Variation in the root microbiome of Populus trees was explained by season and soil properties (Shakya et al., 2013). Therefore, interannual variation seems to have less importance compared with seasonal variation on timescales shorter than successional scales, e.g., for ECM fungi (Bahram et al., 2015, and references therein). However, over longer timescales, the fungal communities undergo significant successional changes. Indeed, we found an increase in both fungal diversity and proportion of ECM fungi over time. An increase in the proportion of ECM fungi, after establishing mycorrhizal Salix on arable land, would be expected (Dickie et al., 2013), and is analogous to forest regrowth after clearcut (see, e.g., Kyaschenko et al., 2017). As plants grow, they supply the surrounding soil with C, influencing soil chemistry and biota (Leake and Read, 2017), and larger trees have the capacity to support more symbiotic fungi. Since tree biomass was larger in Rostock compared with Uppsala (mean shoot biomass production over the whole cycle after the first cutting, 8.7 Mg ha −1 and 5.7 Mg ha −1 in Rostock and Uppsala, respectively), this difference in tree biomass likely explains the higher proportion of ECM fungi found in Rostock.

Methodological Considerations
The potential problem with differences in sequencing depths due to sequencing chemistry (amplicon pool 1: 276,111 sequences, 25% passing QC; amplicon pools 2-4: mean number of sequences 672,661, 63% passing QC; see Materials and Methods) was handled by rarefying data. AM fungi were only detected infrequently in our study, even though Salix was planted on a former agricultural field where they should be present and even promoted by the no-till management of Salix. Low AM fungal colonization of Salix has been described by several authors (e.g., Püttsepp et al., 2004;Hrynkiewicz et al., 2012). Furthermore, low detection of AM fungi may also be due to low abundance (in terms of DNA) relative to other fungal taxa, or because the analytical methods used to quantify fungi (ITS2 primers) were not optimal for this fungal group (e.g., Delavaux et al., 2020). Although these primers work well for amplifying the fungal ITS2 region for Ascomycota and Basidiomycota (Ihrmark et al., 2012), they are less frequently used in AM community studies (Glomeromycota) due to poor amplification (Tedersoo et al., 2018). Furthermore, sampling soil instead of mycorrhizal roots may be masking direct genotype effects on root-associated fungi. It should also be stressed that our results are correlative and that manipulative experiments would be required to consolidate our conclusions and confirm causal relationships.

CONCLUSIONS
Salix spp. in short rotation coppice was here successfully used as model system to explain environmental changes in soil ecological properties, especially fungal diversity. Fungal communities were largely independent of plant community diversity but were affected mainly by site-specific conditions. This significant site specificity underlines the need of consideration of diverse sites to draw general conclusions on temporal variations and functioning of fungal communities. A significant increase in ECM colonization of soil under the pioneer tree Salix on agricultural soils was evident and points to a changed litter decomposition and soil C dynamics during Salix growth, where symbiotic nutrient exchange might be an important driver to consider.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
MW initiated the field experiment together with CB. SH and PF carried out field sampling. SH extracted soil DNA and prepared one sequencing pool. PF performed most sequence analyses and multivariate analyses. SH carried out other statistical analyses, using R. SM contributed to regression analysis. PF drafted the manuscript. SH, CB, MW, and SM contributed to data interpretations and article development, and approved the submitted version. All authors contributed to the article and approved the submitted version.

FUNDING
The research was partly funded by grants from the Swedish Energy Agency (36654-1 and 36654-2) to MW and CB. SM and MW acknowledge support from the Swedish Research Agency Formas (2016-00998).