We selected six species of Suillus based on their occurrence in forests in Minnesota, USA, the availability of cultures for genome-sequenced strains, and the availability of historical herbarium data for both the fungi and their associated hosts. These species included S. weaverae EM37, S. spraguei EM44, S. americanus EM31, S. clintonianus FC179, S. ampliporus FC55, and S. luteus UH-Slu-Lm8-n1. All genomes were originally sequenced as part of Lofgren et al., 2021 (Lofgren et al., 2021), except for S. luteus UH-Slu-Lm8-n1, which was sequenced as part of Kohler et al. (2015). Whereas, S. americanus, S. weaverae, and S. spraguei specialize on white pine (Pinus sg. Strobus), S. luteus associates with red pine (Pinus sg. Pinus), and S. clintonianus and S. ampliporus associate with tamarack (Larix) (Nguyen et al., 2016).

The phylogenetic relationship between the six Suillus species was estimated using PHYling v.2 (Stajich and Tsai, 2023), employing the BUSCO conserved single-copy marker set (fungi_odb10). Reconstruction was accomplished using IQTree v.1.6.12, with model selection using ModelFinder (Kalyaanamoorthy et al., 2017). The optimal model based on BIC score was VT+F+R4, which was run with 1000 replicates of Ultrafast Bootstrapping, on a total of 35,550 parsimony-informative sites. AMT gene-tree reconstruction was accomplished by aligning all AMT genes using MAFFT v.7 (Katoh and Standley, 2013), with subsequent alignment trimming using ClipKIT (Steenwyk et al., 2020), removing a total of 34 sites. The best model for AMT gene reconstruction was determined to be JTT+G4 according to ModelFinder, and IQTree was run on a total of 312 parsimony informative sites.

To investigate the genomic capacity of Suillus species to use inorganic N sources, we conducted genome mining and comparative genomic analysis for genes involved in the metabolism of both ammonium and nitrate. Genes involved in ammonium assimilation were identified via homology to the three characterized ammonium transporter (AMT) genes in Hebeloma cylindrosporum (Javelle et al., 2001, 2003) HcAmt1 (GenBank ID AAM21926), HcAmt2 (AAK82416), HcAmt3 (AAK82417). BLASTP searches were conducted as above, first to identify the relevant orthogroups, and then to assign and quality filter individual orthogroup members. After one orthogroup was found to be conspicuously missing a gene from S. ampliporus, an additional whole-proteome search was conducted to determine the loss status of this gene in the S. ampliporus assembly, with BLASTP search parameters set as above.

Genes involved in nitrate assimilation were identified via homology to characterized genes in Laccaria bicolor, including those encoding nitrate transporters (Lbnrt, protein ID 254042), nitrate reductase (Lbnr, ID 254066), and nitrite reductase (Lbnir, ID 291348) (Kemppainen et al., 2010). Initial BLASTP searches were run as above using representative proteins for every orthogroup to identify the above gene families, and then run on all sequences contained in those orthogroups to implement stringent quality filtering using the cutoffs for identifying putative proteases in the MEROPS database (see below). Sequences of each quality filtered orthogroup along with their corresponding L. bicolor reference genes and exons were aligned using MAFFT v.7 (Katoh and Standley, 2013) with the parameters –maxiterate 1000, –localpair, and analyzed using Jalview (Waterhouse et al., 2009).

To investigate the genomic capacity of Suillus species to utilize organic N sources, we conducted genome mining and comparative genomic analysis of genes involved in both and protein and chitin catabolism. Identification of genes involved in protein catabolism was accomplished using orthology to protease type examples housed in the MEROPS database (Rawlings et al., 2018). First, we used OrthoFinder to cluster gene families (Emms and Kelly, 2019), retaining the longest sequence in each family as a representative for each orthogroup. We then removed annotated inhibitors from the MEROPS SCAN library v.12.5 (n = 673 inhibitors) and used the remaining .fasta files (n = 4,335) to create a curated database of non-redundant proteases. We used BLASTP to search for matches between the sequence representatives for each orthogroup and the protease database with parameters e-value = 1e-10 and max_target_seqs = 10. We implemented subsequent QC filtering to retain only results with >40% sequence identity, >70% query coverage, >70% target coverage, and retained only the top hit according to e-value in the case of multiple matches. Annotations were applied to these results by mapping the protein target IDs to MEROPS database annotations.

