Unusual Metabolism and Hypervariation in the Genome of a Gracilibacterium (BD1-5) from an Oil-Degrading Community

Genome assembly and curation reveal a hypervariable gene.The draft Gracilibacteria (BD1-5) genome binned by Hu et al. (6) from sample BD02T64, taken 64 days after the start of the laboratory experiment (see Materials and Methods), was selected for further curation as it comprised just 6 scaffolds. We verified that these scaffolds cluster tightly together on a tetranucleotide emergent self-organizing map (see Fig. S1 in the supplemental material), supporting their derivation from a single genome (8). Protein predictions for all of these six scaffolds required use of an alternative code in which the UGA stop codon is translated as an amino acid. Consistent with prior studies of Gracilibacteria, the genes were predicted using genetic code number 25 (UGA translated as glycine [9]). There have been two main ideas proposed to explain how alternative coding arises. The first relates to the low GC content of some (but not all) of the genomes it occurs in. The currently described genome fits this pattern (28.87% G+C). Alternatively, McCutcheon and Moran (10) invoke loss of peptide chain release factor 2 (encoded by prfB), which recognizes UGA codons, to explain the reassignment of stop UGA to tryptophan (code no. 4) in insect symbionts. Consistent with the hypothesis of McCutcheon and Moran (10), prfB was not detected in the genome of the gracilibacterium studied here, or in any other available BD1-5 genome. However, peptide chain release factor 1 (prfA) was detected, and the gene coding for it is widely identified across the CPR.

Prior to read-based curation, the six scaffolds were tentatively condensed into two based on perfect overlaps at scaffold ends. Local assembly errors were removed by curation, and unplaced paired reads were used to close gaps. Reads mapped to the scaffolds were visualized in Geneious (11). Notably, a region where two scaffolds were joined based on end overlaps was identified as incorrectly assembled based on the absence of perfect read support. Inaccurate read placements were associated with a hypervariable repeat locus. By manual step-by-step repositioning of paired reads (first placing reads anchored into the repeat locus boundaries, then considering paired read distances and repeat composition), it was possible to generate a representation of the sequence through the locus (Fig. 1). Due to the large size of locus compared to the paired read distance, it is impossible to determine the exact number of repeats or if the locus exhibits cell-to-cell variation in repeat number per locus. However, based on the average sequencing depth, the approximated locus is probably of about the correct length and not highly variable in terms of the number of repeated sequences.

• Open in new tab
• Download powerpoint

FIG 1

Repeat locus from the BD1-5 genome. Colored arrows represent repeated sequence blocks, the sequences for which are shown in the “Repeats” insert. Sets of arrows represent reads, and reads linked within this region to paired reads are indicated by a thin connecting line.

We verified the final genome path by calculating the cumulative GC skew of the closed chromosome sequence and identified the pattern expected for normal bacterial bidirectional replication (Fig. 2). The final assembly comprises 1.34 Mb, 1,243 protein coding genes, 33 tRNA genes and one set of rRNAs (Table 1). According to an RAxML tree based on 16S rRNA genes, the closest relative to our organism was sampled from deep sea sediments; other closely related sequences are from marine environments (Fig. 3).

• Open in new tab
• Download powerpoint

FIG 2

Diagram showing the GC skew (gray dots) and calculated cumulative GC skew (green line) across the finished BD1-5 genome. The pattern is typical of a correctly assembled genome of a bacterium that undergoes bidirectional replication from origin to terminus.

• Open in new tab
• Download powerpoint

FIG 3

Phylogenetic placement of the Gracilibacteria genome from sample BD02T64 reported here. The 16S rRNA tree was constructed using the maximum likelihood method RAxML. The small black circles indicate nodes with values of >70% bootstrap support. 16S rRNA genes retrieved from genomes are indicated by green circles. Dotted circles represent published draft genomes, and the full circle indicates the finished and curated genome from this study. Colored circles indicate the type of ecosystem from which sequences were obtained. The full tree file is provided in the Data Set S1.

The predicted amino acid sequence of the BD1-5 gene containing the repeat region is shown in Text S1, part A, in the supplemental material. Some repeat types occur in blocks, and some repeat types alternate, but overall the most striking feature of the locus is the high level of apparent cell-to-cell heterogeneity (Fig. 1). Variant calling in reads mapped to the full-length protein identified 17 synonymous single nucleotide substitutions and zero nonsynonymous substitutions, with the exception of instances occurring only on a single read. In fact, this gene contains the highest proportion of synonymous substitutions in any genes in the BD1-5 genome.

