Genetic Melting Pot in Blacklegged Ticks at the Northern Edge of their Expansion Front.

Blacklegged ticks (Ixodes scapularis) are considered to be the main vector of Lyme disease in eastern North America. They may parasitize a wide range of bird and mammal hosts. Northward dispersal of blacklegged ticks has been attributed largely to movement of hosts to areas outside of the current range of the tick, in conjunction with climate change. To better understand the drivers of range expansion in the blacklegged tick, we need investigations of the genetic connectivity and differentiation of tick populations at a fine spatial scale using appropriate markers. In this study, we investigated genetic connectivity and differentiation in blacklegged ticks, in an area of putatively recent advance in Ontario and Quebec, Canada, using microsatellite markers. Our findings suggest patchy differentiation of alleles, no spatial pattern of genetic structure, and genetic subdivision within sites, which are consistent with the very limited evidence available near the leading edge of range expansion of blacklegged ticks into Canada. These findings are consistent with the prevailing hypothesis, drawn from a variety of fields of study, suggesting that migratory birds from a variety of regions may be bringing hitchhiking ticks northward into Canada. © The American Genetic Association 2020.

Blacklegged tick (Ixodes scapularis) populations are increasingly observed in eastern Canada, particularly in southern regions bordering the United States. Their establishment in new areas is attributable to multiple factors, including changes in climate and the environment that have affected both survival of ticks and movement of tick hosts, namely songbirds and mammals (Steere et al. 2004; Ogden, St-Onge, et al. 2008; Simon et al. 2014; Ferrell and Brinkerhoff 2018; Halsey et al. 2018). Long-distance transport of ticks by the means of migratory songbirds is suggested as a major contributing factor for the dispersal of ticks outside of their current range distribution (Ogden, Lindsay, et al. 2008). Transport by mammals has also been suggested as a means of dispersal for blacklegged ticks on a relatively smaller scale (Leighton et al. 2012; Talbot et al. 2019). In eastern North America, blacklegged ticks are considered to be the main vector of Borrelia burgdorferi sensu stricto (s.s.), one of the main agents of Lyme disease, and may efficiently transmit the pathogen to humans and domestic animals (Steere et al. 1983, 2004; Thompson et al. 2001). An association between northward propagation of B. burgdorferi s.s. and northward range expansion of blacklegged ticks has been observed in the same region of Ontario, Canada, from 2010 to 2016 (Clow et al. 2017; Kulkarni et al. 2019). The bordering cities of Ottawa in Ontario and Gatineau in Quebec are situated near the current northern edge of the distribution of blacklegged ticks (Roy-Dufresne et al. 2013; Soucy et al. 2018), making this an important region for understanding factors related to tick and tick-borne pathogen emergence.

Feeding ticks will stay attached to their hosts for a few days to complete a single blood meal (Sonenshine and Roe 2014; Nadolny 2016). Coupled with the fact that blacklegged ticks may parasitize a wide range of bird and mammal hosts (Steere et al. 2004; Shapiro 2014; Halsey et al. 2018), this creates multiple opportunities for short- and long-ranged dispersal of ticks by their hosts. Recently, increasing climatic suitability of areas north of the historic range of I. scapularis in Canada has allowed establishment of reproducing populations of the species (Ogden, St-Onge, et al. 2008; Soucy et al. 2018). Migratory songbirds usually follow major flyways during their spring and fall migrations (Bingman 1980; DeLong et al. 2005; Kirsch et al. 2013; Pocewicz et al. 2013; Vansteelant et al. 2017). Thus, regions surrounding these flyways may be prone to long-distance colonization by hitchhiking blacklegged ticks. Movement by hitchhiking on mammals, such as deer and rodents, is also prevalent in areas bordering forest and road edges (Allan et al. 2003; Trout Fryxell et al. 2015; Delgado et al. 2017; Hornok et al. 2017; Ehrmann et al. 2018), and areas with high canopy cover (Werden et al. 2014; Ferrell and Brinkerhoff 2018). These areas would therefore potentially be subject to short-distance colonization.

