The sampling sites are under the Indo-Burma biodiversity hotspot, located near the China-Myanmar border (Figure 1) in Tengchong town (Dongshan, Huangpo, Rehai and Qinghai) and Longling town (Banglazhang), Yunnan Province, China. Sampling was conducted from 16 different sites including water, wet sediment, and microbial mats during August 2017, targeting especially cyanobacterial population and their isolation (Figure 1 and Table 1). Limnological features like pH, temperature (°C), and conductivity (μS/cm) were recorded using a portable pH meter equipped with a pH probe, a temperature probe, and a conductivity probe (Thermo Orion 420C-01A), respectively, on the spot. Microbial biofilms from various substrates were scrapped using sterile forceps and spatula into 50-ml Falcon tubes. Water samples were collected in pre-sterilized 1-L polyethylene bottles for chemical analysis. These samples were immediately stored on dry ice for transportation and later kept at a −80°C deep freezer in the laboratory until further analysis. Limnological parameters including concentration of cations (calcium, magnesium, potassium, ammonium, sodium, and lithium) and anions (nitrite, nitrate, phosphate, sulfate, bromide, chloride, and fluoride) were analyzed using Dionex ICS-5000 + ion-exchange chromatography system equipped with Dionex IonPac CS12A and Dionex IonPac CS18 column (250 × 4 mm, 5 μm, ThermoFisher, Waltham, MA, United States), respectively.

The next-generation sequencing library preparations and Illumina MiSeq sequencing were conducted at GENEWIZ, Inc. (Suzhou, China). DNA samples were quantified using a Qubit 2.0 Fluorometer (Invitrogen, Carlsbad, CA, United States). DNA (30–50 ng) was used to generate amplicons using a MetaVx™ Library Preparation kit (GENEWIZ, Inc., South Plainfield, NJ, United States). Primers targeting V3–V4 hypervariable regions of the 16S rRNA gene were proved to be able to capture the bacterial diversity in various environments (Klindworth et al., 2013). Thus, these regions were selected to generate amplicons and follow taxonomy analysis. These regions were amplified using primer pairs V3V4-F (CCTACGGRRBGCASCAGKVRVGAAT) and V3V4-R (GGACTACNVGGGTWTCTAATCC). First-round polymerase chain reaction (PCR) products were used as templates for second-round amplicon enrichment PCR. At the same time, indexed adapters were added to the ends of the 16S rRNA amplicons to generate indexed libraries ready for downstream NGS sequencing on Illumina Miseq. DNA libraries were validated by Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, United States), and quantified by Qubit 2.0 Fluorometer. DNA libraries were multiplexed and loaded on an Illumina MiSeq instrument according to the manufacturer’s instructions (Illumina, San Diego, CA, United States). Sequencing was performed using a 2 × 300 paired-end (PE) configuration; image analysis and base calling were conducted by the MiSeq Control Software (MCS) embedded in the MiSeq instrument. The QIIME data analysis package was used for 16S rRNA data analysis. The forward and reverse reads were joined and assigned to samples based on barcode and truncated by cutting off the barcode and primer sequences. Quality filtering on joined sequences was performed and sequences with length < 200 bp, ambiguous bases, and mean quality score ≥ 20 were discarded. Then, the sequences were compared with the reference database (RDP Gold database) using the UCHIME algorithm to detect chimeric sequences, and then the chimeric sequences were removed. The effective sequences were used in the final analysis. Sequences were grouped into operational taxonomic units (OTUs) using the clustering program VSEARCH (1.9.6) against the Silva 119 database pre-clustered at 97% sequence identity. The Ribosomal Database Program (RDP) classifier was used to assign taxonomic category to all OTUs at a confidence threshold of 0.8. The RDP classifier uses the Silva 123 database, which has taxonomic categories predicted to the species level.

