Tick Collection

The 11 closed-canopy forested sites included in this study are distributed throughout the Lehigh Valley region of eastern PA, mostly within Lehigh County (Fig. 1). Most of the sampling was performed at the same sites as Edwards et al. (2015). Sites were selected in 2015 for nymphal tick density (ability to collect >50 I. scapularis nymphs within 2 hr), public access, and protected status. Each site contains at least 5 hectares of tulip tree and red oak-dominated forest. Ticks were collected by dragging with a 1 m² white corduroy cloth. Over 450 drags, each timed for 30 min, were performed at a slow walking pace. An average of 11 drags (SD = 6), each ∼620 m in length, were performed at each site, each year, except for RP which was dropped from the study in 2018 due to an absence of ticks. GPX files were collected for individual drags with the iPhone Trails Application Version 2.13 (iosphere GmbH, Cologne, Germany). Time while moving was estimated with GPXSee version 5.17 to adjust for time spent collecting ticks from the drag cloth. Sites were visited at least once per year, and GPS data were used to replicate the approximate sampling locations among years at each site. Dragging was generally performed in leaf litter in the close vicinity of trails so that our sampling would represent places that were likely to be encountered by the public and contain questing nymphs, based on our previous sampling experience (Edwards et al. 2015). All sampling was between May 24 and June 24 (Weeks 21–25) of each year in 2015, 2016, and 2017 and between June 6 and July 11 (Weeks 23–28) in 2018. Both periods coincide with the peak of nymphal activity observed in mid-western PA (Simmons et al. 2015). Sampling was generally performed between 10:00 and 15:00 each day and never when the ground was wet. Ticks were transferred into separate vials for each drag, placed into 80% ethanol and stored at -20°C until they were processed for DNA extraction or at -80°C for long-term storage.

DNA Isolation

DNA was extracted from a subset of ∼50 nymphs per site, per summer 2015–2017. Within each site, an equal proportion of nymphs were selected from each of the vials representing separate drags. DNA was extracted from individual nymphs with the bead-beating protocol of Crowder et al. (2010) and the QIAamp DNA Mini Kit (Qiagen, Valencia, CA) as described by Edwards et al. (2015) and validated by Ammazzalorso et al. (2015). Tick DNA yield was estimated with a nuclear ribosomal internal transcribed spacer (ITS2) qPCR assay (Strube et al. 2010). Measures to avoid sample contamination were consistent with Piedmonte (2018). DNA extractions were performed in groups of 15 plus one blank in which a nymph was not included. If on rare occasions, tick DNA was detected in a blank (Ct ≤40), all of the samples in the group were discarded. An average yield of 55 ng total DNA per nymph (Av. Ct = 20.9, SD 1.4) was obtained, as estimated from an ITS2 qPCR standard curve. Our yields were consistent with Ammazzalorso et al. (2015), who used a similar method. Samples with less than ∼3.2 ng DNA (ITS2 qPCR Ct ≥ 25) per nymph were excluded from the study. DNA samples were kept at -20°C or on ice until use and then at -80°C for long-term storage.

Pathogen Testing

Three molecular methods (qPCR, High-Resolution Melting (HRM) assays, and DNA sequencing) were used to identify and differentiate between six microbes (B. burgdorferi, B. miyamotoi, two strain variants of A. phagocytophilum, B. microti, and Babesia odocoilei (Piroplasmida: Babesiidae)) to a level of specificity consistent with CDC recommendations (CDC 2019b). Each sample that was considered ‘positive’ for any microbe was tested in two different ways, either by two distinct qPCR assays, a combination of qPCR and HRM or by qPCR and DNA sequencing. All individual nymphal samples were screened with a Borrelia-specific 23S rDNA (23S) qPCR assay from Courtney et al. (2004). All 23S qPCR-positive samples were re-tested with a B. burgdorferi-specific outer surface protein A (OspA) qPCR assay from Dibernardo et al. (2014). All 23S qPCR-positive and OspA qPCR-negative ticks were re-tested with a B. miyamotoi-specific glycerophosphodiester phosphodiesterase (GlpQ) qPCR assay (Dibernardo et al. 2014). Samples that were 23S-positive, OspA-negative, and GlpQ-negative by qPCR assays were re-tested by an HRM assay to differentiate between B. burgdorferi and B. miyamotoi. All samples were screened with an A. phagocytophilum major outer membrane protein (Msp2) qPCR assay (Courtney et al. 2004) in pools of five samples. If a pool tested positive, then each sample within the pool was tested individually with the same Msp2 qPCR assay. Since the Msp2 qPCR assay does not differentiate between A. phagocytophilum strain variants (Ap-V1 and Ap-ha), DNA sequencing of a 16S rDNA (16S) nested conventional PCR fragment (Massung et al. 1998) was used to distinguish between strain variants for Msp2 qPCR-positive samples. All samples were screened with a B. microti-specific 18S rDNA (18S) qPCR assay (Teal et al. 2012) with pools of five ticks. If a pool tested positive, then each sample within the pool was tested individually with the same 18S qPCR assay. All Babesia 18S qPCR-positive samples were either re-tested with a duplex Babesia 18S HRM assay or by sequencing a region of the 18S rDNA gene that differs between B. microti and B. odocoilei.

