The hot springs of Anhoni and Tattapani are located in the margins of Gondwana coalfields of India and have been reported as important geothermal resources (Pandey and Negi, 1995). The presence of volcanic tuffs and interlayer basic silts, which are flat intrusions of igneous rocks formed between the pre-existing layers of rocks, has also been reported under these regions. The volcanic sediments are known to constitute hydrocarbons and other organic compounds and also release organic gases (Farooqui et al., 2009). The long-term seepage of hydrocarbons, either as macroseepage or as microseepage can bring about a diverse array of mineralogical and chemical changes that are favored by the development of near-surface oxidation or reduction zones (Khan and Jacobson, 2008). Surface geochemical methods are based on the premise that hydrocarbon gas components migrate from sub-surface petroleum reserves through faults and fractures and leave their signatures in near surface soils (Price, 1986; Tedesco, 2012). Petroleum microseepage in the soil causes several chemical reactions modifying the oxidation-reduction potential of the soil which plays an important role in the mobility of elements. It can be inferred that in Anhoni, a lot of igneous bodies cut the sedimentary bodies and water present deep along these joints gets heated from the mantle thus, finding space to come to the surface. As they travel up and cool they start to deposit the elements with lower solubility i.e., the hydrocarbons, progressively. The same water continues to cycle in this process and brings elements associated with hot springs from the mantle to the crust, making it enriched in hydrocarbons and other inflammable compounds.

Anomalous amounts of vanadium, chromium, nickel, cobalt, manganese, mercury, copper, molybdenum, uranium, zinc, lead,

and zirconium are positive indicators of petroleum deposits (Duchscherer, 1983). The migrating hydrocarbons create a reducing environment in the soil and subsurface, which increases the solubility of many trace and major elements (Mongenot et al., 1996). The samples from Anhoni have anomalously high concentrations of V (120.42–267.59 ppb), Cr (53.77– 189.14 ppb), Co (733.48–2882.70), Mn (83.51–255.74 ppb), Hg (333.00–555.00 ppb), Cu (102.9–276.67), and Mo (72.55–171.66) compared to regions (Bowen, 1979). This strongly indicates the presence of microseepage of petroleum deposits through the Narmada-Son lineament fracture zone into the near surface soils at Anhoni. However, this hypothesis needs further studies, which is beyond the scope of the current study.

The Tattapani geothermal field at Chhattisgarh consists of several hot springs having a broad temperature range from 52 to 98◦C and spread over an area of around 0.5 km2 (Sarolkar, 2010). The temperature of different reservoirs at this site has been reported to be as high as 230 ± 40◦C at a depth of 2 km and 112 ± 30◦C at a depth of 1 km (Vaidya et al., 2015). In another study, the chemical analysis of water from this site was found to contain moderate chloride, sulfate, silica, and sodium content, followed by low potassium, calcium, and arsenic content (Sarolkar and Das, 2006). This geothermal site has been reported majorly for its extremely high temperature, due to which a lower microbial diversity and a distinct functional profile was observed in the results of this study.

Phylum abundance results obtained from the amplicon sequencing data showed phylum Proteobacteria to be abundant in almost all the samples. Proteobacteria has also been reported from many studies based on the 16S rRNA analysis of hot springs with moderately high and very high temperatures (44– 110◦C) at various geographical locations, including India (Bowen De Leon et al., 2013; Chan et al., 2015; Ghelani et al., 2015). Since, Proteobacteria have been found in other hot spring studies including Indian hot springs, and was also one of the abundant taxa in this study, it appears to be indigenous to this region. Other phyla, such as Thermi and Thermotogae were also found abundant at both the sites. Phylum Thermotogae observed as abundant in Tattapani is found associated with hot springs of high temperature and is reported majorly for its high heat tolerance (Chan et al., 2015; Kanoksilapatham et al., 2015). Tepidimonas sp., belonging to phylum Proteobacteria observed as abundant in BAN, CAN, TAT-2, TAT-4, and CAP samples, has been isolated and sequenced from Chhoti Anhoni in an earlier study (Dhakan et al., 2016). The draft genome construction and analysis of the Tepidomonas taiwanensis genome revealed the presence of genes for sulfur metabolism, ammonia metabolism, nitrogen fixation, assimilation of organic acids, and a large variety of proteases. Gulbenkiania mobilis, another Proteobacteria was also cultured from the same environment and its first draft genome was reported in another study (Saxena et al., 2015). Though culturable, this genome was not among the most abundant genomes found in this study. The genome analysis of Gulbenkiania mobilis revealed its unique sulfur-metabolizing properties.

The metagenomic analysis of Anhoni region revealed enrichment of hydrocarbon degrading microbes, enzymes, and pathways. The presence of pathways such as benzoate, toluene, and xylene degradation at Anhoni, specifically in CAP sample, underscores the potential of the inherent microbial community to metabolize hydrocarbons as a source of energy. P. stutzeri and Acidovorax sp. which are known to harbor the enzymes for hydrocarbon degradation were found abundant at Anhoni. P. stutzeri is among the important alkane-degrading microorganisms reported for the bioremediation of crude oil, oil derivatives and aliphatic hydrocarbons (Lalucat et al., 2006). Potential degradation of toxic hydrocarbons such as, benzene, xylene and benzoate by P. stutzeri are well-reported in the literature and is also observed in the results of this study. Acidovorax sp. which was among the abundant species in CAP is also reported for its hydrocarbon degrading properties at high temperatures (Singleton et al., 2009). Microbes belonging to phylum Proteobacteria found abundant at this site are known to be highly abundant in petroleum-contaminated terrestrial and aquatic environments (Kimes et al., 2013; Bao et al., 2016). A 16S rRNA-based study from a hydrothermal vent (270– 325◦C) with similar geochemistry at the Gulf of California, where the organic matter in the sediments are converted to aliphatic and aromatic hydrocarbons under high temperature, suggested a resident sulfate-reducing bacterial population, such as Desulfobacter, Desulforhabdus, Thermodesulforhabdus, etc. (Dhillon et al., 2003).

The known high percentage of methane gases at Anhoni corroborates well with the presence of methanotrophic species such as M. capsulatus (Sarolkar, 2015). The higher abundance of this pathway, particularly in CAP, shows the importance of this site for the isolation of novel methanotrophic species. M. capsulatus is known to convert formaldehyde through Ribulose-P pathway to form Acetyl-CoA for subsequent energy generation pathways which were evident in our data (**Supplementary Figure 8**; Ward et al., 2004). Our data also suggests the conversion of formaldehyde to formate which then enters the Acetyl-CoA pathway. However, this pathway is not reported to be carried out by M. capsulatus (Ward et al., 2004), indicating the presence of other methanotrophs in CAP. Meiothermus ruber, another abundant species in CAP is a thermophile observed to be associated with moderately high temperature (Tindall et al., 2010). Methane-producing and sulfur-utilizing thermophilic bacteria belonging to the genus, such as Thermococcus, Acinetobacter, Pseudomonas, and Methylobacterium have also been reported to be associated with high-temperature petroleum reservoirs (50–120◦C) in California (Orphan et al., 2000).

One of the interesting insights from this study is that all enzymes involved in a particular hydrocarbon degradation pathway were not found in a single microbial species; however, as a microbial community, the entire set of genes for the degradation of that hydrocarbon could be completed. Furthermore, the results also hint toward the existence of novel species in these hot springs which might harbor the complete degradation pathways and can be confirmed by their isolation and sequencing, therefore, providing leads for further studies. The species and enzymes identified in this study and future studies from this site could be used as promising bioremediation agents in oil spills and other hydrocarbon-contaminated regions (Kimes et al., 2013; Gaur et al., 2014).

Tattapani hot spring having an extremely high temperature (61.5–98◦C) displayed a diverse thermophilic population and their pathway analysis suggested the functional adaptations of these thermophiles to survive at high temperatures. Extreme environments such as hot springs with very high temperature are commonly found to be inhabited by chemolithoautotrophic thermophilic species majorly from archaebacteria, and from Aquifex and Thermotoga phyla in eubacteria (Stetter, 1996, 1999). The high abundance of P. aerophilum in this environment is an interesting finding as it is a unique hyperthermophilic archaeal species reported extensively for its nitrate reducing properties with the help of a molybdenum-associated nitrate reductase which is distinct from other mesophilic nitrate reductases (Afshar et al., 2001). P. aerophilum also found to be abundant in Manikaran hot spring in India (96◦C), is also known to grow anaerobically by dissimilatory nitrate reduction, which is an exclusive catabolic feature of this species (Afshar et al., 2001; Bhatia et al., 2015). Fervidobacterium thermophilum and Fervidobacterium pennivorans are among the other thermophilic eubacteria observed to be abundant in Tattapani, and are reported to produce thermostable cellulases and keratinases (fervidolysin) respectively, which have potential application in the conversion of biomass to biofuels at high temperature and other biotechnological processes (Kim et al., 2004; Wang et al., 2010). Overall, a high species diversity was observed at the sites with moderately high temperature (40–55◦C) i.e., from the Anhoni site, compared to the site at extremely high temperature i.e., TAT-1 (98◦C; **Supplementary Figure 5**). The presence of hydrocarbons at Anhoni site could be a contributing factor for the increase in species diversity at this site, as it provides a rich carbon source to the microorganisms thus, helping them to survive despite high temperature (Head et al., 2006). Large community diversity is generally observed in petroleum or hydrocarbon-contaminated environments and is also observed in the Anhoni samples (Kimes et al., 2013; Bao et al., 2016).

The presence of genes for replication and DNA repair pathways appear to be important in Tattapani microbial community to survive in adverse physical conditions such as, high temperature, which cause considerable damage to DNA. The enrichment of genes associated with DNA repair systems and homologous recombination have been well-reported in extreme environments (Xie et al., 2011; Jimenez et al., 2012). Some recent studies have shown that the evolutionary rate of microbial communities is governed by the environmental conditions (Gupta and Sharma, 2015), and microbes in extreme habitats evolve faster with extensive DNA repair system and high mutation rates to cope with the deleterious effects of environment on their genomes as compared to those in stable environments (Li et al., 2014). Despite having a higher abundance of pathways such as replication and nucleotide metabolism, lower species richness was observed in Tattapani which indicates that very few microbial species are able to adapt and survive in extreme conditions.

The present metagenomic exploration of the hot springs located in central India is perhaps the largest comprehensive metagenomic study of any geographical location in India including hot springs, carried out using both 16S rRNA amplicon and shotgun sequencing. It has provided novel insights into the taxonomic, functional, and metabolic diversity of the unique microbial community surviving in these hot springs. The presence of chemoorganotrophic thermophiles degrading complex hydrocarbons at the Anhoni hot springs and extensive survival mechanisms at the Tattapani hot springs are novel revelations of the study. The present analysis is however limited by the availability of known species of microbes in reference datasets and opens up new opportunities to decipher the genomes of yet unknown microorganisms in these hot springs.

#### AUTHOR CONTRIBUTIONS

VS and RS conceived the idea. RS and PW designed and performed the experiments and sequencing. AC and AG carried out the elemental analysis and interpreted its results. DD and PM performed the computational analysis. RS, DD, PM, AG, and VS analyzed the data and interpreted the results. RS, DD, and PM prepared all the figures and tables. RS, DD, PM, AG, and VS wrote and reviewed the manuscript.

#### ACKNOWLEDGMENTS

We thank MHRD, Govt of India, funded Centre for Research on Environment and Sustainable Technologies (CREST) at IISER Bhopal for providing financial support. However, the views expressed in this manuscript are that of the authors alone and no approval of the same, explicit or implicit, by MHRD should be assumed. We thank NGS and HPC Facility at IISER Bhopal for carrying out the sequencing runs and computational analysis, respectively. RS and DD acknowledge Department of Science and Technology-INSPIRE and University Grants Commission (UGC), Government of India, respectively, for providing research fellowships. We thank Ankita Roy and Vibhuti Shastri for reading the manuscript and providing valuable suggestions.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fmicb. 2016.02123/full#supplementary-material

Supplementary Figure 1 | Geographical location and photographs of sampling sites of Anhoni and Tattapani hot springs. The given map was retrieved from Google Maps (2016).

Supplementary Figure 2 | Graphical representation of elemental composition at the seven sample sites.

Supplementary Figure 3 | Alpha diversity in the Hot spring samples.

(A) Shannon index and (B) Observed Species were calculated by rarefying from 100 to 2.2 million sequences at a step size of 0.1 million. This analysis was carried out using amplicon reads.

#### REFERENCES

Afshar, S., Johnson, E., de Vries, S., and Schroder, I. (2001). Properties of a thermostable nitrate reductase from the hyperthermophilic archaeon Pyrobaculum aerophilum. J. Bacteriol. 183, 5491–5495. doi: 10.1128/JB.183.19.5491-5495.2001

Supplementary Figure 4 | PCA plot prepared using Hellinger distances based on eggNOG proportions.

Supplementary Figure 5 | Multiple comparisons of number of observed species with moderately high, high and extremely high-temperature locations. Multiple comparisons using Tukey's test were performed between samples grouped on the basis of temperature (moderately high, high, and extremely high). The mean ± *SD* are plotted with significant variations shown as ∗∗∗*p* ≤ 0.001.

Supplementary Figure 6 | Significantly discriminatory genera (p ≤ 0.05) in Anhoni and Tattapani sites.

Supplementary Figure 7 | eggNOG functional classes and their distributions in the two sites. The box plot shows the abundance of eggNOG functional categories in Anhoni and Tattapani.

Supplementary Figure 8 | A pathway map of methane metabolism showing the observed KOs (highlighted in red) in this dataset. Methane is oxidized to methanol via methane monooxygenase and converted to formaldehyde with the help of methanol dehydrogenase. The KOs which are commonly present in all three sites (CAP, CAN, and BAN) of the Anhoni region are shown.

Supplementary Figure 9 | A pathway map of toluene degradation showing the observed KOs (highlighted in red) in this dataset. Toluene is degraded to Benzoate and 3-methylcatechol which enters the benzoate degradation and xylene degradation pathways respectively for further downstream processes. The KOs which are commonly present in all three sites (CAP, CAN, and BAN) of the Anhoni region are shown.

