Genome sequence of Blochmannia pennsylvanicus indicates parallel evolutionary trends among bacterial mutualists of insects

Genome features

The B. pennsylvanicus genome consists of a 791,654-bp circular chromosome that we sequenced to 12× coverage (Table 1; Fig. 1; Supplemental Fig. S1; GenBank entry CP000016). This AT-rich genome has a relatively low percentage of coding DNA (76.7%) and shows a GC skew that distinguishes the leading and the lagging strands of replication [(G–C)/(G+C) averaged across 2.5-kb sliding windows]. We identified predicted ribosomal binding sites (RBSs) for 479 of 610 open reading frames (ORFs), but only 42 had the canonical AGGAG motif. Surprisingly, two pseudogenes (rpmD and uvrD) had detectable RBSs.

Differential gene loss, yet complete stability of genome architecture within Blochmannia

The 86-kb size difference between the Blochmannia genomes (B. pennsylvanicus, 792 kb; B. floridanus, 706 kb) largely reflects differential gene loss between the two lineages, with a greater extent of loss in B. floridanus. Both genomes possess several intact genes that are missing or pseudogenes in the other genome (Fig. 2). Although lateral gene transfer could, in principle, account for genes distinct to one Blochmannia genome, all genes grouped closely with enterobacterial homologs and showed amino acid and codon usage patterns typical of their resident genome. Notably, the order and strand orientation (i.e., genome architecture) is completely conserved for their 635 shared genes, including cases in which a functional gene in one genome is a pseudogene in the second.

Assuming that the common ancestor of the two Blochmannia strains encoded at least their combined set of 615 ORFs, then gene loss or inactivation in the lineage leading to B. floridanus occurred at an approximate rate of one ORF per 0.64–0.8 Myr [(25 ORFs lost)/(16 Myr) to (25 ORFs lost)/(20 Myr)]. Gene loss in the B. pennsylvanicus lineage is apparently ∼6.5 times slower, with loss or inactivation of one ORF per 4.0–5.0 Myr [(4 ORFs lost)/(16 Myr) to (4 ORFs lost)/(20 Myr)]. Estimated rates of gene loss in B. floridanus exceed those in Buchnera, in which one ORF was lost or inactivated per 2.70–3.60 Myr between Buchnera–B. pistaciae and Buchnera–Acyrthosiphon pisum [(111 ORFs lost)/(150 Myr × 2) to (111 ORFs lost)/(200 Myr × 2)] and per 1.70–2.38 Myr between Buchnera–A. pisum and Buchnera–Schizaphis graminum [(59 ORFs lost)/(50 Myr × 2) to (59 ORFs lost)/(70 Myr × 2)]. Rates of loss may be underestimated if the same genes were deleted independently along lineages. Mechanisms that influence rates of gene loss may include changes in the strength and/or efficacy of selection to maintain genes, as well as differences in underlying rates of gene knockouts and deletion (see Discussion).

Truncations and frameshifts in otherwise conserved Blochmannia genes

In all, 13 Blochmannia genes show a significant (20%–40%) truncation compared to Escherichia coli orthologs but apparently encode functional proteins. These genes include aceF, aroK, aroQ, ftsK, ftsY, hfq, mreC, pheA, rpoZ, thrS, trpD, yfcB, and yqeI. In all but one case (hfq, see below), truncations are shared by B. pennsylvanicus and B. floridanus and thus likely occurred before the divergence of these lineages. Certain truncations are also shared between Blochmannia and Wigglesworthia (yfcB, rpoZ, ftsK), Buchnera (aroQ), or among all three mutualist groups (aceF, ftsY).

