Abstract

Abstract

The impacts of climate change on human and wildlife health remain poorly understood. Carlson and colleagues develop a new method for reconstructing these impacts, and show that since 1950, environmental conditions in the western United States have become more favorable for plague (Yersinia pestis), including both its maintenance in wild mammals and spillover risk for humans.

Open in a separate window

1. INTRODUCTION

The distribution and burden of infectious diseases will be entirely reshaped by global environmental change. Scientific consensus suggests that over the next century, the combined effect of climate change, land degradation and transformation, and increasing human–wildlife contact will bring about a massive increase in the spillover of pathogens that originate in wildlife (zoonotic diseases) (Carlson et al., ²⁰²⁰; Estrada‐Pena et al., ²⁰¹⁴) and the burden of infections transmitted by arthropods (vector‐borne diseases) (Caminade et al., ²⁰¹⁹; Rocklöv & Dubrow, ²⁰²⁰). While there are substantial research efforts working to project these future changes (Carlson et al., ²⁰²⁰; Kraemer et al., ²⁰¹⁹; Ryan et al., ²⁰¹⁹), the impacts of current environmental change on infectious disease burden in the world today are underexplored. Based on current evidence, land use change is the best‐supported leading driver of zoonotic emergence (Gibb et al., ²⁰²⁰; Loh et al., ²⁰¹⁵); much less is known about climate change impacts to date. This is due, in large part, to methodological limitations: the “detection and attribution” methods that are best suited to this problem require substantial data on disease prevalence or incidence over extensive periods, as well as complicated model designs (e.g., counterfactual climate scenarios without climate change) (Carleton & Hsiang, ²⁰¹⁶; Ebi et al., ²⁰¹⁷).

Instead, many projections of climate change impacts rely on ecological niche models (also known as species distribution models), a set of regression and machine learning approaches that relate climate to the geographic range of a species (Escobar & Craft, ²⁰¹⁶; Hay et al., ²⁰¹³). Usually, these approaches are an oversimplification of reality, especially for pathogens: for example, a map of anthrax (Bacillus anthracis) may classify west Texas as an endemic zone, even though the system is characterized by epizootics that are sometimes years apart (Blackburn et al., ²⁰⁰⁷). Ecological niche models are therefore an imperfect tool for exploring climate change impacts. These methods work well for mapping current distributions, for projecting single‐time‐slice distributions under future climates, and—in some recent work—for projecting continuous‐time change (Trisos et al., ²⁰²⁰). Retrospective work to reconstruct climate change impacts is much rarer, and is usually restricted to work that builds two species distribution models for contrasting time intervals and compares them (Dobrowski et al., ²⁰¹¹; Kulkarni et al., ²⁰¹⁶; Nakazawa et al., ²⁰⁰⁷). This approach is less than ideal, forcing researchers to violate the assumption that species’ geographic ranges are at equilibrium (Araújo et al., ²⁰⁰⁵; Gallien et al., ²⁰¹²); to aggregate data into somewhat arbitrary time periods; and to compare models trained on nonindependent but nonoverlapping datasets, which will generate different biological response curves simply because of model uncertainty. In this framework, it is also difficult to eliminate alternate hypotheses for why a species’ apparent distribution might change, like noise in the detection process or shifting abundance patterns.

Recently, a growing set of tools have tried to grapple with the temporal variability exhibited by the distribution of infectious diseases. Though most disease maps are treated as the long‐term average of temporally dynamic processes, time‐specific ecological niche modeling has been proposed as an alternative that captures the dynamic nature of transmission. Almost always, though, these methods have been implemented at the finest temporal scales: monthly (Peterson et al., ²⁰⁰⁵) or seasonal (Kaul et al., ²⁰¹⁸; Schmidt et al., ²⁰¹⁷). Yet, this approach has been mostly untested as a way of understanding disease distributions over multiple years—and, ideally, of contextualizing the impacts of environmental change over decades (but see Nakazawa et al. (²⁰⁰⁷)).

