MEEGID 1970
No. of Pages 8, Model 5G
4 June 2014 Infection, Genetics and Evolution xxx (2014) xxx–xxx 1
Contents lists available at ScienceDirect
Infection, Genetics and Evolution journal homepage: www.elsevier.com/locate/meegid 5 6
Population genetic structure of the Lyme disease vector Ixodes scapularis at an apparent spatial expansion front
3 4 7
Q1
8 9 10
Q2
Rebecca R. Kelly a, David Gaines b, Will F. Gilliam a, R. Jory Brinkerhoff a,c,⇑ a
11 12
Department of Biology, University of Richmond, 28 Westhampton Way, Richmond, VA 23173, United States Division of Environmental Epidemiology, Virginia Department of Health, Richmond, VA 23219, United States c School of Life Sciences, University of KwaZulu-Natal, Pietermaritzburg, South Africa b
a r t i c l e
1 2 4 8 15 16 17 18 19
i n f o
Article history: Received 16 January 2014 Received in revised form 18 May 2014 Accepted 21 May 2014 Available online xxxx
20 21 22 23 24 25 26 27
Keywords: Disease ecology Range expansion Acari Zoonosis Tick-borne disease Epidemiology
a b s t r a c t Modeling and empirical evidence suggests that Lyme disease is undergoing geographic expansion from principal foci in the midwestern and northeastern United States. Virginia is at the southern edge of the current expansion zone and has seen dramatic rise in human Lyme disease cases since 2007, potentially owing to a recent increase in vector abundance. Ixodes scapularis is known throughout the eastern US but behavioral or physiological variation between northern and southern lineages might lead northern-variant ticks to more frequently parasitize humans. We hypothesized that recent spatial and numerical increase in Lyme disease cases is associated with demographic and/or spatial expansion of I. scapularis and that signals of these phenomena would be detectable and discernable in population genetic signals. In summer and fall 2011, we collected nymphal I. scapularis by drag sampling and adult I. scapularis from deer carcasses at hunting check stations at nine sites arranged along an east–west transect through central Virginia. We analyzed 16S mtDNA sequences data from up to 24 I. scapularis individuals collected from each site and detected a total of 24 haplotypes containing 29 segregating sites. We found no evidence for population genetic structure among these sites but we did find strong signals of both demographic and spatial expansion throughout our study system. We found two haplotypes (one individual each) representing a lineage of ticks that is only found in the southeastern United States, with the remaining individuals representing a less genetically diverse clade that is typical of the northern United States, but that has also been detected in the American South. Taken together, these results lead us to conclude that I. scapularis populations in Virginia are expanding and that this expansion may account for recent observed increases in Lyme disease. Ó 2014 Published by Elsevier B.V.
29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
51 52
1. Introduction Vector-borne zoonoses are inherently complex disease systems that are keenly sensitive to anthropogenic and natural environmental change (Mills et al., 2010; Thompson, 2013). The geographic distributions of tick-borne pathogens, as well as phenology of host seeking ticks, have been linked to climate change (reviewed in Gage et al., 2008) and there have been numerous reports of tick-borne diseases increasing in incidence and spatial extent (Danielova et al., 2006; Materna et al., 2008; Medlock et al., 2013). Although large-scale surveillance efforts to generate risk maps for tick-borne pathogens have been undertaken, such efforts are highly resource-intensive and the static models generated by these efforts provide only a ‘snapshot’ image of dynamic
53 54 55 56 57 58 59 60 61 62 63 64
⇑ Corresponding author at: Department of Biology, University of Richmond, 28
Q3
Westhampton Way, Richmond, VA 23173, United States. Tel.: +1 18044841592. E-mail address: [email protected] (R. Jory Brinkerhoff).
processes (Diuk-Wasser et al., 2010, 2012). Molecular approaches, however, can lend insights into evolutionary dynamics and can be used to infer population histories and pathogen transmission dynamics (Guhl and Ramirez, 2013). The spread of the Lyme disease agent, Borrelia burgdorferi, and its principal vector in eastern North America (Ixodes scapularis), into Canada has justifiably been the focus of a great deal of attention (e.g. Koffi et al., 2012; Leighton et al., 2012; Ogden et al., 2009, 2013; Wu et al., 2013). The population genetic structure of I. scapularis in Canada suggests recent demographic expansion (Krakowetz et al., 2011; Mechai et al., 2013) and the genotypic heterogeneity among Canadian B. burgdorferi samples revealed the same suite of genotypes in both southeastern Canada and the eastern United States, albeit with Canadian populations showing signals of founder events (Ogden et al., 2011). There has been less attention paid, however, to other potential Lyme disease expansion fronts. In a large scale study of I. scapularis occurrence throughout the eastern United States, Diuk-Wasser et al. (2010, 2012)
http://dx.doi.org/10.1016/j.meegid.2014.05.022 1567-1348/Ó 2014 Published by Elsevier B.V.
Please cite this article in press as: Kelly, R.R., et al. Population genetic structure of the Lyme disease vector Ixodes scapularis at an apparent spatial expansion front. Infect. Genet. Evol. (2014), http://dx.doi.org/10.1016/j.meegid.2014.05.022
65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82
MEEGID 1970
No. of Pages 8, Model 5G
4 June 2014 2 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125
R.R. Kelly et al. / Infection, Genetics and Evolution xxx (2014) xxx–xxx
developed environmental models to estimate tick population density and Lyme disease risk and found evidence consistent with overall spatial expansion of tick populations. Hamer et al., 2010 demonstrated apparent northward movement of I. scapularis along the eastern shore of Lake Michigan over a period of four years (Hamer et al., 2010). Lee et al. (2013) compared tick occurrence on hunter-shot deer from 1981 to 2009 and concluded that I. scapularis expanding from west to east throughout Wisconsin. Brinkerhoff et al. (submitted for publication) showed data consistent with a shifting geographic distribution of Lyme disease cases in Virginia and suggested, based on field sampling and vector phylogenetic analysis, that tick populations in the western part of the state may be larger than they had been in previous decades. The presence of I. scapularis and incidence of Lyme disease in the southern United States is highly contentious; the occurrence of Lymelike disease in areas where I. scapularis is uncommonly encountered and the importance of other putative disease agents such as Borrelia lonestari contribute to the confusion about tick-borne disease in southern and southeastern states (Stromdahl and Hickling, 2012; Lantos et al., 2013). As a result, it is critical to explore the potential southward expansion of Lyme disease into states from which it has historically been rare or absent. The evolutionary history of I. scapularis, inferred from variation in mitochondrial rDNA and other markers, suggests that this species was relegated to southern North America during Pleistocene glaciation events, and that northern North America was recolonized by founding I. scapularis populations that moved north and diversified following glacial retreat (Norris et al., 1996; Qiu et al., 2002; Humphrey et al., 2010; van Zee et al., 2013). The purported result of this movement is higher genetic diversity of I. scapularis in the southern United States, including a more genetically diverse lineage (hereafter referred to as ‘‘southern clade’’) that is only found in the southern US, and a second lineage (hereafter referred to as ‘‘American clade’’), that is found in both northern and southern states (Norris et al., 1996; Mixson et al., 2004). Potential differences between these two lineages that may affect host seeking or other behaviors (e.g. Durden et al., 2002) could account for differences in infection prevalence and risk to Lyme disease between northern and southern states (Diuk-Wasser et al., 2012; Stromdahl and Hickling, 2012). Since 2007, increasing Lyme disease cases have been reported in Virginia and Brinkerhoff et al. (submitted for publication) demonstrated that Virginia counties with recent increases in human
Lyme disease cases are associated with relatively high I. scapularis densities. Our goal was to use population genetic analysis of Ixodes scapularis sampled throughout central Virginia to identify signals of recent spatial or demographic population expansion. We expected to find relatively lower genetic diversity in ticks associated with recent increases in human LD (i.e. western VA) and higher diversity in locations that have historically been associated with I. scapularis populations (eastern VA). Similarly, we expected to see population genetic structure at the state-wide level and evidence of recent spatial and/or demographic expansion at western sites.
126
2. Materials and methods
137
2.1. Tick collection
138
Tick collection sites were selected along an east–west transect from the Atlantic coast to high-elevation counties in west-central Virginia (Fig. 1). These sites represent variation in recent LD emergence in humans with the westernmost sites having seen dramatic increases in incidence since the surveillance case definition for Virginia was modified to require two-tier testing (i.e. enzymelinked immunosorbent assay and immunoglobulin M-positive Western blot) in 2008 (Fig. 1, Brinkerhoff et al. (submitted for publication)). Host-seeking I. scapularis nymphs were collected by drag sampling from three sites summer 2011 (Fig. 1; Brinkerhoff et al., submitted for publication). At each of these sites, we haphazardly selected and established five 100 m transects through deciduous or mixed evergreen-deciduous forest and collected ticks by dragging a 1 1 m corduroy cloth along both sides of each transect for a total of 1000 square meters of sampling area (Falco and Fish, 1992). For each transect sampling event, we checked the drag cloth at 20-meter intervals and placed all observed ticks in uniquelynumbered vials containing a 70% ethanol solution. Each site was visited four times between May and July 2011 (Brinkerhoff et al. (submitted for publication)). Adult I. scapularis were collected from hunter-shot deer at six additional sites in October and November 2011 (Fig. 1). During these collection events, we carefully removed as many adult I. scapularis as possible while hunters reported to check-in stations. Host-seeking ticks were stored in 70% ethanol and adult ticks were stored in 70% ethanol at 4 °C until identification by light microscopy using dichotomous keys (Sonenshine,
139
Fig. 1. Locations of field sites sampled in this study. Triangles represent sites from which adult I. scapularis were collected at deer hunts and diamonds represent sites where I. scapularis nymphs were collected by drag sampling. Summary statistics and sample sizes associated with each site are indicated in Table 1. County shading represents the mean annual change in Lyme disease incidence between 2008, when the surveillance case definition was made more rigorous in Virginia, and 2011 (incidence data from Brinkerhoff et al. (submitted for publication)).
Please cite this article in press as: Kelly, R.R., et al. Population genetic structure of the Lyme disease vector Ixodes scapularis at an apparent spatial expansion front. Infect. Genet. Evol. (2014), http://dx.doi.org/10.1016/j.meegid.2014.05.022
127 128 129 130 131 132 133 134 135 136
140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164
MEEGID 1970
No. of Pages 8, Model 5G
4 June 2014 3
R.R. Kelly et al. / Infection, Genetics and Evolution xxx (2014) xxx–xxx
170
1979). Based on availability, adult ticks including males and the least engorged females were selected from each hunt for B. burgdorferi testing (when possible, approximately equal numbers of male and female ticks were selected from each site). We subsequently used a random number generator to select a subset of up to 25 ticks per site for genotypic analysis.
171
2.2. Molecular methods
172
194
Following flash-freezing in liquid nitrogen, we pulverized ticks in individual microcentrifuge tubes with a sterile plastic pestle. We extracted DNA from individual I. scapularis using Qiagen DNeasy Blood and Tissue Kit following the manufacturer’s protocols except we performed two 50 ul DNA elutions per sample for a total of 100 ul DNA solution. We amplified a 411-bp portion of the I. scapularis 16S rRNA gene by PCR using published primers (Norris et al., 1996) and verified amplification by visualization on agarose gels stained with ethidium bromide. Successfully amplified PCR products were purified using a Qiagen PCR Clean-up Kit before quantification of DNA concentration. I. scapularis 16S amplicons were then bi-directionally sequenced on an ABI 3730xl DNA Analyzer using Big DyeÒ Terminator chemistry. We tested all samples for presence of B. burgdorferi by PCR amplification of a portion of the 16S–23S intergenic spacer region (Liveris et al., 1999; Bunikis et al., 2004) and the outer surface protein C (ospC) gene (Bunikis et al., 2004) with resulting amplicon identity confirmed by sequencing using methods described above. Chromatograms (I. scapularis and B. burgdorferi) were assembled using Sequencher 5.1 and contigs were aligned with reference sequences using MEGA 5.2 (Tamura et al., 2011). We used MacClade 4.08 to identify redundant haplotypes and used Clustal W implemented in MEGA 5.2 to align all unique haplotypes.
195
2.3. Quantitative analysis
165 166 167 168 169
173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193
196 197 198 199 200 201 202 203 204 205 206 207
We used MEGA 5.10 to select among molecular evolutionary models and to construct a maximum likelihood phylogeny (2000 bootstrap replicates) representing all unique haplotypes as well as those haplotypes detected in previous studies (Qiu et al., Q4 2002; Trout et al., 2009; Krakowetz et al., 2011; Table S1). We used Arlequin 3.5 (Excoffier and Lischer, 2010) to explore patterns of population genetic structure. We tested for population structure as both significant differences in allele frequencies among sites and by pairwise FST values that differed significantly from zero. We tested mismatch distributions within individual sites as well as among all samples for departures from models representing demographic and spatial expansion. Finally, we tested for depar-
Q9
tures from selective neutrality using Tajima’s D and Fu’s FS, both of which assess the difference between observed and expected numbers of alleles based on the observed amount of nucleotide diversity. Of these, Fu’s FS is particularly sensitive to departures from neutrality as a result of demographic expansion (Fu, 1997). We constructed bootstrapped (1000 replicates) minimum spanning trees using MSTGold 2.2 (Salipante and Hall, 2011) and visualized the network with Graphvis 1.0 (Gansner and North, 1999).
208
3. Results
216
We collected a total of 303 I. scapularis nymphs by drag sampling in spring and summer 2011 (Brinkerhoff et al. (submitted for publication)). We also collected 1195 adult I. scapularis from hunter-killed deer in fall 2011. Of these, 273 adults were selected for B. burgdorferi screening (Table 1). We were able to amplify 16S rDNA from up to 24 individual ticks per site and all successfully amplified DNA provided unambiguous bi-directional 16S sequence data at 411 nucleotide sites, 29 of which were polymorphic with an average of 2.08 alleles per locus. Haplotype diversity ranged from 0.324 to 0.905 and nucleotide diversity (per locus) ranged from 0.0021 to 0.0083 (Table 1). We detected a total of 24 unique 16S haplotypes (NCBI accession numbers in Table S1, sequence alignment found in Table S2), of which, 13 were previously unreported (Table S1). All but two haplotypes detected in this study (16S14 and 16S18) represented American clade I. scapularis with two individuals, one collected at site YR and the other at site PSP, representing southern clade I. scapularis (Fig. 2, Table S2). Haplotype 16S7 was the most frequently-encountered sequence at all sites and constituted between 23% and 64% of haplotypes at a given site (Table 2). We detected no significant population structure among study sites (Fisher’s exact test of global differentiation, P = 0.24 after 100,000 Markov steps). Moreover, only three out of 36 pairwise FST values were significantly greater than zero at an unadjusted alpha of 0.05: YR:UR (FST = 0.031, P = 0.045), YR:CSF (FST = 0.082, P < 0.001), CSF:UR (FST = 0.083, P = 0.018). We did detect significant departures from selective neutrality at most sampling sites: Fu’s FS was negative at eight of nine sites and significantly lower than zero (alpha = 0.02) at six of nine sites (Table 1). Tajima’s D reflected slightly different patterns with significant differences from zero (alpha = 0.05) at four sites, with two sites having significantly negative values for both indices (Table 1). In accordance with observed values FS values, mismatch analysis of pairwise nucleotide differences revealed departure from a demographic expansion model at only two sites (YR and CSF; Table 1). When all samples were
217
Table 1 Numbers of deer sampled (where relevant), I. scapularis collected and included in the study, along with observed B. burgdorferi infection prevalence. Molecular data include numbers of haplotypes observed, measures of nucleotide diversity, and estimates of departure from selective neutrality and mismatch distributions characteristic of demographic expansion.
* à
Site
Deer sampled (avg. number ticks/deer)
Borrelia burgdorferi infection prevalence
Samples genotyped (haplotypes detected)
Nucleotide diversity (pi) per locus
Segregating sites (S)
Eastern shore NWR (ES)* York river SP (YR)* Crawford SF (CSF) Pocahontas SP (PSP)* University of Richmond (UR) Sabot hill (SH)* Buckingham Co. (BCO)* Nelson Co. (NCO)* Lesesne SF (LSF)
4 (15.0) 8 (53.1) N/A 11 (23.6) N/A 11(12.2) 19 (2.9) 4 (60.3) N/A
2/38 (5.26%) 2/50 (4%) 0/34 (0%)à 8/50 (16%) 2/48 (4.2%)à 6/39 (15.4%) 1/27 (3.7%) 8/69 (11.6%) 43/216 (19.9%)à
22 (7) 23 (10) 17 (3) 24 (9) 19 (6) 24 (10) 24 (10) 22 (7) 24 (10)
0.0025 0.0083 0.0026 0.0058 0.0021 0.0045 0.0046 0.0023 0.0055
6 22 3 18 5 9 9 6 11
Tajima’s D (p-value)
1.31 1.76 1.38 2.06 1.23 1.19 1.59 1.48 1.38
(0.07) (0.03) (0.08) (0.01) (0.11) (0.12) (0.05) (0.04) (0.08)
Fu’s FS (p-value)
Mismatch p-values (demographic, spatial expansion)
3.19 (0.006) 1.76 (0.20) 1.22 (0.744) 2.04 (0.124) 2.81 (0.014) 4.25 (0.009) 4.11 (0.006) 3.54 (0.003) 5.37 (0.001)
0.795, 0.849 0.014, 0.675 0.047, 0.417 0.235, 0.136 0.095, 0.042 0.405, 0.328 0.945, 0.984 0.824, 0.789 0.630, 0.549
Indicate sites from which ticks were collected from deer carcasses. Previously reported in Brinkerhoff et al. (submitted for publication).
Please cite this article in press as: Kelly, R.R., et al. Population genetic structure of the Lyme disease vector Ixodes scapularis at an apparent spatial expansion front. Infect. Genet. Evol. (2014), http://dx.doi.org/10.1016/j.meegid.2014.05.022
209 210 211 212 213 214 215
218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251
MEEGID 1970
No. of Pages 8, Model 5G
4 June 2014 4
R.R. Kelly et al. / Infection, Genetics and Evolution xxx (2014) xxx–xxx
Fig. 2. Maximum likelihood phylogenetic tree of I. scapularis 16S sequences detected in this study, inferred using Tamura’s 3-parameter model (Tamura, 1992) with gammadistributed rates of nucleotide substitution among sites. Haplotypes identified in this study are indicated by numbers 16S1 through 16S24. Reference sequences include NCBI accession number as well as two-letter abbreviation of the state from which samples were collected. Bootstrap values are based on 2000 replicates. American and southern clade designations follow Norris et al., 1996. Accession numbers for sequences presented in this paper as well as source citations for reference sequences are found in Table S1.
252 253 254 255 256 257 258 259
analyzed as a single panmictic population, mismatch analysis revealed a pattern consistent with spatial (goodness-of-fit test P = 0.612), but not demographic (goodness-of-fit test P = 0.0005), population expansion (Fig. 3). We also detected significantly negative D (D = 1.96; P = 0.002) and FS values (FS = 13.4; P < 0.0001), driven by the difference between number of observed (24) and expected (9.14) alleles in the study. To explore the possibility that host-seeking nymphs and adult ticks collected from deer might
represent separate subpopulations, potentially differentiated by host-seeking behavior, we divided our samples into two groups representing collection method and life stage (i.e. dragging for nymphs and collection of adults from deer carcasses). Overall, there was no detectable structure within or among these groups (AMOVA FCT = 0.002, P = 0.65; FST = 0.008, P = 0.13; 5040 permutations) with most (99.25%) of the haplotypic variation coming from within populations. MSTGold 2.2 estimated 2991 possible
Please cite this article in press as: Kelly, R.R., et al. Population genetic structure of the Lyme disease vector Ixodes scapularis at an apparent spatial expansion front. Infect. Genet. Evol. (2014), http://dx.doi.org/10.1016/j.meegid.2014.05.022
260 261 262 263 264 265 266 267
MEEGID 1970
No. of Pages 8, Model 5G
4 June 2014 R.R. Kelly et al. / Infection, Genetics and Evolution xxx (2014) xxx–xxx Table 2 Haplotype frequencies at each of nine field sites. Haplotype accession numbers are found in Table S1 and sequences are found in Table S2. Haplotype
16S1 16S2 16S3 16S4 16S5 16S6 16S7 16S8 16S9 16S10 16S11 16S12 16S13 16S14 16S15 16S16 16S17 16S18 16S19 16S20 16S21 16S22 16S23 16S24
Site ES
YR
CSF
PSP
UR
SH
BCO
NCO
LSF
2 2 1 0 2 1 13 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
2 0 0 0 6 0 8 0 0 0 1 1 1 1 1 1 1 0 0 0 0 0 0 0
0 0 0 0 0 0 14 0 0 0 0 0 0 0 0 0 0 0 0 0 1 2 0 0
0 1 0 1 3 1 13 1 2 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0
0 1 2 1 4 0 10 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 2 2 0 2 0 9 2 2 2 0 0 0 0 0 0 0 0 0 1 1 0 0 0
3 0 1 0 2 0 12 1 1 0 0 0 0 0 0 0 0 0 1 0 0 1 1 1
3 1 1 1 1 1 14 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
2 1 2 0 5 1 6 2 3 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0
Fig. 3. Frequencies of observed pairwise mismatches among haplotypes (gray bars) compared with mismatch distributions expected under demographic expansion (filled triangles, solid black line) or spatial expansion (open diamonds, dashed black line) for all haplotypes at all sites combined. 268 269 270 271 272 273 274 275 276 277 278 279 280 281
minimum spanning trees among the 22 American clade haplotypes we detected; the minimum spanning tree with highest bootstrap support (39% overall) is presented in Fig. 4. As previously reported (Brinkerhoff et al. (submitted for publication)), we detected highest infection prevalence among nymphal samples at site LSF (Table 1). Among adult ticks collected from deer, we detected B. burgdorferi DNA at all sites, albeit mostly at low prevalence (Table 1). Because we would expect higher infection prevalence among adult than nymphal ticks and our dataset included specimens of both life stages, we did not perform any analyses of spatial variation in infection prevalence. We also decided not to analyze spatial patterns of tick parasitism because of generally low sample sizes at each site and high levels of observed variation in parasite burdens among individual deer.
5
4. Discussion
282
Uncommonness or absence of I. scapularis at inland sites in Maryland, Virginia, and North Carolina in previous decades (Sonenshine et al., 1965; Aperson et al., 1990; Amerasinghe et al., 1992, 1993; Diuk-Wasser et al., 2006), combined with observed increase in human LD cases in Virginia from 2008 to 2011 (Fig. 1, Brinkerhoff et al. (submitted for publication)), suggests that I. scapularis populations in western Virginia have undergone recent expansion, either from small relictual populations of endemic ticks, or of ticks that have recently invaded from adjacent states where I. scapularis is more abundant. The population genetic analysis of I. scapularis in Virginia we present here is consistent with recent demographic and spatial expansion throughout most of the state, with the most apparently stable populations occurring in areas with a longer-term association with I. scapularis (Sonenshine et al., 1965, 1995) and Lyme disease cases in humans (Heimberger et al., 1990). Specifically, we found evidence of demographic population expansion at most of our field sites, including all of the sites in the piedmont and mountains, as well as haplotype frequency patterns consistent with spatial population expansion at eight of nine sites and at the state-wide level. We did not find that southeastern Virginia I. scapularis populations were dominated by southern clade haplotypes, suggesting that any potential genetic bases for behavioral differentiation among populations are not well-resolved with the 16S rDNA variation. Overall, our findings support the hypotheses that (1) Virginia populations of I. scapularis tend to be undergoing demographic expansion and that there is spatial expansion occurring throughout much of the state (2) Virginia I. scapularis populations at inland and coastal sites are dominated by American clade tick populations and that genetic variation in this locus is not associated with observed behavioral differences between eastern and western populations. In contrast with a population genetic analysis of I. scapularis 16S sequences at a spatial expansion front in Canada (Krakowetz et al., 2011), we failed to detect significant genetic structure in and among Virginia I. scapularis populations, likely due to the much smaller spatial scale of this study. However, as in previous studies (Krakowetz et al., 2011; Mechai et al., 2013), we did find evidence of recent spatial and demographic population expansion and the haplotype and nucleotide diversity we observed were consistent with these and other reports (i.e. Norris et al., 1996). The concordance between Fu’s FS, which is particularly sensitive to demographic expansion (Fu, 1997), and nucleotide mismatch analysis suggests that tick populations in the central and western – but not eastern – parts of the state are growing. The one exception to this pattern is the easternmost tick population (ES) which showed signs of demographic expansion. This result may be explained by stochastic extinctions in this flood- and storm-surge-prone coastal location, with recolonization by ticks that drop from migratory birds or are imported from deer at nearby coastal locations. Our failure to detect differences in genetic or nucleotide diversity among sites suggests that ticks may not be dispersal-limited. We expected that eastern populations, assumed to be older based on tick surveillance and human Lyme disease case data through the 1990s (Sonenshine et al., 1965; Aperson et al., 1990; Heimberger et al., 1990; Amerasinghe et al., 1992; Amerasinghe et al., 1993; Sonenshine et al., 1995; Diuk-Wasser et al., 2006, 2010), would have accumulated more genotypic diversity than recently established western populations, but it appears that the same overall suite of haplotypes occurs throughout the state and there are no clear spatial patterns in haplotype or nucleotide diversity (Table 1). We note, however, that the two southern clade I. scapularis we detected were found at eastern sites (YR and PSP) associated with the
283
Please cite this article in press as: Kelly, R.R., et al. Population genetic structure of the Lyme disease vector Ixodes scapularis at an apparent spatial expansion front. Infect. Genet. Evol. (2014), http://dx.doi.org/10.1016/j.meegid.2014.05.022
284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345
MEEGID 1970
No. of Pages 8, Model 5G
4 June 2014 6
R.R. Kelly et al. / Infection, Genetics and Evolution xxx (2014) xxx–xxx
Fig. 4. Bootstrapped (1000 replicates) minimum spanning tree containing all American clade haplotypes detected in this study. Bold lines represent >90% bootstrap support, solid lines represent >70% bootstrap support, and dashed lines represent
No. of Pages 8, Model 5G
4 June 2014 Infection, Genetics and Evolution xxx (2014) xxx–xxx 1
Contents lists available at ScienceDirect
Infection, Genetics and Evolution journal homepage: www.elsevier.com/locate/meegid 5 6
Population genetic structure of the Lyme disease vector Ixodes scapularis at an apparent spatial expansion front
3 4 7
Q1
8 9 10
Q2
Rebecca R. Kelly a, David Gaines b, Will F. Gilliam a, R. Jory Brinkerhoff a,c,⇑ a
11 12
Department of Biology, University of Richmond, 28 Westhampton Way, Richmond, VA 23173, United States Division of Environmental Epidemiology, Virginia Department of Health, Richmond, VA 23219, United States c School of Life Sciences, University of KwaZulu-Natal, Pietermaritzburg, South Africa b
a r t i c l e
1 2 4 8 15 16 17 18 19
i n f o
Article history: Received 16 January 2014 Received in revised form 18 May 2014 Accepted 21 May 2014 Available online xxxx
20 21 22 23 24 25 26 27
Keywords: Disease ecology Range expansion Acari Zoonosis Tick-borne disease Epidemiology
a b s t r a c t Modeling and empirical evidence suggests that Lyme disease is undergoing geographic expansion from principal foci in the midwestern and northeastern United States. Virginia is at the southern edge of the current expansion zone and has seen dramatic rise in human Lyme disease cases since 2007, potentially owing to a recent increase in vector abundance. Ixodes scapularis is known throughout the eastern US but behavioral or physiological variation between northern and southern lineages might lead northern-variant ticks to more frequently parasitize humans. We hypothesized that recent spatial and numerical increase in Lyme disease cases is associated with demographic and/or spatial expansion of I. scapularis and that signals of these phenomena would be detectable and discernable in population genetic signals. In summer and fall 2011, we collected nymphal I. scapularis by drag sampling and adult I. scapularis from deer carcasses at hunting check stations at nine sites arranged along an east–west transect through central Virginia. We analyzed 16S mtDNA sequences data from up to 24 I. scapularis individuals collected from each site and detected a total of 24 haplotypes containing 29 segregating sites. We found no evidence for population genetic structure among these sites but we did find strong signals of both demographic and spatial expansion throughout our study system. We found two haplotypes (one individual each) representing a lineage of ticks that is only found in the southeastern United States, with the remaining individuals representing a less genetically diverse clade that is typical of the northern United States, but that has also been detected in the American South. Taken together, these results lead us to conclude that I. scapularis populations in Virginia are expanding and that this expansion may account for recent observed increases in Lyme disease. Ó 2014 Published by Elsevier B.V.
29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
51 52
1. Introduction Vector-borne zoonoses are inherently complex disease systems that are keenly sensitive to anthropogenic and natural environmental change (Mills et al., 2010; Thompson, 2013). The geographic distributions of tick-borne pathogens, as well as phenology of host seeking ticks, have been linked to climate change (reviewed in Gage et al., 2008) and there have been numerous reports of tick-borne diseases increasing in incidence and spatial extent (Danielova et al., 2006; Materna et al., 2008; Medlock et al., 2013). Although large-scale surveillance efforts to generate risk maps for tick-borne pathogens have been undertaken, such efforts are highly resource-intensive and the static models generated by these efforts provide only a ‘snapshot’ image of dynamic
53 54 55 56 57 58 59 60 61 62 63 64
⇑ Corresponding author at: Department of Biology, University of Richmond, 28
Q3
Westhampton Way, Richmond, VA 23173, United States. Tel.: +1 18044841592. E-mail address: [email protected] (R. Jory Brinkerhoff).
processes (Diuk-Wasser et al., 2010, 2012). Molecular approaches, however, can lend insights into evolutionary dynamics and can be used to infer population histories and pathogen transmission dynamics (Guhl and Ramirez, 2013). The spread of the Lyme disease agent, Borrelia burgdorferi, and its principal vector in eastern North America (Ixodes scapularis), into Canada has justifiably been the focus of a great deal of attention (e.g. Koffi et al., 2012; Leighton et al., 2012; Ogden et al., 2009, 2013; Wu et al., 2013). The population genetic structure of I. scapularis in Canada suggests recent demographic expansion (Krakowetz et al., 2011; Mechai et al., 2013) and the genotypic heterogeneity among Canadian B. burgdorferi samples revealed the same suite of genotypes in both southeastern Canada and the eastern United States, albeit with Canadian populations showing signals of founder events (Ogden et al., 2011). There has been less attention paid, however, to other potential Lyme disease expansion fronts. In a large scale study of I. scapularis occurrence throughout the eastern United States, Diuk-Wasser et al. (2010, 2012)
http://dx.doi.org/10.1016/j.meegid.2014.05.022 1567-1348/Ó 2014 Published by Elsevier B.V.
Please cite this article in press as: Kelly, R.R., et al. Population genetic structure of the Lyme disease vector Ixodes scapularis at an apparent spatial expansion front. Infect. Genet. Evol. (2014), http://dx.doi.org/10.1016/j.meegid.2014.05.022
65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82
MEEGID 1970
No. of Pages 8, Model 5G
4 June 2014 2 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125
R.R. Kelly et al. / Infection, Genetics and Evolution xxx (2014) xxx–xxx
developed environmental models to estimate tick population density and Lyme disease risk and found evidence consistent with overall spatial expansion of tick populations. Hamer et al., 2010 demonstrated apparent northward movement of I. scapularis along the eastern shore of Lake Michigan over a period of four years (Hamer et al., 2010). Lee et al. (2013) compared tick occurrence on hunter-shot deer from 1981 to 2009 and concluded that I. scapularis expanding from west to east throughout Wisconsin. Brinkerhoff et al. (submitted for publication) showed data consistent with a shifting geographic distribution of Lyme disease cases in Virginia and suggested, based on field sampling and vector phylogenetic analysis, that tick populations in the western part of the state may be larger than they had been in previous decades. The presence of I. scapularis and incidence of Lyme disease in the southern United States is highly contentious; the occurrence of Lymelike disease in areas where I. scapularis is uncommonly encountered and the importance of other putative disease agents such as Borrelia lonestari contribute to the confusion about tick-borne disease in southern and southeastern states (Stromdahl and Hickling, 2012; Lantos et al., 2013). As a result, it is critical to explore the potential southward expansion of Lyme disease into states from which it has historically been rare or absent. The evolutionary history of I. scapularis, inferred from variation in mitochondrial rDNA and other markers, suggests that this species was relegated to southern North America during Pleistocene glaciation events, and that northern North America was recolonized by founding I. scapularis populations that moved north and diversified following glacial retreat (Norris et al., 1996; Qiu et al., 2002; Humphrey et al., 2010; van Zee et al., 2013). The purported result of this movement is higher genetic diversity of I. scapularis in the southern United States, including a more genetically diverse lineage (hereafter referred to as ‘‘southern clade’’) that is only found in the southern US, and a second lineage (hereafter referred to as ‘‘American clade’’), that is found in both northern and southern states (Norris et al., 1996; Mixson et al., 2004). Potential differences between these two lineages that may affect host seeking or other behaviors (e.g. Durden et al., 2002) could account for differences in infection prevalence and risk to Lyme disease between northern and southern states (Diuk-Wasser et al., 2012; Stromdahl and Hickling, 2012). Since 2007, increasing Lyme disease cases have been reported in Virginia and Brinkerhoff et al. (submitted for publication) demonstrated that Virginia counties with recent increases in human
Lyme disease cases are associated with relatively high I. scapularis densities. Our goal was to use population genetic analysis of Ixodes scapularis sampled throughout central Virginia to identify signals of recent spatial or demographic population expansion. We expected to find relatively lower genetic diversity in ticks associated with recent increases in human LD (i.e. western VA) and higher diversity in locations that have historically been associated with I. scapularis populations (eastern VA). Similarly, we expected to see population genetic structure at the state-wide level and evidence of recent spatial and/or demographic expansion at western sites.
126
2. Materials and methods
137
2.1. Tick collection
138
Tick collection sites were selected along an east–west transect from the Atlantic coast to high-elevation counties in west-central Virginia (Fig. 1). These sites represent variation in recent LD emergence in humans with the westernmost sites having seen dramatic increases in incidence since the surveillance case definition for Virginia was modified to require two-tier testing (i.e. enzymelinked immunosorbent assay and immunoglobulin M-positive Western blot) in 2008 (Fig. 1, Brinkerhoff et al. (submitted for publication)). Host-seeking I. scapularis nymphs were collected by drag sampling from three sites summer 2011 (Fig. 1; Brinkerhoff et al., submitted for publication). At each of these sites, we haphazardly selected and established five 100 m transects through deciduous or mixed evergreen-deciduous forest and collected ticks by dragging a 1 1 m corduroy cloth along both sides of each transect for a total of 1000 square meters of sampling area (Falco and Fish, 1992). For each transect sampling event, we checked the drag cloth at 20-meter intervals and placed all observed ticks in uniquelynumbered vials containing a 70% ethanol solution. Each site was visited four times between May and July 2011 (Brinkerhoff et al. (submitted for publication)). Adult I. scapularis were collected from hunter-shot deer at six additional sites in October and November 2011 (Fig. 1). During these collection events, we carefully removed as many adult I. scapularis as possible while hunters reported to check-in stations. Host-seeking ticks were stored in 70% ethanol and adult ticks were stored in 70% ethanol at 4 °C until identification by light microscopy using dichotomous keys (Sonenshine,
139
Fig. 1. Locations of field sites sampled in this study. Triangles represent sites from which adult I. scapularis were collected at deer hunts and diamonds represent sites where I. scapularis nymphs were collected by drag sampling. Summary statistics and sample sizes associated with each site are indicated in Table 1. County shading represents the mean annual change in Lyme disease incidence between 2008, when the surveillance case definition was made more rigorous in Virginia, and 2011 (incidence data from Brinkerhoff et al. (submitted for publication)).
Please cite this article in press as: Kelly, R.R., et al. Population genetic structure of the Lyme disease vector Ixodes scapularis at an apparent spatial expansion front. Infect. Genet. Evol. (2014), http://dx.doi.org/10.1016/j.meegid.2014.05.022
127 128 129 130 131 132 133 134 135 136
140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164
MEEGID 1970
No. of Pages 8, Model 5G
4 June 2014 3
R.R. Kelly et al. / Infection, Genetics and Evolution xxx (2014) xxx–xxx
170
1979). Based on availability, adult ticks including males and the least engorged females were selected from each hunt for B. burgdorferi testing (when possible, approximately equal numbers of male and female ticks were selected from each site). We subsequently used a random number generator to select a subset of up to 25 ticks per site for genotypic analysis.
171
2.2. Molecular methods
172
194
Following flash-freezing in liquid nitrogen, we pulverized ticks in individual microcentrifuge tubes with a sterile plastic pestle. We extracted DNA from individual I. scapularis using Qiagen DNeasy Blood and Tissue Kit following the manufacturer’s protocols except we performed two 50 ul DNA elutions per sample for a total of 100 ul DNA solution. We amplified a 411-bp portion of the I. scapularis 16S rRNA gene by PCR using published primers (Norris et al., 1996) and verified amplification by visualization on agarose gels stained with ethidium bromide. Successfully amplified PCR products were purified using a Qiagen PCR Clean-up Kit before quantification of DNA concentration. I. scapularis 16S amplicons were then bi-directionally sequenced on an ABI 3730xl DNA Analyzer using Big DyeÒ Terminator chemistry. We tested all samples for presence of B. burgdorferi by PCR amplification of a portion of the 16S–23S intergenic spacer region (Liveris et al., 1999; Bunikis et al., 2004) and the outer surface protein C (ospC) gene (Bunikis et al., 2004) with resulting amplicon identity confirmed by sequencing using methods described above. Chromatograms (I. scapularis and B. burgdorferi) were assembled using Sequencher 5.1 and contigs were aligned with reference sequences using MEGA 5.2 (Tamura et al., 2011). We used MacClade 4.08 to identify redundant haplotypes and used Clustal W implemented in MEGA 5.2 to align all unique haplotypes.
195
2.3. Quantitative analysis
165 166 167 168 169
173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193
196 197 198 199 200 201 202 203 204 205 206 207
We used MEGA 5.10 to select among molecular evolutionary models and to construct a maximum likelihood phylogeny (2000 bootstrap replicates) representing all unique haplotypes as well as those haplotypes detected in previous studies (Qiu et al., Q4 2002; Trout et al., 2009; Krakowetz et al., 2011; Table S1). We used Arlequin 3.5 (Excoffier and Lischer, 2010) to explore patterns of population genetic structure. We tested for population structure as both significant differences in allele frequencies among sites and by pairwise FST values that differed significantly from zero. We tested mismatch distributions within individual sites as well as among all samples for departures from models representing demographic and spatial expansion. Finally, we tested for depar-
Q9
tures from selective neutrality using Tajima’s D and Fu’s FS, both of which assess the difference between observed and expected numbers of alleles based on the observed amount of nucleotide diversity. Of these, Fu’s FS is particularly sensitive to departures from neutrality as a result of demographic expansion (Fu, 1997). We constructed bootstrapped (1000 replicates) minimum spanning trees using MSTGold 2.2 (Salipante and Hall, 2011) and visualized the network with Graphvis 1.0 (Gansner and North, 1999).
208
3. Results
216
We collected a total of 303 I. scapularis nymphs by drag sampling in spring and summer 2011 (Brinkerhoff et al. (submitted for publication)). We also collected 1195 adult I. scapularis from hunter-killed deer in fall 2011. Of these, 273 adults were selected for B. burgdorferi screening (Table 1). We were able to amplify 16S rDNA from up to 24 individual ticks per site and all successfully amplified DNA provided unambiguous bi-directional 16S sequence data at 411 nucleotide sites, 29 of which were polymorphic with an average of 2.08 alleles per locus. Haplotype diversity ranged from 0.324 to 0.905 and nucleotide diversity (per locus) ranged from 0.0021 to 0.0083 (Table 1). We detected a total of 24 unique 16S haplotypes (NCBI accession numbers in Table S1, sequence alignment found in Table S2), of which, 13 were previously unreported (Table S1). All but two haplotypes detected in this study (16S14 and 16S18) represented American clade I. scapularis with two individuals, one collected at site YR and the other at site PSP, representing southern clade I. scapularis (Fig. 2, Table S2). Haplotype 16S7 was the most frequently-encountered sequence at all sites and constituted between 23% and 64% of haplotypes at a given site (Table 2). We detected no significant population structure among study sites (Fisher’s exact test of global differentiation, P = 0.24 after 100,000 Markov steps). Moreover, only three out of 36 pairwise FST values were significantly greater than zero at an unadjusted alpha of 0.05: YR:UR (FST = 0.031, P = 0.045), YR:CSF (FST = 0.082, P < 0.001), CSF:UR (FST = 0.083, P = 0.018). We did detect significant departures from selective neutrality at most sampling sites: Fu’s FS was negative at eight of nine sites and significantly lower than zero (alpha = 0.02) at six of nine sites (Table 1). Tajima’s D reflected slightly different patterns with significant differences from zero (alpha = 0.05) at four sites, with two sites having significantly negative values for both indices (Table 1). In accordance with observed values FS values, mismatch analysis of pairwise nucleotide differences revealed departure from a demographic expansion model at only two sites (YR and CSF; Table 1). When all samples were
217
Table 1 Numbers of deer sampled (where relevant), I. scapularis collected and included in the study, along with observed B. burgdorferi infection prevalence. Molecular data include numbers of haplotypes observed, measures of nucleotide diversity, and estimates of departure from selective neutrality and mismatch distributions characteristic of demographic expansion.
* à
Site
Deer sampled (avg. number ticks/deer)
Borrelia burgdorferi infection prevalence
Samples genotyped (haplotypes detected)
Nucleotide diversity (pi) per locus
Segregating sites (S)
Eastern shore NWR (ES)* York river SP (YR)* Crawford SF (CSF) Pocahontas SP (PSP)* University of Richmond (UR) Sabot hill (SH)* Buckingham Co. (BCO)* Nelson Co. (NCO)* Lesesne SF (LSF)
4 (15.0) 8 (53.1) N/A 11 (23.6) N/A 11(12.2) 19 (2.9) 4 (60.3) N/A
2/38 (5.26%) 2/50 (4%) 0/34 (0%)à 8/50 (16%) 2/48 (4.2%)à 6/39 (15.4%) 1/27 (3.7%) 8/69 (11.6%) 43/216 (19.9%)à
22 (7) 23 (10) 17 (3) 24 (9) 19 (6) 24 (10) 24 (10) 22 (7) 24 (10)
0.0025 0.0083 0.0026 0.0058 0.0021 0.0045 0.0046 0.0023 0.0055
6 22 3 18 5 9 9 6 11
Tajima’s D (p-value)
1.31 1.76 1.38 2.06 1.23 1.19 1.59 1.48 1.38
(0.07) (0.03) (0.08) (0.01) (0.11) (0.12) (0.05) (0.04) (0.08)
Fu’s FS (p-value)
Mismatch p-values (demographic, spatial expansion)
3.19 (0.006) 1.76 (0.20) 1.22 (0.744) 2.04 (0.124) 2.81 (0.014) 4.25 (0.009) 4.11 (0.006) 3.54 (0.003) 5.37 (0.001)
0.795, 0.849 0.014, 0.675 0.047, 0.417 0.235, 0.136 0.095, 0.042 0.405, 0.328 0.945, 0.984 0.824, 0.789 0.630, 0.549
Indicate sites from which ticks were collected from deer carcasses. Previously reported in Brinkerhoff et al. (submitted for publication).
Please cite this article in press as: Kelly, R.R., et al. Population genetic structure of the Lyme disease vector Ixodes scapularis at an apparent spatial expansion front. Infect. Genet. Evol. (2014), http://dx.doi.org/10.1016/j.meegid.2014.05.022
209 210 211 212 213 214 215
218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251
MEEGID 1970
No. of Pages 8, Model 5G
4 June 2014 4
R.R. Kelly et al. / Infection, Genetics and Evolution xxx (2014) xxx–xxx
Fig. 2. Maximum likelihood phylogenetic tree of I. scapularis 16S sequences detected in this study, inferred using Tamura’s 3-parameter model (Tamura, 1992) with gammadistributed rates of nucleotide substitution among sites. Haplotypes identified in this study are indicated by numbers 16S1 through 16S24. Reference sequences include NCBI accession number as well as two-letter abbreviation of the state from which samples were collected. Bootstrap values are based on 2000 replicates. American and southern clade designations follow Norris et al., 1996. Accession numbers for sequences presented in this paper as well as source citations for reference sequences are found in Table S1.
252 253 254 255 256 257 258 259
analyzed as a single panmictic population, mismatch analysis revealed a pattern consistent with spatial (goodness-of-fit test P = 0.612), but not demographic (goodness-of-fit test P = 0.0005), population expansion (Fig. 3). We also detected significantly negative D (D = 1.96; P = 0.002) and FS values (FS = 13.4; P < 0.0001), driven by the difference between number of observed (24) and expected (9.14) alleles in the study. To explore the possibility that host-seeking nymphs and adult ticks collected from deer might
represent separate subpopulations, potentially differentiated by host-seeking behavior, we divided our samples into two groups representing collection method and life stage (i.e. dragging for nymphs and collection of adults from deer carcasses). Overall, there was no detectable structure within or among these groups (AMOVA FCT = 0.002, P = 0.65; FST = 0.008, P = 0.13; 5040 permutations) with most (99.25%) of the haplotypic variation coming from within populations. MSTGold 2.2 estimated 2991 possible
Please cite this article in press as: Kelly, R.R., et al. Population genetic structure of the Lyme disease vector Ixodes scapularis at an apparent spatial expansion front. Infect. Genet. Evol. (2014), http://dx.doi.org/10.1016/j.meegid.2014.05.022
260 261 262 263 264 265 266 267
MEEGID 1970
No. of Pages 8, Model 5G
4 June 2014 R.R. Kelly et al. / Infection, Genetics and Evolution xxx (2014) xxx–xxx Table 2 Haplotype frequencies at each of nine field sites. Haplotype accession numbers are found in Table S1 and sequences are found in Table S2. Haplotype
16S1 16S2 16S3 16S4 16S5 16S6 16S7 16S8 16S9 16S10 16S11 16S12 16S13 16S14 16S15 16S16 16S17 16S18 16S19 16S20 16S21 16S22 16S23 16S24
Site ES
YR
CSF
PSP
UR
SH
BCO
NCO
LSF
2 2 1 0 2 1 13 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
2 0 0 0 6 0 8 0 0 0 1 1 1 1 1 1 1 0 0 0 0 0 0 0
0 0 0 0 0 0 14 0 0 0 0 0 0 0 0 0 0 0 0 0 1 2 0 0
0 1 0 1 3 1 13 1 2 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0
0 1 2 1 4 0 10 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 2 2 0 2 0 9 2 2 2 0 0 0 0 0 0 0 0 0 1 1 0 0 0
3 0 1 0 2 0 12 1 1 0 0 0 0 0 0 0 0 0 1 0 0 1 1 1
3 1 1 1 1 1 14 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
2 1 2 0 5 1 6 2 3 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0
Fig. 3. Frequencies of observed pairwise mismatches among haplotypes (gray bars) compared with mismatch distributions expected under demographic expansion (filled triangles, solid black line) or spatial expansion (open diamonds, dashed black line) for all haplotypes at all sites combined. 268 269 270 271 272 273 274 275 276 277 278 279 280 281
minimum spanning trees among the 22 American clade haplotypes we detected; the minimum spanning tree with highest bootstrap support (39% overall) is presented in Fig. 4. As previously reported (Brinkerhoff et al. (submitted for publication)), we detected highest infection prevalence among nymphal samples at site LSF (Table 1). Among adult ticks collected from deer, we detected B. burgdorferi DNA at all sites, albeit mostly at low prevalence (Table 1). Because we would expect higher infection prevalence among adult than nymphal ticks and our dataset included specimens of both life stages, we did not perform any analyses of spatial variation in infection prevalence. We also decided not to analyze spatial patterns of tick parasitism because of generally low sample sizes at each site and high levels of observed variation in parasite burdens among individual deer.
5
4. Discussion
282
Uncommonness or absence of I. scapularis at inland sites in Maryland, Virginia, and North Carolina in previous decades (Sonenshine et al., 1965; Aperson et al., 1990; Amerasinghe et al., 1992, 1993; Diuk-Wasser et al., 2006), combined with observed increase in human LD cases in Virginia from 2008 to 2011 (Fig. 1, Brinkerhoff et al. (submitted for publication)), suggests that I. scapularis populations in western Virginia have undergone recent expansion, either from small relictual populations of endemic ticks, or of ticks that have recently invaded from adjacent states where I. scapularis is more abundant. The population genetic analysis of I. scapularis in Virginia we present here is consistent with recent demographic and spatial expansion throughout most of the state, with the most apparently stable populations occurring in areas with a longer-term association with I. scapularis (Sonenshine et al., 1965, 1995) and Lyme disease cases in humans (Heimberger et al., 1990). Specifically, we found evidence of demographic population expansion at most of our field sites, including all of the sites in the piedmont and mountains, as well as haplotype frequency patterns consistent with spatial population expansion at eight of nine sites and at the state-wide level. We did not find that southeastern Virginia I. scapularis populations were dominated by southern clade haplotypes, suggesting that any potential genetic bases for behavioral differentiation among populations are not well-resolved with the 16S rDNA variation. Overall, our findings support the hypotheses that (1) Virginia populations of I. scapularis tend to be undergoing demographic expansion and that there is spatial expansion occurring throughout much of the state (2) Virginia I. scapularis populations at inland and coastal sites are dominated by American clade tick populations and that genetic variation in this locus is not associated with observed behavioral differences between eastern and western populations. In contrast with a population genetic analysis of I. scapularis 16S sequences at a spatial expansion front in Canada (Krakowetz et al., 2011), we failed to detect significant genetic structure in and among Virginia I. scapularis populations, likely due to the much smaller spatial scale of this study. However, as in previous studies (Krakowetz et al., 2011; Mechai et al., 2013), we did find evidence of recent spatial and demographic population expansion and the haplotype and nucleotide diversity we observed were consistent with these and other reports (i.e. Norris et al., 1996). The concordance between Fu’s FS, which is particularly sensitive to demographic expansion (Fu, 1997), and nucleotide mismatch analysis suggests that tick populations in the central and western – but not eastern – parts of the state are growing. The one exception to this pattern is the easternmost tick population (ES) which showed signs of demographic expansion. This result may be explained by stochastic extinctions in this flood- and storm-surge-prone coastal location, with recolonization by ticks that drop from migratory birds or are imported from deer at nearby coastal locations. Our failure to detect differences in genetic or nucleotide diversity among sites suggests that ticks may not be dispersal-limited. We expected that eastern populations, assumed to be older based on tick surveillance and human Lyme disease case data through the 1990s (Sonenshine et al., 1965; Aperson et al., 1990; Heimberger et al., 1990; Amerasinghe et al., 1992; Amerasinghe et al., 1993; Sonenshine et al., 1995; Diuk-Wasser et al., 2006, 2010), would have accumulated more genotypic diversity than recently established western populations, but it appears that the same overall suite of haplotypes occurs throughout the state and there are no clear spatial patterns in haplotype or nucleotide diversity (Table 1). We note, however, that the two southern clade I. scapularis we detected were found at eastern sites (YR and PSP) associated with the
283
Please cite this article in press as: Kelly, R.R., et al. Population genetic structure of the Lyme disease vector Ixodes scapularis at an apparent spatial expansion front. Infect. Genet. Evol. (2014), http://dx.doi.org/10.1016/j.meegid.2014.05.022
284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345
MEEGID 1970
No. of Pages 8, Model 5G
4 June 2014 6
R.R. Kelly et al. / Infection, Genetics and Evolution xxx (2014) xxx–xxx
Fig. 4. Bootstrapped (1000 replicates) minimum spanning tree containing all American clade haplotypes detected in this study. Bold lines represent >90% bootstrap support, solid lines represent >70% bootstrap support, and dashed lines represent
