The fungal collaboration gradient drives root trait distribution and ecosystem processes in a tropical montane forest

Plant roots have a large diversity of form and function, which is also related to their degree of mycorrhizal association. This is known as the fungal collaboration gradient, where thinner roots acquire resources by themselves, and thicker roots depend on mycorrhizae. In this study, we, for the first time, implement the fungal collaboration gradient in a trait-based dynamic vegetation model (DVM, LPJ-GUESS-NTD). We test if the DVM can predict fine-root-trait distributions and estimate the effects of arbuscular-mycorrhiza-fungus (AMF)-mediated nutrient uptake on ecosystem processes along an elevation gradient in a tropical montane forest in southern Ecuador. The model reproduces the observed fine-root traits of specific root length (SRL) and AMF colonization along the elevation gradient, which ranges from low AMF colonization at 1000 m (25 %) to high AMF colonization at 3000 m (61 %). When AMF-mediated nutrient uptake is deactivated, site average biomass values are reduced by up to 80 %. Accounting for AMF-related belowground traits also affects simulated community leaf traits, suggesting linkages between below- and aboveground traits as AMF promotes more leaf-acquisitive traits. In addition, deactivation of AMF uptake reduced simulated soil C stocks by up to 68 %. The model suggests that the collaboration gradient has a substantial influence on vegetation diversity and functioning as well as soil carbon in the study system. We thus advocate more explicit treatment of fine-root traits and mycorrhizae in DVMs. The model scheme here is based on general trade-offs and could be implemented in other DVMs and be tested for other study regions.

1 Introduction

Belowground processes are becoming an increasingly important research topic in ecology, due in part to their overwhelming role in regulating biogeochemical cycles (Beillouin et al., 2023). Soils at depths of up to 200 cm store an estimated 2400 Pg C globally, highlighting their importance in the carbon cycle (Batjes, 1996). This is nearly 9 times the amount stored in global forests (Santoro et al., 2021). Fine roots, the plant's interface with soil (as opposed to coarse roots, which have a more structural role), play an important role in driving ecosystem processes (Bardgett et al., 2014; Weigelt et al., 2021). In the tropics, fine roots store up to 50 Mg C ha⁻¹ (Jackson et al., 1997), with a productivity of around 6 $\mathrm{Mg}\mathrm{C}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ (Finér et al., 2011), and are a major input to soil C stocks (Rasse et al., 2005).

Fine-root morphological traits play a critical role in nutrient and water uptake. These traits may also shape species coexistence and thus community composition in specific environments (Nie et al., 2013). Biotic interactions of fine roots and other organisms are widespread, and conservative estimates suggest that ca. 20 000 plant species could rely on soil biota to persist in natural and especially nutrient-poor environments (Van Der Heijden et al., 2008). For instance, mycorrhizae support nutrient acquisition in exchange for C (Bardgett et al., 2014; Bennett and Groten, 2022), with many plant species exhibiting fine-root traits that maximize this interaction, such as increased diameter and cortex area (Gu et al., 2014; Kong et al., 2014; Shi et al., 2023; Valverde-Barrantes et al., 2021). Mycorrhizal fungi are a major player in the global C cycle, drawing an average of 3 %–13 % but up to 50 % of the plant partner's net primary production (NPP) (Hawkins et al., 2023). This means that fine-root traits may influence larger-scale ecosystem processes.

By analysing the co-occurrence of plant traits, researchers have identified relationships between them and developed the concept of the global spectrum of plant form and function (Díaz et al., 2015; Guerrero-Ramírez et al., 2020; Kattge et al., 2020). This includes the conservation trade-off axis, where morphological and stoichiometric traits (e.g. leaf C : N, tissue densities) are related to high-productivity or high-longevity strategies (Chave et al., 2009; Díaz et al., 2015; Wright et al., 2013). This concept highlights how trade-offs in physiological and morphological traits influence species coexistence (Shipley et al., 2006). Advances in belowground trait measurements have enabled researchers to synthesize fine-root traits within the global spectrum of plant form and function (Weemstra et al., 2016, 2023; Weigelt et al., 2021). Fine-root stoichiometry traits and root tissue density were also observed to produce a similar conservation gradient; however, other fine-root morphological traits such as diameter and specific root length did not seem to align with the existing conservation axis (Carmona et al., 2021).

A promising plant form and function trait spectrum framework for understanding fine-root morphological trait variation is the fungal collaboration gradient, one of the main components of the root economics space framework proposed by Bergmann et al. (2020). In addition to concentrating nutrients produced from their own metabolism in their tissues, fine plant roots and their morphological traits significantly affect their capacity to forage for additional nutrients and water. For instance, thinner and longer fine roots have a higher total soil absorptive area per mass when compared to thicker and shorter fine roots (McCormack et al., 2015). The conserved presence of thicker fine roots in several plant species, however, shows that individuals with this trait may also receive benefits. Thicker, large-diameter fine roots that are usually associated with greater cortex : stele ratios have a higher anatomical capacity to be colonized by mycorrhizae, which can mediate the transfer of nutrients and water (Gu et al., 2014; Kong et al., 2014; Shi et al., 2023; Valverde-Barrantes et al., 2021). Although plants must transfer carbon to fungi as part of their partnership, the fungus's extensive hyphae networks can significantly boost nutrient and water absorption. This collaboration offers thicker-rooted plants an alternative strategy to relying solely on roots for the uptake of water and nutrients (Kakouridis et al., 2022). A trade-off emerges between a cheaper do-it-yourself strategy, with fine roots of high specific root length (SRL) that provide sufficient absorptive area without the need for fungal collaboration, and “outsourcing” to the mycorrhizae partner, whereby a plant with fine roots of relatively thick diameter and low surface area “pays” (i.e. transfers) carbohydrates to the fungus to benefit from the high absorptive capacities of hyphae (Bergmann et al., 2020). A spectrum of strategies thus emerges: at one end, plants invest heavily in fine-root traits that enhance independent nutrient uptake, while at the other end, they rely more on fungal partners to exchange nutrients for carbon.

Empirically based concepts connecting patterns to processes such as the spectrum of form and function (Díaz et al., 2015) are very useful for dynamic vegetation models (DVMs), which depend on generalized representations of ecology. DVMs are software simulators of plant community, physiology and edaphic processes used to conduct experiments that are difficult or impossible to produce in field settings, such as characterizations over large temporal and spatial scales from complex treatments (Prentice et al., 2007; Quillet et al., 2009; Sitch et al., 2008). Aboveground plant trait variation has been implemented in DVMs using the leaf economics spectrum, a framework that classifies leaves based on a trade-off between rapid growth and resource conservation (Wright et al., 2004). Including such frameworks in DVMs has provided insights into nutrient dynamics and community resilience (i.e. Sakschewski et al., 2015, 2016; Dantas de Paula et al., 2021). In contrast, belowground processes, despite their critical ecological roles, remain underrepresented in these models (Langan et al., 2017; Sakschewski et al., 2021). In particular, despite mycorrhizae being fundamental to plant growth, not many attempts have been made to explicitly include them in DVMs (Dantas de Paula et al., 2019, 2021; He et al., 2021; Kou-Giesbrecht et al., 2021; Thurner et al., 2023). None of these approaches used trait-based models, meaning that the coexistence of species with lower and higher mycorrhiza colonization rates was not present, and results were highly parameter dependent. Since data on mycorrhiza physiology and structure are scarce and are complex to measure (Godbold et al., 2006; Van Der Heijden et al., 2008), these models may be producing results that are unrealistic. Nevertheless, including variation in fine-root traits and mycorrhizal colonization in DVMs holds the potential to improve predictions of the vegetation response to climate change, as was the case with leaf and wood traits (Sakschewski et al., 2016). Due to their increased absorption capabilities, mycorrhiza can effectively support plant survival in harsher environments, increasing resilience to disturbances such as droughts (Das and Sarkar, 2024).

