Metatranscriptomics captures dynamic shifts in mycorrhizal coordination in boreal forests.

Carbon storage and cycling in boreal forests-the largest terrestrial carbon store-is moderated by complex interactions between trees and soil microorganisms. However, existing methods limit our ability to predict how changes in environmental conditions will alter these associations and the essential ecosystem services they provide. To address this, we developed a metatranscriptomic approach to analyze the impact of nutrient enrichment on Norway spruce fine roots and the community structure, function, and tree-microbe coordination of over 350 root-associated fungal species. In response to altered nutrient status, host trees redefined their relationship with the fungal community by reducing sugar efflux carriers and enhancing defense processes. This resulted in a profound restructuring of the fungal community and a collapse in functional coordination between the tree and the dominant Basidiomycete species, and an increase in functional coordination with versatile Ascomycete species. As such, there was a functional shift in community dominance from Basidiomycetes species, with important roles in enzymatically cycling recalcitrant carbon, to Ascomycete species that have melanized cell walls that are highly resistant to degradation. These changes were accompanied by prominent shifts in transcriptional coordination between over 60 predicted fungal effectors, with more than 5,000 Norway spruce transcripts, providing mechanistic insight into the complex molecular dialogue coordinating host trees and their fungal partners. The host-microbe dynamics captured by this study functionally inform how these complex and sensitive biological relationships may mediate the carbon storage potential of boreal soils under changing nutrient conditions.

Complex microbial communities (microbiomes) are intrinsically associated with all higher organisms and can contribute functions spanning a continuum, from essential mutualistic services to parasitism (1). DNA-sequencing approaches have greatly increased our understanding of the complexity of microbiomes, from environments as diverse as the human gut (2), to oceans (3) and soils (4). Recent sequencing approaches to measure gene expression within microbial communities (metatranscriptomics) have enabled researchers to complement these taxonomic surveys with real-time functional insights in a range of clinical (5, 6) and ecological (7, 8) settings. Metatranscriptomics can also capture the host transcriptome, providing insight into both community functions and their interaction with active host processes (9, 10). To date, the technical difficulties associated with metatranscriptomics have largely limited its application to controlled laboratory experiments or to interactions between a host and select symbiotic partners. These experiments are essential to generate detailed molecular-level mechanistic understanding of the interaction between single microbial species and their host, but fail to provide understanding of how those species will interact within the complex communities present in nature.

Soil microbes, such as fungi, are primary mediators of carbon (C) and nitrogen (N) cycling in terrestrial ecosystems, and thus have a critical influence on nutrient availability and soil C sequestration and storage (11, 12). In addition, fungi directly interact with over 85% of all vascular plant species via mycorrhizal symbiosis (13), a form of mutualism that has shaped the evolution and expansion of terrestrial plant life. The success of this symbiotic association has resulted in the frequent adoption of mycorrhizal strategies throughout fungal evolution, with current estimates suggesting 40,000 to 50,000 fungal species can form mycorrhiza with plant hosts (14). In doing so, these fungi have experienced varying degrees of diminished autonomy, as they have successively lost biological functions connected to their pathogenic or saprotrophic origins (15). For these kingdom-spanning associations to establish, a complex molecular dialogue must occur to structurally and physiologically coordinate the host with its prospective fungal partners. A crucial component of these interactions are microbial effectors: small secreted proteins deployed by microbes (both mutualists and pathogens) that promote colonization via modulation of host defense mechanisms (16), alteration of cellular architecture (17), or appropriation of host physiology (18).

In the nutrient-limited (NL) soils of boreal forests, ectomycorrhizal (ECM) fungi predominate, and play vital roles in mobilizing mineral nutrients, such as N and phosphorus (P), exchanging them with plant roots for energy-rich C compounds, such as simple sugars (19). Currently, nutrient availability is a critical limitation of boreal plant productivity and C storage, and a restraining force maintaining host-plant dependency on their ECM partners (20). Anthropogenic N deposition and higher rates of mineralization linked to global warming are predicted to diminish this dependency and reduce growth limitations (21), leading to decreases in the proportion of photosynthates allocated belowground and increases in aboveground vegetative C storage, ranging from 5 to 30 parts of C sequestered per unit of additional N (22–24). Until recently, aboveground plant litter was regarded as the primary origin of soil-sequestered C. However, a recent study using bomb-14C modeling showed that at least half of the stored C in boreal soils is derived from plant roots and fungal mycelium (25). As these forests comprise the largest terrestrial biome and harbor more stored C than temperate and tropical forests combined (26, 27), alterations to the diversity, composition, and function of soil- and root-associated microbial communities resulting from increased nutrient availability, particularly N, could be expected to impact ecosystem productivity and stability, with global consequences for C storage and release (28–30).

