Autotrophic carbon fixation pathways along the redox gradient in oxygen-depleted oceanic waters.

Anoxic marine zones (AMZs), also known as 'oxygen-deficient zones', contribute to the loss of fixed nitrogen from the ocean by anaerobic microbial processes. While these microbial processes associated with the nitrogen cycle have been extensively studied, those linked to the carbon cycle in AMZs have received much less attention, particularly the autotrophic carbon fixation - a crucial component of the carbon cycle. Using metagenomic and metatranscriptomic data from major AMZs, we report an explicit partitioning of the marker genes associated with different autotrophic carbon fixation pathways along the redox gradient (from oxic to anoxic conditions) present in the water column of AMZs. Sequences related to the Calvin-Benson-Bassham cycle were found along the entire gradient, while those related to the reductive Acetyl-CoA pathway were restricted to suboxic and anoxic waters. Sequences putatively associated with the 3-hydroxypropionate/4-hydroxybutyrate cycle dominated in the upper and lower oxyclines. Genes related to the reductive tricarboxylic acid cycle were represented from dysoxic to anoxic waters. The taxonomic affiliation of the sequences is consistent with the presence of microorganisms involved in crucial steps of biogeochemical cycles in AMZs, such as the gamma-proteobacteria sulfur oxidisers, the anammox bacteria Candidatus Scalindua and the thaumarcheota ammonia oxidisers of the Marine Group I.

Introduction

Anoxic marine zones (AMZs) are open oceanic ecosystems in which dissolved oxygen concentrations are undetectable (≤1–10 nM), and nitrite accumulates to levels of ≥0.5 μM at intermediate depths of the water column (from ~100 to ~800 m depth; Thamdrup et al., 2012; Ulloa et al., 2012). Major permanent AMZs can be found in the Arabian Sea, the Eastern Tropical North Pacific (ETNP) and Eastern Tropical South Pacific Oceans (ETSP). Anaerobic microbial processes occurring in these systems impact the oceanic bioavailability of nitrogen at the planetary scale (Codispoti et al., 2001). While microbial processes associated with the nitrogen cycle have been extensively studied (Lam and Kuypers, 2011), those linked to the carbon cycle in AMZs have received much less attention, particularly the autotrophic carbon fixation, which is a crucial component of the carbon cycle. To date, the following six pathways of autotrophic carbon fixation have been identified (Berg et al., 2010; Hügler and Sievert, 2011): the Calvin‐Benson‐Bassham cycle (CBB), the reductive Acetyl‐CoA or Wood‐Ljungdahl pathway (rAcCoA), the reductive tricarboxylic acid cycle (rTCA), the 3‐hydroxypropionate bicycle (3HP), the 3‐hydroxypropionate/4‐hydroxybutyrate cycle (3HP/4HB), and the dicarboxylate/4‐hydroxybutyrate cycle (DC/4HB). Among them, the CBB cycle of photosynthesizers contributes most to the global synthesis of organic matter (Hügler and Sievert, 2011). Nevertheless, many microorganisms, including the anammox (anaerobic ammonium oxidation) bacteria (Woebken et al., 2008) and the abundant sulfur oxidisers (Walsh et al., 2009; Schunck et al., 2013; Léniz et al., 2017) in AMZs, have the genetic potential for carrying out autotrophic carbon fixation in dark marine ecosystems (Swan et al., 2011). Previous studies that have measured carbon fixation rates in the ETSP AMZ have reported a peak of values ranging from 9 to 130 nmol C L−1 d−1 below the photic zone (100–500 m), which has been attributed to the dark autotrophic carbon fixation (Löscher et al., 2016; Callbeck et al., 2018). However, no global study of these six pathways of autotrophic carbon fixation has been made along the redox gradient, from oxic to anoxic conditions, present in the water column of AMZs. Since some of the pathways utilize oxygen‐sensitive enzymes (Berg et al., 2010; Hügler and Sievert, 2011), a partitioning among them along this gradient is expected. To test this hypothesis, we analysed 42 metagenomes and 29 metatranscriptomes collected from geographically diverse AMZs (Fig. 1 and Supporting Information Table S1). The results show a clear vertical distribution (Figs. 2 and 3 and Fig. S1) of the relative abundances of sequences encoding the marker genes related to the six known pathways of autotrophic carbon fixation (Supporting Information Table S2; Kanehisa, 2000; Caspi et al., 2017).