For the chitin assessment, we annotated CAZymes and Auxiliary Activity (AA) enzymes based on orthology to all currently accepted CazyDB CAZyme classes accessed via the MycoCosm web portal (Levasseur et al., 2013; Grigoriev et al., 2014; Drula et al., 2022). Recent advances have revealed that chitin degradation involves enzymes beyond classically studied families, with members such as GH2, GH9, GH16, GH23, and GH35 demonstrating activity on chitin-derived substrates (Viborg et al., 2019; Vaaje-Kolstad et al., 2019; Zeng et al., 2025) Identification of relevant CAZyme families was accomplished via the manual curation of CazyDB and CAZypedia annotations of all functionally characterized enzymes and associated carbohydrate-binding modules active on chitin, chitosan, and chitin degradation derivatives (CAZypedia Consortium, 2018; Drula et al., 2022). This comprehensive approach captures chitin-active enzymes beyond canonical chitinases and ensures we do not overlook potentially significant enzymatic activities contributing to fungal chitin metabolism. These included the endo-acting enzymes Chitinase (EC 3.2.1.14), Chitin deacetylase (EC 3.5.1.41), Chitin oligosaccharide deacetylase (EC 3.1.1.-), N-acetylglucosamine-6-phosphate deacetylase (EC 3.5.1.25), and the chitin active lytic polysaccharide monooxygenases (LPMOs) encoding Lytic chitin monooxygenase (EC 1.14.99.53). To assess exo-acting enzymes relevant to chitin metabolism, we included Chitin exo-β-1,4-N-acetylglucosaminidase (EC 3.2.1.-), Chitin exo-β-1,4-N-acetylglucobiosaminidase (EC 3.2.1.200), Di-N-acetylchitobiase/reducing-end exo-hexosaminidase (EC 3.2.1.-), [reducing end] exo-chitinase (EC 3.2.1.201), β-1,6-N-acetylglucosaminidase (EC 3.2.1.-), Chitin β-1,3/1,6-glucanosyltransferase (EC 2.4.1.-) and β-N-acetylhexosaminidase (EC 3.2.1.52). To assess chitosan metabolism we included Chitosanase (EC 3.2.1.132), and Exo-β-1,4-glucosaminidase (EC 3.2.1.165).

To evaluate how well our genome-based predictions of N use matched metabolic activity in vitro, we measured biomass production of the same six genome-sequenced Suillus speices when grown on diverse N sources. These sources included inorganic N in the form of ammonium, and organic N in the form of protein, protein-tannin complexes, protein-mineral complexes, and chitin. Cultures of each species were first grown for six weeks on a sterile cellophane membrane covering solid modified Fries media at 21 °C in the dark. The modified Fries media contained 2.72 mM (NH₄)2•tartrate, 2.05 mM MgSO₄·7H₂O, 5.88 mM KH₂PO₄, 0.18 mM CaCl₂·2H₂O, 0.34 mM NaCl, 1, 0.24 mM H₃BO₃, 20 μM ZnSO₄·7H₂O, 5.01 μM CuSO₄·5H₂O, 50.29 μM MnSO₄·H2O, 0.16 μM (NH₄)₆Mo₇O₂₄·4H₂O, 73.99 μM FeCl₃·6H₂O, 55.5 mM d-glucose, 55.51 μM myo-inositol, 0.30 μM thiamine·HCl, 0.10 μM biotin, 0.59 μM pyridoxine·HCl, 0.27 μM riboflavin, 0.82 μM nicotinamide, 0.73 μM p-aminobenzoic acid, and 0.46 μM Ca-pantothenate, supplemented with 0.23 g/L of casein enzymatic hydrolysate and 0.5 g/L of yeast extract (Fries et al., 1987; Op De Beeck et al., 2018). The pH of the medium was adjusted to 4.8 using HCl, and agar was added at a concentration of 12 g/L.