Within the repeat region, the incidence of synonymous versus nonsynonymous substitutions is shown in Fig. 1 (insert). The four main repeat nucleotide sequence variants are indicated in orange, green, yellow, and blue, along with their translated sequences. Single-incidence sequences are indicated by white bars. Notably, the nucleotide sequences of the four major repeat variations all translate to the tripeptide amino acid motif PTD. Given that the Gracilibacteria population cells share near-identical nucleotide sequences genome-wide, except within this specific locus, we infer that the repeat-bearing protein may be under strong pressure to evolve at the nucleotide level. If cells acquire nonsynonymous substitutions in the repeat protein, they are apparently strongly selected against.

The PTD repeat motif is found in hypothetical proteins and predicted surface proteins of a few other organisms, including some that are eukaryotic sporozoite surface protein 2-like. A secondary structure prediction of the BD1-5 protein suggests only β sheet and coils, with the repeat motif in a coil region. We predict a single N-terminal transmembrane domain and extracellular localization of the remainder of the protein sequence, including the repeat region (see Fig. S2 in the supplemental material).

We investigated codon usage in the repeat gene. The codons for D are GAC and GAT, with usage of 6.2:6.7 in the repeat gene, whereas the expected genome-wide incidence is GAC:GAT = 0.83:4.79. Synonymous substitutions within the repeats could cause ribosome pausing and modulate rates of protein folding (12, 13). While we considered that atypical codon use in this region might indicate selection for translational pausing, the corresponding 5′ tRNA anticodon position (G) enables a wobble pair to recruit the same tRNA to either GAC or GAT. The codons for P are coded for by CCA, CCG, CCT, and CCC, with an expected incidence of 1.3:0.13:0.89:0.13, whereas the repeat gene has an incidence of 11.82:0.25:1.23:0. Thus, there is evidence for strong selection for the CCA codon (considered further below). Inosine is the only tRNA-proline 5′ anticodon base that could accommodate all synonymous variants. The more prevalent codons, CCA and CCT would recruit two different tRNAs. Intriguingly the tRNA-proline gene carried by the gracilibacterium genome corresponds to CCA or CCG. In the third position of the tripeptide, T can be ACT, ACC, ACA, or ACG with an expected incidence of 2.1:0.41:2.4:0.36. However, within the repeat gene, ACT, ACC, ACA, and ACG occur with an incidence of 4.19:0:6.9:6.4. Again, with the exception of anticodon inosine wobble pairing, reliance on rare codons may indicate selection for translation pausing.

If it is advantageous for coexisting cells to have highly variable rates of translation, one might expect that the sequences would make maximal use of the available codons. Counter to this, we see reduced codon diversity. Thus, we considered that variation in the secondary structure of the RNA sequence in the repeat array may be selected for. In the secondary structure prediction for the repeat region, we note the periodic alternation of stems, comprising mostly Watson-Crick base pairs, and loops (see Fig. S3 in the supplemental material). Notably, the CCA codon (specifically the first C) is at the base of the bubbles and closes them, paired to G’s from either the first base of the first codon or the last base of the third codon. Stem-loops impact RNA folding, can stabilize mRNA, and provide recognition sites for RNA binding proteins. We speculate that nucleotide variation may impact the translation rate of this gene and lead to variation in the fitness of different population members.

We searched the genomic region flanking this gene but did not identify a known mechanism for site-directed mutagenesis within the repeat locus. Genes with functions linked to DNA repair and recombination are found in close downstream proximity (uvrC excinuclease [5,798 bp downstream], an exodeoxyribonuclease III gene [8,336 bp downstream], and DNA recombination-mediator protein gene, dprA [30,798 bp downstream]). Perhaps this organism possesses a DNA mutator, which mediates targeted diversification in the repeat locus. It is unlikely that the organism is deficient in repair enzymes, as sequence variation is not elevated elsewhere in the genome. Perhaps nucleotide heterogeneity arose due to suppressed proofreading in this region, but we have no explanation for how this might have occurred.

In the current study, it is difficult to evaluate locus length variation because read lengths are short compared to the length of the repeat arrays. Locus length variation is expected, given the presence of perfect repeat arrays. In bacterial genomes, repeat regions may expand and contract due to either recombination or slipped-strand mispairing (SSM [14, 15]), resulting in population variability in terms of tripeptide motifs that may impact three-dimensional protein structure and ligand binding. The relationship between microsatellite length and point mutation has been described elsewhere and generally predicts that as a locus expands, base substitutions accumulate and suppress further SSM (16, 17). If SSM is undesirable, it is advantageous to include nucleotide variants that offset repeat pairing and thus prevent slippage.