Here, we derived parameters from the literature and field measurements to implement fine-root-trait variations and the fungal collaboration gradient (FCG) in a trait-based DVM (LPJ-GUESS-NTD). This DVM includes N and P cycles, the latter of which is particularly important for tropical forests (Du et al., 2020). Previously, an LPJ-GUESS-NTD version that included aboveground trait diversity, N and P cycling, and a simple mycorrhizal representation was applied to a well-researched elevation gradient at tropical montane forest (TMF) biodiversity hotspot sites in the southern Ecuadorian Andes. The model reproduced an observed gradient of aboveground vegetation traits and C processes along the elevation gradient, showing that nutrient dynamics might play a very important role in vegetation changes and C cycling along this gradient (Dantas de Paula et al., 2021). However, variation in fine-root traits was not accounted for, and mycorrhizal symbioses were highly simplified and prescribed. Such an approach can only yield limited site-specific insights about belowground traits and mycorrhiza. To address these gaps, this study develops a dynamic approach that integrates detailed fine-root-trait and fungal interaction data into DVMs for broader ecological applications. It builds upon more detailed ecophysiological processes whose understanding is currently not entirely generalizable to DVM applications but which could be essential for linking above- and belowground plant traits. We applied the new model version to the Andes TMF and tested it against field observations of fine-root and arbuscular mycorrhizal fungus (AMF) measurements along the altitudinal gradient (Homeier and Leuschner, 2021; Pierick et al., 2021). This local setting was chosen not only to build on previous results but also because we were able to apply the FCG using only one mycorrhiza type (AMF). We aim here to test with our model implementation the general hypothesis that the FCG is an important factor behind the observed fine-root trait distribution, forest biomass and productivity. More specifically, we hypothesize that (1) in line with the mycorrhizal colonization gradient and field measurements of fine-root traits, as available nutrients decrease with elevation, simulated community average values of SRL decrease, diameters increase and colonization rates by AMF increase when the FCG is active. Next, we removed the mycorrhizal fungi in a simulated exclusion experiment and expect that (2) in the absence of AMF, plant biomass and productivity would differ from observations, and SRL between the different sites of the elevation gradient would be similar. In other words, this latter result would imply that AMF drives morphological fine-root diversity.

2 Materials and methods

2.1 Site description

To test our hypotheses and to drive and evaluate our simulations, we used site-level data from three different elevations (1000, 2000 and 3000 m a.s.l.) in the Cordillera Real in the eastern range of the south Ecuadorian Andes. These three sites harbour old-growth forests and are located within or close to Podocarpus National Park, where a wealth of biotic and abiotic measurements have been carried out since the early 2000s (Beck et al., 2008; Bendix et al., 2013, 2021). The terrain is predominantly steep and rugged, characterized by ridges and valleys. The primary parent materials for soil formation include Paleozoic metamorphosed schists and sandstones interspersed with quartz veins. Cambisols are the most prevalent soil type at lower elevations, while Planosols and Histosols dominate at higher elevations around 2450 m (Wilcke et al., 2008). The main recurring disturbance in the area is landslides (Wilcke et al., 2003). The annual mean temperature at these elevations decreases from around 20 °C (1000 m) to ∼10 °C (3000 m), and annual precipitation increases from 2230 mm (1000 m) and 1950 mm (2000 m) to 4500 mm (3000 m) (Bendix et al., 2021; Dietrich et al., 2016; Moser et al., 2007). Due to several abiotic and biotic factors, of which temperature is the most relevant, decomposition and mineralization of organic matter decline with elevation along this range, resulting in an increase in the thickness of the organic layer from around 5 cm at the 1000 m site to more than 50 cm at the 3000 m site. Conversely, inorganic nutrient stocks, which were studied as the plants' main nutrition source, decrease with elevation. In particular, nitrogen (both inorganic forms of N and ${\mathrm{N}}_{\mathrm{i}}={\mathrm{NH}}_{\mathrm{4}}-\mathrm{N}+{\mathrm{NO}}_{\mathrm{3}}-\mathrm{N}$) stocks vary from 5.14 kg N_i ha⁻¹ at 1000 m to 13.02 kg N_i ha⁻¹ at 2000 m and 0.64 kg N_i ha⁻¹ at 3000 m on average (Dantas de Paula et al., 2021; Velescu and Wilcke, 2020). Bray extractable inorganic phosphorus (P_i) measurements show no significant differences between elevation sites (Dietrich et al., 2016), suggesting that P limitation may play a minor or co-limiting role relative to N on this gradient (Homeier et al., 2012). Nutrient availability and the resulting limitation of plant growth along the elevation gradient are considered important drivers of the diversity of ecosystem processes and species functional traits in this region: for example, with increasing elevation, forest productivity and biomass stocks decrease (Homeier and Leuschner, 2021), and the community trait distribution becomes more conservative (Homeier et al., 2021; Pierick et al., 2023, 2024) (e.g. lower specific leaf area, lower fine-root nutrients).

Extensive floristic inventories in the area allowed us to distinguish between three main forest types: the forest types studied are classified as evergreen pre-montane forest (1000 m), evergreen lower-montane forest (2000 m) and evergreen upper-montane forest (3000 m). The floristic composition changes rapidly with elevation, as most tree species in the study area are only found in narrow elevational belts (Homeier et al., 2008). The evergreen pre-montane rainforest at the lowermost study site reaches 40 m in height, with common tree families being Fabaceae, Moraceae, Myristicaceae, Rubiaceae and Sapotaceae. It is replaced at 1300–2100 m by a smaller-statured lower-montane rainforest, with Euphorbiaceae, Lauraceae, Melastomataceae and Rubiaceae as characteristic tree families, and above 2100 m by upper-montane rainforest with a canopy height rarely exceeding 8–10 m. Common tree families of the latter forest type are Aquifoliaceae, Clusiaceae, Cunoniaceae and Melastomataceae. With increasing elevation, forest biomass and productivity decrease, whereas the root–shoot ratio increases (Homeier and Leuschner, 2021). Concerning the belowground perspective, which is the main focus of this work, it is known at these sites that average SRL decreases with increasing elevation, from around 3000 cm g⁻¹ at 1000 m to 1500 cm g⁻¹ at 3000 (Pierick et al., 2023). The percentage of root length colonized by arbuscular mycorrhizal fungus (AMF), which interacts with almost all species in the study area (Kottke and Haug, 2004), increases from an average of 25 % at the lowest elevation to 61 % at the highest (Camenzind et al., 2016). This altitudinal trend, as well as nutrient addition experiments at the sites (Camenzind et al., 2014), suggests that AMF has an important role in plant N uptake in our studied TMF. This is in contrast to the common view that AMF is mostly relevant for P, which has been challenged (Hodge and Storer, 2015).

