Redox traits characterize the organization of global microbial communities.

The structure of biological communities is conventionally described as profiles of taxonomic units, whose ecological functions are assumed to be known or, at least, predictable. In environmental microbiology, however, the functions of a majority of microorganisms are unknown and expected to be highly dynamic and collectively redundant, obscuring the link between taxonomic structure and ecosystem functioning. Although genetic trait-based approaches at the community level might overcome this problem, no obvious choice of gene categories can be identified as appropriate descriptive units in a general ecological context. We used 247 microbial metagenomes from 18 biomes to determine which set of genes better characterizes the differences among biomes on the global scale. We show that profiles of oxidoreductase genes support the highest biome differentiation compared with profiles of other categories of enzymes, general protein-coding genes, transporter genes, and taxonomic gene markers. Based on oxidoreductases' description of microbial communities, the role of energetics in differentiation and particular ecosystem function of different biomes become readily apparent. We also show that taxonomic diversity is decoupled from functional diversity, e.g., grasslands and rhizospheres were the most diverse biomes in oxidoreductases but not in taxonomy. Considering that microbes underpin biogeochemical processes and nutrient recycling through oxidoreductases, this functional diversity should be relevant for a better understanding of the stability and conservation of biomes. Consequently, this approach might help to quantify the impact of environmental stressors on microbial ecosystems in the context of the global-scale biome crisis that our planet currently faces.

Biological communities are conventionally described as assemblages of species whose ecological roles are known or predictable from their observable morphological characteristics. In the early twentieth century, Lotka and Volterra pioneered the development of theoretical ecology using species numbers as the master variable in differential equations that describe the interactions and complexity of ecological systems (1). Since then, most theoretical ecologists have used species numbers as the ecological unit for developing an extensive body of theory, which includes elaborate mathematical models to explain the dynamics of populations and communities (1). In practice, this approach requires the categorization of every observed individual into a taxonomic unit —which is not a trivial task in some cases (2), and it is definitively a problem in microbial ecology (3–5). In the latter context, microbial ecologists face three main problems. First, observable morphological attributes do not provide sufficient discriminatory or functional characterization. Second, the isolation of microbial species to assess their physiology and ecological function is rarely possible, a phenomenon that is related to the so-called Great Plate Count Anomaly (6). Third, prokaryotic genomes are highly dynamic, mainly due to pervasive horizontal gene transfers and the effect of mobile DNA elements and phages (7). Microbial ecologists have employed molecular taxonomic markers, primarily the small subunit ribosomal RNA gene, to address the first and second problems, thereby operationally defining species and estimating their abundances and taxonomic diversity (8). This taxonomic approach has been used to explain and predict the microbial dynamics in diverse environments (9, 10). In such a context, the Earth Microbiome Project initiative has recently reported microbial taxonomic diversity per biome on a global scale with the use of standardized protocols to provide an organized and complete catalog of microbes (11). However, several studies have reported inconsistent taxonomical correlations under apparently similar ecological scenarios, finding better consistency only when using multiple protein-coding genes as traits and when the whole community is analyzed as the ecological unit (12–16). This has been performed in an attempt to address the third above-mentioned problem. After all, it is the function, not the taxonomic information, that has the actual ecological relevance (17). Unfortunately, the selection of the ecologically relevant categories of protein-coding genes for use is not evident in the broad context of planetary biomes (6, 18, 19). We analyzed 247 metagenomes from 18 biomes (Fig. 1) to tackle this issue and to determine under which specific nonexclusive set of genes the differences between biomes are the highest. These gene sets included protein-coding genes with associated orthology in the KEGG database (a typical approach in trait-based analyses), enzyme-coding genes, transporter-associated genes, and taxonomic marker genes (Fig. 1). We found that the set of genes that were encoding enzymes better differentiated the biomes than the other gene categories. In particular, the profiles of genes that were encoding oxidoreductases composed the set with the highest cohesion and separation of biome groups, suggesting that they can better describe the association of the microbial communities to their respective biomes. In addition, we found no correspondence in biome maximum diversity between the functional and the taxonomic approaches. An oxidoreductase-based description of microbial communities also serves as a convenient proxy for an energetic description of ecosystems as these proteins are responsible for redox reactions, which are the processes by which every living organism uses energy from and modify the chemical characteristics of the environment (20).

Fig. 1.

