A set of genes (739 singletons) from our prior report lacked significant similarity to genes in other Rickettsia genomes (Gillespie et al., 2012a); accordingly, these were used as queries in Blastn searches against the NCBI ‘Rickettsia’ databases (taxid 780). Genes with significant similarity to genes from rickettsiae present within the geographic distribution of I. scapularis were removed. Further, we excluded genes encoded on plasmids, those present within chromosomal RAGES, and those encoding transposases or related elements. The remaining genes were then evaluated using bioinformatics analysis to determine the likelihood that they encode functional proteins.

For toxins, proximal (~100 aa) C-terminal sequences of R. buchneri protein REIS_1424 (EER22217) and another rickettsial species “Candidatus Jidaibacter acanthamoeba” str. UWC36 protein NF27_IC00050 (KIE04387) were used as queries in Blastp searches to compile and analyze diverse proteins harboring significant similarity across complementary proximal C-terminal sequences. For antidotes, entire sequences for R. buchneri protein REIS_1423 (EER22217 with an adjusted start site adding 41 aa at the N-terminus) and “Candidatus Jidaibacter acanthamoeba” str. UWC36 protein NF27_IC00040 (KIE04386) were used as queries in Blastp searches to compile and analyze diverse proteins harboring significant similarity across the entire lengths of the queries. Analyses utilized our HaloBlast method, which is a combinatorial Blastp-based approach for interrogating proteins for LGT (Driscoll et al., 2013). Individual Blastp searches were conducted against five distinct taxonomic databases: 1) “Rickettsia” (NCBI taxid 780)”, 2) “Rickettsiales” (taxid: 766) excluding “Rickettsia”, 3) “Alphaproteobacteria” (taxid: 28211) excluding “Rickettsiales”, 4) “Proteobacteria” (taxid: 1224) excluding “Alphaproteobacteria”, 5) “Bacteria” (taxid: 2) excluding Proteobacteria”, and 6) “minus bacteria”). All subjects from each search were ranked by Sm score (= b * I * Q, where b is the bitscore of the match, I is the percent identity, and Q is the percent length of the query that aligned), a comparative sequence similarity score designed to de-emphasize highly significant matches to short stretches of the query in favor of longer stretches of similarity (Driscoll et al., 2013). The “halos” or separate database searches were then compared to one another to determine the taxon with the strongest similarity to the query sequences.

HaloBlast subjects from the searches with REIS_1424 and NF27_IC00050 as queries were analyzed in two ways. First, only sequences matching the proximal (~100 aa) C-terminal sequences of the query were compiled and aligned with MUSCLE using default parameters (Edgar, 2004). The entire alignment was then visualized as sequence logos using WebLogo (Crooks et al., 2004). Second, two representative sequences per halo were selected for domain predictions across the entire protein. EMBL’s Simple Modular Architecture Research Tool (SMART) (Letunic and Bork, 2017) and/or the Protein Homology/analogY Recognition Engine V 2.0 (Phyre2) (Kelley and Sternberg, 2009) were used to predict and evaluate the following domains: UBA (ubiquitin-associated) (Mueller and Feigon, 2002); haemagglutination activity site (Kajava et al., 2001); hemagglutinin repeats (Pfam ID PF13332); Peptidase M43 domain (Rawlings and Barrett, 1995); endonuclease III (Bruner et al., 2000); RHS repeat (Busby et al., 2013); VENN motif (Aoki et al., 2010; Zhang et al., 2011); DUF637: hemagglutinin-/hemolysin-associated domain (PF04830); alanine-rich-conserved phenylalanine (ALF) motif (Yeats et al., 2003); Laminin_G_3 (PF13385); LamG-like jellyroll fold domain (Liu et al., 2007; Weyer et al., 2007); HintN domain (Perler, 1998). Individual protein schemas were generated using Illustrator of Biological Sequences (Liu et al., 2015) with manual adjustment.

