Skip to main content

Range-wide genetic analysis of Dermacentor variabilis and its Francisella-like endosymbionts demonstrates phylogeographic concordance between both taxa

Abstract

Background

The American dog tick, Dermacentor variabilis, is an important vector of pathogens to humans, wildlife and domestic animals in North America. Although this tick species is widely distributed in the USA and Canada, knowledge of its range-wide phylogeographic patterns remains incomplete.

Methods

We carried out a phylogenetic analysis of D. variabilis using samples collected from 26 USA states and five Canadian provinces. Tick samples (n = 1053 in total) originated from two main sources: existing archives (2000–2011), and new collections made from 2012 to 2013. We sequenced a 691 bp fragment of the cox1 gene from a subset (n = 332) of geographically diverse D. variabilis. DNA extracted from individual ticks (n = 1053) was also screened for a Francisella-like endosymbiont, using a targeted 16S rRNA sequencing approach, and important pathogens (Rickettsia spp. and Coxiella burnetii), using species-specific quantitative PCR assays.

Results

Maximum parsimony analysis of cox1 sequences revealed two major groups within D. variabilis with distinct geographical distributions: one from the eastern USA/Canada (Group 1) and one from the west coast states of the USA (California and Washington; Group 2). However, genetic subdivisions within both of these two major groups were weak to moderate and not tightly correlated with geography. We found molecular signatures consistent with Francisella-like endosymbionts in 257 of the DNA extracts from the 1053 individual ticks, as well as Rickettsia spp. and Coxiella burnetii in a small number of ticks (n = 29 and 2, respectively). Phylogenetic patterns for Francisella-like endosymbionts, constructed using sequence data from the bacterial 16S rRNA locus, were similar to those for D. variabilis, with two major groups that had a nearly perfect one-to-one correlation with the two major groups within D. variabilis.

Conclusions

Our findings reveal a distinct phylogenetic split between the two major D. variabilis populations. However, high levels of genetic mixture among widely separated geographical localities occur within each of these two major groups. Furthermore, our phylogenetic analyses provide evidence of long-term tick-symbiont co-evolution. This work has implications for understanding the dispersal and evolutionary ecology of D. variabilis and associated vector-borne diseases.

Background

The American dog tick (Dermacentor variabilis) is the most widely distributed native three-host tick in North America and is an important vector of pathogens to humans, wildlife, and domestic animals across the continent [1, 2]. It plays a critical role in the natural maintenance and spread of multiple tick-borne pathogens, such as Spotted Fever group Rickettsia spp. (causing Rocky Mountain spotted fever and other spotted fever group rickettsioses), Francisella tularensis (causing tularemia), Coxiella burnetii (causing Q-fever), and Anaplasma spp. (causing anaplasmosis) [3]. Dermacentor variabilis is abundant in the eastern USA and southeastern Canada, with a more restricted western population that occurs in the Pacific coast states of the USA; additional disjunct populations have also become established in restricted localities in the intermountain USA and Canada (Alberta, British Columbia, Idaho and Montana) [2, 4,5,6,7]. Overall, regional climatic and geographical differences (e.g. temperature and elevation) appear to be the prevailing factors that affect the distributional limits of D. variabilis and it has been hypothesized that these factors, along with associated physiological limitations of the tick, could result in genetic divergence of populations [8]. Interestingly, a westward expansion of D. variabilis has been documented in recent decades, highlighting the potential for range expansion due to climatic shifts across this region [6, 7, 9].

In the tick gut, endosymbionts and other non-pathogenic microbes significantly outnumber pathogenic organisms [10]. The nature of the relationship between ticks and their symbionts remains understudied and, therefore, poorly understood. Non-pathogenic endosymbiotic bacteria are described from D. variabilis and the abundance and composition of this microbiome can vary considerably among populations, even in those that are separated by relatively short geographical distances [5, 11]. The presence of these tick-associated symbionts has been shown to affect the acquisition of some pathogenic bacterial species [12] and, therefore, influence pathogen distributional patterns. A recent study in a very closely related species, D. andersoni, demonstrated that there was a negative correlation between the presence of the rickettsial symbiont R. bellii and the level to which the tick-borne pathogen A. marginale will replicate [13]. Furthermore, the “Rickettsial Exclusion Hypothesis” described the prevention of transovarial transmission of R. rickettsii when non-pathogenic endosymbiotic rickettsiae were prevalent [14, 15]; additional studies have experimentally confirmed this rickettsial interference phenomenon with other spotted fever group rickettsiae [16, 17]. The possibility that symbionts might influence competence for tick-borne pathogens has led to an increased interest in the study of endosymbiont prevalence and distribution [18,19,20].

Genetic analyses of D. variabilis are important for understanding its natural history as a host and vector, and to reconstruct its evolution and dispersal mechanisms. This important disease vector boasts an extensive and expanding range, which emphasizes the need for comprehensive phylogeographic studies to better understand its distributional patterns and characterize the relationship between these ticks and their microbiota. Multiple population genetic studies of D. variabilis have identified important phylogeographic trends for this species [5, 7, 21,22,23] and the bacterial pathogens that they vector; however, the breadth of these studies has been limited by geographical scope (1 to 3 states or provinces).

Here, we address this gap in knowledge with a phylogeographic study of range-wide samples using mitochondrial sequencing. The mitochondrial 16S rRNA and cytochrome c oxidase subunit 1 (cox1) genes have been used for phylogenetic and population analyses of many tick species in the past [24,25,26,27,28,29,30,31] and have been shown to be phylogenetically informative in discerning relationships from the subfamily level down to that of population. The aims of this study were to: (i) determine if D. variabilis exhibits broad population structure throughout its range using sequence data from a large number of individuals and populations representing the known species distribution, and (ii) examine whether there is an association between the presence and population structure of Francisella-like endosymbionts and the patterns observed within D. variabilis. We also screened all D. variabilis samples (n = 1053) for Rickettsia spp. and Coxiella burnetii to better assess the relationship between these two infectious agents and this important arthropod vector.

Methods

Sample collection and DNA extraction