Biomes and categories of genes. (A) Sketch of the biomes from which metagenomes (as proxies for microbial communities) were included in this paper. The animal-associated biome includes metagenomes from terrestrial animals only. A complete list and origin of these metagenomes can be found in the . (B) Organized list of the biomes illustrated in A. The number of metagenomes per biome is shown in parentheses besides the biome name, which is displayed in the color code utilized in the rest of the figures. (C) Categories of gene profiles considered in the analyses. All protein orthologies refer to the protein orthologies present in the KEGG protein database. The fourth rank taxonomy typically corresponds to a phylum in the prokaryotic taxonomy (see for details).

Results and Discussion

From an ecological point of view, the functions of communities represent the most relevant information about an ecosystem. In microbial ecology, when these functions are fine grained to molecular processes through functional genes, it is natural to ask whether all of them have the same ecological relevance to differentiate one biome from another (Fig. 1). Our results show that redox functions supported the highest statistical differentiation among biomes when taxonomic and functional sets of genes were compared (Table 1 and ). The discriminatory power of oxidoreductase genes for grouping biomes can be visualized in networks of correlations using different gene categories with metagenomes as nodes (microbial communities, colored according to biome origin) and correlations as edges (Fig. 2 and ). Metagenomes from different biomes were more separated in the networks of oxidoreductases than in the network of taxonomic markers, which is the visual expression of the better cohesion and separation results as shown in Table 1. Hierarchical clusterings of these profiles (Fig. 2 and ) revealed the following three main groups of biomes: a group of apparent anoxic or suboxic biomes (animal associated, some hot springs, subterranean ecosystems, marine sediments, subseafloor, and mangrove sediments), a group of aquatic biomes (freshwater and different types of marine ecosystems), and a group of soil-associated biomes (grassland, forest, deserts, and rhizosphere). Note that environments associated with oxygen minimum zones did not cluster with the first above-mentioned group. The oxygen-limited condition shared by these ecosystems is not reflected in this clustering because the microorganisms in the pelagic low-oxygen environments mainly exploit chemolitoautotrophic metabolisms instead of the anaerobic degradation of organic matter that normally occurs in, for example, anoxic sediments. This analysis also showed that metagenomes from extreme ecosystems, such as acidic cave biofilms, some hot spring systems, and hypersaline environments, cluster outside of these three main groups.

Table 1.

Different profiles of genes separating metagenomes into biomes groups

Category	PERMANOVA F statistic	Cohesion	Separation	Cohesion + separation	Profile size
All KEGG protein orthologies	10.66858	0.59033	−0.33468	0.25565	6,789
All enzymes	10.64720	0.69027	−0.43105	0.25922	1,826
Oxidoreductases	12.47645	0.69076	−0.42962	0.26113	484
Transferases	9.64637	0.76038	−0.52017	0.24021	541
Hydrolases	12.16139	0.70443	−0.45038	0.25405	423
Lyases	11.21940	0.75581	−0.51344	0.24237	211
Isomerases	7.80187	0.80306	−0.57657	0.22649	103
Ligases	9.97776	0.86361	−0.69899	0.16463	94
Transporters	7.23900	0.44584	−0.25863	0.18721	1,869
Taxonomy (species)	1.83457	0.13495	−0.06620	0.06875	4,011
Taxonomy (fourth rank)	2.81049	0.36668	−0.24392	0.12275	365

Different sets of profiles of relative gene abundances (with stabilized variances) were evaluated to determine under which of them the separation of metagenomes (microbial communities) into biome groups was most significant. The PERMANOVA statistical test (all P values < 0.001) indicates that gene profiles of oxidoreductases were the set of genes with higher statistically supported differences between the biomes (the higher the F statistic, the more likely to reject the null hypothesis of no differences among groups). The cohesion/separation of biomes groups confirmed that result (third, fourth, and fifth columns). The last column of this table indicates the total number of subcategories that each category displayed in this paper. For example, considering all of the datasets analyzed, we found 484 oxidoreductase genes and 4,011 species. To resolve the low differences among the F statistics of some enzyme categories, additional analyses were required, and their results can be found in .

Fig. 2.

Association of microbial metagenomes and biomes. (A) Network representation of the microbial metagenomes by profiles of oxidoreductase and taxonomic gene ranks. The nodes correspond to metagenomes, colored according to their biome of origin (Fig. 1). The edges represent maximal information coefficients (MICs). In the network associated with oxidoreductases (Left), all MIC ≥ 0.5 are shown. The network associated with taxonomic profiles (Right) was drawn with all MIC ≥ 0.1 as in this case, and these values were significantly lower. These taxonomic edge weights were increased by 0.4 to give visual balance to the plot. The differential clustering of biomes (colors) is explained by the values of cohesion and separation from Table 1 (see also ). (B) A simplified version (topology only and grouped per biomes) of the hierarchical clustering of the metagenomes based on oxidoreductase gene profiles (). Support values higher than 90% are shown in the plot.