Few systems provide a better opportunity to test this approach than plague, a globally cosmopolitan zoonotic infection caused by the bacterium Yersinia pestis. The bacterium is typically maintained by rodent (and sometimes shrew) hosts and the fleas that live on them or in their nests (Anisimov et al., ²⁰⁰⁴; Krasnov et al., ²⁰⁰⁶; Mahmoudi et al., ²⁰²⁰). Suitable plague reservoir host species are often thought to exhibit a standing variation in their susceptibility to the bacterium, such that some hosts generate a high enough (and often lethal) bacteremia in the blood for the efficient transmission of the bacterium by flea vectors, while the host population as a whole survives (Anisimov et al., ²⁰⁰⁴; Lowell et al., ²⁰¹⁵; Nilsson et al., ²⁰²¹). Which local flea species are important for the spread and persistence of the pathogen depends on their typical host's abundance, and the taxonomic range of hosts that they feed upon; both of these factors are correlated with intrinsic variation in their vector competence (Krasnov et al., ²⁰⁰⁶), and have important implications for plague dynamics, including spillover risk for humans. While other transmission modes are sometimes relevant to both spillover and epidemic transmission, including consumption of infected meat (Malek et al., ²⁰¹⁶) and droplet transmission (Randremanana et al., ²⁰¹⁹), respectively, the relationship between rodents, fleas, and bacteria gives shape to most salient ecological questions.

The global distribution of plague has been far from stable over the past two centuries; the Third Pandemic (late 18th to mid–20th Century) (Bramanti et al., ²⁰¹⁹; Jones et al., ²⁰¹⁹; Neerinckx et al., ^2010a; Seal, ¹⁹⁶⁰; Xu et al., ²⁰¹⁴) in particular was responsible for the introduction of Y. pestis into many new regions that were environmentally suitable but otherwise uncolonized, particularly the Americas (Link, ¹⁹⁵⁵; Morelli et al., ²⁰¹⁰; Schneider et al., ²⁰¹⁴). In some of these regions, outbreaks have faded over time, while in others, plague foci have persisted and the pathogen has become endemic, maintained by a sylvatic cycle in rodent reservoirs and flea vectors (Schneider et al., ²⁰¹⁴). Rodent biodiversity hotspots may be particularly conducive to the formation of these reservoirs (Sun et al., ²⁰¹⁹), a possible case of biodiversity amplification effects (Halliday & Rohr, ²⁰¹⁹; Luis et al., ²⁰¹⁸), where the diversity of competent hosts allows a virulent pathogen to be maintained at more stable levels. Though underexplored, emerging evidence also suggests that plague may persist in the soil, possibly by acting as an endosymbiont with amoebas (Benavides‐Montaño & Vadyvaloo, ²⁰¹⁷), from which it sporadically can reinfect burrowing rodents (Boegler et al., ²⁰¹²). Soil conditions may therefore further constrain the distribution of plague reservoirs (Ayyadurai et al., ²⁰⁰⁸; Baltazard et al., ¹⁹⁶⁴; Eisen et al., ²⁰⁰⁸; Karimi, ¹⁹⁶³). Like other pathogens that can persist in the soil (Carlson et al., ²⁰¹⁹; Limmathurotsakul et al., ²⁰¹⁶), provisionary evidence suggests that plague may be limited by soil salinity (Barbieri et al., ²⁰²⁰; Malek et al., ²⁰¹⁷), soil organic carbon, and alkalinity (Neerinckx et al., ²⁰⁰⁸). Though these factors may have limited influence in the short‐term dynamics of plague in any one location, at continental scales, they could reasonably be expected to shape where plague foci have become established.

Both experimental and ecological analyses suggest that plague dynamics are also highly sensitive to climatic conditions. The disease's sensitivity to bioclimatic conditions has been documented throughout its life cycle, but is particularly pronounced on the arthropod level, where temperature (and to a lesser degree humidity) influence the rate at which various flea species move through their life cycle (Bacot, ¹⁹¹⁴). Flea species differ in their temperature sensitivities (Bacot, ¹⁹¹⁴), making the local composition of flea communities an important consideration. Temperature also directly influences biochemical aspects of the transmission efficiency of the plague bacterium, particularly when temperatures rise above 27°C (Eisen et al., ²⁰⁰⁹), presumably by negatively influencing the stability of the biofilm that the bacterium forms in fleas. Temperature also influences rodent populations, including through a mechanism generally referred to as a trophic cascade: climatic anomalies influence primary productivity, driving changes in rodent density, which in turn change the density and biting preferences of fleas (Parmenter et al., ¹⁹⁹⁹; Xu et al., ²⁰¹⁵). The combination of these environmental sensitivities, when playing out across the scale of ecosystems, can lead to widespread synchronicities in plague epizootic periods (Kausrud et al., ²⁰⁰⁷).

