We used the MG-RAST server, because the metagenomes data deposited in the server are open and accessible publicly. In addition, the server contains functional analysis and can directly export functional and taxonomic composition matrices from metagenomic data. Instead of directly obtaining available shotgun metagenomic data from MG-RAST server, we only selected soil metagenomes that have been published in peer-reviewed journals to ensure the quality and completeness of the metagenomic data, because we believed that the peer-review process in publication papers can confer a tighter control over the quality of the metagenomic data that have been included in our study. We searched peer-reviewed publications from 2012 to 2018 using the Web of Science database for our literature survey on December 10, 2018. Search keyword combinations included “soil, metagenome, MG-RAST,” “soil, shotgun sequencing, MG-RAST,” and “soil, Illumina, MG-RAST.” The following criteria were applied: (1) studies that analyzed and/or directly deposited soil metagenomes (generated using shotgun sequencing that involves randomly breaking the DNA into a collection of small fragments that are sequenced individually to obtain reads) to the MG-RAST database; (2) soil metagenomes in the MG-RAST database that could be publicly accessible. Based on the Study ID and/or the MG-RAST ID reported in the publications, we further extracted data matrix of soil metagenomes from the MG-RAST database (Chen et al., 2021a,b). In total, this study included 845 soil metagenomes across 17 climate zones around the world, extracted from 56 MG-RAST studies published in 51 peer-reviewed papers (Figure 1). Details of each soil metagenome extracted from publications and the MG-RAST database was given in Supplementary Table 1. For example, Study ID, MG-RAST ID, sample name, bp (base pair), and reads were obtained from metadata of each soil metagenome in MG-RAST database. Generally, the ranges (mean ± standard error) of bp (base pair), raw reads, SEED Subsystems hits, and RefSeq hits were 2.41×10⁹ ± 1.34×10⁸, 2.05×10⁷ ± 1.34E×10⁶, 2.01×10⁶ ± 1.22×10⁵, and 5.11×10⁶ ± 2.96×10⁵, respectively.

For each soil metagenome, using the “Analysis” function of the MG-RAST server, we loaded the SEED Subsystems database, which contain functional information structured into three levels, with Level 1 being the most general class while the Level 3 being the most specific class (Function) (Overbeek et al., 2013), as the metagenomic functional potentials, and loaded from the RefSeq database taxonomic information at genus, family, order, class, and phylum levels (Taxonomy) (Tatusova et al., 2013) as taxonomic compositions. Although RefSeq database had taxonomic resolution to species, it was not recommended by the server for most application using shotgun sequence data to infer taxonomic information below the genus level, so we kept the lowest taxonomic information at the genus level in this study.

The analyses were performed using default settings (maximum e-value cutoff was 1e^–5, minimum identity cutoff was 60%, and minimum alignment length was 50) (Meyer et al., 2008). We further merged data matrix of each functional or taxonomic group extracted from different studies together to generate new datasets of microbial profiles of functional potentials annotated by the SEED Subsystems database and taxonomic compositions annotated by the RefSeq database. We chose Subsystems database for functional grouping rather than the KEGG Orthology (KO) (Kanehisa et al., 2015), Clusters of Orthologous Groups (COG) (Galperin et al., 2014), and Non-supervised Orthologous Groups (NOG) (Huerta-Cepas et al., 2015) databases because Subsystems had more than 21 classification at Level 1, comparable to RefSeq taxonomy at phylum levels, allowing us to conduct comparison between taxonomic and functional annotations. We chose RefSeq database rather than the traditional ribosomal RNA databases, such as RDP (Ribosomal Database Project) (Cole et al., 2008), Greengenes (Desantis et al., 2006), or Silva LSU/SSU (Pruesse et al., 2007) databases, because taxonomic hits in the RefSeq database were over 1,000-fold higher than the rRNA databases, rendering the resolution comparable to functional hits for comparison (Supplementary Table 1). Since the assembled and unassembled sequences may result in potential bias in the diversity of taxonomy and function annotations, all the publication included in this study have deposited or submitted their unassembled raw sequences in MG-RAST server. Even though some studies reported their own processing data based on de novo assembled contigs in their publications, as long as they used MG-RAST server to deposit their unassembled raw sequences, they were all included in this study to increase the coverage of our datasets.