To conduct the N use assay, a single layer of autoclaved glass beads was added into 6-well plates as described in Maillard et al (Maillard et al., 2023). Modified Fries medium (made as above) was used as the base medium, omitting the casein enzymatic hydrolysate and the yeast extract, and with a lower (NH₄)2•tartrate concentration of 0.109 mM. A treatment using the base medium without any additional N supplementation served as a control to account for fungal growth using solely (NH₄)2•tartrate and any N potentially contained in the vitamin and micro-nutrient solutions. The other treatments included the base medium supplemented with (NH₄)2•tartrate (0.460 g/L), bovine serum albumin (BSA) (0.437 g/L), and chitin from shrimp shells (1 g/L). (NH₄)2•tartrate was used as an inorganic N source (Chalot and Brun, 1998). BSA was used as a model for protein degradation, having been previously used to demonstrate organic N acquisition from proteins by ECM fungi (Shah et al., 2013; Wang et al., 2020, 2021; Op De Beeck et al., 2018) and chitin from shrimp shells was used as in Maillard et al (Maillard et al., 2018, 2023) to study the degradation of chitin as a source of organic N. In forest soils, proteins occur not only freely but can also form complexes with phenolic molecules or absorb onto mineral surfaces (Adamczyk et al., 2017), which limits their biodegradability (Knicker et al., 1999; Jilling et al., 2018). Thus, to represent the complexity of protein forms in soils, we also included protein-tannin and protein-mineral complexes as organic N sources. Protein-tannin complexes were prepared as described in Bending and Read (1996), and protein-mineral complexes were prepared as described in Wang et al. (2021) using goethite, a naturally abundant iron oxide in soils as a model mineral. For the protein-tannin and protein-mineral treatments, the BSA concentration was the same as in the free-BSA treatment (0.437 g/L), with almost all BSA considered complexed with tannins or absorbed on geothire based on Bending and Read (1996) and Wang et al. (2021). The respective concentrations of (NH₄)2•tartrate, BSA (either free, complexed, or absorbed), and chitin were chosen to contain the same mass of N per treatment, and to be N-limited with a C:N ratio of approximately 28.6, making the additional C incorporated in chitin, protein, tartrate, or tannic acid negligible relative to glucose C. We added 2.4 mL of each of the five media treatments (low NH₄ control, NH₄, free-BSA, BSA-tannin, and chitin) into the glass-bead-covered wells of 6-well plates. For each of the Suillus species, we cut small mycelium disks (3 mm in diameter) from the cellophane membrane close to the hyphal front and placed one disk in each well (n = 3 wells per treatment). The plates were incubated in the dark at 21 °C and harvested at three timepoints (12-, 24-, and 35-days post-inoculation). The mycelium collected was quickly blotted on a paper towel, then dried at 50 °C, and the dry biomass measured on an analytical balance. To control for variation in the amount of starting material, the dry mass of the initial transplant (day 0) was estimated by averaging 10 transplants per species. To account for any free goethite in the BSA-goethite treatment, all Suillus biomass grown with BSA-goethite as an N-source were combusted at 550 °C for 12 hours in a furnace. Goethite alone was combusted to account for goethite mass loss during combustion. Samples from the ammonium treatment (high biomass) were also combusted to evaluate fungal biomass remaining after 12 hours in 550 °C, and determined to be negligible. The goethite contribution in the original fungal biomass samples was then subtracted to correct the fungal biomass values. Three values were considered outliers (all more than 5 times the average) and subsequently excluded (n=1 S. clintonianus, n=1 S. spraguei and n=1 S. weaverae all at the 35 day time point).