All of these lines of evidence suggest that plague should be broadly sensitive to environmental change, and that in systems where trends in plague occurrence have been tracked, an anthropogenic signal might be detectable. The United States is the perfect system to test this approach, as data in this region are particularly abundant; human case data goes back over a century, to plague's first introduction on the Pacific coast in early 20th Century (Adjemian et al., ²⁰⁰⁷; Ben Ari et al., ²⁰⁰⁸; Link, ¹⁹⁵⁵). Moreover, the US Department of Agriculture (USDA) has collected plague seropositivity data from wildlife for multiple decades through the USDA National Wildlife Disease Program (Bevins et al., ²⁰¹²). In total, these national datasets include more records than many global studies of pathogen distributions (Carlson et al., ²⁰¹⁹), making this system a data‐rich testing ground. In addition, these data also cover nearly a century of environmental change, a temporal scope that allows time‐specific ecological niche modeling to be implemented. This also allows us to revisit one of the only previous attempts at this approach, which compared models of plague risk in the western United States based on case data in three multiyear time slices (1965–1969, 1980–1984, and 1995–1999), and concluded that plague risk had expanded since 1950 and would continue to do so in future (Nakazawa et al., ²⁰⁰⁷).

In this study, we revisit this prediction by using two independent data streams (human cases and wildlife serology) in a machine learning model called Bayesian additive regression trees (BART) (Carlson, ²⁰²⁰) (see Methods for a detailed explanation). Section 1 Climatic reconstructions are readily available for the duration of our study (1950–2017), allowing us to use annual climate layers (including long‐term anomalies) corresponding to the year of each plague case. This pairing allows us to improve model precision relative to long‐term averages, to differentiate areas of ephemeral versus persistent risk, and to identify the fingerprint of environmental change in risk trends. We also test whether the distribution of plague in this region is responsive to rodent biodiversity or soil chemistry and macronutrients, offering detailed insights into the factors that maintain plague risk. Finally, we propose a new approach that harnesses BART with random intercepts (riBART) to account for historical variability in detection and sampling, allowing us to confidently identify the signal of changing environmental conditions in plague occurrence over time. In doing so, we propose the first extension of ecological niche modeling that nods toward the ultimate aim of detection and attribution of anthropogenic climate change impacts on the geographic distribution of infectious diseases.

2. METHODS

Despite recent interest in modeling the distribution of major infectious diseases (Hay et al., ²⁰¹³; Kraemer et al., ²⁰¹⁶), there is no definitive global map of plague reservoirs. All existing global plague maps have been derived from expert opinion (Stenseth et al., ²⁰⁰⁸); all modeled products so far have been produced for national or continental scales (see Table S2). Plague ecology is regionally variable enough that this patchwork approach has the advantage of being tailored to relevant local predictors. However, the mix of modeling methods, variables, and spatiotemporal scales makes it nearly impossible to compare these models and develop any consensus on the biological or geological factors that determine where plague reservoirs can exist, and where not. In this study, we adapt predictors that have previously worked in other similar work on plague, and develop novel models of spatiotemporal risk patterns in the western United States based on human and wildlife data spanning from 1950 to 2017.

2.2. Modeling

Dozens of statistical methods have been applied to species distribution modeling in the past few decades, with a wide range of performance (Norberg et al., ²⁰¹⁹). Over the past few years, classification and regression tree (CART) methods—including random forests and boosted regression trees—have become especially popular for mapping the geographic distribution of infectious diseases (Bhatt et al., ²⁰¹³; Carlson et al., ²⁰¹⁹; Hieronimo et al., ²⁰¹⁴; Pigott et al., ²⁰¹⁴; Richards et al., ²⁰²⁰; Shearer et al., ²⁰¹⁸). Here, we use a fairly new method, Bayesian additive regression trees (BARTs), implemented with the R package embarcadero as a species distribution modeling wrapper for the dbarts package (Carlson, ²⁰²⁰; Dorie, ²⁰²⁰). BART is a powerful new method with growing application in computer science, and often performs comparably to other CART methods like random forests and boosted regression trees (Chipman et al., ²⁰¹⁰). In the embarcadero implementation, BARTs have several unique features that make them a powerful tool for disease mapping, such as model‐free variable importance measures, and automated variable selection; posterior distributions on predictions, as a measure of uncertainty; posterior distributions on partial dependence plots; two‐dimensional and spatially projected partial dependence plots; and various extensions, including random intercept models.

Like other CART methods, BART makes predictions by splitting predictor variables with a set of nested decision rules (“trees”) that assign estimated values to terminal nodes (“leaves”). Predictions are generated based on a sum‐of‐trees model, where a set of $n$ trees with leaves $(T_1,M_1),...,(T_n,M_n)$ each make predictions $g(·)$ that are added together, for a total estimate:

$$\widehat{Y}=f(X)=\sum _{j=1}^{m}g(X;T_j,M_j)+\epsilon ;\epsilon :\mathcal{N}(0,{\sigma }^2)$$

For logistic classification problems (like species distribution modeling), BART uses a logit link function:

$$\widehat{Y}=f(X)=\mathrm{\Phi }[\sum _{j=1}^{m}g(X;T_j,M_j)]$$

