Comparative hologenomics of two Ixodes scapularis tick populations in New Jersey

Shotgun sequencing and taxonomic binning

Shotgun sequencing output and homology search results are summarized in Table 1; of the 24.0M and 25.9M respective NWSE and PVIL reads remaining after trimming, 97.1% and 97.2% had BLASTn hits against the ‘nt’ database + I. scapularis genome while 15.4% and 16.0% had BLASTx hits against the ‘nr’ protein database. Of the 751,985 and 810,827 respective assembled NWSE and PVIL metagenome scaffolds, 99.3% and 99.3% had BLASTn hits against the ‘nt’ database + I. scapularis genome while 22.8% and 22.1% had BLASTx hits against the ‘nr’ protein database. Far fewer matches were found when searching protein databases due to the fact that the majority of our genome data presumably represent non-coding host tick DNA, however this step was used as a secondary source of evidence for taxonomic binning as protein homology may remain conserved for metagenome constituents that do not have a sequenced genome or have diverged at the nucleotide level.

Examination of the BLASTn output (Fig. 1, Summarized in Table 2 (raw BLAST output and taxonomic annotations in Data S1–S8)) revealed a majority of the raw reads had top hits to ten genera: the proteobacterial genera Rickettsia, Anaplasma and Wolbachia; the spirochaete genus Borrelia; the Arthropod (e.g., host tick derived) genera Ixodes, Rhipicentor and Ornithodoros; the apicomplexan genus Babesia; and the nematode genera Loa and Brugia. Inspection of reads derived from tick genera other than Ixodes were found to be from conserved genomic loci for which equal-scoring top hits were obtained to multiple species. A similar pattern emerged when analyzing the assembled scaffolds: eight genera (Rickettsia, Anaplasma, Wolbachia, Borrelia, Ixodes, Babesia, Loa and Brugia) represented all but 14 of the 1.55M total scaffolds with BLASTn hits (Table 2, Data S3 & S4). By contrast, DIAMOND translated protein homology (i.e., BLASTX-like) searches of both reads and scaffolds contained hits to a broader diversity of taxa (Table S1; Data S5 & S8) although further inspection revealed this to be artifactual; a majority of non-coding host tick genomic DNA in the sample, when translated, had weak local alignments to proteins and/or mobile genetic elements of unrelated organisms in the ‘nr’ database. For example, a combined 271,697 reads from both locations had top hits to Ehrlichia minasensis (a tick-borne pathogen (Cabezas-Cruz et al., 2016)), however these hits comprised only 39 unique NCBI accessions for which the majority were reverse transcriptase and integrases, and all of which had perfect or near-perfect BLASTn alignments to the Ixodes scapularis nucleotide genome sequence thus making them likely retroelements. For this reason, we focused predominantly on the BLASTn results moving forward.

16S taxonomic diversity

Our 16S amplicon sequencing produced 418,528 paired reads from the NWSE library and 498,534 paired reads from the PVIL library. After quality trimming and primer removal, 376,944 (90.0%) and 461,306 (92.5%) respective NWSE and PVIL read pairs were overlapped into a single fragment prior to analysis (Table 1). All merged amplicons had BLASTn hits to the NCBI ‘nr’ database (Summarized in Fig. 1 and Table S2; raw BLAST output in Table S3). The taxonomic composition of these hits mirrored those of the shotgun metagenome sequencing, as 99.8% of the NWSE 16S fragments were assigned to the genera Rickettsia, Borrelia, and Wolbachia, while 99.90% of the PVIL fragments were assigned to the same three genera with the addition of Anaplasma.