Fig. 1.

Location of study sites. AT: Alburtis Mountain Road Tract (40.506, -75.598); LH: Little Lehigh Headwaters Wildlife Park (40.480, -75.694); BP: B. Leroy and Elizabeth Burkhart Preserve (40.516, -75.487); SP: South Mountain Preserve, Wildlands Conservancy (40.547, -75.476); GA: Graver Arboretum of Muhlenberg College (40.800, -75.359); LM: Lehigh Mountain Park 40.604, -75.423); RM: Riemert Memorial Bird Haven, Wildlands Conservancy (40.502, -75.563); RP: Raker Preserve, Muhlenberg College (40.692, -75.706); SM: South Mountain Park (40.599, -75.371); SW: Scholl Woodlands Preserve (40.564, -75.434); TP: Trexler Preserve (40.658, -75.622). Inset, the three shaded counties of Pennsylvania are from right to left, Berks, Lehigh, and Northampton.

qPCR Assays

All of the qPCR primers, probe sequences, and methods used in this study are the same as Edwards et al. (2015). All qPCR assays were performed in 20 μl reactions with an Applied Biosystems StepOne Plus thermal cycler with a fast cycling program. A 2-μl aliquot was sampled from each tick (∼0.02 tick equivalents), containing ∼1 ng of total DNA. All reactions were performed with TaqMan Fast Advanced Master Mix (Life Technologies, Carlsbad, CA) or an equivalent formulation (PrimeTime, Integrated DNA Technologies, Coralville, IA). Standard TaqMan Gene Expression Assays (Life Technologies) were used according to the manufacturer's specifications. Assays were generally considered negative if they had not reached threshold Ct by 40 cycles. In the few cases where typically sigmoid-shaped curves were observed after 40 cycles, the samples were re-tested with nested PCR followed by sequencing, or with an HRM assay. The positive control for B. burgdorferi was genomic DNA from strain B31 (35210D) purchased from the ATCC (Manassas, VA). The B. microti positive control was total DNA purified from B. microti infected mouse blood (ATCC 30222). The Anaplasma positive control was total DNA from an A. phagocytophilum-infected tick (Ap-variant 1), kindly provided by Michael Hutchinson of the PA Department of Environmental Protection, Bureau of Labs-Vector Management, Harrisburg, PA. Our positive control for B. miyamotoi was DNA from a B. miyamotoi infected tick, kindly provided by the Northeast Wildlife DNA Laboratory at East Stroudsburg University, East Stroudsburg, PA.

High-Resolution Melt Analysis

A duplex Borrelia spp. integral outer surface membrane protein p66 (p66) HRM assay was developed to distinguish between B. burgdorferi and B. miyamotoi. This HRM assay was used for samples that were Borrelia 23S qPCR positive but OspA-negative and GlpQ negative in qPCR assays. The Tm of the respective Borrelia spp. p66 amplicons are predicted by uMelt (2019) to differ by ∼3°C, as they share only 45/123 identities and differ in total GC content by 6%. Five GlpQ qPCR-positive samples were used to validate the p66 HRM assay for B. miyamotoi, with an average Tm = 73.9°C (SD 0.5). Six OspA qPCR-positive samples were used to validate the p66 HRM assay for B. burgdorferi, with an average melting temperature (Tm) = 77.3°C (SD 0.2). A Babesia spp. 18S HRM assay with a single primer pair was developed to distinguish B. odocoilei from B. microti. The Tm of the respective Babesia spp. 18S rDNA amplicons are predicted by uMelt (2019) to differ by ∼4°C as they share only 95/119 identities and differ in total GC content by 6%. A B. odocoilei sample that was identified by DNA sequencing was used to validate the assay, with a reproducible Tm value of 81.6°C. Seven B. microti samples that were identified by DNA sequencing had an average Tm value of 79.2°C (SD 0.1). Primers used in both HRM assays are in Table 1. HRM assays were performed with PowerUp TM SYBR Green MasterMix (Applied Biosystems, Foster City, CA) using a StepOnePlus RealTime PCR system (Applied Biosystems) in a 20 μl volume with 0.5 mM of each primer. Forty cycles of amplification were followed by a melt curve analysis, in which the temperature was increased from 65 to 95°C by 0.2°C, at 10s intervals. Dissociation curves and their derivatives were generated by StepOne software.

DNA Sequencing

PCR products from B. miyamotoi and B. microti positive nymphs were amplified and sequenced with methods and primers that were described by Edwards et al. (2015). Anaplasma samples were identified to strain with nested conventional PCR reactions that amplified a fragment of the 16S rDNA (16S) sequence (Massung et al. 1998) and were sequenced with an optimized primer (Table 1) that was designed using the Primer3 program (Untergasser et al. 2012) by Genewiz, South Plainfield, NJ. Variations at two positions were analyzed to differentiate between A. phagocytophilum variants (Massung et al. 1998).