The metagenomic data deposited in the MG-RAST server were either assembled or unassembled sequences, both of which were included in this study to increase the coverage of our datasets.

The geographic coordinates of latitudes (LAT) and longitudes (LONG) of each soil metagenome were directly obtained from publications. Based on LAT and LONG, climate data of mean annual temperature (MAT) and precipitation (MAP) of the study sites for each soil metagenome were extracted from the WorldClim version 2 dataset (Fick and Hijmans, 2017) using the R package “raster.” Because not all study provided edaphic property data, we obtained soil properties, including C (soil total carbon content, g C kg-1 soil), pH (soil pH measured in a 1:5 w/v suspension of deionized water, no unit), sand (percentage of soil particles ranging from 2.0 to 0.05 mm in diameter, %), silt (percentage of soil particles ranging from 0.05 to 0.002 mm in diameter, %), and clay (percentage of soil particles smaller than 0.002 mm in diameter, %), from the SoilGrids database¹ (Hengl et al., 2017). Typically, the model data for climate and edaphic properties were used to compare their correlation with functional and taxonomic diversities of each soil metagenome tested in our study. To examine how microbial taxonomic and functional diversities differ globally, soil metagenomic data was classified into 17 climate zones based on the main classification of Koeppen-Geiger Climatic Zones (Kottek et al., 2006) using the R package “kgc.”

To minimize bias caused by different sequencing depths and read lengths among studies, we standardize the hit of each functional or taxonomic category in the data matrix to relative abundance through dividing by total hits to calculate beta-diversity of all tested samples. We found the functional and taxonomic diversities among different samples did not depend on the sequencing depths and read lengths of the tested samples. Even if the size of the dataset might affect the direct comparison among soil metagenomes, this study mainly focused on assessment of the relationships between functional and taxonomic diversities as well as climate and soil distances. Thus, the size of the dataset should have similar overall effects on both functional and taxonomic compositions in the same soil metagenomes, which would not lead to significant bias as a result. Based on the relative functional abundance at level 3 (Function, 1,135 categories) and the taxonomic abundance at genus level (Taxonomy, 1,113 categories), we calculated Bray-Curtis similarity and constructed triangular pairwise Bray-Curtis similarity matrices, which were further transformed to lists of pairwise Bray-Curtis similarities ordered by sample names in PRIMER 7 (Plymouth Routines in Multivariate Ecological Research Statistical Software, v7.0.13, PRIMER-E Ltd., United States) (Clarke and Gorley, 2015). Pearson’s correlation between these Bray-Curtis matrices of functional and taxonomic similarity were constructed to examine the relationship between functional and taxonomic diversities. Based on the MAT and MAP data (Climate) and soil properties of C, pH, sand, silt, and clay (Soil), triangular pairwise Euclidean distance was calculated to construct triangular pairwise Euclidean distance matrix, which were further transformed to a list of pairwise Euclidean distance ordered by sample names in PRIMER 7. To evaluate how functional and taxonomic similarity decays with climate and soil distance, Pearson’s correlations were constructed between the pairwise Euclidean distance of climate and soil data and Bray-Curtis similarity of function and taxonomy. To analyze functional and taxonomic composition structures of soil metagenomes annotated in Subsystems database at Level 3 (Function) and in RefSeq database at genus level (Taxonomy), PCoA (principal coordinates analysis), and PERMANOVA (Permutational multivariate analysis of variance) were conducted using the triangular pairwise Bray-Curtis similarity matrices in PRIMER 7. To evaluate how latitude affects climate, soil properties, and functional and taxonomic compositions, principal coordinate (PCoA) 1 scores of Euclidean distance of climate and soil data, and Bray-Curtis similarity of function and taxonomy were used to construct Pearson’s correlations with the absolute latitude of soil metagenome locations grouped by 17 climate zones. To assess the associations of functional and taxonomic diversities with climate and soil properties, DistLM (distance-based linear model) (Anderson, 2004) was conducted in PRIMER 7 using forward procedure and adjusted R². Pearson’s correlations were further constructed to assess the relationships of specific functional and taxonomic compositions with climate and soil properties. One-way analysis of variance (ANOVA) with P-values adjusted by Bonferroni-correction for multiple comparisons was conducted using SPSS 22.0 software (Chicago, IL, United States) to evaluate the statistical differences in the relative abundance of dominant functional and taxonomic compositions (mean > 1%) among climate zones after the normality of residues and homogeneity of variance were checked using Shapiro-Wilk and Levene test, respectively. The significance level was set at α = 0.05 unless otherwise stated. Heat maps were constructed using HeatMapper (Babicki et al., 2016).

