ARTICLE https://doi.org/10.1038/s41467-019-13019-2 OPEN Global mycorrhizal plant distribution linked to terrestrial carbon stocks Nadejda A. Soudzilovskaia 1*, Peter M. van Bodegom 1,10, César Terrer 2,3,10, Maarten van’t Zelfde1, Ian McCallum4,10, M. Luke McCormack 5, Joshua B. Fisher6,7, Mark C. Brundrett 8, Nuno César de Sá 1 & Leho Tedersoo9,10 1234567890():,; Vegetation impacts on ecosystem functioning are mediated by mycorrhizas, plant–fungal associations formed by most plant species. Ecosystems dominated by distinct mycorrhizal types differ strongly in their biogeochemistry. Quantitative analyses of mycorrhizal impacts on ecosystem functioning are hindered by the scarcity of information on mycorrhizal dis- tributions. Here we present global, high-resolution maps of vegetation biomass distribution by dominant mycorrhizal associations. Arbuscular, ectomycorrhizal, and ericoid mycorrhizal vegetation store, respectively, 241 ± 15, 100 ± 17, and 7 ± 1.8 GT carbon in aboveground biomass, whereas non-mycorrhizal vegetation stores 29 ± 5.5 GT carbon. Soil carbon stocks in both topsoil and subsoil are positively related to the community-level biomass fraction of ectomycorrhizal plants, though the strength of this relationship varies across biomes. We show that human-induced transformations of Earth’s ecosystems have reduced ectomycor- rhizal vegetation, with potential ramifications to terrestrial carbon stocks. Our work provides a benchmark for spatially explicit and globally quantitative assessments of mycorrhizal impacts on ecosystem functioning and biogeochemical cycling. 1 Institute of Environmental Sciences, Leiden University, 2333 CC Leiden, the Netherlands. 2 Institut de Ciència i Tecnologia Ambientals (ICTA) Universitat Autonoma de Barcelona, Barcelona, Spain. 3 Physical and Life Sciences Directorate, Lawrence Livermore National Laboratory, Livermore, CA 94550, USA. 4 Ecosystems Services and Management Program, International Institute for Applied Systems Analysis, Schlossplatz 1, A-2361 Laxenburg, Austria. 5 Center for Tree Science, The Morton Arboretum, 4100 Illinois Route 53, Lisle, IL 60532, USA. 6 Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Dr., Pasadena, CA 91109, USA. 7 Joint Institute for Regional Earth System Science and Engineering, University of California, Los Angeles, CA, USA. 8 School of Biological Sciences, Faculty of Science, University of Western Australia, Crawley 6009 WA, Australia. 9 Natural History Museum and Institute of Ecology and Earth Sciences, University of Tartu, 14a Ravila, 50411 Tartu, Estonia. 10These authors contributed equally: Peter M. van Bodegom, César Terrer, Ian McCallum, Leho Tedersoo *email: n.a.soudzilovskaia@cml.leidenuniv.nl NATURE COMMUNICATIONS | (2019)10:5077 | https://doi.org/10.1038/s41467-019-13019-2 | www.nature.com/naturecommunications 1 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-13019-2 M ycorrhizas are mutualistic relationships between plants an ecosystem and soil carbon content in topsoil (0–20 cm), and fungi, in which fungi supply plants with nutrients medium (20–60) and deep (60–100 cm) subsoil layers. and plants provide carbon to fungi1. Among mycor- rhizal types, arbuscular mycorrhiza (AM), ectomycorrhiza (EcM) and ericoid mycorrhiza (ErM) are geographically the most Results and Discussion widespread, colonizing over 85% of vascular plants across vege- Assembly of mycorrhizal vegetation maps. To generate global tated terrestrial biomes1–4. Due to the facilitation of plant maps of mycorrhizal vegetation, we estimated biomass fractions nutrient acquisition1 and the large biomass of fungal networks in of AM, EcM, ErM and NM plants within each combination of soil5, the presence and type of mycorrhiza are among the key continent × ecoregion × land cover type. Supplementary Fig. 1 determinants of ecosystem functioning6–9 and biogeochemical illustrates the data assembly processes for the maps. Ecoregions cycling10–13. Thus, the types of mycorrhizal associations present follow Bailey32 (Supplementary Data 1), and land cover types likely also affect the global distribution of soil carbon stocks. were retrieved from the ESA CCI land cover map33, which spe- There is growing evidence that ecosystems dominated by EcM cifies cover and biomass fractions of trees, shrubs and herbaceous and ErM vegetation exhibit higher topsoil carbon to nitrogen plants (Supplementary Data 2). For each combination, we ratios (C/N) compared with ecosystems dominated by AM determined the dominant species or group of species from 1568 plants11,12,14,15, although in temperate forests the pattern may be vegetation surveys (Supplementary Data 3). For these species, we reversed in deeper soil layers16. The mechanisms driving these determined mycorrhizal type using the FungalRoot database differences are heavily debated in current literature, with distinct v1.034 (see Supplementary Data 4 for data sources). Integrating physiological traits of mycorrhizal fungi most likely playing a these data, we obtained mycorrhizal plant biomass fractions of critical role16–20. AM, EcM, ErM and NM plants for each combination of Bailey Although it can be argued that high abundance of EcM plants ecoregion, continent, and land cover type (Supplementary Data 5 is a consequence rather than a driver of high soil C stocks, a large and 6). These fractions were overlain on a global grid. body of recent findings provides evidence that EcM symbionts Our maps (Fig. 1) provide quantitative estimates of the may be the key drivers of topsoil carbon accumulation through distribution of aboveground biomass fractions among AM, EcM two interacting mechanisms. First, EcM fungi produce greater and ErM plants within areal units of 10 arcmin. The use of a biomass of more recalcitrant mycelium compared to AM fungi5. detailed map of ecoregions32 provides much greater resolution Second, while EcM fungi are more efficient in taking up N in N- compared with the biome-based patterns of mycorrhizal poor soils than AM fungi or roots10,21, EcM fungi immobilize distributions reported by Read3 > 25 years ago, whereas the land most of the N in their own biomass. This suppresses saprotrophic cover map33 enabled us to provide accurate spatial positioning of decomposition process22 and reinforces the competitive advan- ecosystem boundaries based on satellite-derived data, explicitly tage of EcM and ErM plants via enhanced N limitation22,23. taking into account human-driven transformations of vegetation. A full understanding of global carbon and nitrogen stocks We validated the map data using four independent datasets: (i) requires quantitative models on the distribution of mycorrhizal forest biomass structure for Eurasia35, (ii) a global dataset of forest types in ecosystems18. Despite the existence of regional maps of biomass structure used for an analysis of mycorrhizal impacts on current24,25 and past26 mycorrhizal vegetation, and on the dis- carbon vs nitrogen dynamics19, (iii) estimates of mycorrhizal tribution of mycorrhizal fungal species27,28 we still lack global associations in the USA based on satellite remote sensing36, and information on the distribution of biomass of mycorrhizal plants, (iv) West Australian map of mycorrhizal root abundance24 which is a much better proxy for mycorrhizal impacts on eco- (Supplementary Fig. 2). This validation revealed that the vast system functioning than the biodiversity of mycorrhizal sym- majority of the data (87% of the AM data points and 89% of the bionts. While current terrestrial biosphere models simulate EcM data points) deviate by < 25% from the measurements19,35,36, feedbacks between the carbon cycle and vegetation distribution29, when excluding ESA land use classes comprising poorly resolved most models ignore mycorrhizal types and their effects on combinations (i.e. mixed classes of land cover types33, such as nutrient cycling. Integration of such information is expected to “Tree cover, broadleaved, evergreen, closed to open (>15%)”, see provide a more realistic simulation of carbon and nutrient fluxes Methods for details) that were difficult to couple to our associated with plant nutrition17,18,21 and soil carbon cycles15,21. classification scheme. The relationship between the validation Quantitative models of the distribution of mycorrhizal vegetation data and our estimates is shown in Supplementary Fig. 3. constitute an important missing link between the known effects of Our maps of mycorrhizal vegetation were assembled based on mycorrhizas in biogeochemical cycles and their global impacts30. multiple published datasets, using a number of conversion factors Human activities such as forest logging, urbanization and to obtain per pixel values of mycorrhizal plants biomass fractions. agricultural practices have altered 50–75% of the Earth’s terres- These conversions as well as the fact that the plant species trial ecosystems31, transforming areas with previously natural distribution data (Supplementary Data 4) originates from multi- EcM and ErM vegetation into AM and non-mycorrhizal (NM) ple sources constitute important uncertainty sources in our vegetation. However, the impact of anthropogenic land use shifts dataset. We examined the uncertainty of our maps based on on biogeochemical cycles associated with mycorrhiza have uncertainties of tree, shrubs and herbaceous plant fractions remained poorly known due to the lack of appropriate spatial within the land cover types37, and the number of data sources information. used to assess mycorrhizal fractions of plant biomass within each Based on a comprehensive quantitative evaluation of plant- combination of Bailey ecoregion × continent; see Methods for mycorrhizal associations and the distribution of vascular plant details. Supplementary Fig. 4 shows spatial distribution of species across biomes and continents, we assembled high- uncertainties. The mean uncertainties of AM, EcM, ER and resolution digital maps of the global distribution of biomass NM maps are 19.6, 17.6, 14.6 and 15.0% at the 90% confidence fractions of AM, EcM, ErM and NM plants. Building on these interval. Overall, tropical areas have the highest uncertainties of maps, we assessed: (i) the amount of aboveground biomass car- the mycorrhizal fraction data, reaching 50% (AM) in the Amazon bon currently stored in each type of mycorrhizal vegetation; (ii) region. Therefore, our maps should be used with caution for these the impact of conversion of natural ecosystems to croplands on areas. Future sampling efforts of mycorrhizal vegetation distribu- the distribution of mycorrhizal types globally; and (iii) the rela- tion should be more focussed on tropical areas of Asia, Africa and tionships between relative abundance of AM and EcM plants in South America. 2 NATURE COMMUNICATIONS | (2019)10:5077 | https://doi.org/10.1038/s41467-019-13019-2 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-13019-2 ARTICLE a b c d 0 25 50 75 100 Fig. 1 Percentage of aboveground plant biomass of mycorrhizal vegetation. a Arbuscular mycorrhizal plants, b ectomycorrhizal plants, c ericoid mycorrhizal plants, and d non-mycorrhizal plants. The map resolution is 10 arcmin. See Supplementary Fig. 4 for associated uncertainty values. Source data are provided as a Source Data file Mycorrhizal vegetation and aboveground carbon stocks. By fraction of AM, EcM, ErM and NM plants in a cropland-free linking our maps of mycorrhizal vegetation to satellite observa- world (Supplementary Figs. 6 and 7). The current biomass frac- tions of global aboveground biomass carbon38, we estimated the tions of AM plants have increased in Europe, parts of Asia and amount of aboveground biomass carbon stored in arbuscular, North America, but declined in Africa, Asia (mostly India) and ecto-, ericoid and non-mycorrhizal vegetation as 241 ± 15, 100 ± South America, coinciding with increase in non-mycorrhizal 17, 7 ± 1.8 and 29 ± 5.5 GT (mean values ± uncertainty at 90% vegetation (Fig. 3, Supplementary Fig. 8). Our analysis suggests confidence interval; Fig. 2). In this analysis, the data were scaled that EcM biomass has declined in all continents, primarily due to to a resolution of 15 arcmin to match biomass estimates38. Most a replacement of natural forests by agricultural lands, whereas of the aboveground carbon stock stored in arbuscular mycorrhizal ErM biomass has remained unchanged. vegetation is situated in tropical forests (Fig. 2a). Supplementary Table 1 shows per-biome distribution of the carbon stocks among mycorrhizal types. Biomass fractions of mycorrhizal types and soil C stocks. Recent field research in the US temperate forests suggests that soil carbon content increases with increasing EcM abundance in Impacts of land transformations on mycorrhizal vegetation. topsoil layers; but, depending on forest type, this relationship may Agricultural practices drive the replacement of natural vegetation be reversed in deeper soil16,42. This is in agreement with the by facultatively AM crops1,39, which could also be de facto non- Microbial Efficiency-Matrix Stabilization hypothesis, which pre- mycorrhizal due to destruction of hyphal networks by ploughing dicts that ecosystems with rapid decomposition, such as most and excess fertilisation40,41. Using past vegetation estimates, AM-dominated forests19, enhance soil organic matter (SOM) Swaty et al.26 showed that across conterminous USA, agriculture stabilization by accelerating the production and deposition of has reduced the relative abundance of ectomycorrhizal plants microbial residues43–45. compared with other mycorrhizal types. However, global quan- To separate the effects of biome and mycorrhizal type on soil C tifications of agricultural impacts on distribution of mycorrhizas on a global scale, we modelled the relationships between soil have not been possible until now. Based on the current land use carbon content, biome type46, and biomass fractions of AM and data underlying our maps (Supplementary Data 4), we assessed EcM plants. We did not analyse the relationship between ErM mycorrhizal distributions on Earth in the absence of croplands. plant cover and soil C content, due to a small proportion of the For each ecoregion-continent-land cover combination that con- ErM plant biomass in the majority of ecosystems. tained croplands, we replaced current biomass fractions by esti- We conducted separate analyses for the topsoil (uppermost mates of per-grid cell biomass fractions in AM, EcM, ErM and 20 cm soil layer) and subsoil (20–60 and 60–100 cm soil layers), NM plants that would be expected at these locations based on as obtained from the ISRIC-WISE Soil Property Databases, at a natural vegetation types (see Methods for details, and Supple- resolution of 30 arcsec47. The data sources used for these analyses mentary Data 7–8 for data). Based on these data, we generated are independent: the ecoregion classification, and hence mycor- maps presenting potential natural distributions of biomass rhizal type distribution, does not account for edaphic parameters, NATURE COMMUNICATIONS | (2019)10:5077 | https://doi.org/10.1038/s41467-019-13019-2 | www.nature.com/naturecommunications 3 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-13019-2 a b c d 0 3 6 9 12 Fig. 2 Amount of carbon stored in plant biomass in vegetation of different mycorrhizal types (Mt C per-grid cell of 15 arcmin). a Arbuscular mycorrhizal plants, b ectomycorrhizal plants, c ericoid mycorrhizal plants, d non-mycorrhizal plants. The amount of aboveground biomass carbon stored in arbuscular, ecto-, ericoid and non-mycorrhizal vegetation is 241 ± 15, 100 ± 17, 7 ± 1.8 and 29 ± 5.5 GT (mean values ± uncertainty at 90% confidence interval), respectively whereas soil C data are unrelated to that of vegetation. Model Our analysis of the relationships between mycorrhizal type and comparisons were based on the Akaike information criterion soil carbon is based on ancillary maps, which feature large (AIC). The relative importance of each predictor was examined uncertainties. Analyses of relationships between ISRIC-WISE using the Lindemann-Merenda-Gold (LMG) metric, providing predicted soil carbon and the original data that were used to the fraction of variance explained by each predictor, within the generate the ISRIC soil map yield R2-values in the range of total variance explained by the model. 0.4–0.649,50. This uncertainty adds ambiguity to our analysis, and Our global analysis revealed a moderately strong positive reduces the reliability of quantitative estimates of the relation- relationship between the biomass fractions of EcM plants and ships between EcM plant biomass fractions and soil C. However, both topsoil and subsoil carbon (Fig. 4, Supplementary Figs. 8 and these uncertainties equally apply to the analysis of AM vs soil C as 9). Consistent with the current paradigm45, cf. ref. 11, our analysis well as to that of the EcM vs soil C, due to the fact that the revealed that biome is the main predictor of soil carbon stocks. analysis is based on the same geographical data points. Therefore, Even so, the grid cell biomass fraction of EcM plants still we consider that the high uncertainty of the ISRIC-WISE soil data accounted for one third of the explained variation in both topsoil is unlikely to affect the qualitative nature of our conclusion that carbon and in subsoil (Table 1). The interaction between biome AM and EcM vegetation differently relate to soil C. and EcM biomass fractions was significant (P < 0.001) but only marginally important (LMG = 1%), suggesting that the increase Implications. Our work provides the quantitative estimates of the in topsoil carbon along with an increase in EcM plant biomass is global biomass distribution of arbuscular, ecto-, ericoid and non- mostly independent from the environment (Table 1, Fig. 4, mycorrhizal plants, accounting for human-induced transforma- Supplementary Figs. 8, 9, Supplementary Table 2). The total tion of habitats. Previous research has shed light onto the dis- aboveground carbon stock stored in the EcM plants was a tribution patterns of mycorrhizal plant and fungal species relatively worse predictor of soil carbon stocks compared with the richness27,28, and onto low-resolution (1 arcdegree) distribution EcM relative biomass fraction, showing in all cases higher AIC patterns of mycorrhizal trees51. In contrast, our maps directly and lower R2. In contrast to EcM, AM biomass fractions per-grid reflect the global distribution of biomass fractions of mycorrhizal cell showed negative but inconsistent relationships to soil carbon, plants across all biomes and all main vegetation types. Availability with a contrasting positive trend in tundra (Fig. 4, Supplementary of such data at high resolution of 10 arcmin provides an oppor- Fig. 10). The latter relationship explained only a small amount of tunity for multiple potential analyses aimed at unravelling variance in tundra soil carbon stocks (12% in top 0–20 cm, 1.7% mycorrhizal impacts on ecosystem functioning at large- in 20–60 cm layer and 0.2% in 60–100 cm layer), and could arise geographical scales. from less accurate data of soil C in tundra47 due to local Our maps were derived at spatial resolutions allowing landscape heterogeneity and prevalence of facultative AM plants identification of the global patterns of mycorrhizal distributions with no or low mycorrhizal colonization48. and are most appropriate for global and large-geographical scale 4 NATURE COMMUNICATIONS | (2019)10:5077 | https://doi.org/10.1038/s41467-019-13019-2 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-13019-2 ARTICLE a b c d −100 0 100 Fig. 3 Changes in biomass fractions of mycorrhizal vegetation induced by crop cultivation and pastures. a Arbuscular mycorrhizal plants, b ectomycorrhizal plants, c ericoid mycorrhizal plants, d non-mycorrhizal plants. Purple colours indicate losses, green colours indicate gains. Uncertainties are shown in Supplementary Fig. 7. Source data are provided as a Source Data file analyses of mycorrhizal impacts on ecosystem functioning and topsoil and subsoil, while AM vegetation has more variable, global drivers thereof. Recent estimates suggest that the total soil weaker and mostly negative relationships. This analysis does not carbon loss due to agricultural practices accounts for 133 GT, provide evidence of causality of this relationship as multiple with great acceleration of losses during the past 200 years52. Our environmental variables such as climate, soil nutrients, especially analyses accounting for pre-agricultural patterns in EcM plant nitrogen availability, and soil texture may affect both soil carbon distribution point to large-scale losses of ectomycorrhizal and mycorrhizal plant distributions. Nonetheless, our study vegetation, with potentially strong effect on the amount of C establishes a quantitative framework to test the relation between stored in soils. Analyses of agricultural impacts presented in this the dominance of mycorrhizal types and soil C stocks. Complete paper are based on the assumption that these impacts are limited and directional understanding of the complex nature of the to shifts in plant species composition, and do not encompass hierarchy of environmental drivers controlling soil carbon shifts in soil water and nutrient availability, which could affect patterns requires further detailed (experimental) investigations activity of mycorrhizal fungi, depth distribution of mycorrhizas in of the hierarchy of different predictors and importance of local the soil, and shifts among mycorrhizal fungal species composition edaphic variables. due to the introduction of exotic species. While such simplifica- Our estimates of the carbon stocks in AM and EcM above- tions are necessary for the analyses reported in this paper, they ground biomass together with the quantitative relationships should be considered when interpreting our results. Furthermore, between soil carbon stocks and AM or EcM plant dominance in our analyses do not address other human impacts that can lead to ecosystems provide qualitative insights into the global carbon shifts among AM and EcM vegetation, such as climate change, cycle, highlighting the substantial role of mycorrhizas therein. introduction of invasive species and nitrogen deposition. The Tropical forests, usually dominated by AM symbiosis (Fig. 2), latter is known to be an especially important driver of contain 162 GT (44%) of global aboveground biomass38, whereas mycorrhizal vegetation shifts42,53, as it negatively affects abun- the predominately EcM temperate and boreal forests altogether dance of ectomycorrhizal fungi in soil54. Given that nitrogen store only 21% of global aboveground biomass carbon38, deposition leads to replacement of ectomycorrhizal plants by indicating that contribution of EcM vegetation to the above- arbuscular and non-mycorrhizal vegetation, it further enhances ground biomass carbon is relatively small. In contrast, below- agricultural impacts on soil carbon losses. ground carbon stocks are positively correlated to the proportion The question whether increased domination of EcM plants in of EcM plant biomass, suggesting that mycorrhizal contribution an ecosystem is associated with higher soil carbon content across to large carbon stocks in these regions occurs primarily through both top- and subsoil is heavily debated11,12,16,18–20,30. However the carbon supply to belowground organs and mycorrhizal fungi, previous studies have been based on a limited number of which is further emphasized with slowed decomposition observations11 or regional-scale analyses11,12,16,53,55. Our analysis processes22. Furthermore, our analyses revealed a relatively shows that across large geographical scales, higher cover of EcM stronger relationship between soil C and EcM plant biomass vegetation is broadly associated with greater soil C stocks in both fraction (%) than between soil C and the amount of carbon stored NATURE COMMUNICATIONS | (2019)10:5077 | https://doi.org/10.1038/s41467-019-13019-2 | www.nature.com/naturecommunications 5 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-13019-2 a ecosystem functioning, without introducing a circular reasoning 40 caused by the use of common environmental variables. Inclusion of mycorrhizal distribution into vegetation models would provide a benchmark for testing hypotheses about mycorrhizal impacts Topsoil C content (kg/m2) 30 on ecosystem functioning and related ecosystem services. Our maps enable quantifying relationships between mycorrhizal abundances in ecosystems as well as soil and vegetation carbon 20 content in global-scale analyses of biogeochemical cycles. In particular, the results of our study suggest that restoration of native vegetation especially in abandoned agricultural and barren 10 land may help alleviate anthropogenic soil carbon losses and ameliorate increases in atmospheric greenhouse gases. 0 0 25 50 75 Methods Biomass fraction of EcM plants Data used for construction of mycorrhizal vegetation maps. To construct the maps of mycorrhizal biomass fraction distribution, we integrated (1) data con- b cerning dominant plant species and their growth forms within each continent × 40 ecoregions × land cover type combination, (2) data describing mycorrhizal type of these species, and (3) data estimating cover and biomass fractions of trees, shrubs and herbaceous plants within individual land cover types. Here we defined a Topsoil C content (kg/m2) 30 species or a set of species with a similar mycorrhizal strategy and growth form as ‘dominant’ if it constituted > 20% of vegetation biomass. Supplementary Fig. 1 shows a flowchart of the map assembly process. 20 We selected the global ecoregion map of Bailey32 with 98 ecoregions (Supplementary Data 1), provided by the Oak Ridge National Laboratory Distributed Active Archive Center57 (spatial resolution 10 arcmin), as a basis for mapping of global-scale distribution of mycorrhizal types. This map was preferred 10 over that of biomes provided by Olson and co-workers46 and others because of higher level of detail for ecoregions and because the boundaries of ecoregions were more strongly related to the distribution of mycorrhizal types. Ecoregions spanning 0 across multiple continents were considered separately for each continent. We used 0 25 50 75 100 the continent division based upon the FAO Global Administrative Unit Layers (http://www.fao.org/geonetwork/srv/en/). Biomass fraction of AM plants We determined the land cover types present in each ecoregion using a satellite observation-based map for the year 2015 generated by the European Space Biome Agency33. This map includes 38 land cover categories such as croplands, urban areas, grasslands and forests of various types, with a spatial resolution of 300 m 1Tropical broadleaf forests (Supplementary Data 2). 2Subtropical broadleaf forests Based on vegetation surveys assigned to ecoregions (1568 data sources; 3Temperate broadleaf and mixed forests Supplementary Data 3) we determined the dominant plant species and their growth form for each continent × ecoregion × land cover type combination. We 4Coniferous forests used relative abundances of plant species averaged across different data points, as 5Boreal forests/taiga available across vegetation surveys, with equal weight to all observations. 6Montane grasslands and shrublands We assigned mycorrhizal type (AM, EcM, ErM, NM) to each dominant species using the FungalRoot database34. Species with dual EcM-AM and AM-NM 7Grasslands categories were allocated to both types equally (50% weight). Mycorrhizal status of 8Tundra species with no empirical records were extrapolated from congeneric and confamilial species. Therefore, all Diapensiaceae and Ericaceae species were Fig. 4 Quantitative relationships between topsoil (0–20 cm) C and biomass considered ErM58, except for Enkianthus (AM), Arbuteae, Pyroleae, Monotropeae fraction of mycorrhizal vegetation in natural ecosystems. a EcM plants and and Pterosporeae (all subtypes of EcM). Because of multiple incorrect reports and b AM plants. The outcomes of individual models are presented in the alternative definitions for EcM, we took a conservative approach by considering plants to be EcM only when this was supported by multiple independent studies Supplementary Table 2. Croplands were excluded from the analysis. Per- and the proportion of conflicting reports was < 50%59. Although most crop plants biome predictions are shown in different colours. Source data are provided are able to form arbuscular mycorrhizas, intensive agricultural practices and as a Source Data file breeding may lead to reduction or loss of mycorrhizal infection40,41. Therefore, rain fed and flooded croplands were considered to feature partly AM and partly NM vegetation (Supplementary Data 5, 7), unless data indicating dominance of NM vegetation was available. Crop species that belong to Brassicaceae family were in EcM plants. These findings suggest that belowground carbon considered NM. allocation by plants through mycorrhizal pathways is not directly Combining the prevailing dominant plant species, their growth form and proportional to the aboveground plant biomass, supporting the mycorrhizal type, we estimated the biomass proportions of EcM, AM, ErM and view of the importance of mutualism-parasitism trade-off in NM vegetation in each ecoregion by continent by land cover combination ectomycorrhizal associations56 across biomes. (Supplementary Data 5 and 6). We considered that in forests with a sparse understorey, trees contribute 90–95% of the biomass and the understorey accounts Taken together, this study provides a benchmark for relating for 5–10% of biomass36,60–63. In forests with a dense layer of shrubs, trees, shrubs ecosystem processes to the functioning of distinct types of and herbs/dwarf shrubs contribute 70 ± 15%, 20 ± 10% and 10 ± 5% to plant mycorrhizas on a global scale. So far, quantitative global biomass36,60–63, respectively. In shrublands, we consider shrubs and herbs to information about mycorrhizal distribution was virtually absent account for 90 ± 5% and 10 ± 5% of biomass, respectively35,60–63. We considered that woodlands harbour 30 ± 10% of biomass in trees, 30 ± 10% of biomass in despite the high demand for such data12,13,30. In spite of some shrubs, and the remaining biomass in herbaceous vegetation35,60–63. This resulted uncertainty, our mycorrhizal distribution maps provide an in biomass proportions of each mycorrhizal type in continent × ecoregion × land essential source for systematic analyses of mycorrhizal biogeo- cover type (Supplementary Data 5, 6). graphy and environmental drivers. Because our maps are based As we focused on the biomass of mycorrhizal plants and not on species diversity; we did not attempt to map the distribution of orchid mycorrhiza. Orchid on field data, and not on a machine-learning model trained with species are never abundant in ecosystems in terms of biomass, and they are environmental variables, they provide independent data for therefore unlikely to play an essential role in biogeochemical cycles at large regional examining the relationships between mycorrhizal status and scales. 6 NATURE COMMUNICATIONS | (2019)10:5077 | https://doi.org/10.1038/s41467-019-13019-2 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-13019-2 ARTICLE Table 1 Summary of generalized linear models (glm) predicting soil carbon stocks Predicted variable Model R2 Predictor P-value LMG (%) Topsoil C 0–20 cm EcM + Biome + EcM × Biome 0.53 EcM <0.001 42 Biome <0.001 57 EcM × Biome <0.001 1 Subsoil C 20–60 cm EcM + Biome + EcM × Biome 0.38 EcM <0.001 39 Biome <0.001 60 EcM × Biome <0.001 1 Subsoil C 60–100 cm EcM + Biome + EcM × Biome 0.33 EcM <0.001 35 Biome <0.001 64 EcM × Biome <0.001 1 Topsoil C 0–20 cm AM + Biome + AM × Biome 0.54 AM <0.001 38 Biome <0.001 56 AM × Biome <0.001 6 Subsoil C 20–60 cm AM + Biome + AM × Biome 0.33 AM <0.001 29 Biome <0.001 67 AM × Biome <0.001 2 Subsoil C 60–100 cm AM + Biome + AM × Biome 0.32 AM <0.001 31 Biome <0.001 67 AM × Biome <0.001 2 Predictions are made for C at 0–20, 20–60, and 60–100 cm depth and are based on biome and fraction of EcM or AM plants in vegetation biomass. R2—Cragg and Uhler’s pseudo R2. LMG—relative importance of individual predictors in a model examined through the Lindemann-Merenda-Gold metric. The LMG shows the percentage of variance explained by each of model predictors within the entire variance explained. The P values show the outcome of ANOVA type I models (n = 78883 in all models). Source data are provided as a Source Data file Assembly of raster maps of mycorrhizal vegetation. We generated raster maps The dataset of Fisher et al.36 provides the relative cover of AM and EcM plants based on the proportional mycorrhizal type biomass data. We overlaid the raster from Landsat scenes centred on four sites in USA: Lilly-Dickey Woods (Indiana), map of Bailey ecoregions (10 arcmin resolution)57 with the raster of ESA CCI land long-term research site of Smithsonian Conservation Biology Institute (Virginia), cover data (300 m resolution)33, which we converted to 10 arcmin using a nearest Tyson Research Center Plot (Missouri), and a long-term research site of Wabikon neighbour approach. The resulting raster was overlain with the polygon map of Forest Dynamics (Wisconsin). Given that the dataset comprises forested areas only, continents, rasterized at 10 arcmin. To each pixel, we assigned the corresponding we considered areal coverage of AM and EcM plants provides a good estimate of mycorrhizal type proportions, considering the prevailing combination of Bailey AM and EcM plant biomass. Using this dataset, we directly compared the AM and ecoregion × land cover in each continent. Because some ecoregions covered mul- EcM coverage per pixel with the data of our maps. tiple isolated parts of continents and proportions of mycorrhizal type distribution In the datasets19,35,36 biomass and areal fractions of AM and EcM plants are differed in these regional areas by > 2-fold, we split these ecoregions into two or always considered to sum up to 100%. As these datasets do not provide more subregions using the ArcGIS 10.2.2 Raster Calculator. information about non-mycorrhizal and ericoid plants, we estimated these as 5–10% and 0–20%, respectively, depending on the dominant tree association, and accordingly reduced the values of AM and EcM biomass fractions in the validation calculations. While we considered this approach to be acceptable for the validation Impact of croplands. We estimated the effect of agriculture on mycorrhizal type of AM and EcM data, the data quality is not high enough to validate the NM and distribution by substituting the proportions of mycorrhizal types in croplands ErM maps. (land-cover types 10, 11, 12, 20 and 30; see Supplementary Data 2) with estimates As these three datasets19,35,36 represent forest data, we evaluated whether all of mycorrhizal type proportions from natural land cover types in these ecoregions. data points or raster cells36 were indeed located in forest areas. This was done using For ecoregions naturally harbouring more than one vegetation type (e.g., grass- the ESA land cover categories data33. All data points that were located out of the lands, shrublands and forests), we considered that the expected vegetation repre- current areas registers by ESA33 as forests were excluded from the analysis. sents a mixture of these land cover types. This resulted in an additional dataset The West Australian map of mycorrhizal root abundance24 provides describing combinations of vegetation as defined by ecoregion × continent × land information about the percentage of total biomass of plant roots featuring AM and cover without croplands (Supplementary Data 7 for continent data and Supple- EcM root colonization and about the percentage of non-mycorrhizal root biomass. mentary Data 8 for island data). Maps of potential mycorrhizal type distribution in We used this data as a proxy for biomass fractions of AM, EcM and NM plants. In a cropland-free world (Supplementary Fig. 4) were created based on this dataset, order to quantify the differences between the validation datasets and our maps, we following the aforementioned procedures. have calculated the Mean Averaged Error (MAE) of the difference between our In this analysis, we did not consider forest plantations to be croplands. maps and validation datasets. MAE expresses the deviation between two spatial Therefore, changes in land cover induced by forest restoration through, for datasets from the 1:1 line. For AM, EcM and NM vegetation fractions MAE is instance pine plantations or eucalypt plantations, are not addressed as vegetation 18.7%, for EcM it is 13.6%, for NM it is 4.7%, respectively. Due to a virtual absence changes induced by cropland cultivation. The total area occupied by AM tree of information about distribution of solely ErM vegetation, direct validation of the plantations exceeds that of EcM tree plantations64,65 (Supplementary Table 3), ErM maps was impossible. suggesting that exclusion of tree plantations leads to conservative estimates of the To further assess the uncertainties in the maps, we examined which land use global reduction of EcM vegetation. classes represent the data points that deviate from the observed data by more than 25% units. Our analysis showed that the large proportion of those deviations (60% Map validation. The maps of the current distributions of mycorrhizal biomass for EcM and 40% for AM) fall into those land use classes that represent a poorly fractions were validated using the datasets of forest biomass structure for Eurasia35, described mixture of evergreen or mixed forests and grasslands, i.e. ESA classes global analysis of impacts on mycorrhizas on carbon vs nitrogen dynamics19, the described as various forms of “closed to open (>15%) forest” (Supplementary USA-based analysis of mycorrhizal associations conducted with remote sensing Data 2). Further improvement of the ESA classification data will provide a techniques36, and the map of mycorrhizal root biomass in West Australia by possibility to improve precision of our maps. Brundrett24 (see Supplementary Fig. 2). Due to the rarity of datasets on field-examined mycorrhizal vegetation The data of forest biomass structure for Eurasia35 provide information on per- distributions at large special scale we had to validate our datasets used available plot tree species abundances for a large number of European sites. As the data data on plant species distribution and to accept a number of assumptions and/or contain all records obtained since the 19th century, we used only the data recorded recalculations in order to make the data comparable. These adjustments may have after 1999. Using our database of plant species and associated mycorrhizal types we affected the quality of the validation dataset and therefore the validation. assigned every tree species with its mycorrhizal type (1344 data points, Supplementary Fig. 2). This provided us with a per-site data of the relative biomass of AM and EcM trees. We used these data as proxies for AM and EcM biomass Map spatial uncertainty analysis. We quantified the uncertainty in our maps of fractions to compare with the data in our maps. We used the same approach for the mycorrhizal vegetation fractions by applying the error propagation rules to the data of Lin and co-workers19, which represent plot-based records of vegetation formulas used to calculate the biomass fractions of mycorrhizal plants per-grid cell. structure for 100 sites across the globe accompanied with data about plant- For this we used the data provided by refs. 35,60–63 to estimate the uncertainty mycorrhizal associations. associated with relative biomass of trees, shrubs, and herbaceous vegetation in each NATURE COMMUNICATIONS | (2019)10:5077 | https://doi.org/10.1038/s41467-019-13019-2 | www.nature.com/naturecommunications 7 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-13019-2 CCI land cover class. The uncertainty in the proportion of each mycorrhizal type in a All models were examined independently and evaluated based on the Akaike given Bailey × continent combination was set to 1/√n, where n = b + c/3 + g/20. In Information Criterion. this formula, b is a number of literature sources describing the vegetation compo- In all models, we assessed the total variance explained by the model by R- sition in a given Bailey ecoregion × continent combination (Supplementary Data 3), squared metric and the relative importance of each predictor using the Lindemann, while c and g are the numbers of literature sources describing vegetation composition Merenda and Gold (LMG) metric. This metric shows the proportion of variance at continent level and global levels, respectively. These latter sources were given less explained by each of model predictors within the entire variance explained. weight, because of their lower spatial explicitness, though these sources were used Given significant interactions between biome factor and biomass fractions of only if they were providing information relevant for the combination of Bailey mycorrhizal plants (Table 1), we fitted generalized linear models with EcM or AM ecoregion × continent under consideration. This procedure was applied to the maps biomass fractions as predictors and topsoil and subsoil C stocks as response of the current distribution of mycorrhizal fractions (Fig. 1) and for the distribution of variables to the data within individual biomes. The outcome of this analysis is mycorrhizal fractions in the cropland-free world (Supplementary Fig. 5). The reported in the Supplementary Table 2. As this analysis encompasses considerable resulting maps of uncertainties are shown in Supplementary Figs. 4 and 6, respec- errors on both axes of this regression, we have additionally checked if a model II tively. We calculated the uncertainties in the estimations of changes in biomass regression would yield qualitatively similar results and confirmed that this was the fractions of mycorrhizal vegetation induced by crop cultivation and pastures (Fig. 3) case (Supplementary Table 4). by applying the error propagation rule for the mathematical operation of subtraction. Additionally, we examined whether the amount of carbon stored in AM and The resulting uncertainty map is shown in Supplementary Fig. 7. EcM plants per unit area was a better predictor of soil C stocks than the fractions of AM or EcM biomass per unit area. This analysis was performed in exactly the same manner as the analysis of the impacts of AM and EcM biomass. We assessed the Carbon stocks in the aboveground mycorrhizal biomass. In order to estimate resulting models using the Akaike Information Criterion and R-squared, and the amount of carbon stored globally in the biomass of plants that belong to detected that these models were worse than the those based on the fractions of AM different mycorrhizal types, we multiplied the mycorrhizal type biomass fractions or EcM biomass per unit area. We also examined how fractions of AM or EcM by the data of the global distribution of carbon stored in aboveground plant bio- biomass per unit area, were related to soil C-to-N ratio, and detected that these mass obtained from passive microwave-based satellite observations38. As these data relationships had a pattern very similar to that of soil C. have resolution of 15 arcmin, we converted our data of mycorrhizal biomass fractions to this resolution. To calculate the total amount of carbon stored in AM, EcM, ErM and NM plants (illustrated in Fig. 2, and reported per biome in the Dealing with spatial autocorrelations. The above-described analyses were Supplementary Table 1), we summed the data on carbon stocks per unit area checked for spatial autocorrelation using a Moran’s I metric. These tests revealed a globally and per biome46. spatial autocorrelation of residuals, which is expected given the nature of our soil To assess uncertainty of estimations of carbon storage in mycorrhizal vegetation and vegetation data with spatial points nested within biomes. We examined how we used the rule of uncertainty propagation through the operation of incorporation of a spatial dependence structure affected our models. Accordingly, multiplication. For this we used the uncertainty of mycorrhizal biomass fraction we ran generalized least square (GLS) models of the same sets of predictors as data and the per biome uncertainty of biomass carbon data, as provided by Liu described in the section “Statistical analysis” of the relationship between soil carbon et al.38. Given that the data of carbon stored in vegetation of different mycorrhizal content and biomass distribution of mycorrhizal types, but included spatial cor- types is a product of two datasets, each featuring a high uncertainty38, the values of relation structures. Because calculation of the global spatial autocorrelation matix is aboveground biomass carbon stored in arbuscular, ecto-, ericoid and non- computationally intensive, we chose a Monte-Carlo type approach. For each mycorrhizal vegetation are locally associated with large uncertainties (e.g. analysis we drew a random sample of 2% of the data points (1624 points) 100 Supplementary Table 1). Therefore, we recommend to use these data exclusively for times, taking care that all biomes are included into each sample, and subsequently large-scale estimates of biomass carbon storage. However, following the central ran a GLS accounting for autocorrelation structure of the data. For this we used the limit theorem, the uncertainties of the total amount of carbon stored in vegetation algorithms proposed by Zuur et al.67. The averaged outcomes of these analyses for of each mycorrhizal type are much lower: 15, 17, 1.8 and 5.5 GT for AM, EcM, ErM the relationship between topsoil C, AM, and EcM fractions of vegetation biomass and NM vegetation, respectively. are shown in the Supplementary Table 5. The analyses accounting for auto- correlations yielded the same conclusions as the analyses that did not account for autocorrelations: there is a positive relationship between biomass fraction of EcM Statistical analysis. We examined whether soil carbon content was related to the plants in vegetation and soil C, and this relationship is biome-independent (non- biomass fractions of AM and EcM plants, using generalized linear model regres- significant Biome × EcM interaction); in contrast, the relationship between biomass sions, with topsoil (0–20 cm) or subsoil (20–60 cm and 60–100 cm) C content per fraction of AM plants in vegetation and soil C was idiosyncratic and biome- m2 as a response variable, and biome and ecosystem biomass fraction of AM or dependent. EcM plants as predictors. The data for soil carbon content in the top 20 cm of soil Including spatial coordinates explicitly in the model, or explicit accounting for were obtained from the ISRIC-WISE Soil Property Databases at a resolution of 30 autocorrelation matrix, can be problematic to interpret because covariation arcsec47. As we were interested in relationships between AM or EcM coverage and between spatial coordinates and environmental variables can obscure the soil carbon in natural vegetated environments, we excluded urban and agricultural interpretation of the relative importance of the predictors68. Given that the main areas, lakes and “Rock and Ice” areas according to the ESA land cover categories33, goal of our analysis was to detect the global links between AM and EcM dominance from the analysis. We also excluded wetlands, inundated areas, and extremely dry and soil C, but not to predict this relationship in new areas or under future climatic regions (deserts and Mediterranean areas), because we considered that harsh scenarios we prioritized the interpretation of the models rather than their abiotic conditions instead of biotic interactions are likely to shape biogeochemical predictive power, and therfore report the outcomes of the models that do not cycles in these areas; and we excluded five land cover categories for which our maps account for autocorrelations (Table 1). showed higher uncertainties, i.e. those where the extent of forest vs grassland cover was unclear. Those land cover categories included “Tree cover, broadleaved, evergreen, closed to open (>15%)”, “Tree cover, broadleaved, deciduous, closed to Reporting summary. Further information on research design is available in open (>15%)”, “Tree cover, needleleaved, evergreen, closed to open (>15%)”, “Tree the Nature Research Reporting Summary linked to this article. cover, needleleaved, deciduous, closed to open (>15%)” and “Sparse vegetation (tree, shrub, herbaceous cover)”. Data for biome types were obtained from the map of terrestrial biomes46. To Data availability All the authors’ data used to create the maps of mycorrhizal vegetation biomass, and all the create a more balanced dataset, we combined all natural grasslands into one single authors’ data underlying figures and tables in the manuscript text and supplementary category and all coniferous forests into another single category. To minimize the material are available as Source Date files. All maps presented in this paper have been impacts of imprecise estimations of biome borders we excluded from forest biome areas that we recognized as grasslands according to the ESA land cover data33. deposited into the DRYAD Digital Repository at https://doi.org/10.5061/dryad.47d7wm387. Similarly, we excluded grassland biome areas that were fully covered by forests and All the links to the publicly available datasets used within the process of assembly of the shrubs according to ref. 33. However, given that AM-EcM-NM interactions are maps of mycorrhizal vegetation, in addition to the authors’ data, are provided in the likely to relate to C stocks in the areas featuring forest-grassland and forest-tundra reference list of the paper. mosaics, we kept such areas in the dataset. As soil carbon content data is known to feature high uncertainties47, we opted to run the analyses at a resolution of 1 Code availability arcdegree, which allow us to grasp large-scale tendencies, while reducing the All the codes and necessary data files used to the generate maps of mycorrhizal problem of P-value fallacy66. vegetation biomass are available at the GitHub repository at https://github.com/ Our models comprised topsoil and subsoil carbon content as a response nasoudzilovskaia/Soudzilovskaia_NatureComm_MycoMaps. variable. As predictors, we examined the following combinations of data: (i) biome, biomass fraction of EcM plants, and their interactions, (ii) biome, biomass fraction of AM plants, and their interactions, (iii) biome, amount of carbon stored in the Received: 3 November 2018; Accepted: 3 October 2019; aboveground biomass of EcM plants (product of our EcM map and ref. 38), and their interactions, and (iv) biome, amount of carbon stored in the aboveground biomass of AM plants (product of our AM map and ref. 38), and their interactions. 8 NATURE COMMUNICATIONS | (2019)10:5077 | https://doi.org/10.1038/s41467-019-13019-2 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-13019-2 ARTICLE References 31. Ellis, E. C., Goldewijk, K. K., Siebert, S., Lightman, D. & Ramankutty, N. 1. Smith, S. E. & Read, D. J. Mycorrhizal Symbiosis (Academic Press, 2008). Anthropogenic transformation of the biomes, 1700 to 2000. Glob. Ecol. 2. Brundrett, M. C. Mycorrhizal associations and other means of nutrition of Biogeogr. 19, 589–606 (2010). vascular plants: understanding the global diversity of host plants by resolving 32. Bailey, R. G. Ecoregions. The Ecosystem Geography of the Oceans and conflicting information and developing reliable means of diagnosis. Plant Soil Continents (Springer New York, 2014). 320, 37–77 (2009). 33. ESA. CCI Land cover map 2015 (https://www.esa-landcover-cci.org/) (2017). 3. Read, D. J. Mycorrhizas in ecosystems. Experientia 47, 376–391 (1991). 34. Soudzilovskaia, N. et al. FungalRoot: global online database of plant 4. Brundrett, M. C. Mycorrhizas in natural ecosystems. Adv. Ecol. Res. 21, mycorrhizal associations. Preprint at https://www.biorxiv.org/content/ 171–313 (1991). 10.1101/717488v1 5. Leake, J. R. et al. Networks of power and influence: the role of mycorrhizal 35. Schepaschenko, D. et al. A dataset of forest biomass structure for Eurasia. Sci. mycelium in controlling plant communities and agroecosystem functioning. Data 4, 170070 (2017). Can. J. Bot. 82, 1016–1045 (2004). 36. Fisher, J. B. et al. Tree–mycorrhizal associations detected remotely from 6. Bennett, J. A. et al. Plant-soil feedbacks and mycorrhizal type influence canopy spectral properties. Global Change Biol. 22, 2596–2607 (2016). temperate forest population dynamics. Science 355, 181–184 (2017). 37. Hartley, A. J. et al. Uncertainty in plant functional type distributions and its 7. Teste, F. P. et al. Plant-soil feedback and the maintenance of diversity in impact on land surface models. Remote Sens. Environ. 203, 71–89 (2017). Mediterranean-climate shrublands. Science 355, 173–176 (2017). 38. Liu, Y. Y. et al. Recent reversal in loss of global terrestrial biomass. Nature 8. Shi, Z., Wang, F. & Liu, Y. Response of soil respiration under different Climate change. Nat. Clim. Change 5, 470–474 (2015). mycorrhizal strategies to precipitation and temperature. J. Soil Sci. Plant Nutr. 39. Wang, B. & Qiu, Y. L. Phylogenetic distribution and evolution of mycorrhizas 12, 411–420 (2012). in land plants. Mycorrhiza 16, 299–363 (2006). 9. Wurzburger, N., Brookshire, E. N. J., McCormack, M. & Lankau, R. 40. Thompson, J. Decline of vesicular-arbuscular mycorrhizae in long fallow Mycorrhizal fungi as drivers and modulators of terrestrial ecosystem disorder of field crops and its expression in phosphorus deficiency of processes. New Phytologist 213, 996–999 (2017). sunflower. Aust. J. Agric. Res. 38, 847–867 (1987). 10. Tedersoo, L. & Bahram, M. Mycorrhizal types differ in ecophysiology and 41. Mbuthia, L. W. et al. Long term tillage, cover crop, and fertilization effects on alter plant nutrition and soil processes. Biol. Rev. 94, 1857–1880 (2019). microbial community structure, activity: Implications for soil quality. Soil Biol. 11. Averill, C., Turner, B. L. & Finzi, A. C. Mycorrhiza-mediated competition Biochem. 89, 24–34 (2015). between plants and decomposers drives soil carbon storage. Nature 505, 42. Jo, I. et al. Neglecting plant–microbe symbioses leads to underestimation of 543–545 (2014). modeled climate impacts. Sci. Adv. (in the press). 12. Phillips, R. P., Brzostek, E. & Midgley, M. G. The mycorrhizal-associated 43. Cotrufo, M. F., Wallenstein, M. D., Boot, C. M., Denef, K. & Paul, E. The nutrient economy: a new framework for predicting carbon-nutrient couplings Microbial Efficiency-Matrix Stabilization (MEMS) framework integrates plant in temperate forests. New Phytologist 199, 41–51 (2013). litter decomposition with soil organic matter stabilization: do labile plant 13. Terrer, C., Vicca, S., Hungate, B. A., Phillips, R. P. & Prentice, I. C. inputs form stable soil organic matter? Glob. Change Biol. 19, 988–995 (2013). Mycorrhizal association as a primary control of the CO2 fertilization effect. 44. Cotrufo, M. F. et al. Formation of soil organic matter via biochemical and Science 353, 72–74 (2016). physical pathways of litter mass loss. Nat. Geosci. 8, 776–779 (2015). 14. Read, D. J., Leake, J. R. & Perez-Moreno, J. Mycorrhizal fungi as drivers of 45. Poirier, V., Roumet, C. & Munson, A. D. The root of the matter: linking root ecosystem processes in heathland and boreal forest biomes. Can. J. Bot. 82, traits and soil organic matter stabilization processes. Soil Biol. Biochem. 120, 1243–1263 (2004). 246–259 (2018). 15. Terrer, C. et al. Nitrogen and phosphorus constrain the CO2 fertilization of 46. Olson, D. M. et al. Terrestrial Ecoregions of the World: A New Map of Life on global plant biomass. Nat. Clim. Change 9, 684–689 (2019). Earth: a new global map of terrestrial ecoregions provides an innovative tool 16. Craig, M. E. et al. Tree mycorrhizal type predicts within-site variability in the for conserving biodiversity. Bioscience 51, 933–938 (2001). storage and distribution of soil organic matter. Global Change Biol. 24, 47. Batjes, N. Harmonized soil property values for broad-scale modelling 3317–3330 (2018). (WISE30sec) with estimates of global soil carbon stocks. Geoderma 269, 61–68 17. Brzostek, E. R., Fisher, J. B. & Phillips, R. P. Modeling the carbon cost of plant (2016). nitrogen acquisition: Mycorrhizal trade‐offs and multipath resistance uptake 48. Newsham, K., Upson, R. & Read, D. Mycorrhizas and dark septate root improve predictions of retranslocation. J. Geophys. Res. 119, 1684–1697 (2014). endophytes in polar regions. Fungal Ecol. 2, 10–20 (2009). 18. Shi, M. et al. Neglecting plant–microbe symbioses leads to underestimation of 49. Hengl, T. et al. SoilGrids250m: Global gridded soil information based on modeled climate impacts. Biogeosciences 16, 457–465 (2019). machine learning. PLoS ONE 12, e0169748 (2017). 19. Lin, G., McCormack, M. L., Ma, C. & Guo, D. Similar below-ground carbon 50. Batjes, N. H. et al. WoSIS: providing standardised soil profile data for the cycling dynamics but contrasting modes of nitrogen cycling between world. Earth System Science. Data 9, 1–14 (2017). arbuscular mycorrhizal and ectomycorrhizal forests. New Phytologist 213, 51. Steidinger, B. S. et al. Climatic controls of decomposition drive the global 1440–1451 (2017). biogeography of forest-tree symbioses. Nature 569, 404–408 (2019). 20. Averill, C. & Hawkes, C. V. Ectomycorrhizal fungi slow soil carbon cycling. 52. Sanderman, J., Hengl, T. & Fiske, G. J. Soil carbon debt of 12,000 years of Ecol. Lett. 19, 937–947 (2016). human land use. Proc. Natl Acad. Sci. USA 36, 9575–9580 (2017). 21. Terrer, C. et al. Ecosystem responses to elevated CO2 governed by plant–soil 53. Averill, C., Dietze, M. C. & Bhatnagar, J. M. Continental-scale nitrogen interactions and the cost of nitrogen acquisition. New Phytologist. 217, pollution is shifting forest mycorrhizal associations and soil carbon stocks. 507–522 (2018). Global Change Biol. 24, 4544–4553 (2018). 22. Fernandez, C. W. & Kennedy, P. G. Revisiting the ‘Gadgil effect’: do interguild 54. Ekblad, A. et al. Production and turnover of ectomycorrhizal extramatrical fungal interactions control carbon cycling in forest soils? New Phytologist 209, mycelial biomass and necromass under elevated CO2 and nitrogen 1382–1394 (2016). fertilization. New Phytologist 211, 874–885 (2016). 23. Näsholm, T. et al. Are ectomycorrhizal fungi alleviating or aggravating 55. Zhu, K., McCormack, M. L., Lankau, R. A., Egan, J. F. & Wurzburger, N. nitrogen limitation of tree growth in boreal forests? New Phytologist 198, Association of ectomycorrhizal trees with high carbon-to-nitrogen ratio soils 214–221 (2013). across temperate forests is driven by smaller nitrogen not larger carbon stocks. 24. Brundrett, M. C. Distribution and evolution of mycorrhizal types and other J. Ecol. 106, 524–535 (2018). specialised roots in Australia. Ecol. Stud. 230, 361–394 (2017). 56. Karst, J., Marczak, L., Jones, M. D. & Turkington, R. The mutualism- 25. Menzel, A. et al. Distribution patterns of arbuscular mycorrhizal and non- parasitism continuum in ectomycorrhizas: a quantitative assessment using mycorrhizalplant species in Germany. Perspect. Plant Ecol. Evolution Syst. 21, meta-analysis. Ecology 89, 1032–1042 (2008). 78–88 (2016). 57. ORNL DAAC. Dataset “Bailey Ecoregions of the Continents (Province)” 26. Swaty, R., Michael, H. M., Deckert, R. & Gehring, C. A. Mapping the potential https://doi.org/10.3334/ORNLDAAC/1388 (2017). mycorrhizal associations of the conterminous United States of America. 58. Kohout, P. Biogeography of ericoid mycorrhiza. Ecol. Stud. 230, 179–193 Fungal Ecol. 24, 139–147 (2016). (2017). 27. Tedersoo, L. et al. Global diversity and geography of soil fungi. Science 346, 59. Tedersoo, L. & Brundrett, M. C. Evolution of ectomycorrhizal symbiosis in 1078–1088 (2014). plants. Ecol. Stud. 230, 407–467 (2017). 28. Davison, J. et al. Global assessment of arbuscular mycorrhizal fungus diversity 60. Pan, Y., Birdsey, R. A., Phillips, O. L. & Jackson, R. B. The structure, reveals very low endemism. Science 349, 970–973 (2015). distribution, and biomass of the world’s forests. Annu. Rev. Ecol. Evolution 29. Fisher, J. B. et al. Modeling the terrestrial biosphere. Annu. Rev. Environ. Syst. 44, 593–622 (2013). Resour. 39, 91–123 (2014). 61. Lamotte, M. in Tropical Ecological Systems 179–222 (Springer, 1975). 30. Soudzilovskaia, N. A. et al. Quantitative assessment of the differential impacts 62. Ruiz-Benito, P. et al. Stand structure and recent climate change constrain of arbuscular and ectomycorrhiza on soil carbon cycling. New Phytologist 208, stand basal area change in European forests: a comparison across boreal, 280–293 (2015). temperate, and Mediterranean biomes. Ecosystems 17, 1439–1454 (2014). NATURE COMMUNICATIONS | (2019)10:5077 | https://doi.org/10.1038/s41467-019-13019-2 | www.nature.com/naturecommunications 9 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-13019-2 63. Poulter, B. et al. Plant functional type classification for earth system models: M.C.B. provided data for map validation. All authors contributed to preparation of results from the European Space Agency’s Land Cover Climate Change manuscript drafts. Initiative. Geosci. Model Dev. 8, 2315–2328 (2015). 64. FAO. Global Forest Resource Assessement (2000). 65. FAO, OECD. Oilseeds and Oilseed Products. OECD-FAO Agricultural Outlook Competing interests The authors declare no competing interests. (OECD Publishing, 2013). 66. Hoenig, J. M. & Heisey, D. M. The abuse of power: the pervasive fallacy of power calculations for data analysis. Am. Statistician 55, 19–24 (2001). Additional information 67. Zuur, A. F., Ieno, E. N., Walker, N. J., Saveliev, A. A. & Smith, G. M. Mixed Supplementary information is available for this paper at https://doi.org/10.1038/s41467- Effects Models and Extentsions in Ecology with R (Springer, 2009). 019-13019-2. 68. Dormann, C. F. (2007) Effects of incorporating spatial autocorrelation into the analysis of species distribution data. Glob. Ecol. Biogeogr. 16, 129–138 (2007). Correspondence and requests for materials should be addressed to N.A.S. Reprints and permission information is available at http://www.nature.com/reprints Acknowledgements This study was supported by the vidi grant 016.161.318 issued to NAS by The Nether- Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in lands Organization for Scientific research and grants EcolChange, MOBERC and published maps and institutional affiliations. 1399PUT issued to L.T. by the Estonian Research Council. I.M. was supported by the H2020 projects GeoEssential (689443) and MAGIC (727698). C.T. was supported by a Lawrence Fellow award through the Lawrence Livermore National Laboratory (LLNL). This work was performed under the auspices of the US Department of Energy by LLNL under contract DE-AC52-07NA27344 and was supported by the LLNL-LDRD Program Open Access This article is licensed under a Creative Commons under Project No. 20-ERD-055. J.B.F. was supported by funding from DOE BER TES and Attribution 4.0 International License, which permits use, sharing, NASA IDS. The authors thank Daisy Brickhill for editing the paper. Michala Phillips adaptation, distribution and reproduction in any medium or format, as long as you give assisted J.B.F. in data transfer. J.B.F. contributed to this research, in part, from the Jet appropriate credit to the original author(s) and the source, provide a link to the Creative Propulsion Laboratory, California Institute of Technology, under a contract with the Commons license, and indicate if changes were made. The images or other third party National Aeronautics and Space Administration, California Institute of Technology. Government sponsorship is acknowledged. The authors also thank Milagros Barcelo, material in this article are included in the article’s Creative Commons license, unless Stijn Vaessen and Jinhong He who helped in compiling data on plant–mycorrhizal status. indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from Author contributions the copyright holder. To view a copy of this license, visit http://creativecommons.org/ N.A.S. and C.T. generated the idea. N.A.S., P.M.v.B., L.T. and C.T. designed the study, licenses/by/4.0/. and analysed the data. N.A.S. wrote the paper. N.A.S. and I.M. prepared the figures, with I.M. leading the role. L.T. provided mycorrhizal type assignments to vegetation. C.T., I.M., N.C.d.S. and M.v.Z. performed geoinformatics data processing. M.L.M., J.B.F. and © The Author(s) 2019 10 NATURE COMMUNICATIONS | (2019)10:5077 | https://doi.org/10.1038/s41467-019-13019-2 | www.nature.com/naturecommunications