Statistical Analysis of DON

A Bayesian negative binomial regression analysis was used to model I. scapularis DON. Negative binomial regression was selected in place of Poisson regression to account for overdispersion in the nymph counts. Two regression models were analyzed. The number of nymphs/100 m² per day was the response variable in both models. For the first model, only year was used as a categorical predictor. Bayesian 95% highest posterior density intervals for each parameter were calculated. For each of the regression coefficients, an independent normal distribution with mean 0 and standard deviation 10 was used for the prior. For the overdispersion parameter, a gamma distribution with shape parameter 0.5 and rate parameter 0.1 was selected as the prior. Posterior distributions for all nymphal density parameters were estimated using Markov chain Monte Carlo Methods implemented in R (R Team Core 2018) with the bayesm (Rossi 2018) package. Convergence of the chains was determined by monitoring the Monte Carlo standard errors with the batch means procedures described in Flegal et al. (2008). Several starting values were used to initiate the Markov chain for each parameter to ensure convergence of the chains. After an initial burn-in period of 500,000 samples, 4,500,000 samples of the posterior were recorded. The size of the Markov chain was increased to ensure there were at least 100,000 effective samples from the posterior of each parameter. Inference was performed again with more disperse prior distributions for the regression coefficients, but there were no significant changes in the results. A second Bayesian negative binomial regression model was used to assess the impact of weather variables. This model considered the year, day of the year, and the following additional weather predictors: accumulated degree-days above 10°C since 1 March, precipitation for the past 30 d, daytime relative humidity for the past 7 d (RH7) and 8–14 d prior to the sampling date (RH14). Meteorological data were from the National Weather Service station at Lehigh Valley International Airport (KABE).

Table 1.

Primers unique to this study

Statistical Analysis of Nymphal Infection Prevalence

Bayesian 95% credible intervals (95 CI) for nymphal infection prevalence were constructed in aggregate across years at every observed site. A Jeffrey's prior distribution was utilized for the proportion of infected nymphs. Credible intervals were constructed by calculating the highest density interval for the posterior distribution of the nymphal infection prevalence for each site. Credible intervals were also generated separately for each year and at every site incorporating temporal dependence, but the results are not presented here as there was no significant change in the infection prevalence estimates for each of the pathogens from year to year. Posterior distributions for all infection parameters were estimated using Markov chain Monte Carlo methods implemented in R (R Team Core 2018) with the rjags (Plummer 2018) package. As in the DON analysis, batch means procedures and multiple starting values were used to confirm convergence of the Markov chains. After an initial burn-in period of 10,000 samples, 1,000,000 samples were generated from the joint posterior distribution. The size of the Markov chain was increased until there were at least 100,000 effective samples from the posterior of each parameter. Informative priors were also considered, but some of the results for the individual sites were sensitive to the more informative priors. The noninformative (Jeffrey’s) prior results are presented here since the results for the nymphal infection prevalence estimates for each of the pathogens did not change significantly as the dispersion of the priors was increased for this prior, which indicates that the nymphal infection prevalence results presented here are driven by the data rather than the prior distribution. When reporting prevalence data from other published studies (i.e., Barbour et al. 2009, Hersh et al. 2014, Keesing et al. 2014, Prusinsky et al. 2014, Hutchinson et al. 2015, Farone et al. 2018), 95CI were calculated with the EpiTools package from Ausvet (2019) using the Jeffrey's prior (Brown et al. 2001).

A permutation test was developed to test for dependence between the prevalence of different pathogens in the same nymph, as performed in Hersh et al. (2014). To test for dependence, infection outcomes (infected or not) were simulated every year at each site for each tested nymph assuming the prevalence of each pathogen was independent. The nymphal infection prevalence of each pathogen at each site every year was tallied from the observed data, and then infections were simulated at each site for every year. In this approach, for the simulated infections, being infected with one type of pathogen has no impact on the probability that a nymph would be infected with a different type of infection. Once the simulated infection status of each nymph was determined for each of the pathogens, the total number of nymphs co-infected with more than one type of pathogen were counted. The process was repeated 200,000 times, and the median number of coinfections in the simulations was calculated. Also, a 95% central interval of the number of coinfected ticks in the simulations was computed by calculating the 2.5th and 97.5th percentiles from the number of co-infected nymphs in the simulations. Finally, the simulations were compared with the observed number of co-infected ticks to calculate a permutation test P-value.

Fig. 2.

DON at individual field sites. DON is estimated as the average number I. scapularis nymphs/100 m² in drags that were performed at each site, each year. Abbreviations and locations of sites are shown in Fig. 1. Error bars are represented as SEM.

Tick Collection and DON

