We collected the root trait data of plants in China from the global Fine-Root Ecology Database (FRED¹) (Iversen et al., 2017) and the published literature (listed in Note 1 in Supplementary Material). We focused on first- and second-order roots that are defined by branching order, where the most distal roots are numbered as first order and where second-order roots begin at the junction of two first-order roots (Pregitzer et al., 2002). The first- and second-order roots generally belong to absorptive roots (Guo et al., 2008; Freschet and Roumet, 2017), which have high physiological activity and resource acquisition efficiency (McCormack et al., 2015). Thus, the studies that used diameter-based fine roots (e.g., ≤1 mm, ≤2 mm, and ≤0.5 mm) and absorptive roots without a clear root branching order were excluded. To ensure data quality and homogeneity, root trait data were obtained according to the following criteria: (1) studies conducted in fields were included and those conducted in croplands, aquatic ecosystems, greenhouses, and laboratories were excluded, in order to minimize the effects of management disturbance; (2) root samples were collected from mature and healthy plant individuals to minimize the effects of ontogeny (Alvarez-Flores et al., 2014); (3) root samples were collected from live roots and data from dead roots were excluded to reduce the confounding effects of root vitality; and (4) root samples were collected from native species and non-native species were excluded. When a species occurred at multiple sampling sites, all site–species trait values were recorded.

Species name and taxonomic nomenclature were standardized and corrected according to the Plant List² using the “plantlist” package. A total of 407 site × species observations of 306 species from 72 families and 174 genera were collected. The species were classified into seven phylogenetic clades according to APG IV classification (APG, 2016): gymnosperms, chloranthales, monocots, magnoliids, basal eudicots, asterids, and rosids. Basal eudicots include species that are not included in the clades of rosids and asterids. The information on plant growth form (herbaceous and woody) was obtained from the original literature and the TRY–Categorical Traits Dataset³ (Kattge et al., 2011). Our dataset included 53 herbaceous and 253 woody species, respectively.

The information regarding the type of mycorrhizal association was collected and corrected according to the original literature and previously published databases (Wang and Qiu, 2006; Hempel et al., 2013; Phillips et al., 2013; Valverde-Barrantes et al., 2017; Guerrero-Ramírez et al., 2021). In this study, mycorrhizal types were classified into five groups based on the Global root traits (GRooT) database (Guerrero-Ramírez et al., 2021): AM (221), ECM (68), ericoid mycorrhizal (ERM, 2), AM-ECM (8), and AM-non-mycorrhizal (AM-NM, 7).

The species included in this study accounted for almost all major biomes in China, including forests, grasslands, and alpine tundra. The forests were divided into four biome types: tropical forest (latitude < 23.5°N), subtropical forest (latitude 23.5°–34°N), temperate forest (latitude 34°–50°N), and boreal forest (latitude > 50°N). Species with multiple biome entries were categorized into the biome in which they had the most observations. In this dataset, the tropical forest biome contained 31 species, subtropical forest 144 species, temperate forest 77 species, boreal forest 10 species, grassland 15 species, and the alpine tundra 29 species.

The phylogenetic tree was constructed based on the comprehensive angiosperm species-level phylogeny of Zanne et al. (2014) updated by Qian and Jin (2016). The phylogenetic tree was constructed using the “S.PhyloMaker” function in the “phytools” package (Revell, 2012). All 306 species were included in this phylogenetic tree.

Before data analyses, all root trait data were log₁₀-transformed to meet the assumption of normality. In this study, all data analyses were conducted at the species level, species mean trait values were thus used. We only focus on interspecific trait variation in our analysis.

Blomberg’s K (Blomberg et al., 2003) and Pagel’s λ (Pagel, 1999) values were calculated to assess the strength of phylogenetic signals of root traits using the “phytools” package (Revell, 2012). Blomberg’s K values > 1 suggest higher phylogenetic conservatism than expected owing to the Brownian motion, and values <1 suggest weaker phylogenetic conservatism. Pagel’s λ values close to 0 indicate phylogenetic independence, and values closer to 1 indicate that the trait distribution perfectly complies with Brownian motion (Münkemüller et al., 2012).

Correlations among root traits were examined using Pearson’s correlation analysis using the “Hmisc” package and phylogenetic independent contrasts (PICs) that account for the phylogenetic relatedness among species using the “picante” package (Felsenstein, 1985). Linear and non-linear regressions were also performed to evaluate the pairwise relationships among root traits. These analyses were conducted for all species and within different plant growth forms, mycorrhizal types, and biomes.

