Mycorrhizal Distributions Impact Global Patterns of Carbon and Nutrient Cycling

Most tree species predominantly associate with a single type of mycorrhizal fungi, which can differentially affect plant nutrient acquisition and biogeochemical cycling. Uncertainties in mycorrhizal distributions are non‐trivial, and current estimates disagree in up to 50% over 40% of the land area, including tropical forests. Remote sensing capabilities for mycorrhizal detection show promise for refining these estimates further. Here, we address for the first time the impact of mycorrhizal distributions on global carbon and nutrient cycling. Using the state‐of‐the‐art carbon‐nitrogen economics within the Community Land Model version 5, we found Net Primary Productivity (NPP) increased throughout the 21st century by 20%; however, as soil nitrogen has progressively become limiting, the costs to NPP for nitrogen acquisition—that is, to mycorrhizae—have increased at a faster rate by 60%. This suggests that nutrient acquisition will increasingly demand a higher portion of assimilated carbon to support the same productivity.

such mechanisms through Earth System models (ESMs), which allow comprehensive and spatially explicit assessment of the impacts of future climate on biogeochemical cycles in terrestrial ecosystems.
A large part of plant nitrogen and phosphorus is provided by fungal root symbionts (van der Heijden et al., 2015); thus, it is likely that mycorrhizal associations explain a large fraction of the variance in plant response to eCO 2 (Drake et al., 2011;Kivlin et al., 2013;Orwin et al., 2011;Sulman et al., 2017;Terrer et al., 2016Terrer et al., , 2018Terrer et al., , 2021. However, the global spatial distributions of these mechanisms as well as their potential impacts are still uncertain (Norby et al., 2017;Sulman et al., 2019). Only a handful of ESMs consider mycorrhizal nutrient acquisition when calculating carbon assimilation and allocation (Wang et al., 2010;Zaehle et al., 2015;Goll et al., 2017). The Community Land Model version 5 (CLM5) within the Community Earth System Model (CESM) currently enables an explicit representation of the functional differences between distinct types of plant symbiotic associations (Brzostek et al., 2014;J. B. Fisher et al., 2010;Lawrence et al., 2019;Shi et al., 2016). However, until recently, one of the major challenges in generating global estimates of nutrient limitation on the global carbon cycle is related to a lack of information of the spatial distribution of nutrient-acquiring plant-microbe symbioses. Despite the availability of regional maps of present and past plant symbiotic status (Brundrett, 2017;Jo et al., 2019;Menzel et al., 2016;Swaty et al., 2016), scientists have only recently begun to develop explicit global data about mycorrhizal and nitrogen fixing associations (Davies-Barnard et al., 2020).
Recently, scientists developed methods for extrapolating spatially sparse measurements into large-scale maps suitable for applications within ESMs (Shi et al., 2016;Soudzilovskaia et al., 2019;Steidinger et al., 2019;Sulman et al., 2019). These developments for the first time enable examination of how mycorrhizal spatial distributions are related to global carbon and nitrogen cycles. In this study, we seek a better understanding of mycorrhizae on global carbon and nitrogen cycles through incorporation of multiple state-of-the-art spatial distributions of mycorrhizal associations in a global land surface model. We first compare four global maps of mycorrhizal associations. Second, we perform transient global runs of CLM5 with climate change, nitrogen deposition, land use change, and eCO 2 through the 20th and 21st centuries in order to understand the impact of global changes associated with different spatially variable mycorrhizal representations. Finally, we evaluate possible feedback effects that spatial changes in mycorrhizal association due to climate change may have on the global carbon cycle in the future following a projected map presented in Steidinger et al. (2019).

Land Surface Model Description: CLM5
CLM5 includes the Fixation and Uptake of Nitrogen (FUN) module for calculating the carbon costs for each pathway of plant nitrogen uptake: symbiotic biological nitrogen fixation, non-mycorrhizal and mycorrhizal uptake of soil nitrogen, and nitrogen retranslocation from leaves (Allen et al., 2020;Brzostek et al., 2014;Fisher et al., 2010;Shi et al., 2016). Plants shift uptake pathways to minimize the carbon costs of nitrogen uptake. FUN simulates uptake from the two major types of fungi that plants associate with: arbuscular mycorrhizae (AM) and ectomycorrhizae (ECM) fungi.
To generate the trade-offs between AM, ECM, and non-mycorrhizal root uptake, FUN within CLM5 uses an estimate of the percentage of aboveground biomass per grid cell that associates with each mycorrhizal type (Brzostek et al., 2014;Shi et al., 2016). The carbon cost of nitrogen uptake (N cost ) from soil by mycorrhizae, for each soil layer j, is controlled by two uptake parameters (k n,pathway and k c,pathway ) that pertain respectively to the relationship between soil nitrogen and nitrogen uptake, and between fine root carbon biomass and nitrogen uptake, given as: ,pathway ,pathway cost,pathway, min, root, where k n,pathway (kgC m −2 ) and k c,pathway (kgC m −2 ) varies according to whether the pathway considered is referring to ECM or AM uptake. The uptake parameters were chosen so that the thresholds in nitrogen availability and root biomass mirror empirical shifts in the abundance of AM and ECM fungi and plants across fertility and latitudinal gradients as defined and calibrated by Brzostek et al. (2014). For ECM, k n = 0.05 kgC m −2 and k c = 0.3 kgC m −2 ; while for AM, k n = 0.1 kgC m −2 and k c = 0.05 kgC m −2 . N smin,j and c root,j are the soil nitrogen content (gN m −2 ) and fine root carbon biomass (gC m −2 ), respectively (see Supporting Information S1 for details).

Coupling Mycorrhizae Spatial Distribution Into CLM5
Plant functional types (PFTs) are used to classify plants according to their physical, phylogenetic, and phenological characteristics. The parameters associated with each PFT is determined or inferred from observable characteristics of the land surface following different land types. A spatial data product can be added as a 2D variable varying as function of latitude and longitude, but because land surface models usually work with the concept of PFTs, adding a third spatial dimension (i.e., PFT) can represent within grid cell heterogeneity and improve accuracy of land surface processes, as well as reduce model uncertainty (Braghiere et al., 2019). Here, given new spatial distributions of mycorrhizal associations based on observations at different spatial resolutions, we modified CLM5 and added mycorrhizal association types per PFT within each gridcell. The total carbon cost of nitrogen uptake (N cost ) from soil by mycorrhizae is calculated as: where perECM is the fraction of ECM association present in a PFT within a model grid cell.  Figure 1 and Supporting Information S1 for details about each map). Map A presents binary ECM values per PFT with most of the tropics having 100% AM and the boreal regions having 100% ECM. The other maps were derived from observations following different methodologies.