Only three species of ticks were collected by dragging in closed-canopy forests, with I. scapularis being the most common (Table 2). Not a single specimen of any stage of H. longicornis or A. americanum was observed in any drag or collected among the 6,398 ticks reported in this study. I. scapularis DON at each site each year is presented as nymphs/100 m² per drag (Fig. 2). By analyzing GPX files with GPXSee v5.17 software, DON was estimated at each site as a function of distance dragged and time while moving. There was a strong linear association (R² = 0.950) between nymphs/h per drag and nymphs/100 m² per drag. The DON estimates for nymphs/30 min drag correlated with nymphs/100 m² per drag with an R² value of 0.737. Thus, in cases where distance measurements were not available for particular drags (e.g., when a phone was shared between team members), the average distance of each drag recorded at the same site on the same day was used to estimate nymphal density. Modeling was performed with our nymphs/100 m² data as recommended by (CDC 2019b) guidelines. Among almost all sites, we observed consistent variations in DON as a function of year (Fig. 2). Thus, we first analyzed the impact of year on DON. The initial model presented here is a Bayesian negative binomial regression model with year as a categorical predictor. The 95CI for the overdispersion parameter was estimated to be (1.80, 12.67), so there is overdispersion in the nymphal count data, and a Poisson regression model would not be appropriate. The resulting Bayesian credible intervals for the expected nymphs/100 m² each year are presented in (Fig. 3). Since there is no overlap in the 95CI, the years 2016 and 2018 have significantly lower DON than the years 2015 and 2017. With respect to the potential role of weather in understanding the differences in DON, 95CI for the coefficients of each of the weather covariates included zero. Thus, none of the weather variables considered in the regression model were statistically significant predictors of DON, nor was the day of the year.

Table 2.

Total ticks collected by dragging at 11 forested sites in the Lehigh Valley region

Pathogen Testing

Summary results of pathogen testing with samples prepared from I. scapularis nymphs collected in 2015–2017 and expressed as nymphal infection prevalence are shown in Table 3. Of the 1,721 samples that were screened with a Borrelia 23S qPCR assay (Courtney et al. 2004), 435 were positive (Av. Ct = 33.1 SD 2.2). All 23S qPCR-positive samples were tested with a B. burgdorferispecific OspA qPCR assay (Dibernardo et al. 2014). Of these 435 samples, 423 were OspA qPCR-positive (Av. Ct = 34.8 SD 2.2). Of the 12 23S-positive and OspA-negative qPCR samples, six were B. miyamotoi-specific GlpQ qPCR-positive (Av. Ct = 33.5 SD 1.9), and six were glpQ qPCR-negative. To identify these remaining six GlpQ qPCR–negative samples, we developed a duplex Borrelia p66 HRM assay. Four of the six samples were identified as B. burgdorferi by p66 HRM with Tm values between 77.2 and 77.8°C. One of the six samples had a Tm of 76.1°C which was out of our observed p66 HRM Tm range for either Borrelia species, and another sample could not be amplified. Thus, two of our 435 Borrelia 23S qPCR-positive samples remain unidentified.

Fig. 3.

Negative Binomial Regression model for DON based on year. The response variable is nymphs/100 m2 per day, as estimated from time while moving from individual GPX files using GPXSee version 5.17. Each dot on the graph represents the expected nymphs/100 m2 collected using the negative binomial regression model in a particular year. The lower and upper bounds represent the 95% credible interval.

Of the 1,721 samples that were screened with the A. phagocytophilum Msp2 qPCR assay, 60 were positive (Av. Ct = 27.9 SD 2.3). This qPCR assay, developed by Courtney et al. (2004), does not differentiate between the ‘human active’ (Ap-ha) strain and the strain that has not been shown to infect humans, known as Ap-variant 1 (Ap-v1). To differentiate between Anapalasma strain variants, Keesing et al. (2014) developed a 16S HRM assay that was reported to have nonoverlapping Tm values for each strain following the amplification of a 58 bp fragment containing the same divergent 16S sequence described by Massung et al. (1998). Using their 16S HRM method, we obtained indistinguishable Tm values with both Anaplasma strains. We could not resolve this disparity, but note that the Ap-ha and Ap-v1 strain reference sequences (Massung et al. 1998) differ reciprocally at only two positions (27A>G followed by 37G>A) in the 58 bp fragment used for their 16S HRM analysis. We were able to amplify and sequence the divergent 16S sequence (Massung et al. 1998) in 59 of our 60 Msp2 qPCR-positive samples, finding that 45 samples were Ap-V1 and 14 samples were Ap-ha (Table 3). One of our Msp2 qPCR-positive samples could not be amplified and remains unidentified.

Table 3.

Overall nymphal infection prevalence at 11 sites in the Lehigh Valley region of eastern Pennsylvania (2015–2017)