Both shotgun and 16S amplicon sequencing identified the same set of bacterial taxa (Rickettsia, Borrelia, Wolbachia, and Anaplasma) as the major constituents, indicating (1) that the shotgun sequencing approach identified bacterial constituents despite the much greater abundance of tick DNA, and (2) that the two I. scapularis pools examined here contain few other microbial genera. Recent experiments employing surface sterilization, tissue dissection and/or rigorous statistical analyses have found that I. scapularis and I. pacificus (as well as species of Amblyomma and Dermacentor) tend to harbor only a limited number of permanent bacterial species, primarily endosymbionts (Lado et al., 2020; Pavanelo et al., 2020; Couper et al., 2019; Ross et al., 2018). This is reflected in the fact that genomes of many tick-associated endosymbionts and pathogens have lost suites of effector-immunity genes that mediate inter-bacterial competition, suggesting that these microbes encounter a limited community diversity within the tick (Ross et al., 2018).

Given the agreement between taxonomic affinity of shotgun genome data generated and 16S taxonomic assignment using BLASTn (the Wolbachia 16S BLASTn hits for example were >99.5% identical spanning a local alignment of >400 nt (Table S3)), we chose to present representative data from only the BLAST analyses herein.

Rickettsia endosymbiont

In each sample, the majority of non-tick reads (90.1% and 88.6% from NWSE and PVIL, respectively) had top hits to the proteobacterial genus Rickettsia (Table 2) and most likely derive from the known I. scapularis endosymbiont Rickettsia buchneri (Kurtti et al., 2015). The greatest number of Rickettsia reads hit the draft genome sequence of Rickettsia monacensis (NCBI accession LN794217), a spotted fever pathogen, while the majority of DIAMOND hits (167,949/365,340) were to Rickettsia endosymbiont of Ixodes scapularis (REIS; a synonym of Rickettsia buchneri (Kurtti et al., 2015; discussed below) Data S5 & S6), followed by R. amblyommatis and R. buchneri. We attribute this to exclusion of the R. buchneri symbiont genome from the ‘nt’ nucleotide database coupled with the presence of a predicted proteome in ‘nr’. To test this, both NWSE and PVIL libraries were mapped to the genomes of R. monacensis and R. buchneri str. ISO7 (NCBI accession GCF_000696365.2) concurrently, with non-specific hits (those mapping equally well to both genomes) discarded. Of the 200,410 and 298,528 reads aligned for each respective library, only 3,089 (1.5%) and 4,183 (1.4%) were to the R. monacensis genome. The reads that mapped to the R. buchneri genome had an average pairwise identity of 99.2%, while those mapping to R. monacensis had a pairwise identity of only 89.4% and are unlikely to be derived from this species. Taken together, these results allow us to conclude that the Rickettsia strain recovered in our Ixodes populations is R. buchneri.

Rickettsia endosymbionts are common among Ixodes spp. and there is evidence that they share a coevolutionary past (Díaz-Sánchez et al., 2019). Rickettsia buchneri has been documented at a high prevalence in many I. scapularis populations from the Northeastern US (e.g., 46.0% in Maine, 64.9% in Pennsylvania, 65%–100% in New York, and 93.3% in Connecticut) (Moreno et al., 2006; Pokutnaya et al., 2020; Steiner et al., 2014; Tokarz et al., 2019). Its occurrence is much higher in female I. scapularis than in males, localizing primarily within the ovaries where it is passed vertically to larvae (Kurtti et al., 2015; Zolnik et al., 2016). While endosymbionts are known to have a range of impacts on their tick hosts (Bonnet et al., 2017) information on the impact of R. buchneri infection on I. scapularis is sparse. However, some studies have indicated that it may provide them with folate (Hunter et al., 2015) and alter their motility (Kagemann & Clay, 2013). It is also possible R. buchneri could affect the susceptibility of I. scapularis to colonization by pathogens, as has been observed for R. bellii in Dermacentor andersoni (Gall et al., 2016).