The biological mats were washed thoroughly by double-distilled water to remove inorganic debris, followed by homogenization by a glass homogenizer (1 ml) and then inoculated into Petri plates containing solidified BG11 medium (Rippka et al., 1979) at a slightly higher temperature of up to 40 ± 2°C under 30 μmol photon m^–2 s^–1 continuous light intensity in a temperature-controlled culture room. Axenic cultures of cyanobacteria were obtained by repeated sub-culturing on Petri plates.

To verify the thermotolerant features of these isolated cyanobacterial species, they were incubated in 10-ml test tubes at high temperature (50°C) with initial OD₇₃₀ ∼0.06 under a continuous light intensity of 30 μmol photon m^–2 s^–1 without shaking.

Morphological features of the cyanobacteria were studied extensively under an Olympus CX31 microscope equipped with an Olympus DP-72 digital camera (Olympus, Tokyo, Japan) and software DP2-BSW (version 2.2) at 20 ms exposure time. The enumeration of phenotypic features was later considered for the identification of the taxa.

Genomic DNA of all the isolates was extracted by using HiPurA™ plant genomic DNA Miniprep Purification Spin Kit. 16S rRNA gene from the isolated DNA was amplified using cyanobacteria-specific primers CYA106F and CYA781R (Nübel et al., 1997). PCR amplification consisted of an initial denaturation at 95°C for 5 min, 32 cycles of denaturation at 95°C for 30 s, annealing at 55°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 MgCl₂, dNTPs (100 μM each), 0.25 μM each primer, 2.5 U of DNA polymerase (Ex-Taq) (TaKaRa, Japan), and ∼50 ng of total DNA. PCR products were purified using an E.Z.N.A. Gel Extraction Kit (Omega Bio-Tek, United States) according to the manufacturer’s instructions. The obtained sequences of 16S rRNA gene were sequenced and the ambiguities were removed by visual editing. The sequences were compared with the NCBI GenBank database using the BLASTn¹ and sequences of uncultured organisms were excluded.

The paired-end nucleotide sequences have been submitted in NCBI SRA database with the accession IDs SRR11796781–SRR11796799 and SAMN14918520–SAMN14918537 for BioSample accessions under the BioProject PRJNA628565. Sequences of 16S rRNA genes obtained during the present study were deposited in GenBank with accession numbers (MK625304–MK625325).

All OTUs belonging to cyanobacteria and 16S rRNA gene sequence of cyanobacterial isolates were compared with the NCBI (National Center for Biotechnology Information) database using BLAST. An initial BLAST search of the NCBI GenBank database provided candidate sequences for comparing relationships of the cyanobacterial isolates with previously characterized cyanobacterial species. Different sets of sequences were combined together to construct 3 types of phylogenetic trees in the present study, which are as follows: (1) the sequences of major OTUs of cyanobacteria/bacteria of corresponding sampling spots and 16S rRNA gene sequence of cyanobacterial isolates; (2) the sequences of all reported OTUs of cyanobacteria and 16S rRNA gene sequence of cyanobacterial isolates; and (3) the 16S rRNA gene sequences of other hot spring cyanobacteria reported worldwide and 16S rRNA gene sequence of cyanobacterial isolates of the present study. Multiple alignments were created with reference to the selected GenBank sequences and Neighbor-joining phylogenetic trees were constructed from dissimilar distance and pairwise comparisons with the Jukes–Cantor distance model using the MEGA (molecular evolutionary genetics analysis) program, version 5 (Tamura et al., 2011). Bootstrap replications of 1,000 were assessed.