Of the 1,721 nymphal samples that were screened with the B. microti 18S qPCR assay (Teal et al. 2012), 62 samples had typical sigmoid-shaped amplification curves. Edwards et al. (2015) found that this 18S qPCR assay (Teal et al. 2012) rarely, but unexpectedly yielded a positive reaction with B. odocoilei, which infects deer and is not known to cause human disease. B. odocoilei is likely to be more abundant than B. microti in some regions (Armstrong et al. 1998), and a failure to distinguish between them can lead to overestimates of B. microti infection. Accordingly, each Babesia 18S qPCR-positive sample was re-tested by 18S DNA sequencing or with an HRM assay that we developed for a divergent region within the respective 18S sequences. Of these, 17 samples were identified as B. microti by 18S DNA sequencing and 31 samples were identified by 18S HRM (Total B. microti = 48, Av Ct = 36.8 SD 2.6). The average Tm of the 18S HRM assay was for these 31 B. microti samples was 79.3°C (SD 0.1). Three Babesia 18S qPCR-positive samples were identified as B. odocoilei by DNA sequencing and two by our 18S HRM assay, with Tm values of 81.6 and 81.7°C (total B. odocoilei = 5). Six Babesia 18S qPCR-positive samples with late amplification curves (Ct ∼ 42) had an average 18S HRM value of 78.5°C (SD 0.2), which diverged from the expected Tm range for either Babesia species. Three other late-amplifying 18S qPCR-positive samples (Av Ct = 41.2) could not be amplified by HRM or conventional PCR. Thus a total of 9 of 62 samples with typical sigmoid-shaped 18S qPCR-curves remain unidentified.

Nymphal Infection Prevalence and Coinfections

Overall (i.e., all sites combined) annual and 3-yr nymphal infection prevalence values for all pathogens are in Table 3. All of the 11 sites of this study had B. burgdorferi-infected nymphs in each of the 3 yr (2015–2017) that samples were screened for pathogens. Individual sites were generally consistent in their annual B. burgdorferi prevalence, so expected prevalence values were aggregated and 95CI were generated (Fig. 4). Although 2017 had the highest overall B. burgdorferi prevalence (Table 3), there was no evidence of significant variation in overall B. burgdorferi prevalence among the 3 yr, as the 95CI for each year overlapped. There was significant variation in the 3-yr prevalence among the different sites (Fig. 4). For example, SM was the site with the highest prevalence (36.7%; 95CI 29.5–44.4), and SP had the lowest prevalence (11.9%; 95CI 7.6–17.7). The sample size in this study was too small to assess the site-by-site variation in other pathogens besides B. burgdorferi. However, two sites (LM and SM) accounted for 38 of the 48 B. microti positive samples, and no B. microti-positive samples were found at six of the 11 sites. The 14 A. phagocytophilum variant-ha positive samples were distributed among five sites.

The results for the permutation test for association between infection types are presented in Table 4. The results indicate that a nymph is more likely to be co-infected with both B. burgdorferi and B. microti than would be expected if the infections occurred independently (P = 0.002). Similarly, a nymph is more likely to be infected with both B. burgdorferi and A. phagocytophilum variant-ha (human-active) (Ap-ha) than expected if these pathogens occurred independently (P < 0.001). Conversely, there is not a significant positive interaction between B. burgdorferi and A. phagocytophilum strain variant 1 (Ap-V1), which is not known to infect humans (P = 0.999).

Tick Species

Recent changes in the spatial distribution of ticks and their associated pathogens have been observed at regional and global scales (Sonenshine 2018). The primary goal of this study was to estimate the density of I. scapularis nymphs and to survey a subset of these ticks for the most common nonviral human pathogens for which they are known vectors (Nelder et al. 2016). However, our active surveillance also provided some information on the presence of other ticks that could be collected by dragging in forested habitats preferred by I. scapularis. With the caveat that our collection strategy was optimized for I. scapularis, these baseline data may be useful in documenting future changes in the distribution of other species of ticks. Very few American dog ticks were collected in the drags included in this study, but our sampling method in forested areas minimized the likelihood of collecting this tick which prefers the old-field-forest ecotone (Sonenshine and Levy 1972). D. variabilis were highly abundant in the grassy areas immediately adjacent to a forested site (Trexler Preserve) that we surveyed. The invasive and generally parthenogenetic longhorned tick was reported in adjacent Hunterdon County NJ and was found to overwinter at the site (Rainey et al. 2018). It was not collected by dragging in the four summer seasons of this study, suggesting that the veterinary pest and potential vector of human disease has not yet become widely established in forested regions of the Lehigh Valley region. In places where H. longicornis has invaded, it prefers pasture to forested habitats (Heath 2016), but it has been recently collected along with I. scapularis nymphs by dragging in forested habitats of the Hudson Valley region (Dr. Matthew Frye, NYS IPM, Cornell University, personal communication).

Fig. 4.

B. burgdorferi nymphal infection prevalence at individual sites. Each dot on the graph represents the expected infection prevalence at each site aggregated across the years 2015–2017. The lower and upper bounds represent the 95CI. Abbreviations and locations of sites are shown in Fig. 1.

Table 4.

Permutation test for nymphal coinfection

