Metabolic and genomic analysis elucidates strain-level variation in Microbacterium spp. isolated from chromate contaminated sediment

Sample collection and isolation

Soil was collected in Seymour, IN from a Department of Transportation site known to have contamination from chromium, lead, and organic solvents (Kourtev, Nakatsu & Konopka, 2006; Nakatsu et al., 2005). Bill Jervis from the Indiana Department of Transport provided site access (no permit was required) and the project did not involve endangered or protected species. Bacterial strains were isolated from the contaminated soil as described in Kourtev, Nakatsu & Konopka (2009). Briefly, isolates were enriched on 50% trypic soy agar (TSA) and 0.25 mM Cr(VI). Colonies were picked from the TSA plates and placed in 50% TSB with .25 mM Cr(VI). Individual isolates were grown at varying amounts of chromate enriched TSA and tryptic soy broth (TSB). The isolates were allowed to grow and then selected based on their varying abilities to resist and reduce Cr(VI). Isolates were maintained on 50% TSA with 2 mM Cr(VI) as well as in 30% glycerol stocks stored at −80 °C to help preserve active pure cultures.

DNA extraction

Bacterial isolates used in this study (Cr-K1W, Cr-K20, Cr-K29, and Cr-K32) were grown in 250 ml of 50% TSB amended with 2 mM potassium chromate (K₂CrO₄). Following inoculation, isolates were incubated at 30 °C, 225 rpm for 24–72 h. Cells were harvested by centrifuging the cultures in 250 ml Nalgene bottles at 8,000 × g for 18 min. The supernatant was removed and the pellet was washed with autoclaved nanopure water following the same procedure. Cell pellets were extracted using the FastDNA^® Spin Kit (MP Biomedical, Santa Ana, CA) with one modification. Before extraction, cell pellets were resuspended in autoclaved nanopure water before 200 µl was transferred to the Lysis Matrix tube provided by the kit. DNA samples were recovered in 100 µl of DES solution and stored at −20 °C.

Chromate resistance and reduction experiments

All growth and reduction experiments were conducted in 50% TSB at 30 °C while shaking at 225 rpm. Chromium was added when needed in the form of K₂CrO₄. Growth was measured on a UV-Vis spectrophotometer (Varian Cary 50, Agilent) at 600 nm. To test the ability of isolates to reduce Cr(VI), chromate was added to the media to a final concentration of 2 mM. Cultures, inoculated from −80 °C stocks, were grown overnight and then used to subculture fresh medium to an initial OD₆₀₀ of 0.004. Cultures were prepared in triplicates for each growth and reduction experiment. OD₆₀₀ and Cr(VI) reduction for each culture were measured at 6, 12, 24, 48, 72, 96, and 120 h. A Cr(VI) reduction assay was performed as described by Urone (1955), with minor modifications. Briefly, one mL of each culture was centrifuged at 7,000 rpm for seven minutes to remove biomass. Ten µL of the supernatant were added to 966 µL of deionized water to which 3.40 µL of sulfuric acid was added. 1,5-Diphenylcarbazide (20 µL) was added to the sample and allowed to set for 10 min for full color development. Readings were taken on a Cary UV-Vis spectrophotometer (Agilent technologies, Santa Clara, CA) at 540 nm. Negative (uninoculated) controls of 50% TSB with 2 mM K₂CrO₄ were used to determine TSB-specific abiotic reduction of Cr(VI). In addition, relative Cr(VI) reduction rates were calculated following Zhu et al. (2008).

Resistance determinations were performed in 5 ml of 50% TSB to which K₂CrO₄ was added to final concentrations of 5, 10, 20, 40, 60, 80, and 100 mM. Cultures were prepared in triplicate for each resistance experiment and the OD₆₀₀ was taken at 0 and 96 h. Growth was determined against an abiotic control.

Metabolic screening