Our capacity to predict the impact of climate change and anthropogenic influences on a range of ecosystems—and the effect this will have on the global C budget—requires rigorous understanding of the complex microbial communities involved and their interactions with host species. An essential component of this is to understand the diverse biological processes these microbial communities contribute to the global C budget, and how these processes might be altered by changes in symbiotic coordination. Variations in the host–soil microbiome dynamic can be driven by environmental perturbations, either via direct anthropogenic inputs, such as fertilization, or indirect inputs, such as increased atmospheric N deposition or enhanced mineralization rates in warming soils, and can vary along natural nutrient or hydraulic gradients. To begin to develop mechanistic insight into how changes in symbiotic coordination may occur in environments as complex as forest soils, we generated high temporal resolution metatranscriptomes to profile Norway spruce [Picea abies L. (Karst.)] fine roots and over 350 root-associated fungal species over a growing season from NL) and long-term nutrient-enriched (NE) plots in the boreal region of northern Sweden (Materials and Methods and ) and developed robust and reproducible data-analysis pipelines that can be applied broadly to study host–microbiome coordination. This resource is publicly available as the Boreal Rhizospheric Atlas (BRA) online tool to facilitate community utilization (https://www.boreal-atlas.info/). Our findings link the mechanisms by which trees limit C flow and augment defense processes as nutrient availability increases with the profound effect this can have on fungal community structure, function, and host–microbe coordination.

Results

Alterations to Fungal Community Structure.

Field-sampled fine roots of Norway spruce contained RNA originating from the host tree and associated fungal communities. Sequencing this RNA (RNA-seq) revealed a significantly greater relative proportion of fungal RNA-seq reads in samples from NL compared to NE treatment (Fig. 1) (P < 0.001; Wilcoxon signed-rank test), and supporting microbial phospholipid fatty acids data collected from soils at this site, which showed decreased fungal abundance in NE soils (31). NE conditions had a profound effect on the fungal community transcriptome at all time points throughout the growing season (349 up-regulated and 573 down-regulated differentially expressed [DE] Kyoto Encyclopedia of Genes and Genomes orthologs across the season; adjusted P < 0.01) while there was no apparent signal of seasonality (Fig. 1). Taxonomic profiling of the data revealed shifts in community structure that were confirmed in the subsequent year using metatranscriptomics and DNA ITS amplicon data (2012) ().

Fig. 1.

Transcriptomic overview of root-associated fungal community structure in samples from NL and NE plots from the same forest site over a growing season. (A) Relative distribution of fungal reads observed per time point in NL and NE conditions. An asterisk (*) indicates significant difference between NL and NE (P < 0.001; two-tailed Wilcoxon signed-rank test). (B) Vector plot of the PCA of a fungal metatranscriptome based on samples collected at 19 time points over a growing season, in NL or NE soils. Each vector denotes the distance and orientation between the average of the NL replicates for a given time point (dot at the origin of the arrow), and those of the NE samples (point of the arrow). Weeks with active fertilization occurring have been indicated by a dotted vector. (C) The number of reads assigned to each fungal family were normalized to the total number of fungal reads present and the top 10 most represented fungal families in each treatment were identified (11 families in total), accounting for 98.7% of all fungal reads. The normalized read count of these 10 families were plotted across the sampling time course for both NL and NE conditions. The period of active fertilization has been indicated by dotted shading. (D) Shannon diversity index of fungal families in NL and NE forest plots. An asterisk (*) indicates significant difference between NL and NE (P < 0.001; two-tailed Wilcoxon signed-rank test). (E) PCA of fungal reads mapped to specific reference genomes (fungal species) in NL and NE conditions. The nine fungal species with the highest read representation are indicated. PC1 was largely attributable to the number of reads assigned to a specific species (86% variance explained) and was visualized using node size. Node color represents species identity and the colored line around each node indicates treatment, with green = NL, gold = NE. (F) Relative hierarchy of the 11 most represented fungal species in NL and NE, representing 49.9% of all fungal reads. An asterisk (*) indicates this species is not in the corresponding list in the opposing treatment.

Briefly, 150 fungal families were identified associated with NL roots, with the most abundant belonging to the phylum Basidiomycota, of which the families Cortinariaceae (27.6%) and Atheliaceae (22.6%) predominated, followed by the Ascomycetes Hyaloscyphaceae (3.4%), and Gloniaceae (1.6%) (Fig. 1). The majority of reads within each of these families were assigned to a single reference genome: 1) Cortinarius glaucopus (Cortinariaceae), 2) Piloderma olivaceum (Atheliaceae), and 3) Cenococcum geophilum (the only species of the Cenococcum genus; Gloniaceae). For Hyaloscyphaceae, reads were primarily assigned to two reference genomes: Meliniomyces bicolor and Meliniomyces variabilis (synonymy: Hyaloscypha). These taxonomic assignments were conditional to a ≥85% similarity to the available reference genomes () and may represent a number of closely related species belonging to the same genus. For example, amplicon data () identified 10 swarm operational taxonomic units assigned as Piloderma and over 50 as Cortinarius. Despite this caveat, for the remainder of this article the full species names of the best-matching reference genomes are used to avoid confusion.

Fungal family diversity increased significantly in NE conditions (P < 0.001; Wilcoxon rank sum test) (Fig. 1), accompanied by decreased read representation of formerly dominant Co. glaucopus (from 27.6 to 7.0%) (Fig. 1). Similarly, NE conditions resulted in a ∼50% reduction in read abundance assigned to P. olivaceum; however, this remained the most dominant species with the reduction reflecting a decrease in overall fungal abundance rather than a complete alteration in fungal abundance hierarchy (Fig. 1 ). While reads assigned to these Basidiomycota decreased in NE conditions, read abundance and species dominance increased for the Ascomycota Ce. geophilum (from 1.6 to 6.8%) and Meliniomyces spp. (from 3.4 to 4%) (Fig. 1).

Remodeling of Fungal Functions.