Supplementary Figure 10 | A pathway map of xylene degradation showing the KOs (highlighted in red) observed in our dataset. Xylene is degraded and converted to Acetyl-CoA by the community microbes and enters the TCA cycle for energy generation. The KOs which are commonly present in all three sites (CAP, CAN, and BAN) of the Anhoni region are shown.

Supplementary Figure 11 | A pathway map of benzoate degradation showing the observed KOs (highlighted in red) in this dataset. Benzoate is degraded and converted to intermediates such as Pyruvate and Succinyl-CoA, finally entering into TCA cycle for energy generation. The KOs which are commonly present in all three sites (CAP, CAN, and BAN) of the Anhoni region are shown.

Supplementary Figure 12 | Bar charts showing the proportion of key enzymes of benzoate pathway observed in this dataset along with their taxonomic annotations.

Supplementary Figure 13 | Utilization of hydrocarbons by the microbial community for generation of energy. The hydrocarbons are utilized as a substrate by the chemoorganotrophic thermophiles to produce energy. The degradation pathways which were observed to be completely present in Anhoni samples are shown with some of the intermediates and their end-products.

Supplementary Table 1 | Elemental analysis of the hot spring samples.

Supplementary Table 2 | Number of 16S rRNA (V3 hypervariable region) amplicon reads obtained per sample.

Supplementary Table 3 | Sequencing and assembly statistics of the hot spring samples obtained after de novo assembly. The size of contigs ranged from 1817 to 189,789 bp indicating a taxonomic diversity across different samples.

Supplementary Table 4 | The alpha diversity metrics calculated for 16SrRNA and metagenomics reads are shown for both sites (mean ± SD).

Supplementary Table 5 | Discriminatory taxa of Anhoni and Tattapani identified using Random Forest Analysis.

Altschul, S. F., Gish, W., Miller, W., Myers, E. W., and Lipman, D. J. (1990). Basic local alignment search tool. J. Mol. Biol. 215, 403–410. doi: 10.1016/S0022-2836(05)80360-2

Asnicar, F., Weingart, G., Tickle, T. L., Huttenhower, C., and Segata, N. (2015). Compact graphical representation of phylogenetic data and metadata with GraPhlAn. PeerJ 3:e1029. doi: 10.7717/peerj.1029

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2017 Saxena, Dhakan, Mittal, Waiker, Chowdhury, Ghatak and Sharma. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Aerobic Lineage of the Oxidative Stress Response Protein Rubrerythrin Emerged in an Ancient Microaerobic, (Hyper)Thermophilic Environment

#### Juan P. Cardenas1,2† , Raquel Quatrini3 and David S. Holmes1,2 \*

#### Edited by:

Kian Mau Goh, Universiti Teknologi Malaysia, Malaysia

#### Reviewed by:

William P. Inskeep, Montana State University, USA Amy Michele Grunden, North Carolina State University, USA Xi-Ying Zhang, Shandong University, China

> \*Correspondence: David S. Holmes<EMAIL_ADDRESS>
†Present address: Juan P. Cardenas, uBiome, Inc., San Francisco, CA, USA

#### Specialty section:

This article was submitted to Extreme Microbiology, a section of the journal Frontiers in Microbiology

Received: 03 August 2016 Accepted: 31 October 2016 Published: 18 November 2016

#### Citation:

Cardenas JP, Quatrini R and Holmes DS (2016) Aerobic Lineage of the Oxidative Stress Response Protein Rubrerythrin Emerged in an Ancient Microaerobic, (Hyper)Thermophilic Environment. Front. Microbiol. 7:1822. doi: 10.3389/fmicb.2016.01822 1 Center for Bioinformatics and Genome Biology, Fundacion Ciencia & Vida, Santiago, Chile, 2 Facultad de Ciencias Biologicas, Universidad Andres Bello, Santiago, Chile, 3 Laboratory of Microbial Ecophysiology, Fundación Ciencia & Vida, Santiago, Chile

Rubrerythrins (RBRs) are non-heme di-iron proteins belonging to the ferritin-like superfamily. They are involved in oxidative stress defense as peroxide scavengers in a wide range of organisms. The vast majority of RBRs, including classical forms of this protein, contain a C-terminal rubredoxin-like domain involved in electron transport that is used during catalysis in anaerobic conditions. Rubredoxin is an ancient and large protein family of short length (<100 residues) that contains a Fe-S center involved in electron transfer. However, functional forms of the enzyme lacking the rubredoxin-like domain have been reported (e.g., sulerythrin and ferriperoxin). In this study, phylogenomic evidence is presented that suggests that a complete lineage of rubrerythrins, lacking the rubredoxin-like domain, arose in an ancient microaerobic and (hyper)thermophilic environments in the ancestors of the Archaea Thermoproteales and Sulfolobales. This lineage (termed the "aerobic-type" lineage) subsequently evolved to become adapted to environments with progressively lower temperatures and higher oxygen concentrations via the acquisition of two co-localized genes, termed DUF3501 and RFO, encoding a conserved protein of unknown function and a predicted Fe-S oxidoreductase, respectively. Proposed Horizontal Gene Transfer events from these archaeal ancestors to Bacteria expanded the opportunities for further evolution of this RBR including adaption to lower temperatures. The second lineage (termed the cyanobacterial lineage) is proposed to have evolved in cyanobacterial ancestors, maybe in direct response to the production of oxygen via oxygenic photosynthesis during the Great Oxygen Event (GOE). It is hypothesized that both lineages of RBR emerged in a largely anaerobic world with "whiffs" of oxygen and that their subsequent independent evolutionary trajectories allowed microorganisms to transition from this anaerobic world to an aerobic one.

Keywords: rubrerythrin, evolution, phylogeny, comparative genomics, microaerophilic, hyperthermophiles, GOE, cyanobacteria

**84**

## INTRODUCTION

fmicb-07-01822 November 16, 2016 Time: 14:5 # 2

The ability to combat oxidative stress is a widespread feature found in most organisms including many obligate anaerobic and microaerophilic organisms. The Great Oxygenation Event (GOE), that has been hypothesized to occur approximately 2.3 billion years ago (Lyons et al., 2014), would likely have initiated an adaptation process to attenuate the threatening exposure to increased levels of reactive oxygen species (ROS). The development of efficient ROS scavenging mechanisms would have facilitated the co-evolution of redox proteins that could take advantage of the energetically advantageous use of oxygen as a terminal electron acceptor.

Reactive oxygen species are partially reduced oxygen compounds that are produced as byproducts of oxygen reduction and, in many conditions, their exacerbated production can lead to severe' stress and even cellular death (Imlay, 2003). Several mechanisms have evolved to mitigate the oxidative stress caused by ROS including direct mechanisms and indirect scavenging of ROS (Farr and Kogoma, 1991; Maaty et al., 2009; Zuber, 2009). The most stable ROS is hydrogen peroxide (H2O2) and a variety of enzymes have evolved to remove it, such as catalases, peroxiredoxins and other peroxidases (Mishra and Imlay, 2012). Among the H2O2–scavenging enzymes, one variant with an important role in stress survival in several microorganisms, is rubrerythrin.

Rubrerythrin (RBR) is a member of the Ferritin-like superfamily (FLSF), consisting of proteins with a variety of different functions such as iron storage and/or iron detoxification (ferritins, Dps proteins, bacterioferritins), ubiquinone biosynthesis (COQ7 proteins), and radical scavenging (Mn-catalases, RBRs; Andrews, 2010). These proteins possess a four-helical bundle, configured in two anti-parallel helix pairs, forming a di-iron center, mediated by the coordination of (at least) six highly conserved residues among the fourhelix structure (Andrews, 2010). In RBR, the di-iron center is responsible for the reduction of H2O2 and organic hydroperoxide (Dillard et al., 2011). The role of RBR in the reduction of H2O2 has been experimentally verified in several different organisms including aerobes (Sato et al., 2012), cyanobacteria (Zhao et al., 2007), and also obligate anaerobes (Weinberg et al., 2004). The classic RBR, found in the obligate anaerobe Desulfovibrio vulgaris contains, in addition to the FLSF domain, a C-terminal domain related to the rubredoxin family, putatively involved in electron transfer during catalysis (Andrews, 2010). This feature is conserved in most representatives of this family (Iyer et al., 2005), including a version of RBR with its rubredoxin domain in the opposite orientation with respect to the classic version called "reverse RBR" (Andrews, 2010). This suggests that the short rubredoxin-like domain has an important role in RBR function.

Despite the strong conservation of the rubredoxin-like domain in RBR proteins, there are also functional versions of this enzyme without the rubredoxin-like domain such as sulerythrin and ferriperoxin from Sulfolobus tokodaii (Wakagi, 2003) and Hydrogenobacter thermophilus (Sato et al., 2012). This onedomain RBR, composed only of the FLSF domain, was found to be widely distributed in microbial organisms, raising questions about how it evolved and whether the lack of the C-terminal rubredoxin-family domain had functional implications for the catalysis of H2O2. However, in obligate anaerobes, a functional and physical association between RBR and other proteins was demonstrated, suggesting that some other proteins might supply the missing catalytic function to the one-domain RBR (Weinberg et al., 2004; Le Fourn et al., 2011).

In this report, we have undertaken a phylogenomic analysis of RBRs using a number of techniques including conventional phylogenetic approaches, sequence similarity networks and comparisons of genomic neighborhoods. The objective was to derive a plausible trajectory of the evolution of RBRs and, if possible, to link this trajectory with postulated changes of temperature and atmospheric oxygen concentrations during the early stages of the earth over 3 billion years ago.

### MATERIALS AND METHODS

#### Compilation of RBR Sequences

Rubrerythrin sequences were obtained from NCBI nonredundant (NR) database using a two-step filter: first, the NR database was analyzed using HMMsearch (HMMer version 3.0) against PF02915 (PFAM domain for Rubrerythrin, E-value < 10−6 ) followed by a RPS-BLAST search against COG, recovering all proteins with significant similarity to COG1592 (E-value < 10−10); all sequences with lengths less than 100 residues and/or with COG coverage values less than 70% were discarded. **Supplementary Table S1** lists the sequences found using this strategy.

### Sequence Similarity Network

The sequence similarity network elaboration was as described by Atkinson et al. (2009). A total of 4527 RBR sequences found by the aforementioned method were clustered using CD-HIT software (Fu et al., 2012), resulting in 2631 representatives comprising groups defined by 90% identity. These filtered sequences were analyzed in a BLASTp-all-versus-all round (nofilter, default parameters) with a threshold E-value of 10−35. The pairwise bit scores were used as measure of distance for the network, visualized in Cytoscape 3.0.2 using Organic layout.

#### Phylogenetic Analyses

Protein sequences were aligned using MAFFT (Katoh and Standley, 2014) with the L-INS strategy. The phylogenetic analyses were performed using either maximum likelihood (ML) or Bayesian inferences (BI). For ML, the trees were elaborated in PhyML, version 3.1 (Guindon et al., 2009) with different sets of parameters for each case. For BI, the trees were elaborated in MrBayes version 3.2.6 (Huelsenbeck and Ronquist, 2001), with different parameters for each case (see "Results"). The selected substitution models for both ML and BI analyses were selected using ModelGenerator (Keane et al., 2006).

For the "aerobic-type" RBR tree, the alignment containing 127 ungapped and unambiguous sites for 335 taxa was analyzed by MrBayes with five substitution rate categories and gammadistributed rate variation, using LG as the prior model. Bayesian

analysis was run for five million generations (in two independent runs, using four chains and a heating parameter of 0.1); trees were saved every 200 generations and posterior probabilities calculated after discarding the first 33% of trees.

For the DUF3501 and the rubrerythrin-associated Fe-S oxidoreductase (RFO) protein families (see below), the trees were constructed using PhyML. For DUF3501, the alignment containing 118 ungapped and unambiguous sites for 258 taxa was used; the ML-based tree was performed using the LG model with four substitution rate categories and gamma-distributed rate variation, with a proportion of invariable sites. The tree topology search was performed by the combination of NNIs and SPRs strategies and the approximate likelihood ratio was computed by the SH-like branch support test (Anisimova et al., 2011). The RFO tree was computed using similar parameters, from an alignment containing 289 ungapped and unambiguous sites for 206 taxa. Additionally, a tree of the "cyanobacterial group" RBR orthologs was computed from an alignment containing 289 ungapped and unambiguous sites for 206 taxa, using the same configuration used for the two latter phylogenetic analyses.

#### Other Analyses

Genomic contexts for the genes encoding rubrerythrin were obtained from IMG-JGI (Markowitz et al., 2012), using a set of 4461 complete prokaryotic genomes. Protein domain analyses were carried out using InterProScan (Quevillon et al., 2005). Co-localized genes with RBR were searched in the complete prokaryotic genome collection using BLAST (Altschul et al., 1997).

#### RESULTS AND DISCUSSION

#### Identification of Different Rubrerythrin Groups

In order to study the evolution of RBR forms, the first step was the preparation of a trustworthy set of these proteins. The search of putative RBR homologs was made using a dual criterion selection that combines the sensibility of Hidden Markov Model (HMM) based search (PF02915) with a confirmative application of RPS-BLAST against the COG1592 family (corresponding to the RBR family). This double filter allowed a discrimination to be made between RBRs and other members of the FLSF, including members of the Mg+2 -protoporphyrin IX monomethyl ester cyclase subfamily, which were also included in the HMM profile for PF02915.

A sequence similarity network approach was implemented in order to identify subfamilies inside the RBR family. Similarity networks are acceptable approximations for the study of huge protein families, although they are not a replacement for phylogenetic studies (Atkinson et al., 2009). The use of a given E-value (or score) threshold in an "all-vs-all" BLAST permits the recovery of clusters that can be visualized in a network fashion, where each node is an homolog or family member, and each edge is the measure of the pairwise score obtained by BLAST. The distances between different node groups are inversely proportional to the value of their pairwise scores: higher BLAST pairwise scores between multiple sequences promote their clustering, and vice versa.