Two studies on spatial distribution of blacklegged ticks (Talbot et al. 2019) and B. burgdorferi strains (Mechai et al. 2018), in regions of recent advance in Canada, found an effect of landscape elements that are more commonly associated with mammal movement, which tends to be local in scale, such as proximity to roads and forest connectivity. Currently, although migratory birds are known to carry B. burgdorferi-infected blacklegged ticks into their breeding ranges in Canada (Ogden, Lindsay, et al. 2008, 2015), sparse evidence demonstrated an effect of avian movement, which tends to be long-distance in scale, on spatial distribution of blacklegged ticks and B. burgdorferi in Canada. However, a predictive model-based study suggested a significant effect of both types of movements on tick invasion and establishment across Canada (Leighton et al. 2012). This could have to do with how avian movement operates, which tends to not largely be affected by landscape elements or geographic distance, and can therefore be difficult to track with conventional methods.

A considerable body of work has investigated various population genetic and phylogeographic aspects of the blacklegged tick across North America (Araya-Anchetta et al. 2015). An interesting outstanding pattern from these studies is the consistently high genetic structure, at the scale of the species range in eastern North America, despite the use of migratory songbirds for dispersal (Ogden, Lindsay, et al. 2008). Local gene flow limitations have been reported in several studies, particularly near the northern range limit of the species (Humphrey et al. 2010; Krakowetz et al. 2011; Araya-Anchetta et al. 2015). These local-scale studies have focused on nuclear and mitochondrial markers with low mutation rates. A more recent study, also near the northern edge of the distribution of blacklegged ticks in Canada, focused on microsatellite DNA (Leo et al. 2016). A similar pattern of local genetic structure was observed, whereby patchy differentiation of alleles was present around the study area, without any clear spatial pattern in distribution of these genetic differences (Leo et al. 2016). The authors noted high inbreeding coefficients, which may be caused by genetic subdivision within sites. The authors concluded that the observed patterns were consistent with a scenario of repeated long-distance movement in the species, for example, through bird migration, as opposed to a scenario of a slowly advancing colonization front, for example, through mammal dispersal (Leo et al. 2016). Here, we investigated the consistency of these patterns across the northern leading edge of the distribution of blacklegged ticks in Canada using microsatellite DNA.

In this study, we investigated genetic connectivity and differentiation in blacklegged ticks, in an area of putatively recent advance in Ontario and Quebec, Canada, using microsatellite markers. Based on previous studies in a similar setting, we predicted patchy genetic differentiation across sites, weak or nonexistent spatial pattern of genetic structure, and evidence of genetic substructuring within sites. Our study is one of a limited number of studies performed at the northern edge of the distribution of blacklegged ticks that investigated patterns of genetic connectivity and differentiation using microsatellite markers. These studies are useful in our understanding of patterns in which these ectoparasites are propagating the agent of Lyme disease northward into Canada.

Materials and Methods

Data Collection

We used a set of adult blacklegged ticks (Ixodes scapularis) previously collected from 2017 to 2018 in 2 rounds of sampling: one in late spring/early summer and one in the fall. Details of the collection methods can be found in Kulkarni et al. (2018). Sampling sites are situated near Gatineau, Quebec, and Ottawa, Ontario, in Canada. We used only adult blacklegged ticks to prevent sampling bias arising because of spatial aggregation of individuals from the same family group (i.e., associated either by a parent–offspring or a full-sibling relationship), which has been observed to influence results of population genetic analyses (Goldberg and Waits 2010). We included ticks from sites where at least 3 adult ticks were available for use in this study, to minimize sampling bias of nonresident ticks, and we selected a maximum of 10 ticks in any one site to limit sampling bias arising from large differences in sample size among sites.

Genetic Analyses