While conventional amplicon-based approaches assayed community structure (), metatranscriptomics revealed the functional impact of community restructuring within a coexpression network. Broadly speaking, the network divided into two hemispheres: 1) genes with higher average transcript abundance in NL conditions and 2) genes with higher average transcript abundance in NE conditions (Fig. 2 and ). Hemisphere 1 was overrepresented in reads assigned to Co. glaucopus or P. olivaceum (P < 0.001; one-proportion z-test), with functional enrichment of cell cycle processes and meiosis, amino acid and secondary metabolite biosynthesis, and signaling pathways (e.g., AMPK, MAPK, TGF-β, ErbB, and calcium signaling pathways; adjusted P < 0.05 in all cases; Fisher exact test) (). In contrast, hemisphere 2 was overrepresented in reads assigned to Ce. geophilum or Meliniomyces spp., and enriched in functions indicative of carbohydrate starvation in fungal symbionts (32), such as autophagy (e.g., lysosome, glycan degradation, protein processing in the endoplasmic reticulum) and alternate energy pathways (e.g., fatty acid metabolism, glyoxylate cycle, and gluconeogenesis). Functions linked to primary energy metabolism (e.g., pyruvate metabolism, tricarboxylic acid cycle, and oxidative phosphorylation) were also significantly enriched. Taken together, these changes indicate that although exhibiting transcriptional hallmarks of starvation, the metabolic plasticity of species within hemisphere 2 allows them to effectively utilize alternative energy sources (e.g., lipids) and persist in NE soils.

Fig. 2.

Remodeling of fungal community metatranscriptome in samples from NL and NE plots from the same forest site over a growing season. (A) Visualization of the root-associated fungal community network of samples collected from NL and NE forest plots over a growing season. Average expression profiles of the eight largest network modules are displayed, with NL in green and NE in gold, and the SD indicated by gray shading. Significantly enriched functional categories (adjusted P < 0.05; Fisher exact test) for each network module can be found in . The network modules were grouped into two hemispheres, with hemisphere 1 (indicated in green) displaying higher average expression in NL conditions, and hemisphere 2 (indicated in gold) displaying higher average expression in NE conditions. Fungal species that are significantly enriched in each hemisphere (P < 0.0001; one-proportion z-test) are indicated. (B) Low-dimensional visualization of all fungal genes in both NL and NE conditions. The gene nodes have been colored based on the fold-change differences between NL and NE conditions. Where possible, the branching arms have been assigned to distinct taxonomic assignments, with species names reflecting the reference genome to which RNA-seq reads were assigned. (C) The genes of the three fungal species with the highest assignment of RNA-seq reads (mapped to a specific reference genome) were isolated and separated based on three criteria: 1) fold-change was <−0.5 (higher abundance in NL), indicated in green; 2) fold-change was >0.5 (higher abundance in NE), indicated in gold; and 3) the fold-change was < 0.5 but >−0.5. The average expression profiles of the genes meeting these three criteria have been visualized and corresponding functional enrichment analysis can be found in .

Taxon-specific expression profiles (Fig. 2 ) and associated functional modifications () were identified for the three species with the highest representation in the fungal read pool (P. olivaceum, Co. glaucopus, and Ce. geophilum). In response to NE conditions, P. olivaceum and Co. glaucopus underwent extensive transcriptional modulation, with decreased transcript abundance of functions linked to cell growth, nuclear pore complexes, and proteolytic processes, reflecting a decrease of growth in these species. In contrast, the impact of NE conditions on the transcriptome of Ce. geophilum was comparatively subtle, indicating that the observed increase in community dominance did not result from a profoundly altered lifestyle strategy. Functional enrichment analysis of genes that remained stable between treatments (70% of Ce. geophilum genes) revealed a high level of similarity with functions that were prevalent in P. olivaceum during NL but lost in NE. These functions were associated with metabolism and growth: for example, the degradation of polycyclic aromatic hydrocarbons (such as lignin), aromatic amino acid metabolism, and cell growth. Unique to Ce. geophilum in both conditions was the production of bioactive molecules with structural functions, such as melanin biosynthesis (tyrosinase and polyketide synthase); intracellular receptors of fungal innate immunity, which play a role in the detection of nonself and mediating symbiotic cross-talk (signal transducing ATPases containing WD, ANK, and TPR repeat domains) (33); and antimicrobial functions (nonribosomal peptide [NRP] synthetase component F and NRP synthetase/α-amidoadipate reductase, isopenicillin N synthase, and phospholipase A2). This reflects Ce. geophilum navigating interactions with competing microbes (34, 35) and suggests its success under NE conditions is a combination of metabolic versatility, environmental resilience, and chemical suppression of microbial competitors.

Norway Spruce Root Transcriptome.

The temporally resolved Norway spruce root transcriptome revealed extensive seasonal remodeling common to both NL and NE conditions (Fig. 3) (4,671 up-regulated and 4,999 down-regulated DE genes across the season; Benjamini–Hochberg adjusted P < 0.01; absolute log2 fold-change 0.5). Functional enrichment tests for sets of coexpressed genes within the inferred coexpression network (network modules) were carried out using two complementary ontology-driven approaches: 1) gene ontology (GO) (36), with adjusted P < 0.05; parent–child adjusted significance test (37); and 2) Mapman (38), with adjusted P < 0.05, Fisher exact test. This analysis revealed seasonal transcriptomic alterations primarily indicative of growth, including functions associated with: 1) cell cycle, cell proliferation, and cytoskeleton microtubular networks; and 2) cell wall biosynthesis and remodeling, which coincided with warming soil temperatures, increases in water availability, and a maxima in daylight hours () (network modules 3 and 5). In contrast, the impact of NE had more subtle effects on the root transcriptome than the effect of seasonality (1,689 up-regulated and 956 down-regulated DE genes; NE vs. NL).