where $\mathrm{\Phi }$ is the standard normal cumulative distribution. An initial set of $n$ trees is fit, and then altered in an MCMC process based on a set of random changes to the sum‐of‐trees model (e.g., new splits added, levels rearranged, or leaves pruned). An initial burn‐in period is discarded, and then a set of posterior draws ($f^{\ast }$) create the posterior distribution for $p(f|y)\equiv p(\mathrm{trees}|\mathrm{data})$.

BART is easily implemented out‐of‐the‐box, even with a full Bayesian MCMC component. Three priors control the ways decision trees change: the probability each variable is drawn for a split, the probability of splitting values tested, and the probability a tree stops at a certain depth. In the simplest form, the first two can be set as uniform distributions, while the latter is usually set as a negative power distribution; they can also be adjusted using a full cross‐validation approach. This is handled automatically in the dbarts package, for which embarcadero is a wrapper. More advanced implementations with complex prior design are sometimes appropriate; for example, a Dirichlet distribution on the variable importance prior can help identify informative predictors in high dimensionality datasets (dozens or hundreds of covariates). However, in our case, we had confidence all variables were biologically plausible based on expert opinion.

2.2.1. The base models

We ran two separate baseline models, the first using human data from 1950 to 2005, and the second using the wildlife data from 2000 to 2017. For the human model, we used the number of cases recorded each year in each county to generate a set of random georeferenced pseudopresence points. We then generated seven pseudoabsence points in each year to create a roughly balanced design, for a total of n = 430 pseudopresence points and n = 392 pseudoabsence points. For the wildlife model, we balanced the design by subsampling seronegative animals in equal number to seropositive ones, for a final n = 5,002 true presence points and n = 4,759 true absence points.

Both models were run with the full predictor set, followed by an automated variable set reduction procedure implemented in embarcadero that formalizes the recommendations of Chipman et al. (²⁰¹⁰). In BART, variable importance is “model free,” measured as the number of splitting rules involving a given variable (but incorporating no information on the proportional effect on the outcome variable, or proportional improvement in the model predictions). In models with fewer trees (small n), informative variables tend to be selected more often, while uninformative variables are selected rarely or drop out entirely. This property of BART establishes a rubric that can be used to identify an informative variable set. First, an initial model is fit with all variables 100 times each for six different settings of ensemble size (n = 10, 20, 50, 100, 150, and 200 trees). Plotting the average importance of variables at each level offers a qualitative diagnostic of how informative each predictor is. Next, an initial set of 200 models with n = 10 trees are run, and variable importance is recorded and averaged across models. Models are run again (200 times) without the least informative variable from the first fit, and this is performed iteratively until only three variables remain; the variable set with the lowest average model root mean square error (RMSE), and therefore highest accuracy on the training data, is selected. Finally, we plot variable importance (including standard deviations based on model permutations).

Final models were run with the reduced variable set, with recommended BART model settings (200 trees, 1000 posterior draws with a burn‐in of 100 draws) and hyperparameters (power = 2.0, base = 0.95 for the tree regularization prior, which limits tree depth). We then used the retune function in embarcadero to run a full cross‐validation panel on the three prior parameters. retune runs a full cross‐validation across the k hyperprior (values of 1, 2, and 3), the base parameter (0.75 to 0.95 in increments of 0.05), and the exponent parameter (1.5 to 2 in increments of 0.1), and returns the model with the parameter combination that generates the minimum root mean squared error.

For the wildlife model, the final variable set included: temperature mean, maximum, and anomaly; rodent richness; elevation; and five soil traits (calcium, sodium, iron, clay, and sand). The model validated well on training data (AUC = 0.836). For the human model, the final variable set included a similar subset: precipitation anomaly; temperature mean, maximum, and anomaly; rodent richness; elevation; and four soil traits (sodium, iron, clay, and calcium). The model also validated well on training data (AUC = 0.909).

2.2.2. Alternate formulations

As a final check of model performance, we ran a separate model with the same predictor sets that withheld the years 2000–2005 from both. On the test dataset for humans (n = 64), the model performed very well by the standards of external cross‐validation (AUC = 0.820); on the test data for wildlife (n = 796), the model also performed well (AUC = 0.775). This indicated that both models were performing adequately.