Principal component patterns for SRL, fine-root diameter, tissue density and N content (Pierick et al., 2021, 2024) are very similar to the global ones described in Bergmann et al. (2020), where fine-root tissue density and N form one axis (i.e. fast-slow gradient) and SRL with fine-root diameter another, in line with the FCG.

2.2 Model description

To simulate the trait distributions in the Andes tropical montane forest (TMF), we use the LPJ-GUESS DVM (Smith et al., 2014) with the NTD (nutrient-trait dynamics) implementation (Dantas de Paula et al., 2021). This model includes individual representations of each tree, a Farquar-based photosynthesis implementation, detailed tree population dynamics (establishment, mortality, disturbances) and abiotic competition processes (light, water, nutrients, space; Smith et al., 2001, 2014). In the NTD version of LPJ-GUESS (see Dantas de Paula et al., 2021, for the full description), plant diversity is included through trait variation as random values of specific leaf area (SLA) and wood specific gravity (WSG) when new tree saplings are established. These traits are related to further traits, and the trait–trait relationships define the major trade-off axes, such as SLA to tissue C : N ratios. For the model in this study, most trait–trait relationships and tissue stoichiometry were not parameterized from global data but from measurements at the study site (Dantas de Paula et al., 2021). The trait combinations each plant receives upon establishment impact its competitive success compared to other individuals in any specific environment, with less-adapted individuals suffering, e.g. having lower growth than competitors and, therefore, higher growth efficiency mortality (the mortality component in the model related to negative net primary productivity (NPP) values).

Soil organic matter (SOM) dynamics were adopted from the CENTURY model (Parton et al., 1988, 1993, 2010), with organic matter pools and the N and P cycles also included. Organic matter enters the SOM model through vegetation litter input, and total soil organic C for a particular point in time can be estimated by summing up all C in the SOM pools. The C : N and C : P ratios of litter, its size (e.g. leaf versus coarse woody debris), soil temperature, soil humidity, and available soil nutrients influence SOM decomposition rates, soil C accumulation and nutrient dynamics. Tissue C : N and C : P ratios also determine the nutrient demand that must be met by root uptake; otherwise, photosynthesis and growth become limited by N and/or P. Nutrient limitation also drives increased C allocation to roots, and higher nutrient uptake than demand drives increased C allocation to leaves. LPJ-GUESS-NTD includes a simple implementation of AMF-mediated plant nutrient uptake (Dantas de Paula et al., 2021), which was required to reconcile measured available soil nutrient concentrations with observed vegetation structure and trait distributions. This representation was based on Kirschbaum and Paul (2002), in which vegetation demand for N and P are met by additional uptake from the surface microbial pool. However, this approach had several limitations, among which were (1) no representation of individual AMF mass per plant, i.e. a “big mushroom” approach (analogous to the “big leaf” approach; Luo et al., 2018xs); (2) a fixed AMF colonization rate for all individuals and (3) no C costs for AMF-mediated nutrient acquisition. Also, fine-root traits were fixed among individuals, meaning that the total fine-root surface area for nutrient uptake depended only on fine-root biomass, and fine-root morphological traits such as SRL were not considered. This previous implementation therefore did not provide a realistic relationship between the C invested in nutrient acquisition (as fine roots or mycorrhizae) and effective acquisition capacity. Although fine-root biomass in general terms translates into higher nutrient or water uptake capacity, there is high variation in their biomass values, and this is thought to result from fine-root morphological trait variation (Kokko et al., 1993). Even though fine-root architecture (i.e. root distribution along different soil layers) and rooting depth are considered in our model and in other models (Langan et al., 2017; Sakschewski et al., 2021), to our knowledge, morphological traits such as fine-root diameter and SRL have not been implemented.

2.3 Implementation of the fungal collaboration gradient

We consider fine-root surface area to be a better proxy for uptake capacity than fine-root biomass. We consider fine roots as those having a diameter of less than 2 mm, as defined in our reference field samples (Pierick et al., 2021, 2024). The inclusion of SRL as an individual trait permits the calculation of the total fine-root surface area using the following equation:

$$\begin{array}{}\text{(1)}& A_{\text{root}}=\text{SRL}{d}_{\text{root}}\mathit{\pi }{C}_{\text{root}},\end{array}$$

where A_root is the total (absorptive) fine-root surface area (m²), SRL is in m (kg C)⁻¹, d_root is the fine-root diameter in metres calculated from SRL (see below) and C_root is the fine-root C biomass (kg C m⁻²). Individuals with higher values of SRL and thinner roots (low diameter) have larger values of A_root (since SRL and d_root have a nonlinear relationship) per fine-root mass. The relationship between SRL and d_root was produced using species data measured at our study site (Pierick et al., 2021).

We included SRL as a randomized trait in the tree's individual establishment in addition to the model's existing traits of SLA and WSG, as well as the AMF C mass as a mass pool for each tree. In this updated model version, AMF biomass and area (C_AMF and A_AMF) are implemented as an extension of the root system for each individual tree. For herbaceous vegetation, individuals are not distinguished in LPJ-GUESS. Accordingly, C_AMF and A_AMF are not individual-based for herbaceous vegetation, with only one pool for C₃ and one pool for C₄ herbaceous vegetation. The total hyphal surface area for each individual's AMF is calculated as in Eq. (1), substituting A_root, d_root and C_root for A_AMF, d_AMF and C_AMF, respectively, and SRL for specific hyphal length (SHL). Values for SHL and d_AMF (Table A1) are fixed for all individuals (no fungal trait variation) and are based on Raven et al. (2018). Due to their very thin hyphae, AMFs typically have specific surface areas that are several times larger than those of fine roots with the same C mass (Raven et al., 2018). The total soil area a plant individual can explore for nutrients is thus equivalent to A_root+A_AMF. Cooperation with AMF benefits an individual plant, increasing the surface area for its nutrient uptake to levels that cannot be reached by its fine roots alone. However, cooperation with AMF implies fungal C_AMF costs that must be borne by the plant. We expect then that in more-nutrient-limited environments, such as in higher-elevation areas, mycorrhiza colonization will be high. On the other hand, at lower-elevation sites where light competition plays a large role, investment in mycorrhiza will not lead to higher fitness, leading to lower colonization rates.

AMF C dynamics for each individual's C_AMF were implemented similarly to the structural plant C compartments, i.e. leaf, wood and fine roots. AMF respiration is included and depends on the current daily C_AMF and a fixed AMF C : N ratio value (Orwin et al., 2011). C is allocated to AMF biomass (C_AMF) at each time step and is calculated as

$$\begin{array}{}\text{(2)}& \mathrm{\Delta }{C}_{\text{AMF}}=\mathrm{\Delta }{C}_{\text{inc}}\cdot r_{\text{max}}\cdot m_{\text{col}},\end{array}$$

