Table 1 provides information about the genomes.

Predicted protein sequences corresponding to all Acidithiobacilli proteomes were sorted using an all-vs.-all BLASTP script based on Best Bidirectional BLAST Hit (BBBH; Altschul et al., 1997) with an E-value of 1e-5. Protein families were constructed based on 50% of identity and 50% of coverage in the alignments (Altschul et al., 1997), assigning each protein to one protein family. The families with predicted proteins shared by all strains were selected and denominated the core-genome (Williams and Kelly, 2013; Hudson et al., 2014). The Acidithiobacillus core-genome was compared using BLASTP version 2.2.26 (Altschul et al., 1997) against NCBI non-redundant (NR) database in August of 2015, using a minimal E-value of 1e-5. Core families with exclusive similarity with Acidithiobacillus members, and not associated with any other microorganism, were selected and denominated unique (orphan) core genes. The selected unique protein families were checked manually using BLASTP, Psi-BLAST (Altschul et al., 1997) and HMMer version 3.0 (Eddy, 1998) against NR database with an E-value of 1e-4 to confirm their exclusive association with the Acidithiobacillus genus. The locus tags of the respective genes are provided in Table 2.

Collinear blocks between the genomes and conservation of gene neighbors were determined by MAUVE (Darling et al., 2010), RAST (Aziz et al., 2008; Overbeek et al., 2014; Markowitz et al., 2014a) and IMG-JGI (Markowitz et al., 2014b; Dhillon et al., 2015). Genomic contexts were visualized using Artemis of Sanger (Brettin et al., 2015).

IslandViewer (Rutherford et al., 2000) was used to predict genomic islands.

Protein coding sequences were annotated using an integrated pipeline consisting of BLASTP (Altschul et al., 1997) searches against NR database of NCBI with an E-value cutoff of 1e-3, Pfam (Punta et al., 2012), TigrFAM (Consortium, 2014), and Uniprot (Hofmann and Stoffel, 1993) database comparisons. Transmembrane regions in protein sequences were predicted with TMHMM (Haft et al., 2003) and TMPRED (Krogh et al., 2001). Computation of isoelectric point and molecular weight were made with ExPASy web tool (Bjellqvist et al., 1993; Nakai and Horton, 1999; Gasteiger et al., 2005).

Synonymous and non-synonymous substitution rates were calculated as follows: amino acid alignments of unique (orphan) core genes were constructed using MUSCLE (Edgar, 2004), and used as input for PAL2NAL (Suyama et al., 2006) with the nucleotide sequences to create the codon alignments of gene core families. The ratio of non-synonymous (K_a) to synonymous (K_s) nucleotide substitution rates (K_a/K_s ratios) were calculated using SeqinR package of R project (Charif and Lobry, 2007). Mean K_a/K_s ratios were assigned for individual unique (orphan) core genes (families I–V) by averaging all pairwise ratios within each family.

A combination of computational prediction tools PSORTb (Nakai and Horton, 1999; Yu et al., 2010), CELLO (Yu et al., 2006) and ProtCompB¹ (Yu et al., 2004) were used to perform whole genome analysis of unique core protein subcellular localization via the Sec Mechanism and Tat signal prediction (Natale et al., 2008; Bagos et al., 2010). The results derived from three prediction algorithms tools were combined according to majority to obtain a more accurate protein subcellular localization prediction.

Prediction of lipoproteins signals was made with LipoP Server (Juncker et al., 2003).

16S rRNA sequences from Acidithiobacillus genomes were identified by BLASTN-based script using an E-value threshold of 1e-5 and the databases GREENGENES (DeSantis et al., 2006), RDP (Cole et al., 2009) and SILVA (Pruesse et al., 2007) and were aligned using MAFFT (Katoh et al., 2002, 2005) alignment tool with L-INS strategy. Phylogenetic trees were constructed with MrBayes (Huelsenbeck and Ronquist, 2001; Ronquist and Huelsenbeck, 2003) and PHYML (Guindon et al., 2010), using the substitution model predicted for jModelTest2 (Guindon and Gascuel, 2003; Darriba et al., 2012).