We extracted DNA from whole male ticks and from one half of female ticks, cut longitudinally, using the DNeasy Blood & Tissue Kit (QIAGEN, Germantown, MD), according to the manufacturer’s instructions. We used a 2720 Thermal Cycler (Applied Biosystems, Foster City, CA; ABI) to execute the polymerase chain reaction (PCR) amplification of 10 microsatellite loci developed specifically for I. scapularis by Fagerberg et al. (2001), using their original protocol. We visualized PCR products by 1.5% agarose gel electrophoresis using SYBR Safe (Invitrogen, Waltham, MA) on a Chemidoc MP Imaging System (Bio-Rad, Hercules, CA) to check the quality and size of amplified fragments. We then sized products on a 3730xl DNA Analyzer (ABI). We called all microsatellite genotypes for each species using Geneious Software v.2019.2.1 (Biomatters, Ltd., Auckland, New Zealand), and we checked all calls manually.

Statistical Analyses

We looked for evidence of genetic differentiation in the study area using 2 genetic clustering approaches. For both analyses, we used site of collection as prior population information. First, we used a Bayesian clustering algorithm implemented using Structure v2.3.4 software (Pritchard et al. 2000). We computed probability of data to estimate the number of genetic clusters (K), using 100 000 iterations, and a burn-in period of 10 000 iterations. We executed 10 runs for each K value between 1 and 10. We used the admixture model, individual alpha for each population, and allele frequencies correlated among populations and inferred a separate lambda for each population. We kept all other parameters at default value. We used the Evanno et al.’s (2005) method to determine the most probable K. We executed the analysis a second time, through 10 runs, using the most probable K determined in the previous analysis, to calculate membership probability of each individual to each identified genetic cluster. We used the Kopelman et al.’s (2015) method to average membership probability across major modes. Second, we conducted a discriminant analysis of principal components clustering analysis using the “find.clusters” function of the “adegenet” package (Jombart 2008) in R v3.6.0 (R Development Core Team, Vienna, Austria). We used the “kmeans” method and used 95% as the minimal percentage of the total variance of the data to be expressed by the retained axes of the principal component analysis. We kept all other parameters at default value. We used Bayesian information criterion values to determine the most likely K value. For both analyses, we assigned each individual to its highest probability genetic cluster for mapping purposes. We expected K values higher than 1 for both analyses, with no clear association with sites of collection, which would be consistent with patterns of patchy genetic differentiation.

We looked for evidence of spatial patterns of genetic structure across the study area, using 2 matrix-based linear multiple regression approaches, one at the site level and one at the individual level, and 2 global indices. In both regression approaches, we measured the effect of geographic distance on genetic distance. At the site level, we used the index of differentiation Hedrick’s G’ ST (Hedrick 2005). We calculated G’ ST for all pairs of sites, using the “fastDivPart” function of the “diveRsity” package (Keenan et al. 2013) in R v3.6.0. We calculated geographic distance (in kilometer) between sampling sites, corrected for sphericity of the earth, using the “rdist.earth” function from the “fields” package (Fields Development Team 2006) in R v3.6.0. We then fit pairwise genetic distance to geographic distance using a multiple regression on distance matrices, using the “MRM” function from the “ecodist” package (Goslee and Urban 2007) in R v3.6.0, which uses a Mantel test derived linear regression model. We assessed significance through a permutation procedure (9999 replicates) that takes into account nonindependence of data points in distance matrices (Legendre et al. 1994; Lichstein 2007). At the individual level, in addition to testing the effect of geographic distance, we looked for evidence of temporal and sex-biased effects, which could obscure subtle patterns of spatial genetic structure. Therefore, we tested the effect of geographic distance, year of collection (2017 and 2018), season of collection (spring and fall), and sex, on genetic distance, in a multiple regression on distance matrices model. We calculated Smouse and Peakall’s (1999) individual genetic distance for each pair of individuals, using GenoDive v2.0. For each pair of individuals collected in the same year, we assigned a value of 0, and for those collected in different years, we assigned a value of 1. We proceeded the same way for season of collection and sex. For visualization purposes, we conducted a standard principal coordinate analysis of Smouse and Peakall’s individual genetic distances, using GenAlEx v6.51 (Peakall and Smouse 2012). Last, we computed Hedrick’s global G’ ST among sites (Hedrick 2005), and we conducted an analysis of molecular variance, using GenoDive v2.0 (Meirmans 2012). We expected no strong effect of geographic distances on genetic distances, which would be consistent with weak or nonexistent patterns of spatial genetic structure.