Fig. 3.

Transcriptomic overview of Norway spruce fine roots sampled from NL and NE plots over a growing season. (A) Vector plot of the PCA of Norway spruce fine roots collected at 19 time points over a growing season. Each vector denotes the distance and orientation between the average of the NL replicates for a given time point (dot at the origin of the arrow), and those of the NE samples (point of the arrow). Weeks with active nutrient enrichment are indicated by a dotted vector. (B) Visualization of Norway spruce fine root coexpression network of samples grown on NL and NE forest sites over a growing season. The network has been clustered into distinct network modules and key modules mentioned in the text are indicated. (C) The average expression profiles of the network modules mentioned in the text are indicated for NL (green) and NE (gold), and the SD is indicated by gray shading. A dot plot illustrates significantly (adjusted P < 0.05) enriched GO and MapMan functional categories for these network modules. Average expression profiles and significantly enriched functional categories for all network modules can be found in .

Broadly speaking, NE conditions did not fundamentally alter seasonal expression, impacting instead the magnitude of the expression dynamics. For example, in both NL and NE samples, genes belonging to network modules 2 and 12 were characterized by a series of transient bursts in transcript abundance across the season, with the most extreme occurring at week 26 (June 27 to July 3). Although temporally coincident, the magnitude of this induction was greater in NE samples (Fig. 3) (network modules 2 and 12). These modules were significantly enriched in functions associated with signaling and mounting defense responses, such as effector-triggered immunity signaling (nucleotide-binding leucine-rich repeat [NLR] effector receptors and protein kinases), cell wall biogenesis and thickening, oxidative stress, and jasmonic acid biosynthesis. Dehydration-responsive element binding (DREB)-type transcription factors were also significantly enriched in these modules, and overexpression of this class of transcription factor has been linked to enhanced biotic and abiotic stress tolerance in plants (39). Notably, the transient burst in expression of genes encoding defense responses at week 26 coincided with conditions that support novel fine root production and ECM colonization, such as the first substantial rainfall of the year and a period of rapid soil warming (). In contrast to the elevated expression of stress-associated functions, a significant decrease in transcript abundance was observed in NE for a suite of sugar efflux transporters, primarily the Sugars Will Eventually be Exported Transporters (SWEET) family of transporters (P < 0.001; two-tailed Wilcoxon signed-rank test) (). These uniporters regulate the passage of sucrose, glucose, and fructose across cell membranes and are one of the primary exporters of C into the rhizosphere (18, 40). Furthermore, their promoters are targets for microbial transcription activator-like effectors, leading to the appropriation of plant sugar transport during infection (18, 41, 42). The lowered expression of sugar transporters coupled with enhanced defense and effector-triggered immunity responses observed under NE conditions (Fig. 3 and ) (network modules 2 and 12) represent converging mechanisms by which host trees redefine their economic relationship with the fungal community in response to increased nutrient availability.

Changes in Symbiotic Coordination.

As the extensive seasonal dynamics observed in the tree did not appear to drive community-level fungal processes (Fig. 1), we carried out a targeted analysis of the three most dominant fungal species to uncover specific signals of symbiotic coordination masked by global analyses. For both P. olivaceum and Co. glaucopus, all of the fungal network modules that were strongly correlated with Norway spruce network modules in NL conditions lost correlation in NE conditions (17 and 18 module comsbinations, respectively; Spearman’s rank correlation coefficient [R] > 0.7 or < −0.7) (Fig. 4). In both genera, these fungal network modules were enriched in functions observed in hemisphere 1, such as signaling, cell cycle processes, and the metabolism of amino acids and secondary metabolites. The most connected Norway spruce module that was decoupled from Co. glaucopus in NE conditions (Norway spruce network module 10) was significantly enriched in plant solute transport functions: for example, NITRATE TRANSPORTER 2 (NRT2), PHOSPHATE TRANSPORTER 1 (PHT1), AMMONIUM TRANSPORTER (AMT2/3-type), multiple phospholipid and fatty acid transporters, and a range of cation transporters associated with Fe, Cu, Ca, and K transport. Notably, this Norway spruce module became coordinated with two P. olivaceum network modules in NE (P. olivaceum network modules 2.1.3 and 2.1.5), which were enriched in fungal functions linked with glycolysis, nucleobase transport, peptidase activity, and the glyoxylate cycle (). The contrasting responses observed for these two dominant ECM genera could explain P. olivaceum’s persistence in relative community dominance in NE (Fig. 1). Alternatively, the observed fitness of P. olivaceum in both nutrient conditions could signify the mapping of several closely related Piloderma spp. to the P. olivaceum reference genome. In this scenario, differences in physiology or exploration strategies among these species could explain the perceived fitness of the single reference species.

Fig. 4.

