nifH gene expression and diversity in geothermal springs of Tengchong, China.

Terrestrial hot springs have been suggested to harbor diverse diazotrophic lineages by using DNA-based nifH gene phylogenetic analysis. However, only a small amount of diazotrophs were ever confirmed to perform nitrogen fixation. In order to explore the compositions of active diazotrophic populations in hot springs, the in situ expression and diversity of nifH and 16S rRNA genes were investigated in the sediments of hot springs (pH 4.3-9.1; temperature 34-84°C) in Tengchong, China, by using high-throughput sequencing. The results showed that active diazotrophs were diverse in the studied Tengchong hot springs. The main active diazotrophs in high-temperature hot springs were affiliated with Aquificae, while those in low-temperature hot springs belonged to Cyanobacteria and Nitrospirae. Such dominance of Aquificae and Nitrospirae of diazotrophs has not been reported in other ecosystems. This suggests that hot springs may harbor unique active diazotrophs in comparison with other type of ecosystems. Furthermore, there were significant differences in the phylogenetic lineages of diazotrophs between hot springs of Tengchong and other regions, indicating that diazotrophs have geographical distribution patterns. Statistical analysis suggests that the expression and distribution of nifH gene were influenced by temperature and concentrations of ammonia and sulfur seem in Tengchong hot springs. These findings avail us to understand element cycling mediated by diazotrophs in hot spring ecosystems.

Introduction

Biological dinitrogen (N2) fixation, as an important nitrogen source of nitrogen-limited ecosystems, provides bioavailable fixed nitrogen for the synthesis of biological macromolecules, metabolism, and cell growth of microbial community (Falkowski, 1997; Zehr et al., 2003). Nitrogen is usually limited in terrestrial geothermal springs due to the rapid assimilation and transformation of inorganic nitrogen (Estrella-Alcamán et al., 2015; Lin et al., 2015). Therefore, nitrogen fixation processes may exist to support microbial demand for nitrogen in hot springs.

With the use of DNA-based methods, nifH gene coding for the ferric protein subunit of the nitrogenase complex of diazotrophs has been proved to be widespread in hot springs (temperature: 35-84°C; pH: 2-9) of Yellowstone National Park (YNP, United States) (Miller et al., 2006; Hall et al., 2008; Hamilton et al., 2011a). Subsequently, nifH gene was also detected in hot springs in Porcelana (Chile; Estrella-Alcamán et al., 2015), Nakabusa (Japan; Nishihara et al., 2018a) and Iceland (Cousins et al., 2018). These previous studies showed that hot spring diazotrophs are diverse and mainly belong to Cyanobacteria, Proteobacteria, Firmicutes, Chloroflexi, Aquificae, and Nitrospirae. However, the DNA-based methods can only indicate the presence of diazotrophs, but cannot reveal their activity and whether they promote nitrogen fixation process in hot springs. To confirm the presence of active diazotrophs, in situ expression of nifH gene and(or) nitrogenase has been investigated in some YNP hot springs (Steunou et al., 2006, 2008; Hamilton et al., 2011b,2014; Loiacono et al., 2012) and Porcelana (Estrella-Alcamán et al., 2015; Alcamán-Arias et al., 2018). Generally, two cyanobacterial genera, i.e., Mastigocladus and Synechococcus, were confirmed as active diazotrophic taxa in hot springs with photosynthesis (Steunou et al., 2006, 2008). In contrast, in chemotrophic hot springs, most nifH gene transcripts were related to the Aquificae genera of Thermocrinis and Hydrogenobacter, while only a few nifH gene transcripts were affiliated with Nitrospirae (genus Thermodesulfovibrio) and Chloroflexi (genus Roseiflexus). These previous results indicated that active diazotrophs of hot springs were just limited to a few taxa. However, considering the high diversity of diazotrophs revealed by DNA-based studies, it is reasonable to hypothesize that more undetected groups of diazotrophs may play an important role in nitrogen fixation process in hot springs.

In addition, one interesting finding of previous in situ studies is that Aquificae-related nifH gene transcripts dominate in the YNP high-temperature hot springs, but have never been reported to appear in other types of environments (Loiacono et al., 2012; Hamilton et al., 2014). This suggests that the composition of diazotrophs is very unique in terrestrial hot springs. Indeed, Aquificae-related bacteria widely exists in terrestrial hot springs with high temperature chemosynthetic function and play important roles in the cycling of carbon, sulfur and hydrogen (Hedlund et al., 2015). These above findings suggest Aquificae may also be the main executors of nitrogen fixation in global hot springs with chemosynthesis. This speculation was supported by the retrieval of two Aquificae strains capable of nitrogen fixation isolated from one hot spring in Nakabusa, Japan (Nishihara et al., 2018b). But more evidence is needed to prove the expression of nifH gene by Aquificae.

