Comparative genomic analysis of a new tellurite-resistant Psychrobacter strain isolated from the Antarctic Peninsula

Strain isolation and culture conditions

P. glacincola BNF20 was isolated from a sediment sample collected at King George Island, Antarctica (S62°11′37.6″; W58°56′14.9″) during the ECA-48 Chilean Antarctic Expedition (January 2012). Bacteria were grown at 25 °C as described previously (Arenas et al., 2014) in Lysogenic Broth (LB) medium (Sambrook & Russell, 2001) supplemented with tellurite (200 µg/ml). Strains were identified by sequencing the 16S rRNA gene (accession MF806171) and determining the fatty acid profile. The 16S rRNA gene was sequenced at Pontificia Universidad Católica de Chile using Sanger sequencing with the primers 8F (5′-AGAGTTTGATCCTGGCTCAG-3′) (Turner et al., 1999) and 1492R (5′-ACGGCTACCTTGTTACGACTT-3′) (Lane, 1991). Fatty acid analyses were carried out at DSMZ, Braunschweig, Germany (Kämpfer & Kroppenstedt, 1996). The strain was deposited at DMSZ (Germany), accession # 102806.

To determine the minimal inhibitory concentration (MIC), bacteria were grown overnight in LB medium with shaking at 25 °C. Subsequently, saturated cultures were diluted 1:100 with fresh medium and grown to OD₆₀₀ ∼ 0.4–0.5. Then, 10 µl were added to 990 µl of LB medium containing serial dilutions of defined toxicants in 48-well culture plates. The plates were incubated with constant shaking for 24 h at 25 °C. Assayed toxicants included K₂TeO₃, K₂CrO₄, CdCl₂, ZnCl₂, CuSO₄, HAuCl₄, AgNO₃, NiSO₄, NaAsO₂ and Na₂HAsO₄.

16S rRNA gene phylogenetic analysis

A phylogenetic tree of P. glacincola BNF20—based on the partial 16S rRNA gene sequence—was constructed with bootstrap values based on 1,000 replications (Felsenstein, 1985). The nearly complete 16S rRNA gene sequence (1,516 nt) was obtained by merging the PCR sequenced amplicons (accession MF806171) and the sequence obtained by whole genome shotgun sequencing (accession AMK37_RS07000). Sequence alignments, assembly and comparisons, along with best model calculation and construction of the phylogenetic tree were carried out using the MEGA software version 6.0 (Tamura et al., 2013). As outcome, Jukes Cantor model and Pairwise Deletions for gaps treatment was the best fitting model for these sequence data. Nucleotide sequence positions from 16 to 1,535 were considered, according to the E. coli K-12 16S rRNA gene sequence numbering (accession AP012306). Scale bar represents 0.01 substitutions per-nucleotide positions. Moraxella osioensis DSM 6998^T (JN175341) was used as outgroup. The following 16S rRNA sequences from Psychrobacter strains were collected from GenBank (accession numbers are given in parentheses): P. glacincola DSM 12194T (AJ312213); P. adeliensis DSM 15333T (AJ539105); P. urotivorans DSM 14009T (AJ609555); P. arcticus DSM 17307T (AY444822); P. cibarius DSM 16327T (AY639871); P. cryohalolentis DSM 17306T (AY660685); P. frigidicola DSM 12411T (AJ609556); P. fozii NF23T (AJ430827); P. inmobilis DSM 7229T (U39399); P. namhaensis DSM 16330T (AY722805); P. aquimaris DSM 16329T (AY722804); P. luti NF11T (AJ430828); P. alimentarius DSM 16065T (AY513645); P. fulvigenes JCM 15525 (AB438958); P. piscatorii JCM 15603 (AB453700); P. jeotgalli JCM 11463T (AF441201); P. arenosus DSM 15389T (AJ609273) and P. okhotskensis JCM 11840 (AB094794).

Preparation of genomic DNA

