Evolution of Wolbachia mutualism and reproductive parasitism: insight from two novel strains that co-infect cat fleas

Assembly and annotation of wCfeT and wCfeJ

Our recent report detailed (1) the assembly of Wolbachia reads generated from C. felis genome sequencing, (2) genome circularization and naming of two novel strains (wCfeT and wCfeJ), and (3) complete genome annotation (Driscoll et al., 2020). All relevant information pertaining to wCfeT and wCfeJ is available at NCBI under Bioproject PRJNA622233. Complete genome sequences are available in the NCBI RefSeq database at accession identifiers NZ_CP051156.1 (wCfeT) and NZ_CP051157.1 (wCfeJ).

Wolbachiae phylogenomics

Protein sequences (n = 84,836) for 66 sequenced Wolbachia genomes plus five additional Anaplasmataceae (Neorickettsia helminthoeca str. Oregon, Anaplasma centrale Israel, A. marginale Florida, Ehrlichia chaffeensis Arkansas, and E. ruminantium Gardel) were either downloaded directly from NCBI (n = 48), retrieved as genome sequences from the NCBI Assembly database (n = 13), contributed via personal communication (n = 8; Michael Gerth, Oxford Brookes University), or sequenced here (n = 2) (Table S1). For genomes lacking functional annotations (n = 15), gene models were predicted using the RAST v2.0 server (n = 12) or GeneMarkS-2 v1.10_1.07 (n = 3; (Lomsadze et al., 2018)). Ortholog groups (n = 3,038) were subsequently constructed using FastOrtho, an in-house version of OrthoMCL (Li, Stoeckert & Roos, 2003), using previously established criteria (i.e., an expect threshold of 0.01, percent identity threshold of 30%, and percent match length threshold of 50% for ortholog inclusion (Driscoll et al., 2017, 2020; Smith et al., 2015; Hagen et al., 2018)). A subset of single-copy families (n = 12) conserved across at least 60 of the 66 genomes were independently aligned with MUSCLE v3.8.31 (Edgar, 2004) using default parameters. Regions of poor alignment were masked with trimal v1.4.rev15 (Capella-Gutiérrez, Silla-Martínez & Gabaldón, 2009) using the “automated1” option. All modified alignments were concatenated into a single data set (2,803 positions) for phylogeny estimation using RAxML v8.2.4 (Stamatakis, 2014). The gamma model of rate heterogeneity and estimation of the proportion of invariant sites were used. Branch support was assessed with 1,000 pseudo-replications. Final ML optimization likelihood was −52350.085098.

BOOM characterization

CI gene characterization

Screening colonies and wild populations of cat fleas for wCfeT and wCfeJ

Upon receipt from collaborators, cat fleas were stored at −20 °C in EtOH until processing for DNA extraction and qPCR. Fleas were sexed under a stereomicroscope and surface sterilized for 5 min with 1% bleach solution, followed by 5 min with 70% EtOH, and finally with three washes using molecular grade water. Individual fleas were flash-frozen in liquid N₂ and ground to powder with sterile pestle in a 1.5 ml microcentrifuge tube. DNA was extracted from dissected tissues using the Zymo Quick-DNA Miniprep plus kit following the solid tissue protocol (Zymo Research; Irvine, CA, USA). The presence of wCfeT, wCfeJ, and C. felis DNA was determined via qPCR using the following primer sets: Cfe#18S|179, 5′-TGCTCACCGTTTGACTTGG-3′ and 5′-GTTTCTCAGGCTCCCTCTCC-3′ (Reif et al., 2008); wCfeT#ApaG|75: 5′-GCCGTCACTGGCAGGTAATA-3′ and 5′-GCTGTTCTCCAATAACGCCA-3′, wCfeJ#CinA|76, 5′-AGCAACACCAACATGCGATT-3′ and 5′-GAACCCCAGAGTTGGAAGGG-3′. Each primer set amplicon was sequenced before use to ensure specificity. qPCR was performed using a QuantStudio 3 with the PowerUp SYBR green mastermix in a 20 μl reaction volume containing 2 μl of sample DNA with 400 nM of each primer. Cycling parameters were consistent with the “fast cycling mode” including a 2 min cycle at 50 °C, a 2 min cycle at 95 °C, 40 cycles with 3 s at 95 °C and 30 s at 60 °C followed by a melt curve analysis. Standard curves were generated to quantify the number of DNA copies present in each reaction. Results were expressed as the ratio of Wolbachia amplicons to C. felis 18S rDNA copies to normalize for varying efficiency of DNA extractions. Fleas were considered positive for wolbachiae if the flea 18S rDNA Ct value is 30 or lower and the Wolbachia Ct value is 37 or lower.

