Genomic data of Bartonella and other bacterial organisms related to this study were downloaded from the National Center for Biotechnology Information (NCBI) GenBank (http://www.ncbi.nlm.nih.gov/, last accessed August 15, 2014) or from the website of the Bartonella Group Sequencing Project, Broad Institute of Harvard and the Massachusetts Institute of Technology (http://www.broadinstitute.org/, last accessed August 15, 2014). Bartonella species were grouped into four lineages (L1–L4) plus B. tamiae and B. australis, following the current taxonomy (Engel et al. 2011; Pulliainen and Dehio 2012; Guy et al. 2013). For the purpose of this study, we will refer to all bartonellae except B. tamiae as eubartonellae. This is based on the recognition that B. tamiae has been described as clearly distinct from all other currently known bartonellae lineages (Kosoy et al. 2008; Guy et al. 2013). A total of 28 Bartonella species were examined in this study (table 1). Table 1

Basic Information of Bartonella Genomes and Corresponding Samples Assessed in This Study

Initial discovery analysis of putative HGT events in metabolic pathways was assisted by an automated pipeline (available in the Dittmar Lab: https://github.com/DittmarLab/HGTector, last accessed August 15, 2014). This pipeline is based on a computational method of rapid, exhaustive, and genome-wide detection of HGT, featuring the systematic analysis of BLAST hit distribution patterns combined with a priori defined hierarchical evolutionary categories (Zhu et al. 2014). Batch BLASTP of Bartonella protein-coding genes was performed against the entire NCBI nr database (E value cutoff = 1 × 10⁻⁵, other parameters remain default). Genes that have less than a statistically relevant threshold of the expected number of hits based on known close relatives of Bartonella, but meanwhile show multiple top hits from taxonomically distant organisms (non-Rhizobiales groups), were considered to be candidates of HGT-derived genes and were subject to further phylogenetic analyses (see below) (see Zhu et al. 2014 for details on pipeline). Particular attention was paid to genes involved in the core central intermediate metabolism and cell wall formation, which have been identified in previous studies on bacterial metabolism (Zientz et al. 2004).

Phylogenetic analyses were employed to validate the putative horizontal and vertical histories of the genes identified in the initial discovery screen. Phylogenetic patterns nesting a Bartonella gene within a homologous gene clade of a candidate donor group, or as strongly supported sister group of a candidate donor group, were considered significant evidence supporting the horizontal transfer from this particular donor to Bartonella (Koonin et al. 2001; Nelson-Sathi et al. 2012; Husnik et al. 2013; Schonknecht et al. 2013). Nucleotide sequences of metabolic genes of interest (i.e., phospholipid pathway) were extracted from Bartonella genomes as well as genomes of selected organisms that represent the putative donor group and its sister groups. Sequences were aligned in MAFFT version 7 (Katoh and Standley 2013), using the L-INS-i algorithm. The MAFFT program was called from the “Translational Align” panel of Geneious 6.1 (Biomatters 2013). Alignment edges were trimmed manually, if needed. The phylogenies of single-gene families were reconstructed based on nucleotide and amino acid sequence alignments (to check for congruence) in a Bayesian Markov chain Monte Carlo (MCMC) statistical framework using MrBayes 3.2 (Ronquist et al. 2012), as well as a maximum likelihood (ML) method implemented in RAxML 7.7 (Stamatakis 2006). The Bayesian MCMC runs had a chain length of 20 million generations, with the sample frequency set as 1,000. The optimal nucleotide substitution models for all three codon positions were computed in PartitionFinder 1.1 (Lanfear et al. 2012). Three independent runs were performed for each data set to ensure consistency among runs. Trace files were analyzed in Tracer 1.5 (Drummond and Rambaut 2007) to check for convergence in order to determine a proper burn-in value for each analysis. A consensus tree was built from the retained tree-space, and posterior probabilities are reported per clade. The ML was run implementing the GTR + G model (for all codon positions) and a bootstrap analysis was performed to gauge clade support.

In order to determine the frequency, components, and boundaries of the putatively horizontally transferred genetic material, genomic environments were manually examined in Geneious 6.1 (Biomatters 2013). Our assumptions are that multiple independent transfers of a gene would likely result in different gene environments being affected. Likewise, if different Bartonella species share the same gene environment adjacent to horizontally transferred genetic material, and the transferred genes follow the previously detected vertical evolutionary pattern for bartonellae, presumably a single ancestral HGT event can be inferred for all species in that lineage. Putatively HGT-derived genes and their adjacent genomic elements were identified in recipient and donor genomes and compared across species and within lineages. Results from this analysis were mapped onto the Bartonella species tree (see below).

Selection analyses were carried out to gage selective pressures operating on all genes in the phospholipid pathway. Selection was assessed using the ML method in the Codeml program of the Phylogenetic Analysis by Maximum Likelihood (PAML) 4.7 package (Yang 2007). As the first step, an analysis under the one-ratio model (M0) was performed to estimate a global ω value (dN/dS ratio) across the phylogenetic tree. Global selective pressures were assessed using the site models (M1a, M2a, M8, and M8a). Evolutionary rates of particular branches of interest (ω₁) versus the background ratio (ω₀) were computed using the branch model (model = 2). Selective pressures operating on subsets of sites of these branches were calculated using the branch-site models (model A and A1). The significance of change of ω value and evidence of positive selection was assessed using the likelihood ratio test. Positive sites were identified using the Bayes Empirical Bayes (BEB) analysis (Yang et al. 2005).

The tertiary structure of the GpsA (NAD(P)H-dependent glycerol-3-phosphate dehydrogenase) protein in Coxiella burnetii (Gammaproteobacteria: Legionellales) was used to model the position of the identified sites with positive selection in horizontally transferred gpsA genes (Seshadri et al. 2003; Minasov et al. 2009).

The possibility that the horizontally acquired gpsA genes underwent convergent evolution in Bartonella, relative to their ancestors was explored. Potential ancestral states of the gpsA genes before HGT were reconstructed under model A in Codeml of PAML (see above). The sequence was then compared with the consensus sequence of extant Bartonella species. A statistical approach recently introduced by Parker et al. (2013) was applied to identify the signatures of convergent evolution of gpsA versions after horizontal acquisition. In brief, this method is based on the significance of differences in sitewise log-likelihood supports among a commonly accepted species tree and given alternative convergent topologies under the same substitution model.

In order to parsimoniously map HGT events to the evolutionary history of bartonellae, phylogenetic relationships among Bartonella species were inferred using standard phylogenetic and phylogenomic approaches as follows:

Phylogenomic analysis: The proteomes of 23 annotated Bartonella genomes were downloaded, from which orthologous groups (OGs) were identified using OrthoMCL 2.0 (Li et al. 2003) with default parameters (BLASTP E value cutoff = 1 × 10⁻⁵, percent match cutoff = 80%, MCL inflation parameter = 1.5). OGs that have exactly one member in each and every genome were isolated, resulting in 516 OGs. Members of each of these OGs were aligned in MAFFT (Katoh and Standley 2013) and refined in Gblocks 0.91b (Castresana 2000) to remove problematic regions. An optimal amino acid substitution model for each OG was computed in ProtTest 3.3 (Darriba et al. 2011) using the Bayesian information criterion. The 516 alignments were concatenated into one data set, based on which a phylogenetic tree was reconstructed using the ML method as implemented in RAxML (Stamatakis 2006) with 100 fast bootstrap replicates. Bartonella tamiae was used as an outgroup, with all other bartonellae treated as ingroup.

In order to further explore the distribution of Bartonella clades with HGT-derived metabolic genes in blood-feeding insects, we screened a global sampling of 21 species of Siphonaptera and Hippoboscoidea for gpsA sequences (table 2). All of these samples had been positive for Bartonella gltA gene and 16S rRNA detection by polymerase chain reaction (PCR) in previous analyses (Morse et al. 2012, 2013). Genomic DNA was extracted from each individual specimen, using the DNeasy Blood & Tissue Kit (Qiagen Sciences Inc., Germantown, MD), following the animal tissue protocol. The quality and concentrations of DNA were assessed with a NanoDrop spectrometer (Thermo Fisher Scientific, Wilmington, DE). Bacterial gpsA diversity was assessed by amplification of gpsA genes from each sample using specific primers and reaction conditions: Helicobacter-derived gpsA (He) (see Results) forward: 5′-ATG AAA ATA ACA RTT TTT GGW GGY GG-3′, reverse: 5′-TTA ATA CCT TCW GCY ACT TCG CC-3′; Enterobacteriales-derived gpsA (Ar/Se) forward: 5′-GGT TCT TAT GGY ACY GCW TTA GC-3′, reverse: 5′-TAR ATT TGY TCG GYA ATT GGC ATT TC-3′. Subsequent TA cloning (if applicable) was performed to isolate amplicons. Based on previous studies of the microbial diversity of bat flies, we expect a subset of species to harbor Arsenophonus and like organisms (ALOs) as endosymbionts (Morse et al. 2013; Duron et al. 2014). In these species, we specifically targeted Arsenophonus-type gpsA for comparative purposes. Sequence analysis and phylogenetic analysis followed the standard protocols described above. Table 2

PCR-Verified gpsA Types in Bartonella-Positive Insect Samples

Given the conserved nature of the bacterial phospholipid pathway, the ancestral loss of glpK and the Glp system in stem bartonellae followed by a gpsA loss at the base of eubartonellae likely created an ancestral population of bartonellae unable to use either glycerol or glucose metabolites (fig. 1). Our analyses therefore suggest that the ancestors of extant eubartonellae relied directly on G3P, which can be imported into the cytoplasm by the Ugp system that remains intact in all bartonellae. G3P is known to be an intermediate metabolite of a strictly intracellular biochemical pathway in prokaryotic and eukaryotic cells, and does not occur stably in the extracellular environment (e.g., blood) (Cronan 2003; Spoering et al. 2006). Therefore, G3P capture and utilization by ancestral bartonellae were likely accomplished by cytoplasmic associations to a living cell. Previous research has shown that in the absence of readily available G3P, the ability of gpsA mutant bacteria to form cell membranes is severely compromised, resulting in the cessation of cell growth (Bell 1974; Cronan and Bell 1974). This functional peculiarity may explain the slow extracellular growth of the gpsA-less lineage 3 bartonellae on blood agar and their more successful isolation in living cell lines (Heller et al. 1997; Podsiadly et al. 2007). Specifically, the four extant representatives of lineage 3 bartonellae (B. rochalimae, B. clarridgeiae, B. sp. 1-1C, and B. sp. AR 15-3) are possibly the surviving representatives of an ancient lineage, as they are still lacking glpK, the Glp system, as well as gpsA. The ubiquitous loss of these important genes in the phospholipid pathway prior to the evolution of extant eubartonellae certainly suggests that bartonellae had an early intracellular beginning.

Our results provide strong evidence that bartonellae gpsA was acquired from three independent prokaryotic sources outside of alpha-proteobacteria after a single initial loss at the base of the eubartonellae lineage. The repeated retentions of HGT-derived gpsA in the Bartonella genomes confirm the functional importance of gpsA, in the context of the loss-and-regain hypothesis of Doolittle et al. (2003). Based on an array of studies related to HGT, it has been hypothesized that genes that are selectively advantageous in the new organisms have a higher probability of being retained (Kuo and Ochman 2009). Results suggest vertical inheritance and global stabilizing selection after gpsA transfer in Bartonella lineages, as well as the maintenance of ORFs in all transferred genes. Taken together with the previously confirmed expression of gpsA in bartonellae (Saenz et al. 2007; Omasits et al. 2013), the above evidence supports the functionality of the gpsA genes after transfer. Furthermore, the protein sequence alignment shows that all functionally important sites are conserved among the HGT-derived gpsA versions (fig. 3). This, combined with the overall stabilizing selection on horizontally transferred gpsA genes (see Results), implies that the acquired genes are likely to have inherited their original functional role in the biosynthesis of bacterial membrane lipids. However, the significant elevation of evolutionary rates in all three acquired genes and the detection of specific sites under positive selection suggest that amid the overall stabilizing selection, the genes may still undergo functional evolution to adapt to bartonellae-specific metabolic pathways. Moreover, although different gpsA genes were inserted into distinct genomic loci, it is notable that their genomic contexts typically include clusters of genes that are involved in bacterial lipid metabolism (see Results; fig. 4). This suggests that they have been integrated into the existing transcriptional regulation machinery of lipid metabolic genes, which may have facilitated their retention (Lercher and Pal 2008). Bartonellae gpsA mutants are known to result in an abacteremic phenotype (Saenz et al. 2007), pointing to the importance of functional gpsA in pathogenicity and hematogenous spread. Therefore, it is possible that the horizontal reacquisitions of functional gpsA facilitated the hematogenous spread of bartonellae to diverse hosts through blood-feeding vectors, as confirmed for the majority of extant bartonellae.

For two prokaryotic organisms to exchange genetic material likely requires interactions in micro- or nanospace, either directly between bacteria (e.g., conjugation) or with transfer agents (e.g., bacteriophages) in an appropriate environment (Frost et al. 2005; Polz et al. 2013). Based on the high likelihood of their intracellular beginnings (as inferred from the ancestral losses in the phospholipid pathway), the identified cases of HGT provide not only clues to infer past shared ecological connections between Bartonella and other bacteria but also point to putative ancestral hosts.

In line with this argumentation, we suggest that the two transfers from Gammaproteobacteria are more likely to have occurred in an arthropod. Specifically, the Arsenophonus-derived transfer likely stems from endosymbionts exclusively associated with arthropods. ALOs, the immediate well-supported sister group of the bartonellae gpsA, is ecologically versatile and widely distributed among arthropods, including blood-feeders, such as ticks (Ixodidae), keds (Hippoboscidae), and bat flies (Nycteribiidae) (Trowbridge et al. 2006; Novakova et al. 2009; Morse et al. 2013; Duron et al. 2014). Among blood-feeding parasites, Arsenophonus and like species are still primary endosymbionts of extant hippoboscoid flies, which are among the confirmed insect vectors of bartonellae of lineages 4 and 2 (Halos et al. 2004; Morse et al. 2013). Therefore, we propose that hippoboscoid flies may have already been among the ancestral blood-feeding arthropod hosts of bartonellae providing a biological reservoir conducive for horizontal transfer of genes. Interestingly, the gammaproteobacterial ALOs have never been detected in extant Siphonaptera, the insect order with the most common and speciose Bartonella vectors (Chomel et al. 2009; Tsai et al. 2011; Pulliainen and Dehio 2012). Yet all lineage 4 bartonellae transmitted by fleas carry the Arsenophonus-derived gpsA (Ar) (fig. 3C). These facts imply that the Arsenophonus-derived transfer of gpsA (Ar) to stem L4 Bartonella likely did not occur in fleas, but that use of Siphonaptera are likely evolutionarily derived vectors of lineage 4 bartonellae.

The gpsA transfer from the enterobacterial Serratia involves B. australis, which at present seems to be the only representative of its lineage, although this may change with better sampling coverage. Serratia species colonize diverse habitats, including plants, insects, and vertebrates (Grimont F and Grimont PAD 2006). Given this wide, and largely underexplored host range, it is difficult to pinpoint a specific ancestral host for the HGT exchange of the Serratia-derived gpsA in Bartonella. However, the Serratia-derived gpsA of B. australis strongly nests within a monophyletic clade containing Serratia symbiotica, a known endosymbiont of aphids (Aphidoidea) (fig. 3C). In insects Serratia may function as pathogen or symbiont, and have been shown to invade the hemocoel and the intestinal tract (Grimont F and Grimont PAD 2006). In mammals, infection often is opportunistic and rarely systemic (unless previous immunosuppression exists) (Grimont F and Grimont PAD 2006; Mahlen 2011). Therefore, it is conceivable that an insect host was involved in this transfer too.

In contrast, we suggest that the horizontal integration of Helicobacter-derived gpsA into a Bartonella genome likely took place in a mammalian host, especially given that bartonellae gpsA (He) is firmly nested within the Helicobacter clade. Helicobacter are predominantly mammalian pathogens (Whary and Fox 2004; Rogers 2012), whose hosts overlap well with the known host range of Bartonella species. In their hosts helicobacter-type bacteria typically colonize the gastrointestinal tract and liver, causing peptic ulcers, chronic gastritis, cancer, and other diseases. Meanwhile, blood is a secondary site for some species, where they adhere to erythrocytes, which is also the dominant infection site of bartonellae (Whary and Fox 2004; Dubois and Boren 2007).

It is important to note that the timing of each inferred horizontal transfer coincides with the subsequent speciation of extant bartonellae lineages, yet it is difficult to assess exactly when these events happened on an evolutionary scale. However, the transferred gpsA genes exclusively affect eubartonellae (=every lineage after B. tamiae), and the transferred genes are still closely related to extant prokaryotic groups, allowing donor identification. This could be in part due to strong functional constraints, but may also point to an evolutionarily recent HGT relative to the total lineage age. This would support a picture of a more recent diversification with invertebrates and mammals, as suggested by the hypothesis of “explosive radiation” of bartonellae by other authors (Engel et al. 2011; Guy et al. 2013).

Experimental Bartonella genotyping and subsequent phylogenetic analyses recover expected lineages given currently known host and vector ranges (Halos et al. 2004; Morse et al. 2012) (table 2), and confirm predictions of gpsA origin based on phylogenomic analyses. Specifically, bartonellae of deer keds (L. cervi, Hippoboscidae), which are known vectors of lineage 2 bartonellae, contain Helicobacter-derived gpsA (He) (fig. 3B). A diverse sampling of flea and bat fly species shows Arsenophonus-derived gpsA (Ar), as is expected for vectors of L4 Bartonella species (fig. 3D). Some samples (e.g., Megabothris walkeri) closely cluster with known L4 Bartonella species (e.g., B. doshiae). Others, such as bartonellae from Trichobius species, appear to be phylogenetically distant from any established Bartonella subclades, suggesting putative novel species in bat flies, which has been proposed previously (Billeter et al. 2012; Morse et al. 2012). These findings call for more in-depth studies to characterize extant Bartonella diversity.

The overall topology of the L4 gpsA (Ar) tree (fig. 3D) does not mirror the phylogeny of either fleas (Whiting et al. 2008) or bat flies (Dittmar et al. 2006; Petersen et al. 2007), suggesting the absence of Bartonella-insect coevolution within these two groups. For bat flies, this is in contrast to their primary endosymbionts, which exhibit notable coevolutionary patterns with their fly hosts (Hosokawa et al. 2012; Morse et al. 2013). Flea and bat fly bartonellae are interspersed among each other in the tree, implying frequent horizontal transmission of L4 bartonellae between insect vectors, and low host specificity.

More detailed coevolutionary analyses of mammal–Bartonella and insect–Bartonella relationships are warranted given these findings.