Differences in fungal biomass production across treatments and species were analyzed using a two-way ANOVA with Type III sum of squares, implemented in R (version 4.4.1). To assess relative biomass production across treatments, we first normalized the measurements by subtracting the dry weight of the initial transplant. Then, for each species at each time point, we averaged the biomass generated on the low ammonium control media across n=3 replicates, and subtracted it from each of the treatment biomass measurements for the corresponding time point. Normal distribution of residuals and homogeneity of variance were confirmed using visual inspection. Effect sizes were calculated as eta-squared values (proportion of total variance explained by each factor) by dividing the sum of squares for each factor by the total sum of squares. Post-hoc comparisons were conducted using Tukey’s HSD to identify significant differences between species’ performances within each treatment and between treatments for each species.

In addition to ¹⁵N, the analysis of ¹³C isotopes can provide indirect evidence of N source utilization. The progressive depletion of atmospheric ¹³C compared to ¹²C, due to the release of anthropogenic emissions into the atmosphere is known as the Suess effect (Bassett et al., 2023). These ratios are reflected in plants and other carbon reservoirs, including sporocarps (mushrooms). Historically, plant δ¹³C declined slowly from 1850–1957 during the early industrial revolution, and began declining rapidly starting in 1957 coinciding with globalization, industrialization, and massive increases in fossil fuel consumption (Belmecheri and Lavergne, 2020). By comparing δ¹³C composition of the same species across time, we can also infer N use, where more negative δ¹³C values indicate the assimilation of older, organic, forms of N. We therefore considered both δ¹³C and δ¹⁵N in our isotopic evaluation.

A total of 112 Suillus mushrooms and 131 host plant specimens were selected from the collections at the Bell Herbarium at the University of Minnesota ( Supplementary Table S1 ). All specimens were collected in Minnesota between 1961-2021. Fungal caps and stipes were sampled and processed separately. To conduct elemental and isotopic analysis, air-dried sporocarp and needle tissues were ground and a mass of 2 ± 0.05 mg to the nearest 0.001 mg was analyzed via an Elementar Cube elemental analyzer interfaced to a GV Instruments Isoprime isotope ratio mass spectrometer (Manchester, UK) at Boston University Stable Isotope Laboratory.

Due to normal isotopic fractionation during N metabolism and transport, cap tissue should exhibit lower isotopic values than stipe tissue (Hobbie et al., 2012). As such, samples where δ¹⁵N or δ¹³C were lower in stipe than cap were considered aberrant, potentially indicating contamination, abnormal physiological processes, or sample preservation issues. Sporocarps failing this quality control criteria were therefore excluded from our analysis, resulting in the removal of 5 samples (n=1 S. americanus, and n=2 of each S. clintonianus and S. weaverae). After removing these samples, the final dataset included n = 20 samples of S. americanus, n = 11 of S. ampliporus, n = 24 of S. clintonianus, n = 9 of S. luteus, n = 15 of S. spraguei, and n = 28 of S. weaverae. The host dataset consisted of n = 55 samples of L. laricina, n = 34 of P. resinosa, and n = 42 P. strobus.

For all downstream modeling, values obtained for caps and stipes were averaged for each sporocarp. Mixed effects models were constructed in R to model sporocarp and needle N concentrations, δ¹⁵N, and δ¹³C for both plant and fungal tissues. For all Suillus, to account to non-independence arising from location and seasonality, county and month collected in were included as random effects. Year was included as a fixed effect in all models to test for temporal trends, as well as an interaction between year and species, and species as a fixed effect to account for species-level variation in baseline values. Backwards AIC stepwise selection was performed on all models. Month information was missing for four Suillus collections, which were removed prior to stepwise selection, and reincluded if month was eliminated as a random effect. In many models, stepwise selection eliminated county and month as random effects. These models were then run using simple linear regression via the “lm” function in base R with the preserved variables. To model the δ¹⁵N difference between cap and stipe, sporocarp average N concentration and the difference between cap and stipe N concentrations were included as fixed effects to test the effect of overall N status, and the amount of protein in the cap, respectively.

