Phylogenetically and catabolically diverse diazotrophs reside in deep-sea cold seep sediments

Abstract

Microbially mediated nitrogen cycling in carbon-dominated cold seep environments remains poorly understood. So far anaerobic methanotrophic archaea (ANME-2) and their sulfate-reducing bacterial partners (SEEP-SRB1 clade) have been identified as diazotrophs in deep sea cold seep sediments. However, it is unclear whether other microbial groups can perform nitrogen fixation in such ecosystems. To fill this gap, we analyzed 61 metagenomes, 1428 metagenome-assembled genomes, and six metatranscriptomes derived from 11 globally distributed cold seeps. These sediments contain phylogenetically diverse nitrogenase genes corresponding to an expanded diversity of diazotrophic lineages. Diverse catabolic pathways were predicted to provide ATP for nitrogen fixation, suggesting diazotrophy in cold seeps is not necessarily associated with sulfate-dependent anaerobic oxidation of methane. Nitrogen fixation genes among various diazotrophic groups in cold seeps were inferred to be genetically mobile and subject to purifying selection. Our findings extend the capacity for diazotrophy to five candidate phyla (Altarchaeia, Omnitrophota, FCPU426, Caldatribacteriota and UBA6262), and suggest that cold seep diazotrophs might contribute substantially to the global nitrogen balance.

Fig. 1. Geographic distribution of 11 cold seep sites where metagenomic and metatranscriptomic data were collected.

These samples were originated from five types of cold seeps: gas hydrates, mud volcanoes, asphalt volcanoes, oil and gas seeps and methane seeps. Sites with red asterisks denote that both metagenomes and metatranscriptomes were collected, sites with blue asterisks denote that only metatranscriptomes were collected, and sites without asterisks denote that only metagenomes were collected. Also see details in Supplementary Data 1. The world map was drawn using the ggplot2 package in R v4.0.3.

Fig. 2. Maximum-likelihood phylogenetic tree of non-redundant nitrogenase subunit NifH identified from cold seep metagenomic assemblies.

Homologues of nifH were classified into canonical groups I to III, nitrogen fixation-like groups IV to VI, and newly assigned groups including groups of Methanosarcina-like, VII and VIII. NflH denotes NifH-like sequences. Scale bar indicates the mean number of substitutions per site.

Fig. 3. Relative abundance patterns of 202 nifH genes.

a Relative abundances of 202 nifH genes from different cold seep sediments, shown as RPKM (reads per kilobase per million mapped reads). b Comparison of nifH gene abundances in different types of cold seep ecosystems. n values refer to the number of biologically independent samples for statistics analysis. Asterisks indicate statistically significant differences between groups of mud volcanoes, oil and gas seeps, and methane seeps (determined by two-sided Wilcoxon Rank Sum test; ^* for P < 0.05, ^** for P < 0.01 and ^*** for P < 0.001). Boxplot components: center line, median values; box limits, upper and lower quartiles; whiskers, 1.5× interquartile range; points, outliers. c Significant Spearman correlation between relative abundances of nifH and the oxidative mcrA gene. Percentages were calculated by dividing the RPKM value of nifH genes by the mean of RPKM values estimated from 14 single-copy marker genes. The gray shadow indicates the 95% confidence interval. The abbreviations of the sites are shown in Fig. 1. Source data are provided as a Source Data file.

Fig. 4. Maximum-likelihood phylogenetic trees of nitrogen-fixing MAGs and their NifH protein sequences.

a Phylogenomic analysis of 35 MAGs containing nitrogen fixation genes. This maximum-likelihood tree is based on concatenation of 43 conserved protein sequences. MAGs are colored based on their phylogenetic affiliation at the phylum level. b Phylogenetic analysis of identified NifH protein sequences and genomic context of corresponding nif genes in the same 35 MAGs with nitrogen fixation genes. The scale bar represents one amino acid substitution per sequence position. For both trees, bootstrap values of >70% are indicated as black circles at the nodes, and scale bars indicate the mean number of substitutions per site.

Fig. 5. Metabolic reconstruction of core pathways for nitrogen-fixing MAGs.

a Anaerobic archaeal oxidation of methane; b dissimilatory sulfate reduction; c archaeal methanogenesis; d anaerobic degradation of alkanes by bacteria. Red font indicates that not all MAGs retrieved include the gene (numbers of MAGs with the corresponding gene indicated in parentheses). The percentages between brackets indicate the estimated completeness of the corresponding MAGs. Mtr N⁵-methyltetrahydromethanopterin–coenzyme M–methyltransferase, Mer 5,10-methylenetetrahydromethanopterin reductase, Mtd methylenetetrahydromethanopterin dehydrogenase, Mch methenyltetrahydromethanopterin cyclohydrolase, Ftr formylmethanofuran-tetrahydromethanopterin N-formyltransferase, Fwd formylmethanofuran dehydrogenase, Hdr heterodisulfide reductase, APS adenosine phosphosulfate, Apr adenylylsulfate reductase; Sat sulfate adenylyltransferase, GS glutamine synthetase, GOGAT glutamate synthase, GDH NADH-glutamate dehydrogenase, Ass alkylsuccinate synthase, AssK CoA-ligase, Mcm methylmalonyl-CoA mutase, Pcc propionyl-CoA carboxylase, Acd acetate-CoA ligase (ADP-forming). Detailed enzyme annotation is presented in Supplementary Data 7.

Fig. 6. Genomic context of nitrogen fixation genes.

Gene neighborhoods of nifHDK include retrotransposable or transposon elements, regulatory nitrogen fixation genes, nitrogenase metal cofactor biosynthesis genes and transporter genes. Source data are provided as a Source Data file.

Fig. 7. Evolutionary metrics of nitrogen fixation genes.

a Nucleotide diversity (π) of nifHDK genes at different types of cold seeps; b pN/pS ratio of nifHDK genes at different types of cold seeps. Nucleotide diversity is used to measure genetic diversity within a population (microdiversity), which is calculated using the formula: 1 − [(frequency of A)² + (frequency of C)² + (frequency of G)² + (frequency of T)²]. pN/pS is the ratio of non-synonymous to synonymous polymorphism rates within a population. n values refer to the number of biologically independent samples for statistics analysis. The significances were analyzed by two-sided Kruskal–Wallis Rank Sum test. Boxplot components: center line, median values; box limits, upper and lower quartiles; whiskers, 1.5× interquartile range; points, outliers. Source data are provided as a Source Data file.