To qualify our R. buchneri reads further, we mapped short-reads from both locations (NWSE and PVIL) to the Rickettsia endosymbiont of Ixodes scapularis str. Wikel genome sequence + 3 associated plasmids (NCBI BioProject PRJNA33979). We used this strain as opposed to str. ISO7 (above) as it was assembled to a single chromosome. Results indicate that 199,787 reads (ca. 21x avg. coverage of the reference) from the NWSE library and 299,315 reads (ca. 32x avg. coverage) from the PVIL library align to the symbiont genome. Using this short-read mapping, we detected and manually curated 19 single-nucleotide variants (SNVs; Table S6) present in the symbiont genome data; six variants occurred at >90% frequency in the NWSE sample (three were fixed at 100%) while present in <25% of the reads mapping to the equivalent locus in the PVIL library (three were absent entirely at 0%). This analysis also uncovered two diverged genomic regions in the symbiont genome of each sample, both of which occur in genes encoding conjugative transfer proteins. One region spans nucleotides 706,940–706,952 (RAGE-Be) of the symbiont chromosome scaffold CM000770 in which the PVIL read data support seven linked polymorphisms at high coverage that are absent from the NWSE data (Fig. S2). All seven variants are synonymous and occur within the conjugal transfer pilus assembly protein TraU gene. The second region spans nucleotides 760,409–760,434 (RAGE-F) of the genome scaffold and contains six linked SNVs that are fixed in the NWSE data but occur only at low coverage in the PVIL library (Fig. S3). Four polymorphisms are non-synonymous and result in four amino acid changes in the conjugative transfer protein TraA gene.

Both TraU and TraA are members of a 21-gene integrative conjugative element known as the Rickettsiales amplified genetic element (RAGE; Gillespie et al., 2012). RAGE elements together encode an F-like type IV secretion system allowing translocation of genetic material between diverse intracellular bacteria, and have been discovered on plasmid or chromosomal DNA of tick-borne Rickettsia bellii, R. massiliae and R. tamurae species, as well as R. peacockii and R. felis (Hagen et al., 2018). Rickettsia buchneri however is the only known species found to encode multiple RAGE elements (7 chromosomal and two encoded on plasmids) where they account for 9% of the genome by length (Gillespie et al., 2012; Hagen et al., 2018). The variants described above occur in two such elements, the TraU gene of RAGE-Be and TraA gene of RAGE-F, using the naming scheme of Gillespie et al. (2012). Both variation within and loss of the genes that constitute the RAGE-T4SS module have been reported between the two R. buchneri strains with currently sequenced genomes, as well as within natural tick populations (Hagen et al., 2018). To test whether the variants described here were mapping artifacts caused by the lack of a particular TraA/U element within the reference resulting in secondary read alignment to a sub-optimal location, we repeated the analysis using the R. buchneri strain ISO7 genome reference (Kurtti et al., 2015) with identical results (not shown). These data indicate concerted evolution acting on individual components of the RAGE-T4SS in I. scapularis symbiont genomes at the host (tick) population level.

Pathogens and other bacteria

In addition to the rickettsial endosymbiont, our metagenome sequencing also identified 2,140 NWSE and 3,887 PVIL reads with top hits to the spirochaete genus Borrelia (Table 2), with a majority (38.2% and 44.3% respectively) hitting B. burgdorferi strain MM1, originally isolated in Minnesota (Jabbari et al., 2018; Loken et al., 1985) and associated plasmids (Data S1 & S2). Short-read mapping aligned 1,841 and 3,248 reads to the B. burgdorferi strain MM1 genome. This causative agent of Lyme disease is known to be locally common in these tick populations from prior studies at the same sites, with prevalence ranging from 43.5–52.7% of adults infected (Dolan et al., 2016; Schulze et al., 2013). We did not detect B. miyamotoi in either pool, which is known to be rare across NJ (∼2.7% infected) (Egizi et al., 2018), nor B. mayonii, which has so far only been found in the Midwest (Kingry et al., 2017).