Analysis of the coordination between the host-tree and the three most dominant root-associated fungi. (A) Distribution of coordinated network modules (Rs > 0.7 or < −0.7) between Norway spruce and specific fungal species in NL and NE conditions. The number of coordinated network module combinations between Norway spruce and a given fungi are tabulated in the top left of the figure. Those occurring in NL conditions only are shown in green, those occurring in NE conditions only are shown in gold, and those occurring in both treatments are shown in black. (B) Visualization of coordinated Norway spruce and fungal network modules. Edge coloring is consistent with that used in the tabulated section above. The size of each node represents the number of undirected edges to which it is connected. The term “undirected” means the edges represent a two-way relationship between nodes. Network modules that are mentioned in the text have been indicated. Inclusion of an ‘i’ in a network module name indicates the module has the inverse expression profile of the named module. The difference in tree height and aboveground biomass illustrated is modified from Sigurdsson et al. (54) and Majdi and Andersson (117). Functional enrichment analysis for network modules for these species can be found in .

Ce. geophilum displayed a fundamentally different response from the two Basidiomycetes, with 89% of the combinations with Norway spruce network modules that were strongly coordinated in NL conditions maintaining this coordination in NE (16 module combinations), while more than double this number became coordinated only in NE (Fig. 4) (Ce. geophilum). The modules that maintained or increased coordination (primarily Ce. geophilum network modules 4, 14, 21, and 28) were significantly enriched in functions associated with metabolism (fatty acid desaturases, glycogen synthase, and pyrimidine biosynthesis), cell wall remodeling (glucan synthase and endoglucanases), nutrient transport (ammonia and nucleoside permeases), and signal transduction (histidine kinases) (). Notably, many of the novel module combinations formed in NE conditions included Norway spruce modules that were significantly enriched in plant functions associated with cell wall remodeling, and the biosynthesis and perception of presymbiotic molecular signals: specifically, Norway spruce network modules 14 and 15i, which were coordinated with P. olivaceum in NL but lost coordination in NE conditions, and Norway spruce network module 16. These functions promote reciprocal physical and physiological changes that facilitate symbiosis and included the biosynthesis of terpenes, diterpenes, and strigolactone signal transduction, in addition to a nodulation signaling pathway 1 (NSP1) component, a GRAS transcription factor previously thought to be exclusively linked to nodulation but more recently implicated in mycorrhizal establishment (43, 44).

Coordination between Fungal Effectors and Norway Spruce Genes.

In light of these findings, we targeted the coordination between Norway spruce genes and predicted effectors of the three most dominant fungal species (Fig. 5 and ). Using a combination of predictive tools to establish protein secretion (45, 46) and a machine-learning classifier trained on experimentally supported fungal effectors (47), 33 predicted effectors were identified in P. olivaceum (representing 0.3% of all reads mapped to P. olivaceum), 8 in Co. glaucopus (0.5%), and 23 in Ce. geophilum (0.8%). Putative functions were assigned to the predicted fungal effectors based on the identification of protein domains and transit peptides, which were primarily associated with the modulation of gene expression (transcription factors and subunits of the Mediator complex); redox sensing and signaling; and protein folding, posttranslational modification, and degradation (). Intrinsic disorder, a property that causes proteins to lack a fixed three-dimensional structure while retaining biological activity (48), was the most frequently predicted domain identified in the effectors of all three fungal species.

Fig. 5.

Functional analysis of Norway spruce genes displaying positive and negative coordination with predicted fungal effectors. (A) The present set of genes assigned to Co. glaucopus, P. olivaceum, and Ce. geophilum reference genomes were assessed using the fungal effector predictive tool EffectorP (Materials and Methods). The hierarchically clustered transcript abundance of the predicted fungal effectors is displayed for each species. (B) Spearman’s rank correlation was carried out between the expression profiles of the predicted fungal effectors and all Norway spruce transcripts, in both the NL and NE seasonal time courses. The number of Norway spruce genes that satisfied the given Rs (r > 0.7 or < −0.7) and adjusted P thresholds (<0.001) with a given fungal effector are displayed. Genes positively correlated in NL conditions are displayed in light green, negatively correlated in NL are displayed in dark green, positively correlated in NE are displayed in light gold, and negatively correlated in NE are displayed in dark gold. Significant (adjusted P < 0.05) functional enrichments for these gene groups are provided (). Specific coordinated effectors discussed in Results are marked in bold.

A number of significantly overrepresented Norway spruce functions were strongly coordinated with the predicted effectors of all three fungal species, regardless of soil nutrient status. These included transcriptional regulation, vesicle trafficking, pattern-triggered immunity (e.g., fungal and bacterial elicitors), and pathogen effector-triggered immunity (e.g., NLR effector receptors, enhanced disease susceptibility 1 [EDS1], and phytoalexin-deficient 4 [PAD4]) (Fig. 5). In addition, many Norway spruce functions were specifically coordinated with predicted effectors from a single fungal species and showed highly treatment-specific responses. For example, four predicted Co. glaucopus effectors (including putative subunits of the Mediator complex and trafficking protein particle complex) were positively correlated with Norway spruce genes that were significantly enriched in SWEET-family transporters and response-to-symbiont processes (e.g., reduced arbuscular mycorrhization [RAM]2 and Nodulin 26-like intrinsic protein [NIP] transporters) during NL conditions. Notably, for all four predicted effectors, this coordination with sugar efflux and symbiosis responses were lost in NE conditions. A similar loss of coordination in NE soil was observed in a predicted effector of P. olivaceum (Poli84736) containing a cyclin-like subunit (Ssn8) domain of the kinase module of the Mediator complex, which was coordinated with over 3,000 Norway spruce genes in NL conditions—the most of any predicted effector—but only 175 genes in NE soils.