To examine the potential interactions of microbial functional and taxonomic compositions across the globe, co-occurrence network analysis was performed using the Molecular Ecological Network Analyses Pipeline² (Zhou et al., 2011; Deng et al., 2012). The data matrix of standardized relative abundance, multiplied by 10⁶ to satisfy the requirements of the pipeline, was uploaded to construct the network with default settings, including (1) only keeping the species present in more than a half of all samples; (2) only filling with 0.01 in blanks with paired valid values; (3) taking logarithm with recommended similarity matrix of Pearson’s correlation coefficient; and (4) calculation ordered to decrease the cutoff from top using regress poisson distribution only. A default cutoff value (similarity threshold, S_t) for the similarity matrix was generated to assign a link between the pair of species. Then, the global network properties, the individual nodes’ centrality, and the module separation and modularity were analyzed based on default settings using greedy modularity optimization. Network files were exported and visualized using Cytoscape software (Shannon et al., 2003). The scatter plots of within-module connectivity (Zi) and among-module connectivity (Pi) were constructed to show the network node distribution of module-based topological roles of functional and taxonomic compositions. The threshold values of Zi and Pi for categorizing were 2.5 and 0.62, respectively (Guimerà and Nunes Amaral, 2005; Olesen et al., 2006; Guimerà et al., 2007).

In total, we analyzed 845 soil metagenomes that resulted in 356,090 pairwise comparisons of Bray-curtis similarity in functional (Subsystems Level 3) and taxonomic (RefSeq genus) diversities. Globally, soil taxonomic compositions were more diversely changed than functional potentials (Figures 2A,B). Pearson’s correlation showed that pairwise similarity of function was positively correlated to that of taxonomy as indicated by the linear regression (r² = 0.2831) (Figure 2A), which had higher coefficient than the logarithmic regression (r² = 0.2752). Paired-samples t-test showed that functional potentials had higher similarity (44–99%) than taxonomy (9–100%) (t = 67.5, P < 0.0001), regardless of the association between functional and taxonomic diversities revealed by Pearson’s correlations (Figure 2A). Functional diversity was more conserved at both classification Level 1 and 2, with higher similarity than taxonomy at all five phylogenetic levels examined (Figure 2B).

Several environmental metagenomic studies have reported linear relationships between functional and taxonomic diversities (Fierer et al., 2012, 2013; Leff et al., 2015), indicating a certain degree of dependency of microbial functional potentials on taxonomic compositions. This significant association, however, does not necessarily imply a lower level of microbial functional redundancy, especially in the terrestrial ecosystems. On the contrary, lower levels of diversity in microbial functional potentials, indicated by the lower dissimilarity in metagenomic composition of functional genes than in taxonomic composition, still supports that the relative abundance of microbial functions based on metagenomic genes are more stable than taxonomy in response to ecological and environmental perturbations. It should be noted that the taxonomic database of RefSeq in our analysis was more diverse than traditional ribosomal RNA databases used in the above mentioned studies. In agreement with our results, some prior studies (Fierer et al., 2012, 2013; Leff et al., 2015) also showed a lower extent of metagenomic diversity in functional potentials compared to taxonomic compositions despite their linear correlations. Our previous study has found a significantly positive correlation between pairwise similarity of function and taxonomy in water metagenomes (r² = 0.3344, P < 0.0001), though taxonomy had higher global dissimilarity than function (Chen et al., 2021a), suggesting different degrees of microbial function redundancy between soil and water biomes.