where ΔC_inc is the total biomass increment for this individual at this time step (NPP – respiration), and r_max is a calibrated parameter defining the maximum fraction of ΔC_inc (0–1) transferred to AMF for growth. The individual's mycorrhizal colonization rate m_col (0–1) was defined from a correlation with d_root based on global trait data from the GROOT database (Guerrero-Ramírez et al., 2020) since our study area had no species-level (only site-level) measurements of mycorrhizal colonization rates. Based on the global m_col−d_root relationship, we assume that increasing colonization rates will result in a larger C transfer from plant to AMF, in other words, increased C costs. The parameter r_max is therefore calibrated to establish the slope of the assumed linear relationship between AMF colonization and C increment. We explore its values from 0 to 1 to determine values where maximum benefit for the plant occurs. A graphic visualization of the relationship between m_col and r_max can be seen in Fig. A1.

In the new model, ecological filtering determines the optimal AMF colonization rates, leading to the highest biomass growth rates under given environmental conditions. The individual trees with the optimal AMF colonization then outcompete other individuals. Competition occurs for light, water, N and P. This approach to C transfer between plants and AMF differs from those in which N and P are exchanged for fixed or varying costs of C (Allen et al., 2020; Fisher et al., 2010; Reichert et al., 2023); instead we use a more explicit representation of AMF physiology. Although fungal C demand is represented here, we did not include N or P demand and the uptake of the mycorrhiza themselves. Thus, all nutrients absorbed by the AMF hyphae in the model are transferred to the plant (Fig. 1).

Figure 1 Full schematic of LPJ-GUESS-NTD including the new fine-root and mycorrhiza implementations in relation to Dantas de Paula et al. (2021). For more details, refer to that publication. In the present study, arbuscular mycorrhiza fungus (AMF) was added as an individual component of a plant. A_root and its analogue A_AMF, as described by Eq. (1) and ΔC_AMF from Eq. (2), are indicated. Dark green and red elements represent soil organic matter (SOM) pools, while the purple hexagons represent mineral (plant-available) pools.

Relevant to our model approach for the FCG implementation is the advantage posed by AMF's higher nutrient absorption area and uptake capacity. Research in the literature has yielded some data on the necessary Michaelis–Menten equation parameters for model runs: the maximum uptake rate of N and P (Vmax) and the half-saturation constant (Km) (Pérez-Tienda et al., 2012; Silveira and Cardoso, 2004). These values represent the amount of nutrients taken up per fine-root mass for each time step and how well the uptake occurs under low concentrations, respectively. Values for AMF indicate that due to different N- and P-form transport systems, fungi differ in their fine-root uptake dynamics (Silveira and Cardoso, 2004; Wu et al., 2020) and particularly in the Km values, as AMF can absorb N and P under much lower concentrations than fine roots can (Table A1).

Another relevant parameter in the simulations was the turnover rate of AMF extra-radical mycelium. This has been recognized as a complex field measurement, with estimations as low as 5 to 6 d (Staddon et al., 2003) and 9 d (Godbold et al., 2006). However, here we consider the extra-radical mycelium longevity of 80 d proposed by Raven et al. (2018) based on the elongation rate of up to 3 mm d⁻¹ for AMF hyphae (Olsson and Johnson, 2005; Smith and Read, 2010).

Graphs of the data used in the trade-offs and the equations used in the trait–trait relationships can be seen in Fig. A2 in the Supplement.

2.4 Model drivers and evaluation data

Each of our three elevation sites was simulated with distinct daily data for temperature, precipitation and radiation, interpolated from measured monthly average values from climate station data at each elevation site between 1999 and 2018 (Peters and Richter, 2009; Rollenbeck et al., 2015; Bendix, 2020). Soil data and parameters were taken from the World Soil database (FAO/IIASA/ISRIC/ISS-CAS/JRC, 2012); the N and P deposition rates used were from Dantas de Paula et al. (2021), which were measured weekly during the same 1999–2018 period. Details on how driver data influence ecosystem processes can be found at Smith et al. (2014) and Dantas de Paula et al. (2021).

Minimum and maximum SLA, WSG and SRL values for the initialization of individual trees at establishment were defined using values measured in the field (Báez and Homeier, 2018; Homeier et al., 2021; Pierick et al., 2023). These minimum and maximum observed trait values consider all three elevation sites as one since we do not include dispersal or establishment limitation along this elevational range. Regarding the trade-off relationship between the traits SLA, WSG and SRL, we use for the first two the same correlations defined in Dantas de Paula et al. (2021). For SRL, we derive d_root from it, using the dataset of Pierick et al. (2023), with AMF colonization rates from d_root and m_col values from the GROOT database (Guerrero-Ramírez et al., 2020). Equations, R² values and graphs can be seen in Fig. 2. Model parameterizations can be found in Table A1.

2.5 Modelling protocol and scenarios

Following the standard LPJ-GUESS procedure (Smith et al., 2014), the LPJ-GUESS-NTD model was initialized from the bare ground and was allowed to spin up for 500 years using the driving data to build vegetation and soil C, N and P stocks. Following the spin-up period, a further 200 years was simulated to generate the modelled results by repeating the observed climate data to represent long-term mean conditions. We repeated this process for two defined scenarios: (1) no AMF collaboration (AMF-off), in which the parameter for maximum C transfer to AMF r_max is set to zero and A_total equals A_root, and (2) AMF collaboration (AMF-on), considering the new implementations defined in the current study and where uptake of nutrients depends on A_root+A_amf (Fig. 2). Each simulation was set to represent a total area of 10 ha, and the whole procedure was replicated 30 times to average random effects.

Figure 2 Scenarios used in this study to evaluate the role of the fungal collaboration gradient in shaping specific root length (SRL) distributions. In the top scenario (AMF-off), transfer of C to arbuscular mycorrhiza fungus (AMF) is deactivated (r_max=0), and plants depend only on roots to acquire nutrients. In the scenario below (AMF-on), transfer of C is active (r_max=0.5), and AMF participates in nutrient acquisition.

2.6 Model evaluation

To assess model results, we compared them to the trait and forest structure measurements also used in Dantas de Paula et al. (2021) in addition to comparisons to SRL, fine-root diameter and AMF colonization data reported in Pierick et al. (2021) and Camenzind et al. (2016). The trait distributions from these same sources as well as mean values for SLA, WSG, SRL and fine-root AMF colonization were also used to evaluate the simulated trait distribution results for each of the three elevation sites.

2.7 Sensitivity analysis to infer the maximum carbon fraction transferred to AMF

We included a sensitivity analysis of the parameter r_max (the fraction of plant NPP allocated to growth in each time step that is transferred to AMF when the colonization rate is 100 % or a rate of 1), which is the only parameter that is calibrated and not based on observations. The sensitivity analysis allowed us to assess the full range of possible r_max values (0–1) and identify optima of this parameter that maximize plant productivity and plant–fungi interactions. Given the lack of experimental data and relationships, we assumed here a linear response between the AMF colonization rate and r_max (see Eq. 2).

We executed 30 simulations for each elevation site with the r_max parameter varying from 0 to 1 and examined the effects on traits and ecosystem processes, as well as the total rate of transfer of C from plant to fungi and the C cost per N or P uptake.

3 Results

3.1 Sensitivity analysis of the r_max parameter