HaloBlast subjects from the searches with REIS_1423 and NF27_IC00040 as queries were compiled and aligned with MUSCLE (default parameters), with the entire alignment visualized as sequence logos using WebLogo. Additionally, CDS flanking certain HaloBlast subjects (i.e. those with NCBI reference protein accession numbers) from the searches with REIS_1424 and NF27_IC00050 were evaluated for their size and potential for encoding an N-terminal sequence motif (LS/ADXE/DXQXXXW) determined to be highly conserved in subjects retrieved in Blastp searches with REIS_1423 and NF27_IC00040 as queries. Finally, HMMER (Finn et al., 2011) searches using NF27_IC00040 or NF27_IC00050 were utilized to evaluate our Blastp-based identification and compilation of both toxin and antidote datasets.

For R. buchneri analysis, ticks from New York, New Hampshire, and Pennsylvania (see Figure 6 for locality information) were stored in 100% ethanol at -20°C until isolation. DNA was extracted from I. scapularis adults and nymphs, as well as rickettsiae infecting ISE6 cells (both kindly provided by Drs. Munderloh and Kurtii, University of Minnesota) and using the DNeasy kit (Qiagen) as per manufacturer’s protocols for cell culture and tissue extraction, respectively. Briefly, ticks were surface sterilized with 5 min washes (1% bleach, 70% ethanol, and 1xPBS), cut into quarters with a sterile scalpel blade, incubated with kit-provided digestion buffer with proteinase K at 56°C overnight, and extracted using the tissue protocol with a final elution of 50μl of molecular grade water. Rickettsiae grown in culture were collected in their host cells, and DNA extracted using the cell culture protocol with a 50μl elution with molecular grade water. For analysis of non-target rickettsiae, all bacteria were grown in cell culture prior to DNA extraction and PCR analysis. DNA was qPCR amplified using PowerUp Sybr Mastermix (Thermo) in 20μl reactions containing 400nm of each primer and 50-100ng of DNA. Primers pairs are as follows for R. buchneri-specific targets: Rb-1424-120-F-5’-acaggcgtaaaactagacaatct-3’ with Rb-1424-120-R-5’-aggaaatccaagcttttcaggta-3’ for the amplification of rCRCT and Rb-1423-116-F-5’-gcatagggtttatagcggtgc-3’ with Rb-1423-116-R-5’-ccataagtttcttcctattgtgctt-3’ for the amplification of rCRCA. Rickettsia gltA was amplified for all rickettsiae using the following primers: CSRT-F-5’-tcgcaaatgttcacggtacttt-3’ and CSRT-R-5’-tcgtgcaattctttccattgt-3’ (Stenos et al., 2005). Reactions were amplified under the following conditions: 1 cycle for 2 min at 95°C, 45 cycles at 95°C for 15 sec and 60°C for 30 sec, followed by a melt curve analysis. All primer sets were considered positive if the cycle threshold was 37 cycles or less. All primer sets were validated for range and efficiency of amplification using pCR4-TOPO plasmid standard curves with ligated amplicons. Primer sets described in this manuscript only amplify their intended products as verified by sanger sequencing and melt curve analyses of each reaction. For visualization of qPCR products ( Figure 6 ) gel electrophoresis was performed using a 2% agarose gel with ethidium bromide straining and visualization using a gel imaging station. For transcriptional analysis, R. buchneri growing during log growth in ISE6 cells were collected in 600μl TRIzol (Invitrogen) and RNA extracted using the DirectZol (Zymo) kit using manufacturer’s instructions and on-column DNase treatment. RNA was further DNase treated using the RQ1 DNase (Promega) prior to enzyme removal and concentration with the Zymo Clean and Concentrator-5 kit. The iScript Select cDNA synthesis kit (Bio-Rad) was used for cDNA synthesis using random hexamers in 20μl reactions with 200ng of DNase-free RNA. RNA was determined to be free of DNA with no reverse transcriptase reactions with 200ng of RNA which resulted in no detectable DNA by qPCR. cDNA was analyzed for transcription of rCRCT, rCRCA and I. scapularis β-actin (primers IsActin-95-F-5’-aatcggcaacgagaggttcc-3’ and IsActin-95-R-5’-agttgtacgtggtctcgtgg-3’) using the same qPCR parameters as above.