Phylogenetic analysis of the six Suillus species revealed well-supported evolutionary relationships with all nodes receiving maximum bootstrap support ( Figure 1 ). S. spraguei, S. americanus, and S. weaverae form a strongly supported monophyletic group, while S. luteus occupies a sister position to this trio. In agreement with previous work, all species exhibit distinct evolutionary lineages that correspond to their documented host associations (Nguyen et al., 2016; Lofgren et al., 2021).

N source was the most significant driver of Suillus biomass production (p < 2e-16, F = 139.135, accounting for approximately 64% of the total variation), followed by species (p < 0.001, F = 6.949, ~4% of variation) and the interaction between N source and species (p < 0.01, F = 2.398, ~5% of variation). Overall, S. spraguei had the highest biomass production producing significantly more biomass than all other species (p < 0.05), followed by S. luteus and S. clintonianus. Suillus grown on ammonium produced significantly more biomass than other treatments (p < 0.0001). Among organic N treatments, protein alone supported significantly higher biomass than the other substrates (p < 0.0001), and protein-mineral complexes produced significantly more biomass than protein-tannin complexes (p < 0.01). All species showed a similar hierarchical response to ammonium treatment with more variable hierarchies among organic N substrates. While BSA generally outperformed the protein-mineral complexes and chitin treatments, the relative performance on protein-mineral, protein-tannin, and chitin varied by species. S. weaverae demonstrated the strongest overall response to N source, while S. americanus showed the least differentiation among treatments.

On ammonium, S. spraguei produced significantly more biomass (5.20 mg ± 2.60 SD) than S. americanus (2.86 mg ± 0.98) (Tukey’s HSD, p = 0.021), but no other significant differences were detected on ammonium, despite numerical differences in mean biomass production ( Figure 4A ). When grown on protein, S. spraguei produced the highest mean biomass (an average of 2.58 mg ± 1.68 SD more than controls), followed by S. clintonianus (1.92 mg ± 1.59) and S. ampliporus (1.73 mg ± 1.77). However, all six species displayed positive growth on protein, with biomass production not being statistically significantly different among species ( Figure 4B ). On protein complexed to tannin, all species showed poor or negative growth compared to low ammonium controls. S. americanus maintained the highest biomass production (0.61 ± 0.76 SD), followed by S. spraguei (0.34 ± 0.41 SD), while other species displayed growth inhibition, with S. weaverae and S. ampliporus exhibiting the strongest negative responses ( Figure 4C ). Species varied considerably in their ability to utilize protein absorbed on mineral surfaces, with S. spraguei (1.35 ± 1.66 SD), S. luteus (1.17 ± 1.07 SD), and S. clintonianus (0.85 ± 0.86 SD) displaying positive growth, S. weaverae (0.07 ± 0.95 SD) minimal positive growth and S. americanus (-0.21 ± 0.45 SD) and S. ampliporus (-0.42 ± 0.56 SD) exhibiting growth inhibition compared to controls ( Figure 4D ). On chitin, S. spraguei, S. luteus, S. americanus, and S. clintonianus all displayed slightly positive growth, and S. ampliporus, S. weaverae displayed slight growth inhibition, with the only statistically significant difference between S. spraguei and S. weaverae (p < 0.01) ( Figure 4E ).