**Figure 1** shows a sequence similarity network for the RBR orthologs retrieved from the NR database following the strategy outlined above. From this network (comprising 2631 representative genes), a set of node clusters or groups can be inferred. Four large groups (1–4) connected in a major network were recovered using an E-value threshold (E < 10−45). In addition, there are nine minor isolated groups (5–13), and several smaller groups with disconnected nodes, indicating that, at the BLAST E-value used, those sequences have no similarity with other members of the RBR protein set. Group 5 is included with groups 1–4 in subsequent analyses. Groups 1–5 comprise ∼93% of the sequences.

Group 1 (**Figure 1**) corresponds to the principal group of canonical rubrerythrins, and is comprised almost completely (98.7 % of 1336 sequences) by members with the rubredoxin-like domain with a relatively short spacing between the two C-x-x-C motifs (10–14 residues, we term "short-spaced"). Taxonomic information (**Supplementary Figure S1**) indicates that this group is dominated by members of Clostridia (44.5%) and δ/ε-Proteobacteria classes (11.0%) and the Bacteroidetes phylum (10.2%). Comparison with PDB database showed that the classical RBR protein from D. vulgaris (Swissprot accession P24931, PDB entry 1B71) is a member of this group. Therefore, we propose that group 1 is termed the "Classical RBR" group. The group is populated predominately by obligate anaerobic microorganisms.

Group 2 (**Figure 1**) consists of two sub-clusters: the first (85.3% of the total sequences) exhibits the rubredoxin domain as do group 1 rubredoxins, but unlike group 1 this rubredoxin domain is located in the N terminal of the protein as previously observed in Clostridium acetobutylicum (Riebe et al., 2009). The other sub-cluster (14.7% of the total sequences) does not have a rubredoxin domain. Group 2 was dominated obligate anaerobes such as members of the Clostridia class (64.5 %) and Bacteroidetes phylum (10.5%; see **Supplementary Figure S1**). We propose that group 2 be termed the "Reverse-type" RBR group.

Group 3 (**Figure 1**) is mainly composed of RBR with the C-terminal rubredoxin-like domain and, according to taxonomic information, is dominated by members of δ/ε-Proteobacteria (25.97%) and Clostridia (20.7%) classes and the Euryarchaeota phylum (17.31%). Comparison with PDB showed that a previously crystallized RBR from Pyrococcus furiosus belongs to this group (Uniprot accession Q9UWP7, PDB entry 1NNQ); this RBR contains canonical hydrogen peroxide reductase activity (Dillard et al., 2011). We propose that this group is termed the "Classical group B." Like group 1, both the "reverse" (group 2) and the "classical-B" types of RBR (group 3) are enriched in taxonomic groups dominated by obligate anaerobes.

Group 4 (**Figure 1**) is composed exclusively of genes lacking the rubredoxin-like domain. This group is dominated by members of α/β/γ-Proteobacteria (61.1%), Actinobacteria (9.0%), and Crenarchaeota (including members of Sulfolobales and Thermoproteales orders, 8.3%), as well as bacterial members from Aquificae (including the ferriperoxin from H. thermophilus), Nitrospirae and some few members of

Firmicutes, Deltaproteobacteria, and archaeal members of Thermoplasmatales (from the Euryarchaeota phylum). Unlike the three aforementioned groups, members of group 4 are dispersed in taxonomic groups dominated by (facultative) aerobic organisms. Sulerythrin from S. tokodaii (Uniprot accession F9VPE5; PDB entry 1J30) and a rubrerythrin from Burkholderia pseudomallei (Uniprot accession Q3JK2; PDB entry 4DI0) have three dimensional structures in the PDB. The presence of RBRs without the rubredoxin-like domain is distributed in different sequence clusters of this network, suggesting that these proteins arose independently more than once in evolution. The detection of a well-defined group composed only of members without the aforementioned domain, and additionally, that those homologs are mostly associated with aerobic organisms, strongly suggest that this group is a complete and distinctive evolutionary lineage with common properties and we term it the "aerobic-type" group of RBR.

Group 5 (**Figure 1**) contains two well-compartmentalized sub-clusters; one of them (comprising 41.8% of the sequences) lacks the rubredoxin-like domain. Members of the other subcluster (comprising the 58.2% of the sequences) is distinguished by having a C-terminal rubredoxin-like domain, with a longer spacing between the cysteines motifs (27–32 residues) compared with the C-terminal domain of the classical RBR (10–14 residues). Taxonomically, the group 5 is dominated by members of Cyanobacteria (59.5%) and α/β/γ-Proteobacteria (36.7%). The "long-spaced" rubredoxin-like domain is present in almost all cyanobacterial RBRs and has a different evolutionary origin from the other RBRs. Due to the prevalence of cyanobacterial RBR genes in group 5, we term this the "Cyanobacterial-type" RBR. Given their evolutionary properties, the phylogenetic analysis of this group, and the origin of the "long" rubredoxin domain will be covered below.

### Phylogenetic and Genomic Analyses of the "Aerobic-Type" RBR

In order to analyze the evolutionary history of this proposed lineage, a phylogenetic tree was constructed using BI (**Figure 2A** in summary and **Supplementary Figure S2** in detail). The "aerobic-type" RBR phylogenetic tree (using an appropriate outgroup) showed that the closest clade to the root of the group comprises a set of orthologs from archaeal members (**Figure 2B** in summary and **Supplementary Figure S2** in detail) from the Vulcanisaeta, Caldivirga, and Thermoproteus genera of the Thermoproteales order (Itoh et al., 1999; Mavromatis et al., 2010; Gumerov et al., 2011; Mardanov et al., 2011). All of the extant members of this group are hyperthermophiles and are slightly acidophilic (pH 3.0–6.5). They grow only in anaerobic or microaerobic conditions. For example, a well-supported clade contains orthologs from members of Actinobacteria (species from Acidithrix, Ferrimicrobium, and other related groups), Nitrospirae (species from Leptospirillum and Nitrospira), and Clostridia (Sulfobacillus), associated with extremely acidic (pH < 3) environments (Cardenas et al., 2016). Another wellsupported clade contains only orthologs from members of the Sulfolobales order (from Sulfolobus, Metallosphaera, and Acidianus genera). It is interesting to note that this clade is paraphyletic with respect to other clades containing other orthologs from Archaea (for example, the clades containing orthologs from Thermoproteales or Thermoplasmatales). This

FIGURE 2 | Phylogenetic trees constructed in this study. The circular tree for representative members of the "aerobic-type" RBR (A) and (B) with a phylogram representation of the closest clade to the root of the tree (B). The figure also includes unrooted versions of the trees made from the protein sequences of DUF3501 (C) and RFO (D). For the tree of the "aerobic-type" RBR group, the sequence of a RBR from Cyanothece sp. PCC 8802 (a cyanobacterial RBR) was used as an outgroup taxon, whereas the tree for RFO used the sequence of GlpC from Escherichia coli as outgroup. For each case, the line length represents the respective value for phylogenetic distance. Organism names and more detailed trees are displayed in Supplementary Figures S2–S5.

suggests that the evolutionary history of this lineage of RBR involved multiple different horizontal transferences within archaeal organisms. Another noticeable observation inferable from this tree is the expansion of this RBR lineage inside in the Alpha-, Beta- and Gammaproteobacteria. Most importantly, the phylogenetic tree of RBR suggests that the most ancient group (the closest to the root) is composed of members from the Thermoproteales; this is consistent with the idea that this form of RBR evolved in a microaerobic and thermophilic environment.

### Gene Context Analyses of Genes/Domains Predicted to Be Related to "Aerobic-Type" RBR Function

We have undertaken an analysis of the gene contexts of the "aerobic-type" of RBR in order to shed light on its evolution. A comparative visualization of RBR gene neighborhoods among complete genomes was carried out using the IMG-JGI database (Markowitz et al., 2012). The gene neighborhood comparison (**Figure 3**) showed that the majority of the RBR genes from

the "aerobic-type" were co-localized in a conserved gene cluster with a gene encoding a protein of unknown function (DUF3501, pfam12007), and in several other cases, also co-localized with a gene predicted to encode a member of the COG0247 family (Fe-S oxidoreductase). The role of the DUF3501 protein is unknown. However, the function of the COG0247 family protein has been suggested to be related to electron transfer inside protein complexes, since it contains a conserved CCG domain (pfam02754) and a cysteine rich domain. The latter is present in several oxidoreductase complexes related to energy metabolism under anaerobic conditions such as subunit GlpC in the anaerobic glycerol-3-phosphate dehydrogenase (Cole et al., 1988) and subunits of the HdrD/E from the heterodisulfide reductase complex found in methanogens (Kunkel et al., 1997). We propose that this COG0247 family protein member associated with RBR be termed RFO (Rubrerythrin-associated Fe-S Oxidoreductase). In some instances, DUF3501- and RFO are fused into one gene (e.g., Nitrosomonas spp., **Figure 3**). In other organisms (Thermocrinis spp. and H. thermophilus), the RFO-encoding gene is separated from the other two genes. It is interesting to note that other members of the deeply rooted Aquificae phylum (such as Hydrogenobaculum) have the more common three-gene arrangement. This difference in gene cluster structure among members of the same phylum may be a reflection of the antiquity of the events that resulted in these different gene arrangements.

Additional comparisons suggest that the genes encoding DUF3501 and RFO are not only strongly co-localized with the "aerobic-type" RBR, but also are co-occurrent. An examination of a set of 4461 completed prokaryotic genomes from the IMG-JGI database showed that DUF3501 and RFO genes are not associated with RBR in only one organism (Burkholderia cepacia AMMD) and only in one other genome (Azoarcus sp. CIB) a DUF3501 encoding gene was detected without the other two, RFO and RBR (**Supplementary Table S2**). The exclusive co-occurrence of RBR – DUF3501 was found in 25 organisms, and the triple cooccurrence of RBR, DUF3501, and RFO was detected in 185 genomes. Interestingly, the co-occurrence of RBR – DUF3501 was detected in the members of Thermoproteales, suggesting that this association could be as ancient as the last common ancestor of this group, supporting the contention that the origin of "aerobic-type" RBR lineage, occurred in a high-temperature and microaerobic environment.

### Phylogenetic Analysis of DUF3501 and RFO Families

In order to investigate the possible co-evolution of the "aerobic-type" RBR with the DUF3501 and RFO genes, phylogenetic trees for DUF3501 (**Figure 2C** in summary and **Supplementary Figure S3** in detail) and RFO (**Figure 2D** in summary and **Supplementary Figure S4** in detail) were constructed. These trees are consistent with the late expansion of the respective protein families into the Alpha-, Beta- and Gammaproteobacteria, and with the observed expansion of the co-localizing "aerobic-type" RBR. The case of DUF3501 also showed the non-monophyletic association between the archaeal orthologs from the three different orders (Thermoproteales, Sulfolobales, and Thermoplasmatales), seen previously in the RBR tree. In the case of the RFO tree, BLAST searches of orthologs for this phylogenetic tree showed that the most closely related proteins are GlpCs that encode the subunit C of the anaerobic glycerol-3-phosphate dehydrogenase (involved in glycerol degradation under anaerobic conditions). The

application of GlpC as an outgroup showed that the most ancient group corresponds to Sulfolobales, an archaeal order composed of aerobic thermoacidophiles (Huber and Stetter, 2015), suggesting that both GlpC as well as the "aerobictype" RBR, arose in high-temperature aerobic environments. Both phylogenetic trees for DUF3501 and RFO supporting the contention of the co-evolution of DUF3501 and RFO with the "aerobic-type" RBR.

### The "Cyanobacterial Group": A Parallel Adaptation of Rubrerythrin to an Aerobic Environment

Information obtained from phylogenomic analyses of the "aerobic-type" RBR and its co-occurrent genes suggests an evolutionary trajectory driven by the adaptation to aerobic environments, initiated from an (hyper)thermophilic ancestor. However, another group in the RBR sequence similarity network was also detected (the "cyanobacterial group," group 5 from **Figure 1**) that was not connected to the other four groups and was composed mostly of members from Cyanobacteria (phototrophic oxygen-producers) and α/β/γ-Proteobacteria (facultative anaerobes), raising the question about how this protein group evolved. A phylogenetic tree for members of the "cyanobacterial RBR" group was made using ML (**Supplementary Figure S5**). This analysis showed two welldefined clades: a clade containing the cyanobacterial protein and a clade containing the proteobacterial orthologs. The members of the cyanobacterial clade contained a fusion of the FLSF of RBR with a "long-spaced" rubredoxin domain (data not shown). Based on the well-defined phylogenetic bifurcation and the coherent pattern of inheritance, the most parsimonious explanation is that the original RBR ancestor for this group was a one-domain RBR (as the "aerobic-type" RBR is), but in Cyanobacteria, the protein was fused with the "long-spaced" rubredoxin domain. We hypothesize that this adaptation had a relationship with the exposure to oxygen associated with the development of oxygenic photosynthesis in Cyanobacterial ancestor during GOE.

### Evolution of the "Aerobic" and "Cyanobacterial" Rubrerythrins in the Early Earth: A Model

We propose a model (**Figure 4**) for the evolution of the "aerobic" and "cyanobacterial" rubrerythrins that is consistent with evidence provided by network information, phylogenetic trees and gene/domain contexts of the different forms of rubrerythrins. It is suggested that the "classical" form of rubrerythrin arose by gene fusion of a FLSF domain and a rubredoxin-like domain in a hot and anaerobic environment. This form is predicted to have existed in LUCA (Weiss et al., 2016). Extant obligate thermophilic and mesophilic anaerobes retain this form today [e.g., Thermotoga maritima (Lakhal et al., 2011) and D. vulgaris (Lumppio et al., 1997), respectively].