Determining wCfeT and wCfeJ localization in C. felis tissues

Elward Laboratories (EL, Soquel, CA, USA) cat fleas stored at −80 °C were used to assay tissue localization of wCfeT and wCfeJ. Ovaries and midguts of individual female fleas were dissected under a stereomicroscope. DNA extraction and qPCR were performed as described above.

Identification of novel Wolbachia parasites

Recently, we assembled and annotated closed genomes from two divergent Wolbachia strains, named wCfeT and wCfeJ, that were captured during genome sequencing of the cat flea (EL Colony) (Driscoll et al., 2020). The 16S rDNA sequences of these two strains are identical to those previously identified in another cat flea colony (Louisiana State University) (Pornwiroon et al., 2007; Sunyakumthorn et al., 2008; Gillespie et al., 2015), which historically has been replenished with EL colony fleas. wCfeT and wCfeJ are broadly similar to other closed Wolbachia genomes in terms of length, GC content, number of protein coding genes, and coding density (Table 1), although it is interesting to note that wCfeT exhibits one of the highest pseudogene counts (second only to wCle) while wCfeJ boasts one of the lowest (second only to wOv). Coupled with other recently sequenced genomes for non-Supergroup A and B strains (e.g., wFol and wVulC), these cat flea-associated wolbachiae blur the lines demarcating genomic traits previously used to distinguish basal and sister Wolbachia lineages from those of Supergroups A and B strains (Lefoulon et al., 2020), particularly regarding genome size and mobile element composition (Lefoulon et al., 2020).

Prior reports have indicated the presence of multiple Wolbachia “types” infecting C. felis, with phylogenies inferred from partial gene sequences placing these unnamed Wolbachia strains in supergroups B or F, or in basal lineages (Gorham, Fang & Durden, 2003; Dittmar & Whiting, 2004; Casiraghi et al., 2005). Our robust genome-based phylogeny estimation reveals that both wCfeT and wCfeJ are divergent from Supergroup A and B wolbachiae (Fig. 1; Table S1). wCfeT is similar to undescribed C. felis-associated strains that branch basally to most other Wolbachia lineages (Gorham, Fang & Durden, 2003; González-Álvarez et al., 2017; Vasconcelos et al., 2018), while wCfeJ is similar to undescribed C. felis-associated strains closely related to Wolbachia supergroups C, D and F (Gorham, Fang & Durden, 2003; Casiraghi et al., 2005). The substantial divergence of wCfeT and wCfeJ from a Wolbachia supergroup B strain associated with C. felis from Germany (wCte) indicates a diversity of wolbachiae that are capable of infecting cat fleas. Importantly, the nature of host-microbe relationships for any of these C. felis-associated wolbachiae is unknown.

wCfeT is equipped with biotin synthesis genes

