Predicting the current and future distribution of the western black-legged tick, Ixodes pacificus, across the Western US using citizen science collections

Abstract

In the twenty-first century, ticks and tick-borne diseases have expanded their ranges and impact across the US. With this spread, it has become vital to monitor vector and disease distributions, as these shifts have public health implications. Typically, tick-borne disease surveillance (e.g., Lyme disease) is passive and relies on case reports, while disease risk is calculated using active surveillance, where researchers collect ticks from the environment. Case reports provide the basis for estimating the number of cases; however, they provide minimal information on vector population or pathogen dynamics. Active surveillance monitors ticks and sylvatic pathogens at local scales, but it is resource-intensive. As a result, data are often sparse and aggregated across time and space to increase statistical power to model or identify range changes. Engaging public participation in surveillance efforts allows spatially and temporally diverse samples to be collected with minimal effort. These citizen-driven tick collections have the potential to provide a powerful tool for tracking vector and pathogen changes. We used MaxEnt species distribution models to predict the current and future distribution of Ixodes pacificus across the Western US through the use of a nationwide citizen science tick collection program. Here, we present niche models produced through citizen science tick collections over two years. Despite obvious limitations with citizen science collections, the models are consistent with previously-predicted species ranges in California that utilized more than thirty years of traditional surveillance data. Additionally, citizen science allows for an expanded understanding of I. pacificus distribution in Oregon and Washington. With the potential for rapid environmental changes instigated by a burgeoning human population and rapid climate change, the development of tools, concepts, and methodologies that provide rapid, current, and accurate assessment of important ecological qualities will be invaluable for monitoring and predicting disease across time and space.

About the Authors:

W. Tanner Porter

Roles Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Visualization, Writing – original draft, Writing – review & editing

* E-mail: [email protected]

Affiliations Department of Biological Sciences, Northern Arizona University, Flagstaff, AZ, United States of America, Translational Genomics Research Institute, Flagstaff, AZ, United States of America

ORCID logo https://orcid.org/0000-0002-2611-0924

Zachary A. Barrand

Roles Investigation, Writing – review & editing

Affiliation: Department of Biological Sciences, Northern Arizona University, Flagstaff, AZ, United States of America

ORCID logo https://orcid.org/0000-0001-9281-7314

Julie Wachara

Roles Investigation, Writing – review & editing

Affiliation: Department of Biological Sciences, Northern Arizona University, Flagstaff, AZ, United States of America

Kaila DaVall

Roles Investigation

Affiliation: Department of Biological Sciences, Northern Arizona University, Flagstaff, AZ, United States of America

ORCID logo https://orcid.org/0000-0001-8320-6618

Joseph R. Mihaljevic

Roles Methodology, Writing – review & editing

Affiliation: School of Informatics, Computing, and Cyber Systems, Northern Arizona University, Flagstaff, AZ, United States of America

ORCID logo https://orcid.org/0000-0003-2320-5773

Talima Pearson

Roles Methodology, Writing – review & editing

Affiliation: Pathogen and Microbiome Institute, Northern Arizona University, Flagstaff, AZ, United States of America

ORCID logo https://orcid.org/0000-0003-3696-041X

Daniel J. Salkeld

Roles Conceptualization, Methodology, Writing – review & editing

Affiliation: Department of Biology, Colorado State University, Fort Collins, CO, United States of America

Nathan C. Nieto

† Deceased.

Roles Conceptualization, Funding acquisition, Methodology

Affiliation: Department of Biological Sciences, Northern Arizona University, Flagstaff, AZ, United States of America

Introduction

Tick-borne disease (TBD) diagnoses have steadily risen over the last 20 years across the US, emphasizing the increasing importance of such zoonotic diseases [1]. The rise in TBD has been attributed to several factors, including changes in climate and land-use patterns that influence vector distribution and densities [2, 3]. Historically, TBD surveillance has been conducted at the county, state, or national scale through the National Notifiable Diseases Surveillance System (NNDSS), which requires providers and public health agencies to report cases of specific diseases to the Centers for Disease Control and Prevention (CDC). Several TBDs (Lyme disease, babesiosis, anaplasmosis, tularemia, and rocky mountain spotted fever) are reportable to the CDC [4]. The NNDSS provides important insight into case burdens; however, there are several pitfalls. First, only diagnosed cases are reported, leading to possible under-reporting or over-reporting [5]. Second, cases are reported based on county of residence and not the location of exposure [6]. Third, the presence and population dynamics of the vectors are not considered. These characteristics limit the ability of NNDSS data to predict and explain TBD risk and vector distribution.

Measures of TBD risk, the likelihood of exposure to infected ticks, have traditionally relied on active vector surveillance, requiring researchers to collect ticks from the environment. These collections are most often conducted by either flagging a piece of cloth across vegetation or dragging a piece of cloth across the ground to collect questing ticks [7]. Active vector surveillance has been extensively used to identify and model tick distribution and densities, pathogen prevalence, nymphal infection prevalence, pathogen diversity in sylvatic cycles, seasonality of vectors, risk of spillover, and suitable habitat [8–13]. Active surveillance has long been considered the “gold-standard” for sylvatic TBD surveillance; however, it is resource-intensive, requiring the majority of spatial modeling approaches to aggregate multiple years of data [10, 11], resulting in a loss of temporal resolution (Table 1). For example, some distribution models have utilized active surveillance, collected over three decades, or historical data, collected over the last century, to model current tick distributions. In a stable environment, these methods would likely not impact the outcomes; however, in today’s dynamic world of changing weather patterns, climate, and land use, aggregation of data across multiple years or decades could impact the validity and accuracy of model predictions [14]. Ideally, vector distribution and movement data should be related as closely as possible to appropriate environmental conditions. To achieve that, surveillance methods may be more successful if they trend towards low cost, continuous, rapid, and spatially diverse sample collection.

[Figure omitted. See PDF.]

Table 1. Previously published models identifying tick distribution through species distribution models.

https://doi.org/10.1371/journal.pone.0244754.t001

Opportunistic citizen science tick collections have the potential to fulfill these characteristics. Citizen science or passive surveillance campaigns have been implemented at both state and national scales by encouraging members of the public to mail ticks to a laboratory for subsequent identification or pathogen testing [32–36]. These programs have been proven to be an efficient method for tick collection and can help to inform on a variety of topics that cannot be easily gathered through a single surveillance technique [32]. Citizen science tick collections have generated large datasets that have provided insights on vector spread, links between tick submission frequencies and reported human TBD cases, seasonality of human-tick exposure, human activity during tick exposure, and pathogen prevalence in submitted ticks [30, 33–35, 37, 38]. Additionally, opportunistic and passive surveillance programs have been used to model or inform on a wide variety of TBD topics. For example, these programs have been used extensively to identify the tick species that people interact with, pathogen prevalence data, and the timing of human-tick exposure at a national level or finer scale [32–34, 36–40]. Recently, these data have also been utilized to accurately model human case data from the number of tick submissions received [35, 38].