Examination of the 5′- and 3′-untranslated regions flanking the repeat gene uncovered two sequences capable of forming stem-loops with notably long stems (14 to 15 bp) and 4- to 6-bp loops (see Fig. S4 in the supplemental material). As DNA or RNA structures, these stem-loops may play a role in recombination or as transcriptional regulation signals for the repeat-containing gene, respectively.

Some reads mapped to the BD1-5 repeat region had paired reads that were not placed in that genome. Comparison of the non-repeat regions of these reads and the sequences of their unplaced paired reads to the genomes of other community members revealed 100% nucleotide matches to a region on BD02T64_scaffold_179, part of a draft Colwellia psychrerythraea genome (BD02T64_Colwellia_psychrerythraea_38_180_partial). Thus, we concluded that a region within the genome of this abundant population has the same PTD repeat as found in the Gracilibacteria protein (Text S1, part B). After curation of the region, the C. psychrerythraea protein is predicted to be 3,459 amino acids in length, with a signal peptide and extracellular localization, possible galactose/carbohydrate-binding domains, pectin lyase/virulence domains and parallel β helix repeats. The repeat occurs within a structure that otherwise consists of a mixture of α helices and β sheets but is in neither of these. Besides the repeat region, the C. psychrerythraea protein does not share any sequence identity (<10%) with the Gracilibacteria protein. Also, the TGA codon, which is repurposed in Gracilibacteria, is not used in either of the repeat regions.

Within the Colwellia protein, reads carry up to 11 repeats (see Fig. S5 in the supplemental material). As for BD1-5, it is impossible to detect variation in repeat number in each cell due to the read length limitation, but one read has only five repeats. In virtually all cases, the nucleotide repeat is encoded by a single 9-mer (yellow), This 9-mer is prominent toward the end of the Gracilibacteria repeat region. The loci in both genomes terminate with the same 9-mer (orange in Fig. S5). The essentially perfect repeated sequence would make this region prone to replication slippage, leading to cell-to-cell variation in the number of tripeptides in the protein.

Interestingly, several of the Gracilibacteria proteins encoded immediately adjacent to the variable PTD protein have the highest similarity to proteins in organisms that are not part of the CPR. One is most similar to a protein from Colwellia psychrerythraea, although the percentage of amino acid identity is low (∼53%). To rule out chimeric assembly of sequence from another bacterium in this genomic region, we confirmed the expected alternative coding throughout (and paired-read placements were verified during the main curation phase). Thus, the region encoding the Gracilibacteria variable repeat gene and adjacent genes may have been acquired from a bacterium related to Colwellia psychrerythraea.

Metabolic analysis.The biosynthetic pathways easily recognizable in the genome are for ribosome-based protein synthesis, nucleic acid synthesis and interconversion, DNA repair, peptidoglycan production, secretion, pilus production, and cell division. However, as for other members of the CPR, this gracilibacterium appears to lack the ability to synthesize lipids needed for construction of the cytoplasmic membrane (and, there is no pathway for synthesis of lipid A required for a Gram-negative cell envelope). Thus, these cells are predicted to be either symbionts or closely dependent on other community members for key building blocks. The genome lacks a CRISPR-Cas system for phage defense, but has a restriction modification system that may serve this purpose (18, 19). Absent are almost all pathways for amino acid synthesis, leading us to conclude that amino acids needed for protein biosynthesis are derived through breakdown of externally derived peptides. Many different types of peptidases and proteases are available for this process.

For nucleic acid synthesis the genome encodes the steps required to interconvert nucleotides. We also identified most of the genes required for biosynthesis of purines and pyrimidines from glutamine and aspartate; these genes are relatively uncommon in CPR. IMP can be converted to AMP, ADP, and ATP and incorporated into RNA and DNA. Enzymes were also identified to interconvert forms of GDP and GTP. Genes of the one carbon pool by the folate pathway were identified, enabling transfer of C1 groups during nucleotide metabolism, but genes for folate biosynthesis were not identified.

Perhaps the most surprising feature of this bacterium is the complete lack of genes for glycolysis and the pentose phosphate pathway, which makes this genome distinct from other Gracilibacteria, and possibly even from all other bacteria. At least partial pathways are present in other Gracilibacteria, and the first reported genomes have full pathways to convert glucose to pyruvate and fermentation-based metabolisms (3). More broadly, at least parts of these pathways are present in the most minimal CPR genomes. However, this is the first genome from a major subgroup within Gracilibacteria (Fig. 3), so it remains to be seen whether this is a common trait. The absence of these pathways raises two questions: (i) the nature of central carbon metabolism in these organisms and (ii) how ATP, NADH, NADPH, and ferredoxin are reduced and recycled.