Conversely, the putative effectors of Ce. geophilum demonstrated the greatest level of coordination with Norway spruce genes in NE conditions. For example, a predicted effector (Cgeo64882), containing domains of intrinsic disorder and orthologous to a zinc finger C2H2-type protein in Macrophomina phaseolina, gained 10-times the level of coordination in NE compared to NL conditions. The Norway spruce genes coordinated with this predicted fungal effector were significantly enriched in a large suite of DREB-type transcription factors with homology to specific functionally characterized Arabidopsis thaliana orthologs, which have been linked to the regulation of jasmonic acid biosynthesis (49), repression of disease resistance and mediated cell death (50), and reducing the lignin content of cell walls (51). Taken together, the combination of predictive analytics and high-resolution metatranscriptomics used here provides a strategic approach to identify top candidates for experimental validation, which will dramatically improve our understanding of both the generic and species-specific molecular dialogue occurring between microbial communities and their hosts.

Discussion

Nutrient availability is a critical factor limiting plant productivity and C storage in boreal forest soils. An evolutionary solution to nutrient limitation was the establishment of mutualistic associations, linking plant fitness with root-associated fungal activity. Like most economies, this exchange is sensitive to market forces that can alter the equilibrium of supply and demand (52, 53). Increased nutrient availability, due to anthropogenic N deposition or higher rates of mineralization linked to global warming, is predicted to destabilize this trading status (21), although the severity and lasting consequences of this altered paradigm in boreal symbiosis remains unknown. At the study site used here, 25 y of NE conditions led to a fourfold increase in aboveground tree biomass compared to NL grown trees (54), which coincided with decreases in the flux of photosynthates directed belowground (55), and decreases in autotrophic and heterotrophic soil respiration (56) (). We deployed a high-resolution metatranscriptomic analysis to consolidate these environmental observations with the underlying mechanisms impacted by this symbiotic restructuring. We observed functional evidence of host trees redefining the economic status of their relationship with the fungal community via reduction of sugar transporters that are integral to mycorrhizal establishment (57, 58) (), coupled with an enhanced defense response involving effector-triggered immunity, downstream signal propagation, and stress-associated transcription factors (Fig. 3, modules 2 and 12). The consequence of these alterations on the fungal community was profound, resulting in a reduced relative abundance of sequencing reads of fungal origin (Fig. 1), with a notable decrease in the proportion of reads assigned to ECM Basidiomycota in contrast to an increased proportion of melanotic ECM Ascomycota species (Fig. 1). The sensitivity and resolution of this approach revealed that these changes in community structure were accompanied by, and likely the result of, prominent shifts in broad-scale functional coordination between specific ECM fungi and the host tree (Fig. 4). Perhaps most strikingly, we identified dramatic alterations in the coordination between predicted effectors of the three most prevalent fungal species and over 5,000 Norway spruce transcripts. This finding provides a mechanistic link between the broad structural and functional changes observed in the fungal community, with alterations in the complex molecular dialogue that coordinates the host with its prospective fungal partners.

The response of a single tree to a confluence of effectors from competing fungal species is necessarily an emergent property (59), as individual microbes employ a repertoire of effectors (the “effectome”) that can function pleiotropically (60), redundantly (61), or in concert with each other (62) to affect generic or species-specific processes (63). Global analyses of environmental samples, such as the one described here, can capture these tree-level emergent properties in action, providing crucial insight into their complex in vivo activities and sufficient resolution to identify individual species within the ECM community driving these emergent responses. For Co. glaucopus and P. olivaceum, we observed the transcript abundance of a number of predicted effectors decrease drastically when the resource trading conditions were less favorable in NE conditions, particularly during the period of active fertilization (weeks 23 to 34). In contrast, although the number of Ce. geophilum effectors expressed in NE conditions was the same in NL conditions, the magnitude of their expression was greater in NE conditions (Fig. 5, Ce. geophilum heatmap). Taken together, this suggests the decoupling of P. olivaceum and Co. glaucopus from their host becomes most pronounced following fertilizer application, while the absence of such a response in Ce. geophilum effectors indicates this species is particularly well adapted to N-rich environments.

To date, what little is known of fungal effector function is the result of detailed experimental characterization of individual effectors under controlled conditions. These approaches have provided insights into their diverse modes of action (64), their spatial and temporal kinetics (65), and the host processes commonly targeted (66). Many of the candidate effectors identified in this study possessed qualities linked to interference with protein posttranslational modification (Tyr phosphatase dual domain, protein kinase domain, and Tyr kinase active site), protein stability/folding (DnaJ, calreticulin/calnexin, and peptidyl-prolyl cis-trans isomerase FKBP2/11), or protein degradation (proteasome subunit-β5 and UBQ-conjugating E2). Such mechanisms have been widely described during phytopathogenic interactions (67), for example, during infection the Phytophthora sojae effector Avr1d competitively binds with a soybean E3 ubiquitin ligase, preventing its self-ubiquitination and degradation, leading to increased plant susceptibility to infection (68). Similarly, in Pseudomonas syringae, the effector HopBF1 phosphorylates plant HSP90 to inhibit the chaperone’s ATPase activity, which blocks its initiation of the plant hypersensitive response (69).