Finally, we looked for evidence of genetic substructuring within study sites, otherwise known as the Wahlund effect (Garnier-Géré and Chikhi 2013). The Wahlund effect may be caused by sampling of multiple nonrandom mating populations within each study site. Sites affected by a Wahlund effect would show significant Hardy–Weinberg and/or linkage disequilibrium similarly across all markers, and for each marker, we would observe no clear association between differentiation and inbreeding coefficients. We tested all microsatellite markers for Hardy–Weinberg and linkage disequilibrium, using FSTAT 4.2 (Raymond and Rousset 1994). For each type of test, we corrected for multiple tests using Bonferroni correction, with a threshold α of 0.05. In case of widespread Hardy–Weinberg and/or linkage disequilibrium across loci, we used the determination key proposed by De Meeûs (2018) to determine what demographic or technical processes are most likely driving our data. In these regards, we calculated number of alleles, observed and expected heterozygosities, inbreeding coefficient GIS, and global Hedrick’s G’ ST, using GenoDive v2.0 (Meirmans 2012), for each microsatellite marker separately and overall. We used Pearson’s correlation coefficient, using the “cor.test” function in R v3.6.0, to test correlation between global Hedrick’s G’ ST and GIS coefficients and between proportion of missing data and GIS coefficients, across microsatellite markers. We expected that genetic substructuring within sites, or Wahlund effect, would be the most likely process affecting our data.

Results

We successfully amplified a total 107 individual blacklegged ticks (Ixodes scapularis); of which, 46 were males and 61 were females. The majority of ticks (103; 96%) were collected from 14 sites in areas in and around the city of Ottawa (site identifiers starting with O; Figure 1), while 4 ticks (4%) were collected from 1 site near Gatineau (Site G2; Figure 1). Number of genotyped adult ticks per site varied between 2 and 10, reflecting local densities of blacklegged ticks. One marker (CAG12) was dropped from the analyses due to its monomorphism across all samples. The resulting 9 markers were used for all analyses. Microsatellite genotypic data can be found in Supplementary Table S1.

Figure 1.

Map of study sites where Ixodes scapularis ticks were collected around Ottawa and Gatineau in Canada in 2017–2018, showing genetic cluster membership calculated through (A) Structure (Pritchard et al. 2000) and (B) Discriminant analysis of principal components (Jombart 2008) approaches. Pie charts represent number of samples (size; from 2 to 10) and proportion of genetic cluster membership at each study site (proportion of individuals with highest probability of membership to a genetic cluster; shadings as in legend). Lines connecting sites represent a pairwise Hedrick’s G’ ST differentiation index value of 0.05 or less. Built-up land is shown in dark gray shading and water bodies in lighter shading. Map was created using ArcGIS 10.5 (ESRI, Redlands, CA).

Genetic Differentiation

Our genetic clustering analysis using the Structure approach showed highest probability of a total of 6 clusters in our dataset, with a less likely possibility of 3 clusters (Supplementary Figures S1 and S2). Interestingly, 13 sites of 15 contain individuals that were assigned to a total of at least 3 genetic clusters (Figure 1). Additionally, individuals assigned to a total of 4 of 6 genetic clusters are distributed in at least 7 sites spanning most of the study area (Figure 1). These results suggest no spatial population genetic structure, but rather co-occurrence of genetic clusters within sites around the study area. Our genetic clustering analysis using the discriminant analysis of principal components approach identified a total of 4 clusters in our dataset (Supplementary Figure S3), with no evidence of spatial population genetic structure. There was no concordance in genetic cluster assignment of individuals between the 2 approaches, and both approaches show spatially random patterns of genetic differentiation across the study area (Figure 1). Altogether these results suggest that genetic ancestry of lineages present in the study area may be resolved in different ways depending on which statistical method is used, and both approaches suggest patchy differentiation of alleles across the study area.