A promising avenue of TBD research uses citizen science to collect data on locations of ticks and employs species distribution models (SDMs) to identify the distribution or niche of specific tick species, e.g., this technique created SDMs for Ixodes scapularis in Canada and validated predictions with active surveillance [30]. SDMs are commonly used in vector ecology and conservation biology, as they explore variables that may dictate distribution and identify areas potentially supporting current populations. Additionally, SDMs can be projected across future climate scenarios to help predict species distribution changes in response to climate change [41]. These models often incorporate both abiotic (elevation, climate, etc.) and biotic (land cover, other species, etc.) variables as predictors [41]. SDMs have been previously used to model tick distributions of several tick species in North America (Table 1). The majority of these models rely on presence records derived from active surveillance that have been aggregated across time (e.g., over several decades) and/or space (e.g., county-level presence/absence data), inhibiting the ability to identify accurate distributions of changing populations.

Monitoring the current and future distribution of medically important ticks has significant public health implications and can assist in identifying high disease risk areas. On the US West Coast (California, Oregon, and Washington), Ixodes pacificus, the western black-legged tick, is the most medically important tick vector and is responsible or implicated in the transmission of Borrelia burgdorferi sensu stricto (Lyme disease), B. miyamotoi (hard-tick relapsing fever), and Anaplasma phagocytophilum (anaplasmosis) [42–47]. Here we investigate the application of the dataset from our free national citizen science tick collection program to model the current and future distribution of I. pacificus across the Western US. We hypothesize that the ability to rapidly collect large datasets with citizen science will provide comparable results to previous SDMs and will supplement traditional surveillance programs.

Materials and methods

Citizen science tick collection