Temperature is regarded as one important factor shaping the microbial community of hot springs (Meyer-Dombard et al., 2005; Hou et al., 2013). In normal environments, multiple factors are related to the distribution of diazotrophs. For example, pH and available carbon and nitrogen shape the diazotrophic composition in agricultural soils of China (Hu et al., 2021) and Argentina (Calderoli et al., 2017), while oxygen and iron affect the active diazotrophic lineages in some types of aquatic habitats (Moisander et al., 2014; Larson et al., 2018). However, the effect of temperature and other environmental factors on active diazotrophs in hot spring remains unclear.

There are thousands of hot springs with various hydrothermal characteristics in the Tengchong geothermal zone in Yunnan Province of southwestern China, one of the most active geothermal areas in the world (Du et al., 2005). Previous 16S rRNA gene-based studies showed that Tengchong hot springs were inhabited by highly diverse prokaryotic communities and harbored many potential diazotrophic groups (e.g., Synechococcus, Hydrogenobacter). Many of these groups are unique lineages that have never been detected in hot springs of other regions (Song et al., 2012; Jiao et al., 2020, 2021; Xian et al., 2020). One DNA-based study revealed that the nifH gene composition in two Thengchong hot springs has some common and unique phylogenetic lineages (Feng et al., 2018). However, there is still a lack of comprehensive census on the difference of active diazotrophs among different geothermal regions. Such knowledge gap limits our understanding of the distribution and adaptability of diazotrophs in global hot springs.

In this study, the expression and diversity of nifH and 16S rRNA genes were investigated in geothermal springs with a wide range of temperature and pH in Tengchong. An integrated approach including high-throughput sequencing and geochemical analyses were employed. The main aims of this study are to investigate (1) the composition and distribution of active diazotrophs in Tengchong geothermal springs; (2) the differences of diazotrophic lineages among different geothermal regions; and (3) the influences of environmental variables on the active diazotrophic communities.

Materials and methods

Field measurements and sampling