The group of biomes with apparent anoxic conditions shared distinctive oxidoreductase genes related to methanogenesis, sulfide oxidation, denitrification, hydrogen oxidation, nitrogen fixation, and aromatic aldehydes oxidation (Fig. 3). The animal-associated metagenomes analyzed here were highly diverse, but most of them were related to the digestive systems of animals, making this group slightly biased toward the functional genes that are represented more in these microbial communities. Thus, the functions associated with these diverse biomes should be interpreted with care as it is unlikely that, for example, the human tongue dorsum supports microbial communities exploiting hydrogen oxidation processes. Indeed, hierarchical clusterings separated the microbial communities associated with the parts at the end of the digestive system of animals (cecum, gut, and stool) from other animal-associated metagenomes (human oral mucosa, tongue dorsum, supragingival plaque, anterior nares, and posterior fornix; ). Although the latter subgroup of microbial communities can also be associated with potentially anoxic microhabitats, the former subgroup was found to be functionally closer to the communities from the marine sediments and subsea-floor ecosystems, mainly because of the shared redox functionalities for the degradation of organic matter under anoxic conditions. Notably, gut-associated microbiomes displayed nitrogen fixation capabilities too (Fig. 3), which is consistent with the recent observations (21).

Fig. 3.

Distinctive oxidoreductase genes associated with biomes. These genes were determined by statistically testing that the average of rankings of each oxidoreductase gene within each biome was significantly different from the average ranking in other biomes. Dark and light shades in this figure refer to relative abundances, high and low, respectively. Thus, the rankings for this figure were reversed as a low rank indicates high relative abundance. These values were scaled for better visualization, which means that color shades can only be compared horizontally. Some of these distinctive genes encode oxidoreductases considered associated with important biogeochemical and biochemical processes. For example, CoB-CoM heterodisulfide reductase (methanogenesis), sulfide:quinone reductase (sulfide oxidation), nitrite reductase NADH (denitrification), hydrogenase (hydrogen oxidation), nitrogenase (nitrogen fixation), and aldehyde ferredoxin oxidoreductase (aromatic aldehydes oxidation). Hierarchical clusterings using these values were calculated for convenient grouping of both biomes and oxidoreductase genes.

Marine microbial communities were best characterized by a group of oxidoreductases that includes dimethylglycine dehydrogenase, sarcosine oxidase, and choline dehydrogenase (Fig. 3). These enzymes are involved in the synthesis and degradation of glycine betaine, which is an effective and widely used compatible solute for coping with saline stress (22). Indeed, most algae and some invertebrates produce and accumulate glycine betaine as an intracellular osmolyte (22). Thus, marine microorganisms might take advantage of the availability of this substrate in seawater and can convert it to formate, which can then be used as an energy source or directed to one-carbon metabolism for biosynthesis (23). A direct precursor of glycine betaine is choline, which is also abundant in seawater, as it can represent up to 0.39% of the dry weight of algae (24). A distinctive oxidoreductase gene present in marine microbial communities was 3-hydroxyisobutyrate dehydrogenase, which has been found to play a role in amino acid catabolism (25), as a source of alternative substrates for respiration under metabolic stress situations. Another representative of oxidoreductase encoded in the metagenomes of these microbial communities is aldehyde dehydrogenase NAD+. Polyunsaturated aldehydes are commonly produced by diatoms as a chemical defense against grazers, and their concentrations in seawater can potentially affect the bacterial community structure and diversity (26).

Microbial communities associated with soil were mainly characterized by oxidoreductase genes related to the degradation of aromatic compounds for the carbon source [alcohol dehydrogenase cytochrome c, isoquinoline 1-oxidoreductase, catechol 2,3-dioxygenase, homogentisate 1,2-dioxygenase (27) and phenylacetyl-CoA 1,2 epoxidase (28) (Fig. 3)]. This representation might be explained by the fact that most primary production in soils is returned to the environment as detritus (29), which can be rich in aromatics as they constitute a significant part of lignin in higher plants (27). Genes encoding betaine aldehyde dehydrogenase were also distinctive in soil-associated microbial communities. This enzyme is involved in the biosynthesis of glycine betaine as a compatible solute for alkaline-saline stress (30). In fact, reports indicate that many soil environments are highly alkaline, and transient conditions, such as drought, can significantly increase the alkalinity within cells (31). Additionally, plant root exudates can change the soil chemistry, sometimes creating microhabitats of increased alkalinity (30). Thus, soil microbial communities seem to be genetically prepared to resist saline-alkaline stress by synthesizing their cellular defenses, unlike marine microbial communities that apparently rely more on the environmental availability of glycine betaine or its direct precursors, such as choline or sarcosine. Despite freshwater biome grouping with the marine biomes, its associated microbial communities still share similarities in the abundances of some oxidoreductase genes with the soil biomes, such as in the case of betaine aldehyde dehydrogenase, carbon monoxide dehydrogenase (acceptor), and stearoyl-CoA 9-desaturase (Fig. 3). This observation might be related to the results of a recent study that suggest that freshwater ecosystems might connect the otherwise separated microbial communities (32).