Microbial mats in the sampled thermal springs were both submerged and floating on the surface (Figure 1). Sometimes they appear as gelatinous films blanketed over sandy bottoms of a few hot springs. The textures of the mats were thin slimy gelatinous films to thick leathery mats with color ranging from bluish to yellowish green. The hot springs displayed a wide range of limnological characteristics, i.e., moderate to high temperature (38–90°C), acidic to slightly alkaline pH (4.0–9.4), and a variable conductivity (4–1,953 μS/cm), though the springs are in close proximity (Table 1). Based on the water temperature, the sampling sites were segregated into four groups representing different temperature gradients, i.e., Groups 1 (30–50°C), 2 (50–70°C), 3 (70–80°C), and 4 (80–90°C). Thermal springs Ds001, Hp003, and Qs010 (Group 1) were marked with the lowest temperature range (30–50°C), while the highest temperature range (80–90°C) was documented for the springs Re006, Re007, Re008, Blz012, and Blz016 (Group 4). Springs with temperature range 50–70°C were listed for Re004, Re005, Blz011, and Blz014 (Group 2). Only two springs Blz013 and Blz015 (Group 3) had a temperature range 70–80°C.

A comparative ionic composition and conductivity of water samples are presented in supporting information (Supplementary Figure 1). The abundance of sodium ion (Na⁺), among the cations in all the water samples, was distinctly observed followed by the concentration of potassium (K⁺). Similarly, chloride and sulfate were the dominant anions. The variation in ionic concentrations and conductivity in Longling town hot springs was negligible, but significant in the hot springs of Tengchong town. Total ion analysis (cations and anions) showed its richness in water sample of Re008 and minimum in Hp003. It was also revealed that the conductivity was directly proportional to the concentration of ions present in water samples. However, the derivation of relatedness was not made for the sample Re008, due to insufficient water sampling.

A total of 584,226 reads (1787 OTUs) were identified from these samples. Operational taxonomic unit Venn or Petal diagram was also generated to illustrate the distribution of OTUs among these four groups (Supplementary Figure 2). (1) Group 1/30–50°C (1,186 OTUs): Ds001, Hp003, and Qs010; (2) Group 2/50–70°C (761 OTUs): Re004, Re005, Blz011, and Blz014; (3) Group 3/70–80°C (322 OTUs): Blz013 and Blz015, and (4) Group 4/80–90°C (713 OTUs): Re006, Re007, Re008, Blz012, and Blz016. The OTUs restricted to any individual group were 833, 103, 15, and 60 for Group 1, Group 2, Group 3, and Group 4, respectively. It also showed the information on the distribution of shared OTUs across these four groups, which were 260, 80, 253, 261, 541, and 268 between Groups 1 and 2, Groups 1 and 3, Groups 1 and 4, Groups 2 and 3, Groups 2 and 4, and Groups 3 and 4, respectively. Larger proportions of OTUs (179) were shared between Groups 2 and 4. There was relatively lower overlap between Groups 1 and 3. Common OTUs shared by all the groups was 58. Chao1 was used to compare species richness between different groups (Supplementary Figure 3), which supported that OTU estimate was significantly higher in Group 1, i.e., the lowest temperature gradient in all gradients.

Thirty major OTUs were represented as clustering heat maps (Supplementary Figure 4), which showed the color of different shades based on the relative abundance of each OTU. Among these, OTUs 1, 2, 5, 7, 12, and 21 belong to cyanobacteria, where OTUs 2 and 5 were related to morphological subsection III. Microbial community composition based on pyrosequencing of 16S rRNA gene (Supplementary Figure 5) revealed Cyanobacteria and Proteobacteria representing a larger portion of bacterial community in the hot spring mats of Yunnan Province. Other most common bacterial genera were Firmicutes, Chloroflexi, Bacteroidetes, Acidobacteria, and Deinococcus-thermus. Abundance of cyanobacteria was reported in all the sampling sites except RH007, which was rich with Proteobacteria (Supplementary Figure 5A). Furthermore, it can also be inferred that an extremely high temperature of 85°C and the acidic water (pH 4.0) of RH007 favored the growth of Proteobacteria, Bacteroidetes, and Actinobacteria, but inhibited the growth of cyanobacteria. Another site with the hottest (85°C) but alkaline water (pH 9.3), i.e., RH006, supported the growth of the members of phylum Thermodesulfobacteria and Aquificae rather than Cyanobacteria and Proteobacteria. Our result showed that the bacterial phylum Cyanobacteria and Proteobacteria were dominant in all the studied hot springs. The dominating organism in microbial mat depended on the temperature, pH, and geochemistry of sampling sites (Hou et al., 2013). Species richness of cyanobacteria was higher in the springs HP003, RH009, BLZ014, and BLZ015 than any other spots (Supplementary Figure 5A). The lowest number of OTUs was observed in Group 3, i.e., temperature range of 70–80°C, but interestingly with the highest proportion belonging to cyanobacteria (Supplementary Figure 5B).

