CRISPR arrays were extracted from genomes using CRISPRFinder (Grissa et al., 2007) online server which performs BLAST against dbCRISPR (CRISPR database; last updated on 2017-01-02). CRISPRFinder further classifies the identified CRISPR arrays as true or false based on whether or not they are associated with CRISPR associated genes (Cas) respectively. Cas genes were annotated using CRISPRone (Zhang and Ye, 2017). CRISPRs lacking Cas genes in the vicinity were designated as false/questionable CRISPRs. Only true CRISPRs were selected for analyses. CRISPR arrays have two components: repeats and spacers; both of which were analyzed to study the evolution and probable viral diversity respectively. Classification and clustering of CRISPR repeats and repeat-based Cas gene predictions were undertaken using CRISPRmap, a comprehensive cluster analysis method (based on Markov clustering) which clusters conserved sequence families and potential structure motifs (Lange et al., 2013). Repeats were classified based on 40 conserved sequence families and 33 probable structural motifs. Further, 24 families and 18 structural motifs were considered for the construction of repeat cluster maps. For prediction of potential viruses most frequently associated with Thermus genomes, spacer sequences from all genomes were extracted and BLAST against viral GenBank database of NCBI (Deng et al., 2007) with a threshold e-value of 1. For better stringency, among all matches, only those having 100% identity of more than 20 nucleotides were considered as valid hits.