All four isolates were screened for their metabolic potential using Biolog’s GP2 plates (Hayward, CA) following the manufacturer’s protocol with minor modifications. Specifically, isolates were grown on Biolog Universal Growth (BUG) agar (Biolog, Hayward, CA) at 30 °C after which biomass was transferred to 25 ml of inoculating fluid (IF) (Biolog, Hayward, CA) containing sodium thioglycolate to achieve an OD₆₀₀ between 0.68 and 0.75. An aliquot (150 µl) was placed into each well and then the plates were incubated for 22 h at 30 °C. Plates were read using Biolog’s OmniLog program (version 1.2.01) which corrects for any difference in OD₆₀₀. Absorbance values were converted to 0(-) for negative, 1(+) for positive, and .05(+/ − ) for borderline by the Biolog OmniLog program. Replicates with a minimum of two positives or borderlines or one positive and one borderline were considered positive for the metabolite. A total of 95 metabolites were tested per isolate (Table S6). To ease data analysis, metabolites were grouped into 6 guilds (Zak et al., 1994): (1) carbohydrates, (2) carboxylic acids, (3) amino acids, (4) amines and amides, (5) polymers and (6) miscellaneous. The total number of positives for each guild was then calculated for each isolate and plotted using principal components analysis (PCA) in PAST3 (Harper & Ryan, 2001).

Sequencing, de novo assembly, and analysis