We also recognize that model design can have a substantial effect on machine learning performance, and the downstream biological inference made by using ecological niche models. Given that BART is a relatively new method, it has been comparatively underexplored in this regard, and so a standard panel of “best practices” has not yet been recommended in the literature. However, for transparency about model uncertainty and the influence of subjective decisions on model outputs, we produced four major alternate formulations. First, we produced models that included all variables, rather than using the variable set reduction procedure, for both the human data (Figures S20‐S22) and wildlife data (Figures S23‐S25). We additionally considered two alternate formulations of the wildlife model. In the first, we used pseudoabsences instead of the true absences available in the data (Figure S26). Though this increased model AUC (0.929), and allowed slightly different balancing of the data, it lead to visually apparent overfitting. Finally, we ran an alternate model only using the coyote data in the NWDP dataset, which also performed adequately (AUC = 0.826; Figure S27). Both models were ultimately not selected because they left available, biologically meaningful data unused, and both produced predictions that were slightly less congruous with the human model.

2.2.3. Prediction, delineating foci, and measuring change

Although the models were trained over different intervals, the continuous and standardized set of predictors allowed cross‐prediction over the entire extent of the study (1950–2017). For each layer of annual prediction, we thresholded suitability based on a model‐specific threshold chosen to maximize the true skill statistic on the test data. We mapped areas of “unstable foci” as any region with at least 1 year of suitability, and “stable foci” as any region suitable in every year over the 70‐year interval. This allowed us to compare long‐term spatial patterns between the two models.

2.2.4. Random effects models for interannual variation

Prevalence changes year‐to‐year in both the data and modeled landscapes, but detecting the signal of climate change in that fluctuation can be challenging. There are several reasons prevalence could vary across years: (1) incidence is stochastic but temporally autocorrelated; (2) normal climatic variability (e.g., the Pacific Decadal Oscillation) or other socioecological trends (e.g., rising human populations) might also contribute to interannual variation, including nonlinear trends over time; (3) anthropogenic climate change is directly driving changes in plague risk, or indirectly changing the ecology of the involved species; (4) sampling effort varies between years (for wildlife); or (5) detection rates could change between years, due to testing or surveillance. The last of these is particularly relevant as a possible confounder, given that wildlife diagnostics changed in 2011. A positive trend in plague risk might be generated by increased climatic suitability for plague tranmission, but could also be generated by a consistent increase in plague detection due to improved diagnostics and increased sampling effort, loosely collinear with warming temperatures on the scale of 20 to 70 years.

We propose a new method that uses machine learning approaches (i.e., ecological niche models) to detect the signal of environmental change while adjusting for confounders at a high level. The approach is loosely modeled off the ideas underlying econometric approaches to climate change detection and attribution, which usually use fixed effects panel regression to control for spatiotemporal confounders in climatic signal. By attributing as much variance as possible to spatial, temporal, and other confounders, and then identifying climatic signal in the remaining variance, these approaches can pinpoint the signal of environmental change with a high degree of confidence. So far, no analog to these approaches exists for ecological niche models. Only a handful of studies have even added temporal heterogeneity to ENMs; so far, we know of none that have also independently controlled for interannual variation in detection, sampling effort, or species prevalence.

A solution to temporal confounders is particularly needed in this study, given the challenges of the time‐specific approach. In default settings, BART predictions converge on observed prevalence, that is, $\widehat{Y}=E[Y]=P(Y=1)$. For this reason, we balance the presences and absences, so that the model is as close to $\widehat{Y}\sim \mathcal{N}(0.5,{\sigma }^2)$ as possible. This produces a unique challenge for time‐dependent modeling. Presences are distributed unevenly across years, and consequently, so is positivity. In the human model, this arises artificially, because pseudoabsences are generated evenly across years. We chose this approach to avoid overrepresenting years with more cases in the data, which would introduce an additional collinearity, but as a result relative prevalence varies substantially. This bias also affects the wildlife model more organically; although the number of points per year varies independently of test positivity (cor $=0.113$, $p=0.687$), because most sampling is passive, there is still a wide range in annual prevalence (28% in 2006 versus 80% in 2017, both years with several hundred records), with a net trend toward higher positivity over time. Because prevalence varies between years in both models, the resulting with environmental change could confound the detection of meaningful signals.

Inspired by the econometric approach, we propose a use case for the random intercept BART model (riBART), which has recently been proposed as an extension of the method for clustered outcomes. The approach adds a random intercept term to the model (separate from the tree‐fitting process) based on the identified $K$ clusters

$$\widehat{Y}=f(X)=\mathrm{\Phi }[\sum _{j=1}^{m}g(X_k;T_j,M_j)]+{\alpha }_k$$