Spatial Genetic Structure

Our multiple regression on distance matrices analyses revealed no significant effect of geographic distance on genetic distance, both at the site and individual levels (Table 1; Figure 2). The year of collection showed a weak positive effect on genetic distance at the individual level (Table 1; Figure 3), which suggests possible yearly differences in allelic frequencies in the study area. Global Hedrick’s G’ ST calculated using all markers shows a null value (Table 2). A total of 23 pairs of sites show a Hedrick’s G’ ST value that is lower than 0.05. A total of 17 of these pairs involve either site O13, O20, and/or O23, and connect them with a total of 9 other sites (Figure 1). Analyses of molecular variance show 100% of variation within sites. Altogether, these results suggest no spatial patterns of genetic structure in the study area.

Table 1.

Slope and R2 values for 2 analyses of multiple regression on distance matrices on the effect of geographic distances (km) on genetic distances, for Ixodes scapularis ticks collected in 15 sites around the cities of Ottawa and Gatineau in Canada in 2017–2018

Analysis level	Geographic distance	Season	Year	Sex	R  2
Site	<0.001	—	—	—	0.020
Individual	<0.001	−0.006	0.059	−0.004	0.004

For the site-level analysis, we used pairwise Hedrick’s G’ ST values as genetic distances. For the individual-level analysis, we used pairwise Smouse and Peakall’s individual genetic distances, and conditioned models for the year of collection (0: same year, 1: different year), season of collection (0: same season, 1: different season), and sex (0: same sex, 1: different sex). Bold values are significant at α = 0.05.

Figure 2.

Graph showing the effect of geographic distance (km) on genetic distance (Hedrick’s G’ ST divided by the inverse of Hedrick’s G’ ST), for Ixodes scapularis ticks collected in 15 sites around Ottawa and Gatineau in Canada in 2017–2018.

Figure 3.

Map of the first 2 principal coordinates (with a cumulative explained variance of 10.5%) of Smouse and Peakall’s individual genetic distance for Ixodes scapularis ticks collected in 15 sites around Ottawa and Gatineau in Canada in 2017–2018. Shading of data points is based on the year of sampling.

Table 2.

Population genetic parameters for each analyzed microsatellite locus and averaged across loci, for Ixodes scapularis ticks collected in 15 sites around the cities of Ottawa and Gatineau in Canada in 2017–2018

Locus	N  A	H  O	H  S	G  IS	G’ ST	HWE P
AC4	36	0.54	0.90	0.40	−0.11	<0.01
AG4	10	0.15	0.27	0.45	0.01	<0.01
CTGY17	6	0.09	0.13	0.32	−0.02	<0.01
AC20	4	0.10	0.09	−0.05	0.02	0.21
AC22	5	0.08	0.08	0.02	−0.04	0.73
GATA3	10	0.15	0.47	0.68	−0.03	<0.01
AG25	13	0.56	0.83	0.32	0.16	<0.01
AC8	38	0.81	0.94	0.13	−0.12	<0.01
GATA4	47	0.36	0.98	0.63	0.13	<0.01
All	19	0.32	0.52	0.39	0.00	<0.01

For each locus, we show number of alleles (NA); averaged observed heterozygosity (HO), expected heterozygosity (HS), and inbreeding coefficient (GIS) across study sites; global Hedrick’s differentiation index (G’ ST); and P of Hardy–Weinberg equilibrium (HWE) overall tests. We also show average NA across all loci, and values for all other parameters when calculated using all loci.

Genetic Substructuring