Although a >20% length reduction is often interpreted as evidence for loss of gene function (e.g., Lerat and Ochman 2004), the sequence conservation of truncated Blochmannia genes suggests they encode functional proteins. First, nonsynonymous divergence (dN) between B. pennsylvanicus and B. floridanus is relatively low at the 12 truncated genes they share (average of 0.3012 ± 0.14), whereas synonymous divergence (dS) exceeds 2 for each gene. Although dN/dS is difficult to calculate because of saturation of synonymous sites, this ratio is below 0.13 for each gene. Moreover, apart from RpoZ, protein divergences at truncated genes are generally comparable to other ORFs (average protein divergence of 0.569 ± 0.281, compared to the genome-wide average of 0.565 ± 0.309). RpoZ is quite divergent (2.16), but the fact that dN/dS ≪ 1 suggests the truncation did not eliminate gene function. The one truncated gene (hfq) in B. pennsylvanicus that lacks an ortholog in B. floridanus has a relatively low protein divergence with the closest outgroup, Wigglesworthia (0.576). Thus, in the absence of detailed biochemical and structural information for each of these genes, a 20% truncation appears too conservative a criterion for inferring loss of function. In Blochmannia, truncated genes clearly differ from annotated pseudogenes, in which frameshift and nonsense mutations introduce multiple stop codons throughout the gene.

We also detected a single, short frameshift in each of B. pennsylvanicus ytfM, ybiS, hisH, and ubiF. The first two genes are unclassified proteins, while the latter are required for histidine and ubiquinone biosynthesis, respectively (see Fig. 2 legend). Each indel occurred within a 9–11-bp string of consecutive As (Fig. 3). Reanalysis of B. floridanus pseudogenes revealed a similar pattern at ytfM, ybiS, and gmhB, in which single frameshifts occur within poly(A) or poly(T) tracts. Poly(A) tracts are very common within Blochmannia ORFs (e.g., B. pennsylvanicus has 133 poly(A) or (T) tracts >9 bp long) and likely reflect AT mutational bias and reduced selective pressures. The indels in ytfM, ybiS, hisH, and ubiF each show high coverage in the B. pennsylvanicus shotgun assembly (5–15 independent clones). In the case of hisH, the frameshift occurs across independent PCR products sequenced directly with varying types of sequencing chemistry.

Despite their consistency, we hesitate to interpret these frameshifts as firm evidence of pseudogenes. Notably, apart from a single frameshift, these genes would otherwise encode intact proteins that are relatively conserved between the two Blochmannia genomes, with an average protein divergence (0.67 ± 0.16) that is only slightly above the average for other ORFs (0.565 ± 0.309). The occurrence of frameshifts within homopolymeric tracks is consistent with slippage during transcription (Wagner et al. 1990; Baranov et al. 2005) or during translation (Baranov et al. 2002; Gurvich et al. 2003) that could restore the full-length protein (see Discussion). If the frameshifts have, indeed, disrupted protein functions, these mutations must have occurred very recently since rapid mutation in Blochmannia (Degnan et al. 2004) is expected to erode pseudogenes quickly.

Metabolic similarities of Blochmannia spp.

Analysis of B. pennsylvanicus and reanalysis of other insect mutualist genomes using MultiFun (Serres and Riley 2000) allowed more comprehensive metabolic comparisons across groups (see Fig. 2 and Supplemental Tables S1 and S2). The two Blochmannia genomes share nearly all metabolic pathways thought to contribute to host nutrition, including most biosynthetic functions, sulfate assimilation and metabolism, and hydrolysis of urea (see Gil et al. 2003).

Metabolic differences between Blochmannia spp.

The 30 ORFs that distinguish the two Blochmannia genomes have a range of predicted cellular functions that may alter host–symbiont metabolic exchanges (Figs. 2 and 4). B. pennsylvanicus retains coaADE and dfp for the biosynthesis of coenzyme A, an essential cofactor and substrate for the TCA cycle (Fig. 4A). In contrast, the deletion of these genes from B. floridanus implies that this endosymbiont does not require coenzyme A, uses other enzymes for its synthesis, or imports this cofactor from the host.

