Tree Diversity Reduces Fungal Endophyte Richness and Diversity in a Large-Scale Temperate Forest Experiment

Although decades of research have typically demonstrated a positive correlation between biodiversity of primary producers and associated trophic levels, the ecological drivers of this association are poorly understood. Recent evidence suggests that the plant microbiome, or the fungi and bacteria found on and inside plant hosts, may be cryptic yet important drivers of important processes, including primary production and trophic interactions. Here, using high-throughput sequencing, we characterized foliar fungal community diversity, composition, and function from 15 broadleaved tree species (N = 545) in a recently established, large-scale temperate tree diversity experiment using over 17,000 seedlings. Specifically, we tested whether increases in tree richness and phylogenetic diversity would increase fungal endophyte diversity (the “Diversity Begets Diversity” hypothesis), as well as alter community composition (the “Tree Diversity–Endophyte Community” hypothesis) and function (the “Tree Diversity–Endophyte Function” hypothesis) at different spatial scales. We demonstrated that increasing tree richness and phylogenetic diversity decreased fungal species and functional guild richness and diversity, including pathogens, saprotrophs, and parasites, within the first three years of a forest diversity experiment. These patterns were consistent at the neighborhood and tree plot scale. Our results suggest that fungal endophytes, unlike other trophic levels (e.g., herbivores as well as epiphytic bacteria), respond negatively to increasing plant diversity.


Introduction
Rapid environmental and anthropogenic changes are altering patterns of biodiversity worldwide, degrading ecosystem productivity, function, and services [1,2]. Increased temperatures and human activities over the last century have reduced forest habitats by 50%, which has subsequently diminished plant and animal diversity [3][4][5][6]. Declining biodiversity will likely have important implications for human wellbeing, as decades of research in land and seascapes have demonstrated a positive relationship between biodiversity and critical ecosystem services (e.g., increasing plant productivity, supporting diverse faunal communities, regulating nutrient cycles, and buffering ecosystems against climate change; reviewed by [7,8]). The "Diversity Begets Diversity" Hypothesis (H1): Increasing tree richness (plot and neighborhood level) and phylogenetic diversity (plot and neighborhood level) will cause an increase in tree host foliar fungal endophyte diversity and richness. The "Tree Diversity-Endophyte Community" Hypothesis (H2): Tree species richness (plot and neighborhood level) and phylogenetic diversity (plot and neighborhood level) will explain more variation in endophyte community composition than host identity alone. The "Tree Diversity-Endophyte Function" Hypothesis (H3): Increases in tree richness (plot level) will cause an increase in fungal functional guild richness and diversity.

Site Details
We conducted this study at the BiodiversiTREE experiment in Edgewater, MD, USA (38 • 52 N, 76 • 33 W, Figure 1). The BiodiversiTREE experiment began in 2013 on 13 ha of abandoned maize cropland as a part of the TreeDivNet network, which spans 26 tree diversity experiments and over 1.1 million trees globally [37,38]. The annual mean temperature at the site is 13.2 • C, and the average annual precipitation is 1068 mm. In the spring of 2013, we planted 17,850 one-year-old seedlings from 16 of the most locally abundant tree species in 35 × 35 m experimental plots (N = 70), with 255 equally spaced trees per plot (see Figure 1 for species list and layout). We selected 16 out of the 20 tree species with the highest basal area in Maryland forests [39,40]; G. Parker, unpublished data). Moreover, data from the surrounding CTFS-ForestGEO plots at the site demonstrate that our species were reflective of the surrounding forest matrix [39]. Eight of these tree species associate with arbuscular mycorrhizal fungi (henceforth AM species, Figure 1), and eight species associate with ectomycorrhizal fungi (henceforth ECM species, Figure 1). While some Quercus spp. can concurrently associate with both ectomycorrhizal (ECM) and arbuscular mycorrhizal (AM) fungi, the function of AM is thought to be insignificant for these hosts (e.g., [41,42]). McCormick characterized AM and ECM communities from the same focal species at the BiodiversiTREE site (McCormick, et al., unpublished data); thus, we can confirm that the trees have already established these mycorrhizal communities.
We obtained seedlings as bare-root whips (mean height of 38 cm) from a commercial nursery and planted them directly into the existing soil without fertilizer or mycorrhizal augmentation. Overall, the entire experimental array has relatively minor slope gradients (mean 4.8 • ± 1.93 • , full range = 0.03 • to 13.25 • ) and consistent land-use history (continuously planted with corn Zea mays L. for > 35 year), thus providing relatively homogenous environmental conditions across the site. Cladogram and type of mycorrhizal association (AM = arbuscular mycorrhizal-associated species, ECM = ectomycorrhizal-associated species) for the 16 focal tree species at the site. We used all species but Carya glabra for this experiment. (B) Experimental layout of plots at the site and focal tree sampling. Each plot contains 255 equally spaced trees per plot, which we planted in either monoculture (1 sp.), 4-spp. or 12-spp. plots, totaling 70 plots. We used a well-replicated design across the 15 tree species at the site, totaling 545 trees.