We found no significant case of linkage disequilibrium across all pairs of microsatellite markers, at a threshold of 0.05 after Bonferroni correction. We found significant deviations of Hardy–Weinberg equilibrium at 7 of 9 markers (excluding only AC20 and AC22; Table 2), at a threshold of 0.05 after Bonferroni correction, and all were characterized by homozygote excess. Average number of alleles was 19 (Table 2), and markers CTGY17, AC20, and AC22 showed fewer than 10 alleles. Average expected heterozygosities were higher than average observed heterozygosities, and inbreeding coefficient GIS values were moderately high, for all markers except AC20 and AC22 (Table 2). Pearson’s correlation coefficients between global Hedrick’s G’ ST and GIS coefficients and between proportion of missing data and GIS coefficients, across microsatellite markers, were nonsignificant at a threshold of 0.05 (P = 0.46 and 0.13, respectively). We found small variations in global Hedrick’s G’ ST (most values were between 0.00 and 0.02, and all values were between 0.00 and 0.16), and moderate variations in GIS (most values were between 0.32 and 0.68, and all values were between 0.00 and 0.68; Table 2). Standard error of global Hedrick’s G’ ST calculated using all markers was of 0.01, and standard error of GIS was 0.09.

Discussion

Our results support earlier findings from studies conducted at the northern edge of the spatial distribution of blacklegged ticks in eastern Canada. Indeed, we observed genetic clustering characterized by patchy genetic differentiation that is not distributed in a clear spatial pattern. We did not observe any significant effect of geographic distance on genetic distance in the study area. Our data seem to be affected by a Wahlund effect or genetic subdivisions within study sites. Genetic subdivision may be caused by a variety of biological mechanisms, including behavioral reproductive isolation, geographic barriers to dispersal, and repeated long-distance colonization (Selander 1970; Bohonak 1999; Kempf et al. 2010; Dharmarajan et al. 2011; Garnier-Géré and Chikhi 2013). These findings are similar to those from an earlier study, which analyzed a similar number of sites distributed in an area of similar size, over 4 years, using the same molecular markers (Leo et al. 2016). The authors concluded that the observed patterns were consistent with a scenario of repeated long-distance movement in the species, as opposed to a scenario of a slowly advancing colonization front (Leo et al. 2016). Interestingly, these patterns seem to point to a larger contribution of avian movement compared to terrestrial movement on blacklegged tick genetic connectivity. Our findings and the findings from Leo et al. (2016) may also support earlier findings performed on nuclear and mitochondrial loci, which identified some local gene flow limitations in populations in Canada and northern United States (Humphrey et al. 2010; Krakowetz et al. 2011; Araya-Anchetta et al. 2015), near the leading edge of the expansion front of the species.

The Wahlund effect can be triggered by a number of biological mechanisms that prevent groups of individuals from a putative population from interbreeding with each other (Garnier-Géré and Chikhi 2013). For example, behavioral reproductive isolation (Bohonak 1999) or geographic barriers to gene flow (Selander 1970) would prevent groups of individuals from mating with each other, and genetic drift within each group would lead to genetic substructuring. In blacklegged ticks, reproduction mostly occurs on white-tailed deer (Werden et al. 2014; Halsey et al. 2018), and a number of alternative hosts, such as rodents and birds, may move juvenile ticks across the landscape (Werden et al. 2014; Halsey et al. 2018). Therefore, any potential behavioral or geographic barrier to gene flow would have a limited effect on blacklegged tick populations. An alternative cause of the Wahlund effect is the inadvertent combining of populations with different allele frequencies (Dharmarajan et al. 2011). This appears to be frequent in tick populations, which usually partition according to life-history stages, and further partition according to host species, age, and sex (Dharmarajan et al. 2011). In our study, gene flow patterns of sampled blacklegged ticks are correlated with movement of all the different hosts with which they have been associated at any point in their life cycle. Therefore, there is a high likelihood that our sampling sites are populated by ticks that have been either exclusively transported locally by mammals, exclusively transported over a long a distance by birds, or both processes sequentially. This alternative explanation seems to be the most likely, given that combining of data was done solely based on spatial distribution of study sites, and patchy genetic differentiation and lack of spatial genetic structure were observed in our study area.