We also identified five candidate effectors putatively involved in host transcriptional regulation [e.g., leucine zipper domain, Zn (2)-C6 fungal-type DNA-binding domain, and subunits of the Mediator complex and transcription initiation factor IID]. Compromised immunity arising from effector-mediated interference of host transcription has been described in several plant–microbe interactions, including mutualistic symbioses. For example, the interaction of the symbiosis effectors SP7 of Rhizophagus irregularis with the pathogenesis-related transcription factor ERF19 of Medicago truncatula, and MiSSP7 and MiSSP7.6 of Laccaria bicolor with the jasmonic acid transcriptional corepressor PtJAZ6 and the transcription factors PtTrihelix1 of Populus trichocarpa, respectively (16, 70, 71). Notably, we identified a number of effectors with potential targeting to mitochondria, two of which (one from Co. glaucopus and one from P. olivaceum) possess a domain common to translocase complex proteins of the inner membrane (TIM), which has been linked to mitochondrial-mediated plant immunity in A. thaliana (72, 73). Although mitochondria play a central role in the plant defense response, it is only recently that the effector proteins MoCDIP4 and Avr-Pita from Magnaporthe oryzae were shown to modulate rice susceptibility by targeting a mitochondrial DnaJ protein (74) and a cytochrome c oxidase (75), respectively, revealing a novel pathogen infection strategy targeting host mitochondria. Future efforts to experimentally validate candidate fungal effectors (such as those identified here) (Fig. 5 and ) would be bolstered by the integration of these functional predictions with the foreknowledge of their expression dynamics, not only in regard to treatment and colonization-checkpoints, but in relation to the greater fungal effectome (including self and fungal competitors) and their coexpression with host functions.

The reduction of total fungal reads in NE soils (consistent across growing seasons and methodologies) (Fig. 1 and ) concurs with previous observations that increasing nutrient gradients or active fertilization in boreal soils decreases ECM mycelium and sporocarp production (76, 77). Extramatrical mycelia represent a substantial proportion of total belowground C (78). However, beyond absolute fungal biomass, mycelial taxonomic composition has also been shown to influence the accumulation of C through the production of degradative enzymes (79) and the varying recalcitrance of fungal necromass (25, 80, 81). The Basidiomycota Cortinarius underwent the most dramatic reduction in assigned fungal reads in response to NE conditions (Fig. 1 and ), coincident with a decoupling of broad-scale fungal and tree functions related to growth and solute transport (Fig. 4), as well as highly specific predicted fungal effectors and tree functions associated with supporting mycorrhizal partners (e.g., sugar efflux transporters and symbiont-response processes) (Fig. 5). This genus comprises well characterized ECM fungi with energetically demanding lifestyles that employ medium-distance fringe exploration strategies and secrete oxidative enzymes capable of degrading lignin (79, 82–85). Due to the elaborate foraging strategies of Cortinarius spp., they have some of the highest C demands of all ECM (86) and are thus highly sensitive to reduced C flow from host roots (83, 87–89). The capacity of members of the genus Cortinarius to degrade recalcitrant soil organic matter is unique among ECM fungi and they exert a disproportionate impact on soil C turnover (90). A recent study demonstrated the presence of Cortinarius species-complex members in the organic topsoil of boreal forests resulted in a 33% decrease in local C storage (90). This finding disputes the prevailing consensus that the functional redundancy inherent to complex microbial communities will offset the loss of an individual species. Thus, the collapse in symbiotic coordination and diminished presence of Cortinarius spp. in the fungal community observed in our study suggests that increased atmospheric N deposition and enhanced mineralization rates linked to warming soils will reduce the abundance and activity of these keystone degraders of lignified plant matter and humus in boreal soils.

Contrasting this reduction in Basidiomycete abundance, NE resulted in increases in the Ascomycetes Ce. geophilum (Fig. 1), accompanied by an intensification in transcriptomic coordination with the host tree, at both a broad-scale with fungal processes linked to metabolism, cell wall remodeling, and nutrient transport (Fig. 4 and ), and a fine-scale involving a suite of predicted fungal effectors with putative functions, ranging from transcriptional regulation to protein ubiquitylation (Fig. 5 and ). Ce. geophilum is a globally distributed and highly versatile short-range exploratory ECM fungus with distinctive black hyphae resulting from heavily melanized cell walls (91, 92). This heavy melanization makes these fungi resilient to adverse environmental conditions, such as drought (93, 94), and can act as virulence factors during root colonization (95, 96). Notably, melanized necromass is also highly recalcitrant to decomposition (97, 98), with the mycelia of Ce. geophilum demonstrated to persist in soils up to 10 times longer than that of other ECM species (92, 99). This suggests that while NE resulted in reduced overall fungal biomass (31), it triggered an absolute increase in recalcitrant fungal necromass that is more likely to contribute to long-term C storage. This store of organic C would necessarily benefit microorganisms with the capacity to access it. Functional analysis of Ce. geophilum transcripts highly expressed in both nutrient treatments revealed a significant enrichment of chitinases and functions associated with polycyclic aromatic hydrocarbon degradation () (Ce. geophilum), a cohort of functions that have previously been implicated in the degradation of environmental pollutants and, notably, melanins (100). These enzymes play essential roles in accessing stores of melanized fungal necromass and the self-digestion processes (autolysis) utilized for remodeling of mycelial architecture (80). The expression of these functions in both NL and NE conditions suggests that the capability to oxidatively degrade melanin could be a crucial function for Ce. geophilum, regardless of soil nutrient status, with the corollary that access to these growing reserves of recalcitrant C in NE soil has the potential to further reinforce this species’ position in the rhizospheric community. However, strategic experimentation will be needed to discriminate the overlapping roles these degradative enzymes play in mycelial autolytic remodeling and nutrient scavenging.