The enzymatic steps for bacterial biosynthesis of biotin from malonyl-CoA are largely conserved (Lin & Cronan, 2011) and usually involve six bio enzymes and the fatty acid biosynthesis machinery (Lin, Hanson & Cronan, 2010) (Fig. 2A). We previously identified the first plasmid-borne bio gene operon, which occurs on plasmid pREIS2 of Rickettsia buchneri, an endosymbiont of the black-legged tick (Gillespie et al., 2012). We observed a rare gene order for the pREIS2 bio operon shared only with bio operons encoded on the chromosomes of obligate intracellular pathogens Neorickettsia spp. and Lawsonia intracellularis (Deltaproteobacteria). These unique bio operons formed a clade in estimated phylogenies, indicating LGT of the operon across divergent intracellular species. Subsequent studies, using our same dataset for phylogeny estimation, discovered this bio operon in certain Wolbachia (Nikoh et al., 2014; Gerth & Bleidorn, 2017; Balvín et al., 2018), Cardinium (Penz et al., 2012; Zeng et al., 2018) and Legionella (Ríhová et al., 2017) species. We refer hereafter to this intriguing LGT of bio genes as the Biotin synthesis Operon of Obligate intracellular Microbes (BOOM).

wCfeT carries a complete BOOM with greater similarity to Wolbachia BOOM than those of other bacteria (Fig. 2B). Phylogeny estimation of BOOM and other diverse bio gene sets indicates wCfeT BOOM is basal to other Wolbachia BOOM (Fig. 2C; Fig. S1). However, as the divergence pattern for other Wolbachia BOOM is discordant with genome-based phylogeny (Fig. 1), it is clear that multiple independent BOOM acquisitions have occurred in wolbachiae irrespective of lineage divergences. This is supported by a previous study positing recent BOOM acquisition in wNfla and wNleu (Gerth & Bleidorn, 2017) and a current report that detected partial BOOM loci in wApol, a Supergroup S Wolbachia strain infecting pseudoscorpions (Lefoulon et al., 2020). The nature of BOOM flanking regions in the genomes of wCfeT, wVulC, wCle, wStr, and wLug, which collectively lack synteny and are riddled with transposases, recombinases and phage-related elements (Table S3), also attests to LGT as the origin for BOOM in these genomes. Furthermore, analysis of wCfeT BOOM flanking genes reveals hotspots for additional LGT (Fig. 2D), exemplified by an S-adenosylmethionine importer (EamA) gene (Fig. S1B) that is a signature among many diverse obligate intracellular bacteria (Gillespie et al., 2012; Klinges et al., 2019; Haferkamp et al., 2013; Driscoll et al., 2013).

As some BOOM-containing Wolbachia strains have been shown to provide biotin to their insect hosts (Ju et al., 2019; Nikoh et al., 2014), we posit that wCfeT has established an obligate mutualism with C. felis mediated by biotin-provisioning. Genes involved in CI (and putatively in MK) have been predominantly found in the EAM; the absence of such genes in wCfeT’s WO prophage (Fig. S2) provides no evidence that wCfeT can induce these two forms of RP, implying that its maintenance may indeed rely on nutritional supplementation to the cat flea. Wolbachia strains associated with bedbug (wCle), planthopper (wStr, wLug), and termites (wCtub) with evidence for BOOM (Hellemans et al., 2019) reside in bacteriocytes of their insect hosts, structures associated with nutritional symbiosis in many arthropods (Douglas, 2015). Interestingly, C. felis is not known to contain bacteriocytes, and the ability of wCfeT to synthesize and supplement biotin to cat fleas awaits experimentation.

wCfeJ lends insight on the evolution of reproductive parasitism

wCfeJ carries a putative toxin-antidote (TA) operon that is architecturally similar to the wPip_Pel CinA/B TA operon (Fig. 3A), which was previously characterized as a CI inducer in flies (Chen et al., 2019). Containing at least three distinct variants (LePage et al., 2017; Lindsey et al., 2018), all CinB toxins harbor dual nuclease (NUC) domains in place of the deubiquitinase (DUB) domain of CidB, another CI inducer in flies (LePage et al., 2017; Beckmann, Ronau & Hochstrasser, 2017; Beckmann & Fallon, 2013). Chimeric NUC-DUB toxins, termed CndB toxins (Beckmann et al., 2019), also occur as operons (cndAB) in certain Wolbachia and Rickettsia genomes (Beckmann, Ronau & Hochstrasser, 2017; Gillespie et al., 2018). Although their functions remain unknown, these larger CndB toxins belong to an extraordinarily diverse array of toxins with highly modular architectures that occur in divergent intracellular bacteria (Gillespie et al., 2018). Curiously, despite a size similar to the CinB toxins, wCfeJ’s CinB toxin has much higher sequence identity to larger CndB and NUC-containing toxins found mostly in Wolbachia Supergroups A and B (Fig. 3B). A similar pattern was found for wCfeJ’s CinA antidote, indicating the genes were likely acquired as an operon.

