Spatiotemporal Variations of Plague Risk in the Tibetan Plateau from 1954–2016

Simple Summary Climate variability has influence on plague outbreaks worldwide. Usually, plague cases increase with increasing precipitation. Currently there are many studies on the epidemics of plague in human beings, whereas there are few studies on the dynamic of plague in animal. Nevertheless, animal plague is key in the natural epidemiological cycle of plague. We identified spatiotemporal changes of the plague territories in the Tibetan Plateau only using animal plague records. Our risky plague maps are far superior to the county-based maps used currently and have valuable applications for directly informing conservation and management decisions locally and regionally. Abstract Plague persists in the plague natural foci today. Although previous studies have found climate drives plague dynamics, quantitative analysis on animal plague risk under climate change remains understudied. Here, we analyzed plague dynamics in the Tibetan Plateau (TP) which is a climate-sensitive area and one of the most severe animal plague areas in China to disentangle variations in marmot plague enzootic foci, diffusion patterns, and their possible links with climate and anthropogenic factors. Specifically, we developed a time-sharing ecological niche modelling framework to identify finer potential plague territories and their temporal epidemic trends. Models were conducted by assembling animal records and multi-source ecophysiological variables with actual ecological effects (both climatic predictors and landscape factors) and driven by matching plague strains to periods corresponding to meteorological datasets. The models identified abundant animal plague territories over the TP and suggested the spatial patterns varied spatiotemporal dimension across the years, undergoing repeated spreading and contractions. Plague risk increased in the 1980s and 2000s, with the risk area increasing by 17.7 and 55.5 thousand km2, respectively. The 1990s and 2010s were decades of decreased risk, with reductions of 71.9 and 39.5 thousand km2, respectively. Further factor analysis showed that intrinsic conditions (i.e., elevation, soil, and geochemical landscape) provided fundamental niches. In contrast, climatic conditions, especially precipitation, led to niche differentiation and resulted in varied spatial patterns. Additionally, while increased human interference may temporarily reduce plague risks, there is a strong possibility of recurrence. This study reshaped the plague distribution at multiple time scales in the TP and revealed multifactorial synergistic effects on the spreading and contraction of plague foci, confirming that TP plague is increasingly sensitive to climate change. These findings may facilitate groups to take measures to combat the plague threats and prevent potential future human plague from occurring.


Introduction
In history, plague pandemics have had devastating effects on politics, economics, and demographics [1]. Plague is caused by the bacteria Yersinia pestis, which is directly changing environment. More specifically, We propose a time-sharing strategy which has the following advantages: (1) only plague data isolated from animals were used, which can minimize contingency; (2) environmental variables of ecological significance beyond climate surfaces were selected by considering Y. pestis, hosts, and vectors, which can improve the definition of ecological niche and habitat properties; (3) the plague data were matched by dividing them into five groups according to epidemic characteristics with the meteorological and vegetation data of the same period, which is helpful for matching temporal consistency. From all these data, five separate models from different times were built to define the geographical limits of plague foci as well as their changes. And their predominant influencing factors were identified in each model to clarify the feedbacks of spatial variation under changing climate.

Plague Data Compilation
Epidemic data for Y. pestis were compiled from published literature in the History of Plague Epidemic in China, Plague Prevention and Control in Tibet for 50 Years and Qinghai Plague. There were 1607 Y. pestis isolated in animals and 531 Y. pestis isolated from vectors from 1954-2016 in Qinghai Province and Tibet Autonomous Region. However, we can hardly direct exact geo-location for each Y. pestis from the constrained description, so we used infected regions data which were well documented as the spatial scale of village or hamlet. Then we transformed each occurrence site into spatial points using Google Earth and ArcGIS 10.8 software (ESRI, Redlands, CA, USA). By all the infected regions, we removed duplicates when they occurred in the same hamlet and year. Finally, 623 valid geographic points were compiled as epidemic data that applied in this paper.