Two separate lines of evolution then occurred. In the first, that is postulated to have occurred in the ancestors of the thermophilc Archaeal group Thermoproteales, we hypothesize that the loss of the rubredoxin-like domain was compensated for by the acquisition of DUF3051. This association is proposed to have provided protection in microaerobic environments that might have arisen before GOE due to exposure to low levels of oxygen. These low levels and perhaps transitory increments of oxygen could have originated from biotic or abiotic reactions, such as enzymatic ROS-detoxifying reactions or water UV photolysis (Martin and Sousa, 2016). It is probable that during those "oxygen whiffs," many proteins would be the object of selection pressure to survive in oxygen that could force them to lose or gain new features such as protein domains (Nasir et al., 2014) or physiological activities. In addition, the ancient hot environment of the Archaean eon could helped to accelerate the evolutionary rate of change of these proteins, as postulated to be occurring in current extreme environments (Li et al., 2014). Subsequently, an additional gene, RFO, was added that provided full function in aerobic environments perhaps resulting from GOE. This step is proposed to have occurred in the ancestors of the thermophilc Archaea Sulfolobales. This triple gene form was passed by HGT to the Bacteria and underwent further adaption to work in mesothermic conditions (e.g., <40◦C). It is not known why the additional DUF3051 and RFO genes are not found fused as domains to the FLSF domain but rather they remain as conserved co-localizing but separate genes. Perhaps this architecture provides more opportunities for differential expression of the genes.

The second line, involved the replacement of the rubredoxinlike domain with a separately evolved cyanobacterial rubredoxinlike domain. It is hypothesized that this promoted the function of RBR in high levels of oxygen resulting from oxygenic photosynthesis that is thought to have evolved in the cyanobacterial lineage about the time of GOE.

Oxidative stress is an ancient and widespread phenomenon: several of the extant obligate anaerobes exhibit mechanisms to combat the presence of ROS (Jenney et al., 1999). One of those mechanisms is the use of RBR. This protein has been suggested to be an ancient protein (Andrews, 2010) predicted to be present in LUCA that was itself hypothesized to be an obligate anaerobe and thermophile (Martin and Sousa, 2016). The proposed antiquity of RBR makes it an interesting protein to be analyzed, establishing a case study of the effect of the arise of oxygen in the adaption and origin of new protein (sub)families. It is interesting to note that the vast majority of the known RBRs have the "short-spaced" rubredoxin-like domain, even in a "reverse" form (**Figure 1**). This strongly suggest that, regardless that the very last common ancestor of all rubrerythrins lacked this rubredoxin domain, this fusion was implemented very early in evolution, and maybe even in different independent events (e.g., generating the "classic" and "reverse" RBR separately), and those successful, two-domain RBRs have been inherited by extant aerobic microorganisms.

There is phylogenomic evidence of ancient, massive gene birth events. Large-scale phylogenetic reconstruction of more than 3,000 gene families predicted a massive event of birth and loss of gene families occurring approximately 3.3–2.9 × 109 years ago (David and Alm, 2011), coinciding with or perhaps preceding, with the estimated time of GOE (Bekker, 2014). Another study (using a similar strategy) suggests that aerobic metabolism appeared about 2.9 × 109 years ago (Wang et al.,

Rubrerythrin-associated Fe-S Oxidoreductase.

2011). In that context, it is probable that this "aerobic-type" lineage of RBR, as well as their co-occurrent genes, could arise before GOE or just as GOE was starting. Regarding the case of the "cyanobacterial group" of RBR, the evidence provided by the sequence similarity network (group 5, **Figure 1**) suggests that this group of RBRs has a different evolutionary pathway from the "aerobic-form."

It is interesting to note that the 3-genes form of RBR was most likely to have evolved in the ancestors of two archaeal orders: Thermoproteales and Sulfolobales. Both those two taxonomic orders belong to the Crenarchaeota phylum, and have a proposed ancient point of divergence (Gribaldo and Brochier-Armanet, 2006). Whereas DUF3501 may have arisen in temperatures over 85◦ C (temperature associated with the lifestyle of extant members of Thermoproteales), the RFO gene could have arisen in temperatures between 60–80◦C (associated with the lifestyle of extant Sulfolobales). We hypothesize that the acquisition of DUF3501 was part of an adaptation for function in combined microaerobic and hyperthermic conditions, whereas RFO evolving in direct response to increasing oxygen presence. One observation, consistent with this supposition, is that the RFO gene is not found in microaerobic organisms, such as members of Thermoproteales and some members of Deltaproteobacteria.

We suggest that, since RFO belongs to a group of proteins potentially involved in electron transfer, its linkage with the "aerobic-type" RBR can compensate for the loss of the rubredoxin-like domain, a component of the anaerobic RBRs that has a proposed role in electron transfer during catalysis (Kurtz, 2006). However, this also raises another question: how is the catalysis of electron transfer in the RBR form found in Thermoproteales, where DUF3501 but not RFO is present? This question remains to be experimentally addressed.

The proposed origin and evolution of the "aerobic-type" RBR is consistent with the hypothesis that early life evolved in a hot, anaerobic environment (Martin and Sousa, 2016). It also supports the idea that oxidative stress mechanisms could have been present before GOE. It also provides an exquisite opportunity to study the evolutionary trajectory of proteins as they adapt to increasing oxygen levels and to econiches with lower temperatures.

#### AUTHOR CONTRIBUTIONS

fmicb-07-01822 November 16, 2016 Time: 14:5 # 9

JC, RQ, and DH conceived the project. JC designed and carried out the experiments. All authors analyzed the data. JC wrote the first draft of the manuscript. All authors contributed to other drafts of the manuscript. All authors read and approved the final manuscript.

#### ACKNOWLEDGMENT

This work was supported by Conicyt Basal CCTE PFB16 and Fondecyt 1130683 and 1140048.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fmicb. 2016.01822/full#supplementary-material

#### REFERENCES

FIGURE S1 | Phylogenetic distribution of rubrerythrin (RBR) types belonging to groups 1–5 (see "Text" of article). Sequences were retrieved

from the network shown in Figure 1 and were sorted by taxonomic origin (Phylum or Class) and counted.

FIGURE S2 | Phylogenetic distribution of the "aerobic" group of RBRs. The tree was elaborated as specified in "Material and Methods." Each taxon is tagged with an NCBI Accession number shown in parentheses.

FIGURE S3 | Phylogenetic distribution of members of the DUF3501 protein family. The tree was elaborated as specified in "Material and Methods." The branch support value is symbolized by the color of each node circle (color gradient legend shown). This tree is unrooted. Each taxon is tagged with an NCBI Accession number shown in parentheses.

FIGURE S4 | Phylogenetic distribution of members of the RFO protein family. The tree was elaborated as specified in "Material and Methods." The branch support value is symbolized by the color of each node circle (color gradient legend shown). This tree was rooted using the sequence of GlpC from Escherichia coli as an outgroup. Each taxon is tagged with an NCBI Accession number shown in parentheses.

FIGURE S5 | Phylogenetic distribution of members of the "cyanobacterial

group" of RBRs. The tree was elaborated as specified in "Material and Methods." The branch support value is symbolized by the color of each node circle (color gradient legend shown). In the base of each great clade, the domain fusion/separation is specified. This tree was rooted using the sequence of an "aerobic-type" RBR derived from Vulcanisaeta distributa as an outgroup. Each taxon is tagged with an NCBI Accession number shown in parentheses. Abbreviation: RB, rubredoxin.

#### TABLE S1 | List of rubrerythrin (RBR) sequences found in NR database.

The table contains the NCBI accessions, protein length, and description for each sequence.

TABLE S2 | List of candidate aerobic rubrerythrins, DUF3501 and RFO in prokaryotic finished genomes from IMG-JGI database. This table includes a Venn diagram displaying the co-occurrence of these genes.

from Pyrococcus furiosus. J. Biol. Inorg. Chem. 16, 949–959. doi: 10.1007/s00775- 011-0795-6

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2016 Cardenas, Quatrini and Holmes. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Gene Turnover Contributes to the Evolutionary Adaptation of Acidithiobacillus caldus: Insights from Comparative Genomics

Xian Zhang1,2, Xueduan Liu1,2, Qiang He3 , Weiling Dong1 , Xiaoxia Zhang4 , Fenliang Fan5 , Deliang Peng6 , Wenkun Huang6 and Huaqun Yin1,2 \*

1 School of Minerals Processing and Bioengineering, Central South University, Changsha, China, 2 Key Laboratory of Biometallurgy of Ministry of Education, Central South University, Changsha, China, 3 Department of Civil and Environmental Engineering, the University of Tennessee, Knoxville, TN, USA, 4 Institute of Agricultural Resources and Regional Planning, Chinese Academy of Agricultural Sciences, Beijing, China, 5 Key Laboratory of Plant Nutrition and Fertilizer, Chinese Academy of Agricultural Sciences, Beijing, China, 6 State Key Laboratory for Biology of Plant Diseases and Insect Pests, Institute of Plant Protection, Chinese Academy of Agricultural Sciences, Beijing, China

#### Edited by:

Edgardo Donati, National University of La Plata, Argentina

#### Reviewed by:

Om V. Singh, University of Pittsburgh, USA Jorge Hernan Valdes, Center for Genomics and Bioinformatics/Universidad Mayor, Chile

> \*Correspondence: Huaqun Yin<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Extreme Microbiology, a section of the journal Frontiers in Microbiology

Received: 27 June 2016 Accepted: 22 November 2016 Published: 06 December 2016

#### Citation:

Zhang X, Liu X, He Q, Dong W, Zhang X, Fan F, Peng D, Huang W and Yin H (2016) Gene Turnover Contributes to the Evolutionary Adaptation of Acidithiobacillus caldus: Insights from Comparative Genomics. Front. Microbiol. 7:1960. doi: 10.3389/fmicb.2016.01960 Acidithiobacillus caldus is an extremely acidophilic sulfur-oxidizer with specialized characteristics, such as tolerance to low pH and heavy metal resistance. To gain novel insights into its genetic complexity, we chosen six A. caldus strains for comparative survey. All strains analyzed in this study differ in geographic origins as well as in ecological preferences. Based on phylogenomic analysis, we clustered the six A. caldus strains isolated from various ecological niches into two groups: group 1 strains with smaller genomes and group 2 strains with larger genomes. We found no obvious intraspecific divergence with respect to predicted genes that are related to central metabolism and stress management strategies between these two groups. Although numerous highly homogeneous genes were observed, high genetic diversity was also detected. Preliminary inspection provided a first glimpse of the potential correlation between intraspecific diversity at the genome level and environmental variation, especially geochemical conditions. Evolutionary genetic analyses further showed evidence that the difference in environmental conditions might be a crucial factor to drive the divergent evolution of A. caldus species. We identified a diverse pool of mobile genetic elements including insertion sequences and genomic islands, which suggests a high frequency of genetic exchange in these harsh habitats. Comprehensive analysis revealed that gene gains and losses were both dominant evolutionary forces that directed the genomic diversification of A. caldus species. For instance, horizontal gene transfer and gene duplication events in group 2 strains might contribute to an increase in microbial DNA content and novel functions. Moreover, genomes undergo extensive changes in group 1 strains such as removal of potential non-functional DNA, which results in the formation of compact and streamlined genomes. Taken together, the findings presented herein show highly frequent gene turnover of A. caldus species that inhabit extremely acidic environments, and shed new light on the contribution of gene turnover to the evolutionary adaptation of acidophiles.

Keywords: Acidithiobacillus caldus, comparative genomics, intraspecific diversity, gene turnover, evolutionary adaptation

## INTRODUCTION

fmicb-07-01960 December 2, 2016 Time: 15:36 # 2

Acidithiobacillus caldus (formerly Thiobacillus caldus), a moderately thermophilic, obligately chemolithoautotrophic, and extremely acidophilic sulfur-oxidizing bacterium (Hallberg and Lindström, 1994, 1996), is of interest for its potential role in industrial bioleaching (Rawlings, 1998; Dopson and Lindström, 1999). A. caldus exploits elemental sulfur and a wide range of reduced inorganic sulfur compounds at moderately high temperatures to support autotrophic growth (Mangold et al., 2011; Chen et al., 2012). It is the primary member of a consortium of sulfur oxidizers in different toxic-laden acidic environments, which are termed "extreme environments," including coal pile and spoil, gold-bearing reactor operation, as well as low-grade copper bioleaching heap (Valdes et al., 2009; You et al., 2011; Zhang et al., 2016c). Considering that A. caldus inhabits harsh environments for prolonged periods and accommodates both sudden stress changes and long-term stress conditions in various habitats, gene flow and genetic drift might frequently occur. As such, the flexible gene repertoire generated by gene exchange has imparted A. caldus with extensive genetic material for diversification of function and phenotype. Therefore, research focusing on the correlation between genomic changes and evolutionary adaptation is of great interest.

The accumulation of genomic changes underlying evolutionary adaptation has often been viewed as a complex process, and has been subject to many influences and complications (Barrick et al., 2009). As stated by Carretero-Paulet et al. (2015), homologous genes derived from newly formed subgenomes might undergo asymmetric fractionation via mutational events, which include nucleotide substitutions, gene gains and losses, and changes in genomic structure and organization (Librado et al., 2014). In terms of neutral mutation theory, mutations underlying gene and genome evolution, though not necessarily beneficial, should accumulate at a constant rate by drift (Kimura, 1984). Another view is that the substitution rates for beneficial and deleterious mutations depend on environmental selection, as well as population size and structure (Gillespie, 1991; Ohta, 1992). For many years, the crucial role of gene and genome duplications (namely, neofunctionalization and subfunctionalization) in governing organismal evolution has been acknowledged (Ohno, 1970; Force et al., 1999; Innan and Kondrashov, 2010; Kulmuni et al., 2013). Only in recent decades has great attention been paid to the molecular mechanisms of gene loss (deletion or pseudogenization) as a pervasive source of genetic change, which is believed to be another key evolutionary event that causes adaptive phenotypic diversity (Albalat and Cañestro, 2016). In recent years, a number of analytical methods for population genomics and molecular evolution have provided substantial evidence to determine the relative contribution of diverse evolutionary forces, which shape genome organization, architecture, and diversity in response to environmental perturbations (Librado et al., 2014). In eukaryotes, gene family evolution has often been modeled after a phylogenetic birth-anddeath (BD) process (Nei and Rooney, 2005). This BD model, though suitable to account for single-gene duplications, might not be appropriate for calculating gene turnover rates given that horizontal gene transfer (HGT) events occur in certain organisms (Librado et al., 2014). However, an alternative gain-and-death (GD) stochastic model in a maximum-likelihood statistical framework was applied to circumvent this limitation (Librado et al., 2012). Unlike the birth process, gains in the developed GD model can accommodate all kinds of gene acquisitions, irrespective of their original source, even including HGT (Librado et al., 2014). In this study, we are interested in whether the aforementioned theoretical and analytical approaches can be applied to explain the relationship between genetic change and adaptive evolution of A. caldus inhabiting extraordinarily extreme environments.