The model simulated an increase in average biomass and productivity for all sites, with an increasing maximum C transfer rate from plant to AMF (r_max)(Fig. 3a and b) and an increasing community average AMF colonization rate (Fig. 3c), reflecting the alleviation of N and P limitation (Fig. A3e and f) and the increasing benefit of interacting with fungi for plants. This pattern continues until an r_max of 0.25–0.5, when productivity and biomass reach their highest values. At that point, average community AMF colonization reaches a maximum, having values that are closest to observations. With increasing r_max, average colonization rates decrease, indicating a decoupling of plants and fungi. For this reason, we picked a value of 0.5 r_max for subsequent analysis. We consider an r_max of 0.5 as the value where the largest benefit for the plant occurs and when interaction with AMF occurs; thus, the AMF-on scenario is r_max=0.5. More details and results for the sensitivity analysis of r_max can be seen in Appendix A.

Figure 3 Sensitivity analysis of the parameter r_max, the maximum C allocation to mycorrhizal growth. Each point represents 200-year averages using a particular r_max. Whiskers indicate ± SD from the 200 years of each run. Red lines – 1000 m. Green lines – 2000 m. Blue lines – 3000 m.

3.2 Influence of the fungal collaboration gradient on plant traits

The decrease in SRL and increase in AMF colonization rates observed in the field data along the three elevation gradients were reproduced by the model in the AMF-on scenario. Simulated SRL and AMF colonization rates were within the confidence interval of field measurements for the 1000 and 2000 m a.s.l. sites (metres above sea level; Fig. 4a and b), ranging from 4929.31 to 3339.07 cm g⁻¹ (SRL) and from 0.29 to 48 (AMF colonization), as can be seen in Table 1. When AMF was deactivated, no differences between elevations were found for simulated SRL and AMF colonization (Fig. 4a and b).

Figure 4 Results of simulated and observed C-related values for the elevation gradient: (a) specific root length – SRL; (b) AMF colonization rate; (c) biomass (sum of above- and belowground) and (d) NPP – annual net primary production. Scenarios are AMF-on – C allocation to growth r_max=0.5 – and AMF-off – C allocation to growth r_max=0.0. Grey bars are field measurements, with whiskers indicating confidence intervals.

The distribution graphs of SRL showed clearly how the fine-root traits are shaped by the relationship with AMF in this environment (Fig. 5a). When the FCG was active (AMF-on), the SRL distribution was much closer to the observed distribution than in the AMF-off scenario (Fig. 5a). The resulting AMF-off SRL pattern was very similar to the initial uniform random distribution at the establishment, with a lack of trends and no differences in SRL between the three elevational sites. Viewing SLA in this form as well, it can be seen that the AMF-off community was significantly more conservative (ranging from 47.96 to 75.18 cm² g⁻¹) than the AMF-on one (ranging from 88.79 to 104.97 cm² g⁻¹) (Table 1). A similar pattern was found regarding leaf stoichiometry, as C : N and C : P ratios reflected SLA values and also became more conservative (larger values with a higher proportion of C) when AMF was deactivated (Table 1). WSG, on the other hand, was unaffected by the different AMF scenarios. Model performance thus improved in the AMF-on scenario, reducing median absolute errors for SRL in comparison with the AMF-off scenario. SLA, on the other hand, was slightly worse with the activation of AMF uptake (Table A2 and Fig. A5).

Figure 5 Trait density distributions for (a) specific root length, SRL, and (b) specific leaf area, SLA. The scenarios are AMF-on – C allocation to growth r_max=0.5 – and AMF-off – C allocation to growth r_max=0.0. Vertical dashed lines are the mean values for each scenario.

Table 1 Simulated results for the two AMF treatment scenarios AMF-on (0.5r_max) and AMF-off (0.0r_max) along the three elevation sites (1000, 2000 and 3000 m a.s.l.). T tests (N=30) between AMF-on and AMF-off scenarios presented significant differences for all variables. T tests between elevations were significant as well (1000 compared to 2000 and 1000 compared to 3000) except for SRL, which did not present significant differences between elevations.

3.3 Influence of the fungal collaboration gradient on carbon stocks and productivity

Simulation results with activated AMF (AMF-on scenario; r_max=0.5) exhibited on average higher biomass and productivity values than AMF-off results did. Simulated biomass in the AMF-on simulations was within the observed bounds, whereas with AMF-off, only the 1000 m site produced biomass within the observed bounds (Fig. 4). Biomass differences between the AMF-on and AMF-off scenarios were the largest for the 3000 m elevation site, where the former had $\mathrm{84.9}\mathrm{Mg}\mathrm{C}{\mathrm{ha}}^{-\mathrm{1}}\pm \mathrm{4.85}\text{SD}$ and the latter $\mathrm{16.49}\mathrm{Mg}\mathrm{C}{\mathrm{ha}}^{-\mathrm{1}}\pm \mathrm{2.62}\text{SD}$, an 80.6 % reduction (Fig. 4c and Table 1). Similar differences were found for NPP between both scenarios, and when compared to observations, the AMF-on scenario yielded better NPP estimates, in particular for the 3000 m site (Fig. 4d).

Average soil C stocks increased in the AMF-on scenario from 29.53 at 1000 m to 49.57 Mg ha⁻¹ at 3000 m, while in the AMF-off scenario they decreased from 24.08 to 15.65 Mg ha⁻¹ (Table 1). Although net N mineralization rates were slightly higher in the AMF-off scenario (41.46 $\mathrm{kg}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$) in comparison to the AMF-on one (37.80 $\mathrm{kg}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$) at the 1000 m site, they decreased by more than half at the 3000 m site, from 22.23 to 10.58 $\mathrm{kg}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ (Table 1). Litter stoichiometric properties reflected leaf C : N and C : P ratios, becoming more carbon rich (Table 1).

Total AMF biomass peaked at the 2000 m elevation site at around 0.45 Mg C ha⁻¹ and was lowest at the 1000 m elevation site, 0.35 Mg C ha⁻¹ (Fig. A3b). However, in terms relative to total vegetation biomass, AMF increased with elevation, reaching around 0.4 % of the total plant biomass stocks at the highest-elevation site (Fig. A4).

4 Discussion

Trait-based dynamic vegetation models (DVMs) represent an important tool that can be used to test theoretical frameworks of functional diversity. Our simulated results suggest that the fungal collaboration gradient (FCG) is a major driver of fine-root functional diversity and ecosystem processes. We were able to compare three contrasting environments across altitudes, and the simulations largely followed field observations in magnitude and trend but only when the FCG was activated in the model. These results have important implications for our understanding of the role of fine-root-trait syndromes in driving ecosystem processes, as well as for conservation strategies and ecosystem management.

4.1 The importance and implications of the fungal collaboration gradient for fine-root-trait diversity and ecosystem processes