Considering the phylogenetic relatedness among species, phylogenetic principal component analysis (pPCA) was performed to identify the dominant dimensions of root trait covariation using the “phytools” package (Revell, 2012). pPCA is an extended method of PCA that is a powerful multivariate analysis technique and can summarize a set of data on correlated variables with a few composite, uncorrelated principal components (James and McCulloch, 1990). pPCA has been widely used in the determination of independent axes of functional specialization (e.g., Wang et al., 2018a; Liu C. et al., 2019; Bergmann et al., 2020; Weigelt et al., 2021). The Kaiser’s eigenvalue greater than 1 rule was used to determine the intrinsic dimensionality of root trait covariation, this is, eigenvalues of the principal components greater than 1 were considered significant (Kaiser, 1958; Tabachnik and Fidell, 1996). This rule was used in previous studies that examined the main axes of plant trait covariation (e.g., Comas and Eissenstat, 2009; Chen et al., 2013; Jager et al., 2015; Kramer-Walter et al., 2016). To test whether the pattern of root trait covariation for all species was similar to those for different plant growth forms, mycorrhizal types, and biomes, the pPCA analysis was further repeated for subsets of plant growth forms (herbaceous and woody), mycorrhizal types (AM and ECM), and biomes (tropical, subtropical, temperate forest types, and alpine tundra). ERM, AM-ECM, AM-NM, boreal forest, and grassland did not perform the pPCA analysis as their sample sizes were small (n ≤ 15).

To assess the segregation of plant species by their phylogenetic clades, plant growth forms, mycorrhizal types, and biomes along the dominant axes of root trait covariation, one-way ANOVA was used on the scores of species on the pPCA axes and the significance of pairwise differences was tested using Tukey’s HSD post hoc test.

The variance partitioning analysis was used to determine the relative contributions of plant growth form, mycorrhizal type, and biome according to the pPCA scores of the dominant principal components of root trait covariation for all species. The significances of each factor and their interactions were tested using 999 permutations. The analysis was followed by Legendre and Legendre (2012) using the “vegan” package.

All statistical analyses were conducted using R version 4.0.3 (R Core Team, 2020).

For all species, RD was negatively associated with SRL (Figure 1A and Table 1). RD was not significantly associated with RTD, RNC, RCC, and RCN (Figures 1B–E), but it was negatively associated with RTD after considering the phylogenetic information (Table 1). Regardless of the phylogenetic relatedness among species being accounted for, SRL was always negatively and positively associated with RTD and RNC, respectively (Figures 1F,J). RTD, SRL, and RNC were significantly related to RCN (Figures 1I,L,O). However, RCC did not show a significant relationship with RCN (Figure 1N) but it showed significant relationships with RTD, SRL, and RNC (Figures 1H,K,M).

All root traits (except RCC) exhibited significant phylogenetic signals (Table 2). The pPCA results showed that eigenvalues of the first three axes were greater than 1 (Table 2), indicating that the covariation in root traits was represented by three independent dimensions. The first three primary axes accounted for 85.3% of the total variation of root traits. The first PCA axis (PC1) accounted for 39.0% of the total variation and was mainly related to SRL, RTD, RNC, and RCN (Figure 2A). With increasing scores of the PC1 axis, SRL and RNC increased while RTD and RCN decreased. The second PCA axis (PC2) explained 26.1% of the total variation and showed a negative relationship between RD and SRL. The third PCA axis (PC3) accounted for an additional 20.2% of the total variation and was primarily associated with RCC. The species distribution within the trait space showed that gymnosperms had the lowest PC1 scores with higher RTD and RCN values. Magnoliids species had the highest PC2 values with higher RD and lower SRL values whereas monocots and eudicots species had lower PC2 values with lower RD and higher SRL values (Figure 2A and Supplementary Figure 1).

In terms of plant growth form, RD was negatively related to RTD in woody species, but this relationship was not observed in herbaceous species (Figure 1B and Supplementary Figures 2, 3). Compared to the herbaceous group, the woody group exhibited weaker relationships between SRL and RTD and RNC (Figures 1F,G). The ANOVA results showed that herbaceous and woody species showed significant differences in the scores of the first three axes of root trait covariation (i.e., PC1–PC3) (P < 0.05, Figures 3A,D,G). Herbaceous species had higher scores of the PC1 (e.g., higher SRL and RNC) and the lower scores of the PC2 (e.g., lower RD and higher SRL) compared with woody species. Woody species showed the consistent root trait covariation with the pattern of root trait covariation for all species, representing by three independent dimensions (Supplementary Table 1 and Supplementary Figure 4B). For herbaceous species, the first three axes explained 87.6% of the total variation in root traits (Supplementary Table 1 and Supplementary Figure 4A). Within the first axes, the main difference between herbaceous and woody species was that RD was strongly loaded on PC1 for herbaceous species. In addition, the second dimension was represented by RD, SRL, and RNC; the third dimension was dominated by RCC and RCN for herbaceous species.