A second potential human pathogen was recovered from the PVIL library only, as 529 reads were found to have top hits to the proteobacterial genus Anaplasma (Table 2). Of these, 401 hits were to A. phagocytophilum strain Dog2 (an agent of Anaplasmosis originating from a MN dog (Al-Khedery et al., 2012); Genbank accession GCA_000439795, Data S1 & S2). No DNA corresponding to this pathogen was identified within the NWSE library. While A. phagocytophilum is also known to occur locally, its published infection prevalence in adult I. scapularis is substantially lower than that of B. burgdorferi, which could explain why it was not detected in both samples. Past work has shown that around 20–22% of adult I. scapularis are infected with A. phagocytophilum in Monmouth County (Dolan et al., 2016; Egizi et al., 2018). Two primary strains of A. phagocytophilum have been identified in I. scapularis, one pathogenic to humans (Ap-ha) and one that has not been associated with human illness and primarily circulates in white-tailed deer (Ap-V1), differing by two SNPs in the 16S rRNA gene (Krakowetz et al., 2014). While the single 16S read available in our dataset contained the Ap-ha associated SNPs, we caution that a single read could contain errors and is likely not definitive for strain assignment. However, we were able to verify that the Dog2 strain genome, the closest match in Genbank to our reads, is the Ap-ha strain.

The final bacterium with tangible representation in our data belonged to the proteobacterial genus Wolbachia, the top hit to 24 NWSE and 287 PVIL metagenome reads (Table 2), with 91.7% and 89.2% of respective hits to the Wolbachia symbiont of Cimex lectularius (the common bed bug) strain wCle genome sequence (Data S1 & S2). Wolbachia is commonly an endosymbiont of many hematophagous arthropods (e.g., mosquitoes) estimated to infect up to 60% of known insect species (De Oliveira et al., 2015). To determine the taxonomic placement of the Wolbachia found in this study, we constructed a maximum-likelihood phylogenetic tree using all available Wolbachia hits within the ‘nr’ database to the translated 152 amino acid ORF present on PVIL assembled scaffold NODE_632479, which encodes an elongation factor Tu (EF-Tu) protein that maintains 97.4% amino acid sequence identity over the entire 457 bp scaffold. This tree (Fig. 3) mirrors the branching order of Brown et al. (2016) and places the query sequence sister to Wolbachia strain wCle in a clade of obligate mutualists (clade F of Brown et al. (2016)). Zhang, Norris & Rasgon (2011) reported a novel Wolbachia species from Amblyomma americanum ticks collected in Maryland, and constructed a five-gene phylogeny that also placed this bacterium in an F-group clade containing the Cimex lectularius symbiont. We were unable to place our Wolbachia strain in this tree due to limited genome data recovery of appropriate homologs. At least one prior study posited that the Wolbachia pipientis strain detected in their Ixodes ricinus tick sample was due to cryptic presence of an Ixodiphagus hookeri endoparasitoid wasp (Plantard et al., 2012), however we did not detect hymenopteran DNA in our samples, nor was our Wolbachia strain closely related to W. pipientis.

If residing in the tick itself, the occurrence of Wolbachia is of particular interest due to known impacts of congeners on the reproductive fitness and vector competence of mosquitoes (Caragata, Dutra & Moreira, 2016). These discoveries have led to strategies of controlling mosquito-borne diseases by manipulating Wolbachia infectivity (paratransgenesis). Further characterization of tick-endosymbiont interactions could facilitate the use of paratransgenic strategies to control tick-borne diseases (Díaz-Sánchez et al., 2019).

Eukaryotic sequences

Multiple eukaryote constituents were identified within our metagenome data. Host tick (I. scapularis) DNA expectedly constituted an overwhelming majority of the assembly given the derivation of the tissue sample and large (2.22 Gbp) genome size of the tick with respect to other microbial and symbiont genomes present in the sample; when coupled with a relatively low-output sequencing run as described above, this would result in a fragmented assembly and large number of contigs as witnessed here.