In contrast to wCfeT, wCfeJ harbors a highly degraded WO prophage genome with only a partial EAM, the CinA/B operon, and several other genes encoding hypothetical proteins and transposases (Fig. 3C). This is reminiscent of the CI- and male killer-inducing wRec strain, which lacks a functional WO phage yet possesses TA modules within remnant EAMs (Metcalf et al., 2014).

Synteny analysis with other Wolbachia phages revealed wCfeJ cinB and three downstream genes are colinear with a large (3707 aa) wDacB cndB toxin (Fig. 3D). We posit that gene fission of a large modular CndB toxin created the smaller CinB toxin. This interpretation is consistent with our prior hypothesis stating that more streamlined RP-inducing toxins (i.e., cinB and cidB) originate from larger modular toxins (i.e., cndB and others) that are widespread in the intracellular mobilome, particularly Wolbachia phage and Rickettsia plasmids and integrative and conjugative elements (ICEs) (Gillespie et al., 2018, 2015). We find the alternative hypothesis, that multiple gene fusions create larger modular toxins like CndB proteins, less parsimonious provided that RP toxins tend to degrade and pseudogenize rapidly in rickettsial genomes (Martinez et al., 2020; Gillespie et al., 2018).

wCfeJ is unique among described wolbachiae outside of Supergroups A and B by containing a CI-like TA operon. Some of the Wolbachia strains in closely related Supergroups C, F and D are considered mutualists with their nematode or arthropod hosts, supporting the absence of RP-inducing genes in the genomes of these strains. Still, other strains not included in our analysis, namely those belonging to Supergroups S (Lefoulon et al., 2020; Chafee et al., 2010), P (Glowska et al., 2015), and J (Lefoulon et al., 2016) may eventually be shown to carry RP-inducing genes, especially since some of these strains are not known as mutualists. wFol, another basal Wolbachia lineage (Fig. 1), does carry a myriad of large modular candidate RP-inducing toxins but does not have cin, cid or cnd TA operons (Faddeeva-Vakhrusheva et al., 2017; Kampfraath et al., 2019). This indicates LGT of candidate RP-inducing genes occurs across the breadth of Wolbachia diversity. Collectively, we speculate that wCfeJ acquired a WO prophage, possibly from a Supergroup A or B Wolbachia strain that also infects cat fleas (e.g., wCte), and that viral- or gene-level selection streamlined a cndAB locus into a cinAB operon (Gillespie et al., 2018). Given that the genomes of many Wolbachia reproductive parasites harbor diverse arrays of CinA/B-and CidA/B-like operons (Gillespie et al., 2018; Beckmann et al., 2019), we posit that wCfeJ’s CinA/B TA operon might function in CI or some other form of RP; however, further study is required to evaluate wCfeJ as a cat flea reproductive parasite. Furthermore, an intriguing conjecture is that functions of the larger modular toxins, exemplified by wDacB CndB, whatever they be, are likely the primal functions that gave rise to CI.

Lateral transfer of Wolbachia CI-like antidotes to the cat flea genome

