Genome evolution in an ancient bacteria-ant symbiosis: parallel gene loss among Blochmannia spanning the origin of the ant tribe Camponotini

Preparation of genomic DNA

A single colony of C. obliquus was collected near Morehead City, North Carolina, USA by B Guénard, and P. turneri was collected near Townsville, Australia by SKA Robson (see ‘Acknowledgments’). Voucher specimens were deposited in the Bohart Museum of Entomology, University of California, Davis (UCDC), corresponding to voucher ID numbers CASENT0221021 (C. obliquus) and CASENT0220426 (P. turneri). We used the Qiagen DNeasy Blood and Tissue Kit to prepare genomic DNA from a pooled sample of seven eggs, two larvae, five pupae, eight minor workers, six major workers and five female alates for C. obliquus and a pooled sample of three worker gasters for P. turneri.

Sequencing and assembly of Blochmannia genomes

C. obliquus gDNA was sequenced on an Illumina HiSeq (Illumina, San Diego, California, USA) to generate 100 bp paired end reads. We modified filter_reads.py (https://github.com/nickloman/xbase/blob/master/short-read-assembly/filter_reads.py) for Sanger FASTQ format and used the script to remove paired reads with any bases of quality score <30, which retained 6,184,892 read pairs. We assembled this filtered read dataset using Velvet v1.2.07 (Zerbino & Birney, 2008) with a hash length of 61, exp_cov 200, cov_cutoff 20 and scaffolding turned off. This generated 566 contigs, one of which aligned to B. pennsylvanicus using MAUVE (Darling, Mau & Perna, 2010). We observed an overlap of 116 bp between the contig ends, which suggested that de novo assembly had produced a closed genome. To test this, we used Mosaik (Lee et al., 2014) to align the filtered read dataset against a 1,260 bp sequence encompassing the joined contig ends and flanking regions. This alignment produced no zero coverage regions, which confirmed that the single contig was the closed B. obliquus genome.

To finish the genome, we used a two-step process applying different alignment programs to confirm the majority genotype, or the base represented by the majority of reads at each position. In the first step, we used Mosaik, which allows the user to set the stringency of mismatch tolerance. We aligned the filtered read dataset against the closed genome sequence with a maximum mismatch threshold of 12 and then removed duplicate read pairs with the DupSnoop module. Using Consed (Gordon, Abajian & Green, 1998), we generated a questionable consensus bases report and a highly discrepant indels report for the resulting alignment, which did not identify any positions that needed editing.

The second step in our finishing process invokes the IndelRealigner module of the Genome Analysis Toolkit (GATK) to ensure accurate identification of indels (DePristo et al., 2011). We aligned the filtered read dataset against the closed genome sequence using BWA (Li & Durbin, 2009) and then processed the resulting alignment with RealignerTargetCreator and IndelRealigner from the GATK package. We removed duplicate reads using Picard MarkDuplicates (http://broadinstitute.github.io/picard). Finally, we analyzed the processed alignment using VarScan (Koboldt et al., 2009) to identify positions at which the base in the closed genome sequence differed from the majority of aligned reads. We did not identify any such positions, confirming that the single contig was the closed, finished B. obliquus genome. Genome coverage for this alignment averaged 558x. This sequence is deposited in GenBank as accession number CP010049. The Illumina read dataset is deposited in the NCBI Short Read Archive (SRA) as SRP050154.

P. turneri gDNA was sequenced on an Illumina Genome Analyzer II (GAIIx; Illumina, San Diego, California, USA) to generate 150 bp paired end reads. We used DynamicTrim.pl and LengthSort.pl from the SolexaQA package (Cox, Peterson & Biggs, 2010) to generate trimmed reads of at least 80 bp with quality score >30 for each base. The resulting trimmed read dataset included 9,590,066 read pairs. We assembled this filtered read dataset using Velvet with a hash length of 41, exp_ cov 200, cov_cutoff 20 and scaffolding turned off. This generated 8,275 contigs, four of which aligned to B. pennsylvanicus using MAUVE. To close the four gaps, we generated Sanger sequencing reads and used Phred/Phrap/Consed to assemble and manually examine the sequence.

To finish the B. turneri genome, we followed the two-step process described above for B. obliquus. In the first step, we corrected nine positions based on alignment of the full read dataset against the closed genome sequence using Mosaik. In the second step, we aligned the full read dataset against the corrected sequence using BWA followed by processing with GATK. We analyzed the processed alignment with VarScan and corrected a single position. Genome coverage for this alignment averaged 1223x. The closed, finished genome sequence is deposited in GenBank as accession number CP010048. The Illumina read dataset is deposited in the NCBI Short Read Archive (SRA) as SRP050161.

Annotation of Blochmannia genomes

We used an annotation engine hosted by the Institute for Genome Sciences (IGS) at the University of Maryland School of Medicine to generate an automated annotation of each genome sequence (Galens et al., 2011), which we then manually curated within the MANATEE framework (http://manatee.sourceforge.net/igs/). Protein-coding genes predicted by the annotation engine were removed if they lacked a Blast-Extend-Repraze (BER) alignment score <10⁻⁵ to a protein-coding gene from outside of Blochmannia. We manually examined possible frameshifted genes flagged by the annotation engine. For genes with frameshifts in homopolymer tracts, we included the likely position of the frameshift in the GenBank annotation. Because the frameshifts may be corrected by polymerase slippage (Tamas et al., 2008; Wernegreen, Kauppinen & Degnan, 2010), we consider these to be intact genes.

We curated start sites using BER alignments to Blochmannia and closely related species. When possible, we used the gene name and symbol listed in SwissProt for the homologous gene in E. coli to maintain consistency with existing proteobacterial annotations. For conserved hypothetical proteins or proteins with similarity to a protein family but not a specific family member, we did not assign a gene name and refer to them using the locus tag (for example, BTURN675_020).

After curating the annotations, we analyzed intergenic regions in each genome with RFAM (Burge et al., 2013) and BLASTX (Altschul et al., 1990) to identify uncalled genes and pseudogenes. In both B. obliquus and B. turneri, RFAM identified three RNA-coding genes (ffs, rnpB and tmRNA). To identify protein-coding genes, we aligned intergenic regions to the GenBank non-redundant database using BLASTX with the low complexity filter turned off. We manually examined hits with e-value <10⁻⁵. Pseudogenes were identified by multiple nonsense mutations, frameshifts and/or gaps. Pseudogenes that aligned to intact homologs with >60% coverage had at least two nonsense mutations. We annotated pseudogene coordinates using the boundaries of the BLASTX alignments. Analysis of B. obliquus intergenic regions detected two genes (cyoD and sdhD) and 15 pseudogenes (dnaA, engD, glnA, pdxA, pdxB, pdxJ, ribA, ribB, ribC, ribD, secD, secF, topA, yigB and uvrD). Analysis of B. turneri intergenic regions detected five genes (cyoD, infA, rpmJ, ycaR and yidD), one frameshifted gene (ybeY) and one pseudogene of a hypothetical protein (BTURN675_514).

Phylogenetic analysis

We classified genes into six MultiFun categories (cell processes, cell structure, information transfer, metabolism, regulation, and transport) by searching for the gene name in the EcoCyc database. Some genes are assigned to more than one MultiFun category, and we included all categories listed for each gene. If the gene had no associated MultiFun terms in EcoCyc or if it had no gene name, such as BPEN_040, we considered the gene unclassified.

To construct a phylogeny, we chose Baumannia (NC_007984), Hamiltonella (NC_012751) and Sodalis (NC_007712) as outgroups. We identified orthologs of Blochmannia genes in these genomes using the Reciprocal Smallest Distance (RSD) algorithm (Wall, Fraser & Hirsh, 2003) with default values for divergence (0.8) and e-value (10⁻⁵). For each of the six MultiFun categories, we randomly selected five genes present in all taxa. There were no duplicates in the resulting set of 30 genes. We excluded B. chromaiodes from the phylogeny because its genome sequence is 98.0% identical to that of B. pennsylvanicus (Williams & Wernegreen, 2013). We used TranslatorX (Abascal, Zardoya & Telford, 2010) and MAFFT (Katoh et al., 2005) to construct a multiple sequence alignment for each gene, which we then trimmed with ZORRO (Wu, Chatterji & Eisen, 2012). We concatenated the trimmed amino acid alignments and used MrBayes v3.2.1 (Ronquist & Huelsenbeck, 2003) to construct a majority rule consensus tree. To test whether maximum likelihood methods produced the same topology, we analyzed the same dataset using MEGA v6.0 (Tamura et al., 2013) with a cpREV + G + I + F amino acid substitution model, which had the lowest Bayesian Information Criterion (BIC) score, and assessed branch support using 500 bootstrap replicates.

GC skew and DnaA box motif search

For B. pennsylvanicus (NC_007292), B. obliquus, and B. turneri, we used DNAPlotter (Carver et al., 2009) to construct plots of GC skew with 500 bp window size and 50 bp step size. We used Pattern Locator (Mrazek & Xie, 2006) to search these three genome sequences for the consensus DnaA box motif TTWTNCACA.

BLAST analysis

Because we prepared genomic DNA from whole ants or gasters, the Illumina read datasets include coverage of genomes other than Blochmannia, such as the ant host nuclear and mitochondrial genomes. To determine if the de novo assemblies included contigs from potential Blochmannia plasmids or other bacterial symbionts, we used BLASTN to align all ≥500 bp contigs to the GenBank non-redundant database, limited to bacteria (taxid 2). For contigs with at least one hit of >30% coverage and an evalue of <10⁻⁵, we aligned each contig against the full non-redundant database and examined the top hits.

To test whether the Illumina read datasets contain evidence of riboflavin biosynthesis genes, we constructed BLAST databases of the C. obliquus and P. turneri read datasets. For C. obliquus, we built a BLAST database using the unaligned reads file generated by Mosaik during the initial finishing step (see above), thereby excluding most reads originating from the rib pseudogenes within the B. obliquus genome. We compiled a set of query genes representing orthologs of conspicuously absent genes, including ribABDEF from B. pennsylvanicus (NC_007292), Baumannia (NC_007984) and E. coli MG1655 (NC_000913), ribH from B. pennsylvanicus and Baumannia, ribC from E. coli and the Camponotus floridanus gene for EF-1alpha-F2 (EFN72500). For analysis of C. obliquus reads, we also included ribF and rib pseudogenes from B. obliquus in the query set. We used BLASTN to align the query genes against both BLAST databases and examined the output for alignments with e-values <10⁻⁵.

Certain gene distribution patterns may be explained by either parallel gene loss or acquisition via horizontal gene transfer. To test the hypothesis of horizontal gene transfer, we aligned each of the 55 genes with such patterns against the GenBank non-redundant database using BLAST. We first used BLASTN, and if this search did not return significant hits, we used BLASTX. We examined the taxonomic assignments of the top hits to identify genes with high scoring alignments to bacteria outside of the Enterobacteriaceae, which may support the hypothesis of horizontal gene transfer.

B. obliquus has the fewest protein-coding genes and most pseudogenes

The size and GC content of the B. obliquus and B. turneri genomes are within the ranges observed for Blochmannia of Camponotus sensu stricto (Table 1). Although B. obliquus is on the upper end of the size range, it has the fewest protein-coding genes of the six sequenced Blochmannia. We also detected 15 pseudogenes in B. obliquus, which is an unusually high number for Blochmannia (Table S1). In both B. obliquus and B. turneri, we identified genes that have frameshifts in homopolymer tracts but otherwise are expected to encode a full-length protein. We consider these genes to be functional, because they may be corrected by polymerase slippage during transcription and expressed as full-length proteins (Tamas et al., 2008; Wernegreen, Kauppinen & Degnan, 2010).

B. obliquus is the earliest diverging Blochmannia lineage sequenced

To determine the evolutionary relationships of the Blochmannia lineages, we constructed a phylogeny of 30 randomly selected protein-coding genes shared among Blochmannia genomes and the outgroups Baumannia, Hamiltonella and Sodalis. Using Bayesian methods, the resulting phylogeny has 100% posterior probability at all nodes (Fig. 1). Maximum likelihood methods produced the same topology with ≥99% bootstrap values. In this phylogeny, B. turneri and Blochmannia of Camponotus sensu stricto are more closely related to each other than to B. obliquus. This topology is congruent with recent phylogenies of the ant host taxa that show Colobopsis as a separate lineage rather than a subgenus of Camponotus (Brady et al., 2006). Hereafter, we use the phrase “Blochmannia of Camponotus” to mean Blochmannia of Camponotus sensu stricto, which does not include Colobopsis.

Assemblies show evidence of Wolbachia but no Blochmannia plasmids or other bacterial symbionts

Because we prepared genomic DNA from whole ants or gasters, the Illumina read datasets include coverage of genomes other than Blochmannia. To determine whether the assemblies generated by Velvet included contigs from other symbionts or possible plasmids, we aligned all contigs ≥500 bp against the GenBank non-redundant database using BLASTN. We limited our first search to bacteria (taxid 2). For contigs with hits of evalue <10⁻⁵ and >30% query coverage, we aligned them against the full database. The top hits from this search are shown in Table S2.

Analysis of the C. obliquus assembly identified 15 contigs with BLAST hits that met the criteria outlined above. Three of these are from the ant nuclear genome and one is from the ant mitochondrial genome (Table S2). The best-scoring hits for the remaining 11 contigs are all matches to Wolbachia, an endosymbiont found in many insect species (Werren, Baldo & Clark, 2008), including ants (Russell et al., 2012). We did not detect evidence of any other symbiont or any Blochmannia plasmids.

Analysis of the P. turneri assembly identified eight contigs with BLAST hits that met the criteria outlined above. Four of these are from the ant nuclear genome and three are from the ant mitochondrial genome (Table S2). The best-scoring hit for the remaining contig is to B. pennsylvanicus. To examine this further, we aligned the contig to the B. turneri genome sequence, which produced a much better alignment with 100% coverage and 85% identity. The contig spans the 3′ end of prfB, a 19 bp intergenic region and the 5′ end of lysS. Coverage of this contig averages 54x, whereas coverage of the same region in the B. turneri genome sequence averages 1,380x. We aligned the full read dataset against both the B. turneri genome sequence and the contig and then did a de novo assembly using only the mapped reads. The resulting assembly did not reconstruct the contig, suggesting that it was an artifact of assembly rather than evidence of a plasmid or other symbiont.

Gene content of reconstructed Last Common Ancestor is highly conserved in divergent Blochmannia lineages

Using the six sequenced Blochmannia genomes, we reconstructed the likely gene content of the Last Common Ancestor (LCA) of these lineages. Horizontal gene transfer (HGT) is considered unlikely in Blochmannia due to the isolated intracellular niche of these bacteria. Conservation of gene order in Blochmannia genomes supports the hypothesis that HGT is very rare. Additionally, we used BLAST to identify the best hits for genes with distribution patterns consistent with either independent loss in multiple lineages or acquisition via HGT, and we did not find evidence of HGT (see below). For these reasons, we expect very little, if any, gain of genes in Blochmannia lineages, and our reconstruction of the LCA includes any intact gene found in at least one of the sequenced genomes.

By this definition, the genome of the LCA consists of 690 genes (Fig. 2). Most of these genes are retained in the Blochmannia genomes, with 568 genes, or 82% of the gene content of the LCA, found in all six species. The complete ortholog table can be found as Table S3. Our reconstruction may underestimate LCA gene content, because it is possible that genes were lost independently from all lineages at some point after divergence from the LCA. For example, B. chromaiodes, B. pennsylvanicus and B. obliquus encode a pseudogene of uvrD, but no sequenced Blochmannia has an intact uvrD. This is the only Blochmannia gene found solely as a pseudogene. This gene may have been intact and functional in the LCA, but it is not represented in our reconstruction here. Sequencing of additional taxa, including Blochmannia from other deeply divergent ant host genera such as Opisthopsis, will further refine reconstruction of the ancestral lineage.

Gene content differences among divergent Blochmannia lineages include important cellular functions

DNA replication and repair

The dnaA chromosomal replication initiation protein is intact in B. turneri, detectable as a pseudogene in B. obliquus and missing in all four sequenced Blochmannia of Camponotus (Fig. 2). In other gamma-proteobacteria, DnaA initiates replication by binding to 9-bp sequences called DnaA boxes within the origin of replication (Zakrzewska-Czerwinska et al., 2007). GC skew analysis predicts that the origin of replication is adjacent to mnmG in most Blochmannia species (Fig. 3), so we searched the intergenic regions flanking mnmG in B. obliquus, B. pennsylvanicus and B. turneri for the consensus DnaA box motif TTWTNCACA (Schaper & Messer, 1995).

We found two matches in B. turneri, none in B. obliquus and one match in B. pennsylvanicus. By comparison, Buchnera aphidicola, which retained dnaA, has two DnaA boxes in oriC (Mackiewicz et al., 2004); therefore, it is possible that DnaA can initiate replication in B. turneri with only two consensus DnaA boxes. However, the presence of a DnaA box in B. pennsylvanicus, which lacks dnaA, suggests that these motifs are not necessarily associated with DnaA function in Blochmannia. In fact, when we searched the entire length of each genome, we found multiple matches to the consensus DnaA box motif, which are unlikely to be involved in DnaA function as described for gamma-proteobacteria such as E. coli (Hansen et al., 2006). More work is needed to understand how the differential distribution of DnaA in Blochmannia species affects control of DNA replication and whether DnaA boxes play a role in replication initiation in B. turneri.

Loss of DNA repair mechanisms is a common characteristic of obligate intracellular symbionts of insects. In the six sequenced Blochmannia, two genes involved in base excision repair, mutM and mutY, are differentially distributed. B. obliquus encodes mutM but is missing mutY, whereas B. turneri and the four Blochmannia of Camponotus are missing mutM and encode mutY (Fig. 2). Both protein products act on 8-oxoG, which can mispair with adenine (Michaels et al., 1992). MutM excises 8-oxoG when it is paired with cytosine, thereby initiating base excision repair. If 8-oxoG is not removed prior to replication, MutY removes the mispaired adenine, enabling repair. Inactivation of either gene leads to an increase in GC-to-TA transversions, which may contribute to the reduced genomic GC content observed in Blochmannia and other insect endosymbionts (Lind & Andersson, 2008).

Vitamin biosynthesis

Pyridoxal 5′-phosphate (PLP), the catalytically active form of vitamin B6, is an important enzyme cofactor. With very few exceptions, insects and other animals cannot synthesize vitamin B6 (Tanaka, Tateno & Gojobori, 2005). Biosynthesis of PLP in E. coli and other gamma-proteobacteria occurs via a pathway encoded by seven genes (Mukherjee et al., 2011). This pathway is conserved in B. turneri and Blochmannia of Camponotus, with the exception of epd/gapB, which is missing from all six sequenced Blochmannia (Fig. 4). The function of epd may be fulfilled by gapA, which is present in all six species and can compensate epd mutants in E. coli (Yang et al., 1998). In B. obliquus, only dxs and serC are conserved, whereas pdxABHJ are missing. We detected pseudogenes of three of these genes. Based on these gene losses, B. obliquus appears unable to provision its ant host with vitamin B6.

Riboflavin (vitamin B2) biosynthesis is another vitamin synthesis pathway with differential gene distribution in these Blochmannia lineages. Riboflavin is essential for synthesis of the cofactors flavin mononucleotide (FMN) and flavin adenine dinucleotide (FAD) (Abbas & Sibirny, 2011). Blochmannia of Camponotus encode five riboflavin biosynthesis genes (ribABCDH), which comprise the complete pathway in E. coli. By contrast, these genes are missing in both B. obliquus and B. turneri. We detected pseudogenes of ribABCD, but not ribH, in B. obliquus. The ribF gene, which encodes an enzyme that synthesizes FMN and FAD from riboflavin, is conserved in all six sequenced Blochmannia species.

It is possible that riboflavin biosynthesis genes are encoded on a plasmid, by a secondary symbiont (Lamelas et al., 2011), or even within the ant host nuclear genome (Husnik et al., 2013). As discussed above, we detected Wolbachia but no plasmids or other symbionts in the assemblies. To investigate possible alternative sources of riboflavin, we analyzed the Illumina read datasets. Using BLASTN, we aligned a query set of rib genes from B. pennsylvanicus, Baumannia and E. coli against the reads. We also included the elongation factor alpha F2 (EF-1alpha-F2) gene from Camponotus floridanus in the query set to assess representation of the ant host nuclear genome in the read datasets. Because B. obliquus has pseudogenes of four riboflavin biosynthesis genes, we used the subset of reads that did not align to the B. obliquus genome to build the C. obliquus BLAST database.

We found no evidence of intact riboflavin biosynthesis genes in the read datasets. For reads from P. turneri gasters, we detected no BLASTN hits with e-values <10⁻⁵ to the query riboflavin biosynthesis genes. For reads from C. obliquus ants, some reads aligned to the B. obliquus pseudogenes. Based on examination of the BLAST alignments, these reads likely originated from the pseudogenes, but they have more mismatches than permitted by Mosaik and were included in the unaligned reads file. When we considered only reads that did not align by BLASTN to the rib pseudogenes, we detected only three reads aligning to the query riboflavin biosynthesis genes with evalues <10⁻⁵. A single read, but not its mate, aligned to ribE from B. pennsylvanicus with only 58% coverage and an evalue of ∼10⁻⁶, and a single read pair aligned to ribA from E. coli with 100% identity and 100% coverage, which may be due to contamination or presence of gut-associated bacteria. These single reads contrast with the 558x coverage of the B. obliquus genome. For both read datasets, multiple reads aligned to the C. floridanus EF-1alpha gene with evalues <10⁻²⁰, confirming that the reads include coverage of the host nuclear genome. Our analysis of the Illumina read datasets did not reveal how these symbiotic systems compensate for the loss of riboflavin biosynthesis in the Blochmannia partners.

Phylogenetic framework reveals parallel gene losses

We placed gene losses in a phylogenetic context to identify parallel losses (Fig. 5), which we define as independent loss of the same gene in multiple lineages separated by a lineage that retained the gene. By this definition, we identified 55 genes that were subject to parallel loss in these Blochmannia lineages. An alternative explanation for the distribution patterns of these 55 genes is acquisition by horizontal gene transfer. To test this explanation, we aligned each of the 55 genes to the GenBank non-redundant database using BLAST and identified only two genes (icd and yqiC) for which the top hits were to bacteria outside of the Enterobacteriaceae. The top BLASTN hit for icd was to Candidatus Profftella, a beta-proteobacterial endosymbiont of the Asian citrus psyllid. BLASTN did not return any significant hits for yqiC, but the top BLASTX hit was to Vibrio litoralis. Although it is possible that these genes were acquired by Blochmannia via horizontal gene transfer, for the purposes of our analysis we consider parallel gene loss a more likely explanation for their distribution patterns.

Because of the branching order in the Blochmannia phylogeny (Fig. 1), the 15 genes unique to B. turneri must have been independently lost in both B. obliquus and Blochmannia of Camponotus (Fig. 5), assuming that gene content differences reflect gene loss rather than HGT (see above). Four of these 15 genes are involved in DNA replication and repair, including the chromosomal replication initiation protein dnaA, the DNA mismatch repair protein mutS and DNA topoisomerases topA and topB. Three of the 15 genes (acnB, gltA and icd) encode enzymes of the tricarboxylic acid (TCA) cycle and complete this pathway in B. turneri. By contrast, B. obliquus and Blochmannia of Camponotus retained only TCA genes involved in energy generation. In general, parallel loss of these 15 unique B. turneri genes in B. obliquus and Blochmannia of Camponotus may indicate differing selective pressures on B. turneri and its ant host.

Multiple genes involved in metabolic pathways were subject to parallel loss in the six sequenced Blochmannia. Glutamine synthetase was independently lost in both B. vafer and B. obliquus. In other Blochmannia, this enzyme may play an important role in nitrogen recycling for the ant host (Feldhaar et al., 2007). Parallel loss of glnA suggests that this enzyme is not essential to the symbiosis and its function can be fulfilled by an alternative pathway, such as arginine biosynthesis via carbamoyl phosphate synthase (Williams & Wernegreen, 2010). A few parallel loss events involved genes from the same metabolic pathway, such as coaADE and dfp (coenzyme A biosynthesis), ispADFGH (terpenoid backbone biosynthesis) and ribABCDH (riboflavin biosynthesis). Loss of these genes may reflect relaxed selection on components of a pathway after inactivation of a gene within that pathway.

In addition to genes involved in metabolism, genes encoding outer membrane and secretory proteins, including ompA, oprC, secBDF and tonB, were lost independently in multiple lineages, possibly affecting communication between endosymbiont and host cells. Finally, parallel loss of five hypothetical proteins with unknown function emphasizes the need to better characterize these proteins and understand their contributions to the symbiosis. Maintenance of these hypothetical proteins in particular lineages suggests that they are under selective pressure specific to host lineages.

Blochmannia lineages experienced multiple inversions

Although gene order is conserved among the four sequenced Blochmannia of Camponotus, comparisons including B. obliquus and B. turneri revealed inversions in eight genomic regions involving between two and 34 genes (Figs. 3 and 6). Two of these regions (Figs. 3D and 6D, Figs. 3G and 6G) show evidence of multiple separate inversion events (Table S3), suggesting that they may be inversion “hotspots.” To determine which lineage likely experienced inversion events, we compared gene order in the eight regions to that of close relatives strain HS (CP006569), which is a member of a Sodalis-allied clade, and Baumannia cicadellinicola str. Hc (NC_007984). For some regions, these close relatives provide information on the probable gene order of the ancestral Blochmannia lineage. Based on these comparisons, we hypothesize that two inversions occurred in the lineage leading to B. obliquus (Figs. 3A/6A and Figs. 3B/6B), one occurred in the lineage leading to B. turneri (Figs. 3F/6F), one occurred in the ancestral lineage of B. turneri and Blochmannia of Camponotus (Figs. 3C and 6C), and one occurred in the ancestral lineage of the four sequenced Blochmannia of Camponotus (Figs. 3E and 6E). For the remaining three regions, our comparisons were inconclusive.