Between January 2016 and December 2018, 18,881 ticks were submitted by citizen scientists who contributed to a free tick collection program [32]. The details of the citizen science program can be found in Nieto et al. (2018), but we provide a brief overview here. Citizen scientists were recruited to the program through an initial public relations campaign and a public website (https://www.bayarealyme.org/lyme-disease-prevention/tick-testing/). Additionally, individuals and TBD public awareness groups, not related to the funders or researchers, disseminated information about the program though social media and other advertising platforms. The citizen scientists participated in the program by mailing collected tick(s) and associated data (location of exposure) to the laboratory. The submission form covered information related to locations and characteristics of exposure (date, activity, environment, etc.). Ticks were identified to species using morphological characteristics and tested for the presence of several TBD pathogens.

For the purposes of this study, I. pacificus submissions from the West Coast with GPS points that corresponded to the associated reported county of exposure were included; though further accuracy of the GPS points was not verified by the researchers. GPS points were converted to correspond to the center of each 1km x 1km grid cell and spatially duplicated grid cells were removed from the dataset, so each 1 km x 1 km pixel had a maximum of a single presence point.

Predictor variables

We collected a combination of seasonal climatic and land cover variables to build our species distribution model. We had two goals for this analysis. First, to predict the current distribution of I. pacificus (ecological niche model), we used a model that had both abiotic (climate and elevation) and biotic variables (landcover) to accurately predict the true ecological niche or distribution of I. pacificus across today’s landscape. Second, we sought to predict I. pacificus’ future distribution in response to climate change. To assess these changes and predict the future distribution we created a climate niche model using the current climate (bioclimatic variables only) and predicted it into the future using CMIP5 climate projections. We were limited by the available variables in CMIP5 climate projections, therefore, we had to remove several variables from the ecological niche model (e.g., land cover, seasonal data). Removal of these variables removes important factors that impact tick distribution, however, including these variables without proper future projections would not allow us to accurately model the potential future distribution.

Biotic variables (Climate and elevation)

We collected daily measures of three weather variables (minimum temperature, maximum temperature, and precipitation) from Daymet for the years 2009–2018, which were processed to produce nineteen summary statistics for use in the SDMs [48]. All analyses were performed in the statistical package “R,” version 3.5.2 [49]. Daymet files were downloaded and processed in R through the “FedData” package with a native resolution of 1km x 1km [50]. The maximum temperature, minimum temperature, and cumulative precipitation were calculated across each month and averaged across each year. Bioclimatic variables were then derived through the “biovars” function within the “Dismo” package [51]. In addition to these bioclimatic variables, day length, solar radiation, snow water equivalent, and vapor pressure were also incorporated as predictor variables. Similar to the bioclimatic variables, these were downloaded from “Daymet” through the “FedData” package and were averaged across each month and year between 2009–2018. Additionally, these variables were also averaged across each season (Winter: December-February; Spring: March-May; Summer: June-August; Fall: September-November). Finally, an elevation layer was also included with a native resolution of 1km x 1km [52].

Abiotic variables (Land cover)

Four National Land Cover Database (NLCD) variables were utilized to explain biotic interactions [53]. The native resolution of the NLCD data was 0.03km x 0.03km and was aggregated into four separate rasters with a 1km x 1km resolution, to match the presence data and climate data resolution, by calculating the percentage of each 1km x 1km grid cell that contained (1) forest, (2) scrub, (3) urban low-density, or (4) urban high-density cells. This produced four final variables (% forest, % scrub, % urban low-density, and % urban high density) at a 1km x 1km resolution, while differentiating cells that are heterogeneous. Forest was defined as deciduous forest (NLCD: [41]), evergreen forest (NLCD: [42]), and mixed forest (NLCD: [43]). Scrub was classified using shrub/scrub (NLCD: [52]) and grassland/herbaceous (NLCD: [71]). Urban low density was classified as developed open space (NLCD: [21]) and developed low intensity (NLCD: [22]). Finally, high intensity urban was defined as developed medium intensity (NLCD: [23]) and developed high intensity (NLCD: [24]).

Species distribution modeling

SDMs were created through “R,” utilizing maximum entropy models accessible through the “ENMeval” package [54]. Maxent models were produced using the maximum entropy algorithm [55, 56] and were assessed through the utilization of area under the curve (AUC) calculations where a value of 1.0 indicates a model that perfectly can classify presence vs. absence (e.g. no false positives or false negatives) [56]. Maxent was utilized as it is specifically designed to handle presence only datasets and utilize background points instead of absence locations. A total of 1600 background points were randomly selected from the study area (California, Nevada, Oregon, Washington, and Idaho). In Maxent models, background points are compared to the presence points to identify specific patterns that are indicative of a species niche. This allows the AUC statistic to differentiate inter-model performance, which is described as ability to identify presence points vs. background points. However, this statistic should be interpreted with caution when assessing the overall model performance (e.g. accuracy), which is heavily influenced by the ratio of presence points to background points and the area of interest [56, 57].

Preliminary models were built to select and identify high contribution variables. Within each stack, the presence and background locations were bootstrapped, from the original datasets, 100 times and subsequently modeled. After locations were bootstrapped, locations were split into training and testing data sets utilizing “checkerboard2” algorithm from the ENMeval package. Feature classes were not restricted in the models, and 0.1, 0.25, 0.5, 1, 2, 4, 8, and 10 were used as regularization multipliers accounting for a total of 48 unique models within each bootstrap iteration. In total, 4,800 models were built for each variable stack (suitable habitat model and climate niche model) during the preliminary analysis.

Variable selection was performed with a similar method to Jueterbock et al. (2016). Differing from Jueterbock et al. our methods incorporated bootstrapped model evaluation and do not serially remove variables, instead our methods utilize moderate to high performing models to remove variables that had low maximum contribution values across all models [58]. Preliminary models with a testing AUC greater than 0.8 were used to remove variables that had a maximum contribution across all models of less than 10 percent. Remaining variables were added to a final stack based on average percent contribution (high to low). Variables with high correlation (< -0.8 and > 0.8 Spearman correlation) to previously included variables were removed (S1 Table).

After predictor variable reduction, presence and background points were again bootstrapped 100 times and 4,800 models were generated for each raster stack (ecological niche model and climate niche model). The model with the highest average AUC across 100 bootstrap iterations was selected as the best model and further analyses were conducted with that model. If any variables had a percent contribution of less than 0.5%, it was removed to simplify the model. The ecological niche model of I. pacificus was then predicted over the bioclimatic/climate/elevation/landcover variables to identify current suitable habitat for I. pacificus. Maxent raw output was converted to three sensitivity threshold values by forcing 90%, 95%, or 99% of collected presence points to be predicted as suitable habitat. The 90% and 99% sensitivity thresholds predictions for California were then qualitatively and quantitatively compared to the 90% and 99% ensemble raster from a previous SDM of I. pacificus in California, generated by active collection of ticks by public health biologists [11]. For the comparison, areas that had 1 or more models predicting suitable habitat in the Eisen et al. ensemble raster were considered as suitable habitat. Our predictions were converted to the same spatial extent (California). We then calculated percent agreement of suitable habitat between the models (i.e. percent area that was predicted as suitable in both models) () and overall percent agreement between both projections (i.e. percent area that was predicted as suitable or unsuitable in both models) ().

The climate niche model (bioclimatic variables only) was predicted over future climate predictions using the CMIP5 multi-model ensemble to identify the future climate niche of I. pacificus [59]. Future projections were computed for two Representative Concentration Pathways (RCP) emission scales (RCP 6.0, and 8.5) at 2.5-degree resolution [60]. Suitable habitat estimates in 2050 were generated through the use of 12 (RCP 6.0), and 17 (RCP 8.5) different climate models. Projections were set to a threshold based on a 95% sensitivity and projections from each RCP value were combined to create a final projection. The final projection shows the number of models that predict suitable habitat across the area of interest. A complete list of models and average bioclimatic variables for the area of interest is available in S2 Table.

Results

Citizen science I. pacificus collection

Of the 2,332 I. pacificus submitted, 767 (33%) were submitted with valid GPS position information that corresponded to the reported county of exposure. Multiple submissions within each predictor variable pixel (1km x 1km) were removed, producing a total of 477 unique presence points across the study area: California (n = 397), Oregon (n = 54), and Washington (n = 26) (Fig 1).

[Figure omitted. See PDF.]

Fig 1. Distribution of collected I. pacificus (n = 477) corresponding to unique presence points across the Western US.

https://doi.org/10.1371/journal.pone.0244754.g001

Ecological niche model

The final ecological niche model utilized a total of eight predictor variables: average vapor pressure in the spring, isothermality (BIO3), temperature seasonality (BIO4), land cover- forest, land cover- scrub, land cover- urban low density, average day length (summer), and precipitation of the driest month (BIO14). Average vapor pressure in the spring was the largest contributor (36.7%) to the distribution identifying an increase in the likelihood of suitable habitat with increasing vapor pressure. Isothermality (BIO3) was the second largest contributor (30.4%), predicting increased likelihoods with increased isothermality (Figs 2 and 3; S3 Table). Isothermality (BIO3) had the largest permutation importance with 46.8%, followed by average vapor pressure (spring), and land cover- forest (Fig 3). At a 90% and 95% sensitivity threshold, final predictions showed widespread predicted habitat across the coastal areas of California, western Oregon and western Washington (Fig 4). Additionally, suitable habitat was identified along the Sierra Nevada foothills.

[Figure omitted. See PDF.]

Fig 2. Predicted variable response curves for the suitable habitat model with a linear and quadratic feature class and regularization multiplier of 10.

https://doi.org/10.1371/journal.pone.0244754.g002

[Figure omitted. See PDF.]

Fig 3. Variable contribution (A) and permutation importance (B) in the final ecological niche model.

https://doi.org/10.1371/journal.pone.0244754.g003

[Figure omitted. See PDF.]

Fig 4. Current suitable habitat of I. pacificus across the Western US based on climate and land cover variables.

Normalized raw output (A) and suitable habitat based on a 90% sensitivity threshold (Blue), 95% sensitivity threshold (Green), and 99% sensitivity (Gray) (B).

https://doi.org/10.1371/journal.pone.0244754.g004

Qualitatively, our model predictions of California were similar to an earlier SDM of I. pacificus distribution that used active surveillance by public health biologists (Eisen et al. 2018) (Fig 5A). Quantitatively, when comparing maps that used the 90% sensitivity threshold (i.e. 90% of the citizen science GPS points are predicted as suitable habitat), overall suitable habitat projections had a 51% overlap of suitable habitat between the two projections. When considering overall agreement, the models had an 80% agreement. The majority of the differences in predictions arise from our model predicting suitable habitat within the Sacramento/San Joaquin Valley and urbanized areas (e.g., San Francisco, Los Angeles) (Fig 5A, green areas). Additionally, Eisen et al.’s predictions identify additional suitable habitat in adjacent areas to areas that were predicted as suitable habitat in both models (Fig 5A, red areas).

[Figure omitted. See PDF.]

Fig 5. Comparison of predicted suitable habitat presented here (Porter Suitable) and presented in Eisen et al. 2018 at a 90% sensitivity (A) and a 99% sensitivity (B).

https://doi.org/10.1371/journal.pone.0244754.g005

Qualitatively, at a 99% sensitivity, the predictions produced here begin to diverge from Eisen et al. predicting widespread suitable habitat across the study area. Quantitatively, suitable habitat predictions had a 64% overlap between both projections and an overall agreement of 76% (Fig 5B). These discrepancies arise from a lack of specificity at high sensitivities that are encountered with the citizen science generated model (Fig 5B, green areas).

To build the final model, the preliminary bootstrap analysis and variable selection identified a total of 16 (16/40) variables that were predictive of I. pacificus presence. In decreasing importance, these variables were: isothermality (BIO3), temperature seasonality (BIO4), percent of land cover with low density urban, vapor pressure during the spring, percent of land cover with forest, percent of land cover with scrub, precipitation of driest month (BIO14), average snow water equivalent during winter, percent of land cover with high density urban, precipitation of coldest quarter (BIO19), average day length during the summer, average elevation, average vapor pressure during the summer, precipitation seasonality (BIO15), mean diurnal range (BIO2), and average snow water equivalent in the summer. After variable selection, final bootstrap analysis produced a total of 1,832 (38%) models that had an average testing AUC greater than 0.9 and 3,568 (74%) models had an average testing AUC greater than 0.8. Overall, the models with a linear and quadratic feature classes and a regularization multiplier of 10 had the highest average testing AUC (AUC = 0.95, sd = 0.007). After final model selection, several variables still had minimal percent contribution (< 0.05%) and were removed in the final variable selection step; these included: elevation, percent of land cover with high density urban, precipitation seasonality (BIO15), precipitation of coldest quarter (BIO19), mean diurnal range (BIO2), average snow water equivalent in the summer, average snow water equivalent in the winter, and average vapor pressure in the summer. This created the final model with a total of 8 variables (S3 Table).

Climate niche model

When we restricted the data set to only the climate variables, so that we could use the model to forecast based on future climate projections, the climate niche model of I. pacificus showed similar patterns to the ecological niche model; however, in general, the climate niche model identified a more contiguous niche. Additionally, with a 0.99% sensitivity threshold, the model predicted a widespread climate niche across much of the area of interest (Fig 6). Ensemble future predictions using 2050 RCP 6.0 and RCP 8.5 climate projections identified an overall net loss in the climate niche of I. pacificus. The majority of this loss occurred along the boundaries of the climate niche of I. pacificus (Fig 7). A small amount of range expansion was predicted in the northern ranges of I. pacificus, however, this expansion was predicted in a minority of climate scenarios.

[Figure omitted. See PDF.]

Fig 6. Current climate niche of I. pacificus across the Western US based on averaged bioclimatic variables between 2009–2018.

Normalized raw output (A) and suitable habitat based on a 90% sensitivity threshold (Blue), 95% sensitivity threshold (Green), and 99% sensitivity (Gray) (B).

https://doi.org/10.1371/journal.pone.0244754.g006

[Figure omitted. See PDF.]

Fig 7. Current climate niche of I. pacificus based on a 95% sensitivity (A) and number of future climate scenarios (CMIP5) predicting suitable habitat across the I. pacificus climate niche model in 2050, based on an RCP 6.0 (B) and RCP 8.5 (C) and 95% threshold.

The grey line in panels B and C indicates outline of current climate niche.

https://doi.org/10.1371/journal.pone.0244754.g007

To produce the final climate niche model, the preliminary bootstrap analysis and variable selection identified a total of 7 (7/19) variables that were predictive of I. pacificus presence and included in the final projection. In decreasing percent contribution, these variables were: temperature seasonality (BIO4), isothermality (BIO3), precipitation of driest month (BIO14), minimum temperature of coldest month (BIO6), precipitation seasonality (BIO15), and mean diurnal range (BIO2) (Fig 8). However, permutation importance included precipitation of the coldest quarter (BIO19) and precipitation if the driest month (BIO14) as the most important variables accounting for ~60% of permutation importance (Fig 8). In general, bioclimatic variables with percent contribution greater than 10% were included in the ecological niche model and the climate niche model. Differences in permutation importance and percent contribution between the two presented models are likely due to additional variables (e.g. land cover and vapor pressure) being included in the ecological niche model. During the preliminary variable selection, the final bootstrap analysis produced a total of 2,033 (42%) models that had an average testing AUC greater than 0.9 and 2,454 (51%) models that had an average testing AUC greater than 0.8. Overall, models with a linear and quadratic feature classes and a regularization multiplier of 2 had the highest average testing AUC (AUC = 0.95, sd = 0.006) and the individual variable contribution was above the 0.5% cutoff for all of the variables (Fig 9; S3 Table).

[Figure omitted. See PDF.]

Fig 8. Variable contribution (A) and permutation importance (B) in the climate niche model.

https://doi.org/10.1371/journal.pone.0244754.g008

[Figure omitted. See PDF.]

Fig 9. Predicted variable response curves for the climate niche model with a linear and quadratic feature classes and regularization multiplier of 2.

https://doi.org/10.1371/journal.pone.0244754.g009

Discussion

Traditionally, SDMs identifying suitable habitat for ticks have relied on active surveillance or historical collections that frequently span decades (Table 1). We demonstrate that data from citizen science collections, collected in a two-year timeframe, generated SDMs for I. pacificus that closely reproduced predictions generated through traditional public health surveillance. Data from previous surveillance methods provide a valuable resource for validating new methods, assessing longitudinal changes, and providing long-term average species predictions. Citizen science collections complement these efforts by providing large and spatially diverse sample collections and can be rapidly generated.

In California, citizen science tick collection based suitable habitat and climate niche estimates, presented here, were similar to predictions based on SDMs that were previously generated through active surveillance over a ~30-year period, identifying suitable habitat and climate niches along the Californian coast, the northern bay area to Sacramento and along the Sierra Nevada mountains [11, 23]. Quantitatively, our predictions were extremely similar at a 90% sensitivity (a threshold that forces 90% of collected presence points to be predicted as suitable habitat). Areas that are not consistent could be artifacts of the modeling processes. For example, Eisen et al. restricted their predictions to areas that were classified by land-cover analysis to be forest, grass, scrub-shrub, or riparian, and excluded agricultural and urban areas. Alternatively, our predictions incorporated the percent of landcover types in a 1km x 1km grid cell as variables into the model, allowing suitable habitat to be predicted in areas that may not have been considered by the Eisen et. al model (Fig 5A). An example of this would be areas that were classified as urban but have enough fragments of natural habitat to support tick populations. Through this modeling difference, we include areas within our predictions that are regarded as unlikely habitats for I. pacificus (e.g., metropolitan areas and San Joaquin Valley) that may not have enough fragmented habitat to be suitable for I. pacificus populations. Additionally, Eisen et al. utilized climate data from Daymet aggregated across 35 years, while we used data aggregated across 10 years.

Our 99% sensitivity projections significantly diverged from those presented by Eisen et al., predicting large-scale suitable habitat across California (Fig 5B). This is likely an artifact of extreme sensitivities (> 95%) and citizen science sample collections, which can be prone to inaccurate or mistaken geographical information [32]. With high sensitivities, these samples are forced to be classified as suitable habitat and cause the projections to diverge from previous work. These uncertainties need to be considered when utilizing citizen science collected datasets and can be resolved through using lower sensitivity ranges (< 95%) or through incorporation of targeted active surveillance campaigns to verify tick presence in select locations.

Citizen science approaches have weaknesses: e.g., the impacts of non-structured data collection can produce spatial biases within datasets [61]. These factors are indeed challenging to control for and result in some level of bias. Citizen science tick collections inherently collect ticks that are found by citizens over the collection period; thus, our presence locations reflect where ticks and humans interact, introducing geographical and temporal biases associated with variable sampling effort across space and time, which can over and underestimate distributions. Additionally, our presence points are likely influenced by advertising strategies of the collection program, probably selecting for a high response rate from the San Francisco Bay Area and California. These biases were accounted for within our analysis through the use of random background points and utilizing 90% and 95% sensitivity thresholds, which were more representative of species distributions than continuous predictions. Such spatial biases within presence-only datasets are not specific to citizen science collection methodologies. For example, active surveillance tick collections can suffer from similar biases where ticks are only collected from a few ecological environments (e.g., dense woodland, Eisen et al. 2006), or areas that have established abundant tick populations (for pathogen surveillance); thus, failing to represent tick populations inhabiting different ecological environments (e.g., coastal chaparral), or low-density populations.

Another limitation, specific to citizen science, is that people often do not know precisely where they contacted the tick. Citizen scientists were asked to report where they were exposed to the tick; however, we expect some level of spatial uncertainty. We attempted to account for such uncertainty by using a 1km x 1km spatial resolution and relying on lower sensitivity thresholds. We believe that this variability could drastically impact a fine-scale area of interest modeling campaign (sub-county level) due to the uncertainty of the exact location of the tick exposure, however, at a large-scale (state or regional level), this effect would be relatively minor, and the sample size of citizen science collections likely overpowers this uncertainty spatial uncertainty.

Despite these challenges, citizen science provides a method that supplements other surveillance techniques, providing large-scale sample collection and species distribution monitoring at a spatial and temporal resolution that would be impractical with any other surveillance method. For example, our tick collection program collected specimens from geographical areas (i.e., Oregon and Washington) with sparse active surveillance data. As a result, we could predict an ecological niche of I. pacificus in Oregon and Washington, an area where suitable habitat for I. pacificus, via species distribution model, has not previously been defined. Our predicted current ecological niche in Oregon includes a distribution that extends north to south along the western third of the state, consistent with previous active surveillance campaigns which collected I. pacificus in the western third of the state [62]. Similarly, in Washington, our model predicts patchy suitable habitat along the western third of the state, which is consistent with active surveillance [63].

As tick distributions continue to change due to ecological changes and climate change, predicting the potential future distribution will provide insight into disease/pathogen shifts that could occur, allowing public health agencies to anticipate and combat these shifts. Utilizing future climate predictions (2050, RCP8.5), we identified a 30% loss in the I. pacificus climate niche. However, the accuracy of future model projections relies on shifts in climate and biotic environments to the extent that the environment is no longer habitable for a target population. This model relies on temperature seasonality (BIO4), which is predicted to increase in the future resulting in a potential loss of suitable habitat in the future predictions (Figs 6 and 7). The shrinking climate niche presented seems to counter the current range expansion of other North American tick populations [3, 64], some of which have been attributed to climate change, bird migrations [65–67], and human actions that facilitate tick movement (i.e., livestock transportation) [21].

Conclusions

As we continue to grapple with tick-borne diseases, we need to supplement our processes to match the changing and growing needs. Citizen science tick collections and other collection methods, such as mobile smartphone application-based citizen science programs (e.g., The Tick App [68]) and photographic identification of passively collected ticks (e.g., TickSpotters [69]), have the potential to fulfill this need and allow for a variety of questions to be answered that relate to human exposure, clinical case data, and species distributions. Changing ecological conditions will result in an inevitable shift in tick and disease patterns. Tick surveillance is critical for predicting disease risk and must, in order to be most effective, capture small spatial and temporal variations. Citizen science tick collections have the potential to fulfill this need and allow for a variety of questions to be answered that relate to human exposure, clinical case data, and species distributions. Citizen science collections do have some challenges (i.e., spatial biases, biased sampling) that add complexities into analyses; however, the ability to collect extensive and spatially/temporally diverse data, with limited effort, has significant potential across scientific fields.

Supporting information

[Figure omitted. See PDF.]

S1 Table. Correlation matrix of all considered variables.

https://doi.org/10.1371/journal.pone.0244754.s001

(CSV)

S2 Table. Table of current and future bioclimatic variables from Daymet (2009–2018) and CMIP5 database (Future Scenarios) that have been averaged across the area of interest.

https://doi.org/10.1371/journal.pone.0244754.s002

(CSV)

S3 Table. Predictor variables considered and included in the species distribution model of I. pacificus in the Western US.

https://doi.org/10.1371/journal.pone.0244754.s003

(CSV)

Acknowledgments

We earnestly thank the army of citizen scientists who made these analyses possible. We are grateful to Dr. Brad Butterfield and Dr. David Engelthaler for their valuable discussions and insight; Jason Travis and the entire Nieto Lab for their support and assistance; and Shane Feirer and Kerry Padgett for producing and making available the results from their models presented in Eisen et al. (2018). We thank the anonymous reviewers who provided thoughtful and helpful comments, which improved the quality of the manuscript. We acknowledge the World Climate Research Programme’s Working Group on Coupled Modelling, which is responsible for CMIP, and we thank the climate modeling groups (listed in S2 Table) for producing and making available their model output. For CMIP the U.S. Department of Energy’s Program for Climate Model Diagnosis and Intercomparison provides coordinating support and led development of software infrastructure in partnership with the Global Organization for Earth System Science Portals. Finally, we acknowledge the U.S. Census Bureau for producing and making available county and state shapefiles.

Citation: Porter WT, Barrand ZA, Wachara J, DaVall K, Mihaljevic JR, Pearson T, et al. (2021) Predicting the current and future distribution of the western black-legged tick, Ixodes pacificus, across the Western US using citizen science collections. PLoS ONE 16(1): e0244754. https://doi.org/10.1371/journal.pone.0244754

References

1. Paddock CD, Lane RS, Staples JE, Labruna MB. Global health impacts of vector-borne diseases [Internet]. Mack A, editor. Global Health Impacts of Vector-Borne Diseases: Workshop Summary. Washington, D.C.: National Academies Press; 2016. 221–257 p. Available from: https://www.ncbi.nlm.nih.gov/books/NBK390439/

2. Ostfeld RS, Brunner JL. Climate change and Ixodes tick-borne diseases of humans. Philos Trans R Soc Lond B Biol Sci [Internet]. 2015 Apr 5;370(1665):20140051. Available from: http://www.ncbi.nlm.nih.gov/pubmed/25688022 pmid:25688022

3. Sonenshine D. Range Expansion of tick disease vectors in North America: implications for spread of tick-borne disease. Int J Environ Res Public Health [Internet]. 2018 Mar 9;15(3):478. Available from: http://www.mdpi.com/1660-4601/15/3/478 pmid:29522469

4. Centers for Disease Control and Prevention. National Notifiable Diseases Surveillance System (NNDSS) [Internet]. 2019. Available from: https://wwwn.cdc.gov/nndss/

5. Mead PS. Epidemiology of Lyme disease. Infect Dis Clin North Am [Internet]. 2015 Jun;29(2):187–210. Available from: healthvermont.gov pmid:25999219

6. Beard CB, Visser SN, Petersen LR. The need for a national strategy to address vector-borne disease threats in the United States. J Med Entomol [Internet]. 2019 Sep 3;56(5):1199–203. Available from: https://academic.oup.com/jme/article/56/5/1199/5558470 pmid:31505668

7. Rulison EL, Kuczaj I, Pang G, Hickling GJ, Tsao JI, Ginsberg HS. Flagging versus dragging as sampling methods for nymphal Ixodes scapularis (Acari: Ixodidae). J Vector Ecol [Internet]. 2013 Jun;38(1):163–7. Available from: http://doi.wiley.com/10.1111/j.1948-7134.2013.12022.x pmid:23701621

8. Fedorova N, Kleinjan JE, James D, Hui LT, Peeters H, Lane RS. Remarkable diversity of tick or mammalian-associated Borreliae in the metropolitan San Francisco Bay Area, California. Ticks Tick Borne Dis [Internet]. 2014 Oct;5(6):951–61. Available from: http://dx.doi.org/10.1016/j.ttbdis.2014.07.015 pmid:25129859

9. Salkeld DJ, Nieto NC, Carbajales-Dale P, Carbajales-Dale M, Cinkovich SS, Lambin EF. Disease risk & landscape attributes of tick-borne Borrelia pathogens in the San Francisco Bay Area, California. Stevenson B, editor. PLoS One [Internet]. 2015 Aug 19;10(8):e0134812. Available from: https://dx.plos.org/10.1371/journal.pone.0134812 pmid:26288371

10. Diuk-Wasser MA, Vourc’h G, Cislo P, Hoen AG, Melton F, Hamer SA, et al. Field and climate-based model for predicting the density of host-seeking nymphal Ixodes scapularis, an important vector of tick-borne disease agents in the eastern United States. Glob Ecol Biogeogr [Internet]. 2010 May;22(6):428. Available from: http://doi.wiley.com/10.1111/j.1466-8238.2010.00526.x

11. Eisen RJ, Feirer S, Padgett KA, Hahn MB, Monaghan AJ, Kramer VL, et al. Modeling climate suitability of the western blacklegged tick in California. J Med Entomol [Internet]. 2018 Apr 25;55(5):1133–42. Available from: https://academic.oup.com/jme/advance-article/doi/10.1093/jme/tjy060/4969332 pmid:29697837

12. Gatewood AG, Liebman KA, Vourc’h G, Bunikis J, Hamer SA, Cortinas R, et al. Climate and tick seasonality are predictors of Borrelia burgdorferi genotype distribution. Appl Environ Microbiol [Internet]. 2009 Apr 15;75(8):2476–83. Available from: http://www.ncbi.nlm.nih.gov/pubmed/19251900 pmid:19251900

13. Pepin KM, Eisen RJ, Mead PS, Piesman J, Fish D, Hoen AG, et al. Geographic variation in the relationship between human Lyme disease incidence and density of infected host-seeking Ixodes scapularis nymphs in the Eastern United States. Am J Trop Med Hyg [Internet]. 2012 Jun 1;86(6):1062–71. Available from: http://www.ajtmh.org/content/journals/10.4269/ajtmh.2012.11-0630 pmid:22665620

14. Salkeld DJ, Antolin MF. ecological fallacy and aggregated data: a case study of fried chicken restaurants, obesity and Lyme disease. Ecohealth [Internet]. 2020; Available from: https://doi.org/10.1007/s10393-020-01472-1 pmid:32026056

15. Raghavan RK, Peterson AT, Cobos ME, Ganta R, Foley D. Current and future distribution of the lone star tick, Amblyomma americanum (L.) (Acari: Ixodidae) in North America. Munderloh UG, editor. PLoS One [Internet]. 2019 Jan 2;14(1):e0209082. Available from: https://dx.plos.org/10.1371/journal.pone.0209082 pmid:30601855

16. Springer YP, Eisen L, Beati L, James AM, Eisen RJ. Spatial distribution of counties in the Continental United States with records of occurrence of Amblyomma americanum (Ixodida: Ixodidae). J Med Entomol [Internet]. 2014 Mar 1;51(2):342–51. Available from: http://www.ncbi.nlm.nih.gov/pubmed/24724282 pmid:24724282

17. Springer YP, Jarnevich CS, Monaghan AJ, Eisen RJ, Barnett DT. Modeling the present and future geographic distribution of the lone star tick, Amblyomma americanum (Ixodida: Ixodidae), in the Continental United States. Am J Trop Med Hyg [Internet]. 2015 Oct 7;93(4):875–90. Available from: http://www.ajtmh.org/content/journals/10.4269/ajtmh.15-0330 pmid:26217042

18. Raghavan RK, Goodin DG, Hanzlicek GA, Zolnerowich G, Dryden MW, Anderson GA, et al. Maximum entropy-based ecological niche model and bio-climatic determinants of lone star tick (Amblyomma americanum) niche. Vector-Borne Zoonotic Dis [Internet]. 2016 Mar;16(3):205–11. Available from: http://online.liebertpub.com/doi/10.1089/vbz.2015.1837 pmid:26824880

19. Kessler WH, Blackburn JK, Sayler KA, Glass GE. Estimating the geographic distribution of host-seeking adult Amblyomma americanum (Acari: Ixodidae) in Florida. J Med Entomol [Internet]. 2019 Jan 8;56(1):55–64. Available from: https://academic.oup.com/jme/article/56/1/55/5087684 pmid:30169746

20. James AM, Burdett C, McCool MJ, Fox A, Riggs P. The geographic distribution and ecological preferences of the American dog tick, Dermacentor variabilis (Say), in the U.S.A. Med Vet Entomol [Internet]. 2015 Jun;29(2):178–88. Available from: http://doi.wiley.com/10.1111/mve.12099 pmid:25684582

21. Raghavan RK, Barker SC, Cobos ME, Barker D, Teo EJM, Foley DH, et al. Potential spatial distribution of the newly introduced long-horned tick, Haemaphysalis longicornis in North America. Sci Rep [Internet]. 2019 Dec 24;9(1):498. Available from: http://www.nature.com/articles/s41598-018-37205-2 pmid:30679711

22. Rochlin I. Modeling the asian longhorned tick (Acari: Ixodidae) suitable habitat in North America. J Med Entomol [Internet]. 2019 Feb 25;56(2):384–91. Available from: https://academic.oup.com/jme/article/56/2/384/5213063 pmid:30544234

23. MacDonald AJ, O’Neill C, Yoshimizu MH, Padgett KA, Larsen AE. Tracking seasonal activity of the western blacklegged tick across California. Mosher B, editor. J Appl Ecol [Internet]. 2019 Nov 5;56(11):2562–73. Available from: https://onlinelibrary.wiley.com/doi/abs/10.1111/1365-2664.13490

24. Eisen RJ, Eisen L, Beard CB. County-scale distribution of Ixodes scapularis and Ixodes pacificus (Acari: Ixodidae) in the Continental United States. J Med Entomol [Internet]. 2016 Mar;53(2):349–86. Available from: http://www.ncbi.nlm.nih.gov/pubmed/26783367 pmid:26783367

25. Hahn MB, Jarnevich CS, Monaghan AJ, Eisen RJ. Modeling the geographic distribution of Ixodes scapularis and Ixodes pacificus (Acari: Ixodidae) in the Contiguous United States. J Med Entomol [Internet]. 2016 Sep;53(5):1176–91. Available from: https://academic.oup.com/jme/article-lookup/doi/10.1093/jme/tjw076 pmid:27282813

26. Dennis DT, Nekomoto TS, Victor JC, Paul WS, Piesman J. Forum: reported distribution of Ixodes scapularis and Ixodes pacificus (Acari: Ixodidae) in the United States. J Med Entomol [Internet]. 1998 Sep 1;35(5):629–38. Available from: http://www.ncbi.nlm.nih.gov/pubmed/9775584 pmid:9775584

27. Peterson AT, Raghavan RK. The geographic distribution of ixodes scapularis (Acari: Ixodidae) revisited: the importance of assumptions about error balance. J Med Entomol. 2017;54(4):1080–4. pmid:28591858

28. Hahn MB, Jarnevich CS, Monaghan AJ, Eisen RJ. Response: The Geographic distribution of Ixodes scapularis (Acari: Ixodidae) revisited: the importance of assumptions about error balance. J Med Entomol [Internet]. 2017 Sep 1;54(5):1104–6. Available from: https://academic.oup.com/jme/article/54/5/1104/3861060 pmid:28874013

29. Johnson TL, Bjork JKH, Neitzel DF, Dorr FM, Schiffman EK, Eisen RJ. habitat suitability model for the distribution of Ixodes scapularis (Acari: Ixodidae) in Minnesota. J Med Entomol [Internet]. 2016 May;53(3):598–606. Available from: https://academic.oup.com/jme/article-lookup/doi/10.1093/jme/tjw008 pmid:27026161

30. Soucy J-PR, Slatculescu AM, Nyiraneza C, Ogden NH, Leighton PA, Kerr JT, et al. High-resolution ecological niche modeling of Ixodes scapularis ticks based on passive surveillance data at the northern frontier of Lyme disease emergence in North America. Vector-Borne Zoonotic Dis [Internet]. 2018 May;18(5):235–42. Available from: http://www.liebertpub.com/doi/10.1089/vbz.2017.2234 pmid:29565748

31. Sage KM, Johnson TL, Teglas MB, Nieto NC, Schwan TG. Ecological niche modeling and distribution of Ornithodoros hermsi associated with tick-borne relapsing fever in western North America. Coburn J, editor. PLoS Negl Trop Dis [Internet]. 2017 Oct 30;11(10):e0006047. Available from: https://dx.plos.org/10.1371/journal.pntd.0006047 pmid:29084219

32. Nieto NC, Porter WT, Wachara JC, Lowrey TJ, Martin L, Motyka PJ, et al. Using citizen science to describe the prevalence and distribution of tick bite and exposure to tick-borne diseases in the United States. Stevenson B, editor. PLoS One [Internet]. 2018 Jul 12;13(7):e0199644. Available from: http://dx.plos.org/10.1371/journal.pone.0199644 pmid:30001350

33. Xu G, Mather TN, Hollingsworth CS, Rich SM. Passive surveillance of Ixodes scapularis (say), their biting activity, and associated pathogens in Massachusetts. Vector-Borne Zoonotic Dis [Internet]. 2016;16(8):520–7. Available from: http://online.liebertpub.com/doi/10.1089/vbz.2015.1912 pmid:27248292

34. Xu G, Pearson P, Dykstra E, Andrews ES, Rich SM. Human-biting Ixodes ticks and pathogen prevalence from California, Oregon, and Washington. Vector-Borne Zoonotic Dis [Internet]. 2019 Feb;19(2):106–14. Available from: https://www.liebertpub.com/doi/10.1089/vbz.2018.2323 pmid:30207871

35. Ripoche M, Gasmi S, Adam-Poupart A, Koffi JK, Lindsay LR, Ludwig A, et al. Passive tick surveillance provides an accurate early signal of emerging Lyme disease risk and human cases in Southern Canada. J Med Entomol [Internet]. 2018 Jun 28;55(4):1016–26. Available from: https://academic.oup.com/jme/advance-article/doi/10.1093/jme/tjy030/4911759 pmid:29522180

36. Rand PW, Lacombe EH, Dearborn R, Cahill B, Elias S, Lubelczyk CB, et al. Passive surveillance in Maine, an area emergent for tick-borne diseases. J Med Entomol [Internet]. 2007 Nov 1;44(6):1118–29. Available from: http://www.mmcri.org/deptPages/lyme/downloads/RAND2007.pdf pmid:18047214

37. Salkeld DJ, Porter WT, Loh SM, Nieto NC. Time of year and outdoor recreation affect human exposure to ticks in California, United States. Ticks Tick Borne Dis [Internet]. 2019;(June):0–1. Available from: pmid:31201125

38. Porter WT, Motyka PJ, Wachara J, Barrand ZA, Hmood Z, McLaughlin M, et al. Citizen science informs human-tick exposure in the Northeastern United States. Int J Health Geogr [Internet]. 2019 Dec 7;18(1):9. Available from: https://doi.org/10.1186/s12942-019-0173-0 pmid:31064416

39. Laaksonen M, Sajanti E, Sormunen JJ, Penttinen R, Hänninen J, Ruohomäki K, et al. Crowdsourcing-based nationwide tick collection reveals the distribution of Ixodes ricinus and I. persulcatus and associated pathogens in Finland. Emerg Microbes Infect [Internet]. 2017 Jan 15;6(1):1–7. Available from: https://www.tandfonline.com/doi/full/10.1038/emi.2017.17 pmid:28487561

40. Little EAH, Anderson JF, Stafford KC, Eisen L, Eisen RJ, Molaei G. Predicting spatiotemporal patterns of Lyme disease incidence from passively collected surveillance data for Borrelia burgdorferi sensu lato-infected Ixodes scapularis ticks. Ticks Tick Borne Dis [Internet]. 2019 Aug;10(5):970–80. Available from: https://doi.org/10.1016/j.ttbdis.2019.04.010 pmid:31101553

41. Zimmermann NE, Edwards TC, Graham CH, Pearman PB, Svenning J-C. New trends in species distribution modelling. Ecography (Cop) [Internet]. 2010 Dec;33(6):985–9. Available from: http://doi.wiley.com/10.1111/j.1600-0587.2010.06953.x

42. Eisen RJ, Kugeler KJ, Eisen L, Beard CB, Paddock CD. Tick-borne zoonoses in the United States: persistent and emerging threats to human health. ILAR J [Internet]. 2017 Dec 15;58(3):319–35. Available from: https://academic.oup.com/ilarjournal/article-lookup/doi/10.1093/ilar/ilx005 pmid:28369515

43. Foley JE, Foley P, Brown RN, Lane RS, Dumlers JS, Madigan JE. Ecology of Anaplasma phagocytophilum and Borrelia burgdorferi in the western United States. J Vector Ecol [Internet]. 2004 Jun;29(1):41–50. Available from: http://www.ncbi.nlm.nih.gov/pubmed/15266739 pmid:15266739

44. Lane RS, Brown RN, Piesman J, Peavey CA. Vector competence of Ixodes pacificus and Dermacentor occidentalis (Acari: Ixodidae) for various isolates of Lyme disease spirochetes. J Med Entomol [Internet]. 1994 May 1;31(3):417–24. Available from: https://academic.oup.com/jme/article-lookup/doi/10.1093/jmedent/31.3.417 pmid:8057316

45. Salkeld DJ, Nieto NC, Bonilla DL, Yoshimizu MH, Padgett KA. Borrelia miyamotoi infections in small mammals, California, USA. Emerg Infect Dis [Internet]. 2018 Dec;24(12):2356–9. Available from: http://wwwnc.cdc.gov/eid/article/24/12/17-1632_article.htm pmid:30457525

46. Salkeld DJ, Cinkovich S, Nieto NC. Tick-borne pathogens in Northwestern California, USA. Emerg Infect Dis [Internet]. 2014 Mar;20(3):493–4. Available from: http://wwwnc.cdc.gov/eid/article/20/3/13-0668_article.htm pmid:24565119

47. Nieto NC, Salkeld DJ. epidemiology and genetic diversity of anaplasma phagocytophilum in the San Francisco Bay Area, California. Am J Trop Med Hyg [Internet]. 2016 Jul 6;95(1):50–4. Available from: http://www.ajtmh.org/content/journals/10.4269/ajtmh.15-0707 pmid:27139447

48. Thornton PE, Thornton MM, Mayer BW, Wei Y, Devarakonda R, Vose RS, et al. Daymet: daily surface weather data on a 1-km grid for North America, version 3 [Internet]. ORNL Distributed Active Archive Center; 2017. Available from: https://daac.ornl.gov/cgi-bin/dsviewer.pl?ds_id=1328

49. R Core Team. R: a language and environment for statistical computing [Internet]. Vienna, Austria; 2018. Available from: https://www.r-project.org/

50. Bocinsky RK. FedData: functions to automate downloading geospatial data available from several federated data sources [Internet]. 2016. Available from: http://cran.r-project.org/package=FedData

51. Hijmans R., Phillips S, Leathwick J, Elith J. dismo: species distribution modeling. 2017.

52. U.S. Geological Survey. 100-Meter Resolution Elevation of the Conterminous United States [Internet]. National Atlas of the United States. 2012. Available from: http://nationalatlas.gov/atlasftp-1m.html

53. Yang L, Jin S, Danielson P, Homer C, Gass L, Bender SM, et al. A new generation of the United States National Land Cover Database: requirements, research priorities, design, and implementation strategies. ISPRS J Photogramm Remote Sens [Internet]. 2018 Dec;146(August):108–23. Available from: https://doi.org/10.1016/j.isprsjprs.2018.09.006

54. Muscarella R, Galante PJ, Soley-Guardia M, Boria RA, Kass JM, Uriarte M, et al. ENMeval: An R package for conducting spatially independent evaluations and estimating optimal model complexity for Maxent ecological niche models. McPherson J, editor. Methods Ecol Evol [Internet]. 2014 Nov 18;5(11):1198–205. Available from: http://doi.wiley.com/10.1111/2041-210X.12261

55. Phillips SJ, Dudík M, Schapire RE. Maxent software for modeling species niches and distributions (Version 3.4.1) [Internet]. [cited 2020 May 25]. Available from: http://biodiversityinformatics.amnh.org/open_source/maxent/

56. Phillips SJ, Anderson RP, Schapire RE. Maximum entropy modeling of species geographic distributions. Ecol Modell [Internet]. 2006 Jan;190(3–4):231–59. Available from: https://linkinghub.elsevier.com/retrieve/pii/S030438000500267X

57. Yackulic CB, Chandler R, Zipkin EF, Royle JA, Nichols JD, Campbell Grant EH, et al. Presence-only modelling using MAXENT: when can we trust the inferences? Methods Ecol Evol. 2013;4(3):236–43.

58. Jueterbock A, Smolina I, Coyer JA, Hoarau G. The fate of the Arctic seaweed Fucus distichus under climate change: an ecological niche modeling approach. Ecol Evol [Internet]. 2016 Mar;6(6):1712–24. Available from: pmid:27087933

59. Taylor KE, Stouffer RJ, Meehl GA. An overview of CMIP5 and the experiment design. Bull Am Meteorol Soc. 2012;93(4):485–98.

60. Hijmans RJ. raster: geographic data analysis and modeling [Internet]. 2019. Available from: https://cran.r-project.org/package=raster

61. Zhang G, Zhu A-X. The representativeness and spatial bias of volunteered geographic information: a review. Ann GIS [Internet]. 2018 Jul 3;24(3):151–62. Available from: https://www.tandfonline.com/doi/full/10.1080/19475683.2018.1501607

62. Doggett JS, Kohlhepp S, Gresbrink R, Metz P, Gleaves C, Gilbert D. Lyme disease in Oregon. J Clin Microbiol [Internet]. 2008 Jun 1;46(6):2115–8. Available from: https://jcm.asm.org/content/46/6/2115 pmid:18448697

63. Dykstra EA, Oltean HN, Kangiser D, Marsden-Haug N, Rich SM, Xu G, et al. Ecology and epidemiology of tickborne pathogens, Washington, USA, 2011–2016. Emerg Infect Dis [Internet]. 2020 Apr;26(4):648–57. Available from: http://wwwnc.cdc.gov/eid/article/26/4/19-1382_article.htm pmid:32187009

64. Tufts DM, VanAcker MC, Fernandez MP, DeNicola A, Egizi A, Diuk-Wasser MA. Distribution, host-seeking phenology, and host and habitat associations of Haemaphysalis longicornis ticks, Staten Island, New York, USA. Emerg Infect Dis [Internet]. 2019 Apr;25(4):792–6. Available from: http://wwwnc.cdc.gov/eid/article/25/4/18-1541_article.htm pmid:30681071

65. Newman EA, Eisen L, Eisen RJ, Fedorova N, Hasty JM, Vaughn C, et al. Borrelia burgdorferi sensu lato spirochetes in wild birds in Northwestern California: associations with ecological factors, bird behavior and tick infestation. Stevenson B, editor. PLoS One [Internet]. 2015 Feb 25;10(2):e0118146. Available from: https://dx.plos.org/10.1371/journal.pone.0118146 pmid:25714376

66. Loss SR, Noden BH, Hamer GL, Hamer SA. A quantitative synthesis of the role of birds in carrying ticks and tick-borne pathogens in North America. Oecologia [Internet]. 2016 Dec 26;182(4):947–59. Available from: http://link.springer.com/10.1007/s00442-016-3731-1 pmid:27670413

67. Ogden NH, Lindsay LR, Hanincová K, Barker IK, Bigras-Poulin M, Charron DF, et al. Role of migratory birds in introduction and range expansion of Ixodes scapularis ticks and of Borrelia burgdorferi and Anaplasma phagocytophilum in Canada. Appl Environ Microbiol [Internet]. 2008 Mar 15;74(6):1780–90. Available from: https://aem.asm.org/content/74/6/1780 pmid:18245258

68. Fernandez MP, Bron GM, Kache PA, Larson SR, Maus A, Gustafson D, et al. Usability and feasibility of a smartphone app to assess human behavioral factors associated with tick exposure (the Tick app): Quantitative and qualitative study. JMIR mHealth uHealth. 2019;7. pmid:31651409

69. Kopsco HL, Xu G, Luo C-Y, Rich SM, Mather TN. Crowdsourced Photographs as an Effective Method for Large-Scale Passive Tick Surveillance. J Med Entomol. 2020; 1–9. pmid:31228205