The relationships among root traits differed across mycorrhizal types (Supplementary Figures 2, 3, 5). RD was negatively associated with RTD in AM species but was positively associated with RTD in ECM species (Supplementary Figure 5B). RD was negatively related to RNC in ECM species but showed no correlation in AM species (Supplementary Figure 5C). Notably, most correlations among root traits were stronger in ECM species than those in AM species. The ANOVA results showed that AM and ECM species showed significant differences in the scores of the PC1 but not for the PC2 and PC3 (Figures 3B,E,H). AM species had higher PC1 scores (i.e., higher SRL and RNC values) compared with ECM species. The AM group also showed three independent dimensions of root trait covariation, which was consistent with the pattern for all species (Supplementary Table 1 and Supplementary Figure 4C). ECM species also showed three independent axes of root trait variation. The PC1 axis was dominated by RD, RTD, RCN, SRL, and RNC, the PC2 axis was dominated by RCC, and the PC3 axis was dominated by RD, SRL, and RCN (Supplementary Table 1 and Supplementary Figure 4D).

Moreover, the correlations among root traits differed across biomes (Supplementary Figures 2, 3, 6). RD was negatively correlated with RTD in the tropical and subtropical forests, but not in the other biomes. SRL and RNC values decreased with increasing RTD only in the temperate forest. In addition, there were significant differences in the scores of PC1, PC2, and PC3 among biomes (P < 0.05, Figures 3C,F,I) and had different trait syndromes among different biomes (Figure 2D). Species in the temperate forest had the highest PC1 scores with higher SRL and RNC values, whereas species in the alpine tundra had the lowest PC1 scores with higher RCN and RTD values. Species in the tropical and subtropical forests showed higher variations on the PC2 axis compared with other biomes. The dominant dimensions of the covariation in root traits differed greatly across biomes. In the tropical and subtropical forests, the first three axes explained 90.0 and 86.5% of the total variation in root traits, respectively (Supplementary Table 2 and Supplementary Figures 4E,F). In these two forests, the PC1 axis was heavily loaded by RD, RTD, RNC, and RCN, and the PC2 axis was strongly represented by the negative relationships between RD and SRL. In addition, two independent root trait dimensions – representing root morphological traits (i.e., RD, RTD, and SRL) and nutrient traits (i.e., RNC and RCN) – were detected in the alpine tundra (Supplementary Table 2 and Supplementary Figure 4H).

Plant growth form, mycorrhizal type, and biome together accounted for 52.0, 14.1, and 13.8% of the total variation of PC1, PC2, and PC3, respectively (Figure 4). The PC1 axis was mostly explained by biome alone (33.0%), followed by the interactive effects of biome and mycorrhizal type (14.8%). The PC2 axis was mostly explained by mycorrhizal type alone (7.5%) and biome alone (5.6%), and the PC3 axis was mainly influenced by biome alone (11.8%), followed by plant growth form alone (1.7%).

Our study demonstrated three independent dimensions of the covariation among the six root traits in China. The first dimension was dominated by SRL, RTD, RNC, and RCN that were closely interrelated, which is in accord with the RES expectations. This result is not in agreement with the previous studies reporting that RD and SRL were orthogonal to RTD and RNC (Kramer-Walter et al., 2016; Wang et al., 2018a; McCormack and Iversen, 2019; Bergmann et al., 2020). Within the first axis, in one end, species had high RTD, high RCN, low SRL, and low RNC associated with a resource-conservation strategy, whereas in the other end species had high SRL, high RNC, low RTD, and low RCN associated with a resource-acquisitive strategy (Roumet et al., 2016; de la Riva et al., 2018). As the expected from the RES, we found that SRL was negatively associated with RTD and positively related to RNC, suggesting a trade-off between resource acquisition and construction costs of roots (Makita et al., 2012; Li et al., 2019). According to the RES expectations, fine roots with a higher SRL should be associated with a higher uptake activity and metabolic rate (i.e., higher RNC) and shorter lifespan (i.e., lower RTD) to maximize resource acquisition per investment (Eissenstat et al., 2000; Weemstra et al., 2016; de la Riva et al., 2021b). RD was not related to RNC, RTD, and RCN that strongly loaded on the first axis, thus RD was decoupled from the first axis. In addition, the second axis was dominated by a trade-off between RD and SRL, which did not support the previous studies reporting that the variation of root traits was mostly explained by RD and SRL (Wang et al., 2018a; McCormack and Iversen, 2019; Bergmann et al., 2020). Recent studies have reported that RD and SRL were positively and negatively related to the percentage of mycorrhizal colonization, respectively (Ma et al., 2018; Bergmann et al., 2020), suggesting that roots can enhance resource uptake from soils by constructing either thin-diameter roots with higher SRLs (i.e., “do-it-yourself” strategy), or in contrast thick-diameter roots via reliance on mycorrhizal associations (i.e., “outsourcing” strategy) (McCormack and Iversen, 2019; Bergmann et al., 2020). However, roots with a larger mycorrhizal colonization may increase the resource uptake capacity under resource-limited conditions without necessarily implying an acquisitive strategy (Navarro-Fernández et al., 2016). Therefore, it is needed to uncover how this root–mycorrhizal collaboration gradient links to the differences in resource uptake capacity for species with fast- and slow-traits (de la Riva et al., 2021a). Furthermore, we found that RCC formed an additional axis that was decoupled from the first two dimensions. This lack of correlation may be attributable to two reasons. First, root traits that were related to the first two axes displayed significant phylogenetic conservatism, but RCC did not, suggesting that these trait associations are not tightly coupled. Second, in theory, RCC, RNC, and RCN are mathematically interdependent. However, we observed that RNC was strongly associated with RCN, whereas RCC was not. This indicated that RNC, rather than RCC, was more important in the trade-offs between C investment and resource uptake in roots (An et al., 2021). Our results is not consistent with the global and regional studies (Kramer-Walter et al., 2016; Wang et al., 2018a; McCormack and Iversen, 2019; Bergmann et al., 2020), we found that the RES explained the most variation of root traits and the trade-off between RD and SRL loaded on the second axis in China, suggesting that root trait covariation may not be generalized from global-scale and other regional analyses.