To assess phylogeographic patterns of Dermacentor variabilis ticks throughout the known geographical range of this species, we obtained individual ticks from two separate sources: (i) existing archives maintained by collaborators (2000–2011), and (ii) new collections made from 2012–2013. For the new collections we provided a sampling kit to collaborators, which included a data sheet for recording sample information, such as host, date of collection, and geographical location. Donated ticks were identified using a morphological dichotomous key [1] and stored at -80 °C until further processing.

Genomic DNA was extracted from 1053 individual whole ticks using the 96-well DNeasy® Tissue Kit (Qiagen, Valencia, CA, USA) modified for use with the QIAvac vacuum filtration system. To begin, individual whole ticks (adult ticks or nymphs) were placed into deep-wells (1.2 ml) of a 96-microtube plate (Qiagen) along with 12 Zirconia beads (Catalog # 7361-002000, 2 mm diameter, Glen Mills, Inc., Clifton, NJ, USA). Lysis solution was prepared and added to the wells with the following protocol modifications: 360 μl Buffer ATL, 40 μl Proteinase K, 100 μl Antifoam-A (Sigma, St. Louis, MO, USA), for a total of 500 μl per well. Each 8-well column was then sealed using round-capped collection microtube caps (Catalog # 19566, Qiagen, Valencia, CA, USA) and placed in a secondary containment heat-sealed 9 × 12 inch Kapak bag (Kapak, Minneapolis, MN, USA) for biosafety purposes. The sealed extraction plates were then placed in a Geno/Grinder 2000 mill homogenizer (SPEX CertiPrep, Metuchen, NJ, USA), which was operated at 2000 strokes/min for 10 min [32]. Following bead-milling, plates were incubated overnight at 56 °C and then subjected to 10 min of additional bead milling the following morning. The remainder of the protocol was carried out in a QIAvac filtration system according to the manufacturer’s guidelines, with a final elution volume of 120 μl (60 μl, twice).

Mitochondrial cox1 sequencing

Many of our sampling locations were represented by an abundance of ticks (n > 30 per individual collection) so, to avoid sampling bias, we chose to investigate phylogeographic structure in the mitochondrial cox1 gene using a subset (n = 332) of our 1053 D. variabilis samples from the USA and Canada. To accomplish this, no more than 30 ticks were included from any one sampling location. We designed cox1 primers specific to D. variabilis that amplified a 793 bp fragment of this gene [COI-F2 (5'-CTT AAT TTT CGG CAG TTG AGC A-3') and COI-R1 (5'-GTT CTT TTT TTC CTG TGG AAA AAC-3')]. PCRs were carried out in 20 μl volumes, using 2 μl of extracted genomic DNA (diluted to 4 ng/μl) as template and the following reagents (given in final concentrations): 1× PCR buffer, 2.5 mM MgCl2, 0.2 mM dNTPs, 0.16 U/μl Platinum® Taq polymerase (Invitrogen, Carlsbad, CA, USA), and 0.4 μM of each primer. PCRs were thermocycled according to following conditions: 95 °C for 10 min to release the polymerase antibody, followed by 40 cycles of 94 °C for 60 s, 58 °C for 30 s, and 72 °C for 30 s. PCR products were then treated with ExoSAP-IT (Affymetrix, Santa Clara, CA, USA) using 1 μl of ExoSAP-IT per 5 μl of PCR product under the following conditions: 37 °C for 15 min, followed by 80 °C for 15 min. Treated products were then diluted in the range of 1:2 to 1:20 depending on amplicon intensity (as determined by agarose gel electrophoresis) and sequenced in both directions using BigDye® Terminator v3.1 Ready Reaction Mix (Applied Biosystems, Foster City, CA, USA) with the same forward and reverse primers from the initial PCR. We used 10 μl volumes for sequencing reactions containing the following reagents (given in final volumes): 3 μl of 5× sequencing buffer, 1 μl BigDye® Terminator v3.1 Ready Reaction Mix, 1 μl of a 10 μM primer stock, and 5 μl diluted PCR product. The following thermocycling conditions were used: 96 °C for 20 s, followed by 30 cycles of 96 °C for 10 s, 50 °C for 5 s, and 60 °C for 4 min. An ethanol precipitation technique was used to clean and precipitate the DNA pellet, and Sanger sequencing was carried out using an AB3130xl® automated genetic analyzer (Life Technologies, Grand Island, NY, USA). Sequence chromatograms were edited manually in Sequencher 5.0 (Gene Codes, Ann Arbor, MI, USA); after trimming poor quality sequences from each end we were able to use 691 bp of the cox1 fragment (from 793 bp in the PCR amplicon). Single nucleotide polymorphisms (SNPs) were identified in MEGA version 6 [33] and a maximum parsimony analysis was employed to construct a phylogenetic tree. The tree was rooted with a D. reticulatus cox1 sequence (GenBank: AF132829.1) and a bootstrap analysis of 1000 replicates was used to estimate the robustness of nodes on the tree; bootstrap values were only reported when found to be 50% or greater. To assess the relative amount of homoplasy, the consistency index (CI) and retention index (RI) were calculated following parsimony analysis in MEGA version 6 [33]. The consensus tree was exported into FigTree v1.4.2 (http://tree.bio.ed.ac.uk/software/figtree/), and then into Adobe Illustrator CS4 (Adobe Systems Incorporated, San Jose, CA, USA) for annotation purposes.

We examined the relationship between genetic distance (generated from the 332 cox1 sequences) and geographical distance in D. variabilis (based upon collection locations) by calculating Spearman’s correlation coefficient using the RELATE function in PRIMER-E software version 6 (Plymouth Marine Laboratory, Plymouth, UK). Analyses were carried out for three datasets: (i) the D. variabilis population as a whole (n = 332); (ii) individual D. variabilis (n = 300) assigned to Group 1 in the cox1 phylogeny; and (iii) individual D. variabilis (n = 32) assigned to Group 2 in the cox1 phylogeny. We tested the null hypothesis of no relationship between the genetic and geographical distance matrices (ρ = 0) by generating a simulated distribution of randomly sampled data pairs (one from the geographical distance matrix and one from the genetic distance matrix) taken from 10,000 permutations. Additionally, we calculated both within- and between-group mean genetic distances using a maximum composite likelihood function in MEGA version 6 [33].

Mitochondrial 16S rRNA sequencing

To assess phylogenetic congruence between two commonly used molecular markers and ensure that conclusions drawn from the cox1 gene were supported by an additional mitochondrial locus, we also investigated phylogeographic structure among a small subset (n = 32) of our D. variabilis samples from the USA and Canada by sequencing a 453 bp fragment of the tick mitochondrial 16S rRNA gene using primers 16S+1 (5'-CCG GTC TGA ACT CAG ATC AAG T-3') and 16S-1 (5'-CTG CTC AAT GAT TTT TTA AAT TGC TGT GG-3') [25]. PCR conditions were identical to the cox1 gene, as described above, except an annealing temperature of 58 °C for 60 s was used. Amplicons were visualized, cleaned, and sequenced as described above for cox1, but using the primers 16S+1 and 16S-1. Sequence editing, SNP identification, phylogenetic analysis, and figure processing were also conducted as described above for the cox1 gene, using a D. andersoni sample (GenBank: AY375431) as the outgroup.

Francisella-like endosymbiont and pathogen screening

To understand the relationship between D. variabilis and specific bacterial species, we screened all samples (n = 1053) for the presence of Francisella-like endosymbionts using Sanger sequencing, and additionally, several pathogens using species-specific real-time qPCR assays. First, we targeted the bacterial 16S rRNA locus of Francisella-like endosymbionts using conventional PCR with primers developed specifically for these endosymbionts [20, 34, 35]. Because endosymbiont abundance within a tick can be low, leading to false negative results, all individuals were screened in duplicate to increase the likelihood of detection. Primers F11 (5'-TAC CAG TTG GAA ACG ACT GT-3'; Escherichia coli 16S positions 150–169) [34] and 1227R (5'-CCA TTG TAG CAC GTG T-3'; E. coli 16S positions 1227–1242) [35] were used to generate a ~1093 bp amplicon (1090 bp in Francisella due to two separate deletions). A standard 20 μl PCR was performed with 2 μl (diluted to 20 ng/μl) of the extracted genomic DNA as template and the following reagents (given in final concentrations): 1× PCR buffer, 2.5 mM MgCl2, 0.2 mM dNTPs, 0.16 U/μl Platinum® Taq polymerase, and 0.4 μM of each primer. PCRs were thermocycled using the following conditions: 95 °C for 10 min to release the polymerase antibody, followed by 40 cycles of 94 °C for 60 s, 56°C for 30 s, and 72 °C for 1 min 30 s. In the absence of Francisella-like endosymbionts, this 16S assay occasionally amplified other bacterial species. For this reason, all 16S positive PCR products were then sequenced to confirm the identity and presence of Francisella-like endosymbionts. PCR visualization, cleaning, and sequencing methods were the same as described above for the cox1 gene but used the primers F11 and 1227R. Sequences were trimmed on each end to a final length of 725 bp (E. coli 16S positions 253–972). Sequence editing, SNP identification, phylogenetic analysis, and figure processing was conducted as described above, with a F. tularensis subsp. holarctica sample (GenBank: MG834502) as the outgroup.

All D. variabilis DNA extracts (n = 1053) were also screened for the presence of Rickettsia species and C. burnetii DNA using species-specific qPCR assays. To facilitate rapid analysis of all ticks (n = 1,053), preliminary screening for Rickettsia species was carried out using the Pan-Rickettsia real-time PCR assay [36], a TaqMan® qPCR assay that targets the 23S rRNA gene using “PanR8” oligonucleotides. This real time assay is able to detect a number of Rickettsia species, including R. rickettsii, R. prowazekii and R. typhi, all of which have been reported in North America and cause human infections [37, 38]. Reactions were carried out in 10 μl volumes using 1 μl of extracted genomic DNA (pre-diluted to 4 ng/μl) as template and the following reagents (given in final concentrations): 2× AB TaqMan® Universal master mix (Life Technologies, Grand Island, NY, USA), 1 μM of both the forward and reverse primer, and 0.04 μM probe. Plates were analyzed using an ABI 7900HT Sequence Detection System (Life Technologies, Grand Island, NY, USA); two positive controls (mean CT 19.8) along with 6 NTCs (no amplification) were included with each run. To obtain species level identification of all Rickettsia positive samples from the previous qPCR (CT < 40, n = 29), a 530 bp segment of the OmpA gene was amplified and sequenced. Primers Rr190.70p (5'-ATG GCG AAT ATT TCT CCA AAA-3') and Rr190.602n (5'-AGT GCA GCA TTC GCT CCC CCT-3') [39] were used in a 25 μl PCR, performed with 3 μl of the extracted genomic DNA (pre-diluted to 20 ng/μl) and the following reagents (given in final concentrations): 1× 5 PRIME MasterMix (5 PRIME, Hilden, Germany) and 0.4 μM of each primer. PCRs were thermocycled according to the following conditions: 94 °C for 3 min to release the polymerase antibody, followed by 32 cycles of 94 °C for 45 s, 56 °C for 1 min, and 72 °C for 1 min 30 s. PCR products were visualized, cleaned, and sequenced as previously described, but using the primers Rr190.70p and Rr190.602n. Sequence editing, SNP identification, phylogenetic analysis, and figure processing was conducted as described previously. Additionally, we screened for the presence of C. burnetii among our D. variabilis individuals using the IS1111 TaqMan® real-time qPCR genotyping assay, following the conditions set forth by Loftis et al. 2006 [40]; six positive controls (mean CT 9.2) along with 12 NTCs (no amplification) were included with each run.

Results

Mitochondrial cox1 sequencing

We identified two major groups within D. variabilis (bootstrap support = 100%) that separate according to distinct geographical distributions (Fig. 1 and Additional file 1: Table S1): one occurs in the eastern USA/Canada (Group 1) and the other is restricted to west coast states in the USA (CA, WA; Group 2). This result is in agreement with the major pattern observed in a previous study that investigated the mitochondrial 16S gene in a limited number of D. variabilis collections from the western USA and Canada [5]. We found a large amount of genetic diversity within the 691 bp cox1 fragment: Group 1 was comprised of 118 unique mitochondrial haplotypes, whereas Group 2 had just 12 haplotypes. In total, 110 nucleotide positions were variable (i.e. had at least one SNP): 71 positions were parsimony informative, and 39 positions were not parsimony informative. The vast majority of SNPs in this cox1 segment were synonymous; only 10 nucleotide substitutions were nonsynonymous and each was limited to an individual tick and, therefore, parsimony uninformative. More than 60% of the informative synonymous SNPs (n = 44) separated Group 1 from Group 2. The mean genetic distance between the two major D. variabilis clades was 0.055. In contrast, the mean genetic distance was an order of magnitude lower within samples from Group 1 and Group 2 (0.005 and 0.003, respectively). The consensus tree acquired from three equally parsimonious trees is presented in Fig. 2, with CI of 0.577 and RI of 0.966. This tree is a condensed version of the full maximum parsimony tree with identical genotypes collapsed together; it is labeled with information on collection locations (states and/or provinces) and sample size from each location. An additional expanded version of this tree, including identification and other information on each individual, is presented in Additional file 2: Figure S1.

Fig. 1
figure 1

Sampling distribution of 332 D. variabilis ticks sequenced at the mitochondrial cox1 gene. Circles are placed at the geographical center of each state or province and the relative area of each circle is proportional to sample size. Colors indicate membership to one of the two major phylogenetic groups (blue, Group 1; red, Group 2) within the entire D. variabilis population

Fig. 2
figure 2

Phylogenetic relationships among 332 geographically diverse D. variabilis ticks from the USA and Canada. A maximum parsimony consensus of three equally parsimonious mitochondrial cox1 trees is given, illustrating a deep phylogenetic split within the greater D. variabilis population. D. reticulatus serves as the outgroup species

Several states/provinces (especially KS, ME, ON, SD, VA) had a mixture of dissimilar haplotypes that were spread across the entire phylogeny. The number of unique haplotypes per state or province ranged from two in WA to 20 in IN. In addition, the single most common haplotype (found in 30/332 ticks, or 9%) was widely distributed; it occurred in 11 states and provinces from Alabama to New Brunswick and was relatively basal compared to most other haplotypes. However, a few smaller lineages appeared to display more restricted distributions that were only observed farther north (ND, SD, NE, KS, IA, MN, WI, NY, MA) or restricted to the southeastern USA (VA, SC, GA, FL).

Within the two major groups, we observed very few genetic subdivisions that were associated with geography. The correlation between genetic and geographical distance within the eastern ticks (Group 1) was moderate (ρ = 0.179, P < 0.0001), and samples from Group 2 did not reveal a significant correlation (P > 0.05). The most robust relationship between geographical distance and genetic distance was observed when all D. variabilis samples were used (ρ = 0.453, P < 0.0001, which reflects the large genetic and geographical separation between ticks from Group 1 and Group 2.

Mitochondrial 16S rRNA sequencing

A similar D. variabilis phylogeny, including the major division between Group 1 and Group 2, was obtained using a 453 bp sequence of the 16S rRNA gene from 32 geographically diverse D. variabilis individuals (Additional file 1: Table S1). The consensus maximum parsimony tree acquired from five equally parsimonious trees is shown in Additional file 3: Figure S2, with CI of 0.900 and RI of 0.989. However, the 16S RNA gene yielded less phylogenetic resolution than that of the cox1 gene; this was determined by comparing the percentage of variable nucleotide sites, and additionally, the number of haplotypes represented within each gene fragment [16S = 18 SNPs (3.97%) within a 453 bp segment; resulting in 8 haplotypes vs cox1 = 52 SNPs (7.53%) within a 691 bp segment; resulting in 15 haplotypes] using 25 samples that were analyzed for both genes (Additional file 1: Table S1). As a result, no additional individuals were sequenced at this locus.

Francisella-like endosymbiont and pathogen screening

The main genetic discontinuity that distinguishes D. variabilis in the western USA (Group 2) from all other locations (Group 1) is also reflected in the Francisella-like endosymbionts of these ticks (Fig. 3). Of the 1053 D. variabilis individuals screened, 257 ticks produced robust 16S sequences that were highly similar to the Francisella-like endosymbionts known to be associated with several Dermacentor species [20]. It is important to note that, due to the low abundance of endosymbiont DNA in our samples, our assay likely did not detect all positive individuals. The endosymbiont 16S sequences from Group 1 ticks were primarily a single sequence type labelled “T1” (n = 228). We also found five 16S variants in endosymbionts from eastern ticks that differed by 1–2 SNPs; types T3 through T6 were 97–98% identical. All Group 2 ticks from CA and WA shared a distinct endosymbiont 16S type (T2), with just two exceptions (CA, n = 1; WA, n = 1) that had the T1 signature (Fig. 3). One of these ticks (D.v.0411 CA) was sequenced at the cox1 locus and was assigned to D. variabilis Group 2. This could indicate that T1 endosymbionts have been translocated to the western lineage of D. variabilis at some point. A single SNP (corresponding to E. coli 16S position 258) separates the T2 (A nucleotide) from T1 (G nucleotide) endosymbiont lineages and occurs at the end of a hairpin stem in 16S region 2. All endosymbionts shared a 2-bp insertion (immediately downstream of E. coli 16S position 478) in a side loop located in variable region 3; this insertion appears to be unique to the endosymbiont sequences. The other endosymbiont 16S SNPs in groups T3 through T6 were rare and found only in ticks from the northern distribution of D. variabilis (i.e. 5 ticks from SD, MB, NB and SK). The consensus maximum parsimony tree for the 257 Francisella-like endosymbionts is shown in Fig. 3, with CI of 1.000 and RI of 1.000; an expanded tree provides specific information for individual samples (Additional file 4: Figure S3). Interestingly, we detected relatively high levels of genetic differentiation among endosymbionts found in D. variabilis compared to that within other Francisella species (Additional file 5: Figure S4).

Fig. 3
figure 3

Phylogenetic relationships among 257 geographically diverse Francisella-like endosymbiont sequences from the USA and Canada. A maximum parsimony bacterial 16S rRNA tree is given, illustrating a distinct genetic split within the endosymbiont population, which mirrors that of the tick host

Bacterial pathogens were detected among our 1053 D. variabilis samples but in low frequency. Rickettsia spp. were detected in 29 tick samples (CT < 40; Additional file 1: Table S1). These 29 Rickettsia-positive samples produced robust OmpA sequences that were highly similar to Rickettsia species known to be associated with several Dermacentor ticks in North America [41, 42]. Most notably, amplicons consistent with R. montanensis were detected in 20 individual ticks (1.9% prevalence among all individuals screened), amplicons consistent with R. peacockii in 8 ticks (0.8% prevalence), and one amplicon was consistent with a Rickettsia endosymbiont from Ixodes scapularis (WI tick). Most of the positive ticks (27 out of 29) were sampled in central and eastern Canada from 2000 to 2010. A maximum parsimony tree for rickettsial sequences (Additional file 6: Figure S5) provides location information for individual ticks. We note that, due to the low abundance of Rickettsia DNA in the samples, our assay likely did not detect all positive individuals. Using the IS1111 assay, low levels of C. burnetii were detected in two tick samples: D.v.0451 IN (CT 42.1) and D.v.0239 NB (CT 37.5) (Additional file 1: Table S1).

Discussion

Despite the extensive geographical distribution of D. variabilis in North America, and the critical role that these ticks play in the transmission of pathogenic viruses and bacteria to humans and animals, little work has been done to characterize the diversity within this species across its entire range. Past studies have demonstrated the utility of the mitochondrial 16S rRNA locus in detecting genetic variation among D. variabilis populations [5], but few other markers have been used for phylogeographic studies of the American dog tick. In this study, we generated sequence data for two mitochondrial markers (cox1 and 16S) to improve the understanding of genetic structure across the range of this species and to correlate observed patterns with components of its microbiome. Our analysis of both of these markers suggests that the cox1 locus provided the appropriate genetic resolution for addressing our study questions. Here we present the first range-wide patterns of genetic structure and distribution in D. variabilis ticks across North America while also establishing a correlation in phylogenetic patterns between D. variabilis ticks and one of their bacterial symbionts.

We found strong support for the previously noted genetic divergence between eastern (Group 1) and western (Group 2) populations of D. variabilis [5, 22]. The relatively large number of SNPs that separate these two groups (as compared to the limited diversity within these groups) suggests that geographical isolation has played a role in shaping their genetic divergence. There are two prevailing hypotheses that could explain why these two major groups have distinct geographical distributions. First, it is possible that D. variabilis was recently introduced to the west coast by means of a founder event. This could have been the result of a single event or a small number of dispersal events in which a rare mtDNA type from an eastern population became established in the west and persisted until the present, with no subsequent gene flow back to the east. Dermacentor variabilis does use large mammalian hosts with the ability to travel long distances over a short period of time (e.g. native canids or domestic dogs transported by humans) and these hosts could have been the source for the initial establishment of D. variabilis on the west coast, bridging the continental divide that separates these two populations. This scenario would require gravid female ticks to remain attached to the host for the duration of travel, and then detach in an environment suitable for oviposition and survival of the offspring. But a founder event like this would be unlikely because: (i) the duration of adult attachment to the host is probably too short for long-distance movement on a wild host; and (ii) detachment of a gravid female in a suitable habitat is not likely across great portions of the dry intermountain region of the western USA. In a second scenario, it is possible that D. variabilis individuals on the west coast represent a relict population, whose original range was once more widespread and possibly contiguous with the eastern population during a past geologic epoch. Under this model, we suggest that changes in the availability of refugia during previous glacial episodes could explain the current geographical isolation of these groups [43]. Subsequent habitat discontinuity would then serve to enforce their long-term genetic isolation. A separate study focusing on the North American distribution of the Ixodes ricinus complex [27] found 16 SNPs (4.7%) within a 338 bp segment of the mitochondrial 16S rRNA gene. These SNPs separated northern and southeastern clades and were estimated to represent 35,000 years of evolution. In the current study, we found a similar ratio of SNP variation (6.4%) in the cox1 gene that separated the two major D. variabilis lineages. Although the cox1 appears to mutate at a faster rate than 16S in D. variabilis, the similar ratio of mitochondrial SNPs in D. variabilis and I. ricinus complex may indicate a long-time genetic separation of the two D. variabilis groups, lending support for the glacial refugia scenario.

Within the two main phylogenetic groups, the weak genetic differentiation in D. variabilis populations at a regional scale suggests high levels of genetic mixture, which is probably a consequence of host movements. Ticks are highly vulnerable in the environment when off-host, so their survival and dispersal are intrinsically linked to host behavior and movement [44], which can also influence phylogeographic patterns. Dermacentor variabilis ticks are commonly found on companion animals, such as domesticated dogs, and the potential for anthropogenic influence on host/vector dispersal is quite high. In this study, D. variabilis demonstrated little genetic structure among neighboring populations but strong genetic differentiation among distantly separated populations, indicating high levels of regional dispersal, possibly via domesticated dogs (Fig. 2 and Additional file 2: Figure S1); this observation has also been noted in previous studies [7, 44, 45].

The east-west population divergence observed in both endosymbionts and their tick hosts is suggestive of a long-term co-evolutionary relationship in which the endosymbiont populations have diverged in parallel with their vector hosts. Indeed, the amount of genetic differentiation observed within the Francisella-like endosymbionts from this study is much greater than that observed within other Francisella species, such as F. tularensis, F. piscicida, and F. noatunensis (Additional file 5: Figure S4), suggesting that the endosymbionts have been coevolving with D. variabilis for a very long time. Genetically-isolated tick populations can serve as “islands” where genetic drift can lead to differentiation of both the ticks and the pathogens and symbionts that they harbor. For example, local adaptation in D. variabilis populations has been used to explain a notable regional difference in its ability (and inability) to transmit A. marginale [22]. Similarly, the microbiome of D. andersoni differs from region to region, which may have an impact on competence for A. marginale in this species as well [11, 13, 46]. Because endosymbionts may alter the acquisition of pathogenic species, and a geographical correlation appears to exist among ticks, pathogens, and endosymbionts, these findings highlight that a better understanding of vector-microbe relationship may have implications for managing vector-borne diseases.

Conclusions

Ticks are a highly successful and diverse group of parasites, and studies of their natural history, behavior, and genetic/geographical structure are critically important to improve our understanding of their efficiency as vectors, and for implementing effective strategies for controlling ticks and the pathogens they transmit. Phylogenetic analyses allow us to approximate the magnitude, direction, and timing of genetic dispersal, investigate species microevolution, and make inferences about the future distribution of various tick species. So far, studies of dispersal, colonization, and population genetics demonstrate that tick populations continue to harbor considerable diversity in their behavior and genetic structure, and species-specific biological information is needed to accurately deepen our understanding of ticks and the pathogens they transmit. The present study provides the first data concerning an in-depth phylogeographic analysis of D. variabilis and the Francisella-like endosymbionts associated with this tick species across its entire range and provides evidence of long-term tick/symbiont co-evolution. Using cox1 mtDNA sequencing, we found two distinct genetic clades that suggest present-day populations may have been established during a recent glacial period.

Abbreviations

CI:

Consistency index

cox1:

Cytochrome c oxidase subunit 1 gene

NTC:

Non-template control

PCR:

Polymerase chain reaction

qPCR:

Quantitative PCR

RI:

Retention index

SNP:

Single nucleotide polymorphism

References

  1. Yunker CE, Keirans JE, Clifford CM, Easton ER. Dermacentor ticks (Acari, Ixodoidea, Ixodidae) of the New World - a scanning electron microscope atlas. Proc Entomol Soc Wash. 1986;88:609–27.

    Google Scholar 

  2. Sonenshine DE, Roe RM. Biology of ticks. vol. 2. 2nd ed. New York: Oxford University Press; 2014.

    Google Scholar 

  3. de la Fuente J, Estrada-Pena A, Venzal JM, Kocan KM, Sonenshine DE. Overview: ticks as vectors of pathogens that cause disease in humans and animals. Front Biosci-Landmrk. 2008;13:6938–46.

    Article  Google Scholar 

  4. Dodds DG, Martell AM, Yescott RE. Ecology of American dog tick Dermacentor variabilis (Say) in Nova Scotia. Can J Zool. 1969;47:171–81.

    Article  Google Scholar 

  5. Krakowetz CN, Dergousoff SJ, Chilton NB. Genetic variation in the mitochondrial 16S rRNA gene of the American dog tick, Dermacentor variabilis (Acari: Ixodidae). J Vector Ecol. 2010;35:163–73.

    Article  PubMed  Google Scholar 

  6. Stout IJ, Clifford CM, Keirans JE, Portman RW. Dermacentor variabilis (Say) (Acarina: Ixodidae) established in southeastern Washington and northern Idaho. J Med Entomol. 1971;8:143–7.

    Article  PubMed  CAS  Google Scholar 

  7. Araya-Anchetta A, Scoles GA, Giles J, Busch JD, Wagner DM. Hybridization in natural sympatric populations of Dermacentor ticks in northwestern North America. Ecol Evol. 2013;3(3):714–24.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  8. McEnroe WD. Effect of the temperature regime on Dermacentor variabilis (Say) populations in eastern North America. Acarologia. 1979;20:58–67.

    PubMed  CAS  Google Scholar 

  9. Dergousoff SJ, Galloway TD, Lindsay LR, Curry PS, Chilton NB. Range expansion of Dermacentor variabilis and Dermacentor andersoni (Acari: Ixodidae) near their northern distributional limits. J Med Entomol. 2013;50:510–20.

    Article  PubMed  Google Scholar 

  10. Bonnet SI, Binetruy F, Hernandez-Jarguin AM, Duron O. The tick microbiome: why non-pathogenic microorganisms matter in tick biology and pathogen transmission. Front Cell Infect Mi. 2017;7:236.

    Article  Google Scholar 

  11. Gall CA, Scoles GA, Magori K, Mason KL, Brayton KA. Laboratory colonization stabilizes the naturally dynamic microbiome composition of field collected Dermacentor andersoni ticks. Microbiome. 2017;5:133.

    Article  PubMed  PubMed Central  Google Scholar 

  12. Clay K, Fuqua C. The tick microbiome: diversity, distribution and influence of the internal microbial community for a blood-feeding disease vector. In: Workshop Report: critical needs and gaps in understanding prevention, amelioration, and resolution of Lyme and other tick-borne diseases: the short-term and long-term outcomes. Washington D.C: National Academies Press; 2011.

    Google Scholar 

  13. Gall CA, Reif KE, Scoles GA, Mason KL, Mousel M, Noh SM, et al. The bacterial microbiome of Dermacentor andersoni ticks influences pathogen susceptibility. ISME J. 2016;10:1846–55.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  14. Telford SR 3rd. Status of the “East Side hypothesis” (transovarial interference) 25 years later. Ann NY Acad Sci. 2009;1166:144–50.

    Article  PubMed  PubMed Central  Google Scholar 

  15. Burgdorfer WHS, Mavros AJ. Non-pathogenic rickettsiae in Dermacentor andersoni: a limiting factor for the distribution of Rickettsia rickettsii. In: Burgdorfer WAR, editor. Rickettsiae and Rickettsial Diseases. NY: Academic Press; 1981. p. 585–94.

    Google Scholar 

  16. Macaluso KR, Sonenshine DE, Ceraul SM, Azad AF. Rickettsial infection in Dermacentor variabilis (Acari: Ixodidae) inhibits transovarial transmission of a second Rickettsia. J Med Entomol. 2002;39:809–13.

    Article  PubMed  Google Scholar 

  17. Baldridge GD, Burkhardt NY, Simser JA, Kurtti TJ, Munderloh UG. Sequence and expression analysis of the ompA gene of Rickettsia peacockii, an endosymbiont of the Rocky Mountain wood tick, Dermacentor andersoni. Appl Environ Microb. 2004;70:6628–36.

    Article  CAS  Google Scholar 

  18. Niebylski ML, Peacock MG, Fischer ER, Porcella SF, Schwan TG. Characterization of an endosymbiont infecting wood ticks, Dermacentor andersoni, as a member of the genus Francisella. Appl Environ Microb. 1997;63:3933–40.

    CAS  Google Scholar 

  19. Noda H, Munderloh UG, Kurtti TJ. Endosymbionts of ticks and their relationship to Wolbachia spp. and tick-borne pathogens of humans and animals. Appl Environ Microb. 1997;63:3926–32.

    CAS  Google Scholar 

  20. Scoles GA. Phylogenetic analysis of the Francisella-like endosymbionts of Dermacentor ticks. J Med Entomol. 2004;41:277–86.

    Article  PubMed  CAS  Google Scholar 

  21. Araya-Anchetta A, Busch JD, Scoles GA, Wagner DM. Thirty years of tick population genetics: a comprehensive review. Infect Gen Evol. 2015;29:164–79.

    Article  Google Scholar 

  22. de la Fuente J, Van Den Bussche RA, Kocan KM. Molecular phylogeny and biogeography of North American isolates of Anaplasma marginale (Rickettsiaceae: Ehrlichieae). Vet Parasitol. 2001;97:65–76.

    Article  PubMed  CAS  Google Scholar 

  23. Dharmarajan G, Beasley JC, Rhodes OE Jr. Spatial and temporal factors affecting parasite genotypes encountered by hosts: empirical data from American dog ticks (Dermacentor variabilis) parasitising raccoons (Procyon lotor). Int J Parasitol. 2010;40(7):787–95.

    Article  PubMed  CAS  Google Scholar 

  24. Crosbie PR, Boyce WM, Rodwell TC. DNA sequence variation in Dermacentor hunteri and estimated phylogenies of Dermacentor spp. (Acari: Ixodidae) in the New World. J Med Entomol. 1998;35:277–88.

  25. Norris DE, Klompen JS, Keirans JE, WCt B. Population genetics of Ixodes scapularis (Acari: Ixodidae) based on mitochondrial 16S and 12S genes. J Med Entomol. 1996;33:78–89.

  26. Qiu WG, Dykhuizen DE, Acosta MS, Luft BJ. Geographic uniformity of the Lyme disease spirochete (Borrelia burgdorferi) and its shared history with tick vector (Ixodes scapularis) in the Northeastern United States. Genetics. 2002;160:833–49.

    PubMed  PubMed Central  CAS  Google Scholar 

  27. Rich SM, Caporale DA, Telford SR 3rd, Kocher TD, Hartl DL, Spielman A. Distribution of the Ixodes ricinus-like ticks of eastern North America. Proc Natl Acad Sci USA. 1995;92:6284–8.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  28. de la Fuente J, Almazan C, Van Den Bussche RA, Bowman J, Yoshioka JH, Kocan KM. Characterization of genetic diversity in Dermacentor andersoni (Acari: Ixodidae) with body size and weight polymorphism. Exp Parasitol. 2005;109:16–26.

  29. WCt B, Piesman J. Phylogeny of hard- and soft-tick taxa (Acari: Ixodida) based on mitochondrial 16S rDNA sequences. Proc Natl Acad Sci USA. 1994;91:10034–8.

    Article  Google Scholar 

  30. Caporale DA, Rich SM, Spielman A, Telford SR 3rd, Kocher TD. Discriminating between Ixodes ticks by means of mitochondrial DNA sequences. Mol Phylogenet Evol. 1995;4:361–5.

    Article  PubMed  CAS  Google Scholar 

  31. Mitani H, Takahashi M, Masuyama M, Fukunaga M. Ixodes philipi (Acari: Ixodidae): phylogenetic status inferred from mitochondrial cytochrome oxidase subunit I gene sequence comparison. J Parasitol. 2007;93:719–22.

    Article  PubMed  CAS  Google Scholar 

  32. Allender CJ, Easterday WR, Van Ert MN, Wagner DM, Keim P. High-throughput extraction of arthropod vector and pathogen DNA using bead milling. BioTechniques. 2004;37:730. 2:4.

  33. Tamura K, Stecher G, Peterson D, Filipski A, Kumar S. MEGA6: Molecular Evolutionary Genetics Analysis version 6.0. Mol Biol Evol. 2007;30:2725–9.

    Article  CAS  Google Scholar 

  34. Forsman M, Sandstrom G, Sjostedt A. Analysis of 16S ribosomal DNA sequences of Francisella strains and utilization for determination of the phylogeny of the genus and for identification of strains by PCR. Int J Syst Bacteriol. 1994;44:38–46.

    Article  PubMed  CAS  Google Scholar 

  35. O'Neill SL, Giordano R, Colbert AM, Karr TL, Robertson HM. 16S rRNA phylogenetic analysis of the bacterial endosymbionts associated with cytoplasmic incompatibility in insects. Proc Natl Acad Sci USA. 1992;89:2699–702.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  36. Kato CY, Chung IH, Robinson LK, Austin AL, Dasch GA, Massung RF. Assessment of real-time PCR assay for detection of Rickettsia spp. and Rickettsia rickettsii in banked clinical samples. J Clin Microbiol. 2013;51:314–7.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  37. Azad AF, Beard CB. Rickettsial pathogens and their arthropod vectors. Emerg Infect Dis. 1998;4:179–86.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  38. Fergie JE, Purcell K, Wanat D. Murine typhus in South Texas children. Pediatr Infect Dis J. 2000;19:535–8.

    Article  PubMed  CAS  Google Scholar 

  39. Regnery RL, Spruill CL, Plikaytis BD. Genotypic identification of rickettsiae and estimation of intraspecies sequence divergence for portions of two rickettsial genes. J Bacteriol. 1991;173:1576–89.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  40. Loftis AD, Reeves WK, Szumlas DE, Abbassy MM, Helmy IM, Moriarity JR, et al. Surveillance of Egyptian fleas for agents of public health significance: Anaplasma, Bartonella, Coxiella, Ehrlichia, Rickettsia, and Yersinia pestis. Am J Trop Med Hyg. 2006;75:41–8.

    Article  PubMed  Google Scholar 

  41. Niebylski ML, Schrumpf ME, Burgdorfer W, Fischer ER, Gage KL, Schwan TG. Rickettsia peacockii sp. nov., a new species infecting wood ticks, Dermacentor andersoni, in western Montana. Int J Syst Bacteriol. 1997;47:446–52.

    Article  PubMed  CAS  Google Scholar 

  42. Yunik ME, Galloway TD, Lindsay LR. Assessment of prevalence and distribution of spotted fever group rickettsiae in Manitoba, Canada, in the American dog tick, Dermacentor variabilis (Acari: Ixodidae). Vector-Borne Zoonotic Dis. 2015;15:103–8.

    Article  PubMed  Google Scholar 

  43. James AM, Burdett C, McCool MJ, Fox A, Riggs P. The geographic distribution and ecological preferences of the American dog tick, Dermacentor variabilis (Say), in the U.S.A. Med Vet Entomol. 2015;29:178–88.

    Article  PubMed  CAS  Google Scholar 

  44. Dharmarajan G, Fike JA, Beasley JC, Rhodes OE Jr. Development and characterization of 12 polymorphic microsatellite loci in the American dog tick (Dermacentor variabilis). Mol Ecol Resour. 2009;9:131–3.

    Article  PubMed  CAS  Google Scholar 

  45. Lysyk TJ, Scoles GA. Reproductive compatibility of prairie and montane populations of Dermacentor andersoni. J Med Entomol. 2008;45:1064–70.

    Article  PubMed  CAS  Google Scholar 

  46. Scoles GA, Ueti MW, Palmer GH. Variation among geographically separated populations of Dermacentor andersoni (Acari: Ixodidae) in midgut susceptibility to Anaplasma marginale (Rickettsiales: Anaplasmataceae). J Med Entomol. 2005;42:153–62.

    Article  PubMed  Google Scholar 

Download references

Acknowledgements

We thank the many people who donated tick samples for this project, and United States veterinarians and veterinary staff for assisting with sample collection. We are grateful to Antonia Dibernardo, biologist with the Public Health Agency of Canada, for providing tick samples and epidemiological information from across Canada, and to Kathy Mason, for sending colony reared ticks from Washington and Texas.

Funding

Funding for this study was provided via award 2010-65104-20386 from USDA-NIFA. The funder had no role in the design of the study; collection, analysis, and interpretation of the data; and in the writing of the manuscript.

Availability of data and materials

The data supporting the conclusions of this article are included within the article and its additional files. The sequence datasets generated during this study have been made available in NCBI GenBank database under the accession numbers MG833881-MG834531.

Author information

Authors and Affiliations

Authors

Contributions

ELK generated all molecular data, performed data analysis, and drafted the initial manuscript. NES provided substantial scientific guidance and prepared the final manuscript. GAS provided substantial scientific guidance and critically revised the manuscript. CMH provided guidance in the areas of computational biology and phylogenetic analyses, and critically revised the manuscript. JDB supervised all molecular data generation and data analysis, provided substantial scientific guidance, and critically revised the manuscript. DMW conceived the study, provided substantial scientific guidance, and critically revised the manuscript. All authors revised, read and approved the final manuscript.

Corresponding author

Correspondence to David M. Wagner.

Ethics declarations

Authors’ information

ELK is a former Research Technician at the Northern Arizona University (NAU) Pathogen and Microbiome Institute (PMI), where she performed genetics research on arthropod vectors and pathogens. NES is a Research Specialist Senior at PMI where he performs genetics research on arthropod vectors, vertebrate hosts, and pathogens. GAS is a Research Entomologist with the USDA Animal Disease Research Unit at Washington State University, who conducts basic and applied research on tick-borne pathogens of domestic animals. CMH is an Assistant Professor in the Informatics and Computing Program at NAU. JDB is Adjunct Faculty at NAU and an Associate Director of PMI, where he performs research on population genetics of pathogens, arthropod vectors, and vertebrate hosts. DMW is a Professor at NAU; his research program at PMI is focused on pathogens of humans and animals.

Ethics approval and consent to participate

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1:

Table S1. Dermacentor variabilis (n = 1053) were obtained from 31 states and provinces. We included a subset (n = 491), which represents all D. variabilis that produced data for ≥ 1 analysis (cox1 phylogeny, 16S phylogeny, Francisella-like endosymbiont detection and 16S phylogeny, Rickettsia spp. detection and OmpA phylogeny, and Coxiella burnetii detection). A subset were sequenced at cox1 (cox1D.v., n = 332), and/or 16S (16SD.v., n = 32). We also include ticks that were positive and/or sequenced at three species targets: Francisella-like endosymbiont 16S (16SE, n = 257), Rickettsia OmpA ([R], n = 29), and C. burnetii ([C], n = 2). (XLSX 67 kb)

Additional file 2:

Figure S1. Expanded maximum parsimony phylogeny of 332 Dermacentor variabilis from the USA and Canada based on the cox1 gene. A consensus of three equally parsimonious trees is given. Stars indicate samples for which mitochondrial 16S sequences were also generated. (EPS 1788 kb)

Additional file 3:

Figure S2. Expanded maximum parsimony phylogeny of 32 Dermacentor variabilis from the USA based on the 16S rRNA gene. A consensus of five equally parsimonious trees is given illustrating a phylogenetic split within the greater D. variabilis population, which is congruent with the cox1 phylogeny. (EPS 2121 kb)

Additional file 4:

Figure S3. Maximum parsimony phylogeny of 257 Francisella-like endosymbiont 16S rRNA gene sequences illustrates a distinct split within the endosymbionts, which mirrors the tick host. (EPS 4528 kb)

Additional file 5:

Figure S4. Maximum parsimony phylogenetic analysis among six Francisella species and the Francisella-like endosymbionts from this study based on the 16S rRNA gene illustrates a high level of diversity among endosymbionts when compared to closely related Francisella species. (EPS 1279 kb)

Additional file 6:

Figure S5. Maximum parsimony phylogenetic analysis of the OmpA gene in Rickettsia spp. isolated from Dermacentor variabilis. This phylogeny represents the three Rickettsia species and one Rickettsia-like endosymbiont detected during this study. (EPS 1202 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Kaufman, E.L., Stone, N.E., Scoles, G.A. et al. Range-wide genetic analysis of Dermacentor variabilis and its Francisella-like endosymbionts demonstrates phylogeographic concordance between both taxa. Parasites Vectors 11, 306 (2018). https://doi.org/10.1186/s13071-018-2886-5

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s13071-018-2886-5

Keywords