This study on the microbial composition in these hot springs was intended to explore cyanobacterial diversity in more detail rather than dealing with genera of other microorganisms precisely. Out of the total OTUs, 100 OTUs belonging to cyanobacteria were identified. Ninety-four and 36 cyanobacterial OTUs were found in Tengchong and Longling town, respectively. The OTU Venn diagram was illustrated based on the distribution of cyanobacterial OTUs among the four groups of temperature range and three groups of pH range (Figure 2). It is quite clear that there were diversified cyanobacteria in hot springs with temperature 30–50°C (Figure 2A) and pH > 9 (Figure 2B). This was also supported by the Shannon–Weaver (Figure 3A) and Simpson index (Figure 3B) analysis. When grouped by region, the results show that the diversity of cyanobacteria in Tengchong has changed greatly, while the diversity in Longling town is relatively small (Figures 3C,D).

According to the abundance of cyanobacteria OTUs at each sample point, 30 OTU species with higher abundance were selected to make a heat map (Figure 4). The abundance of different cyanobacteria in each sampling point (Figure 4A) as well as in different temperature ranges was depicted (Figure 4B), which showed that the types of cyanobacteria were quite different. In particular, in the T30–50 group, the cyanobacteria were significantly different from the other three groups, and the abundance of cyanobacteria was higher. The other three groups of cyanobacteria had relatively similar types and lower abundance, which may be due to high water temperature and insufficient nutrients in the water.

The cyanobacterial OTUs found in the present study were assigned to 34 genera and 3 cyanobacterial groups (Synechococcus, Gloeocapsa, Aphanothece, Chroococcus, Cyanobium, Chondrocystis, Thermosynechococcus, Chroococcidiopsis, Pleurocapsa, Xenococcus, Gloeocapsopsis, Uncultured Chroococcales, Leptolyngbya, Oscillatoriales, Phormidium, Microcoleus, Leptolyngbyaceae cyanobacterium, Wilmottia, Arthronema, Pseudanabaena, Pseudanabaenaceae, Planktothricoides, Nodosilinea, Ancylothrix, Jaaginema, Oscillatoria, Microseira, Tolypothrix, Scytonema, Calothrix, Nostoc, Nodosilinea, Chlorogloeopsis, Stigonema, Mastigocladus, Iphinoe, and filamentous thermophilic cyanobacterium) approximately. More than 59% of total morphological subsection was represented by subsection III (59% of total subsections), especially the members of the genus Leptolyngbya. The second most abundant cyanobacterial group belonged to subsection II (18%) followed by subsections I (12%), V (6%), and IV (5%) (Figure 5A).