Simulation Protocols
First, for each ECM map, initial ecosystem carbon and nitrogen stocks for 1850 were generated using a spinup approach using 1850 concentrations of CO 2 (284.7 ppm) and the model's standard climate forcing dataset from the Global Soil Wetness Project Phase 3 version 1 (GSWP3v1) (Kim, 2017) at 1.9° × 2.5° spatial resolution. The Model for Scale Adaptive River Transport (MOSART) was turned on, and ice evolution on land was turned off. Model runs were performed with biogeochemistry mode on (without crops) for 200 years in "accelerated decomposition" mode (see Lawrence et al. (2019) for details) by cycling through the 1901-1920 climate forcing dataset, and for 400 years in regular spin-up mode until soil and plant carbon and nitrogen stocks achieved steady state. Historical simulations were performed from 1850 to 2010 using transient GSWP3 climate, nitrogen deposition, land-use change, and variable atmospheric CO 2 concentration.
Second, in order to illustrate the model sensitivity to changes in global spatial patterns of plant symbiosis due to climate change, we used a projected map of plant symbiotic status for 2070 using a relative concentration pathway (

Calculating Nitrogen Limitation
The risk of nitrogen limitation (NL) can be determined by evaluating if the growth rate of NPP used for nitrogen uptake with time is larger than the growth rate of total NPP with time. If the amount of NPP used for nitrogen uptake increases at a higher rate than the total NPP for a particular grid cell, that grid cell is considered to be at risk of spending too much carbon on nitrogen acquisition, and therefore, NL is closer to 1. On the contrary, if the amount of NPP used for nitrogen uptake increases at a lower rate than the total NPP for a particular grid cell, that area is not considered to be at risk of spending too much carbon on nitrogen acquisition. NL is calculated as: where α 1 is the slope of the linear regression of NPP used for Nitrogen uptake per gridcell (NPP_ NUPTAKE[lat,lon]) with time and α 2 is the slope of the linear regression of NPP (NPP[lat,lon]) plus NPP_ NUPTAKE(lat,lon) with time.

Different Estimates of Plant Symbiotic Status and Impacts on Nitrogen Uptake Pathways
To better visualize the differences from maps presented in Figure 1, the standard deviation of the average ECM fraction (%) across all four maps is shown in Figure 1e. A detailed description of each map can be found in the Supporting Information S1. Although all four maps agree in approximately 60% of the world area, some areas present large standard deviation values (>30%), for example, northern North America, northern and eastern Asia, as well as parts of tropical forests, that is, northwest Amazon, the central part of the Congo Basin, and parts of the maritime continent. These areas would benefit from more measurements of mycorrhizal association and further analysis. Remote sensing capabilities for mycorrhizal detection and hyperspectral Earth System modeling show promise for refining these estimates further (Braghiere et al., 2021;Fisher et al., 2012;Sousa et al., 2021). Map A (the default map in CLM5) presents the closest spatial patterns to map C, indicating an alignment of the assumptions that climate variables are the main drivers of global biogeography of plant-fungi symbiosis and the proposition that fixed values of mycorrhizal associations can be, to an extent, prescribed following PFTs spatial distributions as in map A.
Maps B and D suggest that map A overestimates ECM fraction in boreal regions, as well as drier areas of the world, such as the Atacama, Namibian, Somalian, Mongolian, Sonoran, and Australian deserts. All four maps present an ECM fraction between 30% and 60% in the eastern USA, while the standard deviation among maps over the same area is between 20% and 25%, significantly higher than the standard deviation present in African savannas and grasslands, India, and northeast Brazil (<10%). Comparatively, Europe and southeast Asia also present a standard deviation of 20%-25% in ECM fraction, while map A shows a lower fraction of ECM in Europe (20%-40%), followed by map B (30%-40%), map D (20%-70%), and map C (>90%). In southeast Asia, maps A, B, and D show an ECM fraction between 0% and 40%, while map C indicates higher values (40%-75%). Given that map A is based on PFT values, the biases in particular PFTs are presented in Figure S1.
The ECM-associated (NECM) and AM-associated (NAM) vegetation nitrogen uptake fluxes were the most impacted biogeochemical variables when including spatially explicit mycorrhizal status in CLM5, though the other nitrogen uptake pathways and their associated carbon costs were also impacted. Table S2 shows the average carbon cost per unit of nitrogen (gN kgC −1 ) in the period 2000-2010 for each different nitrogen uptake pathway as predicted by CLM5. On average for the period 2000-2010, the updated total carbon cost per unit of nitrogen globally according to maps B, C, and D increases 2.2% relative to the default map A in CLM5. The main areas where carbon costs of nitrogen uptake became higher are eastern North America, Europe, southeast Asia, and the tropics for mycorrhizal uptake (see Supporting Information S1). Changes in carbon costs of nitrogen acquisition via mycorrhizae uptake are 4.1% higher globally compared to map A because maps B, C, and D present higher fraction of ECM in the tropics and lower fraction of AM in the extra-tropics. While AM fungi act as scavengers for soil nutrients, generally occupying more fertile sites, ECM fungi can mobilize nitrogen from soil organic matter, favorable in sites with lower soil fertility, but also requiring greater carbon expenditure for nutrient acquisition (Chapman et al., 2006;Marschner & Dell, 1994;Read & Perez-Moreno, 2003). Under constant root biomass (0.1 kgC m −2 ), the ECM carbon cost of nitrogen uptake is two times higher than the AM carbon cost of nitrogen uptake. Under constant soil nitrogen (0.025 kgN m −2 ), nitrogen uptake efficiency switches from AM to ECM as fine root biomass increases (Brzostek et al., 2014). Carbon costs devoted to each nitrogen uptake pathway individually are presented in Supporting Information S1.

The Effect of Climate Change and eCO 2 on Nitrogen Limitation
Using the average of all four maps to determine the effect of climate change on NL, Figure 2 shows global NPP (PgC yr −1 ), global carbon cost of nitrogen uptake (NPP_NUPTAKE, PgC yr −1 ), global plant nitrogen demand (PLANT_NDEMAND, TgN yr −1 ), and the global nitrogen uptake (NUPTAKE, TgN yr −1 ). Nitrogen demand is calculated as the total nitrogen that would be required if all assimilated carbon was allocated according to idealized stoichiometric ratios. The CO 2 fertilization effect, with nitrogen deposition, and climate change increased photosynthetic rates across the globe, represented by an increase in NPP by 20% (40 Although the rates of nitrogen uptake systematically increase in response to a higher nitrogen demand by 25%, that is, NUPTAKE of 800 TgN yr −1 in 1850 to 1000 TgN yr −1 in 2010, the associated carbon cost of nitrogen acquisition increased at a faster rate, growing roughly 60% more expensive in 2010 (17.5 PgC yr −1 ) than it was in 1850 (11.2 PgC yr −1 ). In terms of the percentage of NPP spent in nitrogen acquisition, the values increased from 28% of NPP in 1850 to 33% of NPP in 2010. By 2075, it is projected that the NPP used for nitrogen acquisition will reach 35% of total NPP (∼22.5 PgC yr −1 ), suggesting ecosystems will have much less carbon available for allocation and plant growth, possibly becoming more susceptible to extreme events that require extra carbon for re-growth, such as droughts, fires, and insect outbreaks.
All transient runs from 1850 to 2010 with the new maps (B, C, and D) indicated a stronger effect of climate and eCO 2 on nitrogen limitation compared to map A. These findings highlight that as estimated by CLM5, not only has plant demand for nitrogen increased at a faster rate than actual nitrogen uptake, but also that the carbon costs associated with nitrogen acquisition have increased at a faster rate than the extra carbon gained through the CO 2 fertilization effect, that is, plants need to invest more carbon per unit of nitrogen uptaken. This pattern is projected to continue in the future, which means that it is unlikely current plant growth rates will be sustained globally. Figure 3a shows the risk of nitrogen limitation (NL) calculated as described in Equation 3. According to the transient runs from 1850 to 2010 using the default mycorrhizal map in CLM5 (map A), tropical forests have a medium to low risk of being further limited by nitrogen, which is in agreement with studies indicating that intact ancient tropical forests tend to accumulate and recycle large quantities of nitrogen relative to temperate forests (Hedin et al., 2009).
Savannas and forest-grassland transition zones in South America, Africa, and Australia, present a higher risk of NL to plant growth. Parts of the temperate forests in North America, Europe, and Asia, as well as boreal forests present a medium to high risk of nitrogen limitation.

The Feedback Impacts of Mycorrhizal Changes Due to Climate Change
Recent evidence suggests that anthropogenic influences, primarily nitrogen deposition and fire suppression, as well as climate change have increased AM dominance during the past three decades in the eastern United States (Jo et al., 2019). Globally, Steidinger et al. (2019) presented a study using the same environment-mycorrhizae relationships for current climate to project potential changes in the symbiotic status of forests in the future, suggesting that projected climate for 2070 will reduce the abundance of ECM trees by as much as 10%, with major changes in ECM abundance along the boreal-temperate transition zone ( Figure 3b).
Although the magnitude of the time lag between climate change and ecosystem responses is unknown, the predicted decline in ECM trees aligns with previous simulated warming experiments, which have demonstrated that some important ECM hosts decline at the boreal-temperate zones under future climate conditions (Reich et al., 2015). ECM fungi demonstrated increased responses of mycorrhizal fungal biomass under eCO 2 compared to AM fungi (Dong et al., 2018), as the simulated response in the tropics (Figure 3b).
While it has been previously reported that climate change should impact forest symbiosis, no study has ever evaluated the potential feedback of climate change effects on mycorrhizal distribution onto carbon and nitrogen cycles. The difference in NPP for the period of 2016-2075 between the simulations using the future maps of ECM fraction and the simulations using the present-day map C (Steidinger et al., 2019) are shown in Figure 3c.
Large parts of South America, especially areas associated with savannas, present the largest negative feedback effects on NPP due to future climate change impacts on mycorrhizal association, followed by areas with boreal forests. The impact over tropical forests and areas in China seems to benefit from a change in plant symbiotic status in the future. Although, these results should be interpreted carefully due to the limitation of the original forest plot training data in those areas of the globe used in Steidinger et al. (2019), machine learning algorithms indicate more ECM fungi in the tropics in the future, possibly due to the eCO 2 effect on the tropical climate.
In the SSP5-RCP8.5 runs from 2016 to 2075 with present-day plant symbiotic status, the growth rate of nitrogen uptake was 4.8 TgN yr −2 . In terms of carbon costs, NPP is projected to increase at a rate of 265.5 TgC yr −2 , while the carbon cost of nitrogen acquisition is projected to increase at a rate of 130.4 TgC yr −2 . The feedback effect of climate change on the spatial distribution of plant symbiotic status decreases NPP globally (from 58.3 to 58.2 PgC yr −1 ), a negative impact of −23.1 TgC yr −1 . The projected NPP increase rate with the future plant symbiotic status map is 266.2 TgC yr −2 , 0.7 TgC yr −2 faster than the projected NPP without changes in mycorrhizae associations. However, the carbon cost of nitrogen acquisition is projected to increase at a rate of 129.1 TgC yr −2 , versus 130.0 TgC yr −2 in the simulations without changes in the spatial distribution of plant symbiotic status. In terms of global NPP, these changes are predicted to increase carbon costs of nitrogen acquisition by 582.5 TgC.yr −1 , which significantly amplifies the effect of nutrient limitation on plant growth worldwide.  (Steidinger et al., 2019). The projected runs with CLM5 followed the SSP5 scenario in combination with RCP8.5 climate forcing from CESM2, member of CMIP6 simulations.

Conclusions
To overcome the lack of global spatial representations of mycorrhizal associations, recent studies have assembled high-resolution digital maps of the global distribution of biomass fractions of different types of mycorrhizae associations.
In our analysis, we show that differences between these data products have significant impacts on the nitrogen and carbon cycles in CLM5. Nonetheless, this comparison did not aim to determine which map is the most accurate. Rather, we assessed the impact of different mycorrhizal representations in CLM5 to determine signs of changes in the global nitrogen and carbon cycles. We found a negative impact on future NPP due to feedback effects of climate change and eCO 2 on mycorrhizae spatial distribution.
Although the model runs with different spatial representations of plant symbiotic status differ substantially in terms of total nitrogen acquisition and their relative carbon costs, all experiments using the observation based maps do agree that the increasing rate of plant nitrogen demand is higher than the rate of nitrogen uptake as previously reported. Moreover, our simulations found that the carbon costs of nitrogen acquisition also increase at a higher rate than NPP itself, indicating that plants must invest more carbon per unit of nitrogen uptake to sustain growth at current rates globally.
This is the first study using multiple observation-derived global maps of mycorrhizal association within an ESM to estimate the impacts of climate change on mycorrhizae and their feedback on global carbon and nitrogen cycles. Future model developments include better constraining cost parameter functions with model-observation fusion frameworks, as well as the addition of other nutrients and mycorrhizal types (such as ericoid mycorrhizae).

Data Availability Statement
A patch file with the modified version of CLM5 and all python scripts used for analyses and plots are available in https://doi.org/10.6084/m9.figshare.12919385.v1.