Genes involved in imparting competence (16 genes) to T. thermophilus HB27 were used as reference for extracting competence associated genes from individual genomes using BLAST. PilA1-A4 genes were aligned using Hirschberg (KAlign) algorithm (Lassmann and Sonnhammer, 2005). Visual alignment consensus was built at 70% threshold. Relationships among PilA1-A4 genes were inferred by PhyML (Guindon et al., 2009) maximum-likelihood method using HKY85 substitution model (Hasegawa et al., 1985). Median size estimations were made using boxplot function in R (https://cran.r-project.org/mirrors.html).

The genome of T. parvatiensis strain RL was initially assembled into three contigs (totally 2,066,435 bp, 1,886,121 bp, 159 Kbp and 21 Kbp) with G+C content of 68.5%. The 21 Kbp contig, mapped onto the chromosome (BLASTn and mapping with previously generated assembly; Dwivedi et al., 2012). The entire 21 Kbp region, seems to represent an integrated plasmid or a large genomic island incorporated into the genome, based on the annotation of mostly hypothetical genes and transposable elements among genes identified. This was supplemented by differences in mean G+C% of this region (67.3%) as compared to the rest of the genome (68.5%). In an attempt to circularize the largest contig, its ends were aligned against each other and an overlapping region of 13,300 nt was removed. Similarly, in order to circularize the 159 Kbp contig, an overlapping region of 15,853 bp region was removed from the ends of the contig. Finally, two replicons were reconstructed: a chromosome (1,872,821 bp) and a megaplasmid (143,277 bp) (Figure 1) (total size: 2,016,098 bp) (Table 2). The chromosomal origin of replication was located at 158,478–158,778. A total of 12 DNA boxes with consensus sequence of TGTGGATAA were identified spanning the 301 nt OriC region. The total number of predicted coding sequences were 2,383. The genome was found to harbor two rRNA operons and 54 tRNAs and tmRNAs. COG functional category assignment placed a large number of genes to amino acid transport and metabolism (11.04%), general function prediction (12.95%), energy production and conservation (7.54%) and translation, ribosomal structure and biogenesis (7.40%). A number of genes were classified into the unknown function category (7.78%). A large proportion of the genome is strictly attributed to genes needed by the organism for essential cellular processes. T. parvatiensis thrives at a high arsenic concentration (140 ppb). We investigated the presence of arsenic resistance mechanisms in this thermophilic organism. Arsenate reductase gene arsC (1 copy), arsenic efflux pump protein arsB (2 copies) and arsR transcriptional regulator (4 copies) were identified. The above genes belong to the ars (arsenic resistance) operon responsible for the efflux of As(III) out of the cells (Yang and Rosen, 2016). Genes for the oxidation of arsenic (aox operon) were not identified. The mechanism of arsenic detoxification in T. parvatiensis thus involves extrusion of arsenic out of the cells (rather than oxidation).

T. parvatiensis strain RL uniquely harbored three integrated phages in its genome (Supplementary Figure 1)—two on the chromosome and one on the plasmid. This was in contrast to all its close phylogenetic neighbors taken into account in this study namely, T. thermophilus HB8, T. thermophilus HB27, T. thermophilus JL18, and T. thermophilus SG0.5JP17-16 in which no phages could be identified. The first phage (20.3 Kbp) on the chromosome revealed the presence of phage structural proteins such as tail assembly protein and coat protein along with three heat shock proteins which can be directly implicated to the environment, i.e., hot spring water (surface water temperature >95°C). The phage region was associated with two hybrid histidine kinases, which are implicated in two-component regulatory system (Khorchid and Ikura, 2006). The second chromosomal phage (23.2 Kbp) was annotated and revealed a probable integron with attL and attR sites along with integrase encoding gene and flanking tRNA. The gene cassette of this integron had 23 hypothetical proteins and 16 viral proteins. These results also indicate that this integron might denote a super-integron with 37 ORFs captured. Phage regions have been known to play a role in horizontal gene transfer by specialized as well as generalized transduction (Touchon et al., 2017). Genes harbored on the integrated phage regions corresponding to two-component system and heat shock proteins reflect the dispersal of these genes might be an active phenomenon among the population.

Plasmids, including small plasmids, as well as large megaplasmids of up to 440 Kbp are known to be present in Thermus genomes (Table 1). The plasmids of Thermus are known to be the center of plasticity, harboring genes for mobile elements, transposons and a number of biosynthetic clusters (Henne et al., 2004; Bruggemann and Chen, 2006). The megaplasmid pTP143 of T. parvatiensis, contained 181 coding sequences. A large number of the genes harbored on the plasmid were, however, genes belonging to integrases, transposases, mobile elements and hypothetical protein coding genes (81 genes, constituting 44% of the total plasmid genes) (Figure 1) which denote the plastic nature of the plasmid. A low coding density (87.65%) was observed for pTP143, as compared to the chromosome (94.12%). Genetic analysis of T. parvatiensis megaplasmid pTP143 revealed a cobalamin biosynthetic cluster, a denitrification cluster and a carotenoid biosynthesis cluster (responsible for imparting the yellow pigment to the organism). A thermophilic lifestyle demands a robust DNA repair system. Genes required for thermophilic existence most suitably correspond to an elevated number of DNA repair genes (Bruggemann and Chen, 2006). Consequently, a recQ helicase, reverse gyrases, photolyase phrB, sbcC, and sbcD nucleases (implicated in deleting hairpin structures) were found on the megaplasmid. Genes related to stress response, surE (involved in nucleic acid pool maintenance; Proudfoot et al., 2004) and cytochrome P450 (Kelly and Kelly, 2013) were also identified. A number of transcriptional factors known to modulate stress conditions were found on the plasmid. These included a transcriptional regulator IcIR involved in regulation of responses to quorum sensing and toxic stress (Molina-Henares et al., 2005). A transcriptional regulator Crp/Fnr known to be responsive to environmental changes (Körner et al., 2003; Zhou et al., 2012), such as oxidative stress, carbon dioxide concentrations and heavy metal impositions was annotated. Crp/Fnr regulators act by regulating the expression of genes involved in alleviating the respective stress conditions. MerR, a heavy metal modulating transcriptional regulator (Brown et al., 2003), which activates promotors of genes in response to heavy metal influx was annotated. T. parvatiensis harbors a plastic plasmid with mobile and hypothetical gene components. Constituting genes for DNA repair, stress response and transcriptional regulators, we believe that the megaplasmid has an indispensable role to play for the thermophilic survival of T. parvatiensis. Not only essential for survival, the megaplasmid demonstrates a potentially crucial role in communicating and modulating temperature stress via an appropriate response carried out by the transcriptional regulators it harbors.

All species belonging to the genus Thermus, as described in this study, have been isolated from thermophilic environments (mostly hot spring waters) from all over the world, i.e., from USA (Murugapiran et al., 2013), Japan (Henne et al., 2004), India (Dwivedi et al., 2012), South Africa (Gounder et al., 2011), China (Yu et al., 2015; Zhou et al., 2016), Malaysia (Teh et al., 2012), New Zealand (Hudson et al., 1987; Mefferd et al., 2016), and Iceland (Chung et al., 2000). This illustrates that the genus is spread across the globe, thriving at the most extreme environments. Being confined to stressed niches, these G+C rich (64.8–69.4%) organisms possess small genomes, ranging from 2.01 Mb of T. parvatiensis to 2.5 Mb of T. tengchongensis YIM 77401 (mean genome size: 2.25 Mb). Considering the powerful evolutionary forces that have been constantly shaping their genomes, in geographically diverse niches, the genome size of Thermus shows low variability (Table 1). Already known to shed and rearrange genes that are not uniquely essential for survival, the members of the genus have been successful in maintaining their genome sizes close to the average genome size of the genus. A large part of the genomes is comprised of genomic islands, varying from 0 to 5.85% (Supplementary Table 4).

To demonstrate the extent of genomic shuffling, synteny maps of the 10 complete genomes were generated (Figure 3). T. parvatiensis and SG0.5JP17-16 revealed a conserved organizational synteny with each other and lack of inversions and rearrangements. Synteny was conserved among members of T. thermophilus group, however, large blocks of inversions were observed in relation to strains CCB_US3_UF1, T. scotoductus, T. aquaticus, T. oshimai, and T. brockianus. These demonstrate the huge genome wide rearrangements occurring at the genus level (Figure 3). The same observation was reflected in dot-plot comparisons of T. parvatiensis with other members (Supplementary Figure 2).

The genes on Thermus plasmids are known to undergo vast rearrangement events. In some cases, they have been observed to shift from the plasmid to the chromosome and get stabilized there. In other cases, large plasmids with a majority of non-essential but potentially benefit imparting clusters have been discerned. The former case has been observed in T. scotoductus and CCB_US3_UF1, both of which have small plasmids as most of the genes got incorporated on the chromosomes, leaving plasmids with diminished configurations. pTSC8 (8,383 bp), for example, has discarded many non-essential genes like cobalamine biosynthesis pathway genes, plasmid stability genes and chromosome partitioning genes, but retained genes for aerobic and anaerobic respiration to attain a much more compact conformation (Gounder et al., 2011). The latter situation has been observed for plasmid pVV8 (81,151 bp) which was found to be enriched in phosphonate metabolism genes, which are not of common occurrence in Thermus genomes (Ohtani et al., 2012). To evaluate this trend in pTP143 and other plasmids, we surveyed genetic clusters commonly observed on Thermus plasmids. We observed the presence of advantageous traits such as cobalamin biosynthesis on pTP143 and pTHEOS01 (T. oshimai plasmid) and nitrate reduction on pTP143, pTHEOS01 and T. thermophilus JL-18 plasmid. To assess possibly laterally acquired regions on pTP143, we performed a comparative mapping of plasmids based on BLASTn comparisons (identity cut-off: 80, coverage cut-off: 80, e-value: 1e-30) across all the Thermus plasmids. The plasmid pTP143 of T. parvatiensis showed highest identity to the plasmids pTHTHE1601 (identity: 99.3%, query cover: 79%), pTTJL1801 (identity: 97.5%; query cover: 78%), pTT27 of HB27 (identity: 99.0%; query cover: 47%) and pTT27 of HB8 (identity: 98.9%; query cover: 49%). A region on pTP143 that failed to show any homologous regions with any of the plasmids was analyzed as a putatively horizontally transferred region (Supplementary Table 5). This region harbored 3 genes for mobile element proteins, 9 hypothetical protein coding genes, sbcC and sbcD (hairpin structure resolving nucleases) (Figure 4). The genes present on this locus were identified using BLASTp (identity cut-off: 70%, e-value cut-off: 1.00E-15) on the chromosomes of T. aquaticus, T. brockianus, HB8, SG0.5JP17-16, JL-18, CCB_US3_UF1 and even on the chromosome of T. parvatiensis.

The genes on the plasmids of Thermus are under movement and involved in shifting from the plasmid to the chromosome, in order to get stabilized. In case of Thermus species, plasmids perhaps contribute more to the flexibility of genome by either acquiring or shedding of the genes. Incoming genes/pathways may first be incorporated on the plasmid and later be stabilized either on the plasmid itself or through integration into the chromosome. The presence of a genomic island on plasmid pTP143 with multiple mobile elements suggests that it might get mobilized soon. In the process, some other genes may also get shifted from the plasmid to the chromosome. Across Thermus genomes, this process has led to either stabilization of the megaplasmid with megaplasmids acquiring a vast majority of genes (SG0.5JP17-16), or it has led to streamlining of the plasmid (T. scotoductus).

The core genome denotes the conserved functions, whereas, the pan genome denotes the entire genetic potential of a group (Tettelin et al., 2005). The total number of genes constituting the core genome for the genus Thermus were 1177. Functional annotation of genes constituting the core genome placed a high number of genes into categories coding for amino acid metabolism and transport (12.15%), translation (10.07%), energy production and conservation (7.75%) and coenzyme metabolism (6.07%), designating these as the conserved functions specific for the genus. The pan genome was estimated at 5188 and constituted accessory genes and unique genes (singletons). Accessory genes are the ones whose orthologs are present in two or more genomes, but not all the genomes. Singletons are genes that are unique to just one genome out of all those compared. The variability in accessory genome depicts the flexibility of the genome structure. Accessory gene number varied from 661 to 1,166 (mean: 946). High number of accessory genes were observed in T. tengchongensis (1,166), T. brockianus (1,094), JL-18 (1093), T. oshimai (1,053), T. scotoductus (1,032), SG0.5JP17-16 (1027), and T. aquaticus (1003) (Figure 5B). Variability in the number of unique genes (singletons) was observed from 30 to 229 genes (mean 98). A large number of singletons were identified in T. filiformis (229), T. islandicus (212), T. tengchongensis (175), and T. aquaticus (154) (Figure 5B). Out of the genomes having high number of accessory genes (higher than mean), 55.5% of the genomes had high number of singletons too. These genomes were T. aquaticus, T. brockianus, T. oshimai, T. scotoductus, and T. tengchongensis. The accessory and unique components of the genus were enriched in genes belonging to carbohydrate metabolism and transport (8.17 and 8.76% respectively), replication and repair (8.17 and 8.76% respectively), inorganic ion transport and metabolism (5.85 and 5.69% respectively) and signal transduction (5.71 and 5.40% respectively) which reflect the diverse functional counterparts harbored by these organisms. The pan genome of the genus was estimated as an “open” pan genome because a plateau was not observed (Figure 5A) after addition of all genomes to the pan genome plot. Addition of more genomes to the group will lead to expansion of the pan repertoire (Rouli et al., 2015). The genus has thus maintained a conserved core genome, but an expansive and sundry pan genome. The flexibility of the genomes is explained by the high influx of genes into these organisms via horizontal gene transfer. Although, a number of features might get recruited and incorporated in the genome, only those that have a survival benefit for the organism will be retained and the rest will be discarded by the highly active rearrangement events that are continually taking place in these genomes. In due course of time, extensive rearrangement events have led to the establishment of those genomic features that have benefitted the organism for better survival; other dispensable elements were either shed off or transferred to the plasmid. Overall, genome-wide differences and anticipation in accordance with specific gene repertoire required at a niche (niche specialization) can be considered to be driving forces in the evolution of the genus Thermus.

In order to designate strain specific, potentially utilizable attributes of T. parvatiensis, we identified MGIs of T. parvatiensis. The term MGI encompasses those regions of the genome that are identified by mapping metagenomic data from an environment against the genome of an organism isolated from the same niche (Pašić et al., 2009). These regions stand out as “gaps” with little or no reads corresponding to these regions, thus highlighting the strain specific potential that these organisms have accumulated in contrast to the environmental counterparts. Such an analysis was performed for T. parvatiensis by aligning raw metagenomic data from hot spring water (Manikaran, India) onto the replicons (Figure 6). T. parvatiensis chromosome showed a relatively high coverage of reads (4 ×) with five regions with no coverage of reads (MGIs) (Figure 6). T. parvatiensis plasmid pTP143, however had thin coverage of reads (1.9 ×) and five MGIs. Thus, plasmid pTP143 seems to have accumulated more strain specific variations which denote high flexibility of the plasmid. This data further states that the chromosome in case of Thermus is more or less stable in terms of the genes it harbors, however, much influx and rearrangements occur via plasmid. The chromosomal MGIs harbored genes for arginine biosynthesis, iron-sulfur cluster assembly proteins, transcriptional regulators, two component system genes and hypothetical genes (Supplementary Table 6). A MGI (4,811 bp) on the plasmid was found to harbor genes specifically involved in environmental response to stress in general and oxidative stress in particular. These genes included radical SAM (S-adenosylmethionine) domain heme biosynthesis protein (heme is a co-factor for hemoproteins that functions in electron transport chain), cytochrome c552 which is particularly involved in electron transport at low aeration, peptide methionine sulfoxide reductase (MsrA) known to protect against reactive oxygen and nitrogen species (Weissbach et al., 2002). Cobalamin biosynthesis genes were largely detected on the plasmid MGIs, including uroporphyrinogen-III methyltransferase, BluB, adenosylcobalamine-phosphate synthase and CbiG. Other prominently represented genes included plasmid stability genes (ParB, Soj, and StbB) and DNA repair genes (reverse gyrase, sbcC and sbcD). These genes implicate the conservation of low aeration oxidation response, plasmid stability and DNA repair as strain specific features. In order to get insights about the prevalence of these specialized regions in other genomes, we performed BLASTn (identity > 95%, query coverage > 95%, e-value <1.00E-30) of MGI regions with other Thermus genomes. Whereas, genes prevalent on MGI1 of the chromosome did not show significant identity with other members of the genus, MGI2-5 of the chromosome were found to have homologous counterparts in JL-18, SG0.5JP17-16, HB27, and T. brockianus. Plasmid MGIs 1, 3, 4, and 5, similarly could be identified on JL-18, SG0.5JP17-16, HB8, HB27, however, plasmid MGI2 of pTP143 was unique in this respect and significant similarity could not be observed with other members. The genes thus annotated on MGIs, were conserved within the T. thermophilus group indicating close adaptive relatedness. Chromosomal MGI1 and plasmid MGI2 were identified as strain specific MGIs for T. parvatiensis. These attributes suggest the conservation of features that are not directly implicated with the niche, but are retained in the genome as anticipated survival benefits.

Analysis of CRISPRs among Thermus species was performed to get insights into prevalence of viral defense system within the genus. All genomes harbored CRISPR loci, except T. antranikianii. Nine CRISPR arrays were identified in T. filiformis, HB8, HB27 and T. igniterrae, followed by T. aquaticus which had 7 CRISPR arrays (Table 3). T. parvatiensis, on the other hand, was found to harbor only 1 true CRISPR array. One questionable array was detected in T. islandicus. The CRISPR cassette/array number varied from 0-9 with DR length variation from 25 to 37 wide all genomes analyzed (mean CRISPR array count per genome: 5.529). The number of spacers harbored within CRISPR arrays denote the frequency of viral invasions. T. oshimai was found to carry the highest number of spacers (134 spacers). A high frequency of viral attacks in these thermophiles was indicated by a high mean number of spacers harbored by each genome (mean: 72.7 spacers per genome). Five Thermus genomes (29.4%) harbored >100 CRISPR spacers. Within T. oshimai itself, two CRISPR arrays harboring large number of spacers were uncovered. One of the CRISPR arrays harbored 88 copies of a single DR consensus (CGGTCCATCCCCACGGGCGTGGGGACTAC; DR length: 29 bp), with an equally high number of spacers. The other array (DR consensus = CTTTGAACCGTACCTATAAGGGTTTGAAAC; DR length: 30 bp) had acquired 67 spacers. On the contrary, within the same genus, an array containing only 3 spacers was also observed. This indicates high variation among CRISPR elements within the species. Cas genes were extracted from true CRISPRs and annotated (Table 3). The Cas system in Thermus was found to be composed of genes belonging to types I, III, and IV (Supplementary Figure 3). Universal type genes cas1 and cas2 were identified on CRISPR loci. Apart from this, the class I effector cas3 which is responsible for the helicase and DNase activity was also annotated (Table 3). Thus, cleavage of foreign entities entering Thermus genomes is brought about by type I, III, and IV mediated action. Analysis of repeats based on sequence and structure (Figure 7) was performed. Sequence families 1 and 18 were predominant (spotted in 7 and 6 genomes respectively) (Supplementary Table 8). Among structural motifs (based on the classification of motifs into 33 groups), motifs 25 and 5 were most represented (8 and 6 genomes respectively). Motif 6, motif 20, motif 23, and motif 31 were least represented (1 genome each) (Supplementary Tables 7, 8). Additionally, on the basis of repeat-cas binding projections, probable cas genes harbored by Thermus genomes were predicted and found to belong to types I, III, and IV. Separate cluster trees were constructed for structure motifs (based on the classification of structural motifs into 18 groups) to denote the placement of Thermus repeats among all consensus repeats present in the database. A majority of repeats (10 repeats) mapped onto motif 1 and occupied a close phylogenetic position within the cluster tree (Supplementary Figure 4). In order to discern the phages most frequently attacking Thermus genomes, we analyzed spacer matches with viral sequence database (Supplementary Table 9). Positive inferences were based only on the results that satisfied the stringency criteria. A number of phages of different families were detected to infect Thermus species (Figure 7B). Among these, phages of families Sphaerolipoviridae, Siphoviridae, Myoviridae, and Herpesviridae were identified as the most prominent bacteriophages. Earlier Siphoviridae, Myoviridae, Inoviridae, and Tectiviridae have been reported from Thermus species (Yu et al., 2006). However, our analysis revealed the dominance of Sphaerolipoviridae (28.1%), with known thermophilic phages P23-77 and IN93 being prominently detected. A detailed list of viruses, the invasion memory of which is incorporated within Thermus CRISPRs is mentioned in Supplementary Table 9. Our analysis thus reveals the abundant presence of defense mechanism and frequent viral encounters in this genus. Even though a large number of spacers (1,223) were analyzed, only 62 (5.07%) could be assigned significant matches to the viral database.

CRISPRs constitute the characteristic prokaryotic and archaeal adaptive as well as inheritable immune system composed of short repeat sequences (direct repeats/DRs) interspaced with short segments of nucleotides known as spacers. Spacers represent the memory of past invasions by foreign genetic elements like viruses (Barrangou et al., 2007) or plasmids (Marraffini and Sontheimer, 2008). Spacers are incorporated into CRISPR loci whenever a bacteriophage infects the organism. This way, a CRISPR array can be considered as a library of past viral invasions faced by an organism. CRISPR arrays are associated with CRISPR associated genes (Cas genes), which are present in the vicinity of CRISPR arrays. Together with virus specific spacers, Cas genes encode an arsenal of proteins and RNAs, which in conjugation, destroy the foreign element, the next time it invades (Mojica et al., 2005; Barrangou et al., 2007). The CRISPR system is thus, a defense mechanism against bacteriophage invasions on the bacterial genome. CRISPR arrays that lack the requisite Cas genes in their vicinity are known as questionable or false CRISPRs. Predominance of viruses in thermophilic niches and consequently prevalence of CRISPRs in thermophilic genomes is known (Anderson et al., 2011; Weinberger et al., 2012b). We analyzed the CRISPR loci from all Thermus genomes in this study in order to investigate the probable viruses infecting the genus and to uncover the organization of cas genes. Our examination divulged the ubiquity of CRISPR arrays within the genus Thermus, reflecting a resilient viral defense system. The predominance of CRISPRs among the Thermus group suggests presence and activity of phages in thermophilic environments. A wide scenario of phage invasion in T. oshimai was particularly denoted by high number of spacers present in this species. The uneven distribution of CRISPR arrays within a group can be explained by the hypothesis that probably the cost of harboring CRISPR elements in particular bacteria outweighs the benefits harnessed (Weinberger et al., 2012a). Cas genes, which are momentous to the functioning of CRISPRs, are divided into two classes (class I and II), depending on their mechanism of action. Class I CRISPR-Cas systems act by employing a number of Cas proteins which bring about the required action in a cascade of events. Class II systems, on the other hand, rely on single effector proteins for binding and cleavage of the target site. Based on the specific proteins involved, class I is further divided into types I, III and IV and class II is divided into types II and V. Widespread presence of type I (15 genomes) and type III (12 genomes) was observed in Thermus. Type IV systems were detected only in 3 genomes. Type IV systems have been known to be rare (<2%) in both bacteria and archaea (Makarova et al., 2015). Evolutionary relationships were uncovered on the basis of repeat sequences by clustering repeat sequences into conserved sequence families and structural motifs (Lange et al., 2013). CRISPR DRs transcribe repeat RNA sequences which serve as Cas protein binding templates. Repeat sequences show significant conservation in their sequence as well as hairpin structure forming motifs (Lange et al., 2013). Sequence conservation has been used as a criterion for clustering of DRs into families. Similarly, structural motif grouping is based on RNA loop structures. Using this analysis, we were able to identify the distribution of sequence families and structural families within this genus. Interestingly, motifs 20 and 31, which were identified in T. filiformis and HB8, have been known to constitute a mixture of both bacterial and archaeal domains. Viral diversity analysis delineated most probable viral predators for the Thermus group. A high proportion of spacers, however, were left un-assigned, implying the huge viral genosphere that is yet unexplored.

Genes implicated in natural transformability of HB27 include PilA1, PilA2, PilA3, PilA4, ComZ, CinA, DprA, ComEA, ComEC, PilF, PilC, PilM-N, PilN-O-W, ComF, PilQ, and PilD. Apart from HB27, all of the above genes were harbored by T. oshimai, CCB_US3_UF1, T. islandicus, JL-18, and SG0.5JP17-16. In this study, all Thermus genomes were found to show a homogeneous profile with respect to genes except PilA1-PilA4 and ComZ. Among competence genes, PilA gene plays a decisive role in efficient translocation (Schwarzenlander et al., 2009). PilA1-A4 genes and ComZ are present as a cluster in all genomes, however the arrangement and organization of genes show a difference across genotypes. The PilA-ComZ locus was found to be harbored in all Thermus genomes except T. igniterrae and T. antranikianii. Either these strains have not yet acquired the cluster or the cluster has been missed out during sequencing (draft genomes). The PilA-ComZ operon contains PilA1, A2, A3, A4, and ComZ as principal genes. The operon, however, also contains other genes with pilin/putative pilin/pseudopilin domains and genes coding for hypothetical proteins. The PilA-ComZ locus in Thermus was found to be comprised of 9–12 genes out of which 4–7 genes were pseudopilins with significant similarity to the PilA genes of strain HB27 (Figure 8). Apart from pseudopilin genes, the locus is comprised of genes coding for hypothetical proteins, chromosome segregation protein, apolipoprotein and DUF820 superfamily nucleases. PilA1 gene (present in 14 genomes) was found to be duplicated in 7 genomes. The two copies of PilA1 were, however, not identical and showed identity ranging from 67 to 81%. By convention, we have named the PilA1 of HB27 (471 bp) as copy 1. Copy 2 of PilA1 (present in 7 strains) is a smaller gene ranging from 272 to 373 bp (mean = 329 bp). In T. caliditerrae, only a fragmented copy (135 bp) of PilA1 was recovered. PilA2 (582–614 bp) was recovered in 7 strains as a complete copy. In case of T. scotoductus, identity at the C- and N- terminals and non-identity at the middle of the gene sequence was observed. A truncated version of PilA2 was identified in T. tengchongensis. PilA3 was recovered in a complete state (685–746 bp) in nine strains. In T. scotoductus, a truncated PilA3 was detected. The most deviant forms among PilA genes were observed for PilA4. Patterns of alignment among all PilA genes best denote this observation (Supplementary Figure 5). The alignment features are most consistent (>70%) in case of PilA1 through PilA3 genes. In case of PilA4, similarity among genes is observed only at the C- and N- termini of the genes. PilA4 of HB27 is 396 bp in size, however, in the other 13 strains harboring PilA4, only a small part (87–121 bp) showed significant similarity (>80%, e-value: 1E-15) to PilA4 of HB27. Additionally, truncated versions of PilA4 were observed in T. aquaticus, T. caliditerrae, and T. amyloliquefaciens.

A complete PilA-ComZ locus was identified in HB27, T. oshimai, CCB_US3_UF1, JL-18, SG0.5JP17-16, and T. islandicus. However, in other genomes, we could identify the presence of “genetic scars.” Genetic scars are truncated genes or pseudogenes (without start/stop codons), but show significant identity to functional PilA gene regions of HB27 (taken as a reference for all comparisons here). Two genetic scars were identified in T. scotoductus (A3 and A2), two in T. caliditerrae (A1 and A4), one in T. amyloliquefaciens (A4), one in T. aquaticus (A4) and one in T. tengchongensis (A2). Interestingly, an IS4 family transposase, which is known to be found widely across the Thermus genomes was found incorporated in the PilA-ComZ locus of T. aquaticus. This can be an evidence for recent transposition activity at the locus. Another evidence is the presence of a conjugative protein in the PilA-ComZ locus of HB27. Interestingly, a toxin-antitoxin (TA) gene system was found to be incorporated into the competence locus of T. caliditerrae and showed similarity to a TA pair in Rhodothermus marinus. Horizontally transferred regions are generally found to be associated with tRNA genes (Darmon and Leach, 2014). As expected, tRNA genes were located in close proximity of the PilA-ComZ loci of 10 genomes in this study. Out of the ten genomes in which tRNA genes were associated with the competence locus, seven are still not fully formed, thus giving a strong boost to the hypothesis of recent horizontal acquisition of the locus. Genetic relatedness of Thermus strains on the basis of PilA1-A4 genes was reconstructed to infer the evolutionary history of PilA genes on the basis of their sequence development (Supplementary Figure 6). PilA1 copy1 and copy2 gene trees highlighted strain CCB_US3_UF1 and JL-18 as an outlier for both trees and T. parvatiensis as an additional outlier for PilA1 copy1 gene. CCB_US3_UF1 and JL-18 retained their outlier positions in all other trees. In case of PilA2 and PilA3 sequences, T. scotoductus emerged as an outlier. For PilA4, T. parvatiensis, CCB, T. scotoductus and JL-18 lie on the out branches of the dendrogram.

Genes responsible for imparting competence in the genus Thermus can be divided into two groups: the first group of genes are responsible for uptake of DNA and the second group is responsible for transport of DNA into the cell. Some of these genes belong to the T4P family of proteins which form a complex on the membrane of the cell, spanning the S-layer, Outer Membrane (OM), Secondary Cell Wall Polymers (SCWP) and Peptidoglycans (PG) which comprise the periplasm and the inner membrane (IM). Sixteen genes have been known to play a role in natural transformation in HB27, the transformation machinery of which is the most extensively studied. PilQ is a secretin which forms a macromolecular homopolymeric complex on the outer membrane (Burkhardt et al., 2011) and binds to the incoming DNA. Once DNA reaches the periplasmic space, the pseudopilus proteins (pilA1 through pilA4) come into play, forming subunits of the DNA translocator complex present in the periplasm. PilA1-PilA4 pilins form a shaft like structure in the periplasm which is attached to the outer membrane through PilQ (Burkhardt et al., 2011) and to the inner membrane through a motor ATPase PilF. PilF uses the energy of ATP hydrolysis to draw DNA toward the inner membrane (Rose et al., 2011; Collins et al., 2013). Double stranded DNA further binds to the inner membrane protein ComEA. An inner membrane channel protein, ComEC (co-transcribed along with ComEA) has recently been shown to modulate the expression of PilA4 and PilN in relation to environmental cues like nutrient limitation and low temperature (Salzer et al., 2016). The DNA passing through the IM is single stranded DNA. Hence, a nuclease has to be acting on double stranded DNA to make it single stranded. This nuclease has not been identified yet.

A stable arrangement of the PilA-ComZ locus could be observed in T. oshimai, CCB_US3_UF1, JL-18, SG0.5JP17-16, and T. islandicus. The PilA-ComZ locus represents a small horizontally acquired region (genomic islet) on the chromosome. It can be distinguished by high number of hypothetical protein coding genes, gene duplications and their association with tRNA genes (Darmon and Leach, 2014). Gene duplication events observed in case of PilA1 genes demonstrate species radiation forces and amenability of the genomes to evolutionary forces (Roth et al., 2007). The presence of a transposase in T. aquaticus, a conjugative element in HB27, TA element in T. caliditerrae and truncated PilA genes (genetic scars) on the locus indicate recent horizontal origins. Truncated PilA4 genes observed across this genus denote that it has either not developed or has undergone degradation, determined by the presence of genetic scars which show similarity to N-/C-terminals of PilA4 of HB27, but not to the complete gene. The close association of a TA system with PilA-ComZ cluster of T. caliditerrae reflects recent acquisition of this cluster. TA systems have been found to be associated with genomic islands and other mobile elements. Association with a TA system, promotes the maintenance of a horizontally acquired island and stabilization into the host genome (Rowe-Magnus et al., 2003; Iqbal et al., 2015). Generally, horizontally transferred loci are marked with pseudogenes as there is a strong selection pressure against these regions (Hao and Golding, 2010). Pseudogenes are non-functional versions of a previously functional gene, which are in the process of getting lost from the genome (Hao and Golding, 2010). Pseudogenes are known to be associated with recently laterally acquired regions or failed HGT events (Hao and Golding, 2010). The predominance of pseudogenes can be due to the high rates of gene turnover in laterally acquired regions. Some of the truncated genes may get stabilized and some eliminated in due course of time. The evidence thus provided leads us to believe that the competence machinery in Thermus is of horizontal origin and in course of evolution, may get stabilized or eliminated. The presence of a highly efficient transformation system however does not ensure the incorporation of incoming DNA into the host genetic material. Natural transformation in native conditions is activated during environmental challenges such as starvation, wherein DNA is taken up from the environment (Seitz and Blokesch, 2013). Most of the nucleic acid taken up during this process is used to fulfill nutritional requirements. In this process, some amount of DNA may get incorporated into the host genome, thus diversifying the host pan repertoire and expanding the already diverse arsenal of Thermus group.

The pilus structure in Thermus imparted by the T4P genes plays a role not only in twitching motility and natural competence, but also in bacteriophage infection. PilA mutants have been shown to lose not only twitching motility and natural competence, but are also resistant to phage infection in HB27 and HB8 strains (Tamakoshi et al., 2011). On the basis of our comparative analysis for the Thermus group, we propose a link between pilus gene diversification and CRISPR abundance. Our data suggests continued acquisition and evolution of pilus gene structure among the analyzed Thermus genomes. Along with this, a CRISPR system with high number of spacers suggests a robust immune machinery against bacteriophages. In case of T. islandicus, only one questionable CRISPR array was observed, and no CRISPR array was observed in the case of T. antranikianii. Interestingly, a PilA-ComZ locus was also absent in T. antranikianii. The two observations, when coupled together suggest that T. antranikianii might be resistant to a huge proportion of viruses due to lack of pilus system genes which are implicated in phage entry, thus avoiding the need for harboring the CRISPR system. Other Thermus members have, on the other hand, chosen pilus mediated natural transformation as an important evolutionary trait, even though it makes them more susceptible to phage attacks, leading them to harbor more frequent CRISPR defense systems. Thus, natural transformation may be regarded as an overall benefit imparting trait in the small thermophilic genomes of Thermus. Natural transformation has played a role in survival of these organisms since long. Therefore, the dispersal of this system is a rather favorable phenomenon wide this genus even though it imposes an additional cost of harboring CRISPR machinery on them.