Although most biogeochemical processes are widely distributed across different environments (33), some oxidoreductase genes associated with these processes appear to be unimportant for soil and aquatic biomes. This apparent conflict can be explained by the fact that, frequently, the most abundant microbes in these environments are heterotrophs [e.g., members of Acidobacteria in soils (34) and SAR11 clade in the ocean (35)]. Thus, although nitrification, denitrification, sulfur oxidation, and carbon fixation also occur in terrestrial and aquatic ecosystems, their genetic markers are significantly less abundant than the oxidoreductase genes related to heterotrophic metabolisms (). On the other hand, the biomes from the apparently anoxic group (typically harboring fewer heterotrophs) appeared prominently in many of these processes, such as, for example, methanogenesis, hydrogen oxidation, nitrogen fixation, sulfur oxidation, nitrification, and denitrification (Fig. 4). In addition, the oxidative phosphorylation process under suboxic conditions (associated with Cbb3 oxidase, encoded by the ccoN gene, Fig. 4) appeared to be best ranked in these biomes. Despite the pelagic low-oxygen marine biome was not clustered in this group of biomes (Figs. 2 and 3 and ), their metagenomes displayed high genetic representation associated with some of these processes, such as nitrification, denitrification, and sulfur oxidation (Fig. 4). This fact has been described as the beginning of a progressive rerouting of the energy flow into the microbial pathways as oxygen declines in marine ecosystems in detriment of higher trophic levels (36–38). Such progression ends in the extreme situation in which all benthic energy is processed as hydrogen sulfide (36) with concomitant accumulation of nitrite in the intermediate case of the anoxic marine zones (39). Low-oxygen areas in the ocean have rapidly expanded in the past decades, and they are expected to further increase as a consequence of global warming (36, 38). This, in turn, can be affected by the greenhouse gases that are emitted in marine low-oxygen zones as a by-product of anaerobic microbial pathways (36, 38, 39).

Fig. 4.

Biogeochemically relevant processes per biome by oxidoreductase genes relative abundances. Oxidoreductase genes associated with biogeochemical processes and their top five biomes where they were ranked the highest. The biomes per genes are in clockwise order, starting from the biome where the gene was best ranked. For example, the dissimilatory sulfite reductase gene (dsrA; involved in sulfur oxidation and reduction) was found best ranked in the following biomes in this order: hydrothermal vents, subterranean habitats, mangrove sediments, hot springs, and oxygen minimum zones.

The extraordinary dispersal potential of microbes is usually expressed through the old tenet “everything is everywhere, but the environment selects,” which a recent study extends to “every gene is everywhere, but the environment selects” (32). This fact suggests that measures of diversity for conducting large-scale studies of biomes in microbial ecology should include not only the richness, but also the evenness of the distribution of gene categories. By using the inverse Simpson index, we found that microbial taxonomic diversity does not correlate with microbial functional diversity. In our analysis, microbial communities from mangrove sediments were found to be the most taxonomically diverse (Fig. 5). This result is consistent with findings of the recent studies that show that some sediment environments can be more diverse than soils (40), which, in turn, have been traditionally considered to be the ecosystems with the highest microbial diversity (41). However, regarding oxidoreductase genes, grassland soils and rhizospheres were found to be the most diverse biomes (Fig. 5). This finding correlates with observations in plant diversity that suggest that, in the fine grain, grasslands are the most diverse soil biomes, harboring up to ∼90 different plant species per square meter (42). It is noteworthy that the temperate grasslands are currently among the biomes that face the highest ecological risk due to the extensive habitat loss and underprotection (43). To give a quantitative example of the microbial diversity of oxidoreductase genes in grasslands, consider that, on average, ∼130 of their most abundant categories were needed to cover the 70% of the total abundance of these genes. The same coverage percentage needed only ∼40 of the most abundant categories in the subterranean and acidic cave biofilms biomes (Fig. 5).

