Ectomycorrhizal fungal diversity predicted to substantially decline due to climate changes in North American Pinaceae forests

Ectomycorrhizal fungi (ECMF) are partners in a globally distributed tree symbiosis implicated in most major ecosystem functions. However, resilience of ECMF to future climates is uncertain. We forecast these changes over the extent of North American Pinaceae forests.


| INTRODUC TI ON
Forecasting changes in the diversity and composition of microbial communities under anticipated future climates is valuable for concentrating conservation efforts (van der Linde et al., 2018) and predicting changes to ecosystem function (Bissett, Brown, Siciliano, & Thrall, 2013, Duffy et al., 2017, Koide, Fernandez, & Malcolm, 2014.
Loss of host species results in decreased ecosystem productivity and stability across a broad range of taxa (Duffy et al., 2017), including effects on microbes (Duffy et al., 2017, Laforest-Lapointe et al., 2017. Recent advances in molecular methods have made it possible to characterize current continental-scale diversity patterns of soil microbes (van der Linde et al., 2018;Talbot et al., 2014;Tedersoo et al., 2014;. However, continental-scale forecasts under future climates are unavailable for most microbial guilds, making it difficult to predict the consequences of climate change to global biodiversity and ecosystem services. Here we predict how the species richness, relative abundance of hyphal exploration strategies and composition of ectomycorrhizal fungi (ECMF) in North American pine forests will change over the next 50 years.
Ectomycorrhizal fungi are plant symbionts that dominate global temperate and boreal forest soil communities and are implicated in most major ecosystem processes (Phillips, Brzostek, & Midgley, 2013). Some ECMF taxa produce extracellular proteolytic and oxidative enzymes that liberate N from organic complexes (Lindahl & Tunlid, 2015;Talbot, Allison, & Treseder, 2008). For example, relative to ECMF species with short-distance hyphal exploration strategies, long-distance foraging ECMF species are associated with higher activities of organic N mineralizing enzymes (Hobbie & Agerer, 2010;Tedersoo, Naadel, et al., 2012). When ECMF mineralize organic N and transfer it to their host trees under elevated CO 2 , they have the potential to fuel increased photosynthetic rates (Terrer, Vicca, Hungate, Phillips, & Prentice, 2016). Additionally, by removing N from organic complexes, ECMF are hypothesized to competitively inhibit free-living soil microbes that require N to decompose and respire soil organic carbon (Averill, Turner, & Finzi, 2014). However, because not all ECMF species can mineralize organic N (Pellitier & Zak, 2018), diversity losses and shifts in composition could alter ECMF-associated functions.
However, most studies to date that have examined ECMF or whole fungal community responses to simulated climate change have found fairly small effects (Fernandez et al., 2017;Mucha et al., 2018;Parrent, Morris, & Vilgalys, 2006;Tu et al., 2015) relative to natural changes in fungal communities observed along large natural gradients of temperature and precipitation (Jarvis, Woodward, Alexander, & Taylor, 2013;Nottingham et al., 2018;Peay et al., 2017;Talbot et al., 2014;Tedersoo et al., 2014). Yet, few datasets currently exist with spatial resolution necessary to make accurate predictions of ECMF response to climate change across relevant geographic regions (Mohan et al., 2014).
We used next-generation DNA sequencing to determine the ECMF species composition in sites spread across North America (Omernik & Griffith, 2014) (Figure 1, Table S1). To isolate the effect of climate on ECMF communities while minimizing the known effects of vegetation-type variation on microbial community structure and function (Fierer et al., 2012, Tedersoo, Naadel, et al., 2012, we placed all of our sites in forests dominated by single species of trees, all from the family Pinaceae (an obligate ECMF host lineage). The Pinaceae are ideal for exploring environment-community--function relationships across Kingdom Fungi because they have a broad distribution across North America and show low levels of host specificity for mycorrhizal fungi within the family (Ishida, Nara, & Hogetsu, 2007;Rusca, Kennedy, & Bruns, 2006). For example, North American pines readily associate with European pine-associated ECMF (Vellinga, Wolfe, & Pringle, 2009) and co-occurring Pinaceae and angiosperms often share the most common ECMF (Kennedy, Izzo, & Bruns, 2003).
Thus, while our sampling is restricted to pine-dominated systems, the results should be applicable to other forest types.
We fit nonlinear models to the relationships between climate and ECMF species richness, relative abundance, and community dissimilarity. Additionally, to identify ECMF community thresholds to changes along temperature gradients, we performed threshold indicator taxa analysis (Baker & King, 2010). In order to project potential changes to ECMF communities with climate change, we used the climate-envelope approach, which uses statistical models fit along spatial-environmental gradients to predict ECMF responses along temporal-environmental gradients.

| Sampling
We measured ECMF community composition in 68 sites spread across North American forests dominated by a single plant family, the Pinaceae (Table S1). Sampling was carried out in 2011 and 2012 near the period of peak plant biomass production for a given region.
In each plot, 13 soil cores were collected from a 40 m × 40 m grid ( Figure 1d, Figure S1). To look at spatial turnover of community and function at the local to landscape scales, we ensured that each plot had at least one other plot located within a 1-50 km range ( Figure S1b).
At each point in the plot, fresh litter was removed and a 14-cm deep, 7.6-cm diameter soil core was taken and immediately separated into a humic (O) horizon and mineral (M) horizon. This resulted in a total of 26 soil samples collected per plot (13 sample points × 2 horizons).
After removal, soils were kept on ice until processed. Soils were sieved through a 2-mm mesh to remove roots and rocks and homogenized by hand. A ~0.15 to 0.25 g subsample was placed directly into a bead tube from the Powersoil DNA Extraction Kit (MoBio), and the samples were stored at 4°C until DNA extraction. Before extraction, samples were homogenized for 30 s at 75% power using a Mini-Beadbeater (BioSpec). A second subset of soil core x horizon samples was stored in −80°C freezer (within 48 hr of collection) for soil chemical analysis.

| Soil chemistry
Frozen soils were thawed and analysed for pH in a 1:1 water ratio using a glass electrode. Total extractable ammonium and nitrate concentrations were analysed in 2.0 M potassium chloride extracts of each soil sample using a WestCo SmartChem 200 discrete analyzer at Stanford University. For site-level values, we took the average of all soil cores processed for each site. Soil chemical variables were included in statistical models but later dropped during model selection (see Statistical Analysis and Supplemental Materials).

| Molecular methods
To characterize fungal communities, we sequenced the internal transcribed spacer (ITS) region of the nuclear ribosomal RNA genes, the official barcode of life for fungi (Schoch et al., 2012).
Because of improvements in technology during the course of this project, soil samples were sequenced using two different platforms. Soil samples from 25 sites were sequenced via 454 pyrosequencing as per Talbot et al. (2014). The remaining 43 sites were sequenced on an Illumina MiSeq at the Stanford Functional Genomics Facility using the primer constructs and protocols from Smith and Peay (2014). Using a common set of soil samples from this study sequenced on both platforms we have previously demonstrated that both richness and species composition are highly reproducible and strongly correlated between the two platforms ) so that combining samples should not cause any bias in our analyses.

| Bioinformatics
Samples sequenced on the 454 platform were cleaned and denoized in QIIME (Caporaso et al., 2010;Reeder & Knight, 2010), after which we extracted the ITS1 region (Nilsson et al., 2010). For samples sequenced on the Illumina platform, samples were first trimmed using Cutadapt (Martin, 2011) and Trimmomatic (Bolger, Lohse, & Usadel, 2014), and then merged using USEARCH (Edgar, 2010). At this point cleaned 454 and Illumina sequences were merged into a single FASTA file, where sequences were dereplicated, and clustered into specieslevel operational taxonomic units (OTUs) at 97% sequence similarity using USEARCH. We removed all singletons and chimeras, and dropped occurrences <0.025% of relative sequence abundance within a sample to account for tag-swapping (Carlsen et al., 2012). We first used the BLAST tool with the UNITE reference database (Koljalg et F I G U R E 1 (a) Scatter-plots of all pairwise combinations of climatic variables from the = 68 sites spread across North America. Shaded polygons define the range of climates where we project our models of ectomycorrhizal species richness and composition. Points are coloured in red, green and blue according to scaled latitude and longitude, as shown in the map (b), which has a shaded region that corresponds to the geographic distribution of Pine trees from our plot network and the distribution of climate variables from (a). The ectomycorrhizal fungal community of each site was characterized from 13 soil cores in a nested sampling design in each site (c), with samples from the (O)rganic and (M)ineral horizon from each core (d) al., 2013) to eliminate potentially non-fungal taxa; and then used the naïve Bayesian classifier from the Ribosomal Database Project (Wang, Garrity, Tiedje, & Cole, 2007) along with the Warcup ITS reference set (Deshpande et al., 2016) to assign taxonomy. OTUs with sufficiently confident taxonomic assignments were then matched to functional guilds using FUNGuild (Nguyen et al., 2016). The relative abundance of ECMF exploration strategies (% of ECMF OTUs) per site were assigned by matching ECMF genera against published lists (Agerer, 2001(Agerer, , 2006Tedersoo & Smith, 2013), which are available online (www.deemy.de).
Strategies were assigned to the following categories according to Fernandez et al. (2017): contact short (CS), contact medium (CM) and medium long (ML). Because of potential differences in the relationship between DNA copy number and fungal biomass, particularly among functional groups known to differ in their morphology, we restricted our analyses to intraguild comparisons (e.g. differences in relative abundance of ML among sites). Full description of the bioinformatic methods is available in the online supplement.

| Statistical analyses
To correct for variability in DNA sequencing depth between samples from the two sequencing platforms, we first rarefied our sequences We fit generalized additive models (GAMs) of the ECMF species richness and relative abundance of CS, CM and ML exploration strategies of each site as a function of the bioclimatic variables using the 'mgcv' package in R. We choose to use nonparametric GAMs, rather than linear models, so that our statistical models would have sufficient flexibility to capture curvilinear responses of ECMF diversity and abundance to environmental factors. We set the maximum number of knots in all GAM models to 4, which constrains against over-fitting but allows the fitted splines sufficient flexibility to capture potential saturating, sigmoidal, uni-and bi-modal responses of ECMF communities to environmental gradients. Additionally, we performed k-fold cross validation of model results to measure the extent that model performance is stable to data inclusion/exclusion. We divided the data in k = 10 random partitions, each with a different partition of ~10% of data withheld and manually calculated R 2 values for each of the models when compared with the full n = 68 sites using the following equation: R 2 = 1 − sum (actual − predicted response) 2 /sum(actual − mean response) 2 . For models that are 'over-fit', the range of R 2 values should be large, whereas for relatively stable models, the range should be narrow and not substantially less than the R 2 of models fit with data when all k = 10 partitions were present.
In addition to the four climatic variables for which we have recent historical and projected 2070 data, we also considered GAM models with the following covariates: (a) only in situ soil pH and NH4-N; (b) only soil predictors from global rasters of surface horizon chemistry [soil pH in KCl (Hengl et al., 2017) and total N density (Task, 2000)], (c) only estimates of total atmospheric N-deposition (Hember, 2018), and (d) combinations of climate with soil in situ, soil raster, and N-deposition predictors. We selected our climate-only model based on its superior performance relative to model-complexity according to the generalized cross validation statistic (Tables S2-S7, Figures S2-S6).
Model extrapolations are valid only within the range of conditions used to fit those models. To account for this, we projected our models exclusively within the geographic distributions of the 12 most abundant pine tree species from our plots (Table S1) and inside polygons fit around all pairwise combinations of climatic predictors (Figure 1a In order to analyse ECMF community composition, we derived a Sorensen dissimilarity matrix for each pair of sites based on the presence/abence of each OTU. The relationship between geographic distance and difference in climate predictors of each site was analysed using generalized dissimilarity models (GDM) with the package 'gdm' in R. We used the fitted GDM model to project the multivariate ECMF composition in terms of the first three PC axes of predicted community composition from 0 to 255 and plotted the points using three-dimensional colour scaling [PCA 1 (green), PCA 2 (red), PCA 3 (blue)].
In order to identify responses of individual ECMF taxa to changes in temperature gradients, we performed Threshold Indicator Taxa ANalysis (TITAN) using the 'TITAN2' package in R (Baker & King, 2010). First, we aggregated our OTU table to the level of assigned species name and removed all species that occurred in fewer than four sites. Next, we used TITAN to find the individual ECMF species responses to gradients of mean annual temperature, which our analyses identified as having two different ECMF diversity optima.
The analysis returns two metrics for each species: (a) a change point, which splits each species' abundances into two classes along an environmental gradient (above and below a set temperature) in a way that maximizes the fidelity of species' association for one of the two classes; and (b) a standardized z-score, where the magnitude is proportional to the sensitivity of the species to change and the sign indicates that the species' increases or decreases in relative abundance when temperatures exceed the change point (positive and negative z-scores, respectively). For plotting purposes, we display only indicator taxa that are pure and reliable, which are metrics based on the robustness of the sign and magnitude of species response when abundance data are resampled via bootstrapping (e.g. [Baker & King, 2010;van der Linde et al., 2018]). We identified the temperature thresholds for ECMF communities as the peaks in the cumulative distributions of the negative and positive-z values respectively.

| RE SULTS
Climate variables explained 58% of the deviance in ECMF species richness and 41% of variability in species composition among sites (Table 1). K-fold cross validation demonstrates that for both ECMF species richness and community composition that model performance is stable, with coefficients of variation in R 2 values of ~1%.
The most species-rich ECMF communities are associated with sites with high seasonality in temperature and precipitation (Figure 2a).
Seasonality also explains the most variability in ECMF species composition, with the greatest community turnover associated with differences in temperature seasonality (Figure 2b). ECMF respond to historical climates differently by region. The most species-rich ECMF communities in the northern and northwest mountain forests are cold and dry (0°C and <1 m mean annual temperature and precipitation respectively). By contrast, the most

F I G U R E 2
The predictions of partial regression and best-fit splines for general-additive models for species richness (a), with shading along the 95% confidence interval. Total species richness and relative abundance for each site are equal to the sum of the four partial predictions and an intercept value (97.04 for species richness). Points are coloured according to geographic origin as in the map inset on the far right. (b) The expected difference in ECMF community composition, measured as the partial ecological distance, among sites according to generalized dissimilarity models, with lines on the x-axis indicating empirical values, while shading represents the ±1 standard deviation after sampling 70% of sites 10 times. ECMF, Ectomycorrhizal fungi species-rich ECMF communities in eastern temperate forests are hot and wet (>12°C, >1 m precipitation, Figure 2a). Using GAMs, which fit continuous, curvilinear responses along temperature gradients, we found a bimodal relationship between species richness and mean temperature, with separate cold-and hot-diversity optima (Figure 2a).
As a result, our models predict that warming decreases species richness in the relatively cold north/northwest forests and increases species richness in the eastern temperate forests (Figure 4a,b).
Similar to our models of ECMF species richness, the relative abundance of long-distance foraging strategies increases with mean temperature in eastern temperate forests with mean annual temperatures >12°C and declines with mean temperature in southern boreal forests ( Figure S7c). By contrast, both short-and medium-distance foragers increase with rising temperatures (Figure S7a,b). As a result, our models predict that sites that lose species with increased temperatures should also decline in the abundance of long-distance foraging strategies.
TITAN models identify ECMF species that reliably respond either negatively or positively to increases in mean annual temperature ( Figure 3a). Based on the peaks of the cumulative distributions of the change points for ECMF species with negative and positive temperature associations (Figure 3b), we identify two ECMF community thresholds: a cold threshold at 3°C, which is associated with declines among threshold indicator species with negative z-scores; and a hot threshold at 12°C, which is associated with increases among threshold indicator species with positive z-scores. Notably, these ECMF community thresholds occur near the inflection points for the response of overall ECMF species richness to mean annual temperature (Figure 2a).
We compared model predictions using 30-year mean climates , Figure 4a,b) to climate projections for 2070 both with and without mitigation in greenhouse gas emissions (RCP 4.5 and 8.5 respectively). Our models predict that both ECMF diversity and the abundance of long-distance foragers decline at the temperateboreal ecotone (Figure 5a,b). With mitigation in greenhouse gas emissions, median loss of species richness shrinks from 26% to 18%.
By contrast, our models predict increases in ECMF species richness in south/southeastern forests, including the southern extent of the eastern temperate (Figure 5a,b).

| D ISCUSS I ON
We predict climate changes in the next 50 years will result in a median loss of ECMF species richness of 21% in Pinaceae forests throughout a boreal-temperate climate zone comprising more than 3.5 million square kilometres of North America (an area twice that of Alaska state). The boreal-temperate ecotone is already vulnerable to global change due to warming-associated declines in growth and photosynthetic rates of southern boreal trees (Reich & Oleksyn, 2008) and predicted declines in the dominance of ectomycorrhizal trees . Our results suggest that this same transitional region will experience a loss in below-ground diversity and relative abundance of long-distance foraging ECMF species. Moreover, even if greenhouse gas emissions are mitigated-reaching a peak around 2040, followed by a decline into 2100-ECMF species declines may be unavoidable.
Together, these results frame the urgency with which ecologists must fill the knowledge gap separating ECMF community composition and diversity with ecosystem processes.
Our models explain a substantial amount of variability in ECMF diversity (58%) and community structure (48%) using only four climate variables, with soil N-availability and anthropogenic N-deposition being dropped as predictors during model selection due to their low predictive power relative to model complexity (Supplementary F I G U R E 3 (a) The change points (circles) and 95% confidence intervals for ECMF species with negative and positive responses to increasing mean annual temperature (−z and +z respectively). (b) The sum of z-scores (lines with points, left axis) and cumulative distribution of change points (right axis) for ECMF species with negative and positive z-scores. Peaks in the unfiltered sum(z) and sharp increases in cumulative frequency indicate ECMF community thresholds for change along mean annual temperature gradients. ECMF, Ectomycorrhizal fungi Materials). This negative result with respect to N-deposition contrasts with some regional studies of ECMF (Batstone, Dutton, Wang, Yang, & Frederickson, 2017;Jarvis et al., 2013;Pardo et al., 2011;Suz et al., 2014), including a recent continental-scale analysis of ECMF community composition across western Europe (van der Linde et al., 2018). Possibly these differences reflect the relatively steeper climatic (and shallower N-deposition) gradients in North America relative to Europe. While we acknowledge that there are many potential drivers of ECMF community structure, our study design, which focuses exclusively on ECMF in Pinaceae-dominated forest stands, allowed us to isolate the large and regionally divergent responses of ECMF communities to spatial-temporal climate gradients.
The most diverse EMCF forests in our network had highly seasonal temperature and precipitation. Seasonal forests are also associated with ephemeral flushes of nutrients (Voříšková, Brabcová, Cajthaml, & Baldrian, 2014), which ECMF can rapidly absorb, store in networks of soil mycelia and transfer to tree hosts at later times (Read, 1991). The higher predicted ECMF diversity associated with seasonality in precipitation mirror results from experimental manipulations on fungal communities (Hawkes et al., 2011) and suggests F I G U R E 4 Predicted (a) ECMF species richness and (b) ECMF community composition for North American Pinaceae forests using the most recent 30-year climate data. For (a) species richness, the colour legend also shows the cumulative area containing ≤ the same predicted number of species. For community composition (b), pixels with similar colours are predicted to have similar ECMF community compositions, according to the first three axes of PCA of generalized dissimilarity model predictions, which are broken into constituent subplots with red, green and blue colour scaling. ECMF, Ectomycorrhizal fungi that ECMF diversity may be associated with a storage effect of seasonal specialists (Chesson, 2000).
The regionally opposite effects of temperature gradients we observed for ECMF diversity and composition are consistent with regionally opposite effects of climate on host-tree physiology. We found that ECMF diversity and long-distance forager abundance decline with increasing mean temperatures, but only in relatively cold sites from northern and northwest mountain forests. Warming of boreal tree species growing near the boreal-temperate ecotone has also been shown to reduce tree growth and photosynthetic rate (Reich & Oleksyn, 2008;Reich et al., 2015), which causes trees to allocate less C to ECMF (Fernandez et al., 2017). By contrast, simulated warming does not result in declines in photosynthetic rates for temperate tree species, which are adapted to warmer climates, or in boreal tree species growing at their colder, northern range limits (Reich & Oleksyn, 2008 ECMF differ in foraging and dispersal strategy and enzymatic abilities, which can lead to a positive relationship between ECMF diversity and seedling growth from a range of 1-4 ECMF species in experimental tree seedlings (Baxter & Dighton, 2001. Recent work suggests that species loss results in declines in ecosystem productivity and resiliency across a broad range of taxa (Duffy et al., 2017), including plant-microbial symbionts (Laforest-Lapointe et al.,

F I G U R E 5
The projected changes in (a) ECMF species richness and (b) the % in long-distance (ML) foraging ECMF abundance using projected future climates for the year 2070 both with and without mitigations in greenhouse gas emissions (RCP 4.5 and 8.5, respectively). For both (a) and (b), the colour legends to the right also show the cumulative area predicted to change in % by ≤ the amount of the y-axis (dotted and solid lines for RCP 4.5 and 8.5, respectively), with the line segment at the top showing total land area of Alaska (AK) state. ECMF, Ectomycorrhizal fungi; RCP, Relative Concentration Pathways 2017). However, it remains unclear how a system with ~100 ECMF species will respond to a 25%-30% loss or gain in ECMF diversity.
Given that our study and others (van der Linde et al., 2018) detect major shifts in the diversity and composition of ECMF at the continental scale, understanding the contribution of ECMF communities to tree physiology and ecosystem function is a pressing knowledge gap.
Forecasting continental changes in ECMF communities is a first step in projecting the functional consequences of those changes.
Our models identify the boreal-temperate ecotone, which is already identified as a tipping point in the Earth system, as being particularly sensitive to warming-associated diversity declines, as it exists in a region where increases in mean annual temperature can drive ECMF communities into the valley between cold-boreal and warm-eastern temperate diversity optima. Based on our continental-scale models of ECMF-climate relationships, we predict substantial and regionally divergent trajectories for North American ECMF communities under future climates.

ACK N OWLED G EM ENTS
We would like to thank the members of Peay (Stanford University) and Zhu laboratories (University of California Santa Cruz) for constructive comments on this manuscript. This work was supported by NSF BIO 1926438 and NSF DBI 1045658, 1046115, and 1046052.

CO N FLI C T O F I NTE R E S T
The authors declare that they have no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
Data are archived in Dryad: https ://doi.org/10.5061/dryad.h9w0v t4dr (Steidinger, 2019).