Furthermore, our results revealed that the pattern of variation in root traits was closely linked to the phylogenetic structure; that is, species within different phylogenetic clades occupied different locations in the trait space. Gymnosperms generally dominated boreal and subalpine forests that are mainly constrained by low temperatures, as thick roots with high RTD values tend to have high physical robustness to cope with cold environments (Simpson et al., 2016; Wang et al., 2018b; Yahara et al., 2019). Early diverged magnoliids are generally associated with phosphorus-limited tropical and subtropical soils (Ma et al., 2018), resulting in high RD and RNC values (Figure 5). Thick roots maximize the cortex area, supporting more AM colonization as a complementary strategy for nutrient foraging, leading to high RNC (Brundrett, 2002; Comas et al., 2012; Kong et al., 2019). Moreover, recently diverged eudicots species were more concentrated toward lower RD and higher SRL values. This can allow roots to increase their surface area and explore larger soil volumes per unit of C investment – that is, a “do-it-yourself” strategy (Comas et al., 2012; Valverde-Barrantes et al., 2017; Bergmann et al., 2020).

The uncertainties of this study lie in four aspects due to the restrictions of data coverage. First, our study focused on commonly studied mycorrhizal statuses (i.e., AM and ECM), some other mycorrhizal types such as ERM, NM, and dual associations were not considered in our analyses owing to the limited number of root samples. The type of mycorrhizal partnership is an important driver of the variation in root traits and their functions (Valverde-Barrantes et al., 2017; Gao et al., 2021). More detailed studies of mycorrhizal status are needed to elucidate how the evolution of root traits and their mycorrhizal associations affect belowground processes in seed plants. Second, the uneven distribution of data across biomes is likely to impede our understanding of the effect of biomes on root resource uptake. In particular, we had a limited number of species in the boreal forest and grassland, which may explain the weak or insignificant relationships among root traits in these biomes. Further work with a wide representation of species from multiple biomes is crucial to improve our understanding of the role of biomes in resource acquisition strategies in fine-root systems. Third, our analysis focused on the interspecific variation in root traits, thus resulting in some uncertainties regarding the importance of intraspecific variation and plasticity in trait-based studies (Jung et al., 2010; Weemstra et al., 2020). Future studies that incorporate the intraspecific trait plasticity will help to elucidate the trade-offs among root traits related to belowground resource acquisition strategies (Isaac et al., 2017). Finally, we considered six root morphological and chemical traits in our study, however, the root trait covariation may be dependent on the trait variables studied. de la Riva et al. (2021a) demonstrated that SRA (i.e., mass-normalized) was more tightly correlated than SRL with the dry matter content and chemical composition of both roots and leaves along the economics spectrum for Mediterranean vegetation. Therefore, it is important to measure and integrate the root trait data based on the standardized approaches, which can be analogous to traits related to the LES, including root morphology (e.g., SRA and root dry matter content), root chemistry (e.g., RNC and lignin), root anatomy (e.g., cortical area and stele area), percentage of mycorrhizal colonization, and root function (e.g., respiration, decomposition, and resources uptake) (Laliberte, 2017). This would allow us to develop a more general integrated framework related to the trade-offs between root structure and function and their underlying mechanisms.