One hypothesis, proposed by Leo et al. (2016), postulates that migratory birds, originating from a large area across the range of blacklegged ticks, may be responsible for propagating hitchhiking ticks into localized areas in Canada and northern United States, creating melting pots of genetic diversity at the leading edge of their expansion front. Other sources also suggest such a mechanism putatively causing long-distance colonization of blacklegged ticks across North America (Cohen et al. 2015; Khatchikian et al. 2015) and at the leading edge of their distribution in Canada (Morshed et al. 2005; Ogden, Lindsay, et al. 2008). Our results seem to be consistent with this prevailing hypothesis. Two studies in Canada have found a significant effect of proximity to landscape elements that are usually associated with mammal movement, on spatial distribution of blacklegged tick populations and Borrelia burgdorferi strains (Mechai et al. 2018; Talbot et al. 2019). These 2 studies did not suggest any significant effect of landscape elements that are usually associated with avian movement. However, numerous studies have suggested a role of migratory birds in bringing hitchhiking blacklegged ticks in Canada from more southerly areas (Ogden, Lindsay, et al. 2008; 2015; Leighton et al. 2012). We conclude that investigating spatial genetic connectivity and differentiation in blacklegged ticks could lead to insights that may be harder to detect using purely spatial analyses, particularly due to the strong association of these ticks with migratory songbirds.

We found a weak but significant effect of year of sampling on individual genetic distances, which means that different cohorts of ticks may have been sampled between years at the same sites. Combining data from different years may therefore have affected the results of our analyses to a small extent. Our study therefore suffers from small sample sizes within sites and within years of sampling, which is to be expected near the edge of the spatial distribution of most species. Another limitation of our study is the resolution of molecular markers used, which could limit our ability to detect subtle demographic processes, such as tick movement using mammalian hosts. Although microsatellite markers hold more genetic information than almost any other molecular marker (Bowcock et al. 1994), a number of organisms only possess a limited number of these markers across their genomes, including I. scapularis (Fagerberg et al. 2001). On the other hand, single-nucleotide polymorphisms are usually detectable throughout the genome of most organisms with current methods, and large datasets can thus be obtained (Davey et al. 2010; Schmitt et al. 2012; Zimin et al. 2014). Demographic and ancestry predictions using large molecular datasets are more robust and could deepen our understanding of the factors influencing movement of blacklegged ticks in Canada. We suggest that analyzing the subtle effects of landscape structure on gene flow could inform on rate and direction of local-scale transport by mammals, which would best be performed using a large molecular dataset.

Our study adds to the body of work investigating drivers of colonization of blacklegged ticks, the main vector of B. burgdorferi, a bacterial agent of Lyme disease, in areas near the northern edge of the species range in Canada. Our findings further reinforce the prevailing hypothesis suggesting that migratory birds from a variety of regions may be dispersing hitchhiking ticks, which would explain local patchy differentiation that does not follow a clear spatial pattern. Future investigations on the population genetic structure of blacklegged ticks should focus on the detection of subtle gene flow processes to better understand Lyme disease propagation at the local scale. Understanding such subtle processes in detail is crucial to predict the rate and direction of advance of B. burgdorferi and other tick-borne pathogens, in a context of rapid climate and environmental change in Canada and elsewhere.

Supplementary Material

Supplementary material is available at Journal of Heredity online.

Table S1. Site identifier, as in Figure 1, geographic coordinates, and genotypic data at 9 microsatellite markers, for 107 Ixodes scapularis individuals. The first three digits in each string describe the length (in base pairs) of the first allele, and the last three describe the length of the second allele. Missing data is identified with “000000.”

Figure S1. Bar chart showing genetic cluster membership. Probability of membership was calculated through a Structure approach (top: K = 6; bottom: K = 3), for all 107 adult Ixodes scapularis ticks, sorted by sampling site. Site identifiers are the same as in Figure 1.

Figure S2. Graph showing probability (Delta K) for each value of number of genetic clusters (K). Delta K was calculated by the Evanno et al. (Evanno et al. 2005) method following a Structure approach (Pritchard et al. 2000).

Figure S3. Graph showing Bayesian information criterion (BIC) for each value of number of genetic clusters (K). BIC values were calculated through a discriminant analysis of principal components approach (Jombart 2008).

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.