Environmental Variables of Plague Foci
Twenty-two ecophysiologically relevant variables were initially selected as factors that play a significant role in plague occurrence. The explanations behind variable choice and generation of data layers are available in Tables 1 and S1. Generally, as far as the hosts are concerned, their density, distribution and migration are closely related to vegetation conditions. Vegetation growth depends on climatic factors which have indirect effects on hosts. In addition, seasonal variation of vegetation phenology may make hosts emerge earlier in the spring from hibernation. As for fleas, meteorological factors, especially temperature, play a key role on their life cycle, infection rate, bacterial virulence and so on. As for Y. pestis, although it remains controversial, it is certain the strains exist steadily in specific natural foci after a long period of natural selection and evolution. We hence selected soil and topography to delineate finer plague ranges.
Five types of data were collected: topography, vegetation, soil, climate, and river network data. Digital Elevation Model (DEM) data were obtained from the GMTED2010 dataset (Global Multi-resolution Terrain Elevation Data 2010 courtesy of the U.S. Geological Survey) in Google Earth Engine (GEE) [24]. The river network data calculated by DEM can be downloaded directly in RESDC [25], which was then used to calculate Euclidean distances in ArcGIS 10.8. Gravity data were provided by the Technical University of Denmark [26]. The vegetation data were also from GEE in the dataset of LANDSAT3/5/7/8 (GLS1975, LANDSAT SR from courtesy of the U.S. Geological Survey) [27]. Subsequently, normalized difference vegetation index (NDVI) was calculated and the median values of each pixel were exported into five stages: 1958-1979 (hereafter, S1), 1980-1989 (hereafter, S2), 1990-1999 (hereafter, S3), 2000-2009 (hereafter, S4), and 2010-2016 (hereafter, S5). The geochemical landscape data were derived from the Atlas of Plague and its Environments in the People's Republic of China [28]. Detailed classification of the geochemical landscape was in Table S2. The pH in H 2 O at a depth of 200 cm was obtained from GEE in the OpenLandMap dataset [29]. Soil data were obtained from Harmonized World Soil Database (HWSDv1.2), and FAO90 attributes were used. Their detailed list can be found on pages 24 and 25 of the PDF [30]. Soil moisture and climate data were obtained from TerraClimate in GEE [31] and were then processed as mean values for the above five stages. Meanwhile, the minimum monthly temperature values were also obtained between April and October ( Table 1). The reason we chose the above time scale was the consistency of the historical plague prevalence [8].
In addition, global terrestrial Human Footprint maps for 1993 and 2009 [32] were used to clarify the influence of human disturbance by grouping values based on their changes. Using the results from Section 3.2, we combined human footprint maps to determine the relationship between human activity and the spatial variation characteristics of plague areas. We calculated the difference in human footprint between 1993 and 2009 and divided it into six levels based on its changes. Moreover, the average risk and the percentage of areas with risk values greater than 0.5 for different levels of human disturbance were calculated.  [4] Hosts: a prolonged active season [41].
Most noticeably, highly correlated variables may result in over-fitting and hinder obtaining the actual response curves. Hence, we excluded redundant variables that had a high correlation (R ≥ |0.8|) determined by Pearson's correlation analysis [42,43] in R4.0.4 ( Figure S1). The most important variables selected are summarized in Table 1. All environment variables were projected onto Albers-CGCS2000 coordinate system. In particular, to avoid spatial mismatch between the size of organisms and the scale at which climate data were collected, we resampled all the variables to a spatial resolution of 5 km [44], which is the migration distance of marmot. Therefore, we regarded a 5-km spatial resolution as the basic area unit for plague infection among animals.