Members of A. caldus species are ubiquitous throughout many sulfur-rich acidic environments worldwide (**Table 1**), indicating their adaptation to various niches with high concentrations of toxic substrates, such as coal spoil, gold-bearing bioleaching reactor, and copper mine tailing. In recent years, revolutionary technologies and tools have allowed for the rapid characterization of microbial genome sequences (MacLean et al., 2009; Metzker, 2010). Accurate analyses of gene family evolution have been made possible owing to the increasing availability of closely related genomes (Hahn et al., 2007; Sánchez-Gracia et al., 2009; Vieira and Rozas, 2011). Furthermore, acquisition of numerous additional genomes has fuelled a new field termed comparative genomics (Jacobsen et al., 2011), which is useful for investigating microbial genome evolution and even mechanisms for speciation (González et al., 2014; Justice et al., 2014; Ullrich et al., 2016). Comparative surveys based on the available genomes of the two A. caldus strains ATCC 51756 and SM-1 have revealed that both strains harbor a relatively high proportion of unique gene complements (Acuña et al., 2013). These gene complements represent a diverse pool of mobile genetic elements, including insertion sequences (ISs), genomic islands (GIs), and integrative conjugative and mobilizable elements. Yet, limited information is available on the contribution of diverse evolutionary forces to the genomic diversification of A. caldus. Given this knowledge gap, we have isolated and sequenced four new A. caldus strains from different geographic origins (**Table 1**).

In this study, we estimated the phylogenetic relationships of A. caldus strains based on their genomic sequences (four newly sequenced genomes and two existing genomes from a public database), and performed an exhaustive study of the GD dynamics, with special focus on genetic exchange underlying evolutionary adaptation. These findings, to some extent, highlight the role of gene turnover in the evolutionary diversification of A. caldus and adaptation to specific lifestyles and environmental niches.

#### MATERIALS AND METHODS

### DNA Sequencing and Bioinformatics Analysis

Genome sequences for six strains were retrieved in this study, including A. caldus ATCC 51756, SM-1, DX, S1, ZBY, and ZJ.

∗Genome completeness was estimated using the CheckM. Strains SM-1 and ATCC 51756 with complete genome were excluded.

Of these bacteria, the type strain ATCC 51756 was isolated from a coal spoil in Kingsbury, UK (Marsh and Norris, 1983), strain SM-1 was from an industrial reactor used in bioleaching operation (Liu et al., 2007), and the other strains (DX, S1, ZBY, and ZJ) were obtained from the China Center for Type Culture Collection. More details for geographic origins of these four new strains were shown in **Table 1**. Genome sequences of strains ATCC 51756 and SM-1, including chromosomal and plasmid sequences, were downloaded from the GenBank database. For strains DX, S1, ZBY, and ZJ, chromosomal DNA was sequenced by an Illumina MiSeq sequencer (Illumina, Inc., USA), using the paired-end sequencing approach with an average DNA insert size of 300 bp and typical readlength of 150 bp. Subsequently, bioinformatics analysis of raw sequences was performed as described previously (Yin et al., 2014), primarily including quality control, genome assembly, computational prediction of coding sequences (CDS) and other genome features such as rRNA and tRNA, as well as functional assignments against public databases (NCBI-nr and COG). Genome completeness of each strain was also estimated using the program CheckM (Parks et al., 2015). Additionally, circular maps showing chromosome architecture were drawn using the Circos software (Krzywinski et al., 2009).

Intergenomic distance scores were calculated using the web service Genome-to-Genome Distance Calculator (GGDC) 2.1 (Meier-Kolthoff et al., 2013). The distance d(X, Y) between genome X and Y was calculated according to the formula:

$$d\left(X, Y\right) = 1 - \frac{2 \cdot I\_{XY}}{H\_{XY} + H\_{YX}} \tag{1}$$

in which, IXY denotes the sum of identical base pairs over all high-scoring segment pairs (HSPs, which are intergenomic matches), while HXY and/or HYX denote the total length of all HSPs. Heatmap was shown using the software HemI (Deng et al., 2014).

#### 16S Ribosomal RNA (rRNA) Gene-Based and Whole Genome-Based Phylogenetic Tree

Phylogenetic relationship based on 16S rRNA sequences of Acidithiobacillus strains was analyzed using MEGA v5.05 with neighbor-joining method. The robustness of clustering was evaluated by 1,000 bootstrap replicates. Additionally, the phylogenetic relationships between complete and draft genomes from A. caldus strains were estimated. We employed an online platform CVTree3 (Zuo and Hao, 2015) to construct the whole-genome based phylogenetic tree using a composition vector approach. This whole-genome-based and alignment-free prokaryotic phylogeny was validated by directly comparing our result with the taxonomy of these strains, as opposed to performing statistical resampling tests such as bootstrap or jackknife. The genome sequence of Acidithiobacillus ferrooxidans ATCC 23270 was chosen as an outgroup. Subsequently, visualization of phylogenetic tree was executed using the MEGA v5.05 (Tamura et al., 2011).

### Pan-Genome Analysis

fmicb-07-01960 December 2, 2016 Time: 15:36 # 4

Species diversity could be identified by analyzing gene repertoire across all strains of a species, i.e., the pan-genome (Tettelin et al., 2008). PanOCT v3.18 (Fouts et al., 2012) with a BLASTP all-against-all comparison of entire proteins (Evalue ≤ 1e−5 ; sequence identity ≥ 50%) was used to identify shared and unique gene content. Subsequently, annotation of core genome and strain-specific genes was implemented using BLAST against the extended COG database (Franceschini et al., 2013).

### Gene Family Evolution

Groups of orthologous sequences (orthogroups, herein referred to as gene families) in all six A. caldus strains were classified by clustering with OrthoFinder v0.4 (Emms and Kelly, 2015), using a Markov cluster algorithm. Transposable elements were excluded, given that these gene sequences might interfere with our analyses owing to lineage-specific expansions (Carretero-Paulet et al., 2015).

To analyze the evolutionary rates of gene families, we applied the developed computational program BadiRate v1.35 using a GD stochastic model (Librado et al., 2012). The gain (γ) and death (δ) rates of gene families were estimated using a branch-specific rates (GD-BR-ML) model assuming that each phylogenetic branch had its own specific turnover rate.

### Mobile Gene Elements, Insertion Sequence Elements, Transposable Elements, and Genomic Islands

IS family annotation and transposase inspection was done by BLAST comparison (E-value ≤ 1e−5 ) against the ISFinder database with manual detection of the surrounding significant search hits (Siguier et al., 2006). The program SeqWord Genomic Island Sniffer (Bezuidt et al., 2009) was implemented to identify the putative horizontally transferred elements distributed in the chromosome of A. caldus ATCC 51756. Then, the prediction of genes in the putative horizontally transferred elements was performed using the MetaGeneAnnotator (Noguchi et al., 2008). For the other chromosomes, the computational tool IslandViewer 3 (Dhillon et al., 2015), which integrates three different prediction methods including IslandPick (Langille et al., 2008), IslandPath-DIMOB (Hsiao et al., 2003), and SIGI-HMM (Waack et al., 2006), was used to predict GIs. The GC content of GI sequences was calculated using the NGS QC Toolkit (Patel and Jain, 2012). Due to the high number of contigs, A. caldus S1 was excluded from the GI prediction.

### Availability of Supporting Data

The data sets supporting our results in this study are available in the GenBank repository. These Whole Genome Shotgun projects of four newly sequenced A. caldus strains have been deposited at the DDBJ/ENA/GenBank under the accession numbers LZYE00000000 (DX), LZYF00000000 (ZBY), LZYH00000000 (S1), and LZYG00000000 (ZJ). Additionally, the versions described in this paper are version LZYE01000000, LZYF01000000, LZYH01000000, and LZYG01000000, respectively.

## RESULTS AND DISCUSSION

### Overview of the A. caldus Chromosomes

The circular chromosomes of A. caldus strains varied from 2.78 to 3.16 Mb (**Table 1**). A. caldus strains DX, ZBY, and ZJ, which were isolated from a copper mine, possess larger chromosomes than the other strains inhabiting the divergent habitats. Genome-size variations in bacteria correspond to variations in gene number as bacterial genomes are tightly packed, and most sequences are functional protein-coding regions (Mira et al., 2001). Accordingly, strains with larger genome were predicted to harbor more CDSs compared to other strains in this study. Additionally, the evaluation of quality and completeness of genome assemblies supported the reliability of pan-genome analysis, although strain S1 had relatively low genome completeness in comparison with its closely related counterparts (**Table 1**).

In all A. caldus strains, the mean percentage GC content of these chromosomal DNAs (60.90–61.72% for all six strains) was much higher than that observed for other recognized Acidithiobacillus spp., e.g., A. ferrooxidans, A. thiooxidans, and A. ferrivorans. It might be reasonable considering that A. caldus species was known as the only known mesothermophile within the Acidithiobacillales (Acuña et al., 2013), and GC content of prokaryotic genomes was positively correlated with optimal growth temperature (Musto et al., 2004, 2006).

### Evolutionary Relationship of A. caldus Strains

A phylogenetic tree based on 16S rRNA genes of Acidithiobacillus strains preliminarily demonstrated that these four newly sequenced strains in this study were taxonomically affiliated with A. caldus (**Figure 1**). To further identify the evolutionary relationships of A. caldus strains, an whole-genome-based and alignment-free phylogenetic tree was constructed (**Figure 2**). Additionally, GGDC analyses were employed to support the phylogenetic relationship. This phylogenomic tree showed that three strains isolated from the copper mine (namely, ZJ, DX, and ZBY) were clustered together (group 2 in **Figure 2**). Similarly, an earlier study reported that taxonomic clustering of six strains belonging to the genus Novosphingobium was generally influenced by their respective source of isolation (Gan et al., 2013). Further inspection revealed that the geographic distribution of strain ZBY was distinctively differed from those of the other two strains (ZJ and DX), and the genome-contentbased distance matrix implied a slight evolutionary divergence (**Figure 2**). The correlation between intraspecific divergence and geographic distribution was also observed within the closely related A. thiooxidans species by comparative genomic

accession numbers of gene sequences or genomic loci are given in parentheses. Four new strains in this study are highlighted in bold.

analysis (Zhang et al., 2016a). In the group 1 (**Figure 2**), interestingly, A. caldus SM-1 was obtained from a bioleaching reactor used for low grade gold-bearing minerals (Acuña et al., 2013), and the strain ATCC 51756 was isolated from a coal spoil; moreover, phylogenetic analysis revealed that these two strains were more closely related to each other than to the other four strains examined in this study. We therefore suspect that strain SM-1 might originally be isolated from an acidic setting similar to the habitat for ATCC 51756.

Ji et al. (2014) showed that the differences in adaptive evolution were attributable to different econiche by genetically analyzing the marine and freshwater magnetospirilla. Accordingly, we propose that environmental variation, particularly geochemical conditions, might be a determinant of genomic diversity of A. caldus strains. From an alternative perspective, it appears that geographic distribution has less of an influence on hereditary variation in comparison with econiche difference. The findings were consistent with an earlier study showing that environmental heterogeneity has relatively more

influence on microbial biogeography compared to geographic distance (Lin et al., 2013).

#### Gene Contents in A. caldus Strains

Gene prediction showed that the chromosomes of A. caldus strains contained 2,699 (ATCC 51756), 2,833 (SM-1), 2,874 (S1), 2,942 (DX), 3,017 (ZBY), and 2,984 (ZJ) predicted CDS. Functional analysis based on COG categories (Supplementary Table S1) revealed that the four most abundant functional categories within all A. caldus strains were "function unknown [S]," "replication, recombination, and repair [L]," "cell wall/membrane/envelope biogenesis [M]," and "energy production and conversion [C]." As reported by Silver and Phung (1996), high concentrations of toxic substrates such as heavy metals might cause a high rate of DNA damage. Thus, it was expected that CDS involved in COG category [L] would be abundant in A. caldus strains. Additionally, these data can also explain why this finding was distinct from previous studies analyzing the COG classification of other organisms such as marine magnetospirillum Magnetospira sp. QH-2 (Ji et al., 2014), given that the concentrations of potential toxic substrates in the extreme environment were much higher than those in the marine environment.

A previous study based on four genomes of "Ferrovum" strains highlighted the most distinct differences in interspecific metabolisms (Ullrich et al., 2016). However, in our study the assignment of CDS to the COG classification revealed that no significant differences in the number of assigned CDS were observed between the six genomes (Supplementary Table S1), probably suggesting few group-specific metabolic traits.

### Comparison of Inferred Metabolic Traits and Niche Adaptation

#### Comparison of the Central Metabolism

In light of COG assignment aforementioned, we further observed CDS related to the predicted metabolic profiles. Compared with other metabolic models reported in the literature, including carbon metabolism (You et al., 2011; Zhang et al., 2016b), nitrogen uptake (Levicán et al., 2008; Justice et al., 2014), and sulfur oxidation (Mangold et al., 2011; Chen et al., 2012; Yin et al., 2014), all strains in our study were predicted to contain numerous genes involved in central metabolism (Supplementary Table S2). The metabolic potentials of all strains were reconstructed and compared to each other for the identification of shared metabolic features as well as group- or strain-specific traits (**Figure 3**). Comprehensive analysis of these metabolism-related genes focuses on the main differences between the six A. caldus strains. As depicted in **Figure 3**, however, the evidence showed low intraspecific genetic diversity in the predicted metabolic profiles between A. caldus strains. A suite of genes involved in carbon assimilation were found in all strains. A. caldus fixes carbon dioxide via the classical Calvin–Benson–Bassham (CBB) cycle, and harbors a gene cluster predicted to encode carbon dioxide-concentrating protein (CcmK) with various copies, carboxysome shell protein (CsoS), carboxysomal shell carbonic anhydrase (CsoSCA), and ribulose-1,5-bisphosphate carboxylase/oxygenase (RuBisCO; Supplementary Table S2). Moreover, A. caldus operates a complete Embden–Meyerhof pathway (EMP) or glycolysis, pentose phosphate pathway (PPP), and incomplete tricarboxylic acid (TCA) cycle, which lacks the 2-oxoglutarate dehydrogenase complex (Valdés et al., 2008b).