Phylogenetic placement of the six Suillus species agreed with previous analyses (Nguyen et al., 2016; Lofgren et al., 2021) showing clear patterns of host specialization ( Figure 1 ). However, phylogenetic relationships did not consistently align with genomic repertoire or N acquisition strategies. For example, while the three P. strobus associates (S. spraguei, S. americanus, and S. weaverae) formed a well-supported clade, they exhibited markedly different capacities for organic N assimilation ( Figure 3 ). These results were echoed in the genomic analyses of all six species, which revealed notable interspecific variation in the genetic machinery for N acquisition across Suillus that did not consistently follow phylogenetic relationships. In contrast, the genes for inorganic N (ammonium and nitrate) were highly conserved and predominantly single-copy ( Figures 2A-C ). In ECM fungi, ammonium uptake occurs via the ammonium importers AMT1, AMT2, and AMT3, characterized in the ECM species Hebeloma cylindrosporum (Javelle et al., 2001, 2003), Amanita muscaria (Willmann et al., 2007), and Tuber borchii (Montanini et al., 2002). While AMT3 is thought to be a low-affinity ammonium transporter, AMT1 and AMT2 function as high affinity transporters, and are subject to N repression (Javelle et al., 2003; Willmann et al., 2007). Here, we were unable to confidently assign the orthogroups to their respective AMT genes, where both orthogroup’s OG0009316 and OG0006824 appeared to be more closely related to the AMT3 reference ortholog HcAmt3, and OG0006295 was more closely related to HcAmt2 ( Figure 2B ). The inability to confidently map the three gene clusters to their corresponding reference sequences highlights both the high sequence similarity between genes in the AMT superfamily, and the need for functional confirmation to distinguish the identity and impact of high- and low-affinity ammonium transporters. Nitrate uptake into fungal tissues is an active process carried out by N transporters, prior to being reduced to nitrite by nitrate reductases and finally to ammonium by nitrite reductase (Takaya, 2002). Compared to other ECM fungal species Suillus generally grows well on nitrate, producing both significantly more biomass and NO^-2 (indicating reductase activity) than many other ECM genera (Nygren et al., 2008). Despite the conserved nature of inorganic N genes across all species, we observed notable differences in ammonium utilization efficiency ( Figure 4A ). These functional differences suggest that variation in gene expression, post-translational modifications, or transporter kinetics may influence N acquisition and subsequent transfer to host trees, even when core genetic machinery appears conserved.

