Climate change continues to alter the behavior and distribution of species worldwide, with major ramifications for the transmission and risk of infectious diseases, including those caused by zoonotic vector-borne pathogens. This study explores the potential implications of climate change for one such pathogen, Borrelia burgdorferi (the causative agent of human Lyme disease), in Ixodes pacificus ticks of the far-western United States. Nymphal tick infection prevalence and density are compared against several metrics for climate, while also accounting for habitat fragmentation, mammalian species richness, and rodent tick burden to eliminate confounding variables. Findings show that climate extremes, such as those forecast with climate change, correlate with a reduction of B. burgdorferi prevalence in nymphal ticks despite nominal impacts to uninfected tick density, contrasting traditional hypotheses that these changes will increase vector-borne pathogens.

Study system

Data collection for this study took place between April and July from 2016 to 2022. For primary purposes of this study, land patches were designated for field sampling within a four-hour radius of San Francisco State University, San Francisco, California and consisted primarily of oak and scrub vegetation in a variety of settings ranging from urban to protected national forest. Designated patches ranged from 2.5 to 151,123 hectares in size and consisted of 14 different sites - nine designated in 2016, an additional two in 2018, and three more in 2019. Once established, 0.5-hectare sampling grids were marked on each plot with specific attention to canopy cover, slope, and edge proximity, with all sites being established under canopy cover a minimum distance of 20 m from the habitat edge and avoiding north-facing slopes.

Tick collection

Ticks were collected from each plot while questing using a standardized drag-sampling approach twice per season for a total sampled area of 995 m². Within each cohort, ticks were identified for life stage and species using microscopic identification of visual features according to previously established protocols. Once identified, only questing nymphal ticks were used for tick density in this study.

Bacterial diagnostics

The presence or absence of B. burgdorferi in each tick sample was identified using PCR analysis. Whole nymphs were stored at -80 °C until extraction for DNA, which was performed using the Qiagen DNA Easy Blood and Tissue Kit according to manufacturer instructions (Qiagen Redwood City Inc., Redwood City, California, USA; PN 69504). Extracted DNA was then tested using a 5S-23S rDNA PCR protocol to identify the presence of B. burgdorferi sensu lato alongside gel electrophoresis and DNA sequencing, which was used to further distinguish species of pathogen down to B. burgdorferi sensu stricto.

Mammalian species richness and tick burden

Counts for mammalian species richness and rodent tick burden were supported through a combination of mark-recapture and trail camera techniques. Rodent captures were performed for three consecutive days a year at each field site between the months of April and July for all years except 2020 and 2022. Rodents were captured in Sherman live traps (H.B. Sherman Traps, Tallahassee, FL, USA) set at 11.8 m intervals. In total, 49 traps were placed per site. Captured rodents were weighed, sexed, and identified to species using established protocols, which was confirmed with Cytochrome-B molecular analysis, if needed. Attached ticks were collected directly from trapped rodents for rodent tick burden counts. Additionally, two motion detector Bushnell cameras (Bushnell models #119736, #119836, #119836C, Bushnell, KS, USA) were anchored to trees approximately ~1 m above the ground for each field site. These cameras faced in opposite directions and recorded a total of 40 days and nights of wildlife activity for all years except 2020 and 2022. Recordings were used to visually identify all species present at each sample site using morphological characteristics. All mammal trapping was done in accordance with IACUC standards via San Francisco State University (AU16-05R1a and AU19-01R2) and was performed with California Department of Fish & Wildlife approval (collection permits SC-8407 and S-203370009).

Land fragmentation and climate parameters

All landscape measurements were performed manually using ArcGIS software (copyright ESRI) and data layers from the California Department of Forestry and Fire Protection. Climate parameters of interest were provided by VisualCrossing, which is a powerful weather API tool that interpolates between weather stations to fill in missing values and cleanses data for accuracy (Visual Crossing Corporation, Reston, Virginia, USA). As such, it has been validated for measurements at 1 km and 5-minute time point resolution. Captured data included averages for the parameters of interest for the month of April from all sites and years, as this month showed the most extremes for analysis and fell during the peak questing months for tick collection.

Data analysis

All response variables of interest were explored for normality, collinearity, Cook’s Distance, and leverage points using the plot function in the R “stats” package, as well as residual plots and histograms. As all data of interest were considered counts with non-normal distribution and overdispersion, models were run across all years using a Penalized Quasilikelihood GLMM with a Poisson distribution and log link. A random effect was included for site to account for non-independence among the explanatory variables. Stepwise model selection was performed using the “dredge” function in the MuMIn R packages. Year was initially explored as a random effect in all models, but was ultimately dropped due to poor fit. Final model output and all data was back-transformed, and confidence intervals established using the “exp” and “confint” functions in the R “stats” package. R squared values to determine goodness-of-fit were also obtained for all final models using the “R.squared” function in the “aod” package of R.

A PLS-PM was used to assess relationships between climate, landscape, host, and tick parameters of interest for nymphal infection prevalence. This method was chosen due to its ability to explore potentially statistically significant relationships between observational data without strict criteria for data structure. For its construction, latent variables (theoretical outcomes) of interest were described, as well as the path (structural relationship) between each latent variable. Each latent variable was then defined by a second set of manifest variables (measured variables relating to the theoretical concept), and the data structure was explored for statistically significant correlations (i.e., regression coefficients) between manifest variables and their latent variable. Manifest variables with poor relationship to their latent variable (regression coefficient < 0.40; Dillon-Goldstein’s rho < 0.70) were dropped from the final model. The final PLS-PM model was run using the ‘plspm’ package in R and all outputs were explored for goodness-of-fit to determine their influence on the overall model outcome. Bootstrapping was performed to confirm model output and all findings were consistent.