With respect to nitrogen uptake, although A. caldus lacks nitrogenase directing the fixation of molecular nitrogen (Valdés et al., 2008a), assimilation of nitrate, nitrite, and ammonia plays a critical role in meeting nitrogen requirements. A. caldus utilizes nitrate or nitrite via nitrate transporter (NRT) and nitrate/nitrite transporter (Nrt). However, NRT was not present in strain ATCC 51756. Genes associated with dissimilatory nitrate reduction were identified, while a gene involved in assimilatory nitrate reduction (nasA) was absent in strain ATCC 51756. Though absent, it appears that the non-existence of those genes had little

influence on the assimilation of nitrate. Additionally, all strains share the potential to take up extracellular ammonia into the cell via AmtB transporter (**Figure 3**) under low nitrogen levels (Levicán et al., 2008), and to convert it to glutamine via glutamine synthetase.

In recent years, the sulfur oxidation system in A. caldus has been well studied (Mangold et al., 2011; Chen et al., 2012). According to reported sequences, numerous genes related to sulfur oxidation were found. Additionally, all A. caldus strains harbor genes predicted to be involved in sulfate reduction (Supplementary Table S2). Of note, the sor gene encoding sulfur oxygenase reductase, an important enzyme catalyzing a disproportionation reaction of cytoplasmic sulfur (Zhang et al., 2015), was absent in strain SM-1 (**Figure 3**). Group 2 strains lack the gene encoding the putative thiosulfate:quinone oxidoreductase. Thus, whether other alternative genes exist

in these strains needs to be studied further. Similar to the well-studied model for electron transfer of A. ferrooxidans (Valdés et al., 2008a), A. caldus potentially employs the electron transfer pathway from sulfur oxidation to (1) various types of terminal oxidases to generate a proton gradient or (2) to NADH complex to produce reducing power (Chen et al., 2012).

To some extent, investigation of genes involved in central metabolism supported the results of COG assignment that there were no obvious intraspecific differences. In other words, comparison of intraspecific genomes showed that only slight differences were observed in metabolic profiles, at least in central metabolism.

#### Response to Environmental Stress

Microbial response to environmental stresses is always a critical issue in ecological fields (Yin et al., 2015). A long-term experiment with Escherichia coli revealed complex coupling between organismal adaptation and genome evolution, which occurred even in a constant environment (Barrick et al., 2009). In the context of the six A. caldus strains, bacterial adhesion, motility, heavy metal resistance, and organic solvent tolerance were taken into account (**Figure 3**). All strains share a core set of genes potentially related to environmental adaptation (Supplementary Table S2). The presence of genes encoding extracellular polymeric substances precursors and type IV pili in A. caldus suggests a cell adhesion on mineral surface. This trait provides a reaction space between cell and mineral surface, thereby increasing the dissolution of metal sulfides (Watling, 2006; González et al., 2013). Genes assigned to COG category [N] (cell mobility) and [T] (signal transduction) were also observed, but there were few differences between these two groups (Supplementary Table S1). A full suite of genes associated with flagellar assembly were found in all strains, suggesting that A. caldus strains had the capacity to swim across environmental gradients and to colonize new sites.

Extremely acidic environments s, especially bioleaching systems, are regarded as having extremely high concentrations of soluble and potentially toxic substrates such as heavy metals, including arsenic, mercury, copper, and cadmium (Valdés et al., 2008a) and organic extractants, such as Lix984n (Zhou et al., 2012). A series of gene clusters potentially encoding functional enzymes were identified, suggesting that A. caldus has the ability to cope with high concentrations of heavy metal ions. As for organic solvent tolerance, a sixgene cluster, encoding ABC transporter ATP-binding protein, hypothetical protein, toluene tolerance protein, mce-related protein, toluene tolerance protein Ttg2B, and toluene ABC transporter ATP-binding protein, was found in all strains. Additionally, an acrAB-tolC operon potentially encoding AcrB (transporter AcrB/AcrD/AcrF family protein), AcrA (RND family efflux transporter MFP subunit), and TolC (outer membrane efflux protein) in each genome indicated that A. caldus can utilize the pumps associated with resistancenodulation-cell division protein to transfer these organic substrates.

## Pan-Genome Analysis

As shown above, numerous homologous genes associated with metabolic pathways as well as environmental adaptation were observed. To gain a deeper understanding of group- and strainspecific features, pan-genome analysis of A. caldus species was performed. A total of 4,424 CDS acquired from the four newly sequenced chromosomes plus two available chromosomes in the public database were clustered using the PanOCT. Pairwise BLAST comparisons indicated that 1,839 orthologs (41.57%) with a high percentage across all six strains were identified as the A. caldus core genome (**Figure 4**). The remaining variable 1,307 clusters were classified as the A. caldus accessory genome. Furthermore, strain-specific clusters were observed among the six A. caldus strains.

Functional assignment based on the core genome was employed to investigate the proportion of proteins in each COG category. As depicted in **Figure 4**, the core genome in A. caldus strains was commonly enriched in the COG category [M] (cell wall/membrane/envelope biogenesis; 6.36%). Additionally, our results showed that CDSs involving COG categories [C] (energy production and conversion; 6.30%) and [E] (amino acid transport and metabolism; 6.25%) were abundant. The large proportion of these genes indicated that energy utilization and uptake of nutrients in these strains might be more efficient to better adapt to the challenging environment. In other words, these findings were in line with an earlier report detailing that core genes provided functions that were essential to the basic lifestyle of the species (Medini et al., 2005).

Persistent genes encoding essential functions are stably maintained in genomes under constant selection (Nuñez et al., 2013), while dispensable or accessory genes are frequently gained or lost (Medini et al., 2005). Therefore, the accessory genome contributes to intraspecific diversity (Tettelin et al., 2008). Here, we identified many transposases by alignment of accessory genes against the NCBI-nr database (Supplementary Table S3), suggesting roles in shaping the evolution of protein families. Similarly, previously studies based on available genomes revealed that plentiful accessory genes were probably acquired by HGT (Tian et al., 2012; Sugawara et al., 2013). Additionally, it is particularly noteworthy that strain-specific genes were found to be enriched in the COG category [L] (replication, recombination and repair; **Figure 4**), thus supporting the view that the accessory genome confers selective advantages such as niche adaptation.

In particular, a total of 43 and 276 group-specific genes shared by group 1 and group 2 strains, respectively, were detected (**Figure 4**). Functional profiling based on COGs revealed that most of these predicted CDS were assigned to no COG category, probably indicating the existence of many groupspecific CDS with unidentified function (Supplementary Table S4). Further inspection underscored that the abundant genes involved in certain COG categories, including [L] (replication, recombination, and repair), [M] (cell wall/membrane/envelope biogenesis), and [P] (inorganic ion transport and metabolism), might be necessary for the group 2 strains. A reasonable explanation is that copper bioleaching heap, the habitat for group 2 strains, has high concentrations of toxic metals (Zhang et al., 2016c). Microbes in such an extreme environment might

harbor potential strategies to cope with the chemical constraints of their natural functions. Additionally, COG categories [S] (general function prediction only) and [R] (function unknown) were relatively abundant in all groups, further highlighting the role of these unknown functional CDSs in genomic differentiation.

#### Mobile and Transposable Elements

Prediction and classification of transposable elements using ISFinder indicated that a large number of IS elements, which accounted for various proportions of the total CDS in each chromosome (ranging from 1.8 to 5.8%), were randomly distributed over the chromosomes of the A. caldus strains (Supplementary Table S5). Although the types of IS families were similar to each other, their distribution and relative abundance varied with each strain. Among them, some of these IS elements were identified to cluster in flexible chromosomal regions that did not satisfy the criteria of other putative mobile elements such as GIs; these findings were consistent with those from an earlier study (Acuña et al., 2013). As stated by Bentley and Parkhill (2004), the progressive loss of gene order in a prokaryotic genome might be attributed to several events including gene deletion, IS and repeat expansion, as well as recombination or rearrangement. Given this, A. caldus SM-1 as well as ATCC 51756 might have higher genome plasticities compared with other closely related strains, mainly because of the acquisition of IS elements during evolution.

Aside from IS elements, the putative GI elements in all A. caldus strains were also identified. Results showed that several GIs ranging from 4 to 58 kb were widespread in the chromosomes of A. caldus strains (Supplementary Table S6). Additionally, most CDS in the GIs were annotated as hypothetical proteins. Further analyses showed the presence of integrases or mobile genetic elements such as transposase, thereby indicating that various putative GIs might be acquired via HGT. In light of the view that underscores the contribution of horizontal (lateral) gene transfer (HGT) in the expansion of gene repertoires of prokaryotes (Ochman et al., 2000; Gogarten et al., 2002; Treangen and Rocha, 2011), we inferred that the frequency of HGT was high in group 2 strains with larger genomes, conferring a predominant role in shaping their evolution and allowing the acquisition of novel adaptive functions. We emphasized the role of GIs in adaptation to specific lifestyles and environmental niches, considering that many GIs were highly relevant for niche-specific adaptation (Wu et al., 2011). Furthermore, numerous genes in A. caldus species might be obtained by genetic exchange as suggested by the presence of a large load of mobile genetic elements including IS elements, transposases, and GIs. Consequently,

changes in genome structure and gene copy number might provide A. caldus strains with a survival advantage for rapid adaptation and survival in highly acidic and metal-laden environments.

#### Probabilistic Analysis of Gene Family Turnover

Gene families in A. caldus strains were classified as orthogroups using OrthoFinder (**Table 2**). Our classification identified up to 3,109 orthogroups (containing two or more genes in all selected strains), which included 16,470 sequences. There were fewer genes in group 1 strains with small genomes clustered into multigene orthogroups than in group 2 strains. However, these smaller genomes contained more unassigned genes than any orthogroups compared to the others. These results appear to be explained in part by intense fractionation pressure (Carretero-Paulet et al., 2015). In other words, multigene families in smaller genomes might be under continuous deletion pressure and, as a result, these genomes tend to be smaller in comparison with their counterparts in larger genomes.

BadiRate analysis, using a full likelihood method, was applied to examine the evolutionary dynamics of gene families across the A. caldus species, and to characterize the expansion and/or contraction of genomes. The statistical framework not only estimates GD rates in a decoupled manner, using two independent parameters (γ and δ), but also explicitly takes into account certain key features in prokaryotic evolution, such as HGT. A stochastic GD-BR-ML model statistically evaluating the turnover rates demonstrated that a large number of orthologous genes frequently undergoing high gain and/or death events have evolved from ancestral genes (**Figure 2**). Particularly, gene families in A. caldus species rapidly expand through gene gain (duplication) and slowly contract through gene death (deletion or pseudogenization), indicating that the extensive recruitment of genes involved in long-term evolution confers an ecological advantage for survival and proliferation under extremely acidic conditions. Moreover, the phylogenetic branch with the higher death rate (δ = 0.616 and/or 0.808) indicated that group 1 strains with smaller genomes might be derived from free-living ancestors by the genome-reductive evolutionary process (**Figure 2**). Given that genome reduction coincided with the increase in frequency of mobile elements and repeated sequences (Moran, 2003), multiple IS elements identified in strain SM-1 and ATCC 51756 (Supplementary Table S5) might play a key role in mediating intrachromosomal recombination, thereby leading to rearrangements and gene loss. However, the dispensable genes in the above-mentioned microorganisms might suffer extensive loss and non-functionalization. The compact genomes in the given organisms can perform essential functions for cellular survival and replication, as the loss of dispensable genes has little effect on bacterial fitness, at least under certain environmental conditions (Albalat and Cañestro, 2016). Despite their smaller or near-minimal size, all reduced genomes still retain the essential gene set, and are thereby able to support cellular life both in stable and changing circumstances (Moya et al., 2009). Therefore, small genomes in group 1 strains would be more tightly packed by selective reduction, and are thus more streamlined than their larger genome counterparts.

Gene turnover in group 2 strains was also estimated. As illustrated in **Figure 2**, phylogenetic branches showed a lower gene turnover rate in group 2 strains compared to that in group 1 strains. Additionally, we found that the rates of gene death were slightly higher than the gain turnover rates. There were two possible explanations for these results. The number of genes gained from HGT as well as gene duplication events might be significant enough to account for the increase of microbial DNA content and novel functions, and play a key role in evolution (Mira et al., 2001; Navarre et al., 2006). This hypothesis may also be supported by an earlier genetic study on the evolution of Bacillus anthracis virulence, which revealed that key genes that cause anthrax in this bacterium were identified as acquired by HGT (Zwick et al., 2012). However, a conceivable explanation underlying environmentdependent conditional dispensability indicates that genes in a given species would be dispensable if they were related to certain processes that were only required in a specific untested environments (Albalat and Cañestro, 2016). Of note, it is challenging to assess which genes are regarded as dispensable or essential components by coupling genotypes with phenotypes. In view of the complexity of environmental conditions in copper mines, low deletion pressures might provide microbes with a major fitness advantage for growth in adverse environments. Furthermore, large genomes in bacteria correspond to species that have the ability to tackle various environmental stimuli (Schneiker et al., 2007). Accordingly, large bacterial genomes might have an adaptive role in the evolution of group 2 strains.

#### CONCLUSION

Six chromosomes of the extreme acidophile A. caldus were valuable resources for the investigation of genetic diversity and evolutionary adaptation. A phylogenetic tree based on chromosomal sequences of A. caldus species showed a potential correlation between genomic diversity