Potentially addressing the first question, we identified a variety of pathways for production of central carbon currencies. We identified a putative two-subunit ATP citrate (pro-S)-lyase (EC:2.3.3.8) (genes 1051 and 1052), a complex rarely detected in CPR. This annotation (versus citrate synthase) was supported by HMM homology and the presence of the active site residues GHAGA (20). Via this complex, citrate can be converted to acetyl-CoA and oxaloacetate. Citrate may be obtained from external sources via two putative citrate transporters. Intriguingly, both ATP citrate (pro-S)-lyase subunits are most similar to predicted proteins in archaea, suggesting their acquisition via lateral gene transfer. We predict that oxaloacetate derived from breakdown of citrate is converted to pyruvate via a 2-oxoacid ferredoxin oxidoreductase (OFOR). Pyruvate can also be produced from phosphoenolpyruvate via a pyruvate kinase, and from malate (1.1.1.38) and serine (4.3.1.19). Overall, amino acids scavenged from the environment appear to feature prominently in the metabolism of this gracilibacterium, and some are converted into the nitrogen and carbon storage compound cyanophycin.

Addressing the second question, we identified many reactions that oxidize or reduce energy currencies via transformation of small carbon compounds. Specifically, pyruvate conversion to acetyl-CoA via OFOR consumes NADH while reducing ferredoxin. The ferredoxin may be reoxidized via either a cytoplasmic or membrane-bound ferredoxin reductase (FNR) that also converts NADP⁺ to NADPH. NADH may be regenerated in the production of pyruvate from serine or malate. Like citrate, malate may be obtained from external sources. Other reactions, such as those involved in peptidoglycan synthesis and interconversion of tetrahydrofolate compounds, also interconvert energy currencies. Enzymes that respond to oxidative stress response also may provide electron sinks.

Many CPR bacteria generate ATP via substrate-level phosphorylation reactions that produce compounds such as acetate, but genes for production of these short-chain fatty acids were not identified. ATP required for DNA and RNA biosynthesis may be formed via the F-ATP synthase complex (complex V). Given the lack of an electron transport chain that could pump protons, proton motive force (PMF) could be stolen from attached host cells if tight junctions are formed (19, 21). Such close physical associations have been reported for another CPR bacterial group, Saccharibacteria (TM7), which attach to host Actinobacteria cell surfaces (22). Alternatively, proton motive force could be generated by cytoplasmic drawdown of H⁺ via reactions involved in breakdown of amino acids and other compounds, Na⁺/H⁺ antiport, or consumption of H⁺ by superoxide dismutase. The ATP synthase may also be used reversibly to generate proton motive force (as suggested by Wrighton et al. [3]), but no complexes were identified that could make use of the generated PMF. Specifically, there is no indication of hydrogenases, which occur in some other CPR members. Lacking also are other electron transport chain components, such as NADH dehydrogenase, succinate dehydrogenase, and cytochrome c reductase/oxidase, and most steps of the tricarboxylic acid cycle (Fig. 4).

• Open in new tab
• Download powerpoint

FIG 4

Cell cartoon depicting a reconstruction of the metabolism of the gracilibacterium. Bold text indicates prominent functions, blue text indicates resources inferred to be externally derived. * indicates that reactions for biosynthesis of cofactors require a precursor compound. Abbreviations: PEP, phosphoenolpyruvate; UDP-GlcNAc, UDP-N-acetyl-α-d-glucosamine; OFOR, 2-oxoacid ferredoxin oxidoreductase; 3PG, 3-phospho-d-glycerate; 3-PoxyP, 3-phosphonooxypyruvate; 2-oxoglu, 2-oxoglutarate; FNR, ferredoxin reductase; PPPi, PPi, and Pi: phosphate compounds interconverted by inorganic pyrophosphatase; Mk-n, metaquinone; Mk-l, metaquinol; Succ. semiald., succinate semialdehyde; l-Glu, l-glutamate; R-COOH, a carboxylic acid; CPG, cyanophycin; FG, N-formyl-l-glutamate. PTD is a tripeptide repeat.