As boreal forests harbor globally significant amounts of C, understanding the biogeochemistry of boreal soils is essential for informing strategic responses to future climate scenarios. One area where such understanding is needed was revealed by a recent metaanalysis of over 100 elevated CO2 experiments (12) that has challenged the consensus of existing ecosystem modeling. Previous models had predicted increases in terrestrial C sequestration under elevated CO2 conditions, while this metaanalysis suggests instead that the increased biomass of ECM-associated plants will be supported by augmented nutrient uptake and result in reduced soil organic C accumulation. Recent field experiments have shown that trees exposed to elevated CO2 increased belowground C allocation, resulting in enhanced production of fine roots (101) and higher rates of soil CO2 efflux (102). This alteration in belowground C allocation represents a revealing counterpoint to the impact of soil nutrient enrichment examined in the present study.

CO2-driven increases in belowground C-availability could cause reciprocal changes in the dominance and activity of Basidiomycetes (such as Cortinarius and Piloderma species) at the expense of Ascomycetes (such as Ce. geophilum and Meliniomyces species). Potentially, this CO2-driven shift in community composition could result in a gradual increase in oxidative decomposition of soil organic matter (103) and a decrease in recalcitrant fungal necromass, triggering a reduction in the capacity of boreal forest soils to sequester C but an increase in tree biomass. This conflict between observed and modeled C cycling highlights the uncertainty of climate projections that underestimate the importance and scale of rhizospheric processes and emphasizes the need to monitor and integrate the diverse biogeochemical processes contributed by microbial communities and their host plants. Metatranscriptomic approaches, such as the one described here, can aid in functionally resolving these processes and will add to our conceptual framework of global C sequestration and ecosystem stability in future climate conditions. Notably, this approach advances studies of complex microbial communities from being purely taxonomic, using DNA-based methods, to identifying causative process-driven links and alterations in the transcriptional coordination of mycorrhizal associations. Additionally, minor user adjustments to this opensource and reproducible methodology will enable high-resolution monitoring of a range of host–microbe systems, from pathologies arising from dysbiosis in the human body, to applications in a diversity of threatened ecosystems.

Materials and Methods

Samples were collected at the Flakaliden research site (64°07′N, 19°27′E, altitude 310 to 320 m) in the boreal forests of northern Sweden [for details of the site, please refer to Haas et al. (83)]. Treatment of designated experimental plots within the forest stand of an optimal nutrient solution of macro- and micronutrients—including N, P, K, Mg, Ca, S, Fe, B, Mn, Cu, and Zn—commenced in 1987 and was updated yearly based on analysis of the foliar nutrient status at the site (104). Five root sample replicates were harvested from a single NL plot and a single NE plot (, plots 6A and 7A, respectively) at each time point during 19 wk over the 2011 growing season (May to October). Root samples were harvested from 0- to 25-cm-deep soil, cleaned, frozen, and ground. RNA was isolated from 75 mg of tissue using the Ambion Plant RNA Isolation Aid and the Ambion RNAqueous kit, followed by a DNA-free kit DNase treatment. RNA concentration was determined using a Nanodrop NL100 and RNA quality was assessed using an Agilent 2100 Bioanalyzer. The RNA-seq raw data are publicly available in the European nucleotide archive (ENA) with the accession no. PRJEB35805 (105). In 2012, three sample replicates were harvested from each of three replicate NL plots and three replicate NE plots (, plots 4A, 6A, and 14A; and 1A, 7A, 13A, respectively), at four time points over the 2012 growing season. For these samples, DNA and RNA was isolated. The ITS1 amplicon sequencing data generated from these samples were published previously (83) and reanalyzed by Schneider et al. (106). RNA extracted from these Norway spruce root samples was sequenced as described above. This low temporal resolution collection was used to validate the high temporal resolution samples collected in 2011, in reference to 1) year-to-year variation, 2) interplot variation, and 3) a complimentary ITS amplicon sequencing approach.

Norway spruce RNA-seq data were preprocessed and aligned with SortmeRNA (107), Trimmomatic (108), and Salmon (109). Preprocessing and analysis of metatranscriptomic data were implemented in a reproducible snakemake workflow available on Bitbucket (https://bitbucket.org/scilifelab-lts/n_street_1801/). A detailed description of software, including parameters and versions used in the workflow, has been published in a separate article (106). Assembled fungal transcripts were assigned taxonomic annotation based on sequence similarity to available reference genomes, with ≥85% similarity required for species assignment, ≥60% similarity for genus assignment, and ≥45% for phylum assignment (110, 111). If a similarity threshold to a given taxonomic level was not satisfied, the prefix “Unclassified” was assigned before the next highest classification (e.g., “Unclassified Cortinarius” indicates similarity to the available reference genome was >60% but <85%). Functional annotations were assigned using eggnog-mapper (112) and the eggNOG database. For both Norway spruce and fungi, high-dimensional data analysis was performed using principal component analysis (PCA) and extended with PHATE (113). Differential expression analysis and variance stabilization transformation were carried out with DESeq2 (114). Coexpression networks were constructed using Seidr (115) and clustered using Infomap (116). Gene-enrichment tests were carried out with Gofer3, an in-house tool at the Umeå Plant Science Centre (DOI:10.5281/zenodo.3731544). All enrichment tests can be accessed on the BRA web app (https://www.boreal-atlas.info/). All other analyses were executed using in-house tools. All scripts are available on Github (https://github.com/loalon/flakaliden-2011).