where the random intercepts ${\alpha }_k(k\in 1:K)$ are normally distributed around zero (i.e., the $K$ groups are assumed to have normally distributed, independent additive effects on the outcome variable). The error structure of the random effects and the sum‐of‐trees model are assumed independent. Here, we propose that the model can be fit as usual with a random intercept for year, as a way of accounting for temporal heterogeneity as a possible confounder. The yearly random intercept absorbs most of the interannual variation in plague prevalence (i.e., the relative ratio of presences and absences in the data), such that residual variation in prevalence should be roughly constant across years. Identifying the climatic signal in these residual data, and then examining predicted prevalence (without random intercepts) based only on environmental change, allows more confident statements about how environmental change contributes to shifting disease risk.

We revisited the two main models, and used riBART to add an annual random intercept to our model for each year, which we refer to throughout as the “detection” models. Fitting climate–plague response curves after this detrending decouples the possible collinearity between climate trends and coarse interannual signal in the data, which may be driven by natural variation in prevalence or other confounders (e.g., the 2011 change in wildlife testing protocols). We fit both detection models with a random intercept for year, plotted the random effects, and predicted over the 70‐year interval without the random effect included. (All functionality to implement SDMs with riBART is available in embarcadero as an updated release.)

2.2.5. Detecting change over time

To estimate trends for change over time, we fit a linear slope through each pixel‐by‐year. Multiplying by 68 years, we were able to estimate total percent change in suitability since 1950 in a given pixel. We did not limit these to pixels with a significant trend, as any frequentist significance test iterated over millions of pixels would be mostly meaningless. We generated these maps for the two primary models and the two detection models (Figure ​(Figure3),3), as well as (in the supplement) for mean temperature and precipitation (Figures S18 and S19).

Open in a separate window

FIGURE 3

Total percent change in plague suitability, 1950 to present, in the human (top) and wildlife models (bottom), before (left) and after (right) adding random intercepts to control for interannual variation

3. RESULTS

3.1. The distribution of plague

We generated two primary models of plague over time. The first covered 9,761 animals sampled for plague (2000–2017), and performed well (training AUC = 0.836; Figure S1). The second covered a total of 430 human cases of plague (1950–2005), and performed very well (training AUC = 0.909; Figure S2). When both models were rerun with an overlapping “test period” of 2000 to 2005 withheld, they performed adequately, with the human model (AUC = 0.820) performing better than the wildlife model (AUC = 0.775). As both models performed well in temporal cross‐validation, we used both to make annual predictions from 1950 to 2017, and split predictions into binary presence or absence risk maps for each year using the true skill statistic.

Both models found that the majority of plague risk in the western United States is, as expected, found west of the 100th meridian (Figure ​(Figure1).1). The human model mostly predicts risk in the “Four Corners” region (Utah, Colorado, Arizona, and New Mexico), where that risk is relatively stable across years. In contrast, the wildlife model predicts risk fairly expansively west of 100°, including at much higher latitudes than the human risk model. Risk varies much more across years in this model, but several areas are predicted to be environmentally suitable across years from Montana to west Texas. The suitable areas identified by the wildlife model in the southwest are less uniform than the human model, likely reflecting a finer scale differentiation of risk. There are two main reasons the human model might discriminate less in this region: human cases may be reported in different locations than the site of initial spillover, and occurrence points were randomly resampled at the county level (as data have been previously de‐identified).

Open in a separate window

FIGURE 1

Suitability for plague across all years (1950–2017), for humans (left) and wildlife (right). Top panels give mean suitability across all years; bottom panels show areas identified as suitable in no years, at least one, or all 68 years

For the most part, we found that risk areas identified in the human model were a subset of the much broader predictions made by the wildlife model (Figure ​(Figure2),2), with three major exceptions. The human model identified much broader risk in southern Arizona and New Mexico, likely due to how the cases were randomized at county levels. The human model also predicted areas of risk further east, in regions like east Texas or Oklahoma where plague is not known to be endemic (and, in this regard, the wildlife model better captures the known distribution), but conditions may be broadly favorable. Finally, and most notably, the human model predicted plague risk throughout California, in places that have previously been identified as high risk (Holt et al., ²⁰⁰⁹). This likely reflects a deficit of data from Californian sources in our wildlife model, as state wildlife surveillance is curated independently. Together, our findings indicate the value of comprehensive surveillance, and the possibility that zoonotic reservoirs may be more expansive than areas of known spillover.

Open in a separate window

FIGURE 2

The wildlife model's predictions largely encompass the human model's predictions, except in southern Arizona and California (where predictions extend into other areas of suspected plague risk) and west Texas (too far east for plague reservoirs)

3.3. Detecting environmental impacts on change over time