Our prior analysis of the first sequenced R. buchneri genome (Wikel I. scapularis colony) indicated 32% of CDS were absent from other Rickettsia genomes (Gillespie et al., 2012a). Given that dozens of new Rickettsia genomes have been sequenced and assembled since 2012, we revisited this list and allowed for candidate genes to also be present in the R. buchneri str. ISO7 assembly (Kurtti et al., 2015). Further, to increase the likelihood of a stable gene present in all R. buchneri populations, we excluded genes 1) encoded on plasmids, 2) flanked by transposases, 3) containing annotations reflecting an association with MGEs, and 4) containing Blastp profiles indicating pseudogenization (e.g., gene fragments, split genes, etc.). These collective constraints yielded a small list of candidate genes, one of which (REIS_1424, NCBI accession no. EER22217) was selected for further analysis.

REIS_1424 encodes a hypothetical protein of 263 aa; the homolog in the R. buchneri str. ISO7 assembly (KDO03356) is identical yet has a different predicted start site that adds 35 aa to the N-terminus. The only remaining significant Blastp matches are from R. tamurae str. AT-1, with two CDS spanning the entirety REIS_1424 indicating a split gene (WP_215426163 and WP_032138795). This is consistent with the close phylogenetic position of R. tamurae and R. buchneri (Hagen et al., 2018); however, since R. tamurae has not been reported from the Western hemisphere, these CDS are not a concern for utilization of REIS_1424 as a diagnostic for R. buchneri infection.

REIS_1424-based Blastp searches outside of the Rickettsia taxon database yielded only two significant matches: a 2192 aa protein (hypothetical protein NF27_IC00050, KIE04387) from a rickettsial amoeba-associated endosymbiont, “Candidatus Jidaibacter acanthamoeba” str. UWC36 (Rickettsiales: Midichloriaceae) (Schulz et al., 2016), and a 97 aa protein (hypothetical protein E1266_17330, TDB94289) from the actinobacterium Actinomadura sp. 7K534 (Streptosporangiales; Thermomonosporaceae). Both alignments indicated that 1) REIS_1424, NF27_IC00050, and E1266_17330 share over a dozen conserved residues, 2) the REIS_1424 and NF27_IC00050 match aligns the proximal C-terminal sequences of both proteins, and 3) E1266_17330 is truncated and lacks N-terminal sequence outside of the matches. Unlike REIS_1424, NF27_IC00050- or E1266_17330-based Blastp searches yielded many significant matches to diverse bacteria (discussed further below). However, no functional domains for the region shared across these proteins could be predicted with searches against the NCBI Conserved Domains Database or using SMART.

Given that the REIS_1424-NF27_IC00050 match spanned greater sequence in each protein (~137 aa) and “Cand. J. acanthamoeba” is another rickettsial taxon, NF27_IC00050 (aa residues 2068-2192) was used as a proxy for further Blastp searches and in silico characterization. NF27_IC00050-based HaloBlast analysis revealed strongest similarity to certain non-proteobacterial proteins ( Figure 2A ). All obtained sequences matched NF27_IC00050_2068-2192 at their proximal C-termini ( Figure 2B ). Intriguingly, this cohort of proteins (n = 155) varied greatly in size across regions outside of the conserved C-terminal sequence ( Figure 2C ). A wide assortment of domains was predicted for these proteins, with many having modular architectures and other characteristics of contact dependent growth inhibition (CDI) and/or Recombination hotspot (Rhs) toxins ( Figure 2D ). While no functional domain could be predicted for any of the analogous C-terminal regions, sufficient conservation was found to strongly indicate a unifying functional role ( Figures 2E , S1A ). We hereafter refer to these analogous regions as CDI-like/Rhs-like C-terminal toxin (CRCT) domains, and to larger proteins possessing them as CRCT domain-containing proteins (CRCT-DCP).