Our previous work identified CI-like toxin and antidote genes on scaffolds in several arthropod genome assemblies, although antidote genes were found more frequently possibly due to utility in suppressing microbial invaders implementing RP (Gillespie et al., 2018). Many of these genomes contain multiple divergent antidote genes, though none were found to harbor large-scale Wolbachia gene transfer that is known to occur in certain invertebrates (Leclercq et al., 2016; Hotopp et al., 2007; Aikawa et al., 2009; Kondo et al., 2002; Dhaygude et al., 2019; McNulty et al., 2010; Werren et al., 2010). In light of this observation, we originally speculated that CI-like genes are mobilized from wolbachiae to host through WO phage insertion into host nuclear genomes; however, the absence of other WO phage genes in these genomes made this assertion tenuous. Instead, more recent work describing the first WO prophage carrying two sets of CI loci provides compelling evidence for transposon-mediated LGT of CI loci between divergent wolbachiae (Cooper et al., 2019). This phage-independent mechanism for RP gene jumping in wolbachiae may result in the inadvertent capture of these genes by arthropod host genomes.

When we analyzed the cat flea genome for CI-like genes, we discovered that the largest C. felis scaffold (185.5 Mb) was found to harbor two CI-like antidote genes (Fig. 4A). Remarkably, neither antidote has any significant similarity to CinA of wCfeJ (Fig. 4B), indicating that these CI antidotes originated from other, as yet undescribed wolbachiae. The C. felis CinA-like CDS (XP_026474240.1) has greater similarity to CinA proteins than CidA proteins, while the CidA-like CDS (not annotated in the current C. felis assembly) has much greater similarity to CidA antidotes (Beckmann, Ronau & Hochstrasser, 2017; Beckmann & Fallon, 2013; Bonneau et al., 2018; Madhav et al., 2020; Morrow et al., 2020). The C. felis cinA-like gene is found as an exon within a larger flea gene model (LOC113377978), and is well supported (18× coverage, 7.2 TPM) by transcripts from the 1KITE project (Misof et al., 2014). Since these transcripts were generated in fleas from a different colony (Dr. Michael Dryden, Kansas State University), we conclude that this gene transfer occurred during an historical infection predating the establishment of these two colonies. Finally, the C. felis CinA- and CidA-like antidotes each possess the conserved motifs we previously defined for all CI antidotes (Gillespie et al., 2018), with the exception of a truncation of motifs 5 and 6 in the CinA-like protein (Fig. 4C). As other arthropod-encoded antidotes share this C-terminal truncation, it may signify an expendable domain in the host (i.e., a bacterial secretion signal, a motif involved in cognate toxin activation or delivery, etc.).

The incorporation of a Wolbachia antidote into an exon of a host gene provides one explanation for how such a laterally transferred gene would remain functional in the absence of the native regulatory elements encoded on WO prophage genomes. We previously hypothesized that CI antidote genes would be under strong selection in arthropod genomes if their expression could result in rescuing CI induction by reproductive parasites (Gillespie et al., 2018). In support of this supposition, it has been shown that: host genotypes modulate CI (Metcalf et al., 2014; Hornett et al., 2006, 2014; Bordenstein, Uy & Werren, 2003); some hosts show weak (Dmel and D. suzukii (Cooper et al., 2017; Conner et al., 2017; Hamm et al., 2014)) or strong (Culex pipiens and D. simulans (Poinsot, Charlat & Merçot, 2003; Merçot & Charlat, 2004)) CI; and a Wolbachia strain’s CI can be suppressed when transfected across different hosts (Bordenstein, Uy & Werren, 2003). While none of the arthropod-encoded CI antidotes have been characterized, we propose they impart immunity to toxins secreted by intracellular parasites, curtailing chronic parasite drive into populations. At least four independent evolutionary mechanisms might lead to suppression of CI. Gene family expansion of, and/or overexpression of key factors capable of counteracting CI mechanisms is possible; similarly, it was postulated that evolution in Cid interacting proteins like karyopherin-α or P32 could lead to target site resistance (Beckmann et al., 2019). A third possibility contributing to weakening CI, though not suppression per se, is simply the mutational collapse of the operons themselves due to lack of purifying selection (Martinez et al., 2020; Meany et al., 2019). We supply evidence of a fourth possibility whereby CI might be counteracted by germline expression of acquired CI antidotes. Significantly, our finding raises awareness of a potential barrier to Wolbachia-mediated pathogen biocontrol measures.