Fig. 5.

Microbial diversity of biomes. (A) Heat map plot constructed with the inverse Simpson diversity index (true diversity with q = 2) of the taxonomic and functional profiles for the metagenomes, averaged per biome. The dark color shades indicate high diversity. These average values were scaled per profile category for homogenous contrast. Thus, the colors can only be compared along columns, i.e., by biome. For example, regarding oxidoreductase genes, the grassland biome is the most diverse, and the rhizosphere is the second one. On the other hand, the subterranean biome is shown as the less diverse biome in almost every gene category. Note that “All proteins” refer to all proteins with defined orthology in the KEGG database (see for details). (B) Average number, per biome, of oxidoreductase genes (vertical axis) necessary to cover different percentages of total oxidoreductase genes, counted from the most to less abundant. For example, the 60 most abundant oxidoreductase genes in grassland-associated datasets on average covered ca. 45% of the total pool of oxidoreductase genes.

The choice of relevant variables is a critical step in the analysis of any complex system. In microbial ecology, the taxonomic structure of communities has typically been considered a proxy for the microbial ecosystem’s functioning, even though it is often unable to resolve functional genetic traits (44). The need for alternative trait-based approaches has been claimed for years (45), but there has been no agreement on the selection of a relevant set of genes necessary for its practical application (6, 18, 19). In this paper, we evaluated different sets of genes for this purpose, finding that oxidoreductase genes are a convenient choice. The set of transporter genes also has this potential, but its power to differentiate biomes was found to be lower. This is most likely as these genes also suffer from significant redundancy (e.g., there are different transporters for the same substrate, depending on their affinities). Other groups of enzyme genes, such as those associated with hydrolases, also supported a proper separation of biomes (Table 1 and ); however, they are slightly related to biogeochemical processes, mainly through the carbon cycle. In contrast, oxidoreductases are directly involved in most biogeochemical processes and nutrient recycling in every environment. Thus, the diversity of these functions should be relevant to better understand the stability and conservation of biomes, affected by the high disparities between ecosystem conversion and conservation across biomes, which has been recognized as comprising an ongoing biome crisis (43). Indeed, conservation efforts have mainly focused on particular species or local macrocommunities (e.g., polar bears and coral reefs, respectively) but not on the microbial ecological functions that sustain trophic levels, biogeochemical cycles, and the ecosystem services that are derived from them. This omission is likely due to the difficulty of predicting microbial ecosystem dysfunction from environmental stressors using microbial taxonomy information (46). We expect that an oxidoreductase-based description of microbial communities should facilitate this task and help to quantify in future developments the impact of environmental changes on microbial ecosystem functions in the context of the global-scale biome crisis that our planet currently faces.

Materials and Methods

Data Collection and Sequence Analysis.

The metagenomic datasets were collected from metagenomic studies of diverse microbial communities in recent years. The selection of metagenomes was guided by literature search, trying to cover the biomes with at least three “whole genome amplified” metagenomes sequenced with 454 or Illumina technologies. This process resulted in 247 metagenomes, grouped in 18 biomes (Fig. 1). The sources of these datasets are listed in the . The sequences of these datasets were aligned against different protein sequence databases () using the BLASTX algorithm of the DIAMOND software with a bit-score cutoff of 50. With these alignment results, the different profiles listed in Fig. 1 and Table 1 were constructed.

Group Variances Analyses.

The PERMANOVA statistical test was used to assess and compare the degree of separation of metagenomes (microbial communities) into biome groups by using the data profiles (Table 1) with dissimilarity matrices constructed with distances calculated based on nonparametric correlations [MIC and Spearman].

Diversity Estimation.

Each profile of categories, for all of the metagenomic datasets (Fig. 1), was first resampled by a coverage percentage of 95%. True diversity was calculated on the resampled datasets by using the inverse Simpson index. The diversity per biome was calculated as the average of the diversities of all metagenomic datasets from each biome ().

Networks and Clustering.

For each pair of profiles described above, a distance between them was calculated as 1 correlation (correlation as the pairwise maximal information coefficient between the profiles). The networks of metagenomes (Fig. 2) were constructed by writing the graph in the graph exchange XML format and rendered using the Gephi software with the OpenOrd network layout. The hierarchical clustering of biomes was computed with the R package Pvclust with 104 permutations and with a distance based on the Spearman correlation. The genes in Fig. 3 were selected as the top three oxidoreductase genes from each biome whose average ranking was lower than the total average. More details about all these procedures can be found in the .