Experimental Design
The BiodiversiTREE experiment employs a classic BEF manipulation, whereby we planted tree species into either monoculture plots (N = 2 plots per species), four-species plots (N = 19), or 12-species plots (N = 19). Each monoculture plot received 255 individuals per species; each four-species plot received~63 individuals per species; the 12-species polycultures received~21 individuals per species. We generated mixed-species assemblages randomly, but ensured that each species occurred with equivalent frequencies across the entire experiment. We randomly assigned the spatial location of each individual tree in every plot by planting them in an equidistant hexagonal grid with 2.4 m between trees ( Figure 1). We distributed diversity treatments randomly across six fields (which were included in all analyses as random effects), with each plot separated by 3-5 m.

Neighborhood and Plot Phylogenetic Diversity
We calculated phylogenetic diversity of individual tree neighborhoods using the Picante package in R [43], as the standardized effect size of Faith's Phylogenetic Diversity [44,45]. We defined each neighborhood as the six neighbors surrounding each focal tree (see Figure 1), thus indicating the overall relatedness of the neighborhood assemblage (39). We also repeated the analysis of tree phylogenetic diversity at the plot level. We determined the phylogeny of the tree species in the BiodiversiTREE experiment (Smithsonian Environmental Research Center, Edgewater, MD, USA) by pruning the "zanne2014" angiosperm megatree [46] with Phylocom v. 4.2 [47].

Sampling Design
We used a targeted focal tree species design using 545 individuals from 15 tree species (we excluded Carya glabra due to sequencing constraints) to partition the effects of tree richness relative to species-specific differences, while controlling for variation in other factors (e.g., distance to nearby forests/fields and within-plot spatial arrangement; Figure 1b). For each species, we paired monoculture plots with a four-and 12-species plot containing the same species, with the distance to the field/forest edges matched accordingly. Six or in some cases seven individuals from each species were sampled within each of the three plots in the same within-plot spatial arrangement (Figure 1b), ensuring that only the tree neighborhood richness differed across the 1-4-12 species gradient. We repeated this structure for all tree species (N = 15 tree species × 2 plot replicates × 3 planting diversities × 6-7 trees per plot), totaling 545 focal trees sampled ( Figure 1b). Tree height was measured in 2016 to include in models as a covariate, as plant height is known to be an important factor influencing endophyte communities (e.g., [48]). To account for spatial autocorrelation in analyses, we calculated Moran's eigenvector maps (MEMs; [49]) using the adespatial R package v. 0.3-7 [50]. MEMs were calculated by extracting eigenvectors from a spatial matrix that maximizes Moran's I, an index of spatial autocorrelation structures within data. We incorporated the MEMs associated with variation in response variables into linear models as appropriate.