Sampling cruise was performed in Aug 2014. A total of 20 hot springs in Tengchong were selected for field measurements and sampling (Supplementary Table 1 and Supplementary Figure 1). At each hot spring, field measurement and sampling were performed during 9 am-13 pm every day. In the field, water temperature and pH were determined using a Hach pH meter equipped with a pH probe and a temperature probe (PT-10, SARTORIUS, Germany). Nitrite, nitrate, ammonium and sulfide were measured using Hach kits (model CEL 850/product #: 2687900, Hach Chemical Co., Iowa, United States) according to the manufacturer’s instructions. After field measurement, streamers, mat-containing sinters and sediments in different sites of each sampled hot spring were collected and mixed, then were put into 2 ml centrifuge tubes and 15 ml Falcon tubes, respectively. The sample tubes were immediately stored in liquid nitrogen. The samples were kept in liquid nitrogen in the field, during transportation and stored in the laboratory until further RNA extraction.

RNA extraction and cDNA synthesis

Environmental RNA was extracted from 1 to 3 gram of sediment, mat or sinter samples using the PowerSoil RNA Isolation Kit (Mo BIO Laboratories, Carlsbad, CA, United States). To remove the residual DNA, the soluble crude RNA was digested with RNase-free DNase I (Fermentas, United States). The DNase-digested RNA samples were checked for potential genomic DNA contamination by PCR amplification with specific primer sets of 16S rRNA and nifH genes (see below for primer information).

The checked RNA samples were quantified using Nanodrop UV-Vis spectrophotometer (ThermoFisher, United States) and reverse transcribed into cDNA using the Fermentas AMV Reverse Transcriptase (Fermentas, United States). Absence of potential contamination from DNA and chemical reagents was verified by conducting the same reactions without the AMV reverse transcriptase and template, respectively.

Amplification and sequencing of 16S rRNA and nifH genes

The resulting cDNA was employed as the template for amplification of polymerase chain reactions (PCR) of 16S rRNA and nifH genes. The hypervariable V4 region of the 16S rRNA genes were PCR amplified using primer sets of 515F (5′-GTG CCA GCM GCC GCG GTA A-3′) and 806R (5′-GGA CTA CHV GGG TWT CTA AT-3′) (Caporaso et al., 2011). The fragments of nifH genes were PCR amplified using specific primer sets of KAD3 (5′-ATH GTI GGI TGY GAY CCI AAR GCI GA-3′) and DVV (5′-ATI GCR AAI CCI CCR CAI ACI ACR TC-3′) (Gaby and Buckley, 2012). To pool multiple samples for one run of high-throughput sequencing, a sample tagging approach was used (Meyer et al., 2008). Each tag was added to the 5′ end of the reveres primer of 16S rRNA gene and the forward primer KAD3 of nifH gene, respectively. PCR amplification of 16S rRNA gene consisted of an initial denaturation at 95°C for 5 min, 25 cycles of denaturation at 95°C for 45 s, annealing at 54°C for 1 min, and extension at 72°C for 1 min, and a final extension at 72°C for 10 min. Individual reagents and their concentrations were as follows: 1 × PCR buffer with 1.5 mM Mg2Cl, dNTPs (100 M each), 0.25 M each primer, 2.5 U of DNA polymerase (Ex-Taq) (TaKaRa, Dalian, China) and 1ul of cDNA. The PCR amplification of nifH was 35 cycles and annealing temperature was 56°C, and the other protocols were same as the 16S rRNA gene. PCR products were purified using a QIAquick Gel Extraction Kit (QIAGEN, Germany). The purified PCR products were annealed to oligonucleotides that are complementary to an adaptor sequence. Sample libraries were generated from purified PCR products, and subsequently were pooled in equimolar concentrations and added library specific sequencing adapters by NEBNext Ultra (NEB#e7370S/L) assay as followed by instruction, and dual index sequencing of paired-end 250 bp was run on an Illumina Hiseq2500 instrument (Illu- mina, San Diego, CA, United States). The sequencing process was completed in Fixgene Theh Co., Ltd. (Beijing, China).

Forward and reverse reads were merged into full-length sequences using FLASH (Magoč and Salzberg, 2011). The raw sequence data were evaluated and filtered to ensure that >80% of the base calls in a sequence had a Phred quality score of 20 using the FASTQ Quality Filter (q = 20, p = 0.8). Then the barcode and primer sequences were deleted using FLEXBAR (Roehr et al., 2017). Sequences were removed if they were below certain length (16S rRNA gene < 200bp, nifH gene < 250bp) or contained ambiguous bases. Chimeric sequences were identified and discarded using the USEARCH software (Edgar et al., 2011).

The analysis of 16S rRNA gene

A total of 17,000 sequences of 16S rRNA gene were extracted randomly from the high-quality sequence of each sample for the following analysis. The operational taxonomic units (OTUs) of 16S rRNA gene were classified using UCLUST at the 97% similarity level. OTUs that contain less than two sequences were removed. Sequences were taxonomically classified by the Silva database (db128) using RDP algorithm (60% threshold) (Quast et al., 2013). Potential contamination sequences (e.g., unknown, chloroplast, mitochondria, etc.) were removed. The alpha diversity estimates, including coverage (Jiang et al., 2006), Shannon-Wiener’s diversity index (Spellerberg and Fedor, 2003) and Simpson’s diversity index (He and Hu, 2005) were calculated by using Mothur software (Schloss et al., 2009). The nifH gene sequences (a total of 6,243 sequences) that can be identified as the main taxa at the genus level in GenBank were selected, with no terminator in the sequence and a length greater than 600bp as the standard. These selected sequences were employed to construct a local database including potential nitrogen fixation activities, based on the classification of these sequences. By Blasting against the constructed local database, the genus with potential nitrogen fixation activity in 16Sr RNA Gene Data in this study was determined.

The PICRUSt (Langille et al., 2013) and Piohillin (Iwai et al., 2016) (with 99% identity cut-off) were applied to predict the abundance of nifH gene-related 16S rRNA-derived OTUs. For the PICRUSt analysis, the OTU table was normalized by dividing each OTU abundance by the number of the known or predicted 16S rRNA gene copies. Subsequently, the weighted NSTI score was used to measure the phylogenetic distance between OTUs and reference OTUs in each sample, and to evaluate the prediction accuracy of PICRUSt. The final functions for metagenome were predicted with the script predict_ metagenomes.py, generating a table of Kyoto Encyclopedia of Genes and Genomes (KEGG) Orthologs (KOs). Thirdly, the predicted metagenomes were classified into higher categories (e.g., KOs into KEGG pathways). The Piohillin analysis was performed on the website[1] according to the pipeline, and the reference database was used the “BioCyc”.

The analysis of nifH gene

The nifH gene sequences were translated into amino acid sequences by using RDP Framebot[2], and the sequences containing in-frame stop codon were excluded. Subsequently, the remaining sequences were used for OTU analysis with 95% similarity as the cutoff value by using USEARCH software. The alpha diversity estimates, including coverage, Shannon-Wiener’s diversity index, and Simpson’s index of diversity, were calculated by using Mothur software (Schloss et al., 2009). The beta diversity analysis was conducted using QIIME2 (Bolyen et al., 2019), with the script core_diversity_analyses.py. Jackknifed weighted pair group method with arithmetic mean (UPGMA) clustering was performed to compare microbial community similarity among the studied samples based on UniFrac weighted phylogenetic distances.

The obtained nifH gene sequences were classified into OTUs by the phylogenetic analysis combined with representative sequences of known diazotrophic lineages (refer to BLAST results in Non-Redundant Protein Sequence Database and the papers of Zehr et al., 2003 and Raymond et al., 2004). According to the published accession numbers, the nifH gene sequences isolated from hot springs of other regions were downloaded and clustered into OTUs at the 95% similarity level. The sequences of nifH gene OTUs from this study and other regions were combined for the phylogenetic analysis. The Neighbor-joining phylogenetic trees were constructed from dissimilar distance and pairwise comparisons with the Jukes-Cantor distance model using the MEGA7.0 (molecular evolutionary genetics analysis) program. Bootstrap replications of 1,000 were assessed. The phylogenetic tree visualization was created by using the iTOL web server (Letunic and Bork, 2019).

To identify the differences in environmental background of hot springs, principal component analysis (PCA) of physicochemical variables was performed by using R software. Spearman’s correlation and redundancy analyses (RDA) were also performed to explore relationships of microbial community patterns with physicochemical variables. The significance of the RDA models and the explanatory factors were tested by using 999 permutations. The ordination on the x- and y-axis and the length of the corresponding arrows indicate the relative importance.

Nucleotide sequence accession numbers

All the nifH and 16S rRNA gene sequences obtained from this study have been deposited at the National Omics Data Encyclopedia (NODE)[3] under the project OEP003492 with run accession number OER258273 to OER258312.

Results

Sample description

The studied hot springs exhibited a wide range of physical and chemical conditions (Table 1). The temperature of the studied hot springs ranged 34-84°C with pH ranging 4.3-9.1; The concentrations of nitrite, nitrate and ammonium were below 8 mg/L (Table 1). Principal component analysis (PCA) and the Kruskal-Wallis test (p < 0.001) showed that the physicochemical variables differed significantly among the sampled hot springs (Supplementary Figure 2).

TABLE 1

Physicochemical variables of the investigated hot springs.

Sample name	Sample code	Temp (°C)	pH	NO2– (mg/L)	NO3– (mg/L)	NH3+ (mg/L)	S2– (mg/L)
Hashang#1	Hs1	84	8.2	0.005	0.010	0.100	16.000
Shuirebaozha#2	Srbz2	84	8.2	BD	BD	0.080	7.000
Shuirebaozha#3	Srbz3	84	9.1	BD	0.010	0.230	17.000
Hamazui	Hmz	74	7.8	0.003	0.010	0.120	14.000
Hashang#2	Hs2	73	8.2	0.001	0.010	0.370	4.700
Huanxiquan#2	Hxq2	72	6.7	0.017	BD	0.090	21.000
Rehaitiyanqu#5	Rhtyq5	70	4.3	0.004	BD	BD	20.900
Xianrendong	Xrd	68	7.5	0.019	BD	0.020	0.330
Hashang#3	Hs3	66	8.2	0.013	0.017	0.070	0.150
Rehaitiyanqu#2	Rhtyq2	62	5.6	0.004	BD	0.010	24.200
Shizitou#1	Szt1	62	8.4	0.024	BD	5.460	0.490
Jinzequan	Jzq	61	7.2	BD	0.003	0.230	BD
Shuirebaozha#1	Srbz1	61	7	BD	0.013	0.100	0.127
Heinitan#2	Hnt2	56	7	0.007	BD	0.230	0.110
Jiemingquan	Jmq	56	7.8	0.011	0.030	0.160	7.100
Heinitan#3	Hnt3	54	6.8	0.051	0.020	1.260	0.041
Heinitan#1	Hnt1	53	7	0.007	BD	4.030	0.003
Huanxiquan#1	Hxq1	53	7	0.012	BD	7.700	0.045
Hashang#5	Hs5	52	7.2	0.022	BD	0.910	0.081
Shuirebaozha#4	Srbz4	34	6.7	0.003	0.020	0.340	BD

BD: below the detection limit (0.001 mg/L).

Potential diazotrophic taxa based on 16S rRNA gene transcripts

A total of 67 known prokaryotic phyla were identified, and only about 20 of them were dominant (their abundance accounted for above 5% of the total 16S rRNA gene sequences of each sample) in one hot spring. Chloroflex, Proteobacteria, Aquificae, Cyanobacteria, Nitrospirae and, Acetothermia were the most predominant phyla, and they dominated in more than six samples (Supplementary Figure 3). A total of 118 genus-level taxa with potential nitrogen fixation capacity were detected, and the abundance of 19 taxa was greater than 1% in more than one hot spring (Supplementary Figure 4). These identified taxa accounted for 80.7-99.9% of nitrogen-fixing organisms of the studied samples, and belonged to six phyla, i.e., Cyanobacteria, Aquificae, Nitrospirae, Chloroflexi and Proteobacteria (Figure 1A). Aquificae dominated in hot springs Srbz2, Hs1, Hxq2, Rhtyq5, Rhtyq2, Srbz3 (temperature: 62-84°C, pH: 4.3-9.1), accounting for 12-85% and 84-99% of the total 16S rRNA gene sequences and diazotrophic 16S rRNA gene sequences, respectively (Figure 1A). Hydrogenobacter is the only known genus with nitrogen fixation potential detected in Aquificae (Supplementary Figure 4). Cyanobacteria with photosynthetic function are dominant in Hs3, Hnt2, Hnt1, Jmq, Srbz, and Srbz4 hot springs (temperature: 34-66°C, pH: 7.0-8.2), accounting for 12-49% and 55-97% of total 16S rRNA gene sequences and diazotrophic 16S rRNA gene sequences, respectively (Figure 1A). The unclassified genus belonging to “Family I of subsection III” was one dominant cyanobacterial group, accounting for more than 10% of total 16S rRNA gene sequences in Hnt2, Jmq, and Srbz1 hot springs. Unclassified lineage belonging to “Family I of subsection I” was predominant cyanobacterial group in Srbz4 hot spring, accounting for 49% of the total 16S rRNA gene sequences (Supplementary Figure 4). Thermosynechococcus was the main genus of Cyanobacteria in Hs3, Hnt1, and Hnt2 hot springs, accounting for 18, 12, and 14.8% of total 16S rRNA gene sequences, respectively. Nitrospirae was the main potential diazotrophs (46 to 80% of the total diazotrophic 16S rRNA gene sequences) in Hxq1, Jzq, Hs5, Hnt3, and Xrd hot springs. The nitrogen fixation-associated Nitrospirae genus comprised of Nitrospira and Thermodesulfovibrio, with the former accounting for 46 and 71% of total 16S rRNA gene sequences in Hs5 and Hxq1 hot springs, respectively, and the latter accounting for 61, 77, and 80% of total 16S rRNA gene sequences in Hnt3, Jzq and Xrd hot springs, respectively (Supplementary Figure 4). Proteobacteria was detected in all the studied hot springs and covered 54 potential diazotrophic genera. However, only Geobacter, Desulfonema and Desulfobulbus belonging to Deltaproteobacteria and Pseudoalteromonas belonging to Gammaporteobacteria accounted for more than1% prokaryotic communities in more than one hot spring. The genera of Geobacter and Pseudoalteromonas accounted for 42% and 8% of diazotrophic community of Szt1 hot spring, respectively (Supplementary Figure 4).

FIGURE 1

The predictive analysis of nitrogen fixing functions at the level of phylum, enzyme and gene based on the 16S rRNA gene data. Panel (A) showed the compositions of predicted potential diazotrophic taxa; Panel (B) showed the abundance of predicted nitrogenase and nifH gene. Note: Piphillin and PICRUSt were applied for the predictions of nitrogenase and nifH gene, respectively.

Prediction analysis of potential nitrogen fixation function based on 16S rRNA genes

The PICRUSt analysis showed the NSTI score of each sample ranged from 0.05 to 0.28 (Supplementary Table 2). This ensures the predicted results can provide some important insights into microbial functions according to previous studies (Langille et al., 2013; Wang et al., 2016; Fan et al., 2017). The abundance of predicted nitrogenase-related genes was higher than that of nifH gene in most of the studied hot springs, which may be due to different prediction methods. The relative abundance of nifH gene, nitrogenase, and the nitrogen-fixing taxa showed consistent variation trends (Figure 1B). The abundance of nitrogen-fixing taxa was correlated with the abundance of nifH gene (r = 0.69, p < 0.001) and nitrogenase (r = 0.50, p < 0.001), indicating that the analysis of diazotrophic community diversity predicted on the basis of 16S rRNA gene sequence was reliable.

Phylogeny and composition of expressive nifH gene

A total of 1,384,410 high-quality nifH gene sequences were retrieved from the cDNA samples, and they can be classified into 1,587 OTUs at the 95% protein sequence similarity level (Bae et al., 2018). The nifH gene amplification with the DNase-treated RNA yielded negative results, which proved that the obtained nifH gene applications were derived from cDNA transcripts rather than DNA contamination. However, due to the highly sensitivity of high-throughput sequencing approaches (Sogin et al., 2006), the possibility of DNA source cannot be absolutely excluded for some OTUs with very low abundance. Thus, 230 OTUs (accounting for 95.3- 99.8% of the total nifH gene sequences in each hot spring) whose abundance higher than 0.1% in at least one hot spring were remained for the downstream phylogenetic analysis. The obtained nifH gene OTUs were divided into five major groups (I-V) (the phylogenetic types referred to Zehr et al., 2003 and Raymond et al., 2004) that comprised of thirteen lineages. The OTUs belonging to Cyanobacteria, Aquificae, Nitrospirae, and Chloroflexi formed monophyletic clade, respectively, while those affiliated with Proteobacteria, Firmicutes, and Methanobacteria were distributed across on the phylogenetic tree (Figure 2A, see detailed phylogeny in Supplementary Figure 5).

FIGURE 2

Phylogenetic analysis, UPGMA analysis and community composition of nifH gene in the investigated hot springs. (A) Neighbor-joining tree showing the phylogenetic affiliation of the OTUs based on protein sequences. Scale bars indicate the Jukes-Cantor distances. Bootstrap values of 50% (for 1,000 iterations) are shown; (B) UPGMA cluster tree based on Bray-Curtis dissimilarity. (C) Compositions of nifH gene. The numbers in the parentheses are the temperature and pH of each site, respectively.

Jackknifed UPGMA clustering analysis showed that the diazotrophic communities could be divided into the following types: (i) Nitrospira-like, accounting for 67-96% of the sequences in Hxq, Jzq, Hs5, Hnt3, and Szt1 hot springs with temperature ranging 54-62°C; (ii) Aquificae-like, accounting for 38-85% of the sequences in Srbz2, Hs1, Hxq2, Rhyq2, Rhyq5, Xrd, Srbz3, Hs2, and Hs3 hot springs with temperature ranging 62-84°C; (iii) Cyanobacteria-like, accounting for 33-97% of the sequences of Hnt2, Hnt1, Jmq, Srbz1, and Srbz4 hot springs with temperature ranging 34-61°C; (iv) Nitrospira- and Methanobacteria-like (mainly in Hmz hot spring), belonging to group IV (Figures 2B,C). BLAST analyses revealed that the dominated Aquificae-and Nitrospirae-like nifH gene OTUs were closely related to the genera of Hydrogenobacter and Thermodesulfovibrio, respectively (Supplementary Table 3). In addition to the above widely existing dominant groups, there were also unique dominant nifH gene groups in some of the studied hot springs, such as those affiliated with “Beta/gamma-Proteobacteria” in Hxq1 hot spring, the Deltaproteobacteria nifH gene sequences (most of which cannot be assigned to any known genus) in Hnt3, Srbz1 and Hs1 hot springs, and the Chloroflexi-like nifH gene sequences (associated with the Roseiflexus genus) in Hs2 and Jmq hot springs. Furthermore, there were only one or two major OTUs (the relative abundance was above 40%) in most hot springs (Supplementary Figure 6), indicating the main executors of nitrogen fixation process were limited to a few lineages.

Phylogenetic analysis of nifH gene sequences in hot springs of different regions

In order to fully understand the phylogenetic relationships among the nifH gene sequences obtained in this study and those from hot springs of Yellowstone, Japan and Chile, phylogenetic analyses were performed by combining the representative nifH gene sequences from all of up-to-date known researches of diazotrophs in hot springs (more than 20 hot springs with temperature ranging from 35 to 84°C and pH 2-9; Steunou et al., 2006; Hall et al., 2008; Steunou et al., 2008; Hamilton et al., 2011a,b; Loiacono et al., 2012; Hamilton et al., 2014; Estrella-Alcamán et al., 2015; Alcamán-Arias et al., 2018; Nishihara et al., 2018b,c). The nifH gene lineages of hot springs were widely distributed in the phylogenetic tree. Many sequences from the same region were clustered in tight clades within many phylogenetic lineages such as Cyanobacteria, Nitrospirae, Proteobacteria and Aquificae (Figure 3). In addition, some clades represented large phylogenetic branches. For example, the group IV and “Unknown lineage” were composed of Yunnan sequences, and some sequences from Japan formed one monophyletic clade of Aquificae. This indicates that there are many unique nifH gene lineages in the geothermal springs of this region, and the distribution of diazotrophs has biogeographical patterns.

FIGURE 3

A Neighbor-joining tree showing the phylogenetic relationships between the nifH gene sequences obtained in this study and those from Yellowstone, Chile, Japan hot springs.

Influence of environmental factors on the diazotrophic communities

The redundancy analyses indicated that the diazotrophic communities were influenced by temperature, ammonia, nitrite, nitrate and sulfide in the investigated hot springs. Temperature was significantly correlated with the relative abundance of most phylogenetic lineages of nifH gene (Figure 4). Spearman’s correlation analysis showed that the Aquificae was positively correlated with temperature and sulfide, but negatively correlated with ammonia. Cyanobacteria was negatively correlated with temperature. At the OTU level, many of the nifH gene lineages were significantly correlated with temperature, ammonia, and sulfide (Figure 5).

FIGURE 4

RDA of phylogenetic lineages of nifH gene in relation to physicochemical variables. The percent variability explained by each principal component is shown in parentheses in the axis labels.

FIGURE 5

Spearman correlation between the measured environmental variables and the relative abundance of the nifH gene lineages (A), nifH OTUs (B) and the 16S rRNA gene of potential diazotrophic genera (C). ** p < 0.01 and * p < 0.05.

Discussion

Diazotrophic diversity in the Tengchong hot springs and comparative analysis with hot springs in other regions

The results of this study proved that active diazotrophs are widely distributed in hot springs, and many diazotrophic lineages have potential nitrogen fixation activity. Although these diazotrophic lineages have been detected by DNA-based methods, they have not been verified by in situ expression in hot springs. Furthermore, most of the investigated diazotrophic communities were dominated by the lineages of Aquificae and Nitrospirae. However, to our best knowledge, the dominance of these lineages has not been reported as the main diazotrophs in other ecosystems. Therefore, it can be referred that hot springs harbored unique active diazotrophs in comparison with other type of ecosystems. Furthermore, the phylogenetic analysis showed remarkable different diazotrophic genetic lineages between hot springs of Tengchong and other regions.

The nifH gene of Aquificae are mainly derived from Hydrogenobacter and Thermocrinis (Hamilton et al., 2011a,2014; Loiacono et al., 2012), which usually occupy similar niches in high temperature, circumneutral to alkaline terrestrial ecosystems (Eder and Huber, 2002; Reysenbach et al., 2005). YNP hot springs were inhabited by these two genera, but Tengchong hot springs only harbored Hydrogenobacter (Supplementary Table 3 and Supplementary Figure 3). The results about 16S rRNA gene in this and previous studies (Song et al., 2012; Hou et al., 2013) also showed the absence of Thermocrinis species in Tengchong springs, suggesting the distribution of this diazotrophic genus might have geographical patterns. Interestingly, a large number of 16S rRNA gene sequences belonging to the genus Thermocrinis were found in Nakabusa hot springs, Japan, but phylogenetic analysis showed that any nifH gene sequences in those hot springs cannot be assigned to Thermocrinis (Figure 3; Nishihara et al., 2018b). This suggests that Thermocrinis-affiliated bacteria may not have nitrogen fixation ability in some hot springs. This phenomenon can be explained by the fact that some Thermocrinis species living in YNP hot springs may obtain nifH gene from coexisting diazotrophs such as Hydrogenobacter through horizontal gene transfer that was speculated to occur among diazotrophs (Koirala and Brözel, 2021) or prokaryotes in hot springs (Zhu et al., 2019; Li et al., 2021; Qi et al., 2021). In addition, there was one unique Aquificae-related nifH gene lineage in Japan hot springs (Figure 3), but those related nifH gene sequences have stop codons in their amino acid sequences (Nishihara et al., 2018b), indicating they are pseudogenes. So it is reasonable that Thermocrinis did not exist or were not expressed in Tengchong hot springs.

The composition of Nitrospirae diazotrophs differed between Tengchong and YNP hot springs. Nitrospirae diazotrophs were mainly composed of Thermodesulfovibrio genus in YNP hot springs (Loiacono et al., 2012; Hamilton et al., 2014; Thiel et al., 2017). In contrast, besides Thermodesulfovibrio genus, another Nitrospirae diazotrophic lineage belonging to novel genus Dissulfurispira was present in Tengchong hot springs (Figure 3). Dissulfurispira were recently isolated from Nakabusa hot springs of Japan (Umezawa et al., 2021), and have complete nitrogenase genes (nifBDEHKX, genome accession: NZ_AP022873). However, Dissulfurispira was not detected in the YNP hot springs. In addition, it is worth noting that the nifH gene transcripts in five Tengchong hot springs (Hxq, Jzq, Hs5, Hnt3 and Szt) are dominated by Nitrospira, indicating that Nitrospira are major active diazotrophs in Tengchong hot springs. However, this phenomenon has never been reported in other regions (Loiacono et al., 2012; Hamilton et al., 2014; Alcamán-Arias et al., 2018). Such difference may be caused by the different community features of the studied hot springs between Tengchong and other regions (Hou et al., 2013). Tengchong hot springs were predominated by Nitrospira nifH gene, representing low-temperature and non-phototrophic type of communities (the temperature ranged from 54 to 62 °C; 16S rRNA gene analysis showed that the abundance of phototrophic taxa was below 6%, Supplementary Figure 3). In comparison, Yellowstone and Chile hot springs that have been studied represent either high-temperature chemotrophic or low-temperature phototrophic communities (Loiacono et al., 2012; Hamilton et al., 2014; Alcamán-Arias et al., 2018), in which the Nitrospira diazotrophs have weaker adaptability than those coexisting lineages, such as Aquificae or Cyanobacteria.

A certain amount of nifH gene sequences were detected to be related to Roseiflexus, a member of Chloroflexi, in some samples. Combined with the results of in situ studies on YNP hot springs (Loiacono et al., 2012; Klatt et al., 2013), the expression of nifH gene of Roseiflexus seems to be widespread in hot springs. However, previous studies have pointed out that the lack of nifEN gene encoding nitrogenase cofactor assembly protein in Roseiflexus (Dos Santos et al., 2012) leads to the unclear understanding of whether the nifH gene of Roseiflexus spp. encodes functional nitrogenase components. It is believed that “genome downsizing” is one of the important strategies for thermophiles in hot springs to adapt to high temperature (Thomas et al., 2019; Alcorta et al., 2020). So Klatt et al. speculated that Roseiflexus species living in hot springs may directly assemble cofactors on the NifD/NifK heterodimers to offset the loss of nifEN gene function (Klatt et al., 2013). Moreover, nitrogenase activity has been demonstrated in organisms lacking nifEN gene in Elusimicrobia (Zheng et al., 2016). Therefore, the nifH gene expression of Roseiflexus indicates that the species of this genus has nitrogen fixation activity in hot springs.

Environmental factors influencing diazotrophs in the tengchong hot springs

In this study, nifH gene expression was obtained in Tengchong geothermal springs with a wide range of physicochemical parameters. Combined with the in situ expression results of YNP (Steunou et al., 2006; Steunou et al., 2008; Loiacono et al., 2012; Hamilton et al., 2014), Chile (Estrella-Alcamán et al., 2015) and Japanese (Nishihara et al., 2018a) hot springs, it can be seen that the occurrence of biological nitrogen fixation process may not be limited by temperature, pH or geographic location of hot springs. However, some environmental factors may influence the structure of diazotrophic communities in Tengchong hot springs.

Two previous studies demonstrated that temperature influenced the composition of diazotrophs at different locations with the same chemical background in hot spring outflowing channels (Loiacono et al., 2012; Feng et al., 2018). However, no investigation has been performed to assess the contribution of temperature to the shaping of diazotrophic communities among separated hot springs, which generally have different chemistry background. In the investigated samples, the abundance of diazotrophic clades changed regularly with temperature. For example, the Aquificae-related nifH gene transcripts predominated in the hot springs with relatively high temperature (above 62°C), while Nitrospira or Cyanobacteria dominated in the hot springs with temperature below 62°C (Figure 2). BLAST analysis showed that the dominant cyanobacterial nifH gene OTU “Top2” in high temperature samples (hot springs Hnt1 and Hnt2 with 53-56°C; Figure 2) was most closely related (96% identity) to those retrieved from YNP geothermal springs, while OTU “Top9” was predominant in low-temperature samples (the spring Srbz9 with 34°C), which was similar to the sequences from normal temperature lake environments (96% identity). In addition, the RDA (Figure 4) and Spearman’s correlations analysis (Figure 5) also revealed significant statistical correlation between temperature and composition of diazotrophs. Taken together, temperature represents a strong selective pressure driving the distribution of diazotrophic clades in Tengchong hot springs.

In addition, statistical analyses revealed the abundance of many nifH gene OTUs was related to the concentration of ammonia (Figure 5), suggesting the influence of available nitrogen sources on nitrogen-fixating microorganisms (Dixon and Kahn, 2004; Hong et al., 2017; Zilius et al., 2020). Especially, the Aquificae diazotrophs showed negative correlation with ammonia, indicating that the expression of Aquificae-related nifH gene increased with the decrease of ammonia concentration. Furthermore, as far as we known, S2–is rarely reported as a factor affecting the active diazotrophic clade in normal environments. However, in the present study, many diazotrophic lineages, especially the Aquificae, at OTU and phylum levels showed significant correlations with the concentration of S2–(Figure 5), indicating that S2– may affect the distribution of diazotrophs. Such correlations can be explained by the fact that sulfur compounds are common and abundant in hot springs, and they are important electron donors or acceptors and are utilized by many chemolithotrophic organisms to obtain energy (Amenabar and Boyd, 2019). For example, the Aquificae strains isolated from Tengchong hot springs can use S2–, S0, S2O32– as the electron donors (Hedlund et al., 2015), so it is reasonable to find positive correlation between S2–concentration and Aquificae nifH gene abundance (Figure 5). Alternatively, considering that S2– can be easily converted to S0 or S2O32– by biotical and/or abiotic catalysis in hot springs (Amenabar and Boyd, 2019), the positive correlation between S2–concentration and Aquificae nifH gene abundance suggested that the oxidation of S0 or/and S2O32– may be caused by Aquificae diazotrophs in the investigated hot springs. At present, it has been proved that the diazotrophic Aquificae strains isolated from Nakabusa hot spring in Japan can use S2O32– as an electron donor (Nishihara et al., 2018c).

Conclusion

Biological nitrogen fixation process widely exists in terrestrial geothermal springs. The nifH gene transcripts are related to many diazotrophic lineages. The major active diazotrophic lineages belong to Aquificae, Cyanobacteria and Nitrospirae. Aquificae- and Nitrospirae-like diazotrophs dominate in high- and low-temperature hot springs, respectively. In addition, some investigated springs harbored unique predominant clades. The diazotrophs in Tenchong hot springs exhibited significant difference in phylogenetic lineages from those of YNP and Japanese hot springs. Furthermore, temperature and concentrations of ammonia and sulfide may be important environmental factors influencing the expression and distribution of nifH genein hot springs. These results improve the understanding of the distribution and diversity of diazotrophs in hot spring environments.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary material.

Author contributions

Z-QS, W-JL, and HJ conceived the study. Z-QS, LW, and FL performed the on-site measurements and collected the samples. Z-QS, QZ, and DP analyzed the geochemistry and microbiology of the samples. Z-QS analyzed the sequencing data. Z-QS, W-JL, and HJ drafted the manuscript. All authors reviewed results and commented on the manuscript.