The fusion of small, toxin-antidote (TA) pairs to the C-termini of CDI and Rhs toxins has previously been described and is thought to expand the diversity of toxic activities deployed by both CDI and Rhs systems (Aoki et al., 2010; Poole et al., 2011; Zhang et al., 2011; Ruhe et al., 2018). The extreme polymorphic nature of these TA modules indicates bacterial arms races, with selection operating on species- and strain-level recognition that shapes communities (Willett et al., 2015). For instance, many of the CRCT domain-containing proteins we identified in diverse bacteria lack the CRCT domain in closely related strains (data not shown). Furthermore, we found many small proteins containing only the CRCT domain ( Figure 2C and Table S1 ), as well as larger proteins carrying truncated CRCT domains (data not shown), suggesting the mobile nature of these toxins and a rapid “birth and death” process. These observations, combined with small size and limited sequence conservation, collectively challenge computational approaches for identifying these polymorphic toxins. This is evinced by HaloBlast profiles for REIS_1424 that mirror those for NF27_IC00050_2068-2192 once a Rickettsia-specific insertion is removed from the query in Blastp and HMMR searches ( Figures 2E , S1B ).

As CRCTs are often found as TA modules, we interrogated genes up- and downstream of REIS_1424, NF27_IC00050, and the genes encoding the 153 other identified CRCT-DCPs. This revealed probable antidotes, hereafter named CDI-like/Rhs-like C-terminal toxin antidotes (CRCA), adjacent to 37% of the 155 CRCT-DCPs ( Figure 3A ). NF27_IC00040-based HaloBlast analysis mirrored that for NF27_IC00050_2068-2192, revealing strongest similarity to non-proteobacterial proteins ( Figure 3B ). Taxonomic breakdown of these non-proteobacterial proteins for both CRCT-DCPs and CRCAs revealed a majority from Actinomycetia and Cyanobacteria genomes ( Figure 3C ). While all CRCAs are strongly constrained in length (~140 aa), with only a few proteins fused to partial genes with unrelated functions ( Figure 3D ), Blastp profiles for the best scoring matches parsed by taxonomy strongly indicate Rickettsia CRCT/CRCA (rCRCT/CRCA) modules were acquired from distant non-proteobacteria via LGT ( Figure 3E ).

Like analyzed CRCTs, predicted CRCAs possess enough conservation to indicate a common function ( Figures 3F , S2 ). Still, efforts to thread either CRCTs or CRCAs to solved structures of CdiA-CT/CdiI toxin/immunity complexes (Morse et al., 2012; Beck et al., 2014; Johnson et al., 2016; Michalska et al., 2018) were futile. Many of the characterized CdiA-CT toxins are from proteobacterial species and function as Rnases, specifically targeting tRNAs or rRNAs (Willett et al., 2015). Using SMART or searching against the NCBI CDD did not indicate CRCTs harbor nuclease activities, and no similarity of CRCAs to CdiI domains (e.g., cd20694: CdiI_Ct-like) or members of the SUKH superfamily of immunity proteins (Zhang et al., 2011) could be made. CdiI antidotes have been hypothesized to drive CdiA-CT/CdiI module diversification since they evolve faster than CdiA-CT toxins (Nikolakakis et al., 2012). While we observed greater conservation in CRCTs ( Figure 2E ) versus CRCAs ( Figure 3F ), HaloBlast and HMMER (data not shown) searches recovered more CRCA ( Figure 3B ) versus CRCT-DCPs ( Figure 2A ). Still, the presence of CRCT/CRCA modules across diverse bacterial phyla, with some drastic differences in lifestyle (i.e. obligate intracellular versus extracellular, eukaryote-dependent versus environmental, etc.) indicates a common universal cellular target of CRCTs such as membrane, DNA, or RNA previously characterized for all other studied CdiA-CT/CdiI modules (Aoki et al., 2009; Aoki et al., 2010).