P. glacincola BNF20 was grown in LB medium at 25 °C for 24 h with constant shaking. DNA was extracted using the Wizard Genomic® DNA Purification Kit (Promega, Madison, WI, USA). The quality and integrity of gDNA was determined by agarose gel (1%) electrophoresis and by determining the 260/280 nm absorbance ratio in a microplate multireader equipment (Tecan Infinite®PRO; Tecan, Männedorf, Switzerland).

Genome sequencing and annotation

The draft genome sequence of P. glacincola BNF20 was determined by a whole-genome shotgun strategy using the Illumina HiSeq 2000 platform with a mate-pair library of 3 kb (Macrogen®). A total of 10.89 million reads were quality filtered and assembled using the A5 pipeline (Tritt et al., 2012). Open reading frame prediction and annotation was carried out using standard operational procedures (Tanenbaum et al., 2010). Gene models were predicted using Glimmer 3.02 (Salzberg et al., 1998) and predicted coding sequences were annotated by comparison with public databases (COG, PFAM, TIGRFAM, UNIPROT, and NR-NCBI). P. glacincola BNF20 predicted proteome completeness was assessed by the presence/absence of bacterial orthologs according to the OrthoDB database using BUSCO (Simão et al., 2015). The circular genome map was assembled from P. glacincola BNF20 GenBank formatted file (NZ_LIQB00000000.1) using the plotMyGBK wrapper script (https://github.com/microgenomics/plotMyGBK); plotMyGBK uses BioPython and the R platform with the packages rSamTools, OmicCircos, and data.table to produce a vector image of a circular map (Cock et al., 2009; R Development Core Team, 2011; Morgan et al., 2016; Hu et al., 2014; https://github.com/Rdatatable/data.table).

Nucleotide sequence accession and culture collection number

Raw sequence data from P. glacincola BNF20 are available online under the BioProject #PRJNA293364, and Gold ID Gp0145575. The genome project has been deposited at GenBank, accession number NZ_LIQB00000000. The strain was deposited at the DSMZ culture collection, ID number DSM 102806.

Psychrobacter genome dataset

A total of 35 Psychrobacter genomes including P. glacincola BNF20 were retrieved from NCBI’s Genome and JGI GOLD databases (as of February 2017), where 10 and 24 genomes were annotated to the species and genus level, respectively.

Phylogenetic relationships and whole-genome nucleotide identity

The average nucleotide identity (ANI) was calculated for the 35-genome dataset using the pyani Python3 module (Pritchard et al., 2016) and the results were visualized using the data.table and pheatmap R packages (https://github.com/Rdatatable/data.table; https://cran.r-project.org/web/packages/pheatmap/index.html). Thirty-one phylogenetic marker genes corresponding to widespread housekeeping genes dnaG, nusA, rplA, rplD, rplK, rplN, rplT, rpsB, rpsI, rpsM, tsf, frr, pgk, rplB, rplE, rplL, rplP, rpmA, rpsC, rpsJ, rpsS infC, pyrG, rplC, rplF, rplM, rplS, rpoB, rpsE, rpsK, and smpB were identified in each Psychrobacter genome (AMPHORA2; Wu & Eisen, 2008). Each gene was translated under standard genetic code to perform a protein-coding-guided multiple nucleotide sequence alignment, using TranslatorX MUSCLE for the multiple sequence alignment (Abascal, Zardoya & Telford, 2010; Edgar, 2004). Alignments were concatenated using the alignment editor tool Seqotron (Fourment & Holmes, 2016) and the best partition scheme and substitution model was evaluated by PartitionFinder2 (Lanfear et al., 2016). Finally, the software MrBayes v3.2 was used for phylogenetic reconstruction (Ronquist et al., 2012), and the resulting tree was plotted and annotated using FigTree v1.4.3 (http://tree.bio.ed.ac.uk/software/figtree/). Phylogenetic tree annotation was based on the geographic location, according to BioSample database information for each Psychrobacter genome.

Search for metal resistance ortholog genes

To identify metal resistance genes, especially ter genes, a bidirectional Blast analysis was performed using the CRB-BLAST method (https://github.com/cboursnell/crb-blast). The BacMet Metal Resistance database (Pal et al., 2014) was used as target and the 35-genome dataset as query. In addition, each genome was re-annotated using the same methodology to identify syntenic genes based on BacMet and Prokka annotation without the bias of different annotation labels as implemented in Prokka v1.12 (Seemann, 2014). Finally, the results were visualized in their genomic context using the in-house script multiGenomicContext (https://github.com/Sanrrone/multiGenomicContext).

Promoter search

To elucidate if ter genes were under the control of one or more promoters, two promoter prediction tools were used on specific contigs where ter genes were found (LIQB01000002.1, position 189803-204267; LIQB01000003.2, position 286133-301767): PromPredict algorithm and the online program BPROM (Rangannan & Bansal, 2010; http://www.softberry.com/berry.phtml?topic=bprom&group=programs&subgroup=gfindb). Both results were jointly analyzed.

Taxonomic classification of tergenes

To investigate if the presence and number of ter genes was restricted to certain taxonomic levels, we downloaded all the bacterial reference genomes from the NCBI’s RefSeq database (January 2017; ftp://ftp.ncbi.nlm.nih.gov/refseq/release/bacteria/) and performed a bidirectional Blast searches (crb-blast) against the protein sequences of all ter genes in the BactMet database. Then, the in-house script fetchMyLineage (https://github.com/Sanrrone/fetchMyLineage) was employed to obtain the complete lineage of each bacterial genome with at least one ter gene match. The results were finally visualized using the R packages: ggplot2, RColorBrewer, devtools, ggjoy purr andreshape2 packages (R Development Core Team, 2011; http://ggplot2.org; https://cran.r-project.org/web/packages/RColorBrewer/index.html; https://github.com/hadley/devtools; https://cran.r-project.org/web/packages/ggjoy/index.html; Wickham, 2007).

A new Psychrobacter species from the Chilean Antarctic territory

P. glacincola BNF20 was isolated from Antarctic sediments and is a Gram-negative, non-motile, aerobic, oxidase positive, rod-shaped bacterium with an average dimension of 1.66 µm length and 1.09 µm width (Table 1, Fig. 1A). The fatty acid composition was determined at the DSMZ Institute (Germany) and showed that the major fatty acid was cis-9 octadecenoic acid C_18:1 ω9c (63.78%). The morphology description and major fatty acid component agrees with previous studies of Antarctic Psychrobacter isolates (Bozal et al., 2003). Initially, BNF20 was erroneously identified as P. inmobilis, based on a partial 16S rRNA gene sequence (Arenas et al., 2014). However, re-sequencing and a phylogenetic analysis of the partial 16S rRNA gene revealed that it is related to the P. glacincola species, family Pseudomonadaceae from the Gammaproteobacteria class (Fig. 1B). Altogether, morphology (electron microscopy), biochemical properties, partial (Sanger) and full length (NGS) 16 S rRNA gene sequence analysis, and fatty acid composition suggest that isolate BNF20 is member of the P. glacincola species.

P. glacincola BNF20 tolerates high tellurite and chromate concentrations

Several tests were carried out to determine if BNF20 was resistant to multiple metals. Besides tellurite (used in the initial selection), P. glacincola BNF20 was 4 times more resistant to chromate than the sensitive strain E. coli BW25113 under optimal growth conditions (Table 2). However, P. glacincola BNF20 growth was impaired in the presence of all other metal(loid)s tested, including Cu²⁺, Cd²⁺, Hg²⁺, Zn²⁺, ${\mathrm{AuCl}}_4^{1-}$, Ni²⁺, ${\mathrm{AsO}}_4^{2-}$, ${\mathrm{AsO}}_3^{1-}$, and Ag¹⁺.

First draft genome of Psychrobacter glacincola BNF20

Previous studies showed that P. glacincola BNF20 was highly resistant to tellurite (MIC ∼2.3 mM, Arenas et al., 2014). Although tellurite reduction is often accompanied by the formation of black deposits of elemental tellurium in resistant organisms, this phenotype was not observed in P. glacincola BNF20. To further investigate the mechanism(s) of tellurite resistance in P. glacincola BNF20, we sequenced the whole genome in search for genetic determinants implicated in metal(loid) resistance. The assembled genome of P. glacincola BNF20 consisted of 3,490,622 bp, 18 scaffolds, with an average G + C content of 42.76% (Fig. 1C; NCBI Reference Sequence: NZ_LIQB00000000.1). The predicted proteome scored 100% completeness according to the presence of highly conserved ortholog genes in bacteria (BUSCO analysis). A set of 47 tRNA genes and one copy of the rRNA operon were also identified. From a total of 2,968 predicted CDS, 2,872 (96.7%) ORFs matched coding sequences available in public databases, of which 2,515 were assigned (84.7%) or not (352 CDS, 19.31%) to COG categories (NZ_LIQB00000000.1).

P. glacincola BNF20 is evolutionarily divergent from other Antarctic Psychrobacterisolates

The genome sequence of P. glacincola BNF20 was compared to other 34 available genomic sequences by estimating ANI values and performing a multi-locus phylogenetic analysis (Fig. 2). Besides P. glacincola BNF20, the full dataset was composed of 10 named and 24 unnamed Psychrobacter species, respectively. P. glacincola BNF20 exhibited an average nucleotide identity >95% and an alignment fraction of over 80% with 3 isolates designated as Psychrobacter sp. JCM18903 (GCA_000586475.1), Psychrobacter sp. JCM 18902 (GCA_000586455.1) (Kudo et al., 2014) and Psychrobacter sp. P11F6 (GCA_001435295.1) (Moghadam et al., 2016), of which none was isolated from Antarctica. We did not find any genome comparison against BNF20 of >96.5% ANI and >60% alignment fraction, which has been suggested as a “genomic boundary” for bacterial species (Fig. 2A; Varghese et al., 2015). While some of the available genomes come from Antarctic isolates, none of them showed high ANI values (>90%): P. aquaticus (85%; GCA_000471625.1); P. alimentarius (85%; GCA_001606025.1); P. urativorans (85%; GCA_001298525.1); TB15 (84%, GCA_000511655.1), G (86%, GCA_000418305.1); PAMC 21119 (86%, GCA_000247495.2); TB2 (84%, GCA_000508345.1); TB47 (86%, GCA_000511045.1); TB67 (86%, GCA_000511065.1) and AC24 (86%, GCA_000511635.1).

Supporting our previous results, multi-locus phylogenetic analysis showed that P. glacincola BNF20 is more related to P11F6 (isolated from Tunicate ascidians from the Arctic, Moghadam et al., 2016), JCM 18902 and 18903 (isolated from frozen porpoise Neophocaena phocaenoides, Kudo et al., 2014). Antarctic isolates PAMC21119 and G (from King George Island, Moghadam et al., 2016; Che et al., 2013) belong to a polyphyletic group and do not form a monophyletic clade with P. glacincola BNF20, highlighting the heterogeneous nature of the Psychrobacter genus (Fig. 2B). All nodes of the phylogeny were well supported (posterior probability >  0.99).

P. glacincola BNF20 encodes multiple metal resistance determinants

As P. glacincola BNF20 was isolated from King George Island sediments, a place where heavy metal contamination has not been previously reported, we searched for genes known to be involved in metal resistance that could explain the observed tellurite and chromate resistance of strain BNF20 (BacMet database; Pal et al., 2014). Type and gene copy number distribution was not uniform in the 35-genome Psychrobacter dataset (Table S3).

Specifically, ∼100 genes possibly conferring metal resistance were identified in the genome of P. glacincola BNF20, of which some are related to chromate resistance, including chrL (BAC0361; regulatory protein, involved in chromate resistance), chrR (BAC0538; chromate reductase), mdrL/yfmO (BAC0209; multidrug efflux protein yfmO) and ruvB (BAC0355; ATP-dependent DNA helicase), and some to tellurite resistance—the so-called ter genes (Whelan & Colleran, 1992), including terA (BAC0386), terC (BAC0388), terD (BAC0389), terE (BAC0390) and terZ (BAC0392) (Fig. 3, Table S3). Two other genes apparently involved in tellurite resistance, ruvB (BAC0355; ATP-dependent DNA helicase) and pitA (BAC0312; low-affinity inorganic phosphate transporter 1), were also identified.

Organization of ter genes in P. glacincola BNF20

Given that (i) tellurite is by far more toxic for bacteria than other metals (Taylor, 1999) and (ii) it is scarce in the Earth’s crust (Turner, Borghese & Zannoni, 2012), finding tellurite resistance determinants in P. glacincola BNF20 was somewhat unexpected. Since to date the presence of ter genes in Antarctic microorganisms has not been reported, we focused the following analyses our study on them.

The ter genes were originally described as part of an E. coli operon exhibiting the terZABCDE structure (Taylor et al., 2002). P. glacincola BNF20 harbors terA, terZ, terE, terC and terD orthologs, but not terB (Fig. 3A); terA shows the opposite transcriptional orientation than the rest of the ter genes, while terZ is duplicated and is contained in different contigs (Fig. 3B). In addition, the expression of all ter genes in P. glacincola BNF20 seems to be regulated by individual promoters (PromPredict and BPROM analyses), suggesting that they are organized as a gene cluster rather than as an operon.

Three members of the Psychrobacter genus contained one ter gene (P. phenylpyruvicus (terZ, GCA_000685805.1), P. lutiphocae (terZ, GCA_000382145.1) and P. sp. ENNN9 III (terD, GCA001462175.1)), while the rest had different combinations of them (Fig. S1).

In P. glacincola BNF20, the context of the ter gene cluster is similar to other isolates like Psychrobacter sp. JC18902 (GCA_00058655.1), Psychrobacter sp. G (GCA_000418305.1), Psychrobacter sp. TB67 (GCA_000511065.1), Psychrobacter sp. AC24 (GCA_000511635.1), Psychrobacter sp. TB47 (GCA_00051045.1) and P. arcticus 273-4 (GCA_000012305.1). Interestingly, in all analyzed Psychrobacter genomes the ter gene cluster also contains a gene encoding a protein of the TIGR00266 family (unknown function, Fig. S1).

ter genes are distributed over several bacterial Phyla

To determine the frequency of ter genes in known bacterial genomes, their taxonomic distribution was evaluated. In general, ter genes are more commonly found in Gram-positive than in Gram-negative bacteria. Using NCBI’s RefSeq database (>5,000 genomes; accessed January 2017), we found that 48.59% of them contained ter genes (26 out of 30 bacterial Phyla). While, at the genus level, most genomes had one ter gene (67.95%) (Table S4, Fig. 4), others harbor two (2.31%), three (0.69%), four (5.24%), six (4.61%) or seven (1.15%) ter genes. Interestingly, the second most abundant combination of ter genes in genomes was five (18.04%), which could suggest evolutionary constrains.

At the phylum level most Proteobacteria contain one ter gene, with a few exceptions showing up to 7, including Yersiniacee, Morganellaceae, Enterobacteriaceae and Erwiniaceae. A similar pattern is observed in other Phyla, except for Firmicutes where genomes exhibit a defined array of ter genes (Fig. 4). Interestingly, while members belonging to the best represented family in RefSeq, i.e., Streptomycetaceae (149 genomes) exhibit five or six ter genes, in other well-represented families such as Flavobactericidae only 26 out of 114 genomes exhibit five ter genes (23%). Within the Moraxelaceae family, nine out of 45 genomes show five ter genes (20%, including BNF20), which agrees with the complete family database distribution (∼18% with 5 ter genes).