Characterization of wCfeT and wCfeJ infection of cat fleas

Coupled with earlier reports indicating multiple early-branching wolbachiae lineages infecting C. felis (Gorham, Fang & Durden, 2003; Casiraghi et al., 2005), our characterization of two divergent Wolbachia strains infecting cat fleas in the EL colony (Driscoll et al., 2020) as well as fleas maintained for a decade at LSU (Pornwiroon et al., 2007; Sunyakumthorn et al., 2008; Gillespie et al., 2015) prompted us to test individual fleas for double infection. In addition to EL fleas, three other colony strains and three geographically diverse wild cat flea populations were also surveyed (Fig. 5A). Aside from the Modesto Wild strain, which is historically replenished with wild-caught cat fleas (Dr. B. Donahue, 2017, personal communication), strong wCfeJ-wCfeT co-infection was found in colony fleas. Conversely, single wCfeT infection dominated wild cat fleas, with only a few fleas from a large sampling throughout Orange Co. CA harboring wCfeJ. While far greater sampling of both colony and wild populations is necessary, our results hint at a strong selection for wCfeT over wCfeJ in nature, which is concordant with our premise that wCfeT’s BOOM might drive mutualism. The ability of wCfeJ to highly infect colony fleas may arise as a consequence of bottlenecks within colonies, which are known to strengthen CI and maternal artificial transmission (Turelli, 1994). No sex ratio distortions were observed in colony or wild fleas, ruling out other forms of RP as driving factors. Nonetheless, the ability of wCfeT and wCfeJ to co-infect individual cat fleas provides a system to study contrasting forces (nutritional symbiosis and RP) that potentially drive multiple Wolbachia infections in a single invertebrate host.

To gain more insight into the nature of the wCfeT-wCfeJ co-infection, we analyzed the distribution of wCfeT and wCfeJ in EL flea tissues (Fig. 5B). Surmising that both Wolbachia strains are maternally transmitted, we dissected out the ovaries of unfed virgin female fleas and compared Wolbachia density in these tissues vs others (primarily midgut). While wCfeT predominantly localizes in flea ovaries, the cif carrier, wCfeJ, mostly localizes to somatic tissues, although minimal co-localization was observed. The significance of this distribution is unclear. wCfeJ may relocate to ovaries upon a mating cue, or perhaps exert CI via some undefined somatic cell mechanism. Transmission of wCfeJ may also be predominantly horizontal via environmental or direct flea contact (per (Pietri, DeBruhl & Sullivan, 2016)). Sampling co-infected fleas at specific timepoints during development, as well as determining wCfeJ tissue localization in the absence of wCfeT infection, will be required to unravel the intertwined transmission dynamics of these Wolbachia strains in cat fleas.

Evolution of Wolbachia mutualism and reproductive parasitism

Order Rickettsiales is a highly diverse group of obligate intracellular bacteria (Gillespie et al., 2012; Driscoll et al., 2013; Darby et al., 2007). A rapid pace of newly identified lineages (Klinges et al., 2019; Castelli et al., 2019; Tashyreva et al., 2018; Kamm et al., 2019; Gruber-Vodicka et al., 2019; Martijn et al., 2015; Schulz et al., 2016; Yurchenko et al., 2018; Schrallhammer et al., 2013; Montagna et al., 2013; Boscaro et al., 2019; Muñoz-Gómez et al., 2019; Mediannikov et al., 2014; Szokoli et al., 2016; Prokopchuk et al., 2019) provides a rich source for evaluating mechanisms of evolution within constraints of the eukaryotic cell. Phage, plasmids, transposons, ICEs, and other insertion sequences are variably present across rickettsial genomes and mobilize numerous genes to offset reductive genome evolution and provide lifestyle altering traits (Gillespie et al., 2012). Our analyses of Wolbachia BOOM and RP genes prompted us to evaluate the recurrent evolution of nutritional-mediated mutualism and reproductive parasitism across diverse rickettsial lineages (Fig. 6).