The genes encoding families I–V were mapped onto the genomes A. ferrooxidans ATCC 23270, A. ferrivorans SS3, A. caldus ATCC 51756, and A. caldus SM-1 using DNAplotter (Carver et al., 2009). The origin of replication (Ori) of each genome was predicted between dnaN and dnaA as previously described (Valdés et al., 2008) and was used as the zero coordinate to orient the genome maps.

Metagenomic and metatranscriptomic sequences were downloaded from NCBI, JGI (Nordberg et al., 2014), and MG-RAST (Meyer et al., 2008; additional information can be found in Table 4) and were interrogated by BLASTX (Altschul et al., 1997) against the five core protein families with an E-value cut-off of 1e-5. The percent identity and coverage of sequences were analyzed for each alignment.

Figure 1 summarizes the bioinformatics pipeline used to recover five families of proteins and their corresponding genes that are taxonomically restricted to the genus Acidithiobacillus. Using a relaxed cutoff (1e-5) in a BLAST search, they were not detected in the NCBI nr database of more than 90 million proteins that includes the predicted proteins of Thermithiobacillus tepidarius, the nearest extant relative of the Acidithiobacilli.

Since Acidithiobacilli-specific protein families have almost no similarity with known proteins for other non-Acidithiobacilli representatives, we used a collection of bioinformatics resources in order to gain insights into potential protein functions based on hydrophobicity profiles, secondary structure predictions, predicted protein cell localizations and the comparison of consensus and profile sequences to pattern and domain databases (see Section Materials and Methods). Protein function predictions of the five Acidithiobacilli-specific protein families were examined using an analysis of their genomic contexts. Their differential expression was linked to previously published proteomic data derived from cells subjected to changes of pH, which is known to be a major selective pressure for members of the Acidithiobacillus genus (Baker-Austin and Dopson, 2007; see Table 3).

Figure 2 provides an example of the predicted protein properties deduced with bioinformatics tools and comparative genomic analysis for members of family II. Additional information for all five families I–V can be found in Supplemental Files 1, 2. In silico predictions demonstrate the power of integrative genomics approaches to gain insights into gene function. A significant prediction was made for an integral membrane segment with a moderate conservation profile within the family II. From the non-membrane associated portion of the protein, profile sequences were generated that have similarity to a pattern present in periplasmic binding proteins (Dwyer and Hellinga, 2004) and also solute carrier organic anion transporter family member 4A1 (Pizzagalli et al., 2003).

Comparative genome organization data demonstrated that there is conservation of gene neighborhood profiles that include genes predicted for cell division, surface proteins and ABC transport systems (Figure 2 and Supplemental File 3). Table 2 shows a detailed overview of the predicted properties based on amino acid sequences for families I–V.

Information regarding the expression of the genes encoding the five families was extracted from the literature and is presented in Table 3. RNA transcript analysis indicates that all five family genes are expressed in A. ferrivorans SS3 in two different conditions: continuous culture at 20°C (Christel et al., 2016a) and at 8°C (Christel et al., 2016b), adjusted to pH 2.5 with sulfuric acid plus trace elements. A proteomic study of A. ferrooxidans ATCC 23270 on elemental sulfur as electron donor under aerobic and anaerobic conditions (Osorio et al., 2013) showed that family III was expressed in this strain. A proteomic study of A. caldus ATCC 51756 using cells grown at pH 2.5 (optimum growth pH) vs. pH 1 and 4, demonstrated up-regulation of core families I, III, and IV when cells were shifted from pH 2.5 to 1 and that family V was upregulated when cells were shifted from pH 2.5 to 4 (Table 3; Mangold et al., 2013). These data show that the genes for the five families (i) are expressed and thus are unlikely to be mis-annotated open reading frames with no coding capacity and (ii) provide evidence that families I, III, IV, and V could be involved in responses to acid stress at least in A. caldus. It remains to be determined if changes in RNA levels are associated with these genes in the other Acidithiobacilli.

