In California, the western blacklegged tick, Ixodes pacificus Cooley and Kohls, is the principal vector of the Borrelia burgdorferi sensu lato (sl) complex (Spirochaetales: Spirochaetaceae, Johnson et al.), which includes the causative agent of Lyme disease (B. burgdorferi sensu stricto). Ixodes pacificus nymphs were sampled from 2015 to 2017 at one Sierra Nevada foothill site to evaluate our efficiency in collecting this life stage, characterize nymphal seasonality, and identify environmental factors affecting their abundance and infection with B. burgdorferi sl.To assess sampling success, we compared the density and prevalence of I. pacificus nymphs flagged from four questing substrates (logs, rocks, tree trunks, leaf litter). Habitat characteristics (e.g., canopy cover, tree species) were recorded for each sample, and temperature and relative humidity were measured hourly at one location. Generalized linear mixed models were used to assess environmental factors associated with I. pacificus abundance and B. burgdorferi sl infection. In total, 2,033 substrates were sampled, resulting in the collection of 742 I. pacificus nymphs. Seasonal abundance of nymphs was bimodal with peak activity occurring from late March through April and a secondary peak in June. Substrate type, collection year, month, and canopy cover were all significant predictors of nymphal density and prevalence. Logs, rocks, and tree trunks had significantly greater nymphal densities and prevalences than leaf litter. Cumulative annual vapor pressure deficit was the only significant climatic predictor of overall nymphal I. pacificus density and prevalence. No associations were observed between the presence of B. burgdorferi sl in nymphs and environmental variables.
The western blacklegged tick, Ixodes pacificus, Cooley and Kohls, is the primary vector of the Lyme disease bacterium Borrelia burgdorferi sensu stricto (Spirochaetales: Spirochaetaceae, Johnson et al.) in California, and nymphs of this species are responsible for most disease transmission to humans (Clover and Lane 1995). The highest incidence of Lyme disease in western North America occurs in northwestern California; however, areas of increased incidence occur in other regions of the state, including the northern Sierra Nevada foothills east of the Sacramento Valley (Eisen et al. 2006, CDPH 2018).
Despite the broad distribution of human cases of Lyme disease and I. pacificus in California, most abundance and phenology studies of I. pacificus nymphs have focused on northwestern California (Clover and Lane 1995; Eisen et al. 2001, 2002, 2003; Lane et al. 2007; Salkeld et al. 2014). Variation in seasonal activity occurs within this region of the state. For example, in climatically different areas of Mendocino County, peak numbers of nymphs were reached about a month earlier and decreased more sharply in an oak/madrone woodland in comparison to a redwood/tanoak woodland (Eisen et al. 2002). In other coastal counties near Mendocino, Salkeld et al. (2014) found more protracted nymphal activity, beginning in early February in Marin and San Mateo counties, and extending into October in Sonoma County. In contrast, MacDonald and Briggs (2016) observed truncated nymphal activity in central and southern California locations, demonstrating the need to evaluate differences in I. pacificus phenology elsewhere in this ecologically diverse state.
Other Ixodes spp. studies have demonstrated significant effects of temperature, humidity, and saturation vapor pressure on nymph abundance and activity (Vail and Smith 1998, Perret et al. 2000, Diuk-Wasser et al. 2010, Berger et al. 2014). In northwestern California, Eisen et al. (2003) found that I. pacificus nymph density began to decline when mean maximum daily temperatures exceeded 21–23°C and mean maximum daily relative humidity decreased below 83–85%.
Habitat and climatic variables have also been associated with the density or prevalence of B. burgdorferi-infected I. pacificus nymphs in northwestern California (Eisen et al. 2010, Swei et al. 2011). At a microhabitat scale, higher prevalences of B. burgdorferi infection in I. pacificus nymphs were found on logs and tree trunks in dense woodland habitats (Lane et al. 2007). The effects of habitat and climate variables on nymphal I. pacificus activity or their prevalence of infection with B. burgdorferi sensu lato (sl) have not been extensively evaluated in the foothills of the northern Sierra Nevada.
Dragging vegetation over linear transects is an effective way to evaluate adult I. pacificus density and temporal activity (Eisen et al. 2018); however, nymphs of this species do not typically quest on vegetation (Lane et al. 1995, Talleklint-Eisen and Lane 2000). Most nymphs remain on undersides of leaf litter or below the soil surface and are usually sampled by dragging leaf litter (Lane et al. 1995). In Sierra Nevada foothill habitats where I. pacificus adults are often abundant, nymphs have proven difficult to collect by dragging leaf litter in comparison to northwestern California habitats. We have found dragging leaf litter rarely yielded greater than five nymphs per hour of effort in areas where adult I. pacificus were common (CDPH, unpublished data). Additionally, oak woodlands in the Sierra foothills often have a dense understory of shrubs, which minimizes open areas of leaf litter that are easily accessible for sampling. For these reasons, collecting I. pacificus nymphs from leaf litter in the Sierra foothills has been labor intensive and relatively unproductive, making it difficult to accurately determine nymphal distribution, phenology, and prevalence of infection with B. burgdorferi sl.
Lane et al. (2007) observed greater mean densities of I. pacificus nymphs on wood substrates (logs and tree trunks) than in adjoining leaf litter or grass, which suggested a potentially more efficient sampling strategy for this life stage. Greater numbers of nymphs may be collected from these substrates because their larval-stage hosts frequently utilize the same biotopes and drop replete larvae there; the nymphs largely remain on these substrates and utilize an ambush host-seeking strategy to find their next host (Lane et al. 2009).
To improve our nymphal surveillance efficiency and better understand habitat and climatic characteristics that affect nymphal I. pacificus abundance and seasonality in the Sierra Nevada foothills, we systematically sampled four questing substrates over a 3-yr period at one study site. We conducted this surveillance to address three primary objectives. The first objective was to compare our flagging success for I. pacificus nymphs from four questing substrates (leaf litter, tree trunks, rocks, and logs) and compare fine-scale habitat characteristics and temporal associations with nymph abundance. The second objective was to evaluate associations between seasonal nymphal activity and climatic variables. The third objective was to determine whether any associations exist between nymphal B. burgdorferi sl presence and substrate-level habitat characteristics. To accomplish the first objective, we modeled nymphal density (i.e., nymph count per meter flagged) and prevalence (presence or absence of a nymph) at the individual flagging sample level as a function of substrate type, month and year of sampling, and fine-scale habitat variables (i.e., substrate-level models). The second objective was addressed by modeling nymphal density and prevalence (summarized by collection date) as a function of temperature, humidity, precipitation, and vapor pressure deficit transformations (i.e., climate models). To address the third objective, we compared the presence/absence of B. burgdorferi sl in nymphal I. pacificus at the individual sample level with substrate type, month and year of sampling, and fine-scale habitat variables (i.e., infection models).
Materials and Methods
Study Site
The approximately 6.2-hectare (15.3 acre) study site (38.771324, –121.048042) was located near Salmon Falls in the Folsom Lake State Recreation Area in western El Dorado County and ranged in elevation from 150 to 200 m above sea level. The site was selected because of its proximity to a hiking trail where I. pacificus adults are common. From 2015 to 2017, monthly drag sampling along this trail found that adult tick abundance peaked in March at an average of 26 I. pacificus per 100 m sampled (CDPH, unpublished data).
Salmon Falls receives an average of 61 cm of rain annually, with most falling between the months of October and May (Western Regional Climate Center— https://wrcc.dri.edu). The two primary habitat types within our study site were 1) mixed oak woodland, dominated by the tree species Quercus douglasii, Q. kelloggii, Q. wislizeni, and Pinus sabiniana, with an understory of toyon (Heteromeles arbutifolia) and poison oak (Toxicodendron diversilobum) and 2) mixed chaparral, dominated by the shrubs Arctostaphylos spp. and Ceanothus spp. To account for habitat variation, the study site was divided into three sampling zones, corresponding to a gradient of habitat types: zone one was dominated by dense oak woodland, zone three was dominated by chaparral—blue oak woodland, and zone two was an intergrade of zones one and three.
Data Collection
Sampling began in February 2015. Because I. pacificus nymphs were collected during the initial sampling event (17 February 2015), we started sampling in late January in 2016 and 2017. Sampling was conducted approximately every 2 wk thereafter and ended in July each year, when no I. pacificus nymphs were collected. Sampling time of day ranged from 1100 to 1630. Host-seeking nymphs were sampled with 1 m2 white flannel cloths by flagging four substrate types: leaf litter, logs, tree trunks, and rocks. The flannel cloths (i.e., flags) were brushed over or wrapped around accessible surfaces of tree trunks, logs, and rocks (tree trunks were sampled ≤2 m from ground). Leaf litter surface was lightly disturbed (i.e., turned over) with the flag to contact underlying surfaces. The flags were in contact with substrate surfaces for 20–30 s to maximize contact with nymphs. The sampled areas of leaf litter and rocks were estimated by coverage of the flags. The sampled areas of tree trunks and logs were estimated by measuring the height (tree trunks) or length (logs) covered by the flags, recording the average diameter of the sampled substrate, and later calculating the total surface area sampled (described in further detail for log samples below). The estimated area sampled for each substrate ranged from 0.5 to 3 m2 (µ = 1 m2) due to the variable size and accessible areas of selected substrates. After each substrate was sampled, any observed nymphs were removed from the flags and placed in collection vials labeled for individual substrates. Flags were inspected after all nymphs were removed and again before sampling another substrate. In total, six biologists conducted the substrate sampling (range = 2–5 individuals per event). Individuals were randomly assigned to sampling zones for each event and substrates sampled within zones were randomly chosen by the individual collectors. Additional measurements recorded included estimate of canopy cover (low, medium, or high); genus or species (if known) of tree, log, and leaf litter. Canopy cover categorization was based on the estimated proportion of canopy coverage directly above the sampled substrate. Low canopy corresponded with less than 1/3 cover, medium was from 1/3 to 2/3 cover, and high canopy was above 2/3 coverage. Meteorological data (hourly temperature and relative humidity) were recorded at the study site using an EasyLog USB data logger (Lascar Electronics Inc., Erie, PA) located approximately 200 m east of the sampling zones and at a similar elevation and distance from the lake. Temperature and relative humidity data were collected at a single location, not within each zone or at each substrate sampled, to compare with summarized nymphal density and prevalence estimates (averaged for all substrates by collection date). Precipitation data from 2014 to 2017 for Folsom Lake was downloaded from the Western Regional Climate Center ( https://wrcc.dri.edu).
Laboratory Methods
Collected nymphs were identified to species and tested individually at the California Department of Public Health, Vector-Borne Disease Section laboratory for the presence of B. burgdorferi sl—the bacterial complex that includes B. burgdorferi sensu stricto—and Borrelia miyamotoi as described previously (Padgett et al. 2014).
Data Analysis
Data Management and Calculations
To calculate nymph density for an individual sample, the total number of nymphs collected was divided by the total area sampled (e.g., 10 nymphs/2 m flagged = 5 nymphs/meter). The average density for a collection date and substrate type was then calculated for analysis of climate models. Because the lower portions of many smaller logs could not be completely sampled (i.e., the flag did not contact log surfaces adjacent to or touching the ground) an area correction factor was included for logs with diameters < 40 cm. For this correction, we used the calculation for the surface area of a cylinder to represent the area of a log and assumed that an average of 2/3 of a given log's surface was available for sampling. The area calculation was multiplied by 2/3 to account for the underside portion of the log that was unavailable for sampling. This correction was only used on logs with a diameter less than 40 cm because a 1 m2 cloth would typically cover the accessible surface of any log ≥ 40 cm in diameter.
The hourly temperature and relative humidity measurements from the data logger were used to calculate the average daily vapor pressure deficit (DVPD; Murray 1967) and cumulative daily vapor pressure deficit (CVPD; January 1—collection date). The formula for vapor pressure deficit which was used as a basis for all further summary calculations (e.g., DVPD, CVPD) is below.
where RH = relative humidity, SVP = saturation vapor pressure, T = temperature.
Substrate-Level Models
All statistical analyses were conducted using R statistical software (version 4.0.3, R Core Team 2020). To address our first objective, we used negative binomial mixed models (glmmTMB package version 1.0.2.1, Magnusson et al. 2020) to analyze a suite of substrate-level models comparing I. pacificus nymph density (counts with an offset of meters sampled) with substrate type, collection month, collection year, and canopy cover at the individual sample level. We also modeled the data using zero inflated negative binomial mixed models, but they performed poorly when compared with standard negative binomial models. We then used multiple logistic regression to model the prevalence of I. pacificus nymphs (presence or absence) as a function of the same variables used for nymph density. We included the random effect of collection date in all models to account for any issues arising from the repeated measures collection method. However, the variances for the random effect of collection date for the fully parameterized substrate-level models were very small for both density (σ2 = 0.02) and prevalence (σ2 = 3.9 × 10–8), indicating that the effect of collection date was not significant.
Substrate-level negative binomial mixed models of nymph density by species of log and tree (Quercus spp. vs Pinus sabiniana, if known) were conducted separately because of collinearity with the substrate variable in the main model. To assess the broad effect of habitat type and to compare relative density of nymphs among zones, negative binomial mixed models were used on a reduced data set to compare nymph density (counts per meter squared) by sampling zone, using only substrate samples with nymphal I. pacificus present.
Climate Models
To accomplish our second objective, we used negative binomial mixed models to analyze substrate-level count data summarized and averaged for each collection date. Models compared average tick density (for a collection date), collection year, total precipitation for the current year, total precipitation for the prior year, average daily temperature and relative humidity, maximum daily temperature and relative humidity, average and maximum temperature and relative humidity for seven and 14 d prior to collection, average and maximum monthly temperature and relative humidity, average DVPD, average vapor pressure deficit for seven and 14 d prior to collection, and CVPD. All climate models utilized an offset term to account for the total number of meters sampled for each collection date, which allowed analysis of nymph counts per unit effort (i.e., average nymph density). The same climate models analyzed for nymph density were also analyzed for nymph prevalence (i.e., proportion of substrates with at least one nymph collected) averaged by collection date using generalized linear mixed models with a Tweedie distribution. All models included the random effect of collection date; however, the variance associated with this term was always near zero.
Infection Models
Our third objective was accomplished using multiple logistic regression to compare the presence or absence of B. burgdorferi sl in individual nymphs with all the variables described for the substrate-count models, less the meteorological calculations. The random effect of collection date was included for all models.
Akaike's information criterion (AIC) was used for all analyses to select the best-fit model for those models using the same data set. All model coefficients were exponentiated to obtain odds ratios (OR; for logistic regressions) and rate ratios (for count models). Odds ratios and rate ratios were used to compare the variables included in a model with a referential category (e.g., rocks vs leaf litter). A ratio greater than one indicates a positive relationship and a ratio less than one indicates a negative relationship. Because of the nearly identical results of the substrate-level nymph density and prevalence models, graphical representations of the data are only presented for nymphal density. However, model results and general summaries for prevalence are provided in the text and tables.
ArcGIS Pro (v 2.5) was used to map the locations and analyze potential spatial patterns of all nymphs collected throughout the study. The kernel density tool was used to create density maps to compare nymph density among sampling zones.
Results
In total, 2,033 substrates were sampled during 33 surveillance events from 17 February 2015 to 18 July 2017 (Table 1). This total consisted of 519 samples from logs, 506 from rocks, 510 from tree trunks, and 498 from leaf litter. The number of surveillance events and total substrates sampled varied by year, with 9 events (n = 605 substrates) in 2015, 11 events (n = 660 substrates) in 2016, and 13 events (n = 768 substrates) in 2017. The average number of substrates sampled per collection date was 15.6 (range = 5–24) for logs, 15.3 (range = 6–25) for rocks, 15.4 (range = 6–20) for tree trunks, and 15.0 (range = 4–19) for leaf litter. The lowest sample size (total of 21 substrate samples, range 4–6 per substrate) was in January 2016 prior to the first seasonal detection of nymph activity.
In total, 742 I. pacificus nymphs were collected, accounting for 95% of all Ixodid nymphs found. Other nymphs that were collected but not included in further analyses were 25 I. spinipalpis, 3 I. brunneus, 6 Dermacentor occidentalis, 3 Dermacentor spp., and 5 Haemaphysalis leporispalustris.
No I. pacificus nymphs were found during the initial sampling events in January 2016 and 2017. Nymphs were first collected during sampling events conducted from 11 to 17 February for each of the 3 yr of this study. Ixodes pacificus nymphs were found in all subsequent surveillance events until none were collected on 7 July and 6 July, in 2015 and 2016, respectively. Two I. pacificus nymphs were collected on 6 July 2017, but none were found on the next sampling date, 18 July 2017.
Table 1.
Nymphal Ixodes pacificus density, prevalence, and Borrelia burgdorferi sensu lato and Borrelia miyamotoi infection prevalences summarized by substrate type, canopy cover, and year at Salmon Falls, Folsom Lake State Recreation Area
Substrate-Level Analysis and Models
The top ranked substrate-level nymph density (AIC = 2,510.8, Table 3) and prevalence (AIC = 1,580.7) models included the variables collection year, collection month, substrate, and canopy cover. The average nymphal I. pacificus density for the study period was 0.31 nymphs per meter (n/m) sampled (Table 1). Average nymphal density was lowest (0.15 n/m) in 2015 (n = 113 nymphs collected), peaked in 2016 at 0.47 n/m (n = 343), and was 0.31 n/m in 2017 (n = 286). Significantly, more nymphs were collected per meter in both 2016 and 2017, compared with 2015 (Tables 1 and 2). No significant difference was observed in these comparisons when the collection dates with zero nymphs collected were dropped from the analysis. The overall prevalence of I. pacificus nymphs on substrates was 0.19 (Table 1). The between-year patterns of nymphal I. pacificus prevalence were the same as nymphal density, where 2016 and 2017 had significantly more substrates with I. pacificus present, compared with 2015 (Tables 1 and 2).
The density and prevalence of I. pacificus nymphs varied within years (P < 0.001, Tables 2–4). The initiation of I. pacificus nymph activity occurred in mid-February (range = 11–17 February) from 2015 through 2017. Collection dates for peak densities ranged from late March (25 March 2015) to late April (27 April 2017; Fig. 1). A secondary peak in nymphal density was observed in early June of 2016 and 2017. Seasonal peak density periods, defined as sampling dates ≥ 50% peak density, began in mid-late March (14–30 March) for all 3 yr and ended in late May (29 May 2015) or early June (1 June 2016 and 9 June 2017). Peak prevalence occurred in April (11–28 April) each year, with seasonal peak prevalence periods beginning with the first March (10–17 March) samples each year and, identical to peak density periods, ended in late May in 2015 and early June in 2016 and 2017.
Nymphal I. pacificus density varied depending on the substrate type sampled (P < 0.001, Tables 1–3). Compared with leaf litter (0.10 n/m), all other substrates yielded greater average densities of I. pacificus nymphs. Logs (0.62 n/m) and rocks (0.42 n/m) exhibited the greatest difference, while the increase on tree trunks (0.17 n/m) was relatively moderate. This pattern was observed for each year of the study (Fig. 1). Similarly, the prevalence of I. pacificus on logs (0.26), rocks (0.25), and, to a lesser extent, tree trunks (0.14) was greater when compared to leaf litter (0.11; P < 0.01, Tables 1, 2, and 4).
Table 2.
Summary of the coefficients (Coef), rate ratios (RR), odds ratios (OR), standard errors (SE), and P values for all predictor variables in the top ranked substrate-level nymphal models for Ixodes pacificus density and prevalence
Table 3.
All substrate-level Ixodes pacificus density models analyzed for this study
The abundance of I. pacificus nymphs varied significantly by estimated canopy cover (P < 0.01, Tables 1 and 2). Greater densities of nymphs were collected on substrates with medium and high canopy cover compared with substrates with low cover. Nymph prevalence mirrored density estimates and greater prevalences were observed with medium and high canopy cover, compared with substrates with low cover (Tables 1 and 2).
Overall nymphal I. pacificus density varied by sampling zone (Fig. 2). Compared with the chaparral-blue oak zone (1.20 n/m [143/120]), greater densities of nymphs were observed in the dense woodland zone (1.5 n/m [317/211], rate ratio [RR] = 2.00, P < 0.001) and the intergrade zone (1.88 n/m [282/150], RR = 1.80, P = 0.001).
Table 4.
All substrate-level Ixodes pacificus prevalence models analyzed for this study
Tree species (Quercus spp. vs P. sabiniana) also had a significant impact on the density and prevalence of nymphs collected from both tree trunk and log substrates. Pinus sabiniana samples had a greater density of nymphs (0.55 n/m [126/231], RR = 1.74) than Quercus spp. (0.27 n/m [253/922], RR = 0.57, P = 0.005). This pattern was observed for all years of the study. However, prevalence of I. pacificus nymphs exhibited the opposite pattern where the odds of sampling at least one I. pacificus nymph from a Quercus spp. substrate were greater than for a P. sabiniana substrate (OR = 1.61, P = 0.02).
Climate Analysis and Models
Cumulative vapor pressure deficit (CVPD) was the only significant climatic predictor of I. pacificus nymphal density (P = 0.018, RR = 0.99) and prevalence (P = 0.017, RR = 0.99) analyzed in this study. Hereafter, we only present results of nymphal density because of the similarity in results between density and prevalence of I. pacificus nymphs. The relationship between CVPD and average density and prevalence of nymphs was polynomial (CVPD2), not linear, where small and large CVPD values were associated with a lower density and prevalence of nymphs (Fig. 3). Overall, I. pacificus density peaked at an average CVPD of 21.6. In 2015, nymphal density peaked once on March 25 and began to decline at a CVPD of 21.3, whereas nymphal densities peaked twice in 2016 and 2017. The first peak in 2016 (15 April) and 2017 (27 April) was at a CVPD of 21.7 and the second peak was at a CVPD near 52 (1 June 2016 = 52.3, 9 June 2017 = 51.7). In 2016 and 2017, the second peak in nymphal density corresponded with elevated numbers of nymphs collected from rock substrates. No nymphs were collected when the CVPD was greater than 106.
Infection Models
The overall prevalence of B. burgdorferi sl in I. pacificus nymphs was 6.2% (46/742), whereas the prevalence of B. miyamotoi was 1.3% (10/742; Table 1). Borrelia burgdorferi sl prevalence was 11.5% in 2015 (13/113), 2.0% in 2016 (7/343), and 9.1% in 2017 (26/286). Borrelia miyamotoi prevalence was low, did not vary significantly among years, and was not analyzed further. The top ranked model (AIC = 329.9) that evaluated the presence or absence of B. burgdorferi sl in nymphal I. pacificus included the variable collection year. Compared with 2016, the odds of sampling a nymphal I. pacificus positive for B. burgdorferi sl was greater in both 2015 and 2017 (P < 0.001, Tables 1 and 5). None of 24 I. spinipalpis nymphs tested positive for either B. burgdorferi sl or B. miyamotoi. No significant associations were found between the presence of B. burgdorferi sl and substrate type, canopy cover, sampling zone, or tree species.
Discussion
Substrates
Sampling I. pacificus nymphs by flagging rocks, logs, and, to a lesser extent, tree trunks was consistently more productive than flagging leaf litter in this Sierra Nevada foothill location. In areas where I. pacificus nymphs are less abundant or more difficult to routinely collect by dragging leaf litter, sampling other nymphal questing substrates can be a more efficient way to obtain adequate numbers of nymphs to better estimate seasonality, infection prevalence, and Lyme disease risk.
In Mendocino County, Lane et al. (2007) similarly found higher I. pacificus nymph densities on logs and tree trunks than adjoining leaf litter or grass. The authors theorized these nymphs may have evolved a preference for questing on wood versus leaf litter or grass because their principal host, the western fence lizard (Sceloporus occidentalis), utilizes wood substrates for feeding and behavioral activities (basking, courtship), or nymphs may be easier to collect efficiently from the surfaces of logs and trunks. Another possibility, related to the former, is that a range of I. pacificus larval hosts (i.e., lizards, birds, rodents) frequent these substrates for harborage and food, and therefore a greater number of bloodfed larvae are dropped on or near these substrates (Lane et al. 2009). We did not attempt to discern reasons for significantly higher nymph collections from rocks, logs, and tree trunks, but it is likely that a combination of tick and host behaviors, along with the efficiency of collecting nymphs from the surfaces of rock and wood substrates, were responsible for the increased collection of I. pacificus nymphs. For example, although we more often observed the presence of at least one nymph on oak substrates, there may be a difference in the physical structure of pine versus oak bark that resulted in the collection a greater average number nymphs on individual pine substrates. Further study of these and other landscape attributes on nymphal density and prevalence are needed to better evaluate these observations.
Temporal Associations
The observed initial seasonal peak (late March–April) and duration (early February–early July) of I. pacificus nymphal activity at our study site generally falls between the seasonal range observed in northwest and southern California. However, an obvious secondary peak of nymphal activity (early June) was not observed by other studies (Table 6). Eisen et al. (2002) observed peak activity in early May with activity continuing into late July/mid-August in oak-madrone woodland habitat. A shorter duration of total and peak nymphal questing activity, including an earlier peak and faster decline, occurred in warmer, drier habitats (oak woodlands) compared with cooler coniferous habitats (Eisen et al. 2002). In a study that sampled additional sites in Mendocino County, nymphal activity lasted until at least late June and continued into October at three sites (Eisen et al. 2003). Similarly, nymphal activity in Marin, San Mateo, and Sonoma counties peaked from April into June (Salkeld et al. 2014). Nymphs were absent after August in Marin and San Mateo counties, but still present in early October in Sonoma County. In coastal southern California, I. pacificus nymphal activity was truncated in comparison to northwest California studies (MacDonald and Briggs 2016); first activity was detected in late February to early March and ceased by mid-April or mid-May in Santa Barbara County and by end of June in Los Angeles County.
A bimodal temporal distribution of both I. pacificus nymph density and prevalence was observed in all 3 yr at our study site but was much stronger in 2016 and 2017. The first peak in nymphal activity generally corresponded with nymphs collected on logs, whereas the second peak was primarily influenced by nymphs collected on rocks. Although this pattern was strongest in the latter 2 yr of the study, the relatively low numbers of I. pacificus nymphs collected in 2015 may have masked a stronger bimodal pattern of nymphal density and prevalence observed in 2016 and 2017. In 2016, the bimodal seasonal distribution was influenced primarily by a large number of nymphs collected off a few rock substrates (e.g., 35 nymphs collected from one rock in June). However, in 2017, we observed increased late-season density and a more widespread spatial distribution of nymphs on rocks. This, combined with a moderate increase in late-season nymph density on logs, resulted in increased overall nymph density in June 2017 (Fig. 1). If this bimodal activity was not a sampling anomaly, perhaps there is a behavior of vertebrate hosts or some differences in environmental or microhabitat conditions that cause later nymphal activity on rocks compared with other substrates.
Fig. 3.
Nymphal Ixodes pacificus density (nymphs per meter sampled) for all substrates in relation to cumulative vapor pressure deficit. Cumulative vapor pressure deficit is calculated as the sum of the daily average vapor pressure deficits from 1 January to the collection date. Note that y-axis scales vary by substrate type to better visualize the relationship between cumulative vapor pressure deficit and nymphal density.
Table 5.
Summary of the overall prevalence, coefficients (Coef), standard errors (SE), odds ratios (OR), and P values for all predictor variables in the top ranked logistic regression infection mixed model (via Akaike's information criterion) comparing the presence/absence of B. burgdorferi sl and year
Habitat Associations
Eisen et al. (2002) identified variation in nymphal activity based on broad, landscape-level differences in habitat type (e.g., oak-woodland vs conifer). Our study identified significant variation in nymphal density within different oak-woodland habitat types (i.e., sampling zones) in the Sierra foothills region that is also likely influenced by local landscape and climate variables, including differences in habitat type and vegetation composition. The sampling zones dominated by oak woodlands (zones one and two) generally had more canopy cover (e.g., medium to high) and a denser shrub understory, while the chaparral-blue oak habitat (zone three) was more open, drier, and with a less-dense understory. Substrates with more canopy cover may have lower temperatures and higher relative humidity, which may provide better conditions for nymph survival, especially in warmer months. Talleklint-Eisen and Lane (2000) found canopy cover to be highly correlated with the depth of leaf litter; however, their analysis focused on leaf litter depth rather than canopy cover, making it difficult to determine whether observed nymphal I. pacificus density was influenced by canopy cover, leaf litter depth, or a combination of both. The correlation between these two environmental variables suggests a strong link that is undoubtedly associated with fine-scale nymphal activity. Other factors such as host density and distribution, topography, and microclimate variables may also contribute to the variation in density of nymphs observed among habitat types, and these factors likely interact with those we identified to influence the variation in abundance of nymphs at our study site.
Table 6.
Nymphal Ixodes pacificus activity by month at the Salmon Falls, El Dorado County study site, 2015–2017, and as reported in other California studies
Climate Associations
The seasonal activity of nymphs at our study site was influenced by the interaction of local temperature and humidity (i.e., CVPD). Prior studies in coastal northern California showed that the duration and timing of peak I. pacificus questing activity was negatively associated with maximum daily temperature and positively associated with rainfall (Eisen et al. 2002). We observed no significant associations between these and other climate variables and nymphal density; however, CVPD was a significant predictor of overall nymphal density and prevalence in our study. Although the variation in minimum, maximum, and average monthly temperature and relative humidity has been shown to affect peak nymphal activity (Vail and Smith 1998, Eisen et al. 2003), this association can fluctuate greatly. The timing of the initial peak in nymphal density at our study site varied by over a month (33 d) between years. Importantly, regardless of the year or date, the initial peak in nymphal abundance was associated with a CVPD of 21, whereas the mean maximum temperature and relative humidity varied broadly between the annual peak periods. A secondary peak in nymphal activity, primarily associated with rock substrates in early June (1–9 June), corresponded with a CVPD threshold of 52 (prior 20 d mean maximum temperature ranged between 28.3 and 30.9°C and relative humidity was 87.2% each year). The CVPD had a more consistent association with peak nymphal I. pacificus activity than collection date, and nymphal activity peaked earlier in a generally hotter and drier year (2015) and later in years with cooler spring temperatures and increased seasonal precipitation (2016 and 2017). Although we found no statistical association between precipitation and nymphal density, the first significant rainfall since 2011 occurred in early 2016 (Western Regional Climate Center— https://wrcc.dri.edu), which corresponded with an annual increase in nymphal abundance at our study site. The long-term drought in California from 2011 to 2015 may have been a factor in the decreased nymphal density, prevalence, and magnitude of bimodal activity observed in 2015.
The initiation and cessation of I. pacificus nymph activity at our study site typically occurred within the same week each year (except cessation in 2017 was observed two weeks later than 2015 or 2016). Nymph activity started a month earlier than reported from some north coastal California locations (Eisen et al. 2003). The mean maximum temperature and relative humidity for 20 d prior to initiation of nymph activity at our site ranged from 15.6 to 21.1°C and from 94.4 to 97.5%, respectively (CVPD ranged from 3.0 to 9.2). However, CVPD is likely not predictive of the onset of seasonal I. pacificus activity. Daily vapor pressure deficits in January and early February are at a minimum and increasing temperatures are likely responsible for the initiation of nymphal activity (Belozerov 2009, Eisen et al. 2016).
Larger CVPD values (>52), corresponding with increased temperature and decreased relative humidity, typically resulted in lower nymph activity, although there was variation depending on the questing substrate (Fig. 3). However, we only measured temperature and relative humidity at one location and measurements of these variables at the substrate level would provide a clearer picture of their effect on the variation we observed in nymphal activity. Nymphal activity ceased when CVPD ranged from 106.7 to 126.1 and the 20-d prior mean maximum temperature and relative humidity ranged from 37.0 to 42.2°C and from 69.6 to 74.2%, respectively. The overall post-peak decline and cessation of nymphal density and prevalence observed each year was likely influenced by a rapidly drying environment (greater CVPD), thereby increasing the rate of desiccation and subsequent mortality of nymphs as observed by similar studies (Ogden et al. 2004, Diuk-Wasser et al. 2010, Eisen et al. 2010). Eisen et al. 2002 observed that I. pacificus nymphs consistently start to decline when mean maximum daily temperature exceeded 21–23°C and mean maximum daily relative humidity decreased below 83–85%. Similarly, Berger et al. (2014) showed a negative association between the cumulative number of hours below 82% and I. scapularis (a closely related species) nymphal abundance. We observed near perfect collinearity between the cumulative hours below 82% relative humidity and CVPD; therefore, our findings agree with the general patterns observed by previous studies (Eisen et al. 2002, Berger et al. 2014) and show a clear association between the interaction between temperature and humidity, and nymphal I. pacificus density and prevalence.
Public Health Relevance
Our results show that I. pacificus nymphs are more abundant, and therefore present a greater seasonal Lyme disease risk, in this northern Sierra foothill habitat than we previously estimated by dragging only leaf litter. The initial peak period of nymphal exposure risk at our study site occurred earlier in the year (April) and was shorter (March–May) than coastal counties, although nymphal activity associated with shaded rocky habitats resulted in a second peak (early June), effectively extending the observed activity period. However, our study site may not be representative of the region and the presence and abundance of ticks in the Sierra foothills is likely locally variable due to several factors. For example, our study site was relatively low in elevation and I. pacificus can be more abundant in oak woodland habitats at higher elevations, up to approximately 1,500 m, in the Sierra Nevada foothills (CDPH, unpublished data). There is likely a diversity of local nymphal abundance and temporal activity at higher elevation areas, driven partially by canopy cover, overall habitat type, vapor pressure deficit, and the variable effect of its constituent parts (temperature and humidity). For example, at a nearby location, but approximately 400 m higher in elevation, I. pacificus nymphs were found on birds from late April until June (Wright et al. 2000). The interaction among the environmental variables and I. pacificus density and prevalence evaluated at Salmon Falls merits further investigation at other locations in the Sierra foothills to better describe Lyme disease risk and compare nymph abundance in the varied ecological habitats in this region of California. Future studies should further evaluate I. pacificus nymph density and prevalence and the interaction between climate variables and substrates at a finer (microhabitat) level, which may explain some of the observed seasonal variation between substrates.
The prevalence of infected nymphs observed in this study helps to document the risk of exposure to B. burgdorferi sl in this area of California. Although we found no significant environmental associations, it is likely that among year variation in nymphal infection prevalence was due to several environmental factors identified by previous studies (Eisen et al. 2002, 2010). However, the number of samples required to observe any significant associations is likely far larger than were collected for this study. Longer-term studies in the region are needed to determine any fine-scale patterns of B. burgdorferi sl prevalence in I. pacificus nymphs. Regardless, nymphal infection prevalence at Salmon Falls was greater than the state-wide average for both 2015 (4.8%) and 2017 (5.0%), suggesting elevated risk to the public compared with many other areas in California (CDPH 2015, 2017).
Our findings support previous studies showing that contact with wood (trees or logs) may increase the risk of exposure to I. pacificus nymphs (Lane et al. 2004, 2007, 2009), particularly during spring months. We also identified an increased risk of exposure to I. pacificus nymphs associated with contact with rocks, comparable to that of wood substrates. This is of particular public health concern as rocks are often used as a resting place for those recreating in tick-abundant habitats. Given an overall B. burgdorferi sl prevalence of 7.7% of nymphs collected on rocks and the increased likelihood of encountering more than one nymph on a rock, we recommend including information about this potential exposure site in public education materials and evaluating the utility of this substrate type in future broad-scale studies and surveillance protocols.
Acknowledgments
We thank the California State Parks for providing access to Folsom Lake State Recreation Area trails and Joshua Ogawa and Stan Wright for preliminary sampling, which demonstrated the utility of sampling various environmental substrates for I. pacificus nymphs. We also thank Denise Bonilla, who suggested the potential of rocks as a productive nymphal substrate. We would also like to thank interns and staff from the California Department of Public Health, Vector-Borne Disease Section, who assisted with data collection, and especially Kerry Padgett, Vicki Kramer, and Anne Kjemtrup for internal review of the manuscript.