Modelling the Potential Areas of Animal Plague
Maxent quantifies the statistical relationship between predictor variables at locations where a species has been observed versus background locations in the study region. This software has been proven to perform the best in the potential species distribution assessment [45][46][47].
We divided the valid set of plague occurrence points into five subsets according to years (S1-S5). For each period, 10 models were established to guarantee that the results were stable, and every model was created with the occurrence localities by randomly selecting 75% of the occurrence localities as training data with cross-validation, reserving the remaining 25% for testing. The ultimate outcomes were selected from an average of 10 replicates. In addition, we assigned a combination of four types of features to generate models, including L (linear), Q (quadratic), P (product), and H (Hinge) features. We used default values in all runs, that is, threshold = 10 −5 , maximum iterations = 1000 for the algorithm convergence, and regularization value β = 1. The default regularization value was applicable for this study because we did not involve model transfer across either space or time. Each pixel in the study area was assigned a nonnegative probability value using Maxent. We chose the logistics outputs for easier usage and interpretation, which gave an estimate between 0 and 1 for risk of presence. The threshold segmentation was determined by the "maximum training sensitivity plus specific logic threshold" method: all pixels with a value greater than the threshold were classified as risk areas of plague [12,48]. Based on the classified results, the average risks over the TP and the total areas of risk regions were analyzed.
The model evaluation was first verified to perform significantly better than random [45]. Considering that we need to compare performance from different stages, the area under the ROC curve (AUC), which has been proved useful to compare performance between multiple Maxent models, was used to measure performance between the five averaged models. The closer the values are to 1.0, the better is the performance of the model. Models with values above 0.75 were generally considered potentially useful [46]. The importance of variables was evaluated using the Jackknife [49].

Model Performance and Changes of Plague Areas
Each average from the five stages of the 10 replicate models had high AUC scores (>0.9) and consistently performed significantly better than random across the entire spectrum ( Table 2). The results suggested our five-stage models had high performance, and the estimated distributions were a close approximation of the real probability distribution. The predicted risk area increased by 17.7 thousand km 2 from S1 to S2, decreased by 71.9 thousand km 2 from S2 to S3, increased by 55.5 thousand km 2 from S3 to S4, and finally decreased by 39.5 km 2 from S4 to S5 ( Table 2). The average risks showed similar fluctuation trends (Table 2). Consequently, S2 and S4 were regarded as stages of increased plague risk, but the risk decreased in S3 and S5.
Moreover, we compared our results with the officially announced areas of confirmed foci counties. The area of confirmed foci counties across the study period continually increased from 99.79 thousand km 2 in S1 to 687.04 thousand km 2 in S5, which were almost twice as large as our predicted areas in S3 and S4 and 2.6 times as large as that in S5 (Table 2). Figure 1 shows the identified plague territories of the five stages, and each of the five niches estimated almost exactly coincided with their respective occurrence points.

Plague Risk Areas at Different Time
Here, we were mainly interested in pixels that had a higher risk value than 0.5 to analyze spatial variation. The estimated ranges exhibited significant inter-decade variability from 1954-2016. Most noticeably, very high (risk > 0.7) and high (risk > 0.5) risk regions differed significantly between stages in both location and scope.  In the initial stage (S1), very high-risk regions tended to cluster around Qinghai Lake, the northern Qilian area, as well as Gonghe, Xinghai counties in Tibetan Autonomous Prefecture of Hainan and Nagqu area in northeastern Tibet ( Figure 1A). In S2, the high-risk regions spread and were heavily located in the southwest TP; in particular, there were many zones in and near the south of the Tibetan Autonomous Prefecture of Yushu in southern Qinghai ( Figure 1B). Simultaneously, some regions in Lhasa and the surroundings were at higher risk. In S3, plague risk throughout Qinghai Province was reduced ( Figure 1C). However, in S4, a handful of scattered high-risk regions appeared in the south of the Delingha, Wulan county, Dulan county in Haixi Mongolian, and Yushu city in Tibetan Autonomous Prefecture in Qinghai Province ( Figure 1D). In particular, Maxent prediction identified the middle part of the TP as an area of constant high risk for the geographical range of the plague. Specifically, the risk regions were larger in central and southern Tibet plague and were prolonged over a long period from S2 to S5 ( Figure 1B-E).
Overall, we identified abundant discontinuous niches for plague over the entire TP, and large differences in risk patterns were apparent from 1954-2016. The areas underwent repeated spreading and contraction of spatial ranges. Spreading occurred in S2 and S4 while contraction occurred in S3 and S5.