We identified a eukaryotic genome harboring a possible LGT of a rickettsial CRCT-DCP/CRCA system. The genome of the smooth cauliflower coral (Stylophora pistillata) contains a gene encoding a large Rhs-like toxin (YbeQ) that was assembled on an unincorporated scaffold with other bacterial-like genes ( Figure 4 ). YebQ has highest similarity to a few smaller rickettsial proteins that are remnants of degraded Rhs-like toxins (see next section), yet consistent similarity to larger toxins from non-proteobacterial genomes ( Table S3 ). While YebQ does not contain a CRCT domain, smaller S. pistillata genes were found encoding the CRCT and CRCA domains with higher similarity to rickettsial equivalents than other bacteria ( Table S3 ). This attests to the mobile nature of CRCT/CRCA modules and their tendency to incorporate into larger bacterial toxins. It also resonates on our prior work showing another aquatic animal, the placozoan Trichoplax adhaerens, contains LGTs from bacteria (particularly rickettsiae) (Driscoll et al., 2013). While the presence of introns in many of these bacterial-like S. pistillata gene models supports integration, mis-assembly of reads from a rickettsial endosymbiont is also possible.

Inspection of the genomic region where rCRCT/CRCA modules have inserted revealed two interesting findings, both of which further attest to the mobile nature of these polymorphic TA modules and their rapid birth and death process. First, the rCRCT/CRCA loci of R. buchneri and R. tamurae occur in a recombination hotspot adjacent to the SecA gene ( Figure 5A ). This region is highly variable across Rickettsia genomes ( Figure S3A ) and contains small CDS with matches to the S. pistillata YebQ Rhs toxin described above ( Figure S4A ). Despite extraordinary variability in the number and size of CDS in this region across Rickettsia genomes, a conserved tRNA-Ala^TGC locus is always present, corroborating prior observations for CidA-CT/CdiAI modules often inserting near tRNA genes in bacterial genomes.

Second, Blastp searches with a concatenated query (REIS_1427-REIS_1423) revealed matches to additional rCRCA proteins, indicating other rCRCT/CRCA modules elsewhere in Rickettsia genomes ( Figure 5A , dashed box). We designated these divergent rCRCA proteins as components of “rCRCT/CRCA-2” modules, with the above-described system named “rCRCT/CRCA-1” modules. Like the rCRCA-1 protein of R. tamurae, most rCRCA-2 genes are truncated or fragmented, yet a complete protein was found for the Rickettsia endosymbiont of Ixodes pacificus (hereafter REIP) ( Figure S4B ). Despite limited similarity (~30%ID), an estimated phylogeny grouped rCRCA-1 and rCRCA-2 proteins together to the exclusion of other CRCA proteins ( Figures 5B ; S4B ). rCRCA-2 loci all mapped to a second recombination hotspot in Rickettsia genomes adjacent to the conserved BamA and tRNA-Thr^CGT genes ( Figures 5C , S3B ). As with rCRCA-1, we could not identify similarity between rCRCA-2 proteins and domains of CdiI proteins or SUKH immunity proteins.

Searching upstream of rCRCA-2 for potential cognate toxins instead yielded additional genes encoding CRCA antidotes that are duplicated and highly divergent from rCRCA-1 and rCRCA-2 proteins ( Figure 5C ). This arrangement of arrayed immunity antidotes is more characteristic of cdiA-CT/cdiI loci in many proteobacterial genomes (Aoki et al., 2010). Indeed, we were able to identify CdiI-like domains in these proteins using the NCBI Conserved Domains Database (cd20694). Accordingly, we named these antidotes rCRCA-3 proteins. Further inspection of rCRCA-3 genes identified eight genomes with additional copies found in a third recombination hotspot between the genes encoding the cytochrome oxidase subunits 2 (CyoA) and 1 (CyoB) ( Figures 5C , S3C ). Seven Rickettsia genomes have rCRCA-3 genes in both the BamA and CyoA/B recombination hotspots, indicating recent recombination between these loci ( Figure S4C ). We designated rCRCA-3a and rCRCA-3b proteins to distinguish between those located at the BamA or CyoA/B recombination hotspots, respectively.