Reads derived from the Apicomplexan parasite Babesia were identified in both samples (Table 2), with 669 hits in the NWSE library and 124 from the PVIL data. Babesia microti strain RI was the predominant taxon represented in both libraries (71.2% and 98.4% of each respective subset of reads). Like B. burgdorferi and A. phagocytophilum, Babesia microti is also frequently detected in Monmouth County ticks. In a prior field study, 8.2% of I. scapularis adults from PVIL and NWSE were positive for Ba. microti (Schulze et al., 2013). While a single NWSE amplicon (BMY3G:1:2110:24100:13578 (Supplemental File 1)) was identified with 100% homology exclusively to the Ba. microti apicoplast 16S rRNA locus and a further 489 reads had hits to Ba. microti with an average nucleotide similarity score of 99.5%, there was also a moderate proportion of hits to Babesia divergens (n = 148; Data S1) with a much lower average similarity score (86.9% nucleotide homology). A manual BLASTn into the publicly available and fully sequenced Ba. microti genome using these reads as a query returned no hits, and thus they appear to indicate a second species of Babesia is present in our dataset. Proper identification of these reads will likely require follow-up PCR with gene-specific primers, as there were no perfect or near-perfect matches in the NCBI ‘nt’ database and most Babesia species in public gene databases are represented by few (often non-overlapping) loci that prevent a broad phylogenetic analysis with homologs identified in our sequencing. One parsimonious possibility may be Babesia odocoilei, as this protozoal parasite of cervids is known to occur in I. scapularis ticks throughout the Northeastern US (Steiner et al., 2014) and for which a paucity of sequence exists in public databases.

Nematodes represent another eukaryotic component within our sample. We identified 59 and 1,694 reads with top hits to nematodes in the NWSE and PVIL libraries, of which 88.1% and 94.3% were within the family Onchoceridae (Table 2). Three amplicons (six paired reads) were identified with hits to nematode ribosomal RNA, of which one pair maintained 100% nucleotide homology over >300 nt to the small subunit of Filarioidea sp. ex Ixodes scapularis RTS-265 (NCBI accession MK868471.1) reported recently from I. scapularis ticks in New York and Connecticut (Tokarz et al., 2019). The remaining two amplicons were derived from the large subunit rRNA for which no data exist for the latter species, and thus the hits share lesser homology to various species.

Nematodes have previously been detected in Ixodes scapularis from Wisconsin (Cross et al., 2018) and Connecticut (Namrata et al., 2014) and appear to be closely related to the above accession from Tokarz et al. (2019) as well as to nematodes from A. americanum collected in Maryland and Virginia (Henning et al., 2016; Zhang, Norris & Rasgon, 2011). This clade of tick-associated nematodes is usually sister to Acanthocheilonema spp. (Namrata et al., 2014; Tokarz et al., 2019) or Monanema, when included, as in Henning et al. (2016); here their closest match is to Monanema spp. These detections are based on DNA sequencing and not physiological observation, thus their role in the tick (as potential commensals or tick-borne pathogens) remains unknown. In Europe, the brown dog tick Rhipicephalus sanguineus is an intermediate host of canine infectious filaria (Brianti et al., 2012; Olmeda-Garcia, Rodriguez-Rodriguez & Rojo-Vazquez, 1993; Otranto et al., 2012).

Both Tokarz et al. (2019) and Cross et al. (2018) attribute the Wolbachia they detected in their samples to filarial nematodes; while nematodes in the genus Acanthocheilonema have not been found to contain Wolbachia, there is evidence of a historical association: lateral gene transfer from Wolbachia into the Acanthocheilonema genome (Keroack et al., 2016; Mcnulty et al., 2010). Conversely, nematodes in the genus Monanema have not been found to carry Wolbachia (Ferri et al., 2011). Further work is needed to characterize this potentially novel clade of nematodes in US-based Ixodes and Amblyomma ticks and to determine whether the detected Wolbachia is carried by ticks or by their associated nematodes.