and geochemical characteristics. Further analysis revealed that chemical constraint in respective natural habitat might be a determinant contributing to genetic diversification. Apparently, genetic analyses indicated that gene gain and loss were both dominant evolutionary forces in the adaptive evolution of A. caldus species. During adaptation to these adverse environmental conditions, GD rates varied in different settings, resulting in genomic differentiation and speciation. The compact and streamlined genomes might undergo selective deletion pressure, whereas large genomes had been extensively recruited by intraspecific or interspecific genetic exchange. These genome-guided findings in our study, to some extent, provide novel insights into the evolutionary adaptation of A. caldus species.

#### AUTHOR CONTRIBUTIONS

XnZ, XL, and QH conceived and designed the experiments. XnZ and WD performed the experiments. XnZ analyzed the data. XnZ wrote the paper. XoZ, FF, DP, WH, and HY revised the manuscript.

### REFERENCES

### FUNDING

This work was supported by the National Natural Science Foundation of China (No. 31570113 and No. 41573072) and the Fundamental Research Funds for the Central Universities of Central South University (No. 2016zzts102).

#### ACKNOWLEDGMENTS

We thank Dr. Qichao Tu in Zhejiang University and Dr. Guanyun Wei in Nanjing Normal University for helpful discussion and suggestions. Also, we thank the National Center for Biotechnology Information (NCBI) for providing the genomic sequences of A. caldus strains ATCC 51756 and SM-1.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fmicb. 2016.01960/full#supplementary-material

evolution of the extreme acidophile "Ferrovum". Front. Microbiol. 7:797. doi: 10.3389/fmicb.2016.00797

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2016 Zhang, Liu, He, Dong, Zhang, Fan, Peng, Huang and Yin. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Genome Analysis of a New *Rhodothermaceae* Strain Isolated from a Hot Spring

Kian Mau Goh1 \*, Kok-Gan Chan2 , Soon Wee Lim1 , Kok Jun Liew1 , Chia Sing Chan1 , Mohd Shahir Shamsir 1 , Robson Ee2 and Tan-Guan-Sheng Adrian2

*1 Faculty of Biosciences and Medical Engineering, Universiti Teknologi Malaysia, Skudai, Malaysia, 2 Division of Genetics and Molecular Biology, Institute of Biological Sciences, Faculty of Science, University of Malaya, Kuala Lumpur, Malaysia*

A bacterial strain, designated RA, was isolated from water sample of a hot spring on Langkawi Island of Malaysia using marine agar. Strain RA is an aerophilic and thermophilic microorganism that grows optimally at 50–60◦C and is capable of growing in marine broth containing 1–10% (w/v) NaCl. 16S rRNA gene sequence analysis demonstrated that this strain is most closely related (<90% sequence identity) to *Rhodothermaceae*, which currently comprises of six genera: *Rhodothermus* (two species), *Salinibacter* (three species), *Salisaeta* (one species), *Rubricoccus* (one species), *Rubrivirga* (one species), and *Longimonas* (one species). Notably, analysis of average nucleotide identity (ANI) values indicated that strain RA may represent the first member of a novel genus of *Rhodothermaceae*. The draft genome of strain RA is 4,616,094 bp with 3630 protein-coding gene sequences. Its GC content is 68.3%, which is higher than that of most other genomes of *Rhodothermaceae*. Strain RA has genes for sulfate permease and arylsulfatase to withstand the high sulfur and sulfate contents of the hot spring. Putative genes encoding proteins involved in adaptation to osmotic stress were identified which encode proteins namely Na+/H+ antiporters, a sodium/solute symporter, a sodium/glutamate symporter, trehalose synthase, malto-oligosyltrehalose synthase, choline-sulfatase, potassium uptake proteins (TrkA and TrkH), osmotically inducible protein C, and the K+ channel histidine kinase KdpD. Furthermore, genome description of strain RA and comparative genome studies in relation to other related genera provide an overview of the uniqueness of this bacterium.

Keywords: *Rhodothermaceae*, *Rhodothermus*, *Salinibacter*, *Salisaeta*, strain RA, *Rubricoccus*, *Rubrivirga*, *Longimonas*

#### INTRODUCTION

The family Rhodothermaceae is placed under the order Cytophagales of the class Cytophagia within the phylum Bacteroidetes.  Rhodothermaceae is a family of gram-negative, non-sporulating, chemoorganotrophic aerobes shaped as rods or cocci, and most strains are pigmented (Park et al., 2014).

#### *Edited by:*

*Jesse G. Dillon, California State University, Long Beach, USA*

#### *Reviewed by:*

*Bradley Stevenson, University of Oklahoma, USA Xi-Ying Zhang, Shandong University, China*

> *\*Correspondence: Kian Mau Goh<EMAIL_ADDRESS>
#### *Specialty section:*

*This article was submitted to Extreme Microbiology, a section of the journal Frontiers in Microbiology*

> *Received: 13 April 2016 Accepted: 04 July 2016 Published: 14 July 2016*

#### *Citation:*

*Goh KM, Chan K-G, Lim SW, Liew KJ, Chan CS, Shamsir MS, Ee R and Adrian T-G-S (2016) Genome Analysis of a New Rhodothermaceae Strain Isolated from a Hot Spring. Front. Microbiol. 7:1109. doi: 10.3389/fmicb.2016.01109*

Rhodothermus marinus was isolated from a submarine hot spring, at a depth of 2–3 m from the sea surface during low tide (Alfredsson et al., 1988), while Rhodothermus profundi was obtained from a deep-sea (2634 m depth) hydrothermal field (Marteinsson et al., 2010). Rhodothermus obamensis was initially described as a type strain of Rhodothermus (Sako et al., 1996) but was later regarded as a synonym of R. marinus owing to high 16S rRNA gene sequence similarity, high DNA–DNA reassociation values, and similar fatty acid profiles (Silva et al., 2000). For R. marinus DSM 4252T (Nolan et al., 2009) complete genome is available in GenBank (CP001807). The genome of the sister strain R. marinus SG0.5JP17-172 (CP003029) is also available but has not been published. Draft genomes have been deposited for the other strains of Rhodothermus, which include R. profundi DSM 22212T (BioProject ID: PRJNA303571), R. marinus JCM 9785T (PRJDB841), and R. marinus SG0.5JP17-171 (PRJNA52953).

Members of the genus Salinibacter (Salinibacter iranicus, Salinibacter luteus, and Salinibacter ruber) were isolated from a hypersaline evaporating water body (Park et al., 2014). S. ruber requires a high salt concentration, exhibiting optimum growth in media with 15–25% (w/v) total salts (Antón et al., 2002). In addition, sodium chloride and magnesium chloride are vital for the growth of S. ruber (Antón et al., 2002). The complete genome of S. ruber DSM 13855T (CP000159.1) has been assembled into one contig of 3,386,737 bp, while the genome of S. ruber M8 (PRJNA48827) is available in draft form. A species of another genus, the long, rod-shaped bacterium Salisaeta longa S4-4T (= DSM 21114T ; ATTH00000000.1), was isolated from hypersaline water bodies formed from premixed Dead Sea and Red Sea water samples (Vaisman and Oren, 2009). Although few studies have been conducted on strain S4-4T , its genome has been sequenced under the Community Science Program of the Joint Genome Institute (Project ID: 404303).

Rubricoccus marinus SG-29T , a reddish, coccal bacterium isolated from shallow water of the North Pacific Ocean, is the type strain of the type species of Rubricoccus (Park et al., 2011). The average GC content of strain SG-29T (68.9 mol%, determined by HPLC) is the highest value known for Rhodothermaceae. Two years after the report of Rubricoccus marinus, the same group of researchers isolated another bacterium from deep seawater of the western North Pacific Ocean. They proposed a new genus, Rubrivirga, with Rubrivirga marina SAORIC-28T as the sole species (Park et al., 2013). Recently, a new genus of Rhodothermaceae, Longimonas, was proposed (Xia et al., 2015). Longimonas halophila SYD6T was obtained from a marine solar saltern on the coast of Weihai, China (Xia et al., 2015). The 16S rRNA gene sequence similarity between L. halophila SYD6T and S. longa is less than 92% (Xia et al., 2015). Like most other strains of Rhodothermaceae, strain SYD6T is red. At the time of writing, genome information for Rubricoccus, Rubrivirga, and Longimonas were not available.

 Because most Rhodothermus spp. can grow optimally at 65◦C, numerous attempts to mine and examine thermostable proteins from this genus have been reported. Examples of these enzymes include amylase and pullulanase (Gomes et al., 2003), α-L-arabinofuranosidase (Gomes et al., 2000), cellulase (Okano et al., 2014), cellobiose 2-epimerase (Ojima et al., 2011), chitinase (Hobel et al., 2005), and α-galactosidase (Blücher et al., 2000).

In this work, a bacterium designated as strain RA was isolated from a saline hot spring sample. Analyses of the genome, 16S rRNA gene, and housekeeping gene sequences suggested that the strain could represent a new genus of Rhodothermaceae. The purpose of this report is to describe the phenotypic characteristics of this strain and present its genome sequence.

## MATERIALS AND METHODS

#### Water Sample Analyses and Bacterial Characterization

Three water samples were collected from a hot spring near the village of Ayer Hangat (6◦ 25′ 22′′N, 99◦ 48′ 49′′E) and analyzed by Allied Chemists Laboratory Sdn. Bhd. (Malaysia) within 3 days of sampling, in accordance with American Public Health Association and United States Environmental Protection Agency guidelines.

Marine agar (Conda, Torrejón de Ardoz Madrid, Spain) was adjusted to pH 7.6 with 3 M NaOH. A 100 µL aliquot of a water sample was spread on the agar and incubated at 50◦C. After 2 days of incubation, several colonies appeared on the agar. Three distinct colonies were repeatedly streaked on the same medium to obtain pure colonies as confirmed by a single morphology and size when examined directly using a DM300 light microscope (Leica Microsystems, Germany). 16S rRNA gene sequences of these isolates were amplified using the forward primer 27F (5′ -AGAGTTTGATCMTGGCTCAG-3′ ) and the reverse primer 1525R (5′ -AAGGAGGTGWTCCARCC-3′ ) (Lane, 1991; Chai et al., 2012). Sequencing was performed at Malaysian First BASE service provider.

The isolated strains were Gram stained and then examined under a Leica DM300 light microscope. For investigation of endospore formation, endospore staining was conducted on bacterial colony obtained from a fresh culture plate as well as a week-old culture (Schaeffer and Fulton, 1933). The growth temperature range of strain RA was assessed by growing this bacterium in marine broth and were incubated at 4, 10, 37, 45, 50, 55, 60, 65, 70, and 75◦C for up to 3 days. To determine the optimal pH which support the growth of this bacterium, strain RA was grown in non-buffered marine broths with various pH-values (pH 3.0–12 with interval pH-value of 0.5). The salt tolerance of strain RA was determined in marine broth [2% (w/v) NaCl] and half-strength marine broth [1% (w/v) NaCl], at increments of 1.0% (w/v) up to 20% w/v.

Carbohydrate utilization for strain RA was measured using API 50CHB/E test strips (bioMe'rieux, Marcy-l'Étoile, France). Catalase activity and oxidase activity were determined as described by Cowan and Steel (1965). Motility of strain RA was determined using semi-solid medium and bacterial cells were inoculated with a straight wire making a single stab down the center of the tube to about half the depth of the medium, followed by incubation at 50◦C and were examined for up to 3 days. Susceptibility to different antibiotics was determined by the disc diffusion method (Kirby–Bauer antibiotic testing) on marine agar at 50◦C. Bacterial lawns were prepared by spreading the cells on plates with sterile cotton swabs dipped in colony suspensions adjusted to a 2.0 McFarland standard. The following antibiotics were tested: ampicillin, bacitracin, chloramphenicol, erythromycin, gentamicin, kanamycin, nalidixic acid, ciprofloxacin, colistin, penicillin, rifampicin, sulfamethoxazole, and vancomycin.

#### DNA Purification

Cells were scraped from marine agar plates and subjected to genomic DNA extraction using a Qiagen DNeasy Blood and Tissue Kit (Qiagen, Venlo, Netherlands), according to the manufacturer's instructions. A NanoDrop 1000 spectrophotometer (Thermo Scientific, Wilmington, DE, USA) was used to determine the purity (A260/A280 ratio) and concentration of the DNA, which were 1.98 and 1400 ng µL −1 , respectively.

#### Genome Sequencing and Annotation

Whole-genome sequencing libraries were prepared using a Nextera DNA Sample Preparation Kit according to the manufacturer's guidelines (Illumina, Inc., San Diego, CA, USA). Paired-end (2 × 100 bp) sequencing was performed on the HiSeq 2500 platform using the HiSeq Rapid SBS Kit v2 (Illumina, Inc., San Diego, CA, USA). Adapter sequence removal, quality trimming, and de novo genome assembly were conducted using CLC Genomics Workbench version 7.5 (CLC Bio, Aarhus, Denmark). The resulting sequences were then annotated with NCBI Prokaryotic Genome Annotation Pipeline version 2.10, using the "best-placed reference protein set" method, and GeneMarkS+ version 3.1, which integrates information regarding protein alignments, frameshifted genes, non-coding RNA sequences, and DNA-specific statistical patterns typical of protein-coding and non-coding regions into gene predictions (Besemer and Borodovsky, 2005). Annotation and KEGG pathway prediction for strain RA were performed using the online service Pathosystems Resource Integration Center (PATRIC; Mao et al., 2015) as well as Integrated Microbial Genomes/Expert Review (IMG/ER; Markowitz et al., 2009). VirulenceFinder version 1.5 (threshold 85%, minimum length 60%; Joensen et al., 2014) and ResFinder version 2.1 (threshold 60%, minimum length 60%; Zankari et al., 2012) were used to predict the presence of putative virulence factors and antimicrobial resistance genes, respectively.