We found a strong temporal trend in both climate and plague risk since 1950 (Figure ​(Figure3).3). To test whether this signal might be confounded by exogenous factors, we used a novel approach where we trained riBART for each year, and projected the model again without random effects over the 68‐year period (see Methods). The random intercepts identified interannual variation in prevalence (Figures S13 and S14), detrended the data, and allowed the models to identify climate signal minus the confounder without substantially changing overall predictions (Figure S15). Subsequently, we predicted how suitability changed using the same model without random intercepts; this allowed us to be confident that the changing suitability we identified in these “detection models” was the consequence of constant relationships between temperature, precipitation, and plague transmission.

Both detection models identified a meaningful signal of temporal variation. The random intercepts identified a signal of rising prevalence through the wildlife dataset, particularly increasing after 2011 when the diagnostics were changed (Figures S13 and S14). In the human model, we identified a much more subtle long‐term quadratic trend peaking in the 1980s, matching a pattern that has been previously attributed to climate cycles like the Pacific Decadal Oscillation (Ben Ari et al., ²⁰⁰⁸). Surprisingly, the “detection” models identified an even stronger pattern of change over time (Figure S16). In the wildlife model, suitability increased an average of 4.8%, and 4.9% in the detection model, with a much fatter tail to the distribution as well. In the human model, suitability increased an average of 1.7% from 1950 to 2017, and 2.1% in the detection model. In much of the region, we found that plague suitability increased by 30 to 40% over the entire interval. We found that suitability rose most substantially in the wildlife model at high elevations, while spillover risk increased more gently, and peaked at mid‐elevations (Figure ​(Figure4,4, Figure S17). Because the detection models only predict change across years based on the temperature and precipitation variables, we are confident that the increase in plague suitability is driven by the long‐term signal of warming (roughly ${0.8}^{\circ }$ in the region since 1950) and higher anomalous precipitation is responsible for these predicted changes (Figures S18 and S19). We conclude that, even with several confounding factors, environmental change since the 1950s may have helped plague reservoirs become established at higher elevations, and slightly increased the risk of spillover into human populations at mid‐elevations. Over the coming half‐century, previous work suggests that the shift toward higher elevations is likely to continue in this region (Holt et al., ²⁰⁰⁹).

Open in a separate window

FIGURE 4

Environmental suitability for plague has increased substantially at high elevations for wildlife; risk of spillover has increased mildly at mid‐elevations. Lines are given as generalized additive model smooth fits based on the detection models

4. DISCUSSION

Our study shows that human and wildlife data can be used in conjunction to map plague reservoirs and spillover risk in the United States, and to make meaningful inferences about ecological drivers of transmission. We found support for two major hypotheses: the biodiversity amplification effect (higher rodent diversity supports establishment of plague foci) and the trophic cascade hypothesis (anomalously warm and anomalously wet, cold years can increase plague prevalence and spillover risk, respectively, through ripple effects of ecosystem productivity on rodents and flea populations). Support for these patterns has increasingly been found across systems, and points to a view of plague risk where weather conditions (and their impact on flea vectors) in rodent biodiversity hotspots are the primary driver of transmission and spillover. Though the plague “niche” may largely transcend any individual host species (Maher et al., ²⁰¹⁰), in future, our understanding of these mechanisms might be further refined by exploring more granular variation among the species involved. Identifying “true reservoirs” is a nontrivial task in disease ecology (Becker et al., ²⁰²⁰; Viana et al., ²⁰¹⁴); Yersinia pestis can infect hundreds of rodent species (Mahmoudi et al., ²⁰²⁰), and some species once considered key to maintain endemic plague are now known to be spillover hosts from unknown reservoirs (Colman et al., ²⁰²¹; Danforth et al., ²⁰¹⁸). Rodent habitat use adds another dimension to the problem: synanthropic and “wild” reservoirs may have very different richness hotspots (Sun et al., ²⁰¹⁹), which may explain some differences between the geography of maintenance and spillover (as, alternately, could the geography of domestic cat and dog ownership (Campbell et al., ²⁰¹⁹)). These patterns are further complicated by potential variation among fleas, with over two dozen experimentally verified vectors in the United States alone (Eisen et al., ²⁰⁰⁹), which vary in geographic distribution, thermal sensitivity, and affinity for wildlife and human hosts. All of these nuances may point to promising directions for future data synthesis and modeling.

We further found strong evidence that the North American distribution of plague is heavily influenced by soil conditions. The global distributions of soil‐persistent bacteria like anthrax (Bacillus anthracis), tularemia (Francisella tularensis), and botulism (Clostridium botulinum) are known to be constrained by the biochemical properties of soil. Less is known about plague, which is not spore forming, and until recently was mostly thought to behave like a typical vector‐borne zoonosis. It may be that soil properties affect the suitability of burrows for higher flea densities, or determine host homeostasis for minerals that impact the virulence of the infection; plague foci might therefore fall in the narrow range of conditions that can harbor higher densities of fleas, but do not substantially increase the lethality of the infection. However, increasing evidence also suggests that the bacterium can persist in the soil, possibly through symbiotic relationships with amoebas, for weeks to months—and possibly even years (Baltazard et al., ¹⁹⁶⁴; Karimi, ¹⁹⁶³). These complexities underscore the importance of a One Health approach while studying the ecology of plague, which—like anthrax and many other bacterial pathogens—circulates easily among fleas, rodents, other wildlife, humans, and the environment as one interconnected system (Carlson et al., ²⁰¹⁸).

Developing a better understanding of plague in well‐studied systems like the American West will help develop a broader picture of its ecology. At present, all global maps of plague foci have been compiled from expert knowledge; modeled products in the English language are limited to the western United States, China, and Africa (Table S2). In part, this reflects the challenges of sharing, aggregating, and consolidating surveillance data. It may also reflect concerns about model transferability, given that the complex multispecies dynamics of plague reservoirs differ greatly across ecosystems and continents. However, other pathogens with regional host communities and complex environmental persistence have been globally mapped through multinational coordination (Carlson et al., ²⁰¹⁹), and the same synthesis is possible for plague. In the more immediate term, our model also strongly suggests that wildlife reservoirs extend to the US borders with both Mexico and Canada, and could plausibly extend beyond them (recently confirmed for the northern border (Liccioli et al., ²⁰²⁰), but the official World Health Organization map of plague (last updated 2016) includes no reservoirs in Mexico or Canada. Collaborating with national surveillance infrastructure in both countries may help resolve the boundaries of plague transmission more clearly, and reveal foci currently overlooked by global monitoring efforts.

Beyond plague, our study highlights the opportunity for medical geographers to develop new methods that are suited to a rapidly changing world. Here, we proposed two methodological advances that build on existing best practices in infectious disease mapping. First, time‐specific covariates allowed us to train machine learning models on nearly a century's worth of data, improving precision compared to coarsely averaged predictors, and capturing the effects of environmental change. If this approach is integrated with others at finer temporal scales, such as those that consider seasonal aspects of transmission or spillover risk (Kaul et al., ²⁰¹⁸; Schmidt et al., ²⁰¹⁷), this could begin to set the foundation for an early warning system. Second, the use of random intercepts to remove data and detection biases, such as the serology method change in our data sample in 2011, is an important step toward testing climate change impacts using continuous‐time data (similar to how econometric approaches resolve these problems in similar spatiotemporal analyses [Carleton & Hsiang, ²⁰¹⁶]). We propose that when this approach can be taken, it may be used as a first principles method for detecting the signal of environmental change in species’ habitats. This could be a particularly important step toward synthesizing the impacts of climate change on the shifting presence and absence of disease data, especially in cases where prevalence and incidence data are lacking and panel regression approaches cannot be applied. However, this work will still need to be followed by proper “attribution” work that compares predicted patterns to counterfactual scenarios without climate change; at present, all we can conclude with certainty is that weather conditions have changed in a way that trends favorably for plague risk.

Our study also points to a number of gaps in our understanding of environmental change (and consequently, potential methodological limitations). The PRISM data offer a fairly comprehensive view of the recent climate in the United States, and allowed us to identify the role of temperature and precipitation in plague transmission. However, we held both soil and rodent predictor variables constant, and neither are stationary in reality. Soil has changed over the last century due to a combination of climate change and land use change, and unfortunately time‐specific covariates are unavailable; in many cases, our soil layers had to be generated custom to this study, and for the rest of the world these data are even more sparse. Similarly, evidence is strong that most terrestrial species have responded to recent climate change by undergoing range shifts, especially along elevational gradients. If rodents have undergone range shifts, they may have encountered novel vector communities, and the relationship between richness and transmission could change. Similarly, if elevation acts in our models as a proxy for specific rodent‐flea assemblages, range shifts could decouple the observed relationships between elevation and transmission. As other studies have pointed out, these challenges highlight the need to begin integrating zoonotic surveillance and biodiversity monitoring (Carlson et al., ²⁰²⁰).

CONFLICT OF INTEREST

The authors declare no conflicts of interest.

AUTHOR CONTRIBUTIONS

CJC performed the analyses. All authors contributed to the data, and to the conceptualization and writing of the manuscript.

Supporting information

ACKNOWLEDGMENTS

Notes

Carlson C. J., Bevins S. N., & Schmid B. V. (2022). Plague risk in the western United States over seven decades of environmental change. Global Change Biology, 28, 753–769. 10.1111/gcb.15966 [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]

REFERENCES