Absolute latitude of each soil metagenome was used to correlate with climate and soil property data as well as the taxonomic and functional beta-diversities using Pearson’s correlation. Taxonomic compositions showed significant associations with absolute latitudes (P < 0.0001), consistent with the latitudinal response of climate (MAT and MAP) and edaphic (SOC, pH, and soil texture) properties, while functional potentials showed no significant latitudinal gradient (P = 0.7924) (Figure 3). Thus, geographical/climate gradient could influence soil physicochemical characteristics and microbial taxonomic diversity, but due to functional redundancy, the relative abundances of microbial functional potentials were not significantly affected.

To examine how microbial taxonomic and functional diversities differ globally, soil metagenomic data was were compared among the 17 climate zones. However, some climate zones had very few sites, which may affect the analysis of functional and taxonomic compositions in different climate zones. Thus, we combined the climate zones that have numbers of samples less than 10 based on their upper level classifications, leading to 12 climate zones. Across these climate zones, microbial taxonomic compositions were more dissimilar than functional potentials (Figure 4A). The relative abundance of functional potentials remained stable across climate zones (Figure 4B). In contrast, taxonomic compositions significantly varied across climate zones (Figure 4B). Climate and edaphic properties showed a more significant influence on microbial taxonomic compositions than on functional potentials (Supplementary Table 2). DistLM analysis indicated that the cumulative influences of climate and soil characteristics were also greater on taxonomic compositions (16.3%) than on functional potentials (13.5%) (Supplementary Table 3). The decoupling of soil microbial function and taxonomy has been evaluated in our previous study that soil microbial diversities of taxonomy and function responded distinctly to five levels of simulated reduction of microbial species, leading to stable soil microbial functional structure associated with significant microbial species richness decline (Chen et al., 2021b).

Climate data and soil properties were significantly associated with the relative abundance of dominant functional and taxonomic compositions, as revealed by Pearson’s correlations (Figure 5). For functional potentials, nutrient and energy metabolisms were possitively associated with MAT and MAP, but negatively correlated to soil pH, suggesting an enhanced microbial activity with increasing temperature and precipitation, but limited microbial metabolisms in low pH environments. For taxonomic compositions, soil pH was negatively associated with the abundance/occurrence of Acidobacteria and Proteobacteria, which also indicated opposite trends against the effects of MAT and MAP. It has been found that many Acidobacteria are acidophilic, potentially due to their cell specialization and enzyme stability, and thus more adapted to low pH conditions (Kielak et al., 2016). Proteobacteria, a major phylum of bacteria commonly found in the soil, have a wide variety of metabolism types (Spain et al., 2009), particularly some Alphaproteobacteria that can grow at very low levels of nutrients, perhaps making them dominate in acidic soil environments (Feng et al., 2014). However, it has been shown that Proteobacteria increased in relative abundance with increasing soil pH, especially Alphaproteobacteria that were most abundant in the high pH soils (Rousk et al., 2010). Such contradictory conclusions indicate that additional studies are surely necessary to elucidate the interactions of Proteobacteria and soil pH conditions. In addition, both functional and taxonomic diversities gradually decayed with increasing distance of climate and edaphic properties (Supplementary Figure 1).

However, microbial community structure and functional potentials were solely based on the relative abundance of metagenomic genes, which did not directly involve any microbial functional measurement. In our study, functional potentials relied on the presence of functional genes, but not the interaction of the tested gene pathways with microbial activities and the environment. Thus, it was not concluded that microbial functions were the same across all gradients, but rather many individual metagenomic functional genes, which may or may not be active in a given environment, could be observed in the soil samples across the globe. Latitudinal diversity gradient, a decline of biodiversity with latitude (Hillebrand, 2004), has been found in alpha-diversity of microbes in the soil (Zhou et al., 2016) and marine (Fuhrman et al., 2008; Huerta-Cepas et al., 2015) environments due to temperature and pH gradient. Additionally, we further showed this latitudinal correlation of beta-diversity in microbial taxonomic compositions but not in functional potentials in terrestrial ecosystems.