Impacts of Environmental Variables in Different Time
To determine the factors that led to the characterization of distributed differences, we analyzed the alteration of variable importance in each stage. The variable importance dynamics were analyzed using the jackknife test. The quantified results are depicted in Figure 2 where higher values of gain represent a greater contribution when training the models. The changing spatial patterns are profoundly associated with different factor combinations. Intrinsic features had definitive advantages in fitting the distribution of occurrence data in S1. Because these variables contained more useful information, the top three in order were GL (0.93), ST (0.78), and E (0.58). However, in the subsequent stages, ST had displaced GL, and it did not change until PR contributed the highest gain in S5 (Figure 2A). After that, there came to at least one variable feature in the first three rankings instead of merely intrinsic variables since S2. Notably, two variable features appear in S3. More specifically, in S3, the highest gains were ST, followed by PR (0.83), and NDVI (0.52). For S2, S4, and S5, ST, D, and PR had relatively uniform gains (even distribution around 0.5-0.6), but different sorting results.

Impacts of Environmental Variables in Different Time
To determine the factors that led to the characterization of distributed differences, we analyzed the alteration of variable importance in each stage. The variable importance dynamics were analyzed using the jackknife test. The quantified results are depicted in Figure 2 where higher values of gain represent a greater contribution when training the models. The changing spatial patterns are profoundly associated with different factor combinations. Intrinsic features had definitive advantages in fitting the distribution of occurrence data in S1. Because these variables contained more useful information, the top three in order were GL (0.93), ST (0.78), and E (0.58). However, in the subsequent stages, ST had displaced GL, and it did not change until PR contributed the highest gain in S5 ( Figure  2A). After that, there came to at least one variable feature in the first three rankings instead of merely intrinsic variables since S2. Notably, two variable features appear in S3. More specifically, in S3, the highest gains were ST, followed by PR (0.83), and NDVI (0.52). For S2, S4, and S5, ST, D, and PR had relatively uniform gains (even distribution around 0.5-0.6), but different sorting results. In general, the importance of the primary variables and other factors changed continuously across the study period, which resulted in diverse distributions of plague between stages. The main drivers shifted from solely predominant invariable attributes (elevation and geochemical landscape) to include both invariable and variable attributes In general, the importance of the primary variables and other factors changed continuously across the study period, which resulted in diverse distributions of plague between stages. The main drivers shifted from solely predominant invariable attributes (elevation and geochemical landscape) to include both invariable and variable attributes (soil and rainfall). In this conversion, variable importance analyses revealed an increasing explanation of precipitation for plague, with diffusion along a precipitation gradient ( Figure 2B). In addition, relatively uniform gains in the combined variables corresponded to wider spatial distributions ( Figure 2B). In contrast, the spatial ranges tended to shrink when each variable had different gains (Figure 2A). However, on the whole, variables including ST, D, GL, PR, and their combinations were always key factors for plague risk. They not only showed some of the strongest relationships (Figure 2A), but also contained unique information ( Figure 2B). Similarly, NDVI was an important variable that corresponded to the plague risk distribution.

Human Disturbance Evaluation on the Different Distributions
Our modelling environmental variables were natural factors; however, human factors, including urbanization and construction projects, could have altered the environment of some regions during the study period.
The mean values showed a decreasing trend under different degrees of human disturbance (Table 3), indicating that human activities had some restraining effects on plague intensity ( Figure S2). And the differences were highly significant from S2 to S3. Further statistics on the high-risk areas (>0.5) showed that these areas decreased steeply following human disturbance, but then increased after the change in human footprint was greater than 2 ( Table 3). The results revealed that the plague risk has been controlled by human activity to a certain extent since S3. However, with increasing human activity, the plague risk has gradually returned to pre-disturbance conditions.

Discussion
Plague regions of TP were identified using models, which represented the finer and true risky areas better. Because in the current publications, plague risk zonation and management were based on counties. And the risk areas were manually estimated based on the sizes of townships. In addition, the estimates are annual accumulations failing to consider plague dynamics, especially dynamics for decreasing risk. Besides, we found that the most evident indicator of Himalaya plague foci was the existence of independent foci with scattered distributions in the vast geographic regions over the TP. Indeed, the range edges, extents, and clustering of foci exhibited spatiotemporal variations caused by different climate variable combinations and human activity.