Our results show that the FCG is a relevant concept for linking fine-root-trait diversity and ecosystem functioning. The FCG allowed realistic fine-root trait distributions to emerge from the simulations. The FCG also has an important role in vegetation productivity and the organic C cycle. The lower productivity and biomass of AMF-off communities suggest that the lack of mycorrhiza is detrimental to plant nutrition and productivity at the higher elevations in our study regions. The deactivation of C transfer from plant to mycorrhiza represents a long-term fungal suppression experiment from an empirical ecology perspective. This is incredibly valuable in establishing the causality of interactions, but such experiments are understandably difficult to perform under natural conditions in long-standing ecosystems. Field studies that tested either the removal of mycorrhiza or their inoculation show responses proportional to their relevance for the dominant plant species in an ecosystem (Lin et al., 2015). When mycorrhizae are relevant for a community, we can expect that fungal suppression will significantly reduce total biomass as, for example, in O'Connor et al. (2002). In that study, a heavily AMF-colonized herb species had a reduction of 60 % of biomass after the application of fungicide, which is comparable to the 80 % reduction at our highest-elevation site. A meta-analysis of AMF experimental inoculation effects on plant growth found that on average AMF increases biomass by 47 % (Wu et al., 2024). Some small-scale field studies that suppressed fungal colonization in intraspecifically diverse plant populations through molecular rather than chemical methods are perhaps able to more cleanly determine the influence of AMF on biomass and have supported altered biomass and competitive ability of plants in the absence of the ability to associate with AMF (Groten et al., 2023; McGale et al., 2020).

The presence of AMF at the higher-elevation sites was not enough to drive the biomass at the 3000 m sites to the same values observed at 1000 m. As evaluated in a climatic sensitivity analysis in a previous LPJ-GUESS-NTD publication, the gradient of biomass and productivity at our study sites is ultimately driven by temperature impacts on nutrient cycling (Dantas de Paula et al., 2021). In the model, lower decomposition rates result in less nutrient mineralization and availability, which is exacerbated by changes in leaf and litter traits (more C in relation to N and P, tougher leaves that decompose poorly) and limits productivity and biomass. Field studies complement this, indicating that with elevation, more biomass is allocated belowground, suggesting higher competition for nutrients (Leuschner et al., 2013).

Related to this impact on plant growth is the effects of AMF on soil C stocks and nutrient dynamics. This is a complex topic since in the literature the presence of AMF is related to the reduction (Wurzburger and Brookshire, 2017) or increase (Rillig et al., 2001) in soil organic C. Our modelling approach can provide insights into this contrasting effect – AMF may increase plant productivity (by alleviating nutrient limitation), which would tend to increase soil C stocks, while at the same time driving more labile litter and accelerating decomposition rates, which would reduce soil C. Since LPJ-GUESS-NTD includes both litter deposition and stoichiometry (which is a direct consequence of plant trait variation), the balance between increasing soil C input and higher soil C turnover can be effectively disentangled depending on the simulated site. For the particular area simulated in this study, our model suggests that AMF drives higher soil C stocks due to its support of higher plant productivity (Table 1), particularly at the highest-elevation site. This occurs in spite of the AMF-off scenario having lower net N mineralization rates and higher C : N and C : P litter (Table 1), which would promote lower decomposition and thus higher soil C accumulation and stocks in relation to AMF-on. Exploring whether the model could reproduce patterns of higher soil C stocks in the presence of AMF in other regions would be an important model performance test and would help to estimate how global average soil C stocks are affected by AMF. Therefore, the value and necessity of a DVM as implemented here are evident, as they expand our understanding of these results across diverse natural ecosystems and sites. Moreover, this approach is crucial at a scale that is urgently needed in the context of climate change, particularly for conservation efforts.

Our modelling approach also reveals linkages between belowground and aboveground traits. For our simulations, we observed that more-nutrient-limited environments had both low specific leaf area (SLA) and low SRL (high diameter). In addition, deactivating AMF (AMF-off) resulted in communities with significantly lower SLA values. In other words, constraining fine-root diversity (i.e. by excluding AMF) in our model resulted in aboveground changes in leaf traits and a reduction in aboveground trait diversity (SLA range). These results reveal that belowground traits and the FCG can significantly influence aboveground vegetation traits across different altitudes in a tropical system. This is new information for this system and is rarely explored in trait research or DVMs, likely because the complex interactions are hard to untangle in field experiments or observational studies.

Field studies at our elevation gradient sites agree with our modelling results and the hypothesis that plant–microbial competition for N acquisition intensifies with altitude. For instance, experimental additions of N (50 $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$) and P (10 $\mathrm{kg}\mathrm{P}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$) over 5 years resulted in increased CH₄ uptake at the 2000 and 3000 m a.s.l. sites, indicating strong nutrient competition (Martinson et al., 2021). Also, to expand on this, future studies should integrate field and laboratory experiments across different ecosystems and altitudinal gradients. Some methods offer promising avenues to improve our understanding of plant nutritional status and plant–soil competition for nutrients. For instance, stable isotope tracing (e.g. ¹⁵N labelling) can help track N uptake by plants and microbes, while enzyme activity assays can reveal microbial nutrient acquisition strategies under varying temperature conditions (Dunn et al., 2006). Additionally, soil sterilization bioassays can offer insights into the role of soil microbes that affect plant performance (Waring et al., 2016), whilst metagenomics and metatranscriptomics can provide insights into shifts in microbial communities and their functional roles in nutrient cycling (Mendes et al., 2017). Expanding empirical research in these areas will be essential to refine our understanding of plant–microbial interactions in response to changing environmental conditions.

While microbial competition for nutrients can limit plant nutrient acquisition, free-living microbial decomposers play a crucial role in breaking down evolved organic matter and mineral-associated organic matter, thereby releasing soluble nutrients to plants and synchronizing supply and demand (Fontaine et al., 2024). These microbes, including saprotrophic fungi and bacteria, drive nutrient mineralization processes that can enhance plant nutrient uptake, particularly in nutrient-limited environments such as the elevation gradient that we studied. At our sites, the N cycle is closely coupled (i.e. gross N mineralization is equal to NH₄ immobilization, and gross nitrification is equal to NO₃ immobilization), and experimental nutrient additions can alter this equilibrium (Baldos et al., 2015). This is particularly alarming due to the observed increase in anthropogenic N deposition in these areas (Wilcke et al., 2013). Other global environmental changes such as CO₂ and temperature increases may impact plant controls on soil organic matter dynamics that are mediated by microbes through priming effects, leading to the loss of nutrients and ecosystem degradation (Bernard et al., 2022). Balancing the interactions between decomposer activity and plant–microbial competition is, therefore, essential to understand nutrient dynamics along environmental gradients.

4.2 Model advancements and perspectives on plant–fungal interactions

We have presented a relatively simple and effective implementation of the root collaboration gradient, which can be included in other DVMs. A crucial step in the model development was to translate fine-root and hyphal mass into the total surface area through specific root and hyphal length (Eq. 1). Implementing the fine-root surface area calculation based on a varying length-to-mass variable within observed ranges allowed the modelled fine roots to function using observational fine-root physiological parameters for uptake kinetics (Table A1). Fine-root surface area, in particular at the highest-elevation site, is clearly too low to support nutrient uptake and vegetation growth without mycorrhiza. As the model implementation of fine-root nutrient absorption was changed from being based solely on fine-root mass to being based on fine-root surface area, it became clear that mycorrhizal associations (using observed physiological parameters) were necessary to support the observed vegetation biomass under the observed soil nutrient concentrations and nutrient demands, especially for the highest-elevation site in our study.