The global taxonomic co-occurrence network exhibits higher degrees of modularity than the functional network (Figure 6). Almost all nine major modules of taxonomic network were comprised of genera from the same phylum/class, such as Betaproteobacteria in module #1, Bacteroidetes in module #3, Cyanobacteria in module #5, Gammaproteobacteria in module #6 and 8, Actinobacteria in module #7, and Alphaproteobacteria in module #9. Thus, the interaction of bacteria was conserved in taxonomy, especially at the phylum or class levels. Although the numbers of functional categories (1,036 Level 3) were greater than those of taxonomic groups (885 genus), the functional co-occurrence network had significantly lower degree of modularity, with less total nodes and links, indicating less interactions and connectivity than the taxonomic network, which also potentially infers a certain degree of soil microbial functional redundancy from the perspective of co-occurrence networks. Also, other key network indexes of function and taxonomy, such as the average clustering coefficient, geodesic efficiency, and harmonic geodesic distance were significantly different from those of each other and the corresponding random networks with the same network size and average number of links (Supplementary Table 4).

Despites there are numerous connections among modules of functional networks and within each module, five major modules were identified, and functional nodes could be further classified into multiple categories. Functional network had 805 positive and 342 negative links, suggesting that functional potentials had both facilitative and inhibitive interactions, while taxonomic network only had positive interactions, indicating that the response of taxonomic compositions are all cooperative (Faust and Raes, 2012). Thus, soil microbes from the same groups tended to respond in a similar way, while distinct microorganisms rarely interact with each other, suggesting that the interaction of taxonomic compositions are taxa-specific in niches (Faust and Raes, 2012).

Examining closely into the co-occurrence network of functional composition, 21 nodes were identified as either module hubs, which were connected to many nodes within their own modules (high Zi, low Pi), or connectors linking several modules (low Zi, high Pi), as indicated by the Zi-Pi plot (Olesen et al., 2007; Deng et al., 2012; Figure 7). Two functional nodes, namely ammonia assimilation and lactate fermentation, were grouped to network hubs as “supergeneralists,” acting as both module hubs and connectors interacting with the rest of the functional network. Significantly less nodes could be identified as module hubs in the co-occurrence network of taxonomic composition, indicating less interactions among different modules. This is expected given that module was comprised of genera that were mainly from the same phylum/class.

Microbial taxonomic composition and diversity plays a critical role in maintaining ecosystem function, and often shows greater sensitivity to environmental disturbance than microbially-mediated functional potentials (Allison and Martiny, 2008). However, taxonomic information alone provides limited utility in predicting biogeochemical transformations or changes of biogeochemical process rates at ecosystem level (Hall et al., 2018). Because of the decoupling of microbial taxonomic composition and functional potentials as we have discussed in the above sections, understanding which members of the microbial community can provide causative linkages would be critical for accurate prediction of microbial metabolic activity and flexibility across space and time. In this study, by analyzing and generalizing microbial taxonomic compositions and functional potentials based on metagenomic analysis of functional genes, we provide a strong evidence for soil microbial functional redundancy across the globe without any direct measurement of microbial functions. The environmental conditions likely determine the microbial taxonomic composition, and microbial phylotypes sharing similar habitat preferences tend to co-occur (Delgado-Baquerizo et al., 2018; Ramírez-Flandes et al., 2019). Though microbial metabolic functions can be strongly coupled to elemental cycles and certain environmental factors, the decoupling of microbial taxonomic compositions and functional potentials is still inevitable when a low-dimensional functional space is projected to a high-dimensional taxonomic space (Louca et al., 2018).