Basic Factors for Plague Distributions
The variables of intrinsic features are important delineators of the risk distribution of plague. The risk ranges tended to be in a basic niche when intrinsic features were the main factors impacting distribution. This contributed much to the qualitative analysis of foci, because the delineated ranges were closely correlated with the natural focus characteristics of the plague. For example, geochemical landscape and soil type, which had the highest gains in S1 and S3 (Figure S3), showed a wide border with risk regions. Our analysis indicated that high-risk regions were predicted in soils rich in calcium, and this result is compatible with the conclusion of previous studies that plague occurs more easily in calcium soils. In addition, we found that plague tended to occur in leptosols and gleysols, particularly in gelic gleysols [16], where excessive soil moisture accumulates during seasonal melts. The regions with lighter texture and shallower soil layers were ideal for marmots to dig deep burrows and were beneficial for keeping the burrows dry owing to their strong infiltration capacity [50]. Therefore, these soil types are essential and Y. pestis can persist for years in these soil conditions [37]. In that sense, the plague natural foci were concentrated in these areas. Similarly, the distance to a river was a relatively strong and useful predictor across all stages and the risk of plague was higher closer to rivers ( Figure 3A). Habitat preference may explain the close relationship between risk and D as infected hosts often suffer from dehydration as a symptom and need access to water to relieve discomfort; therefore, they are easier to find closer to rivers. 3A). Habitat preference may explain the close relationship between risk and D as infected hosts often suffer from dehydration as a symptom and need access to water to relieve discomfort; therefore, they are easier to find closer to rivers. In addition, the contribution of elevation continuously decreased, reducing from 0.58 to 0.1, indicating that E had less impact on risk distribution ( Figure 3E). The response curves of E showed that the highest risk occurred in S1 at an altitude of 3241.43 m, however, this subsequently shifted to higher but wider altitudes. This finding agrees with studies that found that marmots have widened their ranges and no longer only live on high frigid, subalpine, and alpine meadows at high altitudes [51], and are often found at lower altitudes, around farmland and vegetable fields.

Major Factors That Affect Plague Spatiotemporal Distributions
There was always a strong connection between plague variation and the difference in multifactorial synergy between variables. Climate factors, especially precipitation, are critical for spatial variation. Plague risk increased significantly when PR was between In addition, the contribution of elevation continuously decreased, reducing from 0.58 to 0.1, indicating that E had less impact on risk distribution ( Figure 3E). The response curves of E showed that the highest risk occurred in S1 at an altitude of 3241.43 m, however, this subsequently shifted to higher but wider altitudes. This finding agrees with studies that found that marmots have widened their ranges and no longer only live on high frigid, subalpine, and alpine meadows at high altitudes [51], and are often found at lower altitudes, around farmland and vegetable fields.

Major Factors That Affect Plague Spatiotemporal Distributions
There was always a strong connection between plague variation and the difference in multifactorial synergy between variables. Climate factors, especially precipitation, are critical for spatial variation. Plague risk increased significantly when PR was between 411.54mm and 493.64mm ( Figure 3B). During this study, the cycle changes of plague were in line with precipitation periodical oscillations (10 to 15 years) in the TP over the past several decades [8,52]. In addition, the spatial variation of the plague was consistent with the distribution of precipitation in each stage. For example, regions in Lhasa and the surrounding areas appeared to have more high-risk zones than Qinghai since S2. This was because there was more precipitation and rainy days in Tibet between 1980 and 2013 than in Qinghai [52]. Furthermore, NDVI was also a direct indicator of the diffusion of plague regions. Over the TP, there were many meadow types between northeast and southwest connecting plague niches, which can be easily influenced by precipitation and further cause wider risk distributions. Consequently, in S2, the extensive plague territories were caused by increased participation of NDVI. In fact, these fluctuations in plague and precipitation can be interpreted using the trophic cascade hypothesis. That is, precipitation increases food availability through vegetation-mediated effects, which increases rodent populations and activities, which causes an increase in risk zones [53,54]. However, while these were factors, plague distributions relied on combined effects of a complex synergy of factors rather than just two variables (rainfall and NDVI). For example, in S3, PR rose to 0.83 ( Figure 3A) and NDVI to 0.52 to become the third most important factor. However, large-area foci of spatial connectivity did not appear as a result of significantly varying contributions between the variables.
Simultaneously, other inconstant factors such as solar radiation, temperature, soil moisture, and PDSI also impacted plague risk, generating different risk zones. The SR variable had a high value of 244.41 for very high-risk areas (>0.7). The risk threshold of SR was more similar in S4 than in S1, both of which were stages in which the plague epidemic increased. The peak value T for risk shifted rightward along the axes since S1, exhibiting increasing risks with higher temperatures in the subsequent four stages. Indeed, warmer temperatures impact plague through both flea and bacteria (Y. pestis) survival [55]; higher temperatures advance the date of returning vegetation following winter, leading to the earlier revival of hosts [56,57]. Recent studies show that as temperatures increase, hosts adapted to cooler climates may suffer an increased risk of infectious disease outbreaks [58]. Furthermore, high irradiance values result in surface drying, combined with the finding that summer soil drying was exacerbated by earlier spring greening [59], which may increase migration for local food shortage, leading to more contact with infectious rodents. Moreover, marmots can suffer declines in body condition because of food decline, tending to exhibit easier defenses by infectious fleas [60]. SM and PDSI were equally variables related to water; thus, they were explained through the influencing mechanism the same as precipitation.