Several differences between the Blochmannia genomes involve the biosynthesis, transport, and mediation of cellular wall and membrane components. B. pennsylvanicus retains six distinct ORFs that contribute to the de novo synthesis of peptidoglycan (murein), the major constituent of Gram-negative bacterial cell walls (Park 1996). In addition to MurI, it retains the complete pathway for the biosynthesis of isoprenoids (Fig. 4B), essential substrates for synthesis of peptidoglycan (El Ghachi et al. 2004), and several organic compounds including carotinoids, glycosyl carrier lipids, and the side chains of ubiquinone (Ogura et al. 1997; Kainou et al. 2001). In B. floridanus, the absence of murI and the interruption of isoprenoid biosynthesis implies this bacterium may import both D-glutamate and trans, trans-farnesyl diphosphate (isoprenoid precursor) from the host, might synthesize these components using other enzymes, or produce a cell wall lacking these structural elements (i.e., isoprenoids).

B. pennsylvanicus is the first fully sequenced insect mutualist that retains the entire sec-dependent secretory pathway, including the chaperonin SecB. This general pathway mediates the export and translocation of numerous proteins to the periplasm or inner membrane (Pugsley 1993). The other insect mutualist genomes lack certain components (typically secBDF) of this pathway (Shigenobu et al. 2000; Akman et al. 2002; Tamas et al. 2002; Gil et al. 2003; van Ham et al. 2003), although these losses are not expected to eliminate function (Mushegian and Koonin 1996). B. pennsylvanicus also retains the inner-membrane-bound heat-induced protease HtpX and the periplasmic chaperonins LolB and DsbA. Together, these results suggest that B. pennsylvanicus is better able to respond to cellular stress and ensure proper transport localization, conformation, and decomposition of gene products. Distinct membrane features in B. pennsylvanicus also include its retention of the inner membrane proteins YciC, putatively involved in transport, and BacA, which confers bacitracin resistance.

While insect mutualists have lost many regulatory genes, B. pennsylvanicus encodes three that are missing from B. floridanus (hns, hfq, and mntR). Like B. floridanus and Wigglesworthia, B. pennsylvanicus lacks DnaA, considered important in the initiation of DNA replication (Akman et al. 2002; Gil et al. 2003). In B. floridanus, the HU-like nucleoprotein HlpA might be involved in starting the nucleosome (Gil et al. 2003), and in B. pennsylvanicus, the nucleoprotein Hns might also be recruited to the replication origin for this purpose. Neither Blochmannia genome encodes the θ subunit of the holoenzyme DNA polymerase III, but B. pennsylvanicus retains the τ and γ subunits of DnaX. The τ subunit acts as a molecular tether that couples DnaB (DNA helicase) to the core of DNA polymerase III (α and ε) as the replication fork progresses from the origin (Walker et al. 2000; Gao and McHenry 2001). Thus, the inability of B. floridanus, Buchnera, and Wigglesworthia to transcribe the τ subunit is expected to decrease the efficiency, accuracy, and processivity of the holoenzyme (van Ham et al. 2003). These distinct features of the B. pennsylvanicus replication machinery might contribute to the apparent 29-kb shift observed in B. floridanus GC skew relative to B. pennsylvanicus (see Supplemental Fig. S1).

Long intergenic spacers in B. pennsylvanicus