Lone star ticks are expanding their established range (Stafford et al. 2018) while emerging as a potential vector of several viral and bacterial pathogens (Springer et al. 2014), but not B. burgdorferi (Stromdahl et al. 2018). This tick is also associated with an anaphylactic reaction to certain kinds of meat (Commins et al. 2011). The habitat usage of A. americanum often overlaps with that of I. scapularis, and the two species can be collected by dragging in the same area (Johnson et al. 2017). Maps of the range of this tick often include all of PA (e.g., CDC 2019c). However, Springer et al. (2014) reported very few historical records of the presence of this tick in PA, showing one positive record from the 2010s. Snetsinger (1968) reported that lone star ticks had been collected 14 times from six scattered counties (including nearby Lancaster County), but not in the Lehigh Valley region. A single adult and nymph were reported at Gettysburg National Military Park (Adams County) in 2009 (Han et al. 2014). Johnson et al. (2017) indicate that the tick had not become abundant by 2014 or 2015, as none were detected in 7,500 m² of dragging in an area where 114 I. scapularis nymphs were collected. Springer et al. (2014) noted that the low number of historical records of the lone star tick in PA might be due to unfavorable conditions or a lack of sampling. Thus, as the southeastern species has a range that is generally moving northward (Stafford et al. 2018), continued surveillance for this tick should remain a high priority for all regions of the state.

Density of Nymphs

Our surveillance data provides a snapshot of DON in selected areas of the Lehigh Valley region in four successive years (Fig. 2). Our overall DON estimates with 95CI are shown in Fig. 3. Marked inter-annual fluctuations in DON at the same site have been observed previously (e.g., Hazler et al. 1996), emphasizing the importance of multi-year DON estimates. DON estimates at a given site are influenced by the true population density and the proportion of nymphs that can be captured by dragging. Weather variables that impact tick mortality influence inter-annual variations in population density (Eisen et al. 2016). Temperature and humidity also impact questing behavior (Vail and Smith 2002), which is an important consideration given that the proportion of nymphs of northern origin (such as PA) that quest above the leaf litter at any given time during the peak period of nymphal activity is <10% (Arsnoe et al. 2015). Between temperature extremes, tick survival is positively correlated with humidity, as reviewed by Eisen et al. (2016). Nymphs that do not have access to moisture become desiccated upon exposure to RH <82% (Rodgers et al. 2007). Burtis et al. (2016) showed a negative correlation between DON and hot and dry summer days, which may be explained by nymphal behavior as well as mortality. To assess the impact of weather on our inter-annual variation in DON, it was intuitively appealing to perform a series of Bayesian negative binomial regression models with meteorological data from the Lehigh Valley International Airport weather station (KABE). Using relevant weather parameters from Berger et al. (2014a) and Ostfeld et al. (2018), we did not detect an effect of RH, precipitation, accumulated degree days, or temperature on the striking inter-annual variations of DON (Fig. 3) in our Bayesian negative binomial regression model.

Berger et al. (2014b) noted that coarse-scale weather observations do not capture the daily moisture limiting dynamics that are biologically relevant to tick survival. They also documented a poor linear correlation between RH data recorded at nearest airport weather stations and field data collected within the relevant leaf litter habitat. Even when environmental conditions are meticulously documented in natural settings at the microhabitat scale, the correlation between RH and DON can be difficult to detect. For example, (Berger et al. 2014a) performed multiple observations of the microclimate at the forest-floor level in several repeatedly-sampled forested study plots. They found that RH recorded with 8–14 d before the sampling date (RH14) positively impacts the density of questing nymphs, but an association was not observed between DON and the RH at the time of sampling or even between DON and the average daily RH in the week before sampling.

Table 5.

DON and climate variables

In a 25-yr Connecticut (CT) study, Hayes et al. (2015) found a strong positive correlation between the January Standardized Precipitation Index of the same year and summer DON. January precipitation in CT is mostly in the form of snow, which may protect overwintering ticks from cold temperature extremes. Selected winter temperature and precipitation variables along with the annual DON estimates are presented in Table 5. With only four observation years, we cannot recreate the analysis of Hayes et al. (2015). One pattern that emerges from Table 5 is that the years with the highest DON (2015 and 2017) occurred after winters with the least January precipitation. Our analysis does not suggest that site- and microhabitat-specific weather parameters are unimportant components of the complex set of factors that influence DON in a given season. Rather, it questions the feasibility of predicting summer nymphal activity in the Lehigh Valley region based on low-resolution weather data.

Another factor that could potentially influence our observed annual fluctuations in DON is the regular 2-yr cohort cycle of I. scapularis nymphs in the relatively temperate parts of its range that include PA (Yuval and Spielman 1990). Accordingly, nymphs of the cohort that fed in early summer 2015 could potentially lead to a large cohort of the next generation of nymphs in 2017. Likewise, a relatively small cohort of questing nymphs in early summer 2016 could potentially have contributed to a small cohort of questing nymphs in summer 2018. However, an analysis of 19-years of data from six large field plots in the Hudson Valley region (Ostfeld et al. 2018) found that the density of larval ticks in the previous year had no effect on DON estimates, indicating that ‘demographic forcing’ from larval ticks to questing nymphs is weak.

