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.

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.

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 Mg₂Cl, 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).

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¹ according to the pipeline, and the reference database was used the “BioCyc”.

The nifH gene sequences were translated into amino acid sequences by using RDP Framebot², 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.

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

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).

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).

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.

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).

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.

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.

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).

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.

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, S^2–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 S^2–(Figure 5), indicating that S^2– 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 S^2–, S⁰, S₂O₃^2– as the electron donors (Hedlund et al., 2015), so it is reasonable to find positive correlation between S^2–concentration and Aquificae nifH gene abundance (Figure 5). Alternatively, considering that S^2– can be easily converted to S⁰ or S₂O₃^2– by biotical and/or abiotic catalysis in hot springs (Amenabar and Boyd, 2019), the positive correlation between S^2–concentration and Aquificae nifH gene abundance suggested that the oxidation of S⁰ or/and S₂O₃^2– 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 S₂O₃^2– as an electron donor (Nishihara et al., 2018c).