Climate change influences on the geographic distributional potential of the spotted fever vectors Amblyomma maculatum and Dermacentor andersoni

Methods

Data preparation

We obtained totals of 255 and 586 occurrence points for A. maculatum and D. andersoni, respectively. Those data were obtained from various sources: Global Biodiversity Information Facility (GBIF; http://www.gbif.org), VectorMap (http://vectormap.si.edu/), and BISON (https://bison.usgs.gov) (sources summarized in Supplemental Information 1). We followed Cobos et al. (2018) in cleaning the data to remove errors that clearly fall outside of the known geographic distribution of the species, duplicate records, and localities with missing or meaningless coordinates such as zero degrees latitude and zero degrees longitude (0°N, 0°E), or georeferencing errors (records in ocean and far from coast). We used the spTthin R package to reduce the data spatially based on a 50 km distance filter for several reasons: based on precision of the occurrence points in the area, environmental heterogeneity that present in the area, and to avoid problems with autocorrelation (Aiello-Lammens et al., 2015). In the end, we had 93 and 82 occurrence points for A. maculatum and D. andersoni, respectively. We divided the final occurrence data randomly into two sets: 50% for model calibration and 50% for evaluation steps involved in model calibration. For producing final models, we used the entire cleaned occurrence points.

Delineate the calibration area

The accessible area (M) is the set of places to which the species has had access over relevant time periods, and depends on the dispersal of the species from populations (Barve et al., 2011). Since the movement of tick species is associated with the movement of host species, we assumed ample dispersal abilities for the ticks (Nadolny & Gaff, 2018; Sonenshine, 2018). As such, we created 200 km buffer areas around the known occurrence points for each species (Fig. 1).

Environmental variables

For the current time, bioclimatic variables were downloaded from WorldClim version 1.4, at 10′ spatial resolution (Hijmans et al., 2005) (available at http://www.worldclim.org/). We removed variables 8, 9, 18, and 19 because of known spatial artefacts (Escobar, 2020). The 15 remaining variables were masked to the calibration area (M) for each species. We then used principal component analysis (PCA) to reduce dimensionality and multicollinearity among those variables. After having PCA results, we created 11 sets of environmental variables that represent all possible combinations of the first four principal components to test them with other parameter settings to choose best models during model calibration, following Cobos et al. (2019) (see below).

For future climatic conditions, we used five general circulation models (GCMs) under two representative concentration pathway scenarios (RCP 4.5, and RCP 8.5). Future climate data layers were downloaded from the Climate Change, Agriculture and Food Security (CCAFS) database at 10′ resolution (available at: http://www.ccafs-climate.org/data_spatial_downscaling). GCMs used were (1) National Center for Atmospheric Research (NCAR_CCSM4); (2) Met Office Hadley Centre (HadGEM2); (3) Model for Interdisciplinary Research on Climate (MIROC5); (4) Institut Pierre Simon Laplace (IPSL_CM5A); and (5) Russian Institute for Numerical Mathematics Climate Model Version 4 (INM_CM4). GCM choice was based on frequency of use in other such research applications, and on full availability of scenarios for both RCP scenarios.

Ecological niche modeling and model transfers

The combination of 11 sets of environmental variables, 15 feature classes (all combinations of linear = l, quadratic = q, product = p, hinge = h), and 17 regularization multiplier values (0.1 to 1 at intervals of 0.1, and 2 to 10 at intervals of 1) resulted in 2,805 candidate models for each species. We evaluated candidate models based on statistical significance (partial ROC, P ≤ 0.05; Peterson, Cobos & Jiménez-García, 2018), predictive performance (omission rates, <5%; Anderson, Lew & Peterson, 2003), and a criterion of minimum complexity (Akaike Information Criterion corrected for small sample sizes, AICc; Warren & Seifert, 2011). Specifically, we used differences between particular AICc values and the minimum values (ΔAICc < 2) to select best model parameter settings with which to produce final models.

Final models

For creating final models, we used the complete set of occurrences and the parameterizations selected during model calibration. We created 10 bootstrap replicates, and transferred the models across North America (Mexico, United States and Canada) under current and future scenarios. We calculated medians of all replicate medians from final predictions for each calibration area in which final models were produced to summarize model results. Then, we binarized models using a threshold of allowable omission error rate (E) of 5%, assuming that as a percentage of data may have included errors that misrepresented environments used by the species.

We calculated differences in suitable areas between current and the two future scenarios RCP (4.5, and 8.5). For representing changes of suitable areas, we used the agreement of changes (stable, gain, loss) among the five GCMs per RCP scenario. Simply, for each RCP scenario, we took all projections to future conditions based on distinct GCMs and compared against the current projection, and quantified the agreement of gain and loss of suitable areas, as well as the stability of suitable and unsuitable conditions. All modeling analysis steps were done in R 3.5.1 (R Core Team, 2018) using Maxent 3.4.1 (Phillips et al., 2017), implemented in the kuenm package (Cobos et al., 2019).

Uncertainty in model projections

We used the mobility-oriented parity metric (MOP, considering the nearest 5% of reference cloud) (Owens et al., 2013) to assess strict extrapolation risk. We also calculated variance arising from distinct sources (replicates, parameter settings, GCMs, and RCPs) in our model projections (Peterson, Papeş & Soberón , 2018). Both model variability and strict extrapolation were represented geographically following Owens et al. (2013) and Cobos et al. (2019), respectively.

Model calibration results

From among 2,805 candidate models for each of A. maculatum and D. andersoni, 2,728 and 2,554 were significantly better than random expectations, respectively (pROC test, p ≤ 0.05). Of these models, 2,129 and 761 met the omission rate criteria, (i.e., OR ≤ 0.05) respectively. Based on AICc, 55 and two models were selected as best models for A. maculatum and D. andersoni, respectively. For A. maculatum, models performed better with the variables in set 1 (PC1, PC2, PC3, PC4), set 2 (PC1, PC2, PC3), set 3 (PC1, PC2, PC4), and set 6 (PC1, PC2), whereas for D. andersoni variables in sets 4 (PC1, PC3, PC4) and 7 (PC1, PC3).

Current and future potential distribution

Amblyomma maculatum

Current model predictions for A. maculatum showed suitable areas across the southern United States (Florida, Georgia, South and North Carolina, Virginia, West Virginia, Maryland, Delaware, Kentucky, Tennessee, Arkansas, Alabama, Mississippi, Louisiana, Oklahoma, and Texas), and in the Midwest (Missouri; eastern Kansas; southern Illinois, Indiana, and Ohio), and restricted areas of northeastern states (New Jersey and Pennsylvania). Suitable areas extend to include areas in western states (Arizona, California, Oregon, and Washington), although those areas are not likely accessible to the species (Fig. 2). Our models also showed suitable areas for the species across parts of eastern, western, and southern Mexico (Quintana Roo) (Fig. 2).

Future model transfers showed stable suitable areas (i.e., suitable in current time and in the future time) across the South, Midwest, and the Northeast, in the United States (Fig. 2). Areas of range reduction (loss) were in restricted areas in Kansas, Oklahoma, and Texas. Range expansion (gain) was anticipated in the northeastern (Pennsylvania, New York, Connecticut, Rhode Island, Massachusetts, Vermont, New Hampshire, Maine) and midwestern states (Kansas, Missouri, Nebraska, Iowa, Illinois, Indiana, Ohio, Michigan, Wisconsin) (Fig. 2). In general, we noted greater agreement among models in terms of losses and gains in the RCP 8.5 scenario compared to RCP 4.5.

Dermacentor andersoni

Current-time range predictions for D. andersoni showed broad suitable areas across Washington, Idaho, Oregon, California, Montana, Nevada, Utah, Wyoming, and Colorado, in cases where this species is known to occur in the United States. Climatically suitable areas extended across the Midwest, and Northeast and in some southeastern states (Fig. 3), although these areas are not likely accessible to the species. Currently suitable areas were also observed in parts of central and western Canada (British Columbia, Alberta, Saskatchewan, and restricted areas in Manitoba) (Fig. 3).

Future model transfers showed stable suitable areas across the states listed above, with some degree of reduction in suitable areas in the western states including much of Washington, Oregon, California, Nevada, Arizona, New Mexico, and Utah, and restricted areas in Colorado, Idaho, and Montana (Fig. 3). Predictions for the two RCP scenarios showed closely similar patterns of range stability, expansion, and loss, with more agreement among models in the RCP 8.5 scenario (Fig. 3).

Model uncertainty

MOP results for A. maculatum showed that strict extrapolative areas among future scenarios were concentrated in northern parts in North America, particularly in Canada, and in some restricted areas of the United States and southern Mexico (Fig. 2). Model variability results showed almost no variation coming from replicates and RCPs, but high contribution to variation from GCMs and parameter choice (Fig. S3).

In D. andersoni, we noted high agreement of strictly extrapolative areas in both southern and northern North America, and in lesser degree in the eastern United States and Canada (Fig. 3). High model variability came mainly from parameter choice in the eastern United States; we noted low variation deriving from GCMs, RCPs, and replicates (Fig. S4).