Intergenic spacers in B. pennsylvanicus are significantly longer (average 291 bp) than homologous spacers in B. floridanus (average 180 bp; Wilcoxon Rank sum test, p < 0.0001), are longer than spacers in most other bacteria (Bacillus subtilis, average 121.9 bp; Vibrio cholerae, average 156.4 bp; Escherichia coli K12 (http://ecocyc.org), average 128.8 bp) (Mira et al. 2001) and contribute to its larger genome size compared to B. floridanus (77% vs. 84% coding regions) (Table 1). While spacers of the two Blochmannia genomes are too divergent to align, the strict conservation of gene order allowed us to predict homologous spacers based on their position in the chromosome. A detectable relationship between lengths of spacers in the two genomes suggests a certain degree of conservation of spacer length (Fig. 5A) (Spearman's ρ = 0.7895; p < 0.0001), although some short and medium-length spacers in B. floridanus are quite long (often >1 kb) in B. pennsylvanicus. Multiple gene-finding tools did not detect ORFs or recognizable pseudogenes in any spacers, which have lower GC contents (average of 20.0%) compared to ORFs (average of 32.1%) (Fig. 5B) and may represent eroded pseudogenes and/or regions involved in gene regulation.

Polymorphism within B. pennsylvanicus

Polymorphisms in the pooled symbiont population used for library construction were detectable as well-supported (Phred scores >40) discrepancies in the genome assembly that were represented by at least two independent clones. Nearly all (445/497) polymorphisms are single nucleotide polymorphisms (SNPs). Those located within ORFs occur primarily at third codon positions (Table 2). Although the majority of SNPs are located within ORFs, a disproportionate number of SNPs (35%) and indels (96% of insertions and/or deletions) occur within the intergenic regions (which comprise 23% of the B. pennsylvanicus genome). The two polymorphic indels within ORFs produce amino acid insertions/deletions in sucA and aroE.

Comparison of protein divergences

Wide variation in protein divergence across loci indicates variable functional constraint across Blochmannia proteins (Fig. 6; Supplemental Tables S3 and S5). Amino acid biosynthetic genes, ribosomal proteins, and particular chaperonins are exceptionally conserved (a mean protein divergence of 0.448, with a 99% confidence interval of 0.400–0.504), while surface structures [mean of 0.725 (0.533–0.991)] and hypothetical or unclassified proteins [mean of 0.775 (0.667–0.898)] are quite divergent (Supplemental Table S5). The membrane proteins tolA and tonB are particularly divergent (divergences of 1.74 and 2.64, respectively) (Supplemental Table S3).

Previous studies have demonstrated a negative relationship between GC content of endosymbiont genes and their level of divergence from free-living relatives (Herbeck et al. 2003; Banerjee et al. 2004), suggesting that amino acid changes in proteins under strong functional constraint are less severely affected by AT mutational bias. A strong negative association between GC content and protein divergence in Blochmannia (Supplemental Fig. S3) indicates this relationship also holds when protein divergences are estimated within a mutualist group (rather than to a more distant free-living relative). Protein divergences were, on average, ∼1.88 times faster in B. floridanus lineage compared to B. pennsylvanicus lineage. Genes with particularly elevated rates in B. floridanus include secG, rpsJ, rpsR, and rnpA, each of which evolves more than 10-fold faster in B. floridanus compared to B. pennsylvanicus (Fig. 7; Supplemental Table S4). A 14-fold rate acceleration at ispE may reflect the loss or disruption of other genes involved in isoprenoid biosynthesis in the B. floridanus genome and relaxed selection at this remaining isp gene.

The B. pennsylvanicus genome offered the first opportunity to compare genome-wide patterns of protein evolution in the context of distinct endosymbiotic associations. A previous study showed accelerated evolution at 16S rDNA and at synonymous positions of select Blochmannia genes compared to enteric bacteria and even compared to the rapidly evolving Buchnera (Degnan et al. 2004). Here, we tested whether the Blochmannia genome also undergoes exceptionally fast rates of protein evolution. Strong correlations in divergences at homologous genes indicate parallel functional constraints in Blochmannia compared to Buchnera and E. coli–Photorhabdus luminescens (Fig. 6). However, notably higher protein divergences between Blochmannia spp., despite their relatively recent split, reflects a several-fold acceleration in absolute rates of protein evolution. Among the 302 genes for which RSD detected orthologs across Blochmannia and Buchnera and estimated divergences within each mutualist pair, the average substitution rate in Blochmannia was 0.0132 to 0.0166 amino acid substitutions/site per million years (based on a 20- to 16-Myr divergence, respectively). This is ∼10 times faster than rates in Buchnera for the same gene set (0.00258–0.00362 amino acid substitutions/site per million years, based on a 70–50-Myr divergence). Previous studies have shown that nonsynonymous sites in Buchnera evolve an average of twofold (Canbäck et al. 2004) and up to 10-fold (Clark et al. 1999) faster than in the enterics. Likewise, we found an average ∼50-fold rate acceleration in Blochmannia compared to E. coli–Salmonella typhimurium across the 538 genes for which RSD detected orthologs (0.0142–0.0177 in Blochmannia, and 0.000245–0.000367 for E. coli–S. typhimurium assuming a divergence time of 150–100 Myr). Deviations from this general trend highlight genes that may experience different levels of functional constraint among genomes. For example, relative to other genes in the genome, amino acid biosynthetic genes evolve slightly slower in Blochmannia than in the enterics, while genes involved in translation (comprised largely of ribosomal proteins) evolve slightly faster than expected (Fig. 6).

Comparisons of average protein divergences within and between major functional categories confirmed many of these observations above (Supplemental Table S5). Blochmannia, Buchnera, and E. coli–P. luminescens show similar relative divergences across most functional categories, such as relatively high divergences of unclassified and hypothetical genes, loci for surface structures, cell membrane components; moderate divergences of genes encoding cofactor biosynthesis, information transfer, cell processes, metabolism, information transfer; and relatively low divergences of genes for nucleotide biosynthesis and amino acid biosynthesis. Striking differences among pairs include relatively high divergence of chaperonins and fatty acid biosynthetic genes in Buchnera, the relative conservation of regulation genes in both endosymbionts, and the slightly higher divergence of translation genes in Blochmannia, as noted above.

Blochmannia genome sequencing and assembly

B. pennsylvanicus genomic DNA (gDNA) was prepared from worker and larvae C. pennsylvanicus collected from five colonies at two sites in Falmouth, Massachusetts, USA. The gDNA was either extracted directly from the agarose plugs containing the purified bacterial cells (Charles and Ishikawa 1999) or gelpurified from a chromosomal fragment resolved through Pulsed Field Gel Electrophoresis (PFGE) (Wernegreen et al. 2002). Short (1.5–2.5 kb) insert libraries were generated from hydrosheared DNA using a double adaptor kit (SeqWright Inc.) (Andersson et al. 1996). Plasmid clones were purified and bidirectionally sequenced using BigDye v3.0 chemistry on either an ABI3700 or an ABI3730xl (Applied Biosystems). Detailed methods for library construction and sequencing are provided in the Supplemental material.

Raw sequence data were analyzed by PHRED (Ewing and Green 1998; Ewing et al. 1998) (http://www.phrap.org/phredphrapconsed.html) and screened using BLASTN/X (Altschul et al. 1990) for ant host contamination. Sequence reads that were putatively identified as γ-Proteobacterial (E ≤ 1^–10) were assembled using ARACHNE 2 (Jaffe et al. 2003) (http://www.genome.wi.mit.edu/wga). The resulting contigs were analyzed by hand in CONSED (Gordon et al. 1998) and using BAMBUS (Pop et al. 2004). The B. pennsylvanicus assembly was aligned using LAGAN (Brudno et al. 2003) (http://lagan.stanford.edu/lagan_web/index.shtml) to the published B. floridanus genome (NC005061), which facilitated primer design for gap closure by PCR.

Annotation and metabolic reconstruction

Open reading frames (ORFs) were identified iteratively using GLIMMER v2.10a (Delcher et al. 1999) (http://www.tigr.org/software/glimmer) and gene orthology predictions based on BLASTP sequence similarity to the NR, SWISS-PROT, and ECOLI databases; HMMR searches against Pfam_ls (Bateman et al. 2004); and identification of E. coli orthologs using the Reciprocal Sequence Distance (RSD) program (Wall et al. 2003) (details of the RSD method are noted below under “Comparison of protein divergences”). The few discrepancies among methods were limited to cases of differential loss of one gene from a pair of paralogs (ilvBG, tufAB, argFI), the presence of gene fusions (yidCD) or split genes (trpDG), or failure of RSD to identify a given ortholog because of high sequence divergence.

Three pseudogenes (uvrD, yqiC, rpmD) were detected as regions with similarity to functional ORFs in other genomes, but with multiple indels and missense mutations resulting in stop codons throughout each gene. The two most degraded pseudogenes (rpmD and yqiC) were undetectable by BLASTX and were only identified because of the conservation of gene order between the two Blochmannia genomes. Among the three pseudogenes, uvrD retains the longest (107 amino acids) intact reading frame with similarity to functional orthologs, but even this region is just 15% of the length of UvrD in E. coli and other outgroups. In this sense, annotated pseudogenes clearly differ from truncated B. pennsylvanicus ORFs, which retain at least 60% of the length of orthologous proteins. The three B. pennsylvanicus pseudogenes also differ from genes with single frameshifts within homopolymeric regions (hisH, ubiF, ytfM, ybiS), since the latter would encode intact, relatively conserved proteins if the frame-shift were corrected by slippage during transcription or translation (see Fig. 3, Results, and Discussion).

ORFS and RNAs were manually curated using a Generic Model Organism Database (GMOD) Web browser. ORFs that lacked sequence similarity to any entry in GenBank or the Comprehensive Microbial Resource (Peterson et al. 2001) and lacked any predicted protein domains in Pfam_ls were excluded from the annotation. Functional and pseudo-transfer RNAs (tRNAs) were identified using tRNAscan-SE (http://selab.wustl.edu/cgibin/selab.pl?mode=software), and ribosomal RNAs and structural RNAs were identified by BLASTN searches of the intergenic regions. Blochmannia gene functions and interactions were inferred from orthologs of E. coli K12 MG1655 described in GenProtEC (Serres et al. 2004) (http://genprotec.mbl.edu) and characterized by MultiFun (Serres and Riley 2000), two resources that represent functions of ∼80% of the 4401 genes in E. coli K12. Metabolic pathways were evaluated using the reference pathways available for E. coli at EcoCyc (Karp et al. 2004) and KEGG (Kanehisa and Goto 2000). Genomes of other insect mutualists were reanalyzed in the same manner for a consistent metabolic comparison.

Comparison of protein divergences

The RSD algorithm (Wall et al. 2003) was used to identify the reciprocal best BLAST hits (rbh) between translated ORFs of select bacterial genomes. The program used BLAST to identify potential matches of a given translated gene, aligned all potential matches using CLUSTALW (Thompson et al. 1994), and calculated a maximum likelihood estimation of amino acid substitutions between proteins using PAML (Yang 1997). Protein divergences were based on an empirical amino acid substitution rate matrix (Jones et al. 1992) and accounted for variation in evolutionary rates among protein sites using a γ distribution with shape parameter α = 1.53 (as recommended by Nei et al. 2001). The protein with the lowest divergence was then BLASTed against the first genome, followed by the alignment and divergence calculations. If the protein match with the lowest divergence was the same as the original query sequence, the pair was considered orthologous and the divergence was retained in the output.

Such comparisons were performed within endosymbiont groups: B. pennsylvanicus versus B. floridanus; Buchnera–A. pisum versus Buchnera–S. graminum; and Buchnera–A. pisum versus Buchnera–Baizongia pistaciae. Divergences in endosymbionts were compared to the enterobacterial pairs E. coli versus S. typhimurium, and E. coli versus P. luminescens. Genomes were downloaded from NCBI in June 2004. All genomes were compared to E. coli using RSD for ortholog detection and MultiFun-based functional assignments (detailed in Supplemental material). Rates of protein divergence along the lineages leading to B. pennsylvanicus and B. floridanus were compared using E. coli as an outgroup. Relative rates were calculated as follows, with “B0” representing the common ancestor of the two Blochmannia lineages:

Pairwise dN and dS values reported for certain genes [e.g., truncated genes or ORFs with apparent frameshifts within poly(A) tracts] were calculated using PAML (runmode –2) (Yang 1997). In order to account for the extreme base composition of endosymbiont sequences, we implemented a maximum likelihood model in which codon frequencies were calculated from the average nucleotide frequencies at the three codon positions (CodonFreq = 2). In these cases, translated Blochmannia genes were aligned to each other and to appropriate outgroups using CLUSTALX (Thompson et al. 1997) and back-translated to nucleotide sequences using the program RevTrans (Wernersson and Pedersen 2003).