Our approach to link plant C expenditure to nutrient acquisition differs from other approaches since (1) we explicitly consider fungal C mass and turnover in the model, allowing for their future measurement and validation in the field; (2) plant C investment in mycorrhiza is divided into fungal growth and tissue respiration; and (3) costs of nutrient uptake are not predetermined in the model but are dynamic, varying not only between plant individuals due to fine-root traits (i.e. AMF colonization rates) but also within an individual's life history since C must be allocated to sustain fungal respiration (dependent on fungal C mass, which varies with time). In addition, we consider individual plants in our model, which compete for nutrients that result from a dynamic organic matter submodel, and consider feedbacks from plant tissue stoichiometry. Finally, LPJ-GUESS-NTD is a trait-based DVM not a PFT-based one (plant functional type, where simulated individuals from the same category have the same traits), allowing the exploration of plant functional diversity.

Modelled C costs per N or P uptake (Fig. A3) increase with elevation, which was expected due to the lower nutrient content of the soils at the highest sites. These costs are model emergent and are outcomes of the C investment processes in AMF tissue growth and respiration and of nutrient acquisition through uptake depending on plant and fungi structure and physiology. The approach used in our study is thus different from the ones followed by other models (Allen et al., 2020; Reichert et al., 2023), where the C costs of nutrient uptake or their ranges are defined a priori. The costs defined here are thus connected to other processes within the model and reflect field-measured parameters for AMF dynamics and structure, as defined in Table A1.

4.3 Model limitations and links to observations

Gaps in the current understanding of plant–fungal interactions present a challenge in modelling their processes (Hawkins et al., 2023; Makarov, 2019) and limit the projection of their effects on larger-scale processes such as the C cycle. Key processes that are poorly constrained by data are fine-root and mycorrhiza uptake of nutrients, C transfer between the plant and fungi, their tissue turnover rates and longevity, mycorrhiza C : N ratios (which impact mycorrhiza respiration and total C costs), and the N and P demand of mycorrhiza tissue. In addition to parametrization, field data are also invaluable for model evaluation. With regard to this study, our modelled total AMF mass per hectare reached values of up to 0.45 Mg C ha⁻¹, representing around 0.5 % of total plant biomass, assuming extra-radical mycelium (ERM) longevity of 80 d. Godbold et al. (2006) estimated higher values of around 1.1 Mg C ha⁻¹ for the EuroFACE sites and an ERM longevity of 9 d. These values may, however, be at the upper range of observations, as Parihar et al. (2020) found AMF biomass values between 0.054 and 0.9 Mg C ha⁻¹, which supports our choice of AMF turnover parameter in Table A1. In any case, ERM longevity and total fungal biomass are invariably linked, meaning that increases in the former will result in increases in the latter. Therefore, more measurements of both would be invaluable for models to parameterize these two important components of AMF dynamics. Such data, while in many cases being difficult and costly to measure, could be very important for future simulations and could foster collaborations between modellers and experimental and field ecologists. For example, by exploring the parameter ranges of turnover (as we have done with C investment by varying the r_max parameter), models can estimate the AMF turnover values that produce viable plant communities, which then can be tested against field estimations. This would be particularly important for large spatial- and temporal-scale estimations of increasing atmospheric CO₂ effects on soil C stocks since AMF turnover rates influence soil C (Treseder and Allen, 2000).

The new empirical relationships between SRL, fine-root diameter and AMF colonization included in the LPJ-GUESS-NTD model for this study have led to satisfactory results. However, the relationship between the AMF colonization rate and fine-root diameter as estimated in the field showed considerable variation (Fig. 2). This indicates that although the relationship between fine-root diameter and AMF colonization is considerable, there is still room for improved understanding of this topic. More thorough studies involving fine-root anatomy, physiology and traits may uncover further relationships and axes of variation, which would be invaluable to help improve vegetation models.

How realistic is our modelled relationship between AMF colonization and plant C transfer? In light of the lack of field-based empirical evidence, particularly for the type of environment we simulate, our modelling study represents a theoretical estimate of the link between mycorrhizal colonization rates and C transfer from plant to AMF as a fraction of NPP. The sensitivity analysis suggests a general agreement of Eq. (2) and observations, with an r_max=0.5. The value used for our simulations of r_max=0.5 resulted in around 30 % of NPP (Fig. 3d) or 20 % of gross primary productivity (GPP) being allocated to AMF (including AMF respiration; Fig. A3g). This value is higher than the average values suggested in the literature (Chapin III et al., 2011). Řezáčová et al. (2017) suggested that the average GPP investment of plants in AMF is less than 10 %, implying that less than 20 % of NPP is invested. However, most results exist for crops. Wild plants such as those in our study can allocate substantially more C to AMF (Hawkins et al., 2023). We conclude that the high r_max values in our results might be realistic under severe nutrient limitation and fine-root AMF colonization of around 50 %, as was observed at our highest-elevation site (Pierick et al., 2023). More field research of this maximum transfer rate (having, for example, measurements of NPP, AMF C transfer, and fine-root intra- and extra-radical colonization) would be invaluable to better constrain this important component of NPP, as well as a confirmation of whether the relationship between C transfer and the colonization rate is in fact linear, as was assumed here. A better understanding of this link would be invaluable for management practices to improve estimations of plot-level plant-to-soil C transfer and AMF biomass.

4.4 Future directions for model development, expansion of our approach and model application

Our exploration of fine-root traits and the FCG was limited to three sites in the tropical mountain forests of Ecuador. These areas were chosen due to data availability and are limited in spatial scope, but we argue that our results are relevant in general since (1) the axis of trait variation observed for our sites (Pierick et al., 2021, 2023) followed those estimated globally closely (Bergmann et al., 2020) and (2) AMF is by far the predominant mycorrhizal type phylogenetically (70 % of plant species), in land cover (55 %) and in the global NPP contribution from their plant hosts (63 %; Hawkins et al., 2023).

Other nutrient acquisition strategies may be more relevant regionally. Although at our study sites foraging strategies involving AMF are prevalent in nutrient-poor areas, in strongly P-limited areas of the Amazon lowlands, nutrient acquisition strategies such as the exudation of phosphatases or organic acids might also play a significant role (Reichert et al., 2023). The implementation of exudation C-investment strategies may be necessary to account for increased fitness in nutrient-poor environments and the allocation of C belowground, as up to 16 % of NPP was found to be allocated to the production of organic acids in a P-limited environment (Aoki et al., 2012).