Floristic assessment of the mat-forming cyanobacteria and other associated algae through microscopy identified altogether 39 taxa from these thermal habitats of Yunnan Province (Supplementary Table 1). As stated, the thermal springs of Yunnan provide a wide range of habitats with diverse physicochemical parameters. The cyanobacteria in these extreme habitats thus get ample opportunity to adapt various environmental circumstances. To understand this, we tried to index the ecological specifications of the cyanobacterial species distributed in these thermal habitats. Cyanobacteria reported from Yunnan hot springs were categorized to adapt to various groups, i.e., moderate acidophilic, moderate alkaliphilic, mesothermophilic, thermophilic, and hyperthermophilic, depending on the pH and temperature of the hot springs from where they have reported/isolated (Table 2). Cyanobacteria were identified through extensive phenotypic characterizations (Figure 6 and Supplementary Figure 6). They displayed similar diversities at both species and generic levels. Leptolyngbya was the most diverse genus (8 taxa), followed by Fischerella (3 taxa). For more clarity in understanding the diversity, the taxa were grouped into different sub-sections to estimate their percentage of contribution to the biome (Figure 5B). It was revealed that the cyanobacterial species belonging to subsection III (47% of total cyanobacterial subsection diversity) were the major contributor in mat formation, which included taxa of Leptolyngbya, Thermoleptolyngbya, Phormidium, and Planktothrix, which were observed in 10 sampling sites out of the 16. Leptolyngbya dominated the hot springs of DS001, HP003–HP006, BLZ011, BLZ013, and BLZ014, whereas Thermoleptolyngbya dominated in HP006, HP007, and BLZ011. Likewise, Fischerella colonized the mats in QS010, BLZ012, and BLZ016. Synechococcus, Gloeocapsa, Chlorogloeopsis, and Anabaena were considered as the other associated mat-forming cyanobacterial species. It is noteworthy to mention that the sites with temperature > 80°C had only a few traces of cyanobacterial populations that were detected after extensive microscopy of the natural samples.

Operational taxonomic units (OTUs) belonging to genera Cyanothece, Aphanocapsa, Thermoleptolyngbya, Planktothrix, Limnothrix, Anabaena, and Fischerella were not identified in 16S rRNA sequencing, which may be due to insufficient samples/DNA quantity for high-throughput sequencing study as well as insufficient data in GenBank. The common genera in both 16S rRNA sequencing and phenotypic analysis were Synechococcus, Gloeocapsa, Aphanothece, Chroococcidiopsis, Leptolyngbya, Phormidium, Calothrix, and Chlorogloeopsis. Other cyanobacterial taxa reported in 16S rRNA sequencing analysis were not identified through microscopic observations because of their rarity in occurrence. We believe that the phenotypic characterization of any particular group of microorganisms, especially cyanobacteria, will impart more authenticity along with the 16S rRNA sequencing analysis in understanding the infrastructure of any microbial ecosystem. Further, the available genetic resource databases lack thorough taxonomic information including the identity of a species. Thus, we have adopted both methods together to reveal true cyanobacterial diversity in the hot springs of Yunnan Province, which was represented by a total of 45 cyanobacterial genera. The 16S rRNA sequencing analysis depicted the whole cyanobacterial diversity of the habitats with 100 different types of taxa, but their population and pattern of colonization were deeply elaborated after the phenotypic characterization, where the dominant mat-forming and associated taxa were segregated.

A total of 19 cyanobacterial strains belonging to 10 species and 6 genera were isolated to axenic cultures. These isolates will create the chance to understand the physiology of heat-tolerant cyanobacteria, screen for value-added products, and develop them as model microorganisms. They were closely observed for morphological alterations during their various growth phases, to facilitate more authentic identification. Their closest relatives in database (percentage similarity) and taxonomic assignment are presented in Table 3. Most of the isolates were identified up to species level except Leptolyngbya sp. NK 1-10, Leptolyngbya sp. NK 1-12, Anabaena sp. NK1-14, Fischerella sp. NK 1-16, Fischerella sp. NK 1-20, Fischerella sp. NK 1-21, and Leptolyngbya sp. NK 1-23. All these taxa were not identified up to species level, which could open the scope for identifying them as novel taxa with biotechnological potential in the near future.

With the aim of confirming the thermotolerant feature of these isolated strains and evaluating their survival capacity at higher temperature, they were incubated at OD₇₃₀ ∼0.06 in test tubes and cultivated statically at 50°C under 30 μmol photon m^–2 s^–1 light intensity for 10 days. Compared to their initial OD₇₃₀, most of the cultures, except for NK1-10 and NK1-25, exhibited increased cell densities under these rigorous laboratory culture conditions (Figure 7). The results indicated the survivability of our isolates even at very low cell density under high temperature and further highlighted the significance of hot spring cyanobacterial strain mining using integrative approach.