Sampling Protocol
For each sapling, we took 1.27-cm 2 leaf punches from three randomly selected leaves in July 2016, pooled them in a sterile microcentrifuge tube, and returned them to the lab. We immediately sterilized leaf surfaces by using sequential immersion in 95% ethanol (10 s), 10% chlorine bleach (0.525% NaOCl; 2 min), and 70% ethanol (2 min; [12,14]; adapted from [51]). We confirmed surface sterilization by aseptically rolling the leaf surfaces across a standard growth medium (2% MEA) and incubating for seven days [12,13,22,[52][53][54]. We then flash-froze the leaf samples with liquid nitrogen, pulverized them with sterile pestles, and stored the remains in DNase-free microcentrifuge tubes. We extracted DNA using a Powersoil DNA Isolation kit (formerly Mo Bio Laboratories, now owned by QIAGEN, Carlsbad, CA, USA), and stored the product in 10:1 Tris/EDTA buffer at −20 • C before library preparation.

DNA Sequencing and Data Processing
We prepared amplicon libraries by amplifying the fungal ITS2 region, using the ITS86 and ITS4R universal primer pair [55,56]. We used gel electrophoresis to determine no amplification in PCR negative controls and extraction reagent (we also sequenced negative controls). The sequencing facility at the Dalhousie University Integrated Microbiome Resource (IMR) sequenced samples on the Illumina MiSeq platform (2 × 250 bp). Unprocessed sequence data are hosted by the University of Wyoming at https://doi.org/10.15786/za14-6p32.Scripts and processed data are available at: https: //bitbucket.org/harrisonjg/griffin_biodiversitree/src/master/.
We processed sequences using USEARCH v 10.0240_i86linux32 [57]. We merged paired-end reads (minimum overlap was 10 bp) and removed primer binding regions from the ends of the merged read. We discarded reads with more than a single expected error [58]. We defined exact sequence variants (ESVs) using the UNOISE3 algorithm [59]. ESVs offer many advantages over the traditional 97% similarity threshold used to assign OTUs, including improved resolution, having a precise biological interpretation (a genotype to a locus), and facilitating information sharing among studies (see [60] for further details). We discarded ESVs that were shorter than 64 bp long, matched PhiX, or that had low complexity. We used the Warcup v. 2 [61] and UNITE v. 7.2 [62] training databases to assign taxonomic hypotheses to ESVs using the SINTAX algorithm [63], and furthermore removed host ESVs from consideration. We sequenced negative controls to identify possible contaminants, and we omitted ESVs if > 1% of their total read count was in a negative control. We rarified data to 15,000 reads per sample using the vegan R package [64].

2.7.
Testing the "Diversity Begets Diversity" Hypothesis (H1) and the "Tree Diversity-Endophyte Community" Hypothesis (H2) To estimate the effects of tree species richness on fungal richness and diversity (H1), we used linear mixed effects models with a random intercept for each unique combination of plot replicate and host taxon. We used plot richness, tree height, mycorrhizal association (ECM or AM), and the top three Moran's eigenvector maps (MEMs) as covariates within models. Overall, the model was as follows: Response = plot richness + host tree height + MEM1 + MEM2 + MEM3 + mycorrhizal type + (1|rep:host_taxon).
The reference condition for plot richness was the monoculture treatment. In addition, we calculated linear models to estimate the effects of plot richness on the response variable for each focal tree host, using the same covariates as for the mixed effects models. We did this to extract host species specific effects for plotting purposes. We repeated all models (including those below) with and without outliers; in all cases, results were qualitatively similar to those presented here.
To test for the effects of tree richness and host species on fungal community composition (H2), we conducted permutational multivariate analysis of variance (PERMANOVA) to test whether fungal community composition (Bray-Curtis dissimilarity) differed with tree richness and species. Specifically, PERMANOVA tests whether the centroids of groups of data are different in multivariate space [65]. Because PERMANOVA tests are sensitive to differences in dispersion about the centroid among sampling groups, we consequently tested for homogeneity of variance among groups using the betaDISPER function in the vegan R package [64]. When testing for the effects of tree richness, we accounted for tree species and plot replicate using the "strata" option in vegan. Conversely, when testing for the effects of tree species identity, we accounted for variation due to tree richness and plot replicate. We generated ordinations of fungal community composition using principle coordinate analysis in vegan.
For all models, we used the exponential of Shannon-Wiener diversity ( 1 D) to represent diversity as diversity equivalents [66,67]. We used Benjamini-Hochberg false discovery rate (FDR) p-value corrections to account for multiple comparisons using the multicom R package [68]. We visually assessed residuals from models to ensure minimal departures from normality. When we used Simpson's effective species number ( 2 D) as the response variable, the model residuals departed from normality; however, linear mixed effects models, because they are extensions of linear models, are quite robust to such violations [69].
We computed separate models for each focal tree that tested the effects of tree neighborhood richness and phylogenetic diversity of six surrounding neighbors on fungal richness and diversity (H1), as well as community composition (H2). Lastly, we computed separate models at the plot level to determine the effects of plot-level tree phylogenetic diversity on fungal richness and diversity.

Testing the "Tree Diversity-Endophyte Function" Hypothesis (H3)
To determine whether tree richness had differential impacts on certain guilds of fungal endophytes (H3), we assigned exact sequence variants (ESVs) to functional groups using the FunGuild database, which assigns a functional hypothesis to an ESV based on published ecological information of the individual taxon (e.g., [36,70]). We ran similar models as described above for H1 and H2, using fungal guild richness and diversity as the response variables.
For all fungal species accumulation curves, many samples appear to be somewhat close to approaching an asymptote, suggesting that our sequencing depth is relatively close to encompassing a large extent of endophyte ESV richness in communities among diversity treatments ( Figure S1). Few samples, however, were fully asymptotic, which is typical of sequence-dependent characterizations of microbial communities (e.g., [71][72][73][74] Figure 2, p < 0.01). Specifically, fungal richness decreased by 5% from one-species to four-species plots, and by 10% from one-species plots to 12-species plots.
Results shown are from a linear mixed effects model with random effects of plot replicate by host species combinations. Values in cells are beta coefficient estimates, with 95% confidence intervals in parentheses and p values denoted by asterisks. All coefficients for the effects of richness treatments are relative to monoculture plots. All coefficients for the effects of ectomycorrhizal (ECM) association are relative to those for arbuscular mycorrhizal (AM) association. Height is the height of the tallest stem of the host tree. The effect of spatial autocorrelation was accounted for using Moran's eigenvector mapping (MEM), calculated by extracting eigenvectors from a spatial matrix that maximizes Moran's I, an index of spatial autocorrelation structures within data. We converted all covariates to z scores prior to analysis. *** p < 0.01, ** p < 0.05. Conversely, from one-to four-species plots, fungal Simpson's effective species number increased by 4%; however, from one-to 12-species plots, Simpson's effective species number did not change (Table 1, Figure 2, 1-4 sp.: p < 0.05; 1-12 sp.: p > 0.10). . We used linear mixed effects models with a random intercept for each unique combination of plot replicate and host taxon. The reference condition for plot richness was the monoculture treatment. In addition, we calculated linear models to estimate the effects of plot richness on the response variable for each focal tree host, using the same covariates as for the mixed effects models. Beta coefficients from these linear models are plotted here (some species are very close to zero because there was little change in Simpson's among tree richness treatments). The dotted line is the average effect of diversity treatment across hosts. Overall, fungal richness and diversity equivalents decreased as tree richness increased (richness: p < 0.01; Shannon's: p < 0.01). Conversely, increases in one-to four-species plots increased fungal Simpson's effective species number; however, increases in one-to 12-species plots had no effect on Simpson's (1-4 sp.: p < 0.05; 1-12 sp.: p > 0.10). After correcting for multiple tests with false discovery rates (FDR) values, the effects of tree richness on fungal richness, Shannon's diversity equivalents, and Simpson's effective species number did not differ among tree species (Tables S2-S4; p > 0.05). . We used linear mixed effects models with a random intercept for each unique combination of plot replicate and host taxon. The reference condition for plot richness was the monoculture treatment. In addition, we calculated linear models to estimate the effects of plot richness on the response variable for each focal tree host, using the same covariates as for the mixed effects models. Beta coefficients from these linear models are plotted here (some species are very close to zero because there was little change in Simpson's among tree richness treatments). The dotted line is the average effect of diversity treatment across hosts. Overall, fungal richness and diversity equivalents decreased as tree richness increased (richness: p < 0.01; Shannon's: p < 0.01). Conversely, increases in one-to four-species plots increased fungal Simpson's effective species number; however, increases in one-to 12-species plots had no effect on Simpson's (1-4 sp.: p < 0.05; 1-12 sp.: p > 0.10). After correcting for multiple tests with false discovery rates (FDR) values, the effects of tree richness on fungal richness, Shannon's diversity equivalents, and Simpson's effective species number did not differ among tree species (Tables S2-S4; p > 0.05).

Shannon's Fungal Diversity
Shannon's diversity equivalents decreased as plot richness increased (Table 1, Figure 2, p < 0.01). Specifically, Shannon's diversity decreased by 6% from one-species plots to four-species plots, and by 8% from one-species plots to 12-species plots.
In our models, we found that endophyte richness and diversity were higher in ECM-associated tree species compared to AM-associated species (Table 1, Figure S2, p < 0.05). On average, ECM tree species hosted almost 30 more ESVs compared to AM species, and Shannon's diversity equivalents were on average 55% higher for ECM species compared to AM trees (Table 1, Figure S2).

Tree Neighborhood Richness and Neighborhood and Plot Phylogenetic Diversity
Overall, fungal richness and Shannon's diversity decreased as tree neighborhood phylogenetic diversity of six neighboring trees increased ( Table 2, richness: p < 0.01, diversity: p < 0.01; Figure S3). Similarly, fungal richness and Shannon's diversity decreased as tree richness of six neighboring trees increased (Table S1; richness: p < 0.01, diversity: p < 0.01). Moreover, fungal richness (but not Shannon's diversity) decreased as plot phylogenetic diversity increased ( Table 2, richness: p < 0.05; Figure S4).  Values in cells are beta coefficient estimates, with 95% confidence intervals in parentheses and p values denoted by asterisks. All coefficients for the effects of richness treatments are relative to monoculture plots. Height is the total height of the host tree. SES_PD stands for the standardized effect size of Faith's Phylogenetic Diversity. The effect of spatial autocorrelation was accounted for using Moran's eigenvector mapping (MEM), which were calculated by extracting eigenvectors from a spatial matrix that maximizes Moran's I, an index of spatial autocorrelation structures within data. We converted all covariates to z scores prior to analysis. *** p < 0.01, ** p < 0.05, * p < 0.1.

The "Tree Diversity-Endophyte Community" Hypothesis (H2): Host Tree Identity and Diversity Explained Little Variation in Fungal Community Structure
PERMANOVA models demonstrated that endophyte community composition did not vary significantly among tree species, but did explain~9% of the overall variation in endophyte communities ( Figure S5). By comparison, plot richness was a significant predictor of variation in community composition, but explained little variation overall (<1%) ( Figure S5, p = 0.04). In addition, changes in tree neighborhood and plot phylogenetic diversity only accounted for <1% of variation in fungal endophyte community composition, as determined through redundancy analysis.
Overall, increases in tree richness decreased ESV richness within all guilds as well as guild diversity for three guilds. Specifically, ESV richness within all guilds (not including unidentified ESVs) decreased by 8-11% from one-to 12-species plots (Table 3, Figure 3). Moreover, overall guild diversity decreased by 5% from one-to 12-species plots (t = −2.16, p < 0.05). Pathogen, saprotroph, and parasite diversity decreased by 5-6% from one-to 12-species plots (Table 4, Figure S6; pathogens: p < 0.05; saprotrophs and parasites: p < 0.01).  Results shown are from linear mixed effects models with random effects of plot replicate by host species combinations. Values in cells are beta coefficient estimates, with 95% confidence intervals in parentheses and p values denoted by asterisks. All coefficients for the effects of richness treatments are relative to monoculture plots. Height is the total height of the host tree. The effect of spatial autocorrelation was accounted for using Moran's eigenvector mapping (MEM; see main text). We converted all covariates to z scores prior to analysis. *** p < 0.01, ** p < 0.05, * p < 0.1.  . We used linear mixed effects models with a random intercept for each unique combination of plot replicate and host taxon. The reference condition for plot richness was the monoculture treatment. Increases in tree richness decreased ESV richness within all guilds as well as guild diversity. *** p < 0.001, ** p < 0.01, * p < 0.05.
We identified 42% (1,179) of our ESVs as saprotrophs, 40% (1123) as pathogens, 18% (505) as symbionts, 7% (197) as parasites, and 6% (168) as lichens. We classified 555 ESVs-20% of all ESVsinto more than one guild (note: this is why the percentages do not add up to 100%; see supplementary spreadsheet "Griffin et al. FunGuild results"). The confidence rank assigned to functional guild placements of 2003 ESVs (71%) was "probable," 493 ESVs were "possible," and 310 were "highly probable" (see supplementary spreadsheet "Griffin et al. FunGuild results"). The guilds of 1056 ESVs (27.3%) were unidentifiable.  Figure 3. The effects of tree richness on fungal ESV richness of the most common functional guilds among 15 focal tree species (N = 545 focal trees). We used linear mixed effects models with a random intercept for each unique combination of plot replicate and host taxon. The reference condition for plot richness was the monoculture treatment. Increases in tree richness decreased ESV richness within all guilds as well as guild diversity. *** p < 0.001, ** p < 0.01, * p < 0.05.  Results shown are from linear mixed effects models with random effects of plot replicate by host species combinations. Values in cells are beta coefficient estimates, with 95% confidence intervals in parentheses and p values denoted by asterisks. All coefficients for the effects of richness treatments are relative to monoculture plots. Height is the total height of the host tree. The effect of spatial autocorrelation was accounted for using Moran's eigenvector mapping (MEM; see main text). We converted all covariates to z scores prior to analysis. We calculated the Shannon's diversity equivalents ( 1 D) of all guilds using the relative abundances of the five guilds considered. *** p < 0.01, ** p < 0.05, * p < 0.1.

Discussion
Counter to our hypotheses, we found strong evidence that increasing tree richness and phylogenetic diversity decreased fungal species and functional guild richness and diversity within the first three years of a tree diversity experiment (H1, H2, and H3). These findings differ from studies demonstrating that increases in plant diversity typically increase the diversity of associated trophic levels, including foliar bacterial communities [8,15]. Not only did we expect a positive relationship between tree richness and fungal richness and diversity, we expected that changes in tree neighborhood richness and phylogenetic diversity would strongly structure endophyte community composition (H2); however, this did not occur. In addition, increasing tree richness caused a decrease in guild richness for all five functional guilds and a decrease in diversity for three guilds (the "Tree Diversity-Endophyte Function" hypothesis), suggesting that endophytes may mediate the relationship between tree richness and metrics of ecosystem function (i.e., richness and diversity of pathogens, parasites, and saprotrophs). Overall, our results are somewhat surprising, particularly given that many empirical studies in plant communities fail to detect the effects of plant diversity on ecosystem services in early stages after establishment (e.g., [75][76][77], reviewed by [78]). We discuss our results in detail below.

Counter to Traditional BEF Theory, Tree Richness and Phylogenetic Diversity Decreased Endophyte Diversity
We expected that increases in tree richness would cause an increase in fungal endophyte diversity, but we in fact found the opposite pattern. Indeed, numerous studies have demonstrated that increases in plant richness increase the diversity of other trophic levels, including herbivores and epiphytic bacteria [15,[79][80][81][82]. These patterns are thought to be driven by host specialization, whereby increases in plant richness attract more diverse associate communities (e.g., the "Resource Specialization Hypothesis"; [83,84]). Endophytes, however, especially for smaller seedlings and saplings, may have stronger dispersal limitation compared to plant-associated microbes on larger adult hosts [48,85,86] [48]) demonstrated that increases in grassland biomass among four sites in North America consistently decreased fungal endophyte diversity. They suggested that increases in plant biomass led to higher dispersal rates of competitively dominant fungi, thus decreasing overall fungal diversity. Similarly, increases in tree richness in our experiment may have decreased fungal diversity by increasing the dispersal potential of competitive generalist fungi.
Alternatively, increases in tree richness may result in fewer encounters with suitable hosts for specialist fungi. This idea is consistent with island biogeography theory, whereby more diverse stands represent archipelagos of smaller islands of suitable hosts for specialist endophytes (sensu [87]).
In fact, many studies have demonstrated that pathogenic fungi are less likely to disperse to distant neighbors, particularly in less dense stands (e.g., [88][89][90][91]). In our experiment, monoculture plots contained 255 conspecific trees, whereas 12-species plots only had 21 conspecific trees. Therefore, in monoculture plots, the distance between suitable islands (i.e., among tree hosts) is lower compared to in 12-species plots. For this reason, increases in tree richness may have resulted in less rich and less diverse endophyte communities.
Conversely, fungal endophytes may not be specialized enough for their richness and diversity to increase as tree richness increases. The "Resource Specialization Hypothesis" is predicated on the assumption that faunal communities exhibit a high degree of host specificity (sensu [83]). However, studies have consistently demonstrated that increases in tree richness cause increases in herbivore and bacterial epiphytic communities [15,79], which are typically highly specialized on plant hosts (e.g., [25,26,[92][93][94]). Specifically, over 90% of herbivores exhibit some degree of host preference [84,93,94]), and host identity explains as much as half of bacterial and fungal epiphytic community composition [24][25][26]. Comparatively, our findings, taken with other studies of fungal and bacterial endophytes, demonstrate that host species accounts for as little as 5-9% of the variation of endophyte community composition [14,[95][96][97]. If endophytes were more specialized, monoculture plots would be filled with dominant specialists, as predicted by our dispersal argument. This would lead to an increase in endophyte diversity with tree diversity because tree hosts would primarily host more evenly distributed generalist endophyte taxa.
It is possible that the effects of plant diversity on fungal communities may be caused by indirect trophic cascades rather than dispersal limitations of specialists. For instance, tree diversity may structure below-ground microbiomes or induce changes to leaf chemistry, which can mediate endophyte communities. For example, it has been well documented that soil microbes can trigger systemic chemical responses in plants that increase above-ground growth and host resistance to insects and other enemies ( [98][99][100]; reviewed by [101,102]). In one study, Badri et al. [98] demonstrated that soil microbial communities collected from 11 different sites globally and inoculated into Arabidopsis thaliana roots differentially regulated leaf chemistry (metabolomic) profiles and insect feeding behavior. Specifically, soil microbiome inoculations altered the production of phenolics and sugars, which are important drivers of foliar microbial communities (e.g., [103][104][105][106]). Thus, variations in soil microbial diversity may affect leaf tissue quality (e.g., phenolic content, biomass accumulation, defensive metabolites; [107,108]), which may then directly regulate leaf microbiomes. Ultimately, we suggest that future studies should consider both below-ground and above-ground processes, as well as the leaf metabolome, to fully understand the community-level effects of plant diversity.
We were surprised that the effects of tree richness caused differences in Shannon's diversity equivalents and Simpson's effective species number. Overall, Shannon's equally emphasizes the richness and evenness components of diversity, while Simpson's emphasizes the evenness component, which causes Shannon's to put more weight on rare species and Simpson's to emphasize dominant species [66,67,109,110]. Thus, our results suggest that the most abundant fungal taxa are not responding to neighborhood diversity as much as the rarer fungal taxa. Indeed, the top 20 most abundant ESVs in our study remained relatively constant despite changes in tree richness, while all other ESVs (not top 20) decreased by 3-6% on average as tree richness increased. This suggests that differences in tree communities disproportionately impact rare common species, which suggests that a "core microbiome" may exist among coexisting tree species, independent of tree neighborhoods.
Because plant phylogeny was confounded with the mycorrhizal association (AM vs. ECM) covariate in our analyses, mycorrhizal association in addition to plant taxonomy may be important in structuring above-ground endophyte communities. Indeed, we found that endophyte richness and diversity were higher in ECM-associated tree species compared to AM-associated species ( Figure S2). Specifically, the type of mycorrhizal association (e.g., ectomycorrhizal versus arbuscular) may directly or indirectly influence foliar fungal endophytes, particularly by differentially influencing nitrogen and phosphorus cycling ( [111,112]; reviewed by [113]). As foliar nitrogen and phosphorus concentrations structure microbial communities (e.g., [25,114]), it is possible that host mycorrhizal association type (i.e., ECM vs. AM) may be more important in shaping leaf microbiomes than either host taxonomy or diversity (at plot or neighborhood scales). Therefore, future studies are needed to explicitly test whether mycorrhizal association, rather than host phylogeny, is important in structuring foliar endophytes.
We recognize that the surrounding forest matrix may be a source of fungal inoculum, though we do not believe that there is a disproportionate number of AM vs. ECM mycorrhizal propagules coming from the forest surrounding the site. Indeed, fungal endophytes are largely horizontally transmitted from the surrounding environment [11,16,20] and may even be dispersed via leaf litter [23,115]. Extensive surveying in the surrounding forest (a CTFS-ForestGEO plot is located within three miles of the site; [116,117]) has shown that local forests are comprised of tulip poplar (Liriodendron tulipifera; AM species) in the overstory, as well as several oaks (Quercus spp.; ECM), beech (Fagus grandifolia; ECM) and hickories (Carya spp.; ECM). Red maple (Acer rubrum; AM) and blackgum (Nyssa sylvatica; AM) make up mid-canopies [117]. Thus, because there is a mix of both AM-and ECM-associated trees in the surrounding forest, we are confident that there is no bias in potential AM or ECM propagules from surrounding areas reflected in our results.

Host Tree Richness Mediated the Function of Fungal Endophytes
In addition to fungal richness and diversity, we also found that increases in tree richness decreased fungal functional guild diversity. Specifically, the richness of every functional guild decreased with increasing tree richness, and the diversity of pathogens and saprotrophs decreased with increasing tree richness. We acknowledge that FunGuild results are tentative, particularly given that the database is biased towards soil fungi and pathogens and the impacts of particular fungal taxa may differ among different plants hosts and environmental conditions [36]. Additionally, such classifications may be biased towards pathogens because of a disproportionately high interest in pathogens. Therefore, it is possible that even a larger proportion of foliar endophytes function as saprotrophs than we report here. Overall, however, our results are consistent among all of the metrics tested and the message is clear: increases in tree richness and phylogenetic diversity decrease species and functional richness and diversity of fungal endophytes. We identified almost half (~40%) of the ESVs as putative pathogens. Latz et al. [118] found that greater plant richness indirectly decreased the incidence of soil pathogens through changes in soil bacteria that were antagonistic to pathogens. Thus, it is possible that tree communities cultivate or recruit beneficial bacterial communities to competitively exclude fungal endophytes or produce antimicrobial compounds to ward off fungi (reviewed by [20,33]).

Conclusions
Here, we highlight the importance of plant diversity and its interaction with the plant microbiome and ultimately explore how these interactions may impact plant performance early in forest stand development. Numerous studies have demonstrated that the tree microbiome among forest seedlings establishes early and has significant impacts on plant hosts, potentially determining relative competitive abilities and ultimately shaping plant phenotype and community interactions (e.g., [12][13][14]22,23]). We expand upon these findings to demonstrate that tree community diversity (species and phylogenetic) is critical in driving endophyte community and functional diversity. Specifically, we found that tree neighborhood richness and phylogenetic diversity regulated fungal endophyte assemblages within the first three years of establishment. Future studies should test whether the magnitude of the effects of plant diversity strengthens over time, a trend that is common among plant diversity experiments [76,119,120].
Supplementary Materials: The following are available online at http://www.mdpi.com/1424-2818/11/12/234/s1, Table S1: Effects of tree richness of six neighboring trees and covariates on fungal endophyte diversity and richness; Table S2: The effects of tree (plot) richness on fungal endophyte richness among 15 tree species (N = 545 focal trees); Table S3: The effects of tree taxonomic (plot) richness on fungal endophyte Shannon's diversity equivalents among 15 tree species (N = 545 focal trees); Table S4: The effects of tree taxonomic (plot) richness on fungal endophyte Simpson's effective species number among 15 tree species (N = 545 focal trees); Figure S1: Rarefaction curves for A. all samples, B. 1 sp., C. 4 spp., and D. 12 spp. plot samples compiled together among all host species; Figure S2: The effects of mycorrhizal association (AM vs. ECM) on A. fungal endophyte richness ( 0 D), B. Shannon's diversity equivalents ( 1 D), and C. Simpson's effective species number ( 2 D) among 15 focal tree species (N = 545 focal trees); Figure S3: Partial regression plots from a mixed effects model examining the effects of tree neighborhood phylogenetic diversity on fungal endophyte A. richness and B. Shannon's diversity equivalents; Figure S4: Partial regression plots from a mixed effects model examining the effects of plot-level tree phylogenetic diversity on fungal endophyte A. richness and B. Shannon's diversity equivalents; Figure S5: The effects of A. tree richness and B. 15 focal tree species on fungal endophyte community composition (N = 545 focal trees); Figure S6: The effects of tree richness on fungal ESV diversity (Shannon's diversity equivalents) of the most common functional guilds among 15 focal tree species (N = 545 focal trees). Griffin, et al. FunGuild results. Spreadsheet detailing the functionality of ESVs in study.