Human Disturbance for Plague Distributions
The reduction of plague because of increased human disturbance can be explained by two factors: first, the 1990s saw serious advancements in the prevention and control of plague, enabling comprehensive monitoring and control of plague among animals, which played an active role in weakening plague intensity ( Figure S2) [61]. Second, human disturbance destroys suitable habitats for hosts, resulting in lower abundance thresholds of host populations, thereby decreasing local plague transmission [62]. Unfortunately, other sites with modified habitats (e.g., anthropogenically restored ecological and agricultural sites) may later mediate reservoir abundance and cause changes in rodent and flea community composition, increasing plague risk [35,63]. Additionally, rodents enter sites where the ecology is more susceptible to plague [64]. This explains why plague increased in S5 in our results [65]. Therefore, environmental change should be stressed, and active surveillance programs should be coordinated to prevent the re-emergence of plague [66].
This study had several limitations. There may exist biases in our data as epidemic data were more often collected from general surveillance after human plague was identified. Some sites with adapted environments for plague may thus have not been explored; consequently, our samples may fail to represent all the possible environmental conditions, resulting in a risk underestimation. Additionally, the modelling results were geographical distributions that conformed to the constraint of all the environmental factors based on the epidemic data. However, the barrier restrictions, dispersal disequilibrium, or negative interactions were not considered [67,68], which could result in overestimation of risk [69]. Moreover, we used plague data from the two main hotspots, Qinghai Province and the Tibet Autonomous Region; however, there are risk areas in other provinces, such as the Xinjiang Autonomous Region, which were not considered or modeled (Table 2, Figure 1).
However, these limitations do not invalidate our results and we suggest that our results still demonstrate alterations of the risks and relative influences, which is essential for revealing the long-term effects of climate variation on plague. In the future, projections of plague risk variation under future climate change scenarios are necessary. Homoplastically, when considering human disturbance, it would be invaluable if deeper research could be devoted to understanding the long-term effects on plague risks of eco-relevant exploitation and conservation.

Conclusions
The study demonstrated differences in the ranges of plague for the five different stages over the TP from 1954-2016. For all stages, the modeled distributions were in agreement with the historical data. The modifications of plague risk regions revealed potential synergistic effects between plague and ecological variability. Specifically, we found that the degree of difference in primary controlling factor compositions changed the spatial pattern of the plague. Generally, climate change-induced wetting generates more favorable sites for plague and may further influence niche differentiation of the Himalayan marmot plague. However, despite the strong connection between climate and plague, we detected significant improvements in human plague dynamics. Our results had advantages on the quantitative description of plague dynamics characterization and generalization of occurrence regularity.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/biology11020304/s1, Figure S1: Correlation of environmental variables and occurrence in the five time; Figure S2: Example of risk time series obtained from Maxent results and changes of human footprint; Figure S3: Response curves of soil types and geochemical landscape; Table S1: Environmental variables list; Table S2: The classification of geochemical landscape.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.