Ixodes scapularis microbiome correlates with life stage, not the presence of human pathogens, in ticks submitted for diagnostic testing

Morphological identification of ticks

Ticks used in this study were sent to Connecticut Veterinary Medical Diagnostic Laboratories (CVMDL) by Connecticut residents from across the state and the precise locations where the ticks were acquired by the humans are not available. The ticks used in this study were submitted to CVMDL in spring and summer of 2017 and 2018 (n = 31 and n = 14 respectively). They are removed prior to submission and are either dropped off directly at the laboratory or sent by mail. All ticks in this study were processed by the CVMDL ∼3–4 days post collection. Ticks submitted to CVMDL for diagnostic testing purposes were morphologically identified according to identification keys by Keirans & Litwak (1989) and US National Tick Collection (https://cosm.georgiasouthern.edu/usntc/) using a stereo microscope LAXCO LMSP-Z 230P (Mill Creek, WA) under power of magnification ranging from 6.7 to 45X.

Ticks were then classified based on engorgement (non-engorged, slightly, partial, fully); to control for variation from engorgement status we preferentially chose ticks that were either slightly or partially engorged (Ross et al., 2018). Adult ticks were sexed visually, but the nymphs were not, as nymphs cannot be sexed visually. All samples were extracted on the day they arrived at CVMDL.

Extraction of total DNA from tick specimens

DNA was extracted from individual ticks using a Nucleospin Tissue kit (cat# 740952-250, Macherey-Nagel GmbH & Co. KG, PA). Briefly, ticks were placed in 180 μL of buffer T1or in 360 μL of the same buffer for engorged specimens. Sterile sand was then added to the tube and tick bodies were manually homogenized using sterile disposable wood applicators. After homogenization, 25 µl of proteinase K were added to the mix followed by an incubation step at 56 °C for 4 h. The DNA extraction procedure proceeded then as recommended by the manufacturer. Positive extraction control (PEC) tick was included with each tick DNA extraction, the PEC tick was selected to be free from all tested pathogen. Ticks are processed and tested as they are received; including negative extraction controls on diagnostic samples is not in the standard protocol, as it is not a significant source of error for pathogen detection. To detect potential contamination from extraction, we followed De Goffau et al. (2018) and tested for batch variation based on dates of extraction. We also confirmed via a thorough literature search that the taxa detected in our samples could reasonably have come from ticks (i.e., are ecologically plausible).

qPCR for detection of pathogens

SYBR green-based qPCRs were used to detect genomic DNA of pathogens (see Table 1 for primers used and references). Fast SYBR Green Master Mix (cat# 4444432, Applied Biosystems, Foster City, CA) was used in 25 µl reaction following manufacturer’s recommendations. DNA was amplified using an ABI 7500 Fast Real-Time PCR thermal cycler (Applied Biosystems, Foster City, CA) by an initial denaturation step hold at 95 °C for 3 min followed by 45 cycles at 95 °C for 30 s, 60 °C for 30 s, and 72 °C for 30s. Signal was captured at each cycle at the end of elongation steps. Melting curve analysis was performed at the end of each run for 95 °C for 15s, 60 °C for 1 minute, 95 °C for 15s and 60 °C for 15s. For each real time PCR run, a Non-Template Control (NTC) was used. The NTC produced negative result for each run showing there was no external contamination. A Positive Amplification Control was used in each run as well. A sample is considered positive for a particular pathogen when an amplification curve is observed (Ct value obtained). Each sample was tested for B. burgdorferi, A. phagocytophilum, B. microti, and B. miyamoti, see Table S1 for each sample’s pathogen Ct scores.

Sample selection and 16S rRNA gene sequencing

We selected samples with between zero and three common TBPs (Table 2). One adult tick had three TBPs (A. phagocytophilum, B. burgdorferi, and B. microti, referred to as “[Bb, Ap, Bab]”), seven had two TBPs (three adults: B. burgdorferi, A. phagocytophilum, “[Bb, Ap]”, two adults: B. burgdorferi, B. microti “[Bb, Bab]”, two nymphs: A. phagocytophilum, B. burgdorferi, “[Bb, Ap]”), 25 had one TBP (Table 2). Negative sequencing controls were used to verify the lack of contaminants in reagents used during amplification and sequencing.

All of the adults used for this study were females and the sexes of the nymphs are unknown. We were unable to precisely control the engorgement levels; however, most ticks selected were either partially or slightly engorged. One was engorged (adult) and two were not engorged at all (one nymph and one adult).

DNA extracts were amplified and sequenced at the University of Connecticut Microbial Analysis, Resources and Services center using the standard protocol for amplification and sequencing. Extracts were quantified using the Quant-iT PicoGreen kit (Invitrogen, ThermoFisher Scientific). Partial bacterial 16S rRNA genes (V4, 0.8 picomole each 515F and 806R with Illumina adapters and 8 basepair dual indices (Kozich et al., 2013) were, in triplicate, amplified in 15 ul reactions using GoTaq (Promega) with the addition of 10 µg BSA (New England BioLabs). We added 0.1 femtomole 515F and 806R which does not have the barcodes and adapters to overcome initial primer binding inhibition, because the majority of the primers do not match the template priming site. The PCR reaction was incubated at 95 ° C for 2 minutes, the 30 cycles of 30 s at 95.0 °C, 60 s at 50.0 °C and 60 s at 72.0 °C, followed by final extension as 72.0 °C for 10 minutes. PCR products were pooled, quantified and visualized using the QIAxcel DNA Fast Analysis (Qiagen). PCR products were pooled using QIAgility liquid handling robot after the products were normalized based on the concentration of DNA from 350-420 bp. The pooled PCR products were cleaned using the Mag-Bind RxnPure Plus (Omega Bio-tek) according to the manufacturer’s protocol. The cleaned pool was sequenced on the Illumina MiSeq using v2 2 × 250 base pair kit (Illumina, Inc). Three PCR controls were also sequenced to test for PCR reagent contamination.

Sequence processing

All data analyses and sequence processing were done using R v.3.5.0 (R Core Team, 2019). Sequences were quality controlled, denoised and merged using DADA2 (Callahan et al., 2016) to create a sample by ASV (amplicon sequence variant) matrix. An ASV is an operational taxonomic unit, defined as any unique sequence that passes stringent quality control. Taxonomy of ASVs was assigned using RDP’s Naïve Bayesian Classifier with the Silva reference database v128 (Wang et al., 2007; Quast et al., 2013). Sequences that were identified as mitochondria, chloroplasts, or that were unable to be confidently assigned to any bacterial phylum were removed. We also removed 22 ASVs that were identified as Candidatus Carsonella, an endosymbiont found only in psyllids (Sloan & Moran, 2012). We removed these ASVs because they were nested within the mitochondria clade in our original analysis. We further verified that these 22 ASVs were likely mitochondrial in origin using blastn. We conducted a filtered search, only keeping matches with >90% sequence identity and e-value < 1e-40. We found 15 ASVs were not assigned to Candidatus Carsonella but instead, matched best with I. scapularis and I. pavlovskyi mitochondrial sequences see supplemental Table S2 (Zhang et al., 2000; Morgulis et al., 2008). The remaining 7 ASV sequences could not be assigned to any organism using the previous search parameters.

To remove likely contaminants, we processed the sequences using the Decontam package (Davis et al., 2018), which uses the negative controls to identify likely contaminants.

To calculate phylogenetic diversity metrics, we performed a multiple-alignment of all ASVs using the DECIPHER package in R (Wright, 2015), and constructed a phylogenetic tree with the phangorn package version 2.4.0 (Schliep, 2011).

Data analyses

Data analysis was done using the R packages phyloseq (McMurdie & Holmes, 2013), vegan (Oksanen et al., 2007), and DESeq2 (Love, Huber & Anders, 2014). Code was generated from vegan and phyloseq tutorials available online. We rarefied samples to an even depth of 8,252 reads for alpha and beta diversity analyses, which resulted in the loss of two samples from the dataset for being below 8,252 reads. See Fig. S1 for rarefaction curves.

Shannon diversity index was used to quantify the alpha diversity of the samples. Using a Shapiro test, we determined the data were not normally distributed; therefore, we tested the significance of life stage of ticks using a Kruskal–Wallis test. The alpha diversity was first tested on the rarefied data; however, low biomass samples have the potential to artificially inflate diversity (Salter et al., 2014; Erb-Downward et al., 2020). To address this, we removed ASVs that comprised less than 0.1%, 1%, 5%, and 10% of the overall sequence count in each sample, rarefied again to 8200 reads, and tested alpha diversity using the Kruskal–Wallis test between adults and nymphs (Couper et al., 2019). We conducted multiple pairwise comparisons of alpha diversity looking at different TBP conditions using Wilcoxon rank-sum test.

Beta diversity was assessed by generating pairwise distance matrices of Bray–Curtis, weighted Unifrac, unweighted Unifrac, and Jaccard distances (Lozupone et al., 2007; Xia & Sun, 2017). We plotted NMDS ordinations using phyloseq and tested for significance of the metadata using a Permutational multivariate analysis of variance (PERMANOVA). We tested how the different TBPs correlate to the microbiome. The data were tested by comparing the samples when grouped by their combination of TBPs, refer to Table 2 for TBP combinations. The samples were also categorized by whether they had one of the four pathogens regardless of whether that sample had additional TBPs. We did this in order to see if one TBP correlated to significant variation in the microbiome. Beta diversity was tested using PERMANOVA. We tested life stage by merging the [TBP-] and [Bb] samples when we found no significant difference in those samples’ microbiome, we did this to increase the sample size of both nymphs and adults. Furthermore, we tested for temporal effects on the microbiome and tested changes in beta diversity using the previously mentioned distance metrics. We tested adults and nymphs separately using the whole data set with sequences <1% removed and rarefied to 8200 reads.

To test if TBPs could influence diversity results, we tested using the same methods as before, but we removed all TBP reads and rarefied to a new sequencing depth (3,204 reads). We then tested these results as described previously. To confirm there was no difference, we rarefied the original data set to 3,204 reads with the TBP still present and ran a PERMANOVA on those data again to ensure results were consistent across rarefication depths.

We identified differential abundance of genera from different life stages and different TBPs using DESeq2 (v. 1.24.0) in R. Genera were identified as differentially abundant if the corrected p < 0.05. p-values were corrected with the Benjamini and Hochberg false discovery rate for multiple testing (Benjamini & Hochberg, 1995) .

Relative abundance of taxa in our samples were calculated after the samples had been rarefied (8,252 reads). Low abundance taxa in each sample were grouped if they totaled <1% of the reads at the phylum level and <15% of the genera. The relative abundances for each taxon in each sample were averaged to get the relative abundance for each TBP condition and life stage.

Sequence data and homogeneity check

The initial dataset had 2,449,198 sequences and the final dataset had 2,114,742 high-quality sequences. Sequences per sample ranged from 5,018–126,549 (mean = 46,994). No sequences were identified as likely contaminants.

Homogeneity of variance of the rarefied data was tested using ANOVA and Bray–Curtis distances. We saw no significant difference in life stage or pathogen status (life stage ANOVA: F 1, 4 1 = 0.4685, p = 0.4975; pathogen status ANOVA: F 7, 3 5 = 2.1431, p = 0.06442).

Ixodes scapularis microbiome: bacterial composition

The mean relative abundances of the phyla dominating the [TBP-] samples were Proteobacteria (78.9%), Firmicutes (10.89%), and Spirochaetes (7.4%). [Bb] ticks were similar to [TBP-]: Proteobacteria comprised the majority (77.65%), Firmicutes (14.26%), and Spirochaetes (6.5%). [TBP-] samples had 8 phyla that were >1% abundances in at least one sample. [Bb] had five phyla that were >1% of the relative abundance (Table S3). The genera in both [TBP-] and [Bb] were largely similar and dominated by Rickettsia (43.3% and 43.1% respectively), Fig. 1. [TBP-] and [Bb] had almost equal amounts of rare taxa (relative abundance < 15%): 19.4% and 18.7%, respectively. Complete table relative abundance of genera in [TBP-] and [Bb] samples available, see supplemental Table S4.

Effect of B. burgdorferi infection in I. scapularis microbiome composition

For clarity, the first analyses were between just the [TBP-] and [Bb] ticks. Alpha diversity was not significantly different between [Bb] (n = 13) and [TBP-] (n = 10), Shannon: X² = 0.0038462, df = 1, p = 0.9505, (Fig. 2A). NMDS ordinations showed minimal clustering of samples by the presence of B. burgdorferi (Fig. 2B). Infection status was not significant (p > 0.05) using three of the four distance metrics, although infection status explained 8% of the variation in the unweighted Unifrac NMDS (df = 1, R⁼0.08, p = 0.024, Table 3).

TBPs, single and co-infections effect on the microbiome

We next included all samples, including [TBP-], single infection and coinfected ticks. The most abundant phylum, regardless of TBP status, was Proteobacteria (77.7–98.2%). Firmicutes ranged from 1.5–14.3%, except for in [Bb, Ap, Bab] and [Bb, Bab], where no Firmicutes were present above 1% of the relative abundance. Spirochaetes were in all samples except TBP statuses [Ap] and [Bab]; in the remaining TBP statuses, Spirochaetes ranged from 1.1–13.79% (Table S3). Rickettsia was the most abundant genera in all TBP statuses except for [Ap] and [Bab]. [Ap] had similar levels of Rickettsia and Anaplasma (25.4% and 27.2% respectively). In [Bab], Pseudomonas comprised 53.7% of the microbiome. Anaplasma was also in high abundance in [Bb, Ap] with Rickettsia (48.9%) and Anaplasma (36.2%). Only [Ap] and [Bb, Ap] had Anaplasma with relative abundance greater than 15%. Borrelia average abundance in [Bb] was 4.2%, [Bm] 10.9%, and [TBP-] 7.4% (Fig. 1).

The number of TBPs in a tick did not correlate with beta diversity in I. scapularis. We compared samples that had one, two or three TBPs present and regardless of distance metric, there was no significant correlation between infection status (the number of TBPs) and the microbiome (Bray–Curtis: df = 7, R² = 0.15833, p = 0.586; weighted Unifrac: df = 7, R² = 0.1545, p = 0.598; unweighted Unifrac: df = 7, R² = 0.19906, p = 0.127; Jaccard: df = 7, R² = 0.15858, p = 0.604). To see if the presence of the pathogens themselves were altering the diversity of the microbiome, we removed the TBPs ASVs and beta diversity was still not significantly correlated to TBPs (Table 4). Bray-Curtis NMDS ordinations with and without TBPs ASVs showed minimal clustering (Figs. 3A–3B). No single TBP had a significant effect on the beta diversity in I. scapularis (Table 4), except B. burgdorferi, which was significant when using unweighted Unifrac (R² = 0.06849, p = 0.002) (Table 4).

Shannon diversity indices indicated no significant difference in alpha diversity of any of the TBP statuses (p > .05) with the exception of [Bb, Ana] and [Bm] (Wilcoxon rank-sum: p = 0.032) (Fig. 4). We then separated samples into either [Bb], [Ap], [Bm], or [Bab] regardless of multiple TBPs and we found there were no differentially abundant genera. The only genera that were differentially abundant in the samples were the TBPs in question. Ticks classified as [Bb] or [Bm] had Borrelia significantly more abundant. Similarly, [Ap] ticks had significantly more Anaplasma. [Bab] is a eukaryote and hence there were no Babesia sequences to test.

There was no significant correlation (p >0.05) between the month and year that the tick was collected and beta diversity except when using Jaccard distance with the nymphs (df = 8, p = 0.027, R² = 0.56164).

Microbiome diversity and I. scapularis life stage

Ixodes scapularis life stage correlates with microbiome composition. Alpha diversity differed significantly between adults and nymphs (Kruskal–Wallis, Shannon index: X⁼11.978, df = 1, p = 0.0005384, Fig. 2A). Shannon diversity index in adults was 0.843 ± 0.148 and was 2.241 ± 0.305 in nymphs. Accounting for artificial inflation of alpha diversity in low biomass samples; adults and nymphs were still significantly different when low abundance ASVs (less than 0.1% or 1%) were removed. When ASVs that comprised less than 5%, and 10% of the total ASVs were removed, we no longer saw significant differences between life stages (Table S5). Life stage was significantly correlated with Bray–Curtis, weighted Unifrac, unweighted Unifrac, and Jaccard distance metrics (Table 3). Differences in life stage can be visualized in Bray–Curtis NMDS ordination (Fig. 2B). Life stage explained between 11–15% of the variation observed (Table 3).

Nymphs had eight phyla >1% relative abundance: Proteobacteria (68.0%), Firmicutes (14.8%), and Spirochaetes (13.5%) accounted for most of the communities. Adults had four phyla with >1% relative abundance, with the largest being Proteobacteria (87.5%) and Firmicutes (10.9%) (Table S6). When averaged, adult ticks had seven genera with relative abundance >15%; whereas, nymphs had eight. Rickettsia was more abundant in adults than in nymphs (66.1% vs 16.8%). The averaged relative abundance of genera in nymphs was not dominated by a singular taxa. The most abundant genera identified in nymphs belonged to an unknown genus in Proteobacteria (17.7%), Rickettsia (16.8%), and Borrelia (11.7%) (Fig. 5). Rickettsia was absent in three adult samples; the rest had Rickettsia present at higher than 74.2%. The three adult samples that did not have Rickettsia had a similarly large amount of Lysinibacillus (two samples 46.0% and 64.3%) or Klebsiella (54.3%). Half of the nymphs had Rickettsia; one sample was similar to adults (82.6% relative abundance) the other five ranged from 14.7–28.7%. Five nymphs did not have Rickettsia and have different genera predominating (Table S7).