Figure 1

Map of the origin of the metagenomes and metatranscriptomes from AMZs used in this study. The colour palette represents the oxygen concentrations at 300 m depth and was constructed with the data extracted from the CSIRO Atlas of Regional Seas (CARS2009: http://www.marine.csiro.au/~dunn/cars2009/). The yellow circles represent the location of the metagenomes and metatranscriptomes; details are given in the Supporting Information Table S1.

Figure 2

Relative abundances of a subset of the marker genes of autotrophic carbon fixation pathways, in the metagenomes extracted from AMZs worldwide, and their putative taxonomic affiliations. Metagenomes from the Arabian Sea, the ETNP, and ETSP Oceans were separated by their oxygen concentrations (boxes at the left) according to Thamdrup et al., 2012 and Wright et al., 2012. The sizes of the coloured circles are proportional to the relative abundance of the sequences per metagenome (computed as proportions of the number of reads to the total number of reads per data set and normalized according to the corresponding gene size). The taxonomy used in this figure corresponds to the third taxonomic rank in the NCBI taxonomy, which is normally associated with the class. Supporting Information Table S3 shows the taxonomy at the species level and Table S4 shows the percentages of coverage of the taxa for each marker gene present in the figure. The sequences from each of the datasets were directly aligned to a subset of the KEGG protein sequence database that only contained sequences with a defined KEGG orthology (KO). The algorithm BLASTX from the Diamond software package was used in all these massive sequence alignments. Only the alignments that had a bit‐score equal or higher than 50 were used in the downstream analyses. Afterward, the profiles of KOs and Enzyme Commission (EC) numbers, for each metagenome and metatranscriptome, were created. The mapping file ‘ko2ec’, from the KEGG distribution, was used to convert the KOs to ECs, whenever the KO corresponded to an enzyme. Subsequently, profiles for each metagenomic and metatranscriptomic dataset were created, with the relative abundances of the genes encoding the enzymes listed in the Supporting Information Table S2. The metagenomic and metatranscriptomic sequences associated with the marker genes of Supporting Information Table S2 were extracted from the original nucleotide FASTA files and aligned with BLASTX algorithm from the DIAMOND software package against the NCBI RefSeq, NCBI‐nr and a custom database composed of sequences associated with samples from low‐oxygen marine zones, available in IMG/ER. The sequences in the reference databases with the best bit‐score (whenever it was equal or higher than 50) for each input sequence (aligned to the marker genes only) was used as the putative taxonomic affiliation. When this taxonomic annotation was associated with an ‘uncultured organism’, the second or the third best bit‐scored alignment was used in this process. Singletons were eliminated, and only the fourth most abundant taxa related to each EC were considered in this figure. The results indicated that the distances determined by the correlations of the profiles of proportions of the marker genes (normalized as mentioned before) confirmed the predefined grouping by oxygen concentration with F = 18.10842 (P‐value = 9.999e‐05, Supporting Information Table S5).

Figure 3

Relative abundances and putative taxonomic affiliation of a subset of the marker genes of autotrophic carbon fixation in metatranscriptomes from AMZs. Metatranscriptomes were sampled from the ETNP and ETSP Oceans and were arranged by their oxygen concentrations (boxes at the left). The sizes of the coloured circles are proportional to the relative abundances of sequences, per metatranscriptome (computed as proportions of the number of reads to the total number of reads per meta‐ome and normalized according to the corresponding gene size). The taxonomy used in this figure corresponds to the third taxonomic rank in the NCBI taxonomy, which is normally associated with the class. Supporting Information Fig. S1 and Table S3 show the taxonomy at the species level and Table S4 shows the percentages of coverage of the taxa for each marker gene present in the figure. The results indicated that the distances determined by the correlations of the profiles of proportions of the marker genes (normalized as mentioned before) confirmed the predefined grouping by oxygen concentration with F = 2.78573 (P‐value = 0.02449755, Supporting Information Table S5).

Results and discussion

The analysis of metagenomes (Figs. 2 and S1) showed that marker genes related to the CBB cycle, like RuBisCO (EC4.1.1.39), were ubiquitous along the entire water column, but specially dominated the photic oxic (>90 μM O2), suboxic (<10 μM O2) and, in the case of the metagenomes from the ETSP, also the anoxic waters. In contrast, the sequences associated with the rAcCoA pathway were restricted to the anoxic and suboxic waters, being markedly more abundant in the anoxic metagenomes from the ETNP. Although significantly lower in relative abundances, the sequences related to the rTCA cycle were present in the range from dysoxic to anoxic waters, whereas those associated with the 3HP bicycle were present along the entire redox gradient. The sequences associated with the 3HP/4HB cycle were present along the entire redox gradient, showing more abundance in dysoxic waters (10–90 μM O2) and dominating in the lower oxycline. No sequences unambiguously assignable to the DC/4HB cycle could be detected in this analysis (Supporting Information Fig. S2). Results of the full set of enzymes for each pathway are presented in the Supporting Information Fig. S3.

The most abundant representative taxa for the studied genes in the metagenomes are shown in Figs. 2 and S1. The sequences associated with the CBB cycle were most prominently affiliated to Cyanobacteria (mainly to Prochlorococcus) and photosynthetic protists (Chlorophyta) in the illuminated (photic) waters, as well as to the gamma‐proteobacteria sulfur oxidisers (GSO) along the entire oxygen gradient underneath the photic waters. The ubiquitous presence of GSO in non‐sulfidic oxygen‐deficient marine environments have been associated with a cryptic sulfur cycle (Canfield et al., 2010 and Supporting Information Fig. S4), wherein GSO would be coupling the nitrogen and sulfur cycles, oxidizing reduced sulfur compounds with nitrate impeding the accumulation of sulfuide in these systems (Walsh et al., 2009; Canfield et al., 2010). Other authors have proposed an offshore transport of sulfuide‐rich waters to the open ocean by mesoscale eddies (Callbeck et al., 2018). The sequences associated with the rAcCoA pathway were mainly affiliated to Planctomycetes of the genus Candidatus Scalindua both in the ETSP and the ETNP, and to the genera Candidatus Brocadia and Jettenia in the Arabian Sea (see Supporting Information Fig. S1 and Table S3; Zarzycki et al., 2009). These proposed taxa have been reported to carry out the process of anammox in AMZs (Woebken et al., 2008 and Supporting Information Fig. S4), as well as being able to fix carbon autotrophically (Strous et al., 1999). Except for Chloroflexi, the sequences related to the 3HP bicycle were affiliated to phyla that have not been described to make use of this pathway so far, such as Actinobacteria, Gammaproteobacteria, and Alphaproteobacteria. This unexpected result suggests that this cycle could be more diversified in the tree of life than currently thought, or that an unknown pathway using shared enzymes with the 3HP bicycle could be operating in the ocean. The sequences associated to the 3HP/4HB cycle were found affiliated to Thaumarchaeota of the Marine Group I. Representatives of this phylum have been found abundant in AMZs oxyclines (Belmar et al., 2011), normally associated with the aerobic ammonium oxidation (Supporting Information Fig. S4), and making use of the most energy‐efficient aerobic pathway for CO2 fixation described so far (Könneke et al., 2014). Sequences related to the rTCA cycle (present from dysoxic to anoxic waters) were linked to Nitrospina gracilis, a known aerobic nitrite oxidiser. Therefore, the presence of these sequences in anoxic waters was unexpected. Genomic analyses of these bacteria have shown a versatile metabolism, including adaptations to microaerophilic conditions (Lücker et al., 2013). A recently discovered cryptic oxygen cycle could also support the presence of such aerobic bacteria in anoxic waters (Garcia‐Robledo et al., 2017).

It is evident from Figs. 2 and 3 that the frequency of the studied marker genes in the metatranscriptomes is much lower than in the metagenomes, which may indicate that these genes were not highly transcribed at the particular sampling times, even though the microbial community had the genetic potential for carbon fixation. Some authors have already pointed out that transcribed genes at a particular time represent a very small fraction of the genetic potential of microbial communities (e.g. García‐Ortega and Martínez, 2015). Also, limitations in the technique (e.g. transcripts of rRNA dominate metatranscriptomes, resulting in higher percentages of unassignable sequences than in metagenomes) may also contribute to the differences observed. Thus, the metagenomic data can represent the ecological potential of this low abundance genes (<1%) better than the metatranscriptomic data, whose analysis was only used to confirm those coming from the metagenomic analysis. Furthermore, as we pointed out before, there is quantitative evidence of carbon fixation from rate measurements under the photic zone, that are consistent with our results.

Therefore, in agreement with the metagenomic analysis, the metatranscriptomic data (Figs. 3 and S1) support the ubiquity of the CBB cycle along the entire redox gradient, although with the highest abundance of gene transcripts in the photic layer, related to photosynthetic microorganisms, and GSO in the rest of the AMZ water column with more prominence in the ETSP AMZ. Additionally, transcripts related to the rAcCoA and rTCA pathways, linked to Candidatus Scalindua and Nitrospina gracilis, respectively, were found exclusively in the suboxic and anoxic waters, stressing the point that these carbon fixation pathways thrive when oxygen is low or absent (Berg et al., 2010; Hügler and Sievert, 2011). The patterns of the gene transcription also support the prominent presence of the 3HP/4HB cycle from Thaumarchaeota in oxyclines. The observed pattern for carbon fixation is, thus, consistent with that observed when the metabolic processes of nitrogen and sulfur previously described are considered (Stewart et al., 2011; Hawley et al., 2014 and Supporting Information Fig. S4).

This study provides information about pathway‐specific metagenomic and metatranscriptomic analyses, demonstrating an explicit redox partitioning of the carbon fixation pathways in AMZs (Figs. 2, 3, S1, and Table S5 for statistical analysis). These redox gradients favour growth of chemolithotrophs that uses distinct autotrophy pathways with different sensitivities to oxygen. Thereby favouring, in the absence of light, additional metabolic pathways to the CBB cycle, such as the 3HP/4HB cycle and the rAcCoA pathway, which have been highlighted as the most energy‐efficient pathways of autotrophic carbon fixation in oxic and anoxic environments, respectively (Berg et al., 2010; Könneke et al., 2014). However, a novel variant of the rTCA cycle has been recently discovered (Nunoura et al., 2018, Mall et al., 2018) and estimated to be even more energetically efficient that the rAcCoA pathway (Mall et al., 2018). This novel pathway could not be analysed with our bioinformatic approach, as it uses a conventional citrate synthase for citrate cleavage instead of the rTCA key enzymes, thus lacking, so far, specific marker genes. Therefore, it cannot be ruled out that this rTCA variant, or another still unknown pathway, can be fixing carbon in AMZs.

Picocyanobacteria inhabiting the photic waters of AMZs can be responsible for up to 47% of the organic matter supplied to the anoxic waters of the ETNP and ETSP (Garcia‐Robledo et al., 2017), suggesting a quantitatively important role of the photosynthesis‐based CBB cycle in these oceanic areas. In oxygen‐depleted waters of the ETNP, an autoctonous organic matter production has also been related to cyanophages, in addition to cyanobacterial primary production (Fuchsman et al., 2019). A peak of carbon fixation rates below the photic zone of AMZs has also been observed, which has been attributed to the dark autotrophic carbon fixation by chemolithoautotrophs (Schunck et al., 2013; Löscher et al., 2016; Callbeck et al., 2018). Moreover, in sulfuide plume events that ocurr sporadically in the ETSP AMZ, dark carbon fixation by GSO (CBB cycle) could represent up to ~30% of the surface photoautotrophic carbon fixation (Schunck et al., 2013). These plumes can be exported to the open ocean by mesoscale eddies, enhacing dark carbon fixation offshore (Callbeck et al. 2018). However, in waters not influenced by eddies, only 7% of the total carbon fixation rate measurements would be attributed to the GSO (Callbeck et al., 2018), suggesting that in addition to the CBB cycle, other chemolithoautotrophic pathways are operating in these ecosystems.

Our results suggest the presence of carbon fixation along the entire redox gradient (from oxic to anoxic waters), although through different pathways and except for the CBB cycle, the quantitative contribution of fixed carbon by them remains unknown. Consequently, more measurements of carbon fixation rates along the entire gradient are necessary to better understand the importance of the carbon cycling in AMZs.

Supplementary Table 1 Metadata associated with the metagenomes and metatranscriptomes used in this study. The metagenomic and metatranscriptomic datasets were downloaded from the Sequence Read Archive (SRA), Integrated Microbial Genomes system (IMG/ER), and MG‐RAST databases.

Click here for additional data file.

Supplementary Table 2 EC numbers of enzymes involved in each of the autotrophic carbon fixation pathways of (A) the marker enzymes and (B) all the enzymes described in the six complete pathways for autotrophic carbon fixation. Marker genes encoding enzymes involved in autotrophic carbon fixation were collected from the literature (Berg et al., 2011; Hüger and Sievert, 2011), and checked in the KEGG and Metacyc databases (Kanehisa, 2000; Caspi et al., 2017) to preferably select those genes exclusively implicated in those pathways.

Click here for additional data file.

Supplementary Table 3 Marker genes of autotrophic carbon fixation pathways in metagenomes extracted from AMZs, and their 10 most abundant putative taxonomic association (at the species rank) for all marker genes. Note that sequences putatively related to the marine gamma proteobacterium HTCC2080 are present in euphotic waters, both for the marker genes of the 3HP/4HB and the 3HP pathways. These associations suggest the existence of an unknown pathway in bacteria inhabiting euphotic waters with common reactions to the 3HP/4HB and 3HP pathways. The presence of these genes has previously been reported for the HTCC2080 strain (Zarzycki et al., 2009).

Click here for additional data file.

Supplementary Table 4 Percentages of coverage of the putative taxonomic affiliation for each marker gene of autotrophic carbon fixation considered in Figs. 2, 3, S1 and S2.

Click here for additional data file.

Supplementary Table 5 Permutational analysis of variance (PERMANOVA) for the relative abundance of sequences aligned to all the marker genes listed in Table S2 for metagenomes and metatranscriptomes. To test whether the predetermined groupings of metagenomes and metatranscriptomes (separately) by oxygen concentration, were determined by the different similarity (correlation) matrices created on the basis of the relative abundance of reads aligned to all the marker genes listed in Table S2, we performed a permutational multivariate analysis of variance (PERMANOVA adonis2 implementation from the ‘vegan’ R package, number of permutations = 10000).

Click here for additional data file.

Supplementary Fig. 1 Relative abundance of the complete set of marker genes of autotrophic carbon fixation pathways in metagenomes (A) and metatranscriptomes (B) associated with AMZs, and their putative taxonomic associations. Meta‐omes are from the Arabian Sea, and the Eastern Tropical North Pacific (ETNP) and Eastern Tropical South Pacific (ETSP) Oceans and were arranged according to their oxygen concentrations. The sizes of the coloured circles are proportional to the relative abundance of the sequences per metagenome (computed as proportions of number of reads to the total number of reads per meta‐ome and normalized according to the corresponding gene size). The taxonomy used in this figure corresponds to the species level. Table S4 shows the percentages of coverage of the taxa for each marker gene present in the figure. In the cases of the CBB and acetyl‐CoA pathways, the marker genes show consistency among them, in the partitioning as well as in their taxonomic affiliation.

Click here for additional data file.

Supplementary Fig. 2 Relative abundances of genes encoding enzymes that participate in the thaumarchaeal 3HP/4HB cycle, predicted by Könnene and co‐authors (2014) in Thaumarchaeal genomes. In metagenomes (A) and metatranscriptomes (B) from AMZs (ECs: 6.2.1.‐, 6.4.1.2, 4.2.1.116, 6.4.1.3, 5.1.99.1, 5.4.99.2, 2.3.1.9, and 4.2.1.120). A special procedure was adopted for identifying the sequences associated with the 3HP/4HB and the DC/4HB cycles because they lack exclusive enzymes, and it was not possible to unambiguously assign sequences to them solely based on orthology. For the 3HP/4HB cycle, we chose five enzymes that have been reported as key enzymes by some authors (Berg et al., 2010; Hügler and Sievert, 2011 and Table S2), and that have low participation in other pathways. These enzymes were: 4‐hydroxybutanoyl‐CoA dehydratase (EC:4.2.1.120), 4‐hydroxybutyrate‐CoA ligase (EC:6.2.1.40), Acrylyl‐CoA reductase (EC:1.3.1.84), Malonyl‐CoA reductase (EC:1.2.1.75), and 3‐hydroxypropionyl‐CoA synthase (EC:6.2.1.36). From this set of enzymes, only 4‐hydroxybutanoyl‐CoA dehydratase had representative sequences associated with Thaumarchaeota in the KEGG database. This enzyme catalyses a reaction shared with the dicarboxylate/4‐hydroxybutyrate cycle (DC/4HB), the succinate fermentation to butanoate and the 4‐aminobutanoate degradation V. Thus, we made a complete analysis of all the sequences associated with the enzymes of this cycle. From this analysis, the enzymes that have already been predicted by Könnene and co‐authors in Thaumarchaeal genomes (ECs: 6.2.1.‐, 6.4.1.2, 4.2.1.116, 6.4.1.3, 5.1.99.1, 5.4.99.2, 2.3.1.9 and 4.2.1.120, Könneke et al., 2014) were found in our data and thus affiliated to Thaumarchaeota. Because of this redundancy and for the reasons explained above, even though there are no exclusive marker genes for this cycle, we used the 4‐hydroxybutanoyl‐CoA dehydratase as the marker gene for the 3HP/4HB cycle in Thaumarchaeota. The DC/4HB and 3HP/4HB cycles share the marker enzymes 4‐hydroxybutanoyl‐CoA dehydratase (EC4.2.1.120) and 4‐hydroxybutyrate‐CoA synthase (EC6.2.1.40) (Berg et al., 2010; Hügler and Sievert, 2011) and given that the 3HP/4HB cycle is the only one that has been reported in genomes from Thaumarchaeota, and the DC/4HB cycle has only been reported in the crenarchaeal genomes such as Thermoproteales and Desulfurococcales (Berg et al., 2010), we relied on taxonomy to discriminate between them and discarded the second.

Click here for additional data file.

Supplementary Fig. 3 Complete set of enzymes for each autotrophic carbon fixation pathway for the metagenomes (A) and metatranscriptomes (B). Pathways are represented by pie charts comprising all of their enzymes as slices (Table S2, B), which are either coloured or not coloured according to their presence or absence, respectively. This figure is to analyse the completeness of each of the six metabolic pathways of autotrophic carbon fixation. Marker enzymes are in bold typesetting. These are the key enzymes or marker genes that we use to evaluate the presence of each of the metabolic pathways since they perform exclusive reactions of these pathways, except for the 3HP/4HB and DC/4HB cycles in which we also use the putative taxonomic affiliation of the sequences, since they share the marker genes (see legend of Fig. S2 for more details).

Click here for additional data file.

Supplementary Fig. 4 Diagram of the coupling of the biogeochemical cycles of nitrogen, sulfur and the autotrophic carbon fixation pathways in AMZs.

Click here for additional data file.