The inclusion of other mycorrhizal types could also be important to correctly simulate certain environments. The second-most-important mycorrhizal type, ectomycorrhizal fungi (EMFs), is dominant in boreal and some tropical regions and is dominant in 25 % of the Earth's terrestrial ecosystems (Hawkins et al., 2023). EMFs can receive on average higher fractions of plant productivity compared to C transfers (Hawkins et al., 2023), have different pathways for nutrient uptake (Phillips et al., 2013), and may have a distinct effect on the soil and on vegetation C dynamics (Terrer et al., 2021). Since EMFs can actively extract nutrients from non-labile sources, they may allow higher plant productivity than AMFs in environments where labile sources are extremely poor and may affect simulated soil C stocks and fluxes. For instance, one modelling study showed how the presence of EMF can affect C storage (Moore et al., 2015). Regarding the FCG, EMFs are expected to follow the same SRL trends as AMFs (Bergmann et al., 2020), meaning that the main difference in future EMF model implementations might be accounting for organic nutrient acquisition kinetics. The inclusion of EMF into DVMs should consider important differences from an AMF implementation. First, the access EMF has to organic sources requires the production of specialized enzymes, which are not present in AMF. This means that total cost of nutrient acquisition borne by plants may be higher in EMF symbiosis, as confirmed empirically (Hawkins et al., 2023). Second, kinetic uptake parameters, fungal tissue C : N ratios (which in our model impact mycorrhiza respiration), turnover rates and SHL traits should differ due to distinct EMF physiology and morphology. Within the trait-varying approach, one possible implementation could be to randomize the individual's preference for AMF or EMF (or no mycorrhiza) during establishment (but not in grasslands, which are AMF exclusive). This can be an interesting prospect to test whether AMF or EMF dominance could emerge in a global model following climatic and edaphic conditions, where AMF would dominate in areas with predominantly mineral soils and EMF where soil organic matter has higher stocks or fluxes (Read, 1991). Also, we could test whether AMF and EMF effects on soil C stocks differ under increasing CO₂, as has been suggested (Terrer et al., 2021), and whether mycorrhiza are indeed one of the primary controls of the CO₂ fertilization effect (Terrer et al., 2016).

Large-scale simulations and future projections are the expected next steps for trait-based DVMs, such as the one included in this study. Global or regional maps of key vegetation traits have recently been published using several methods, such as combining machine learning and remote sensing, and represent a new way of showcasing worldwide variation in plant attributes (Moreno-Martínez et al., 2018). These maps can be used in turn to evaluate similar maps from trait-based DVMs, which can provide further ecological information that cannot be sensed from satellites, i.e. on belowground traits and processes. Regarding future projections, an obvious follow-up study to this one would be changes in fine-root traits and fungal collaboration under climate change scenarios. An analysis of data from Free Air CO₂ Enhancement (FACE) experiments showed that the mycorrhizal type may determine the extent of CO₂ fertilization effects (Terrer et al., 2018, 2021), which is a significant uncertainty in future projections of vegetation and organic C cycle changes (Hickler et al., 2015; Walker et al., 2021). More explicit treatments of mycorrhiza in DVMs and land surface models might reduce this uncertainty.

Finally, our study can provide important insights into the recent discussion on the importance of AMF for ecosystem management, in spite of its limitations. First, the relevance of AMF abundance and diversity for agriculture has been argued for (Rillig et al., 2019) and against (Ryan and Graham, 2002), particularly when the criteria for AMF benefits are defined as either yield or sustainability. In this regard, our modelling study confirms the strong and long-term influence of AMF on plant biomass and productivity for environments where nutrients are the most limiting factor but weaker effects where they are not. Second, the shift in both SRL and SLA distributions when AMF is deactivated suggests that plant communities may differ significantly when mycorrhizal symbiosis is absent, particularly in conservative assemblages. Our model is thus in line with suggestions from ecosystem restoration practices: the introduction of AMF inoculation can impact both plant growth and community composition (Lin et al., 2015). In order to better explore the AMF influence on the whole ecosystem from a model perspective, the implementation and evaluation through simulated experiments of several key processes such as soil aggregation, seedling survival, resistance to pathogens and resistance to invasive species, as well as a thorough analysis of AMF effects on soil stocks and fluxes, would be invaluable (Rillig et al., 2019).

5 Conclusions

Even though fine-root traits and AMF are strongly linked and are thought to impact ecosystem functioning, variation in fine-root traits has hardly been addressed in process-based DVMs and C-cycle models. Here, we present a parsimonious, generalized approach to implement fine-root trait diversity and the fungal collaboration gradient in a widely used DVM. With local data from a tropical montane biodiversity hotspot, we show that the model can reproduce fine-root traits, AMF colonization and aboveground vegetation features along a nutrient-limitation gradient. Model results confirm the expected crucial role of AMF for nutrient uptake and vegetation productivity at the high-elevation, strongly nutrient-limited site. The model also reveals potential linkages between fine-root traits, AMF colonization and aboveground vegetation traits. We also show how trait-based DVMs can be a powerful tool for testing ecological hypotheses concerning complex interactions in ecosystems. Future research and model development, however, should focus more on belowground traits and their interaction with mycorrhiza and aboveground traits. Belowground processes and nutrient dynamics might, in the future, become even more important because plants will be even more important due to enhanced photosynthetic activity as atmospheric CO₂ concentrations continue to rise.

Appendix A: Evaluating carbon transfer from plant to fungi

Increasing transfer of C from plants to fungi reduced nutrient limitation (Fig. A3e and f), total vegetation biomass and productivity in our simulations (Fig. 3a and b). The analysis of how the average fine-root traits (SRL and AMF colonization) changed with increasing maximum C costs, r_max (Figs. 3c and A3a), reveals a clear cost–benefit pattern in which too little or too much C transfer, here also related to fungal colonization, is detrimental to plant productivity. Too little C transfer (r_max<0.5) resulted in low total AMF biomass (Fig. A3b). For r_max values above 0.5, individuals with high C transfer and higher AMF colonization rates are outcompeted, driving down the community average AMF colonization rates (Fig. 3c). This occurs since the increasingly large AMF mass these plants would need to support offers no further nutrient limitation alleviation benefit. The resulting AMF simulated biomass from an r_max of 0.5 and above (between around 0.35 and 0.45 Mg C ha⁻¹; Fig. A3b) fell within the 0.054–0.9 Mg C ha⁻¹ of observations (Parihar et al., 2020). Since AMF biomass is an emergent property of the model (i.e. not prescribed), this fit within observations indicates that our assumptions and prescribed parameters (such as AMF turnover) produce a realistic representation of AMF functioning.

Table A1 Model parameters relevant to the new implementations of LPJ-GUESS-NTD. Variables with asterisks were measured at our study sites.

Table A2 Median absolute error measures between the scenarios (AMF-on and AMF-off) and field data for the three elevation sites.

Figure A1 Relationship between the r_max parameter and NPP transfer to mycorrhiza.

Figure A2 Complete scheme of model trait–trait relationships (trade-offs) after the inclusion of the fine-root traits and collaboration gradient. Traits in blue boxes are randomized uniformly at establishment according to a fixed range. Green boxes represent traits whose trade-offs were defined using measurements only from our study site, and red boxes were defined using data from the literature. SLA – specific leaf area; WSG – wood specific gravity; SRL – specific root length.

Figure A3 Sensitivity analysis of the parameter r_max, the maximum C allocation to mycorrhiza growth. Each point represents 200-year averages using a particular r_max. Whiskers indicate ± SD from the 200 years of each run. Red lines – 1000 m. Green lines – 2000 m. Blue lines – 3000 m.

Figure A4 Simulated proportion of AMF vs. total vegetation biomass for the three elevation sites.

Figure A5 Boxplots of (a) specific leaf area (SLA) and (b) specific root length (SRL) for the two simulation scenarios (AMF-on and AMF-off) and field measurements.