The following genomes were downloaded from GenBank: Rhodothermus marinus DSM 4252T (CP001807.1), S. ruber DSM 13855T (CP000159.1), and S. longa DSM 21114T (ATTH00000000.1). Genome comparisons were performed using the average nucleotide identity (ANI) function of IMG/ER. The complete 16S rRNA gene sequence of strain RA was submitted to GenBank under accession number KU517707.

#### Phylogenetic Analysis

The 16S rRNA gene sequences phylogenetic tree was constructed using Neighbor-joining method with 1000 bootstrap replicates. All the methods above were conducted using Molecular Evolutionary Genetics Analysis software (MEGA, Version 6.0; Tamura et al., 2013).

#### Data Access

Sequencing data for strain RA are available online as BioProject PRJNA308615, NCBI taxonomy ID 1779382, Patric genome id 1100069.3, and IMG/ER Taxon ID<PHONE_NUMBER> (GOLD ID Gs0118059). This Whole Genome Shotgun project has been deposited at DDBJ/EMBL/GenBank under the accession number LRDG00000000. The version described in this paper is version LRDG01000000. A Fasta file of the strain RA genome is available in the following figshare link (https://figshare. com/s/60c1a70599d11684aaf7). All Supplementary Materials are available in figshare.

### RESULTS AND DISCUSSION

### Site Description

The hot spring sampled in this study is located near Ayer Hangat (AH) on Langkawi Island, Malaysia. Owing to its relatively low temperature (∼45◦C), the hot spring has been used for balneotherapy for tourists and local residents. The water in this hot spring is trapped within a man-made pit and is therefore nearly stagnant and is slightly yellowish, and there is a thin biomat (<1 cm) on the sides of the pit that is light brown, green, and yellow. Water samples taken from this site contained high concentrations of Na+ (7900 mg L−1 ), Cl− (13,800 mg L−1 ), Mg2+ (390 mg L−1 ), SO2− 4 (950 mg L−1 ), S (480 mg L−1 ), and CaCO3 (5020 mg L−1 ).

#### Cell Morphology and Biochemical Analyses

Three bacteria (designated as strains RA, RB, and RC) were isolated from a water sample taken at the AH hot spring. Based on 16S rRNA gene sequences, strains RB and RC had the closest match with Bacillus spp. and Geobacillus spp., respectively. Since strain RA had low sequence similarity to known Rhodothermaceae (see below), it was selected for subsequent analyses. The strain RA colonies appeared smooth, mucoid, and raised, with ∼1 mm margins. Furthermore, while the majority of Rhodothermaceae strains form orange-reddish colonies (Park et al., 2014), the colonies of strain RA were cream colored (**Table 1**). Meanwhile, light microscopy analysis demonstrated that this isolate is a gram-negative, rod-shaped bacterium that can be arranged in pairs and chains or as lone cells. Rhodothermus spp. are known to occur singly, but not in chains or filaments (Park et al., 2014). Strain RA was motile and catalase positive and it grew well in marine broth and half-strength marine broth but failed to grow on Luria–Bertani, Reasoner's 2A, Mueller–Hinton, or nutrient agar, with or without supplementation with 2% (w/v) NaCl. While the microorganism grew at 37–60◦C and pH 5.0–9.0, optimal growth occurred at 50–60◦C and under circumneutral conditions. Strain RA was unable to grow at 65◦C. Furthermore, strain RA exhibited growth in marine broth supplemented with 1–10% (w/v) NaCl. On marine agar, strain RA was susceptible to erythromycin (zone of inhibition of 25 mm around a 5 µg disc), rifampicin (25 mm TABLE 1 | General features and information regarding the sequencing of the strain RA genome, according to the MIGS (Minimum information about a genome sequence) recommendations.

zone around a 5 µg disc), and penicillin (35 mm zone around a 10 U disc). Lastly, like R. marinus DSM 4252T , strain RA is a non-spore-forming bacterium.

#### Phylogenetic Relationships

For phylogenetic analysis of strain RA, we sequenced the entire 16S rRNA coding region (1516 bp) of the microorganism and compared the resulting sequence against the GenBank and EzTaxon-e databases (Kim et al., 2012). These analyses indicated that the closest relative (99% identity, 98% coverage) of strain RA is the uncultured bacterial metagenome clone KSB113 (JX047086), which was identified in a marine hot spring at Kalianda Island, Indonesia (unpublished report). Moreover, the genera closest to strain RA (86.5–89.3%) are those of the family Rhodothermaceae: Rhodothermus (Alfredsson et al., 1988; Sako et al., 1996; Marteinsson et al., 2010), Salinibacter (Antón et al., 2002; Makhdoumi-Kakhki et al., 2012), Salisaeta (Vaisman and Oren, 2009), Rubricoccus (Park et al., 2011), Rubrivirga (Park et al., 2013), and Longimonas (Xia et al., 2015; **Table 2**, **Figure 1**). The sequence identities of housekeeping genes of strain RA and species of the family Rhodothermaceae are also low, for instance (i) recA, 83% to R. marinus DSM 4252T ; (ii) rpoD, 85% to R. marinus DSM 4252T and 82% to S. ruber DSM 13855T ; and (iii) gyrB, 79% to R. marinus DSM 4252T and 79% to S. ruber DSM 13855T .

We then compared the available genome sequences of several of these strains with the draft genome of strain RA (**Table 2**). In addition to exhibiting low 16S rRNA gene sequence identity (89.3%), strain RA and R. marinus DSM 4252T had a low average ANI-value of 73.28 (**Table 2**). The ANI-value between strain RA and R. marinus SG0.5JP17-171 or R. marinus SG0.5JP17-172 was 73.23, while that between strain RA and R. profundi DSM 22212T was 70.66. Meanwhile, the average ANI-values for strain RA when compared with S. ruber DSM 13855T and S. longa DSM 21114T were 71.85 and 71.47, respectively. Collectively, these data strongly suggest that strain RA is a novel strain of Rhodothermaceae.

### Genomic Features of Strain RA

The draft genome of strain RA is 4,616,094 bp in length (**Table 3**), and the largest contig is 819,789 bp. The N75, N50, and N25-values are 89,274, 152,113, and 242,570 bp, respectively. The average coverage for the 91 contigs obtained was 120 fold. Notably, the genome of strain RA is larger than those of Rhodothermus marinus (3.3 Mbp), S. ruber (3.5–3.8 Mbp), and S. longa (3.1 Mbp; a draft genome with three contigs; **Table 2**). Furthermore, the genome is predicted to contain 3680 coding DNA sequences (CDS) and 50 rRNA genes, and it has an average GC content of 68.28%, which is markedly higher than that of R. marinus DSM 4252T (64.3%; Nolan et al., 2009), S. ruber DSM 13855T (66.2%; Mongodin et al., 2005), and S. longa DSM 21114T (63.5%) but comparable to that of R. marinus KCTC 23197T (68.9%; **Table 2**). Among the CDS, the functions of 2815 (76.49%) sequences could be predicted, while 879 (24.38%) of the total CDS are annotated as enzymes. Protein-coding genes related to clusters of orthologous groups (COGs) are shown in **Figure 2**.

Using the Phylogenetic Distribution of Genes function in IMG/ER, we determined the best hits for protein-coding genes in strain RA at a 60% identity cutoff. A total of 1408 genes (39% of the protein-coding genes) were found to be similar to genes in the phylum Bacteroidetes. The majority of these sequences are affiliated with class Cytophagia and unclassified classes, while only 52 genes are related to Bacteroidia, Flavobacteriia, and Sphingobacteriia. Table S1 shows that 1987 genes (54.7% of the protein-coding genes) were not assigned to any phylum, as they are specific to strain RA. When all protein-coding genes were examined at a 90% sequence identity threshold, 99.75% (3621) of the genes were classified as unassigned. This suggests that the gene sequences of strain RA are more diverse than those of its Rhodothermaceae counterparts (data not shown).

2|GeneraldataandcomparisonofstrainRAandother*Rhodothermaceae*

genera.

TABLE

Consistent with its observed motility, strain RA carries the genes necessary for complete flagellum assembly. Furthermore, genomic analysis using the antiSMASH predictor program (Blin et al., 2013) identified a gene cluster encoding a type III polyketide synthase. The key enzyme within this cluster has a structural domain related to chalcone synthase; thus, strain RA might use this enzyme to produce certain types of aromatic ketones, including chalcones, which possess important properties for industrial use (Singh et al., 2014).

### Comparison of the Genome of Strain RA with those of Other Rhodothermaceae Strains

At the time of writing, the genomes of three R. marinus strains and two S. ruber strains, as well as the draft genome of S. longa DSM 21114T , were publicly available. Unless otherwise specified, we chose to compare the genome of strain RA with the complete genome sequences of R. marinus DSM 4252T and S. ruber DSM 13855T .

RAST analysis identified 215 and 221 strain RA genes that are not carried by R. marinus DSM 4252T and S. ruber DSM 13855T , respectively (Table S2), including genes encoding DnaK-related proteins, cytochrome C oxidase subunit CcoN, superfamily II DNA/RNA helicases (SNF2 family), primosomal protein N′ (replication factor Y) superfamily II helicase, β-ureidopropionase, fucose permease, gluconolactonase, α-L-fucosidase, malto-oligosyltrehalose synthase, α-1,2-mannosidase, endo-1,4-β-xylanase, endogluc anase, choline-sulfatase, putrescine transport ATP-binding protein PotA, arylsulfatase, and sulfate permease. Of these, strain RA likely maintains the genes for sulfate permease and arylsulfatase to withstand the high sulfur and sulfate contents of the AH hot spring. Furthermore, the presence of these genes suggests that strain RA plays a role in the recycling of sulfur within the hot spring. Two cholesterol oxidase genes were also present in strain RA but not in closely related genera. Notably, cholesterol oxidases catalyze the conversion of cholesterol to cholest-4-en-3-one, which has been proposed as a potential drug candidate for treatment of amyotrophic lateral sclerosis (Bordet et al., 2007).

KEGG metabolic pathway analyses indicated that the starch and sucrose metabolism pathways of strain RA are similar to those of Rhodothermus and Salinibacter spp. According to these analyses, strain RA is capable of hydrolyzing starch to maltose, glucose, or trehalose and converting starch to glycogen. Meanwhile, all three genera encode the enzymes (EC <IP_ADDRESS> and <IP_ADDRESS>) necessary for hydrolysis of cellulose to glucose Figure S1.

Likewise, R. marinus DSM 4252T , S. ruber DSM 13855T , and strain RA each encode a type I system (TolC), a Sec-SRP pathway (SecD/F, SecY, SecA, YidC, FstY, Ffh), and a Tat pathway (TatC) for protein secretion. However, while strain RA also carries genes encoding type II (GspD, GspE, GspF) and type VI (VgrG, Hcp, IcmF, DotU, ClpV) secretion systems, virulence factor prediction analyses using the PATRIC database

#### TABLE 3 | Genome statistics of strain RA.

(Wattam et al., 2014), VirulenceFinder 1.5 (Joensen et al., 2014), and ResFinder 2.1 (Zankari et al., 2012) indicated that this strain is likely non-pathogenic.

Certain ABC transporters encoded by strain RA, including those specific for phosphate (PstS, PstA, PstB, PstC), lipoproteins (LolC, LolE, LolD), and iron-siderophores (FhuB, FhuC, FhuD), are potentially produced by R. marinus DSM 4252T and S. ruber DSM 13855T . In addition, there were multiple transporters that were encoded in the genomes of strain RA and R. marinus DSM 4252T but not in that of S. ruber DSM 13855T , including transporters for molybdate (ModA, ModB), iron(III; AfuA, AfuB), zinc/manganese/iron (TroA,

TroB, TroC, TroD), raffinose/stachyose/melibiose (MsmK), sorbitol/mannitol (SmoK), cellobiose (MsiK), chitobiose (MsiK), arabinooligosaccharide (MsmX), and glucose/mannose (MalK; Table S2). Meanwhile, strain RA also encodes an RbsABC transporter for uptake of ribose/xylose. Lastly, strain RA contains genes encoding unique CusS, CusR, CusB, and CusA proteins that are predicted to play a role in copper and/or silver efflux, which are absent from all currently known Rhodothermus and Salinibacter spp.

### Temperature Adaptation Genes of Strain RA

The average water temperature of the AH hot spring during sampling was 45◦C; however, in vitro analysis demonstrated that strain RA is capable of thriving at temperatures up to 60◦C. Meanwhile, analysis of the composition of the spring water detected high levels of NaCl. We therefore endeavored to elucidate the mechanisms by which strain RA withstands thermal and osmotic stresses by characterizing the stress-related genes within the annotated genome.

The linear polyamines putrescine and spermidine have been associated with the thermophilicity of thermophiles (Takahashi and Kakehi, 2010; Goh et al., 2014). Notably, in addition to harboring the genes necessary to synthesize putrescine via Ncarbamoylputrescine amidase, strain RA appears to be capable of taking up putrescine and spermidine from the environment using the PotA and PotC transporters. In contrast, the genes encoding PotA and PotC are not present in the genomes of R. marinus DSM 4252T and S. ruber DSM 13855T . A unique sequence in strain RA was annotated as a primosomal protein N′ (replication factor Y) superfamily II helicase gene, and the corresponding protein sequence exhibited >58% sequence identity to PriA helicases deposited in GenBank. Several of these helicases have been proposed to play important roles in the thermostability of certain thermophiles (Goh et al., 2014). As mentioned earlier, the GC content of the strain RA genome is higher than that of most other genomes of Rhodothermaceae. In some cases, high genome and tRNA GC% levels are associated with increased optimal growth temperatures (Trivedi et al., 2005). Yet, the average tRNA GC content (61.4%) of strain RA is slightly lower than that of the rest of the genomes (68.3%). Additionally, the strain RA genome contains CDS that exhibit sequence similarity to genes encoding chaperones (DnaJ, DnaK, HtrA, metal chaperone) and the known heat shock proteins GroES and GroEL (HSP-60 family co-chaperones), and it contains three copies of the gene encoding the cold shock protein CspA. These CspA proteins exhibited the highest level of sequence identity to counterparts from R. marinus DSM 4252T (67, 78, and 80%, respectively).

### Osmotic Stress Adaptation Genes of Strain RA