At the BamA recombination hotspot, we identified cognate rCRCT-3a toxins with CdiA-CT-like domains (cd20695) in the R. tamurae and REIP genomes ( Figure 5C ). These two toxins have strongest similarity to counterparts in proteobacterial genomes, particularly Pseudomonas and Moraxella species ( Figures S4D, E ). We modeled the R. tamurae rCRCT-3a toxin to the CdiA-CT structure of Cupriavidus taiwanensis (Kryshtafovych et al., 2018) with high confidence, indicating rCRCT-3 toxins are unrelated to rCRCT-1 toxins ( Figures S4F, G ). The rCRCA-3 antidotes could not be modeled to the CdiI structure of C. taiwanensis or any other CdiI structures, making the association of rCRCT-3 with rCRCA-3 supported by genome proximity alone. Collectively, this analysis of rCRCT/CRCA genes in Rickettsia genomes illuminates recurrent genome integration, possibly in larger Rhs toxins that have degraded over time. The presence of complete, yet divergent rCRCT/CRCA modules in different species of the Tamurae/Ixodes Group (TIG) rickettsiae ( Figures 5D, E ) indicates weaponry for interbacterial antagonism may be functional for these species and implicates a previously unrealized mechanism for rickettsial competition in the same arthropod host

During infection of ISE6 cells, which were originally derived from I. scapularis embryos (Munderloh et al., 1999), both genes of the rCRCT/CRCA module are expressed by R. buchneri (data not shown). The same primers did not amplify product in PCR reactions using DNA template from 10 diverse rickettsiae, illustrating their specificity and efficacy for future surveys of blacklegged tick populations ( Figure 6A ). This is demonstrated by testing the primers on a small sampling of blacklegged ticks from three different populations throughout the I. scapularis range ( Figure 6B ). Future sequencing of these loci will provide resolution on the conservation rCRCT/CRCA-1 and whether there is evidence for an arms race between R. buchneri that have diverged throughout the I. scapularis geographic range.

As R. buchneri REIS_1424 lacks a predicted signal sequence, it is likely secreted via one of two Sec-independent pathways: type I secretion system (T1SS) or type IV secretion system (T4SS), which are both conserved in rickettsiae (Gillespie et al., 2015b). Inspection of the neighborhood loci around select CRCT/CRCA modules indicates several diverse secretion pathways for some CRCT-DCPs, including T4SS (Sphingobium chlorophenolicum) and type VI secretion system (T6SS) (Sphaerisporangium krabiense), as well as typical CDI systems (Pseudomonas orientalis and Enterobacter sp. BIGb0383) with nearby cdiB loci, which encode the outer membrane CdiB that translocates the large CdiA protein as a type Vb secretion system (T5bSS) ( Figures S7A-G ). Using CdiB as a query in Blastp searches against the major taxa harboring CRCT/CRCA modules revealed their scarcity in Rickettsiales (none detected in rickettsiae), Actinobacteria, Chloroflexi, and Spirochaetes genomes, but widespread distribution in other proteobacterial and Cyanobacteria genomes ( Figures S7H, I ). This indicates that bacteria employing the CRCT/CRCA modules we describe here for warfare utilize a variety of secretions systems (e.g., T1SS, T4SS, T5bSS, T6SS, and likely others) consistent with the plethora of secretory pathways now characterized for diverse TA modules involved in interbacterial antagonism (Ruhe et al., 2018; Lin et al., 2020). The lack of cdiB genes and evidence that rCRCT toxins were originally appended to larger proteins indicates these modular Rhs toxins were once widespread in Rickettsia genomes. This is reminiscent of the large modular toxins we identified across numerous intracellular bacteria that encode a myriad of eukaryotic-like domains, some of which function in commandeering host reproduction (Gillespie et al., 2018). Like rCRCT/CRCA modules, many of these variable toxins are found adjacent to genes encoding probable antidotes, indicating a recapitulating theme for toxin architecture that persists evolutionarily and drives innovative strategies for colonizing eukaryotic hosts.