Unlike the highly conserved genes for inorganic N access, we identified substantial differences in the abundance of genes encoding both protein and chitin-degrading enzymes ( Figure 3 ). S. spraguei had the highest number of GH18 chitinases (28) compared to only 13 in S. weaverae and S. luteus. This variation is comparable to that observed in previous investigations of intergeneric variation in GH18 abundance across 12 species of ECM fungi which identified 6–22 genes encoding GH18 CAZYmes (Maillard et al., 2023). S. spraguei also possessed an exceptional expansion of aspartic proteases (312 in S. spraguei compared to 34 in S. ampliporus, 22 in S. clintonianus, and as low as 11 in S. weaverae ( Figure 3 ). These aspartic proteases, also known as acid proteases because they are most active under acidic conditions, play diverse roles in both host-interactions and the liberation of organic N, and have previously been shown to be the dominant proteolytic enzyme class enzyme produced by Suillus (Theron and Divol, 2014; Rineau et al., 2016). Additional analysis identified these aspartic proteases as polyporopepsins ( Supplementary Figure S1 , Supplementary Table S2 ), secreted enzymes implicated in the liberation of nutrients from complex organic substrates. In Paxillus involutus (also a member of the order Boletales), polyporopepsins are the most highly upregulated enzymes during organic N assimilation (Shah et al., 2013), and among the most highly upregulated enzymes during late-stage carbon starvation (Ellström et al., 2015).

The interspecific differences in N metabolism identified using genomics corresponded with functional performance in the laboratory bioassays, where all species grew well on ammonium but showed distinct patterns of biomass production on organic N sources that generally aligned with their genomic potential. On protein, S. spraguei produced the greatest amount of biomass, followed by S. clintonianus and S. ampliporus ( Figure 4B ). Most species showed poor or negative growth relative to controls on protein-tannin complexes ( Figure 4C ), suggesting limitations in their ability to access protein bound to phenolic compounds. These results are in agreement with previous work suggesting that phenolics inhibit the mobilization of N for many ECM fungi (Bending and Read, 1996). In contrast, most species maintained positive growth on protein absorbed on to minerals, though individual responses were highly species-specific ( Figure 4D ). These results both confirm the ability of certain ECM fungi to acquire organic N from mineral-associated proteins (Wang et al., 2021), and indicates that while protein bound to minerals may be more accessible than protein bound to phenolic compounds, the ability to access mineral-bound N varies considerably among species. On chitin S. spraguei again demonstrated the greatest biomass production with S. luteus, S. americanus, and S. clintonianus also displaying positive growth relative to controls, and S. ampliporus and S. weaverae displaying growth inhibition ( Figure 4E ). Notably, these trends were largely reflected in the abundance of chitin-active genes, where S. weaverae displayed the lowest number of chitin-active genes overall, and the worst overall growth when grown on chitin- a trend that was reversed for S. spraguei. These results may have been driven by the GH16 family, where S. ampliporus and S. weaverae displayed the lowest number of GH16 chitinases (19 each) in the genomic analyses, and S. spraguei the highest (with 24).

ECM fungi typically exhibit higher δ¹⁵N and lower δ¹³C than saprotrophic species, reflecting their use of ¹³C-depleted host photosynthate and ¹⁵N-enriched N sources. The ¹⁵N enrichment in ECM fungi may result from several processes: preferential transfer of ¹⁴N to host trees, acquisition of ¹⁵N-enriched N from deeper soil profiles, or sequestration of ¹⁴N in belowground mycelia that concentrates ¹⁵N in sporocarp proteins (Hobbie and Ouimette, 2009; Hobbie and Agerer, 2010; Hobbie et al., 2012). In this study, we found Suillus exhibited a range of δ¹⁵N values that fell within those previously reported for ECM fungi (Mayor et al., 2009). The higher concentrations of both total N and δ¹⁵N in cap tissue than in stipe tissue ( Supplementary Figures S2A, B ), matched expectations of tissue-specific protein concentration and isotopic fractionation (Hobbie et al., 2012). This trend was also reflected in the marginally positive relationship between differences in cap-stipe total N concentration and differences in cap-stipe δ¹⁵N ( Supplementary Figure S2C ), indicating that ¹⁵N enrichment corresponds with greater protein allocation to caps. However, the significantly negative relationship between overall sporocarp N status and differences in cap-stipe δ¹⁵N suggests that N-replete fungi exhibit reduced isotopic fractionation ( Supplementary Figure S2D ), possibly reflecting more efficient transport or reduced discrimination against δ¹⁵N when N availability is high.

Our isotopic analyses of historical specimens revealed complex N dynamics that challenge general predictions about terrestrial oligotrophication in coniferous forests. Contrary to expectations, total fungal N concentration increased significantly over time, with sporocarp N rising by approximately 21% from 1961 to 2019, while host plant N concentrations remained relatively stable and species-specific ( Figure 5A ). Strikingly, this finding contrasts sharply with a recent study from the same region that documented declining foliar N in broadleaf trees along with declining ECM sporocarp N over time (Michaud et al., 2024). Such contrasting N trajectories between broadleaf and conifer species have been documented elsewhere, including in Europe (Penuelas et al., 2020), suggesting fundamental differences in how these forest types respond to changing N availability. In particular, conifers may be more responsive to N deposition than broadleaf trees, potentially explaining this difference. Additionally, Suillus species exhibit ecological traits that distinguish them from the genera sampled in Michaud et al. (2024), including greater sensitivity to N deposition, high host-specificity, and unique mycelial morphology thought to be related to organic N acquisition (Lilleskov et al., 2011). The observed increase in Suillus sporocarp N despite proposed oligotrophication therefore suggests that these fungi and their coniferous hosts are responding differently to changing N regimes than broadleaf forest systems, potentially reflecting greater efficiency in organic N acquisition or both responses to and influences on the availability of different N pools across forest types.

Beyond these temporal trends in total N content, our isotopic analyses also provided detailed evidence of species-specific N acquisition strategies, with both host trees and fungi showing stable fractionation patterns over time but distinct species-level differences ( Figures 5B, C ). The carbon isotope data revealed divergent responses to changing atmospheric conditions, with host needles showing a steady decline in δ¹³C over time consistent with the Suess effect ( Figure 5D ), while fungal responses varied markedly by species. In line with the δ¹⁵N results, four species of Suillus (S. luteus, S. ampliporus, S. americanus, and S. clintonianus) exhibited significant declines in δ¹³C over time, mirroring the Suess effect and suggesting reliance on recently fixed carbon sources ( Figure 5E ). In contrast, S. spraguei and S. weaverae maintained stable δ¹³C values despite ongoing atmospheric δ¹³C decline, indicating greater utilization of older soil organic matter pools for both C and N acquisition. While the increased incorporation of δ¹³C likely reflects the use of organic N sources (which contain C-N bonds) as shown in Paxillus involutus (Akroume et al., 2019), it is possible that these patterns are influenced by the direct uptake of C, where stable δ¹³C in S. spraguei and S. weaverae could indicate the some degree of autonomous exploitation of carbon resources in soil organic matter reserves.

While S. weaverae and S. spraguei exhibited high δ¹⁵N values and stable δ¹³C values over time consistent with preferential access to older organic N pools, S. americanus displayed intermediate δ¹⁵N and a decrease in δ¹³C over time. Despite sharing the same host, these species appear to have evolved diverse N acquisition strategies indicating that differences in N use may have evolved in response to local ecological factors independent of host association. Soil N in ECM forests becomes progressively ¹⁵N-enriched with depth, and this functional diversity may reflect resource partitioning strategies that reduce interspecific competition, potentially enabling ECM fungi to specialize on different soil profiles and associated N pools, as demonstrated for Rhizopogon (Mujic et al., 2015). Such specialization on different N sources may have numerous downstream effects on metabolic function, as evidenced by S. americanus being more competitive than S. spraguei (Kennedy et al., 2020), potentially reflecting more efficient resource incorporation or tighter host associations. These results also highlight complex relationships between host and fungal N status that may interact synergistically. For example, while L. laricina maintained higher foliar N concentrations over time compared to P. strobus and P. resinosa, this pattern was reversed for δ¹⁵N ( Figures 5A-C ). These results may reflect fundamental differences in growth environment such as access to recently fixed or mineralized N sources in wetlands, or differences in the physiology, or fractionation patterns of host species. The low δ¹⁵N in the two L. laricina associates, S. ampliporus and S. clintonianus, may not only influence the low δ¹⁵N of their host (L. laricina needle δ¹⁵N was approximately 7.8‰ and 4.0‰ lower in δ¹⁵N than P. resinosa and P. strobus, respectively), but may also be influenced by the host, via the breakdown and assimilation of host litter which contains low δ¹⁵N. The mechanisms underlying these diverse N acquisition strategies and host-fungal interactions likely depend on soil chemical conditions such as pH, which may influence both organic and inorganic N accessibility.

For organic N, pH mediates the activity of acid-dependent enzymes like polyporopepsins, identified in high copy number in S. spraguei. Simultaneously, pH affects inorganic N access, as demonstrated by P. involutus where optimal ammonium uptake occurs between pH 4 and 5.5 (Javelle et al., 1999). Mechanistically, fungal ammonium assimilation via AMT transporters such as those identified in this study, rely on proton motive force through secondary active transport (Williamson et al., 2024). Although the relative contributions of proton gradient versus membrane potential remain unclear, acidic environments likely influence the efficiency of ammonium transport, particularly if proton motive force is primarily driven by proton gradient. Taken together, these patterns indicate that pH has the potential to shape both organic and inorganic N assimilation in ECM fungi (Fernandez and See, 2025), but the mechanisms and ecological consequences of these effects require further investigation.