A variety of transporter types were predicted, presumably addressing the need to acquire compounds from other cells or detritus. There are many hypothetical membrane-associated proteins with multiple transmembrane (TM) domains that also may serve a transport role. Overall, 122 transmembrane proteins (>3 transmembrane domains) and 80 transporter proteins were identified (Table 1; see Table S1 in the supplemental material). The genome encodes an intriguing 990-aa protein predicted to contain 32 transmembrane domains (gene 860). A large-scale analysis of TM-rich proteins in the NCBI nr database revealed that very few have 32 or more TM domains, and only a few related proteins are known (mostly in other Gracilibacteria). The function of this enigmatic protein is uncertain as the only domain predicted is DUF2339 (hypothetical membrane protein).

A notable feature of the Gracilibacteria genome is the prominence of secretion mechanisms and secreted proteins. We identified 66 such proteins using a combination of three methods to predict signal peptide-mediated export, of which only 24 are shorter than 300 amino acids. A further 104 proteins are predicted to be secreted via nonclassical pathways that do not use a signal peptide. Of these, 43 are larger than 300 amino acids. In addition to a sortase (typically found in Gram-positive bacteria and common in CPR), we identified genes of the type II and IV secretion pathways that are generally associated with Gram-negative bacteria, including multiple copies of SecA, -D, -F, -Y, -E, and -G). SecYEG form the central translocase across the inner membrane, SecA guides proteins to the translocase channel and is the ATPase, and SecF promotes release of the mature peptide into the periplasm. Thus, the identified components provide the functions required for secretion in non-Gram-negative bacteria. Intriguingly, 15 general secretion protein G proteins (GspG, alternatively PulG) are predicted, as well as GspE. These are large proteins, on average 528 aa in length. GspG is the major pseudopilin present in a pseudopilus, and GspE is an ATPase involved the assembly of the pseudopili. In addition, we identified around 12 type IV pilus assembly protein subunits, some in multicopy. Type IV pili allow the transfer of genetic material representing PilV, -C, -B, and -W and are involved in twitching motility (the genome also has two pilT genes). PilD (leader peptidase) was also identified. We did not identify PilQ, consistent with lack of outer membrane. Pili may be involved in attachment and interorganism interactions, as well as uptake of DNA. Competence genes were also identified (19, 21, 23).

From the perspective of the cell envelope, the biosynthesis pathway for peptidoglycan is complete, although the requirement for precursor UDP-N-acetylglucosamine from external sources is predicted. Predicted are genes to convert phosphorylated isoprenoid into a precursor for peptidoglycan, but the genome lacks the archaeal mevalonate and bacterial MEP (2-C-methyl-d-erythritol 4-phosphate) pathways. It has geranylgeranyl diphosphate synthase, but the reason is unclear. In addition, we identified three genes that degrade l-lysine and d-glutamate that may feed intermediates into two different steps within the peptidoglycan biosynthesis pathway. The genome contains many genes for polysaccharide synthesis (e.g., no. 444-460) and for proteins with S-layer domains. Thus, we anticipate a cell-wall-containing peptidoglycan with a periodic surface layer, many and potentially diverse pili, and a variety of large extracellular proteins and polymeric substances (Fig. 4). Interestingly, some S-layer proteins may have toxin domains (e.g., 1226, predicted to have polycystin-1, lipoxygenase, and alpha-toxin domains). Other large proteins have annotations suggestive of hostile interactions with other organisms (e.g., insecticidal toxin complex protein [TccC]), and there is a predicted invasin domain in one large protein in the genome.

In terms of the ability to respond to environmental conditions, the genome encodes at least four RelA/SpoT domain proteins, three of them encoded sequentially and one larger multidomain protein encoded elsewhere. These may function in response to nutrient limitation. Also identified are two 8-oxo-dGTP diphosphatase genes to prevent misincorporation of the oxidized purine nucleoside triphosphates into DNA and proteins with antioxidant functions, including superoxide reductase and enzymes to reduce oxidized methionine.

We conclude that the inferred putative symbiotic lifestyle of Gracilibacteria differs in notable ways from those of other obligate host-associated organisms. The genome size is large, compared to those of most obligate host-associated organisms (usually <1 Mbp in length [24]). Host-associated bacteria that have experienced moderate genome reduction retain genes for synthesis of fatty acids and peptidoglycan (but not for lipopolysaccharide [LPS] or phospholipids), whereas those that have undergone extreme genome reduction have essentially no genes for cell envelope biosynthesis (10). In contrast, the gracilibacterium seems to rely entirely on externally derived fatty acids. It retains genes for regulation of gene expression (e.g., two-component systems and various transcriptional regulators), DNA repair, and homologous recombination, whereas these genes are often lost in symbionts (7). Overall, the genomic features of this gracilibacterium only overlap partially with those of host-associated bacteria, which have experienced rapid genome decay.