Vertebrate host and predator diversity and abundance are important drivers of inter-annual variation in Lyme disease risk (Ostfeld et al. 2018). Our lack of site-specific data on hosts or predators precludes any speculation on the biotic factors that influenced the striking inter-annual fluctuations in DON that we observed at the same sites. Within the large ‘envelope’ of climate conditions that allow for tick survival (i.e., lack of mortality from freezing, overheating, dehydration, or drowning), meteorological factors are generally weak predictors of nymphal abundance, compared with other factors such as host availability (Eisen et al. 2016). Regardless of the underlying cause, the importance of observing these fluctuations in DON over 4 yr is that it will allow future studies to employ these data when evaluating inter-annual trends in reported cases of human disease in the Lehigh Valley region.

Borrelia burgdorferi

Our 3-yr overall nymphal infection prevalence of B. burgdorferi was 24.8% (95CI 22.8–26.9) (Table 3). Since the 95CI of each year overlap, there does not appear to have been a significant annual fluctuation in nymphal infection prevalence. Our B. burgdorferi prevalence is similar to the 18.3% (95CI 11.6–26.9) detected in 109 nymphs collected in the Lehigh Valley region in 2013 and 2014 (Edwards et al. 2015). It is also comparable to Hutchinson et al. (2015), who reported 12/24 B. burgdorferi-positive I. scapularis adults in Lehigh County (50%; 95CI 31–69) in 2013. Adult infection prevalence is generally about twofold higher than nymphal infection prevalence as a result of two potentially infective blood meals (Barbour et al. 2009). Our B. burgdorferi prevalence is higher than a study of eight counties in the Hudson Valley region (Prusinksy et al. 2014), reporting a prevalence of 14.4% (95CI 13.2–15.6) in 3,300 nymphs collected in 2003–2006. However, our prevalence is only slightly lower than to a more recent study in Dutchess County NY, also in the Hudson Valley region (Hersh et al. 2014), reporting a prevalence of 29.1% (95CI 27.8–30.5) in 4,368 nymphs collected in 2011–2012. Our study shows a B. burgdorferi nymphal infection prevalence that exceeds 20% for three consecutive years. The local prevalence of ≥20% B. burgdorferi-infected ticks is often cited as a threshold for clinical decisions regarding the administration of prophylactic antibiotics to adults following a tick bite with evidence of at least 36 h of attachment (Wormser et al. 2006). This study reinforces the public health message that personal tick prevention measures are necessary when at work or play in forested areas of the Lehigh Valley.

The local density of questing B. burgdorferi-infected I. scapularis nymphs (DIN) is one factor that influences the risk of contracting Lyme disease. Our DIN estimates for the Lehigh Valley region for 2015–2017 [(overall annual nymphal infection prevalence) × (overall annual DON)] are 1.04, 0.37, and 1.71, respectively. Since overall nymphal infection prevalence was similar each year (overlapping 95CI), then the 2018 DIN was approximated using the overall 3-yr average nymphal infection prevalence (Table 3) and 2018 DON (Fig. 2) as 0.25. Some studies have shown a positive correlation between DIN (also reported as Entomological Risk Index or ERI) and the local Lyme disease burden (Mather et al. 1996, Stafford et al. 1998) while Prusinski et al. (2014) did not find a general correlation. Eisen and Eisen (2016) critically review the concept of ERI and suggest that, within limitations, ERI is most informative as a predictor of infection risk at the community or county scales. Another confounding factor in understanding the relevance of DIN to disease transmission is that different strains of B. burgdorferi vary in their propensity to lead to a local or disseminated infection (Dykhuizen et al. 2008). Lee et al. (2014) also note that the presence of infected ticks in an area does not necessarily imply that human exposure is occurring. In the case of our study, all of the sites that were visited are near a densely populated metropolitan area with ∼820,000 residents and are permeated with hiking trails. We noted wide annual fluctuations in DIN in the Lehigh Valley region as DON differed significantly (Figs. 2 and 3), but B. burgdorferi infection prevalence was roughly constant (Table 3). Data are not yet available to determine whether these fluctuations in DON and DIN are reflected in the number of local cases of tick-borne disease.

Borrelia miyamotoi

Consistent with Edwards et al. (2015), we detected a 3-yr overall nymphal infection prevalence of B. miyamotoi (0.3%; 95CI 0.1–0.7) (Table 3). This prevalence is not significantly different from Barbour et al. (2009) who detected the pathogen in a single sample among 32 nymphs collected in nearby Reading, PA (3.1%; 95CI 0.3–13.7) and one in 100 nymphs collected near Harrisburg, in central PA (1.0%; 95CI 0.1–4.6). Most recently, Farone et al. (2018) estimated a B. miyamotoi infection prevalence of 1.7% (95CI 0.6–3.7) in 304 pooled I. scapularis adults that were collected from deer hunted in a sector of northeastern PA that includes the Lehigh Valley region. Their state-wide estimate based on 1,990 ticks, was ∼0.8% (95CI 0.4–1.2). Our B. miyamotoi nymphal prevalence may be lower than Farrone et al. (2018) because their study only included adult ticks which would have taken two blood meals. B. miyamotoi infected nymphs may be rare in the Lehigh Valley region, but human–tick encounters are common. Our surveillance suggests that some locally transmitted cases of tick-borne relapsing fever should be expected and as Wormser et al. (2019) suggests, the number of cases may be greater than reported.