In addition, RNA transcripts in metatranscriptomes of the Dabaoshan and Yunfu Pond mines in China (Chen et al., 2015) were detected that exhibited sequences similar to families I–V (Table 3, right hand side) from A. ferrooxidans, A. ferrivorans, and A. thiooxidans, although no strain specificity could be determined. This supports the idea that the five families are bona fide genes.

In order to make a comprehensive summary of the potential gene function inferred from all the evidence presented, a schematic summary is presented in Figure 3. Families I, II, and V, have predicted transmembrane segments that, in conjunction with protein sorting signal identification, provide preliminary information about their cellular location. Profile and consensus sequences comparisons against public databases only provided information about family II. Family II sequences have motifs similar to those of periplasmic binding proteins, usually associated with ABC transport for the substrate specific incorporation of nutrients and scarce molecules or beneficial solutes under extreme environmental conditions (Cuneo et al., 2008). We suggest that members of family II could be distant relatives of periplasmic binding proteins whose specific substrate(s) and functional role remains to be investigated.

Families III and IV have predicted protein localizations associated with inner membrane and periplasmic spaces and their strong lipoprotein signatures, in addition to genomic context information, provide clues for their potential role in key physiological processes, such as lipid metabolism. We hypothesize a potential connection between membrane associated lipoproteins, lipid metabolism and membrane stability as a requirement for low pH lifestyle (Baker-Austin and Dopson, 2007; Liljeqvist et al., 2015).

Predicted protein properties of all families I–V, suggest a general involvement in functions associated with membrane processes perhaps involving roles in membrane stability, transport processes, and/or the generation of molecular components to allow the synthesis and incorporation of hydrophobic molecules into the membrane increasing its stability in low pH.

It has been observed in many bacteria that the gene order relative to OriC is highly conserved along the chromosomal replicores (Sobetzko et al., 2012). Also, essential and highly expressed genes tend to be encoded close to oriC (Rocha, 2004). This heightened activity can be attributed to gene dosage effects during chromosome replication especially in rapidly dividing cells, but underlying physical properties of the circular chromosome, including an inferred gradient of DNA superhelical density from the origin to the terminus, are also known to be involved in influencing gene expression (Sobetzko et al., 2012).

In particular, it has been observed that several genes involved in acid stress, including envelope remodeling, are located close to oriC in the gammaproteobacterium Dickeya dadantii (Jiang et al., 2015). Given the possibility that genes of families I–V could be involved in acid stress response and that this response might be associated with chromosome topology, we determine their chromosomal locations on the closed circular chromosomes of A. ferrooxidans ATCC 23270^T, A. ferrivorans SS3^T, A. caldus ATCC 51756^T, and A. caldus SM-1 using DNAplotter (Carver et al., 2009; Figure 4). In all these chromosomes, the five family genes exhibit a tendency to be located nearer Ori rather than the terminus, especially in the cases of A. ferrooxidans and A. ferrivorans. In the latter two chromosomes, the gene order relative to Ori is conserved but is inverted, perhaps due to inter-replicore translocation that is known to be common around Ori in other microorganisms (Eisen et al., 2000; Khedkar and Seshasayee, 2016). Three of the families have genes predicted to DNA handling functions in their gene neighborhoods ordered in tightly clustered associations that could be operons; for example, rmuC (DNA recombination) near family IV, and dnaB and radA (DNA helicase and DNA repair, respectively) near family V. These genes are usually associated with DNA replication and cell division (Figure 4). The juxtaposition of ftsL, an essential cell division protein (Guzman et al., 1992), to the gene encoding family II and its closeness to the family III gene (Figure 4) strongly suggests that family II and III are involved in cell division perhaps through cell envelope remodeling. Their proximity to Ori could enhance the ability of the Acidithiobacilli to respond to changes in environmental acidity at early stages of cell division. Such changes might be more difficult to accomplish during later stages of cell division or at the resting stage.