Phylogenetic relationships among the 16S rRNA gene sequences of isolated cyanobacterial species were analyzed along with the OTUs reported from our next-generation sequencing as well as other 16S rRNA gene sequences of related taxa. The OTUs obtained through next-generation sequencing were analyzed by NCBI BLAST search (one by one) and classified into different subsections (Supplementary Table 2). A total of 100 cyanobacterial OTUs were confirmed, which mainly belonged to subsection III and Leptolyngbya genera. Major OTUs of corresponding sampling spots and 16S rRNA gene sequences of isolated cyanobacteria species from the hot springs of Yunnan Province were compared, which showed that isolated cyanobacterial species belonged to same subsection as major cyanobacterial OTUs of the corresponding sampling site (Supplementary Table 3).

Different dendrograms of isolated cyanobacterial species were constructed using the following information: (1) major OTUs of cyanobacteria and microorganisms reported from 16S rRNA sequencing analysis (Supplementary Figure 7); (2) all OTUs of cyanobacteria reported in 16S rRNA sequencing (Supplementary Figure 8); and (3) already reported cyanobacteria isolated from hot springs (Figure 8).

The sequences of major OTUs of corresponding sampling spots and 16S rRNA gene sequences of isolated cyanobacteria species from the hot springs of Yunnan Province were combined and Neighbor-joining phylogenetic trees were constructed (Supplementary Figure 7). The result showed that OTUs from 16S rRNA sequencing and isolated cyanobacterial species, belonging to the same subsection, clustered together. A major clade belonged to subsection III followed by subsection V. The OTUs 18, 33, 54, and 722 belonged to subsection III and clustered with different species of Leptolyngbya, Thermoleptolyngbya, and Planktothrix (isolated organism of the present study). Chlorogloeopsis fritschii NK1-24 was clustered with OTU1619, which showed the BLAST search similarity with C. fritschii (96%). Our isolated Fischerella sp. NK1-16, 20, and 21 clustered together with OTU298 under subsection V. OTUs 2 and 5 were the dominating cyanobacterial OTUs, which were classified as Uncultured Oscillatoriales and Phormidium autumnale based on the NCBI BLAST search, respectively. They clustered together and separated from the main subsection III clade. Other major OTUs belonging to bacterial taxa (OTU3, OTU4, OTU7, OTU8, OTU20, and OTU27) made a separate clade in the dendrogram.

Phylogenetic relationships between the isolated cyanobacterial species and all OTUs (100 OTUs) obtained through 16S rRNA sequencing of hot springs of Yunnan Province, China showed the clustering of isolated cyanobacteria in two major clades (Supplementary Figure 8). These two clades belong to subsections III, IV, and V. There were many other clades belonging to subsections I, II, and III, but showed separate clustering from our isolated organisms.

The 16S rRNA gene sequences of the isolated cyanobacteria were from BLAST in NCBI and sequences of other related taxa reported from hot springs were retrieved for the construction of a dendrogram (Figure 8). Phylogenetic analyses revealed that subsection III cyanobacteria were clustered in two clades. The first major clade of subsection III comprised different species of genera Leptolyngbya, Oscillatoriales cyanobacterium, Geitlerinema, Thermoleptolyngbya, and Pseudanabaena, while the second clade of subsection III included only Planktothrix and Phormidium. Cyanobacteria under Subsection V, i.e., Fischerella, Mastigocladus, Hapalosiphon, Westiellopsis, and Nostochopsis, were clustered together, while C. fritschii made a separate clade. Cyanobacterial members of subsection IV formed monophyletic clade including Anabaena and Nostoc. Subsection I clustered in different clades comprising mainly Gloeocapsa, Nodularia, and other hot spring cyanobacteria.