Anaplasma phagocytophilum

Our 3-yr overall combined nymphal infection prevalence of A. phagocytophilum (3.4%; 95CI 2.4–3.9) is consistent with the 1.9% (95CI 0.8–3.7) adult and nymphal infection prevalence reported for 423 ticks collected in the Lehigh Valley region by Edwards et al. (2015). Hutchinson et al. (2015) did not detect A. phagocytophilum in 23 adult I. scapularis from Lehigh County but found a 4.0% (95CI 2.0–7.7) prevalence in 227 adults collected in the northeastern region of PA that includes our study area. Prusinsky et al. (2014) found a significantly higher total infection prevalence of 6.5% (95CI 5.7–7.3) in a sample of 3300 nymphs collected in eight Hudson Valley region counties, in 2003–2006.

Trost et al. (2018) identified three distinct clades of this intracellular bacterium by sequence divergence in the ankyrin A gene. The two clades that infect humans are known collectively as the Ap-ha strain and the strain that has not been shown to infect humans is known as Ap-v1. The Ap-v1 strain infects deer (Massung et al. 2005) and several other animals, but rarely P. leucopus (Keesing et al. 2014). Prior to Trost et al. (2018) and in our study, the strains were differentiated by the divergence of two nucleotides in a PCR-amplified 16S rDNA gene fragment (Massung et al. 1998). We found that the Ap-v1 strain was threefold more prevalent (nonoverlapping 95CI) than the Ap-ha strain in I. scapularis nymphs (Table 3). In the Hudson Valley region of NY (Dutchess County), the inverse ratio was detected between the two strains. Keesing et al. (2014) tested 5098 I. scapularis nymphs collected from host animals in 2011–2012 with a novel 16S rDNA HRM assay and documented an infection prevalence of 2.7% (95CI 2.3–3.2) with Ap-ha, and 0.8% (95CI 0.6–1.1) with Ap-v1. They also detected a significantly higher prevalence of the Ap-ha variant in questing nymphs. Hersh et al. (2014) did not observe a significant difference between the expected and observed coinfection frequencies of A. phagocytophilum (variants not distinguished) and B. burgdorferi in 4,368 questing I. scapularis nymphs collected in Dutchess County, NY. We observed a significant difference between the observed and expected coinfection frequencies of the Ap-ha strain and B. burgdorferi in I. scapularis nymphs (Table 4), but found no positive correlation between the Ap-v1 strain and B. burgdorferi.

Our results suggest that there may be differences in the dynamics of A. phagocytophilum between the Lehigh Valley and Hudson Valley regions that warrant further study. From 2008 to 2012, PA had a relatively low human granulocytic anaplasmosis incidence of 0.3 cases per million, compared with neighboring NJ and NY with 10.4 and 23.3 cases per million, respectively (Dahlgren et al. 2015). Surveillance that differentiates between the two A. phagocytophilum variants should continue in the Lehigh Valley region, as a westward spread of the Ap-ha strain into PA could have important clinical consequences.

Babesia microti

Our 3-yr overall nymphal infection prevalence of B. microti (2.8%; 95CI 2.1–3.7) indicates an increased prevalence from our previous study of 423 adult and nymphal ticks (<0.5%; 95CI 0–1.7) (Edwards et al. 2015). The B. microti nymphal infection prevalence increased each year (Table 3) but the result was not statistically significant due to overlapping 95CI. Our B. microti prevalence is similar to a regional study of eight counties in the Hudson Valley region (Prusinksy et al. 2014), reporting a prevalence of 2.7% (95CI 2.2–3.3) in 3,300 nymphs collected in 2003–2006. However, our prevalence is significantly lower than a more recent study in Dutchess County NY, also in the Hudson Valley region (Hersh et al. 2014), reporting a prevalence of 13.7% (95CI 12.7–14.8) in 4,368 nymphs collected in 2011–2012.

As with human granulocytic anaplasmosis, PA has a lower burden of babesiosis than might be expected given its high incidence of Lyme disease (CDC 2019a). Around 94% of the cases of babesiosis in a recent CDC report (CDC 2016) are from seven northeastern and mid-western states that do not include PA. The clinical significance of the positive correlation that we observed between B. microti and B. burgdorferi (Table 4) is that coinfection with both pathogens increases the severity and duration of both diseases (Krause et al. 1996). Previous studies suggest that coinfection enhances the emergence of B. microti in places where B. burgdorferi is already established since B. microti might otherwise have a low level of ecological fitness (Dunn et al. 2014). In their meta-analysis of 21 studies (Diuk-Wasser et al. 2016), more than half showed a greater probability of coinfection with B. burgdorferi and B. microti than predicted by independent infection, whereas none of the studies reported a lower probability of coinfection. Our study contributes additional evidence of a positive association between B. burgdorferi and B. microti in a region where B. burgdorferi-infected nymphs are common. As with the spread of Lyme disease into PA around thirty years ago, it is possible that a westward range expansion of ticks infected with B. microti into the state could manifest itself at an early stage in the Lehigh Valley region.