Taxonomically restricted genes, such as families I–V, are referred to as orphans genes or ORFans (orphan open reading frames; Fischer and Eisenberg, 1999; Pedroso et al., 2008; Tautz and Domazet-Loso, 2011). ORFans can be artifacts of annotation, non-coding RNA genes or protein encoding genes (Prabh and Rodelsperger, 2016). In the case of families I–V, there is evidence that those from A. caldus encode proteins and that families I–V from A. ferrooxidans, A. ferrivorans, and A. thiooxidans express RNA (Table 3). Given the highly conserved sequences similarity between the respective families from the different Acidithiobacillus species, it is reasonable to suggest that all are protein coding genes, as observed for the A. caldus families and are not “merely” RNA genes. However, in order to provide additional evidence for protein coding capacity, selection pressure was measured as the ratio of the synonymous and non-synonymous rates of amino acid substitution (dN/dS), also called omega (ω) for all families. The omega values for families I–V are 0.07, 0.05, 0.03, 0.05, and 0.08 respectively. An ω <1 can be interpreted as evidence for negative selection and most likely such a sequence would correspond to a protein encoding gene (Prabh and Rodelsperger, 2016). The omega values are considerably <1 for all five families providing compelling evidence that they are protein-encoding genes.

The genes encoding families I–V are not found in T. tepidarius that subtends the genus Acidithiobacillus and shares the last common ancestor with it, nor are they found in any other organism that has sequence information in the NCBI nr database. So questions arise as to the origin and evolution of the five families.

We propose three main hypotheses.

Although a lack of definitive evidence leaves all three hypothesis unimpaired, we speculate that the emergence of families I–V could have helped promote by whatever means (direct activity of the encoded proteins, or via sensing or regulatory mechanisms) the ability of the last common ancestor of the Acidithiobacillus genus and T. tepidarius to transition from a neutral pH environment to one that was increasingly acidic and finally to one that was extremely acidic. In this scenario, the transition process could have provided opportunities for the Acidithiobacillus genus to diverge from the T. tepidarius lineage. This hypothesis requires additional evidence, especially experimental evidence, to clearly pinpoint the specific functions and physiological roles of the five families.

In order to evaluate the sensitivity and specificity of the families to discriminate between Acidithiobacillus species, the DNA sequences of families I–V were concatenated for each Acidithiobacillus species and compared by BLASTN against each Acidithiobacillus species. The results are reported as % nucleotide identity between the concatenated probe and each Acidithiobacillus species (Figure 5). The dark blue diagonal indicates high nucleotide identity, as expected, between the concatenated probe and its respective sequences in the corresponding genome. Importantly, the concatenated probes from one species have lower levels of sequence identity when compared to other species. For example, the concatenated probe from A. caldus has only 69% identity (white cell) when compared to sequences present in the genome of A. ferrivorans.

These data indicate that the concatenated families are capable of discriminating between the different Acidithiobacilli species used to build the concatenated probes, but are they capable of phylotyping new genomes that did not contribute to building the probes?

During the course of this investigation four new genomes of A. ferrooxidans (strains BY0502, DLC-5, YQH-1, and Hel18), one A. caldus genome (strain MTH-04) and six genomes of A. thiooxidans were released (Table 1), providing an opportunity to test the discriminatory powers of the family probes on new genomes.

First, the concatenated family probes, described in the previous experiment, were used in BLASTN comparisons with the new genomes. The results are reported as % nucleotide identity between the concatenated probe and each Acidithiobacillus species (leftmost four columns, Figure 6). The concatenated probes clearly have the ability to discriminate between A. caldus MTH-04, A. thiooxidans DMC, A. ferrooxidans BY0502, A. ferrooxidans YQH-1, and A. ferrooxidans Hel18, indicated by the dark blue color (close to 100% sequence identity). However, there is one anomalous identification. A. ferrooxidans BY0502 exhibits the best match with the A. ferrivorans concatenated probe (bottom row), suggesting that this species might not be A. ferrooxidans.

