ABSTRACT
Members of the genus Acidithiobacillus, which can adapt to extremely high concentrations of heavy metals, are universally found at acid mine drainage (AMD) sites. Here, we performed a comparative genomic analysis of 37 strains within the genus Acidithiobacillus to answer the untouched questions as to the mechanisms and the evolutionary history of metal resistance genes in Acidithiobacillus spp. The results showed that the evolutionary history of metal resistance genes in Acidithiobacillus spp. involved a combination of gene gains and losses, horizontal gene transfer (HGT), and gene duplication. Phylogenetic analyses revealed that metal resistance genes in Acidithiobacillus spp. were acquired by early HGT events from species that shared habitats with Acidithiobacillus spp., such as Acidihalobacter, Thiobacillus, Acidiferrobacter, and Thiomonas species. Multicopper oxidase genes involved in copper detoxification were lost in iron-oxidizing Acidithiobacillus ferridurans, Acidithiobacillus ferrivorans, and Acidithiobacillus ferrooxidans and were replaced by rusticyanin genes during evolution. In addition, widespread purifying selection and the predicted high expression levels emphasized the indispensable roles of metal resistance genes in the ability of Acidithiobacillus spp. to adapt to harsh environments. Altogether, the results suggested that Acidithiobacillus spp. recruited and consolidated additional novel functionalities during the adaption to challenging environments via HGT, gene duplication, and purifying selection. This study sheds light on the distribution, organization, functionality, and complex evolutionary history of metal resistance genes in Acidithiobacillus spp.
IMPORTANCE Horizontal gene transfer (HGT), natural selection, and gene duplication are three main engines that drive the adaptive evolution of microbial genomes. Previous studies indicated that HGT was a main adaptive mechanism in acidophiles to cope with heavy-metal-rich environments. However, evidences of HGT in Acidithiobacillus species in response to challenging metal-rich environments and the mechanisms addressing how metal resistance genes originated and evolved in Acidithiobacillus are still lacking. The findings of this study revealed a fascinating phenomenon of putative cross-phylum HGT, suggesting that Acidithiobacillus spp. recruited and consolidated additional novel functionalities during the adaption to challenging environments via HGT, gene duplication, and purifying selection. Altogether, the insights gained in this study have improved our understanding of the metal resistance strategies of Acidithiobacillus spp.
INTRODUCTION
The genus Acidithiobacillus is composed of a group of obligatory acidophilic chemolithotrophic Gram-negative bacteria that conduct energy conservation by oxidizing reduced inorganic sulfur compounds (RISCs), such as thiosulfate (S2O32−), hydrogen sulfide (H2S), and elemental sulfur (S0). This attribute has been exploited in the bioleaching of ores and has also been the main cause of the accelerating formation of polluting acid mine drainage (AMD) water, which brings disruptive consequences for terrestrial and aquatic ecosystems (1–5). Currently, 8 species have been reported in the genus Acidithiobacillus, including Acidithiobacillus albertensis, Acidithiobacillus caldus, Acidithiobacillus cuprithermicus, Acidithiobacillus ferridurans, Acidithiobacillus ferriphilus, Acidithiobacillus ferrivorans, Acidithiobacillus ferrooxidans, and Acidithiobacillus thiooxidans. Acidithiobacillus species are universally found at AMD sites where highly toxic inorganic forms of heavy metals and metalloids such as arsenic, mercury, copper, cadmium, and cobalt often dominate (1, 6–9). As revealed in previous studies, strains of A. ferrooxidans can adapt up to 800 mM Cu2+, 84 mM As3+, 1,071 mM Zn2+, 500 mM Cd2+, 1,000 mM Ni2+, and 50 mM Hg2+ (10, 11), while strains of A. thiooxidans and A. ferrooxidans are able to adapt up to 120 and 100 mM As5+ (12). It is obvious that Acidithiobacillus species surviving in acidic metal-rich leaching environments must possess certain highly efficient and advanced metal resistance mechanisms, rendering them ideal research models to improve our understanding of metal resistance. There are four main strategies that acidophiles generally utilize to increase resistance to metals, including the chelation and precipitation of toxic metals, the efflux of the toxic metals out of the cells by molecular pumps, the enzymatic conversion of metal ions to less toxic forms, and/or enhanced cell walls and membranes that are impermeable to toxic metal ions (13). These strategies are closely related to the genetic properties of the bacteria. The genes conferring toxic metal resistance are prone to clustering together as gene operons in Acidithiobacillus genomes, such as the mer operons that confer mercury resistance (14). However, the mechanisms driving the clustering of genes into operons are far from clear. For example, the mer operon generally contains merA, merD, merR, merP, merC, and merT, although their organization may be different. The reported genes related to mercury, arsenate, cadmium, zinc, cobalt, and copper resistance and their functions in bacteria are summarized in Table 1.
Genes involved in toxic metal resistance and transformation
Horizontal gene transfer (HGT), natural selection, and gene duplication are three main engines that drive the adaptive evolution of microbial genomes, whereas their relative importance is still ambiguous (15, 16). HGT events generally occur via mobile genetic elements (MGEs) carrying functional genes that enhance the adaptability and versatility of the acceptors. Recent studies analyzing gene-transfer webs among a vast number of sequenced genomes revealed that more than 20% of the microbial genes may be acquired through HGT (17, 18). MGEs such as plasmids, phages, genomic islands (GIs), transposons, and insertion sequences are indicative of HGT events. GIs, which are large chromosomal regions usually found adjacent to or integrated within tRNA genes, represent a versatile gene pool (19). Deviant G+C contents, as well as incongruences between phylogenies of functional genes and species, can also be used as indicators of HGT events (20–23). Some previous studies have described the overall metal resistance mechanisms in acidophiles and indicated that HGT was a main adaptive mechanism by which the acidophiles cope with heavy-metal-rich environments (13, 24–27); for example, it was reported that A. ferrooxidans ATCC 53993 harbored GIs, inserted after a tRNA metabolism related gene (MiaB type), containing genes encoding mercury detoxification determinants (merA, merC, and merR) and copper translocating P-type ATPase (13).
However, direct evidence of HGT in other Acidithiobacillus species in response to challenging metal-rich environments and the mechanisms addressing how metal resistance genes originated and evolved in Acidithiobacillus are still lacking. The rapid development of next-generation sequencing technologies has enriched the public database with numerous sequenced genomes. Exploiting these resources, large-scale comparative genomic analyses are improving our understanding of extensive genetic diversity and evolutionary history. To date, 8 species have been reported in genus Acidithiobacillus, and strains of 6 species (A. thiooxidans, A. ferrooxidans, A. caldus, A. ferrivorans, A. albertensis, and A. ferridurans) have been sequenced. To reveal the genetic traits and phylogenetic history of metal resistance in genus Acidithiobacillus, a comparative genomic analysis of 37 strains within the genus Acidithiobacillus was carried out. In addition, molecular phylogenetic analysis, gene contents, molecular conservation, and linear analyses along with selective pressure and codon usage analyses were used to explore the evolution of toxic metal resistance genes in the genus Acidithiobacillus (28–31). Our results provide explicit evidence that HGT played a crucial role in driving the evolution of heavy metal resistance in Acidithiobacillus species surviving extreme environments. The insights gained in this study improve our understanding of the complex evolutionary history of toxic metal resistance genes in Acidithiobacillus spp. and its roles in the biogeochemical cycling of different metals.
RESULTS
Genomic features.A summary of the features of 38 available sequenced genomes of Acidithiobacillus spp. are listed in Table 2. The G+C contents of the 38 genomes range from 52.6% to 61.5%. These genomes vary in size by approximately 2.48 Mb (ranging from 1.70 to 4.18 Mb) with coding sequence (CDS) numbers ranging from 1,617 to 4,359, suggesting substantial strain-to-strain variations.
General features of bacterial genomes used in this study
Phylogenic analyses and pangenome analyses.To associate the distributions of metal resistance genes in Acidithiobacillus spp. with their phylogenetic affiliations, we constructed a phylogenetic tree of the Acidithiobacillus spp. and closely related species on the basis of 16S rRNA gene sequences using neighbor-joining (NJ) methods (see Fig. S1 in the supplemental material). The results showed that Acidithiobacillus species are phylogenetically close to Thermithiobacillus spp. Since no 16S rRNA gene sequence could be retrieved from the genome of Acidithiobacillus sp. strain NORP59, we constructed a phylogenetic tree of the 38 Acidithiobacillus genomes with their phylogenetic affiliations on the basis of whole-genome sequences (see Fig. S2). The results showed that Acidithiobacillus sp. NORP59 forms a clade with Thioglobus sp. strain MED-G25, distant from the standard strains of Acidithiobacillus spp. In addition, Acidithiobacillus sp. NORP59 showed low values of BLASTN-based average nucleotide identity (ANI; average 75%) and tetranucleotide composition regression (Tetra; average 0.758) in comparison with those of other available Acidithiobacillus strains (see Tables S1 and S2), indicating that strain NORP59 was not a member of genus Acidithiobacillus; thus, it was excluded from any further analyses.
To reveal the genomic features specific to each strain, we identified orthologous groups among the 37 Acidithiobacillus genomes using the Bacterial Pan Genome Analysis (BPGA) pipeline (32). Our analysis of the 37 genomes revealed a pangenome containing 141,937 putative protein-coding genes in the genus Acidithiobacillus. Of these 141,937 genes, 16,317 (11.5%) were clustered into the core genome of Acidithiobacillus spp., 95,704 (67.4%) were represented in the accessory genome, and 5,332 (3.8%) were identified as strain-specific genes. The numbers of specific genes ranged from 0 to 811 (Fig. 1A). The results of BacMet database annotation revealed that the accessory genome had a higher proportion (8.2% [7,846 genes]) of genes involved in heavy metal resistance activities (Fig. 1C) than the core genome (5.2% [851 genes]) and a higher proportion of unique genes (0.3% [16 genes]). We constructed a robust phylogenetic tree of the 37 genomes on the basis of the concatenation of the 110 single-copy core genes in each genome using the NJ method (see Fig. S3). The phylogenetic trees, inferred using ML and UPGMA methods (see Fig. S4 and S5), were congruent with the NJ phylogenetic tree, showing that the phylogenetic tree was highly reliable. A chronogram for NJ phylogenetic tree was constructed with a calibration point to explore the divergence times of Acidithiobacillus (Fig. 2). The most recent common ancestor (MRCA) of Acidithiobacillus (i.e., the emergence time of A. caldus) was estimated to have emerged around 800 million years ago (MYA) (Fig. 2), whereas the MRCA of A. thiooxidans, A. ferrooxidans, A. ferrivorans, and A. ferridurans was estimated to have emerged around 439 MYA (Fig. 2). These well-supported core gene trees showed that some strains in different reported species were grouped together, such as A. ferridurans JCM 18981 and A. ferrooxidans IO-2C, Acidithiobacillus sp. strain GGI-221 and strains of A. ferrooxidans, and A. albertensis DSM 14366 and strains of A. thiooxidans. These results suggested that there were mistakes in the classification of Acidithiobacillus spp. The ANI approach was further applied to evaluate the genomic similarity of microorganisms (33). The ANI results (Table S1) justified the conclusion from the phylogenetic analysis. All 37 strains could be classified into 9 species based on an ANI of ≥96%. A comparison of strains JCM 18981 and IO-2C resulted in a high ANI (98.9%), suggesting that they belong to the same species (A. ferridurans). A comparison between strains A. albertensis DSM 14366 and A. thiooxidans spp. resulted in a high ANI (97.3%), suggesting that they belong to the same species (A. thiooxidans). A comparison between strains GGI-221 and A. ferrooxidans spp. resulted in a high ANI (99.5%), suggesting that they belong to the same species (A. ferrooxidans). Thus, in this article, strain A. ferrooxidans IO-2C was renamed A. ferridurans IO-2C, strain A. albertensis DSM 14366 was renamed A. thiooxidans DSM 14366, and Acidithiobacillus sp. strain GGI-221 was renamed A. ferrooxidans GGI-221.
Pangenome analysis of strains in the genus Acidithiobacillus. (A) Petal diagram of the pangenome. Each strain is represented by a colored oval. The center is the number of orthologous coding sequences shared by all strains (i.e., the core genome). Numbers in nonoverlapping portions of each oval show the numbers of CDSs unique to each strain. The total numbers of protein-coding genes within each genome are listed below the strain names. (B) Mathematical modeling of the pangenome and core genome of Acidithiobacillus. (C) Proportions of genes involved in heavy metal resistance in unique genes, accessory genome, core genome, and pangenome according BacMet database annotation.
Chronogram of the 37 Acidithiobacillus strains. The value near each internal branch is the estimated emerge time for that branch. Nodes with fossil record corrections are indicated with a red star.
Model extrapolation of the pangenome of Acidithiobacillus.Vernikos et al. (34) stated that a mathematical extrapolation of the pangenome would be highly reliable provided that sufficient genomes (more than five) are involved. As shown by the deduced power law regression function (Ps [n] = 4451.42n0.36465) (Fig. 1B), the pangenome of Acidithiobacillus had a parameter (γ) of 0.365, falling into the range of 0 < γ < 1, suggesting that the pangenome is open (35). The deduced exponential regression (Fc [n] =2052.65e−0.0714494n) revealed that the extrapolated curve of the core genome followed a steep slope, reaching a minimum of 441 gene families after the 37th genome was added (Fig. 1B). The numbers of core genes were relatively constant, such that an extra genome added would not significantly affect the size of core genome, which was consistent with the notion that the core genes were conserved genes universally present in all strains (36).
Identification of putative mobile genetic elements.Mobile genetic elements (MGEs) are specific genome segments capable of intra- and extracellular movements, carrying gene(s) encoding putative functions, which are considered indicators of HGT events (37). In this study, MGEs including transposases, integrases, genomic islands (GIs), and phage-associated genes were identified and compared in all genomes. Each Acidithiobacillus genome in this study contained a number of transposons (see Table S3). The number of transposon copies per genome ranged from 143 (A. thiooxidans strain CLST) to 334 (A. caldus strain SM-1). Members of the IS3, IS5, IS21, and IS110 families were the most common. The numbers of GI copies per genome ranged from 25 (A. ferrivorans strain SS3) to 66 (A. thiooxidans CLST) (see Table S4). The numbers of prophages and/or prophage remnants ranged from 0 to 6 (see Table S5). The BacMet annotation results revealed that the genomes of Acidithiobacillus spp. contained a vast repertoire of metal resistance genes, with an average 160 BacMet hits in each genome. A visualization of the chromosomes of representative strains of different species of Acidithiobacillus (see Fig. S6) revealed that many of these genes were distributed within genomic island regions or flanking other mobile genetic elements such as transposons, indicating that putative HGT events contributed greatly to the adaptive survival of Acidithiobacillus spp. in toxic metal-rich niches.
Estimation of gene family gains and losses.We further investigated how the metal-resistance-related gene families might have been gained and lost during the evolution of A. thiooxidans, A. ferrooxidans, A. ferrivorans, A. ferridurans, and A. caldus (see Fig. S7). Gene family gains were observed on the branches from A. ferridurans, A. caldus, and the most recent common ancestor (MRCA) of A. ferridurans and A. ferrooxidans. In all, new gene families arose over time probably via HGT, with only a few gene families being lost.
Case study of Hg-related genes.The overall distributions and organization of Hg resistance genes in 37 Acidithiobacillus strains are summarized in Fig. S8. Hg resistance genes merA, merC, and merR were detected in most Acidithiobacillus genomes, while merD, merP, and merT were only detected in A. ferrooxidans, A. thiooxidans, and A. ferrivorans genomes. The gene merB that encodes organomercury lyase was missing in Acidithiobacillus genomes. The mer genes were prone to group together as mer clusters in Acidithiobacillus genomes, which could be categorized into six subgroups (merCAD, merRTPA, merTAR, merACR, merAR, and merA). Some strains in this study possessed two different mer clusters in their genomes (Fig. S8). However, the mer clusters were not detected in Acidithiobacillus sp. strain SH, possibly resulting from incomplete sequencing.
(i) Origin and evolution of mer clusters. The divergent distribution and organization of mer genes in Acidithiobacillus raised a question as to their origin and evolution. The deviant G+C contents can be used to detect HGT (20, 21). Herein, we detected the G+C contents of mer clusters and their corresponding genomes. The results showed that the G+C contents of the merCAD clusters were higher than those of the genomes in A. ferrooxidans strains (62.8% versus 57.6% to 61.0%, respectively); the G+C contents of respective merAPTR clusters were higher than those of the genomes in A. thiooxidans strains (60.5% to 61.1% versus 52.4% to 53.2%), A. ferrooxidans strains (62.6% to 62% versus 58.6%), and A. ferridurans strains (63.7% to 64.1% versus 58.4% to 61.3%). The merTAR and merACR gene clusters also showed variations of G+C contents between mer gene clusters and their corresponding genomes (see Fig. S9). To further elucidate the evolution of the mer gene clusters, we compared the chromosomal regions flanking the mer gene clusters among the 37 Acidithiobacillus strains and found that the genes in the upstream and downstream regions were conserved among strains of the same species except for merACR clusters (see Fig. S10). Several indicators of HGT, such as flanking tRNA genes and genes encoding mobile element proteins, were found near the merRTPA and merACR gene clusters. The results showed that strains of the same species generally shared identical insertion sites, whereas strains from different species had different insertion sites (Fig. S10), suggesting that mer clusters may be acquired more than once in Acidithiobacillus species. merACR gene clusters were found in different insertion sites regardless of whether the species were different or the same, suggesting that merACR clusters may be acquired via different HGT events (Fig. S10).
To elucidate the origin of mer genes clusters in Acidithiobacillus, a NJ phylogenetic tree was constructed on the basis of the MerA protein sequences. As shown in Fig. 3, the strains were separated into three group: group I included merCAD and merRTPA clusters, group II represented merACR clusters, and group III contained merRTAC, merTA, and merAR clusters and stand-alone merA genes. The phylogenetic tree revealed that the group I and mer genes of Burkholderiales members were sister groups, while those in group II were more similar to mer genes from Acidihalobacter spp., and those in group III were more similar to mer genes of Thiobacillus spp. and Acidiferrobacter thiooxydans, implying that the group I mer cluster may have been acquired via HGT from members of Burkholderiales, the group II mer cluster from Acidihalobacter spp., and the group III mer cluster from Thiobacillus spp. or Acidiferrobacter spp. in early evolutionary history. In addition, A. thiooxidans Licanantay containing cluster merRTAC branches near the base of the merTA and merAR clade. We inferred that among group III, gene clusters merRTAC, merTA, and merAR likely originated from A. thiooxidans, followed by gene loss in A. thiooxidans and HGT to A. ferrivorans. Besides, only stand-alone merA genes were detected in strains of A. caldus. The phylogenetic tree revealed that merA1 and merA2 genes of A. caldus form separate groups, in which merA2 genes of A. caldus were more similar to those from A. ferrivorans. The genes upstream and downstream were conserved for merA1 and merA2 genes, respectively (Fig. S10). Several genes encoding mobile element proteins which are indicators of HGT were found near the merA2 genes. All these indicated that merA2 genes of A. caldus arose from early HGT events, most likely from members within genus Acidithiobacillus.
Neighbor-joining (NJ) phylogenetic tree of the MerA protein sequences derived from Acidithiobacillus species strains and other representative species. Bootstrap values indicated at each node are based on a total of 1,000 bootstrap replicates. Branches representing group I mer clusters are in green, those representing group II mer clusters are in blue, and those representing group III mer clusters are in red and pink.
Case study of As-related genes.The overall distribution and organization of As resistance genes in Acidithiobacillus spp. are summarized in Fig. S11. As resistance genes (ars genes), including arsR, arsB, and arsC, were detected in all Acidithiobacillus genomes, while arsA and arsD were only detected in A. ferrivorans genomes and some strains from other species, including A. ferrooxidans strain DLC-5 and A. ferridurans IO-2C (A. ferrooxidans IO-2C), and arsH genes were only detected in genomes of A. ferrooxidans (within arsCRBH cluster) and A. thiooxidans (stand alone) (Fig. S11). Stand-alone arsM genes were also detected in all species. However, aio and arr, involved in arsenite oxidation and respiratory reduction of arsenate, were not found in Acidithiobacillus genomes, suggesting that cytoplasmic As(V) reduction and As(III) extrusion are the main As resistance strategies used in genus Acidithiobacillus spp. The ars genes in Acidithiobacillus genomes grouped together as ars clusters, including four subgroups, arsADCR, arsRB, arsBRC, and arsBRCH. The arsB gene of A. caldus strain MTH-04 was intercepted by an exogenous insertion sequence, which may result in malfunction of this arsB gene, rendering A. caldus MTH-04 vulnerable to toxic arsenic ions (Fig. S11).
(i) Origin and evolution of ars clusters. The distribution and organization of ars genes in Acidithiobacillus raise a question as to their evolution. The results showed that the G+C contents of the arsADCR clusters were higher than those of the genomes in A. ferrivorans strains (60.0% to 63.3% versus 56.4% to 58.6%, respectively); the G+C contents of respective arsBRC clusters were slightly higher than those of the genomes in A. caldus strains (62.3% to 63.1% versus 60.9% to 61.5%) and A. thiooxidans strains (53.1% to 55.3% versus 52.4% to 54.3%), showing the variation of G+C contents between clusters and the corresponding genomes (see Fig. S12). A comparison of the chromosomal regions flanking the ars gene clusters among Acidithiobacillus strains revealed that the genes in the upstream and downstream regions were conserved among strains of the same species (see Fig. S13). There were several genes encoding mobile element proteins (transposases or recombinases) located near the arsADCR and arsBRC clusters, and there were tRNA genes flanking the arsBRC clusters of A. caldus, indicating putative HGT events (Fig. S13). Strains of the same species share the same insertion sites, whereas strains of different species had different insertion sites, suggesting that ars clusters may have been acquired more than once via different HGT events. The flanking regions of the arsBRC gene clusters in strain A. thiooxidans were homologous to the corresponding regions of strains A. ferrivorans SS3 and A. ferrooxidans ATCC 23270 that have no arsBRC; the same phenomenon was found in arsBR gene clusters and arsH genes. These results suggested that ars genes may have been acquired by different HGT events followed by gene loss in some strains (Fig. S13).
To gain insights into the origin of ars genes clusters, NJ phylogenetic trees were constructed on the basis of on the concatenated ArsBRC protein sequences and ArsAD protein sequences. As shown in Fig. 4, the phylogenetic analysis revealed that the strains possessing arsBRC clusters formed a monophyletic clade. Acidithiobacillus spp. likely obtained the arsBRC cluster via HGT from Acidihalobacter, followed by rearrangement and gene loss, since ars gene clusters of Acidihalobacter prosperus branch near the base of the clade of Acidithiobacillus. NJ phylogenetic trees based on concatenated ArsAD protein sequences showed that the ars clusters of Acidithiobacillus and ars clusters of Alcaligenes faecalis, Pseudomonas hunanensis, and Sulfuriferula spp. were sister groups (Fig. 5), implying that Alcaligenes faecalis, Pseudomonas hunanensis, and Sulfuriferula spp. are likely donors of the arsADCR cluster.
Neighbor-joining phylogenetic tree of concatenated ArsBRC protein sequences derived from Acidithiobacillus species strains and other representative species. Bootstrap values indicated at each node are based on a total of 1,000 bootstrap replicates. Branches representing ars clusters of Acidithiobacillus spp. are in blue.
Neighbor-joining phylogenetic tree of concatenated ArsAD protein sequences derived from Acidithiobacillus species strains and other representative species. Bootstrap values indicated at each node are based on a total of 1,000 bootstrap replicates. Branches representing ars clusters of Acidithiobacillus spp. are marked in red.
Stand-alone arsM genes encoding As(III) S-adenosylmethyltransferase were found located distant from canonical ars gene clusters. NJ phylogenetic trees constructed on the basis of the ArsM protein sequences (see Fig. S14) revealed that the Acidithiobacillus strains formed three separate subgroups, of which subgroup I comprises ArsMs from A. thiooxidans, A. ferrooxidans, A. ferridurans, A. ferrivorans, and Acidithiobacillus sp. strain SH, grouping with ArsM from Thiomonas. Subgroup II comprises two deviant ArsMs from A. ferrooxidans and A. thiooxidans (GenBank accession no. WP_113525444.1 and WP_080707752.1, respectively), grouping with ArsMs from Melaminivora and Pseudomonas. Subgroup III comprises ArsMs from A. caldus, grouping with ArsMs from members of Rhodospirillales. Incongruences between the phylogenies of ArsMs and species imply that the subgroup I arsM genes may have been acquired via HGT from members of Thiomonas, subgroup II arsM genes from members of the genera Melaminivora and Pseudomonas, and subgroup III arsM genes from members of Rhodospirillales in early evolutionary history.
Case study of Cd-, Zn-, Co-, and Cu-related genes.Gene clusters encoding Czc/CusABC homologs, which are involved in the detoxification of divalent cations (cadmium and zinc and cobalt) (38, 39) and monovalent cations (copper and silver) (13, 40, 41), were detected in all Acidithiobacillus genomes (see Fig. S15). In addition, genes encoding CzcD and genes encoding copper-translocating P-type ATPase (Cop) were detected in all Acidithiobacillus genomes, most of which located next to czcABC gene clusters. However, genes that encode multicopper oxidase were only detected in strains of A. thiooxidans and A. caldus. The above-mentioned genes represent the dominant genetic traits that contribute resistance to Cd, Zn, Co, and Cu in Acidithiobacillus genomes. Considerable numbers of additional copies of czc-like and cop-like genes were also present. The czc clusters were divided into two subgroups, czcABC and czcCAB (Fig. S15). The gene clusters in the organization of czcABC were observed in strains across all species of Acidithiobacillus, of which strains of A. ferrooxidans and A. ferridurans possess other copies of czcABC located on different sites, suggesting that czc clusters may be acquired more than once. There are also tRNA genes and genes encoding mobile element proteins flanking czcABC clusters (Fig. S15). Each A. caldus strain possesses a czcCAB cluster in addition to an existing czcABC cluster (Fig. S15). The czc gene clusters exhibited more than 90% identity within each subgroup and over 80% identity between two subgroups. Actually, the czc genes were detected in all strains of Acidithiobacillus. The genes in the upstream and downstream regions of a czcABC cluster were conserved among strains of all species of Acidithiobacillus (Fig. S15A), which indicated that this czcABC gene cluster was acquired before the divergence of genus Acidithiobacillus. Notably, cop genes encoding copper-translocating P-type ATPase located upstream of the czcABC clusters are conserved across all species of Acidithiobacillus (Fig. S15A).
(i) Origin and evolution of czc clusters. The G+C contents of the czcABC clusters are higher than those of the genomes in Acidithiobacillus spp. (62.2% to 64.2% versus 52.9% to 61.5%, respectively) except in A. caldus, in which the G+C contents of the clusters are lower than of the A. caldus genomes (59.4% to 59.5% versus 60.9% to 61.4%, respectively); the G+C contents of respective czcBAC clusters are higher than those of the genomes in A. caldus strains (64.9% to 66.0% versus 60.9% to 61.5%, respectively), presenting variations of G+C contents between clusters and the corresponding genomes (see Fig. S16). To gain insights into the origin of czc genes clusters in Acidithiobacillus, a NJ phylogenetic tree was constructed on the basis of concatenated CzcABC protein sequences. As shown in Fig. 6, the strains possessing czcABC and czcBAC clusters formed separate groups. Notably, the phylogeny reveals that the czcABC of Acidithiobacillus and czc clusters of Thiomonas spp. and Acidihalobacter spp. are sister groups, whereas czcBAC of Acidithiobacillus and czc clusters of another clade of Thiomonas spp. are sister groups. These results implied that the czcABC cluster may have been acquired via HGT from members of Acidihalobacter and Thiomonas, while czcBAC may have been acquired from members of Thiomonas.
Neighbor-joining phylogenetic tree of concatenated CzcA, CzcB, and CzcC protein sequences derived from Acidithiobacillus species strains and other representative species. Bootstrap values indicated at each node are based on a total of 1,000 bootstrap replicates. Branches representing czcABC clusters are in green; czcBAC clusters formed in another group are in blue.
(ii) mco genes in Acidithiobacillus spp. Stand-alone mco genes that encode multicopper oxidase, which oxidizes Cu+ to the less toxic Cu2+ form, were found in all strains of A. thiooxidans and A. caldus in this study. An NJ phylogenetic tree based on the multicopper oxidase protein sequences revealed that the multicopper oxidases of Acidithiobacillus were adjacent to those of Acidihalobacter spp., implying that Acidithiobacillus and Acidihalobacter have a common multicopper oxidase ancestor (see Fig. S17). Synteny of the chromosomal regions flanking the mco genes showed that multicopper oxidase was probably lost in strains of A. ferrivorans, A. ferridurans, and A. ferrooxidans (see Fig. S18).
(iii) czc-like and cop-like genes in Acidithiobacillus. Our studies revealed that there are multiple czc-like genes and cop-like genes present in the 37 genomes, including 1 to 5 pairs of czcA-like genes, 1 to 5 pairs of czcB-like genes, 1 to 4 pairs of czcC-like genes, 1 to 4 pairs of czcD-like genes, and 1 to 3 pairs of cop-like genes (see Fig. S20, S21, S22, S23, and S24). An alignment of residues of CzcA and CzcA-like proteins revealed that the “DDE” motifs essential for CzcA function (42) are highly conserved (see Fig. S25). Alignments of residues of Cop and Cop-like proteins revealed that the metal binding domains (CASC …. CASC), translocation domains (CPCAMGLA), and phosphorylation domains (FDKTGTLT) (43) are conserved, indicating that these genes may be functional (see Fig. S26). A phylogenetic analysis revealed that czc-like and cop-like genes from the same species clustered together and that there were incongruences between the phylogenies of functional genes and species, which suggested that multiple czc-like and cop-like genes may have arisen from gene duplication and HGT (Fig. S20, S21, S22, S23, and S24).
Assessment of functionality of metal resistance genes.We also assessed indirectly the functionality of above-mentioned metal-resistance-related genes by assessing the strength of natural selection acting on these genes, as well as the codon adaption indexes (CAIs) of these genes. The results showed that 83.3% of the above-mentioned metal resistance genes had a rate of nonsynonymous substitutions that was lower than the rate of synonymous substitutions (dN/dS < 1), suggesting purifying selection on these genes (Fig. 7 and Table S6). Only a merA gene from Acidithiobacillus thiooxidans strain A02 (dN/dS = 1.51), a cop gene from Acidithiobacillus ferrooxidans strain CCM 4253 (dN/dS = 1.39), and an arsH gene from Acidithiobacillus ferridurans JCM 18981 (dN/dS = 1.82) showed dN/dS ratios of >1, indicating that they might be under positive selection. The lowest dN/dS ratio was observed for arsB genes (average dN/dS = 0.06) and arsM genes (average dN/dS = 0.04), showing strong purifying selection. All these metal resistance genes with dN/dS ratios of <1 were supposed to be most essential for the growth of Acidithiobacillus spp. However, merA genes from A. ferrooxidans, arsR genes from A. thiooxidans, merC genes from A. caldus, merP genes from A. thiooxidans and A. ferridurans, merC genes from A. caldus, mco genes from A. thiooxidans, and some cop-like and czc-like genes showed dN/dS ratios of ≈1, indicating that the selection force had relaxed on them (Fig. 7 and Table S6). Highly expressed genes in bacteria often have a stronger codon bias than genes expressed at lower levels, due to translational selection (44–46). In this study, the expression levels of metal-resistance-related genes in the Acidithiobacillus spp. were predicted using the CAI as a numerical estimator (47, 48). Approximately 53.2% of the metal resistance genes were predicted to be highly expressed genes using CAI cutoff values of greater than 0.73 in A. thiooxidans, greater than 0.75 in A. caldus, and greater than 0.76 in A. ferrooxidans, A. ferrivorans, and A. ferridurans (Fig. 8 and see Table S7). The cutoff values were indicated with average CAI values of Tu genes in each species.
Distributions and ranges of selection pressures on different metal resistance genes of Acidithiobacillus spp.
Ranges of CAI values of different metal resistance genes of Acidithiobacillus spp. with Tu gene as a reference.
DISCUSSION
In the present work, we characterized the genomic features of Acidithiobacillus spp., focusing on the distribution, organization, complex evolutionary history, and functionality of metal resistance genes. With the rapid advances in genomic analysis technology, new measurements such as ANI are being developed to evaluate the genomic similarity between bacteria (33). In this study, we redetermined the phylogenetic status of 38 Acidithiobacillus spp. by using ANI and phylogenetic trees. Our results showed that strains DSM 14366, GGI-221, and IO-2C were misnamed. We therefore reclassified strain A. ferrooxidans IO-2C as A. ferridurans IO-2C, A. albertensis DSM 14366 as A. thiooxidans DSM 14366, and Acidithiobacillus sp. strain GGI-221 as A. ferrooxidans GGI-221. Low ANI values for Acidithiobacillus sp. strain NORP59 in comparison to those of other Acidithiobacillus strains and a distant phylogenetic relation indicated that strain NORP59 is not a member of genus Acidithiobacillus. Our study revealed the complex taxonomy of genus Acidithiobacillus and demonstrated the usefulness and accuracy of advanced bioinformatics measurements such as ANI in rectifying inaccurate identifications from previous studies.
There is generally a positive correlation between the abundance of mobile genetic elements and the frequency of HGT (49). We found that Acidithiobacillus spp. contain 0 to 6 prophages and/or prophage remnants, 143 to 334 insertion elements (IS), and 25 to 66 GIs per genome, suggesting that Acidithiobacillus spp. are rich in mobile genetic elements. The existence of transposable elements and prophages near the heavy metal resistance genes and gene clusters suggested that they may be involved in HGT of heavy metal resistance genes, which is consistent with previous studies illustrating that HGT was a major mechanism that confers upon bacteria crucial traits such as metal resistance for adaptation to toxic metal-containing habitats (15, 50–52). In addition, BacMet database annotation revealed that the accessory genome had a higher proportion (8.2% [7,846 genes]) of genes related to heavy metal resistance activity (Fig. 1C), and models of extrapolation revealed that the pangenome of 37 Acidithiobacillus spp. is “open,” which is in line with the results showing that metal-resistance-associated gene family gain events overcame loss events. We inferred that during evolutionary niche adaption, Acidithiobacillus spp. accessorized their genome armaments with genetic elements transferred from members outside their species or genus, with the help of MGEs. The “accessory genome” was therefore considered a massive resource that provided Acidithiobacillus spp. with unprecedented elasticity to improve their fitness while adapting to harsh conditions (53–55). Our results also showed that the transformation of heavy metals to less toxic forms followed by their efflux out of cells was an often-used metal resistance strategy of Acidithiobacillus spp., which is also a process contributing to the biogeochemical cycling of metals.
We implemented a combination of several approaches to discover putative HGT events, including deviant G+C content detection, phylogenetic analysis, mobile genetic element detection, and chromosomal synteny analysis, since it was difficult to identify HGT events via deviant G+C contents alone if HGT happened to occur between organisms with similar G+C contents (56). The different gene contents and organization of metal resistance genes suggested a complicated evolutionary history of these genes and also suggested different genetic requirements for effective cellular metal detoxification and regulation of metal resistance operons. We supposed Thiomonas spp. to be the donors of the arsM gene and czcABC and czcBAC gene clusters in Acidithiobacillus spp., which was reasonable since members of Thiomonas spp. are ubiquitous in heavy-metal-contaminated environments, especially in AMD sites, the same habitats as Acidithiobacillus species, and usually contained a vast flexible gene pool. For instance, the genome of Thiomonas arsenitoxydans contained a circular chromosome and a plasmid (pTHI) where more than 20 GIs (ThGEI-A to ThGEI-S) resided, which were associated with heavy metal resistance and conjugation (57–59). In addition, the median G+C content of Thiomonas species genomes was 64.0%, which was close to those of czcABC (64.2%) and czcBAC (64.9%) gene clusters in Acidithiobacillus spp. All these indicate that Thiomonas spp. are the most possible sources of multiple heavy metal resistance operons that were acquired by Acidithiobacillus spp. via cross-phylum HGT. Similar traits were also shown in case of Acidihalobacter spp., the putative cross-phylum donors of merACR, arsBRC, and the multicopper oxidase-encoding gene, which are a group of halophilic, mesophilic, heavy-metal-tolerant extreme acidophiles that share the same habitats with Acidithiobacillus spp., with an average genome G+C content (61.9%) similar to those of merACR (60.0%) and arsBRC (61.7%) gene clusters in Acidithiobacillus spp. (60, 61). This could also be applied to Acidiferrobacter spp. that inhabit AMD sites, which were predicted to be the cross-phylum donors of the merAPTR gene cluster and rusticyanin-encoding gene (62). Multiple czc-like genes and cop-like genes were present in the genomes of Acidithiobacillus spp., which may arise from gene duplication. Metal resistance proteins encoded by genes such as czc and cop are in high demand for Acidithiobacillus spp. for coping with harsh metal-rich environments. The presence of duplicate metal resistance genes may be beneficial by enhancing the regulatory elasticity of Acidithiobacillus spp., simply because extra amounts of protein or RNA products are generated (63).
The gene contexts of metal resistance genes in different Acidithiobacillus species were not conserved, which indicated that these genes may be acquired more than once via different HGT events. Synteny of the chromosomal regions flanking the mco genes showed that multicopper oxidases were possibly lost in strains of A. ferridurans, A. ferrivorans, and A. ferrooxidans during evolution, since mco genes in A. thiooxidans experienced relaxed selection pressure (dN/dS ≈ 1) and A. thiooxidans arose before the speciation of A. ferridurans, A. ferrivorans, and A. ferrooxidans in evolutionary history (see Fig. S7 in the supplemental material). We supposed that the function of this enzyme was replaced in A. ferridurans, A. ferrivorans, and A. ferrooxidans by rusticyanin, a periplasmic protein with an affinity for Cu+, forming part of an iron-oxidizing supercomplex and supposed to confer a higher level of copper resistance (64, 65). Rusticyanin was absent in species incapable of oxidizing iron, such as A. thiooxidans and A. caldus, and was acquired via HGT from Acidiferrobacter spp. to A. ferridurans, A. ferrivorans, and A. ferrooxidans, as revealed by a phylogenetic analysis (see Fig. S19). Most nonsynonymous substitutions are deleterious, especially for functionally important genes, since the functionality of important proteins might be substantially affected even if small changes happens in its sequences, resulting in a dN/dS if <1 for these genes considering that mutants of these genes are prone to be eliminated by purifying (negative) selection. On the contrary, when new mutations are beneficial, conferring adaptive advantages, these genes might be under diversifying (directional or positive) selection, reflected by a dN/dS of >1, namely, there are more nonsynonymous than synonymous nucleotide substitutions. Nonsynonymous nucleotide substitutions and synonymous substitutions will accumulate at similar rates (dN/dS ≈ 1) if genes become nonfunctional in the population, since there is no selection pressure acting to maintain the adaptive variant (neutral evolution), and these genes are more prone to loss during evolution (66). The results showed that 83.3% of metal resistance genes in Acidithiobacillus species had a lower rate of nonsynonymous substitution than of synonymous substitution (dN/dS < 1), which is in consistent with the assumption of predominant purifying selection on functionally important metal resistance genes. Besides, approximately 53.2% of the metal resistance genes were predicted to be highly expressed via the CAI (codon adaption index). There is a close correlation between codon usage bias and gene expressivity, which results from the selection for efficient and accurate translation of highly expressed genes (67, 68). The CAI was used to measure the synonymous codon usage bias for a given DNA sequence with an ever-improving algorithm by comparing the similarity between the synonymous codon usage of a gene and the synonymous codon frequency of a reference data set of more than 30 highly expressed genes derived from prokaryotes and lower eukaryotes, which can serve as a reliable indicator of codon usage bias as well as a numerical estimator of gene expressivity for genes of prokaryotes and lower eukaryotes (65, 69). Widespread low dN/dS ratios and relative high CAI values across metal-resistance-related genes emphasized their indispensable functions in supporting the growth of these extremophiles in the face of harsh environments. Altogether, our research revealed that Acidithiobacillus spp. recruited and consolidated additional novel resistant functionalities during adaption to challenging metal-rich environments via HGT, gene duplication, and purifying selection. However, further experimental confirmations of the functionalities and expression levels of metal-resistance-related genes in Acidithiobacillus spp. are still necessary for an advanced understanding of these genes.
Concluding remarks.We performed a comparative genomic analysis of 37 strains within the genus Acidithiobacillus to explore the distribution, organization, functionality, and complex evolutionary history of metal resistance genes in Acidithiobacillus spp. The results indicated that the evolutionary history of toxic metal resistance genes in Acidithiobacillus spp. involved a combination of gene gain and loss, HGT, gene duplication, and gene family expansion. Gene clusters involved in toxic metal resistance, such as mer, ars, and czc-cus, were found to be acquired by early horizontal gene transfer (HGT) events from sources related to Burkholderiales, Rhodospirillales, Acidihalobacter, Thiobacillus, Acidiferrobacter, Alcaligenes, Pseudomonas, Sulfuriferula, Melaminivora, and Thiomonas, many of which share habitats with Acidithiobacillus spp., whereas the multicopper oxidase genes of iron-oxidizing A. ferridurans, A. ferrivorans, and A. ferrooxidans involved in copper detoxification were lost and replaced by rusticyanin genes during evolution. An assessment of the functionality of metal resistance genes showed that 83.3% of metal resistance genes in Acidithiobacillus spp. were under purifying selection (dN/dS < 1) and 53.2% of them were highly expressed. The insights gained in this study have greatly improved our understanding of the metal resistance strategies of Acidithiobacillus spp.
MATERIALS AND METHODS
A total of 38 sequenced strains belonging to genus Acidithiobacillus were collected from the NCBI database, containing one to three sequenced standard strains in each species and four unidentified strains. The general features of the bacterial genomes used in this study are summarized in Table 2.
Average nucleotide identity.Comparisons of average nucleotide identities (ANIs) based on a BLAST algorithm (70) and tetranucleotide frequency correlation coefficient (Tetra) (71) were conducted using the web server JSpeciesWS (33, 72) with default parameters.
Phylogenic analyses.To associate the distribution of metal resistance genes in Acidithiobacillus spp. with their phylogenetic affiliations, we constructed the phylogenetic tree of the Acidithiobacillus spp. with closely related species on the basis of 16S rRNA gene sequences using the neighbor-joining method in MEGA v7.0 (73) and with rRNA genes predicted using the RNAmmer v1.2 (74). Since no 16S rRNA gene sequence could be retrieved from the genome of Acidithiobacillus sp. strain NORP59, we also constructed a phylogenetic tree of the 38 Acidithiobacillus species genomes with their phylogenetic affiliations on the basis of whole-genome sequences with CVTree3 (75). A phylogenetic tree of the 37 Acidithiobacillus species genomes on the basis of the concatenation of the 110 single-copy core genes in a genome was constructed with the maximum likelihood (ML) method using MEGA v7.0 (73) with Jones-Taylor-Thornton (JTT) model (76) and 1,000 bootstraps, rooted by Escherichia coli and Salmonella enterica. The ambiguous areas of alignment were removed by Gblocks 0.91 b (77) before constructing the tree. A neighbor-joining tree and UPGMA tree were also constructed with MEGA v7.0 (73), both under a p-distance model with 1,000 bootstraps of the above-mentioned core genes. A chronogram for the NJ phylogenetic tree of the 37 Acidithiobacillus species with branch lengths reflecting divergence times was inferred using the divergence time between Escherichia coli and Salmonella enterica (140 MYA [millions of years ago]) (78) as calibration points with Timetree Wizard, a built-in program of MEGA v7.0 (73) with the JTT matrix-based model (76) and the real-time method (79).
Genes for toxic metal resistance.Genes involved in heavy metal resistance in Acidithiobacillus genomes were identified by performing BLASTP searches against the BacMet database (80), which contains genes with experimentally confirmed metal resistance functions. Genes meeting the following cutoffs were considered reliable hits: E value, <1e−50; sequence similarity, >35%. The evolutionary history of genes for toxic metal resistance was inferred using MEGA v7.0 (73) with the neighbor-joining method and 1,000 bootstraps.
Pangenome analysis and model extrapolations of Acidithiobacillus.The Bacterial Pan Genome Analysis (BPGA) pipeline (32) was used to identify orthologous groups among Acidithiobacillus testing genomes and to extrapolate the pangenome models of Acidithiobacillus spp. by applying default parameters. The set of genes shared among all testing strains was defined as their core genome, the set of genes shared with more than two but not all testing strains was defined as their accessory genome, the set of genes not shared with other strains in each testing strain was defined as the unique genes, and the set of nonhomologous genes with all testing genomes was defined as the pangenome. The details of all testing strains used are listed in Table 2. Genes involved in heavy metal resistance in the pangenome, core genome, and accessory genome and unique genes were identified by performing BLASTP searches against the BacMet database (80) with cutoffs of an E value of <1e−50 and a sequence similarity of >35%.
In this study, the size of the Acidithiobacillus pangenome was extrapolated by implementing an power law regression function, Ps = κnγ, using a built-in program of the BPGA pipeline (32), in which Ps represents the total number of nonorthologous gene families within its pangenome, n represents the number of tested strains, and both κ and γ are free parameters (35). An exponent γ of <0 suggests the pangenome is “closed,” where the size of the pangenome reaches a constant value as extra genomes are added. Conversely, the species is predicted to harbor an open pangenome for γ values between 0 and 1. In addition, the size of the core genome was extrapolated by fitting into an exponential decay function, Fc = κcexp(-n/τc) + Ω, with a built-in program of the BPGA pipeline (32), where Fc is the number of core gene families, and κc, τc, and Ω are free parameters (81).
Estimation of gene family gains and losses.BadiRate (82) was used with the “BDI-FR-CWP” method to estimate the number of metal-resistance-related gene family gains and losses during the evolution of Acidithiobacillus. The analysis was conducted only on 5 representative Acidithiobacillus strains, including A. thiooxidans CLST, A. ferrooxidans ATCC 53993, A. ferrivorans PQ33, A. ferridurans JCM 18981, and A. caldus SM-1, because the estimation required complete sets of testing gene families, which are only available in species with considerable phylogenetic distance and have high-quality (complete or nearly complete) genome databases. The orthologous groups of metal resistance genes based on reciprocal best hits within each Acidithiobacillus species was inferred using the CD-HIT software (global sequence identity cutoff was 0.9 and bandwidth of alignment cutoff was 20 amino acids [aa]) (83). These data were used to construct a gene family size file (a matrix [column, species; row, number of orthologous groups]). The topology of the core-gene tree was used as the reference tree. The number of gene families related to metal resistance in each internal node was inferred using those numbers in branch species and the phylogenetic branch lengths. Gene family gains and losses in each phylogenetic branch were then estimated using the Wagner parsimony algorithm under the BDI stochastic model and free rates (FR) model, assuming that each branch has its own turnover rates.
Prediction of mobile genetic elements.The ISFinder platform (84) was used to predict and classify insertion sequences (IS) and transposases distributed over Acidithiobacillus genomes, with default criteria. The IslandViewer 4 platform (85), which integrates two prediction methods that make use of sequence composition, i.e., IslandPath-DIMOB (86) and SIGI-HMM (87), and a comparative genomic islands prediction method IslandPick (88), were used to identify putative genomic islands (GIs). PHASTER (Phage Search Tool Enhanced Release) (89) was used for the identification and annotation of prophage sequences within Acidithiobacillus genomes. In addition, tRNA genes were predicted using the tRNAscan-SE v2.0 (90). The genes involved in heavy metal resistance were visualized, along with the tRNAs, insertion sequences, prophage sequences, and putative genomic islands, as well as the core genes, accessory genes, and unique genes, using BRIG-0.95 (91).
Selective pressure analysis and expressivity prediction.The strength of natural selection was assessed by estimating the relative rates of nonsynonymous (dN) to synonymous (dS) nucleotide substitutions of orthologous metal resistance genes by using the HyPhy package (92) with the adaptive branch-site random effects likelihood (aBSREL) method (93) via the datamonkey web server (94). A gene in the node or tip of a given tree was considered under positive (diversifying) selection (dN/dS > 1) or negative (purifying) selection (dN/dS < 1) if the P value was <0.05 using the likelihood ratio test after correcting for multiple testing. The CAI (codon adaption index) was used as a numerical estimator of gene expression level (47, 48), a method previously described by Wu et al. (95). CAIcal (96) was used to calculate the CAI values of the above-mentioned metal resistance genes. CAI values range from 0 to 1.0, with a higher CAI value indicating a higher expression level. Putative highly expressed (PHX) genes related to toxic metal resistance in the Acidithiobacillus spp. were inferred, with Tu genes as the references, which encode elongation factors that are supposed to be highly expressed across most organisms (97).
ACKNOWLEDGMENTS
We thank the NCBI for providing the genome sequences of Acidithiobacillus species strains.
This work was funded by the National Natural Science Foundation of China (NSFC; grants 31570113, 41807332, 41877345, 91851206 [Major Research Plan of "Mechanism Underlying Elemental Cycling on the Earth by Microorganisms in Hydrosphere" launched by NSFC], and 41573072) and open funding of the State Key Laboratory of Comprehensive Utilization of Low-Grade Refractory Gold Ores (Zijin Mining Group Co., Ltd.; no. 738010212).
We also thank the Hunan International Scientific and Technological Cooperation Base of Environmental Microbiology and Application, China. We declare no competing interests.
FOOTNOTES
- Received 5 September 2018.
- Accepted 19 October 2018.
- Accepted manuscript posted online 2 November 2018.
Supplemental material for this article may be found at https://doi.org/10.1128/AEM.02153-18.
- Copyright © 2019 American Society for Microbiology.
REFERENCES
- 1.↵
- 2.↵
- 3.↵
- 4.↵
- 5.↵
- 6.↵
- 7.↵
- 8.↵
- 9.↵
- 10.↵
- 11.↵
- 12.↵
- 13.↵
- 14.↵
- 15.↵
- 16.↵
- 17.↵
- 18.↵
- 19.↵
- 20.↵
- 21.↵
- 22.↵
- 23.↵
- 24.↵
- 25.↵
- 26.↵
- 27.↵
- 28.↵
- 29.↵
- 30.↵
- 31.↵
- 32.↵
- 33.↵
- 34.↵
- 35.↵
- 36.↵
- 37.↵
- 38.↵
- 39.↵
- 40.↵
- 41.↵
- 42.↵
- 43.↵
- 44.↵
- 45.↵
- 46.↵
- 47.↵
- 48.↵
- 49.↵
- 50.↵
- 51.↵
- 52.↵
- 53.↵
- 54.↵
- 55.↵
- 56.↵
- 57.↵
- 58.↵
- 59.↵
- 60.↵
- 61.↵
- 62.↵
- 63.↵
- 64.↵
- 65.↵
- 66.↵
- 67.↵
- 68.↵
- 69.↵
- 70.↵
- 71.↵
- 72.↵
- 73.↵
- 74.↵
- 75.↵
- 76.↵
- 77.↵
- 78.↵
- 79.↵
- 80.↵
- 81.↵
- 82.↵
- 83.↵
- 84.↵
- 85.↵
- 86.↵
- 87.↵
- 88.↵
- 89.↵
- 90.↵
- 91.↵
- 92.↵
- 93.↵
- 94.↵
- 95.↵
- 96.↵
- 97.↵
- 98.
- 99.
- 100.
- 101.
- 102.
- 103.
- 104.
- 105.