Grass species identity shapes communities of root and leaf fungi more than elevation
Study sitesWe sampled foliar fungal endophytes and root fungi (root endophytes and AM fungi) in the Colorado Rockies at the Rocky Mountain Biological Laboratory, Gunnison Co., Colorado, USA (38°57’N, 106°59’W). This region has predictable decreases in air temperature (c. 0.8 °C per 100 m; [40]) and declines in soil nutrients with altitude [41], but increases in precipitation, mainly as snow [42]. The entire region is warming at rates of 0.5–1.0 °C per decade [43].To capture environmental, spatial, and grass-host specific variation in fungal guilds, we sampled 66 sites encompassing 9–13 elevations from each of six altitudinal gradients in July 2014 (Supplementary Table S1, Supplementary Fig. S1). Elevational gradients represented separate mountains in the Gunnison Basin and were located within 20 km of each other. We created a regional climate model to interpolate average number of growing degree days (GDD, base 0 °C), mean annual temperature (MAT), maximum temperature (Tmax), minimum temperature (Tmin), mean annual precipitation (MAP), and mean snow depth (MSD) for each site based on data from 29 local meteorological stations [44]. At each site, soil edaphic parameters were measured on dried soil at the UC Davis soils lab (see [24] for more details) and soil nutrients at Western Ag (Saskatoon, Canada). Soil pH was measured in a 1:1 solution with diH2O, and soil moisture was measured gravimetrically. Physical characteristics of each site (e.g., aspect, soil depth, elevation) were measured as described in Lynn et al. [44]. Environmental variation across sites was large. For example, MAT varied from 7.1 to 13.3 °C, MAP from 563 to 1171 mm, and Total N from 2 to 316 ug/g dry soil (Table S1).Host plant speciesWe focused on grasses because grasslands cover ~20% of Earth’s land surface [45] and dominate subalpine meadows of the Rocky Mountains. In addition, individual grass species spanned the entire elevational range of our study system [46], whereas tree, shrub, and forb species did not. At each location, we sampled nine adult individuals from up to 13 grass species representing five genera (Poaceae, subfamily Pooideae; Supplementary Table S1). Many sites had fewer than 13 grass species present, but all sites, except for two, had at least two grass species. Samples were composited by tissue type (leaves v. roots) and grass species within each site.Fungal compositionCollected root and leaf samples were surface sterilized (1 min in 95% ethanol, 2 min in 1% sodium hypochlorite solution, and 2 min in 70% ethanol) over ice to focus on the endophytic fungal community [34]. Following surface sterilization, samples were rinsed in purified water (Milli-Q Integral Water Purification System, EMD Millipore Corporation, Billerica, MA), stored in RNAlater, and refrigerated. All samples were then frozen in liquid nitrogen and ground using a mortar and pestle. Total DNA was extracted from ~50 mg of ground sample using QIAGEN DNeasy plant extraction kits (QIAGEN Inc., Valencia, CA).Fungal composition was characterized using barcoded primers targeting the ITS2 region for leaf and root endophytes [47], and FLR3-FLR4 primers targeting ~300 bp in the 28S region for AMF [48]. Each PCR contained 5 μL of ~1–10 ng/μL DNA template, 21.5 μL of Platinum PCR SuperMix (Thermo Fisher Scientific Inc., Waltham, MA), 1.25 μL of each primer (10 μM), 1.25 μL of 20 mg/mL BSA, and 0.44 μL of 25 mM MgCl2. For the ITS2 primers, the reactions included an initial denaturing step at 96 °C for 2 min, followed by 24 cycles of 94 °C for 30 sec, 51 °C for 40 s, and 72 °C for 2 min, with a final extension at 72 °C for 10 min. For the 28S primers, reactions started with an initial denaturing step at 93 °C for 5 min, followed by 33 cycles of 93 °C for 1 min, 55 °C for 1 min, and 72 °C for 1 min, with a final extension at 72 °C for 10 min.Three PCR replicates from each sample were pooled and then cleaned and concentrated using a ZR-96 DNA Clean & Concentrator-5 (Zymo Research Corporation, Irvine, CA). PCR was then carried out on all samples to add dual indexes and Illumina sequencing adaptors; each reaction began with an initial denaturing step at 98 °C for 30 s, followed by 7 cycles of 98 °C for 30 s, 62 °C for 30 s, and 72 °C for 30 s, with a final extension at 72 °C for 5 min. Sequencing was performed by the Genomic Sequencing and Analysis Facility at The University of Texas at Austin using paired-end 250 base Illumina MiSeq v.3 chemistry (Illumina, Inc., San Diego, CA). We aimed to obtain a minimum of 30,000 reads/sample for the ITS2 region and 20,000 reads/sample for the 28S region. All sequences are deposited in the NCBI SRA database under accession number (PRJNA639093).BioinformaticsWe processed reads to generate OTUs using commands from USEARCH (v9.2.64). Reads from previous studies [24] and this study were clustered together to improve OTU delineations for a total of 36,754,931 reads. We merged paired-end reads using the fastq_mergepairs from USEARCH with “fastq_maxdiffs” set to 20 and “fastq_maxdiffpct” set to 10 to ensure proper merging at a low error rate. The merged reads and the forward unmerged reads were trimmed at the primer sites using cutadapt with “e” set to 0.2, “m” set to 200, and untrimmed reads were discarded. Merged reads were filtered using fastq_filter from USEARCH with “fastq_maxee” set to 1.0. The forward reads were first trimmed to 230 using fastx_truncate from USEARCH with “trunclen” set to 230 and then filtered by fastq_filter from USEARCH with “fastq_maxee” set to 1.0. We then concatenated the merged and forward reads into one file and de-replicated using fastx_uniques from USEARCH with “minuniquesize” set to 2. After these steps, 11,357,274 sequences remained. We clustered these sequences to form OTUs at 97% similarity [49] using cluster_otus command from UPARSE. The reads (all reads before filtering step) of each sample were mapped to OTUs with usearch_global from USEARCH with “id” set to 0.97. We determined taxonomy for the representative OTUs using sintax from USEARCH with the database set to UNITE all eukaryotes (v. 8.2) “strand” set to both and “sintax_cutoff” set to 0.8 [50]. Representative OTUs were also blasted against Genbank with “perc_identity” set to 80 and “max_target_seqs” set to 50. All OTUs identified as “fungi” were retained, and OTUs labeled as “unknown” or “unidentified” were manually inspected based on blast results to determine retention. Our filtering criteria left between 5 and 418 OTUs per sample (Supplementary Table S2).Due to low fungal abundance in leaves [34], many leaf samples were dominated by plant sequences (average ~78% plant reads). Therefore, fungal sequence numbers in leaf samples were low, despite adequate sequencing depth to capture trends in fungal endophyte communities across sites based on prior analyses [24, 34, 35]. We included only samples that contained at least 50 fungal sequences after data processing (Leaves N = 192, Roots N = 191, AMF N = 251), and most samples had much greater sequencing depth, especially for roots (Supplementary Table S2). Nevertheless, there were no correlations between sequence read depth and richness, alpha diversity, or evenness of our samples (P > 0.05 in all cases), and plant species did not differ in the average sequencing depth for samples (P > 0.05). Data for each fungal OTU were transformed to the proportion of total sequence abundance to minimize any differences in sampling effort [51].Diversity and compositionWe calculated the alpha diversity metrics of richness, Shannon’s Diversity, Inverse Simpson’s Diversity, and Pielou’s Evenness. For each fungal guild, differences among plant species and elevation in alpha diversity were first determined using a general linear mixed effects model with plant species (categorical) and elevation (continuous) as fixed effects and site nested within elevation gradient (e.g., mountain identity, Supplementary Table S1, Supplementary Fig. S1) as random effects to account for the lack of statistical independence among plant species sampled at the same site and among sites located within the same mountain elevation gradient (Supplementary Fig. S1). Models were constructed using the lmer function in R package lme4 [52, 53]. To address, do fungal community patterns along environmental gradients differ among guilds: leaf endophytes, root endophytes, or arbuscular mycorrhizal fungi?, we then compared alpha diversity metrics among fungal guilds using a general linear mixed effects model with fungal guild, plant species, and elevation as fixed effects and site nested within elevation gradient as random effects. In all models, we evaluated parameter fit with analysis of deviance using Wald chi-square tests and corrected for multiple comparisons using a false discovery alpha of 0.05. Differences among grass species were determined using Tukey post-hoc tests.Because elevation is a good proxy for variation in both climate and soil parameters (Supplementary Table S1), in all community analyses, we first ran models with grass species and elevation to parse biotic versus abiotic influences on fungal OTUs, then secondly ran full variance partitioning models with all environmental covariates (Supplementary Table S1, climate, physical, soil) in addition to grass species identity and space (gradient location, Supplementary Fig. S1). Because leaf and root endophytes were sequenced using different primers than AM fungi, we could not compare composition among the three guilds directly. Instead, we compared the relative influence of biotic and abiotic drivers on fungal composition within each guild to compare patterns among guilds. To do so, we first used distance-based redundancy analysis (dbRDA) to analyze the effects of plant host species and elevation on fungal composition for general fungal communities in leaves and roots and separately for AM fungal communities in roots. All models were run on quantitative Jaccard indices of fungal composition for each guild and included site nested within elevation gradient (e.g., mountain side, Supplementary Fig. S1) as random effects. Second, to evaluate which environmental variables most strongly influenced fungal composition, we further partitioned variance in fungal composition due to grass species, climate variables (MAP, MAT, MSD, Tmax, Tmin, and GDD), soil variables (total nitrogen, total phosphorus, nitrate, ammonium, calcium, magnesium, potassium, iron, manganese, sulfur, aluminum, soil pH, soil gravimetric moisture content), physical variables (aspect degree, aspect category (e.g., cardinal direction), slope, soil depth, and elevation) and spatial variables (latitude and longitude) using the varpart function in Vegan v. 2–5.3 [54]. Plots of fungal composition by plant host were also generated using dbRDA separately for each fungal guild. Spatial variables were de-trended and tested for spatial autocorrelation using the ade4 package v. 1.7–16 [55]. When we detected significant spatial autocorrelation eigenvectors, we included these in the spatial variable matrix. To characterize how many fungal taxa occurred in multiple plant taxa and elevations, we used the VennDiagram package v. 1.6.20 [56].Turnover and rewiringTo evaluate whether fungal composition was driven by grasses associating with different fungal taxa or differing relative abundances of the same fungal taxa, we first performed a beta partitioning analysis using betapart v. 1.5.3 [57]. Each fungal guild was analyzed separately. Next, to examine turnover in the abundances of fungal functional groups (pathogens, saprotrophs, mutualists), we defined groups using the FungalTrait database, which merges previous databases into one cohesive framework of 17 functional trait types (referred to here as functional groups; [58]). We recognize that fungal functions are highly environmentally dependent and therefore these functional groups may represent potential function more than actual function. Functional group identity was ascribed to 60% of leaf endophyte and 62% of root endophyte fungal taxa. Then, cumulative abundance of proportionally transformed sequence reads in each functional group was analyzed using a general linear mixed effects model with grass species and elevation as fixed effects and site nested within elevation gradient as random effects, as above. Finally, we defined indicator species within the OTUs that comprised at least 1% of the total abundance of each fungal guild by grass host, gradient, and elevation classes (rounded to the nearest 100 m) using the indicspecies package v. 1.7.9 [59]. Functional group assignments using the FungalTrait database from above were assigned to each indicator taxon [58]. A large percentage of significant indicator taxa out of the total number of OTUs would confirm that turnover in the species identity of fungal associations is stronger than turnover in the relative abundances of the same fungal taxa.Network propertiesTo address does grass-fungal network structure track elevation?, we analyzed four properties that encompass different facets of ecological networks at the site level. First, we calculated network nestedness, or the propensity for specialists to interact with the same plant species as generalists, using the weighted NODF (Nestedness metric based on Overlap and Decreasing Fill; [60]). Second, we calculated complexity as linkage density or the average number of interactions per plant species [61]. Third, to characterize specialization, we used the H2’ Index [62]. Finally, network evenness was calculated as Alatalo’s interaction evenness [63]. In all cases, these network metrics were weighted indices to increase accuracy [64], and calculations were performed in the Bipartite package v. 2.15 [65]. To address, how much do fungal guilds differ in altitudinal variation in network structure?, we compared network-level statistics among fungal guilds using a general linear mixed effects model with fungal guild as a fixed effect, number of grass hosts as a fixed effect, and gradient as a random effect (function lmer in lme4 [52],). We compared relationships with elevation separately for each fungal guild, using general linear mixed effects models with elevation as a fixed, continuous effect, number of grass hosts within the network as a fixed, continuous effect, and gradient identity as a random effect (Supplementary Table S1, Supplementary Fig. S2). We evaluated parameter fit with analysis of deviance using Wald chi-square tests using the car package 3.0–10 in R [66].All data met model assumptions of normality of residuals and homogeneity of variance. All analyses were performed in R v. 3.5.0 [53]. More