In order to determine if this anomaly could be attributed to one (or more) of the families in particular, the experiment was repeated with each individual family (Figure 6). Each family correctly identified the new genomes of A. ferrooxidans, A. thiooxidans and A. caldus with the exception of A. ferrooxidans BY0502. The highest percentage matches of all five families to A. ferrooxidans BY0502 were to the probes built from A. ferrivorans, confirming the results using the concatenated family probe.

Because of the vexing problem of the anomalous A. ferrooxidans BY0502 in which the family I–V probes place it closer to A. ferrivorans than A. ferrooxidans, it was decided to use other approaches to investigate its phylogeny using ANI (Goris et al., 2007) and TETRA (Richter and Rosselló-Móra, 2009). Both approaches indicate that A. ferrooxidans BY0502 is not related to A. ferrooxidans because of the low values of ANI and TETRA, 83.4 and 0.988, respectively, between the two genomes. Nor can it be classified in the A. ferrivorans clade, with low values of 91.7/0.996 (ANI/TETRA values), although it is more closely related to A. ferrivorans than A. ferrooxidans. In order to investigate further the phylogeny of A. ferrooxidans BY0502, 16S rRNA sequence analysis was carried out that placed it in a clade with A. ferriphilus, subtended by the clade A. ferrivorans with a bayesian posterior probability node support of 1 that strongly endorses the proposed phylogeny (Figure 7). Therefore, we suggest that A. ferrooxidans BY0502 is more likely to be an A. ferriphilus-like microorganism; an hypothesis that requires confirmation using other phylogenetic approaches. This example demonstrates the power of the family probes to aid in the identification of the Acidithiobacillus genus with discriminatory powers to suggest species at least for those under interrogation in the present study.

Gaining insight into the structure, organization, and function of microbial communities (microbiomes) has been proposed as one of the major research challenges of the current decade (2020 visions, 2010) and metagenomic and metatranscriptomic approaches present major opportunities for advancing our knowledge in this area. One of the most promising areas of metagenomics research is the use of shotgun methods to sequence random fragments of DNA (or RNA) in an environmental sample. This information can then be analyzed for microbial diversity, prediction of gene functions and biochemical pathway model building. Many bioinformatic approaches have been developed to handle the typically enormous amounts of data generated by metagenomics investigations (e.g., reviewed in Hiraoka et al., 2016).

One of the most straightforward and computationally less demanding approaches to estimate microbial diversity in a microbiome is the use of marker genes (molecular probes; Wu and Eisen, 2008; Liu et al., 2011; Wu and Scott, 2012; Kim et al., 2013; Darling et al., 2014). For example, rRNA sequences from known organisms can be used to computationally search the shotgun sequences for similar sequences or can be coupled with rRNA-PCR to pull out and extend specific sequences. These methods provide an overview of the phylogenetic distribution (phylotyping) of the cell-based life present in a sample but they have their limitations (reviewed in Fabrice and Didier, 2009).

Taxonomically restricted protein encoding genes have been used for phylotyping, including the recombinase A gene family and the RNA polymerase beta subunit (Wu et al., 2011), genes specifically targeting the Acidithiobacilli (Nieto et al., 2009; Nuñez et al., 2014, 2016) and many other examples (Liu et al., 2011; Segata et al., 2011; Wu et al., 2013; Darling et al., 2014). However, such marker genes are subject to HGT and evolutionary rate differences that can exacerbate the interpretation of phylogenies. Since the five families are taxonomically restricted to the Acidithiobacilli and do not appear to be prone to HGT, we decided to examine their ability to identify the Acidithiobacillus genus and to discriminate between different species of the Acidithiobacilli (Figures 6, 7) in environmental metagenomic and metatranscriptomic samples. For the first objective, the amino acid sequence of all five families from all participating Acidithiobacillus species (A. ferrooxidans, A. ferrivorans, A. thiooxidans, and A. caldus) was concatenated (five families × nine species). This was considered as a general probe for the Acidithiobacillus genus (genus-level probe). A second series of probes was constructed where the protein sequences of the five families was concatenated according to species, generating five different probes each one specific for an Acidithiobacillus species (e.g., A. ferrooxidans probe = the concatenation of families I–V of A. ferrooxidans). These probes were then used in a BLASTX searches to interrogate several environmental metagenomes and metatranscriptomes listed in Table 4.