Whole genome shotgun sequencing was performed by multiplexing the genomic DNA onto one lane using the Illumina HiSeq 2000 platform with 100 bp paired end reads using V2 chemistry at Cincinnati Children’s Hospital Medical Center’s Genetic Variation and Gene Discovery Core facility. Raw Illumina genomic reads were trimmed of their adapter sequences using the default setting of the program Trimmomatic (version 0.27) (Lohse et al., 2012). Trimmed reads were checked for quality using FastQC (version 0.10.2; http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and then trimmed for quality using the fastx_ trimmer (-Q33 -l 70) and fastx_ quality_ filter (-Q33 -q 30 -p 50) functions of the FastX toolkit program (version 0.13.2) (http://hannonlab.cshl.edu/fastx_toolkit/). To determine the appropriate range of Kmer length for assembly, the program KmerGenie (version 1.5856) was used (Chikhi & Medvedev, 2013; Zerbino & Birney, 2008). Cleaned pair end reads were assembled using the default setting of the assembly programs Velvet (version 1.2.10) (Zerbino & Birney, 2008) and AbySS (version 1.3.6) (Simpson et al., 2009) at a range of Kmer surrounding the estimated Kmer size from KmerGenie. A third assembly was completed using the raw pair-end reads and the intergrated a5 pipeline assembly (Tritt et al., 2012). This pipeline automates the processes of data cleaning, error correction, contig assembly, scaffolding, and quality control. The draft assemblies were compared and the best assembly picked for each isolate based on their total contigs, N50, genome size, max contig length, and mean contig length. AbySS assemblies for all four isolates were used for all downstream analyses (Table S3). The average coverage of each genome was >200×, with each assembly containing 30–81 contigs (Table S2). The contigs from the selected assembly were annotated and analyzed using the Departments of Energy’s Joint Genome Institutes IMG program (Markowitz et al., 2012). Utilizing IMG’s data annotations, Pfam categories and their broad category gene counts for the respective genomes were extracted and normalized based on the sum of each row. Principal component analysis (PCA) was then performed using PAST3. All protein coding genes from the annotated draft isolate genomes were submitted to the Pacific Northwest National Laboratory’s Species Parallel and Orthology Solver (SPOCS) for analysis to determine pairs of orthologous and paraorthologous proteins between the closely related isolates (Curtis et al., 2013).

Analysis of genomes for the presence of chromate-related genes

To increase the potential for finding chromate related genes, assembled contigs from the four representative isolate genomes were searched against a protein database of chromate related genes acquired from the UniprotKB database (http://www.uniprot.org/). Three databases were downloaded containing the chromate related protein sequences of ChrA (efflux pump), ChrR (reductase), and YieF (reductase). Databases were obtained from text-based queries within the UniprotKB database. Because IMG based annotations did not return any chromate-related genes, assembled contigs were searched manually against individual databases with the BLASTX algorithm using a minimum e-value cutoff of 1e–05. Resultant nucleotide sequences with corresponding hits were translated using the translation tool from ExPAsy (http://expasy.org/) and examined to see if the sequences were found within a larger open reading frame (ORF). For each sequence, the top sequence corresponding to an annotated protein with the highest bit and query coverage was selected. The resultant amino acid sequence was searched using the protein BLAST tool from NCBI to further validate the gene candidate based on domain hits as shown within BLAST. Sequences were further scrutinized by examining the top resultant hit within IMG using their alignment tools and gene neighborhood viewer.

Phylogenetic trees of the putative genes were made by using classic chrR and yieF genes downloaded from UniProtKB as well as from the top hits from homology searches within IMG. The protein sequences were aligned using ClustalW. A phylogenetic tree was then constructed using Maximum likelihood (bootstrap = 50) within the software MEGA (version 6.06) (Tamura et al., 2007).

Bacterial isolation and characterization of chromate resistance and reduction

Four bacterial strains isolated from chromium contaminated soil samples (collected in Seymour, IN) were studied for their ability to resist and reduce Cr(VI). When introduced to Cr(VI), each isolate was able to grow and survive in a minimum of 2 mM chromate but differences in their growth were noted. Overall, three isolates (Cr-K20, Cr-K29, and Cr-K32) grew to relatively dense cultures (OD ≥ 4.0) and had similar growth rates whether grown aerobically with or without Cr(VI) (Figs. 1A–1D and Table S1). In contrast, Cr-K1W doubled in optical density when grown without Cr(VI) and it reached exponential growth earlier than the other isolates (Fig. 1A). When experiments were done with Cr(VI), Cr-K1W was also able to maintain a stable stationary phase for over 100 h, which was not observed for the other isolates under similar conditions. Chromate toxicity appeared to have the greatest effect on Cr-K1W when compared to the other three isolates.

The isolates also had different abilities to reduce and resist chromate stress. Cr-K29 and Cr-K32 were capable of reducing 2 mM chromate within 48 h (Figs. 1C and 1D) and had the relatively fastest reduction rates (Table S1), while resisting up to 100 mM of chromate. Conversely, Cr-K1W and Cr-K20 were only able to resist up to 10 mM and reduced less than 1.2 mM of chromate after 120 h (Figs. 1A and 1B). These two isolates also had the lowest relative Cr(VI) reduction rates (Table S1). Taken together, a connection between reduction and resistance was evident. Chromate resistance has been shown to be connected to the efflux pump gene chrA, part of the CHR superfamily, (e.g., Diaz-Perez et al., 2007; Henne et al., 2009a; Ramirez-Diaz et al., 2008), which can provide a wide range of resistance, 0.3 mM–200 mM (Viti et al., 2013). However, the link between reduction and resistance in these isolates may point to the importance of a Cr(VI) reduction mechanism and other non-chromate specific responses in aiding resistance, instead of Cr(VI) resistance being driven solely by an efflux pump.

Metabolic fingerprinting of the isolates

Metabolic fingerprinting found variation between that isolates that is comparable to their growth differences. A survey of 95 potential metabolites showed a core group of 32 metabolites used among all four isolates (Fig. S1). Furthermore, isolates Cr-K1W and Cr-K20 shared the ability to use an additional seven metabolites, while Cr-K29, Cr-K20, and Cr-K32 shared the ability to use 16 additional metabolites (Fig. S1). Of the four isolates, Cr-K20 had the most unique combination of metabolites utilized (n = 16). Both Cr-K29 and Cr-K32 exhibited the largest amount of overlap in substrate utilization, which may be related to why these isolates had very similar growth patterns (Figs. 1B and 1C). The growth patterns of Cr-K1W and Cr-K20 were unique when compared to the other isolates, which was mirrored by their metabolic fingerprint. Of the six-substrate guilds, Carbohydrates and the Miscellaneous guilds exerted the greatest influence on separating Cr-K1W and Cr-K20 from Cr-K29 and Cr-K32 (Fig. 2). The four remaining substrate guilds (Amines/Amides, Carboxylic acids, Amino Acids, and Polymers) appeared to help separate Cr-K1W away from the other three isolates (Fig. 2).

Genomic characterization of chromate reducing isolates

The four draft genomes obtained in this study had lengths between 3.79 and 3.91 Mbp with a GC content of ∼68% (Tables S2 and S3). Analysis of the 16S rRNA genes indicated that all four isolates had identical (1,405 bp) 16S rRNA sequences, and that they are Microbacterium sp. (Actinobacteria phylum), specifically members of the Microbacterium oxydans clade (Fig. S2). Gene annotation of the four genomes documented between 3,616 and 3,806 predicted protein coding genes (Table S2). Genomes were assessed to be nearly complete with all genomes containing the 35 universal single copy marker genes (Raes et al., 2007) (Table S4).

Analysis of the genomes indicated that the four taxonomically closely-related chromium reducing Microbacterium isolates have a surprising amount of inter-strain genomic variation. A broad based genome comparison of the four isolates showed that they share a large predicted protein core made up of 2,810 proteins (Fig. S3). In addition to the protein core, isolates Cr-K29 and Cr-K32 shared an additional 602 predicted proteins while isolates Cr-K1W and Cr-K20 shared 854 predicted proteins (Fig. S3). When the annotated genomes were analyzed by Principal Component Analysis, using annotated Pfam proteins (Table S5), a similar grouping pattern was observed (Fig. 3). Of the 21 Pfam categories, Carbohydrate Transport and Metabolism, Transcription, and Inorganic Ion Transport and Metabolism, General Function, and Function Unknown were the most important factors separating the annotated draft genomes (Fig. 2). Genome analysis also confirmed the presence of certain metabolic genes that connect to the metabolic fingerprinting. For example, the genome of all four isolates contained genes needed to utilize sucrose (glycosidases, starch and sucrose metabolism KEGG pathway) and all four isolates were capable of doing so. Further, annotated genes for arabinose (L-arabinose isomerase and L-ribulose 5-phosphate 4-epimerase, pentose and glucuronate interconversions KEGG pathway) and xylose metabolism (xylose isomerase, pentose and glucuronate interconversions KEGG pathway) were only found for Cr-K1W and Cr-K20 and these were the only two isolates capable of utilizing those substrates.

Bacterial chromate stress has been shown to impact broad metabolic functions. For example, a proteomics study of Arthobacter FB 24 demonstrated that, when the bacterium was exposed to chromate, a large number of proteins involved in metabolism (e.g., carbohydrate, energy production, and amino acid transport and metabolism) had altered expression (Henne et al., 2009b). In S. oneidensis MR-1, chromate stress response caused variation in iron and sulfur transport genes and also DNA repair genes (Brown et al., 2006). Few genomic differences that related to sulfate transport are evident in the examined Microbacterium genomes. The faster reducers, K-29 and K-32, did have two additional genes that are related to iron transport (ABC-type Fe³⁺ transport system, permease component, and siderophore synthetase component). Thus these may be related to their increased ability to thrive in high chromate conditions.

Genomic annotations related to chromate reduction and resistance

The genomes of the four isolates examined contained putative chromate reductases. An automated IMG annotation of the isolate draft genomes failed to identify chromate reductases or resistant genes. Consequently, custom searches were performed to further examine the genomes for chromate related genes. Manual BLASTX searches were performed using the UniprotKB databases for ChrA, ChrR, and YieF. The search with known chromate transporters, ChrA homologues, returned no hits within any of the four isolate assembled genomes. While past studies have found that the efflux pump gene, chrA, can provide chromate resistance (e.g., Diaz-Perez et al., 2007; Henne et al., 2009a; Ramirez-Diaz et al., 2008), the lack of this gene could suggest that other efflux pumps may be exporting chromate. The genome of each of the isolates did contain two metal associated efflux pumps (arsentite efflux pump and cation/multidrug efflux pump). Alternatively, the physiological link between reduction and resistance could suggest that Cr(VI) reduction is responsible for providing resistance to these isolates.

Curated genome searches were able to identify multiple putative chromate reductases. BLASTX searches against the ChrR UniprotKB database returned a hit in all four isolate genomes. The identified genes were homologues (50–51% identity) to a ChrR reductase in Thermus scotoductus (Opperman, Piater & Van Heerden, 2008). IMG classified these genes as oxidoreductases and part of the Old Yellow Enzyme family. Analysis of Cr-K32’s protein (microk32_01893, 2529442617) using IMG’s homolog toolkit, showed a 100% homology to the protein found in Cr-K29 (microk29_00100, 2529452261) and 89% in Cr-K20 (microk20_02158, 2529450446) but only 35% percent identity to the protein found in Cr-K1W. Since the sequence homology of the protein found in Cr-K1W was low, additional searches within this genome were conducted. The genome of Cr-K1W did have a protein identified as an NADH:FMN oxidoreductase from the Old Yellow Enzyme family (cellulok1w_02457, 2529446915) that shared 100% homology to the chrR-like gene found in Cr-K20. Comparison of the four closest homologs’ gene neighborhoods within IMG showed high similarity of gene composition near the oxidoreductase despite the varying percent identity (Fig. S4). Further, jackhmmer searches (http://www.ebi.ac.uk/Tools/hmmer/search/hmmsearch), of the aligned genes found 51% protein identity, 68% similarity, and with an e-value of 3E–100 to Thermus scotaductus’ chromate reductase.

A similar curated search with YieF (ChrR sequence homologue) was also able to identify putative chromate reductases. The top result from a BLASTX search using the YieF databases identified a predicted flavoprotein in each of the individual isolates. Alignment of the four individual isolate genes within IMG’s homology toolkit revealed high degrees of sequence similarity with Cr-K1W (cellulok1w_02457, 2529446915) and Cr-K20 (microk20_03386, 2529451669) having 100% homology and Cr-K29 (microk29_00694, 2529452855) and Cr-K32 (microk32_00074, 2529440805) having 98% homology to one another. In addition, the genes were also homologues (65% for Cr-K29 and Cr-K32 and 69% for Cr-K1W and Cr-K20) to a predicted chromate reductase (ARUE_c41610) in an Arthrobacter sp. (Niewerth et al., 2012). Gene neighborhoods of the resultant genes did not show as much conservation compared the putative ChrR-like genes, with only small regions around the gene being shared between the fastest reducers (Cr-K32 and Cr-K29) and slower reducers (Cr-K1W and Cr-K20) but no shared genes between all four (Fig. S5).

Automated annotations of the isolate genomes failed to provide any putative chromate reductases, however, custom database searches were more successful. Each genome did contain genes with sequence homology to the chromate reductases, chrR and yieF, of non-model organisms. The putative chrR-like genes found in the isolates are homologous to a chrR gene (GenBank accession number AM902709) found in T. scotoductus, which has been experimentally shown to reduce chromate (Opperman, Piater & Van Heerden, 2008). The putative yieF genes found in the isolates were homologous to an annotated chromate reductase in Arthrobacter sp. RUE61a, a known chromate reducer. Despite being found in all four isolates, the genes did showed sequence variability between the four isolates. Isolates Cr-K29 and Cr-K32 are the faster reducers and their putative genes share more sequence homology when compared to Cr-K1W and Cr-K20, the slower reducers. Phylogenetic trees of both the putative ChrR and YieF genes also displayed this variability with slow reducers or fast reducers being more closely related to one another (Figs. S6 and S7). Neither putative gene was found to group with known chromate reducers, however, this may be an artifact from the lack of known chromate specific genes in Microbacterium. While the reasons for the variable ability for each isolate to reduce chromate are still not well understood, the differences observed between the isolates may be related to sequence similarity, expression, or genome content.