The metagenomes were chosen to include low pH environments such as mining operations and AMD, where Acidithiobacilli have previously been reported, and also environments of intermediate acidity (e.g., Black Smokers, Tui Malila), neutral pH (e.g., Hydrothermal vent, Guaymas Basin), and high pH (e.g., Marine Microbial Communities, Lost City) where Acidithiobacilli have not been detected. Two low pH metatranscriptomes were also included in the analysis. The results of the BLASTX interrogations are shown in Figure 8 and the results are summarized in Table 4.

Inspection of the left hand column of Figure 8 indicates that the genus-level probe detects sequence similarity in all the samples except for the Pink Biofilm from the Richmond mine. This is in agreement with the report that no Acidithiobacilli were detected in the Pink Biofilm but were detected in all the other samples (references provided in Table 4). The absence of Acidithiobacilli in the Pink Biofilm sample could be due to its extremely low pH (pH 0.83) which is thought to be too acidic to support their growth (Tyson et al., 2004). In addition no Acidithiobacilli were detected in samples from the Black Soud Mine, Black Smokers (Tui Malila), Hydrothermal Vent (Guaymas Basin), Marine Microbial Communities (Loihi), Deep Ocean Microbial Communities (Juan de Fuca), Marine Microbial Communities (Lost City), which is also in agreement with the published literature (references found in Table 4). The conclusion is that the Acidithiobacilli genus-level probe appears to have good specificity and sensitivity in detecting Acidithiobacilli in environmental metagenomes but more samples are required to develop statistical support for this assertion.

Table 4 also indicates that the families can be used to interrogate metatranscriptomes and provides additional evidence that the genes of family I–V are transcribed. This evidence was used to construct the right hand column presented earlier in Table 3.

However, caution is required in the interpretation of the use of the species-specific probes. In case A (see Figure 8), both the sequence identity (100%) and sequence coverage (83.5–95.1%) of the A. ferrooxidans probes of families II and III strongly support the contention that sequences corresponding to them are present in the Carnoulès metagenome. However, in case B, although there is good coverage of the A. ferrooxidans family I and V probes (99.3–99.6%), the sequence identity is lower (80–83%). This suggests that these families probably belong to A. ferrooxidans in the metagenome but that they have diverged somewhat from the probe sequences. Recovery of such sequences would expand the number and diversity of such sequences that could be helpful for elucidating their function and shedding light on their evolution. In case C, both the coverage and identity are lower and the hits are to probes developed for A. thiooxidans and A. caldus family III and family IV. This suggests that the Carnoulès metagenome contains A. thiooxidans-like and A. caldus-like organisms that exhibit low sequence similarity to families III and IV, but not to the other families. As in case B, these sequences could be helpful for later studies to help unravel sequence function and evolution. A final case marked by asterisks in Figure 8 illustrates the common finding of sequence similarity to metagenomic reads that are truncated. Truncated sequences that have high similarity to the probes could potentially be extended by PCR using primers designed from the probes and subsequently analyzed.

With these caveats in mind, families I–V satisfy a number of criteria for use as identification markers for Acidithiobacilli in genomic, metagenomic/metatranscriptomic investigations. They are universally present in the genus, not present in other genera and are not subject to HGT. Preliminary evidence also points to association of at least three of the families (Families I, III, and IV) in envelope remodeling and lipid metabolism possibly associated with acid stress response and so could serve as PhyEco (for phylogenetic and phylogenetic ecology; Wu et al., 2013) markers for certain acidic environments including AMD and biomining operations.

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.