## Abstract

In central Asia, the great gerbil (*Rhombomys opimus*) is the main host for the bacterium *Yersinia pestis*, the cause of bubonic plague. In order to prevent plague outbreaks, monitoring of the great gerbil has been carried out in Kazakhstan since the late 1940s. We use the resulting data to demonstrate that climate forcing synchronizes the dynamics of gerbils over large geographical areas. As it is known that gerbil densities need to exceed a threshold level for plague to persist, synchrony in gerbil abundance across large geographical areas is likely to be a condition for plague outbreaks at similar large scales. Here, we substantiate this proposition through autoregressive modelling involving the normalized differentiated vegetation index as a forcing covariate. Based upon predicted climate changes, our study suggests that during the next century, plague epizootics may become more frequent in central Asia.

## 1. Introduction

The great gerbil (*Rhombomys opimus*) is a major host (Pollitzer 1966; Gage & Kosoy 2004) of the plague bacterium (*Yersinia pestis*) in the central Asian deserts. Great gerbils are facultatively social rodents (Randall *et al*. 2005) where one or more family groups inhabit one burrow system. The number of burrow systems and their size remain fairly stable, especially the large, complex burrow systems in high-quality habitat (Naumov & Lobachev 1975). As population densities increase, empty burrow systems of lower quality become reoccupied. While adults may visit colonies up to 400 m away, dispersal movements seem limited to approximately 5 km (Rothschild 1978; Randall & Rogovin 2002). The intensity of reproduction is highest in early summer, with females usually having one but occasionally up to three litters per year (Naumov & Lobachev 1975). The gerbils spend most of the winter in the burrows in a largely inactive state, living off plant material stored during the summer. The climate of the central Asian desert is arid continental, with summer temperatures frequently above 35°C, winter temperatures below −20°C, large diurnal temperature variation and a mean annual precipitation of less than 200 mm. The study region consists mostly of scrubland, with a few more productive regions close to the rivers Ili and Karatal (figure 1*b*).

A number of arthropod parasites inhabit the burrows, and fleas (mostly of the genus *Xenopsylla*) act as the primary vectors of *Y. pestis* infection (Pollitzer 1966; Gage & Kosoy 2004). It has been shown (Davis *et al*. 2004) that gerbil abundance must be above a threshold level for sylvatic plague to spread in the population.

Thus, the spatial structure of gerbil population synchrony over large areas is of special interest since it is of fundamental importance for the frequency, scale and magnitude of the plague epizootics. Our ability to predict plague epizootics, and thereby the risk of transmission to humans, will be improved by increasing our understanding of the dynamics of the great gerbils at large spatial scales. To this end, we use monitoring data on gerbil dynamics covering the period from the late 1940s until the present, together with remotely sensed environmental data, to quantify the role of climate in determining and synchronizing the dynamics of great gerbil populations.

## 2. Material and methods

### (a) Population and prevalence data

Our study area is the PreBalkhash plague focus of southeastern Kazakhstan (74°–78° E and 44°–47° N). For monitoring purposes, the focus has been divided into 10×10 km^{2}, referred to as sectors. Four sectors constitute a 20×20 km *primary square* (PSQ), and four PSQs constitute a *large square* (figure 1*b*). Here, we use data aggregated at the PSQ scale as a trade-off between obtaining spatial resolution without losing temporal continuity. Figure 2 depicts the observed great gerbil density fluctuations.

Great gerbil density estimates (*N*) were calculated as the product of the number of burrows per hectare in the PSQ (*A*), the proportion of inhabited burrows (*O*) and direct count observations of the number of gerbils per burrow from several (10 per sector counted) of the inhabited burrows (*C*; Frigessi *et al*. 2005). Each spring (May–June) and autumn (September–October) during the period 1949–1995, *O* and *C* were recorded for 1–78 (median 54) PSQs.

There are also independent data on plague prevalence, where a mean of 201 gerbils (min 1, max 4734) were trapped during spring and autumn. Those gerbils caught were tested for plague through isolation of *Y. pestis* from blood, spleen or liver smears. These data overlap for about 79% of the population estimates at the PSQ scale. Gerbils infected with plague normally develop a detectable bacteraemia only briefly and sometimes intermittently (Gage & Kosoy 2004), and hence the data will underestimate the prevalence of, or even fail to detect, plague.

### (b) The normalized difference vegetation index data

The normalized difference vegetation index (NDVI; Los *et al*. 2000; Hall *et al*. 2005) is based on the difference between near-infrared and visible light reflected from the ground, thereby giving an index of light absorbed by chlorophyll on the ground from approximately 0 (bare rock) to 0.8 (rainforest). For a discussion of its use in ecological studies, see Pettorelli *et al*. (2005).

As high-quality NDVI data were not available prior to 1983, the April NDVI was modelled from the observed post-1983 spatial mean values and between-year differences in regional air pressure (the latter being available pre-1983; see electronic supplementary material) to obtain proxy data for NDVI for 1948–1982. Such pressure indices provide the direction and strength of the prevailing winds, often reflecting weather conditions. The commonly used NAO and ENSO indices (Stenseth *et al*. 2003) are calculated from such pressure measures, and more local indices are also used for modelling local climatic conditions (Hanssen-Bauer 1999).

### (c) Spatial population structure

The spatial synchrony was measured by the pair-wise correlation between (logged) PSQ population densities, also calculated at time lags of 0–6 seasons (Bjørnstad & Falck 2001). The local synchrony (figure 1*b*) was measured by the Spearman correlation between a PSQ and its eight immediate neighbours (or fewer if data from a neighbouring PSQ were not available).

The complex spatiotemporal autocorrelation structure in the data causes the loss of a number of degrees of freedom for several of the statistical tests performed, and no formal method for determining the correct reduction factor is currently available for a number of the tests and models. Hence, reported *p*-values should be viewed as approximations, but estimated 95% confidence intervals are provided where appropriate. To be conservative, we have based our conclusions only on relationships being significant after halving the degrees of freedom (had there been no spatio-temporal autocorrelation). We based our model selection procedure primarily on the models' ability to predict the first 32 years of data after having been fitted to the last 13 years only. For more details on model selection, see electronic supplementary material.

The selected models were then used for assessing the effect of climate, as indicated by the NDVI, on the seasonal dynamics and predicted densities of the great gerbil (see electronic supplementary material).

### (d) Simulation of population synchrony

In order to investigate the effects of population synchrony on plague prevalence, we simulated plague dynamics at different scales of host synchrony. The host population dynamics were represented by the Moran–Ricker equation (see eqn 7 in the electronic supplementary material) for each cell *k* in an *s*×*s* grid-based simulation. In each simulation run, the grid was divided into a number of sub-areas of size *α*^{2} (areas along the edges being fractions thereof when *s*/*α* was a non-integer). The value of the population growth rate *r* is drawn only once per area per year, only adding a small, uniformly distributed noise term *τ* with mean 0, and is thus highly correlated within each sub-area a, since(2.1)

Thus, when *α*=1, all *r* values are independent and spatial synchrony is nil, and when *α*=*s*, only the noise term differs between cells and spatial synchrony approaches one, over the whole area. *κ* was selected to give values of *r* mainly within the domain corresponding to stable population cycles. As there is evidence for transmission (Davis *et al*. in press) as well as host density affecting plague dynamics at these scales, the prevalence *P* of plague in a cell was determined by population density versus the plague invasion threshold (Davis *et al*. 2004) together with the density of infected cells within transmission range(2.2)where *T* is the epidemiological threshold value for plague to appear spontaneously or spread into the population, and is the ‘infection pressure’ by transmission from infected cells. was estimated by a two-dimensional Gaussian kernel density estimation (Venables & Ripley 2002) on the coordinates of infected cells. Increasing the bandwidth *Φ* of this smoothing kernel thus corresponds to longer transmission ranges, as the presence of plague in one cell gives a probability of infection of cells further away. The simulation was then run for 500 time-steps over a matrix of 50×50 cells for each of 3000 different combinations of *Φ* (transmission range indicator) and *α* (size of highly correlated subareas). Parameters are given in table 5 of the electronic supplementary material.

## 3. Results

### (a) The spatio-temporal population structure

The PSQ gerbil densities are positively spatially correlated across the entire focus. The mean correlation between autumn log densities decreases almost linearly from *ρ*=0.74 [0.60, 0.87] between adjacent PSQs to *ρ*=0.66 [0.35, 0.96] at 250 km. This is large compared to the yearly dispersal distance of the gerbils (Rothschild 1978), but similar to the spatial autocorrelation in between-year variation in NDVI_{April}, which decreases from approximately *ρ*=0.83 [0.75, 0.91] between adjacent PSQs to approximately *ρ*=0.65 [0.50, 0.80] at more than 200 km (refer to figure 2 and electronic supplementary material). April NDVI improves the predictive power of our population model, and may thus be expected to contribute to synchronizing abundances across space (Moran 1953). There were no indications of lagged correlations being greater than non-lagged correlations at any distance (see electronic supplementary material). Nor were they significantly greater in any particular direction.

In the southern part of the focus, there seems to be a breakdown of local synchrony (figure 1*b*), roughly corresponding to the landscape ecological region of the Akdala plain. The gerbil densities of the two other regions (the Bakanas (west) and Saryesikotrau (northeast) plains) seem to be highly synchronized both internally and with each other. Also, mean gerbil density does not vary between the Akdala/Bakanas/Saryesikotrau plains (factorial regression, *r*=0.01, d.f.=81, *p*>0.35). Altogether our results suggest a synchronizing exogenous factor (Moran 1953; Bjørnstad *et al*. 1999).

### (b) Effects of vegetation cover

We ask next whether the mean or variance in gerbil density or population growth rates is related to mean vegetation cover (as measured by the NDVI) or variability in NDVI (s.d.). Exploratory analysis suggested that April was the most informative monthly composite NDVI, i.e. NDVI_{April,} for predicting gerbil dynamics (see electronic supplementary material). Deviations from the mean NDVI_{April} within a PSQ for a particular year correlate weakly with population growth rate over the following winter (*ρ*=0.18, [−0.21, 0.52]). This suggests that more-than-usual vegetation cover in a PSQ gives rise to relatively high population growth. However, the mean gerbil density in a PSQ seems clearly independent of its mean NDVI_{April} (*ρ*=0.02, [−0.36, 0.40]). The standard deviation of the population growth rate (logged) from one season to the next is uncorrelated with mean gerbil density (*ρ*=−0.07, [−0.44, 0.32]), but it is negatively correlated with the mean NDVI_{April} of the PSQ (autumn–spring: *ρ*=−0.56 [−0.75, −0.29], spring–autumn: *ρ*=−0.45 [−0.67, 0.15]) and, less strongly, with the standard deviation of the NDVI_{April} (autumn–spring: *ρ*=−0.42 [−0.65, −0.12], spring–autumn: *ρ*=−0.29 [−0.65, 0.03]). (Note, though, that the standard deviation of the NDVI_{April} is strongly dependent on the mean (*ρ*=0.87 [0.78, 0.92]).) Thus, gerbil density fluctuations are greater in the drier, less productive areas, despite their lesser absolute fluctuations in NDVI_{April}.

The whole-focus mean gerbil density shows between-year autocorrelation for lags of up to 3 years for both spring and autumn, but there is no between-year autocorrelation in the mean NDVI_{April} of the focus (see electronic supplementary material). This autoregressive structure in population densities seems to overshadow the short-term effect of mean NDVI_{April} on the mean gerbil density of the focus (correlation *ρ*<0.05), but this correlation increases consistently as periods of greater length are compared, up to a moving-average period of 15 years, where *ρ*≈0.60, *ρ*<0.05 (see also electronic supplementary material). Thus, over moderately large time-frames, gerbils are more abundant during periods when there is more vegetation.

### (c) The gerbil population model

By model selection (see above and electronic supplementary material), we arrived at a log-linear autoregressive seasonal population model (see eqn 4 of the electronic supplementary material) providing good predictions on the PSQ scale up to two seasons ahead. In addition to one- and two-seasonal lag population densities, it included the mean difference in density from adjacent squares, the mean NDVI_{April} of the PSQ, as well as the yearly deviation from this. The plague prevalence detected in the previous season (averaged over adjacent squares) was also selected. The model coefficients are given in table 1 of the electronic supplementary material.

Combining the one-step predicted values of the spring and autumn into one time-series for each PSQ, they correlate closely (mean *ρ*=0.85, [0.75, 0.93]) with observations before 1983. For the period on which the model is fitted (1983–1995), the correlation is 0.87. In total, the model explains about 73% of the total variation in the log-population density. Its interaction terms provide a nonlinear relationship between climate variables and predicted densities. The model captures the greater variability in the drier (low NDVI) areas, exhibiting a correlation of *ρ*=−0.40 [−0.60, −0.12] at the PSQ scale between mean NDVI_{April} and the standard deviation of predicted log-population density.

Observed versus predicted time-series are shown in figure 1*b* (see also electronic supplementary material).

### (d) The plague dynamics

The chance of plague being observed in a PSQ is significantly increased if it, or its immediate neighbours, were observed to have plague in the previous season (*Χ*^{2}-test *p*<0.001). In addition, when plague is present, prevalence levels depend on the registered prevalence levels in the preceding season (correlation between logit prevalence spring versus previous autumn: *ρ*=0.54 [0.31, 0.71], autumn versus spring same year: *ρ*=0.51 [0.27, 0.69]).

Plague prevalence alone was found to explain a significant but modest part of the logarithmic population growth rate (autumn–spring *R*^{2}=0.03, spring–autumn *R*^{2}=0.04, autumn–autumn *R*^{2}=0.06, all *p*<0.05) at the PSQ scale. Prevalence was thus kept in the population model (as also indicated by the VAIC criterion; see electronic supplementary material) despite the effect of observed prevalence contributing a modest 1% to explained variation in the full model.

We calculated ‘epizootic magnitude’ as the product of plague prevalence and gerbil density for each PSQ at each time-step. Model simulations (see also electronic supplementary material) then predict that the frequency and magnitude of epizootics (which is closely correlated to the rate of human infection; Davis *et al*. 2005) increase with the mean NDVI_{April}. This was mainly due to a dampening of the density fluctuations, resulting in a more continuous presence of plague in the system, not to any increase in mean gerbil densities. In fact, when increasing the NDVI_{April} by less than one standard deviation, on average the mean population density is predicted to decrease with increasing median NDVI.

The relationship between the gerbil density and epizootic magnitude is also evident at the scale of the entire focus. Even though gerbil densities are part of the calculation, the highest correlations between mean epizootic magnitude and mean gerbil density were found using the autumn density 2 years previously in the case of spring magnitude (*ρ*=0.74 [0.38, 0.90]), and using autumn density the preceding year in the case of autumn magnitude (*ρ*=0.65, [0.23, 0.87]). Multiple regressions show that the density at these time lags continued to be highly significant even when controlling for present (same year, same season) density.

As climate signals of different relevance to plague dynamics may be reflected in the NDVI_{April} of dry versus river-watered regions, the area was split into two parts: the dry interior and the well-watered region close to the Ili River. The mean plague prevalence of the entire focus in the autumn was found to correlate with the mean April NDVI in the well-watered region in the preceding year (*ρ*=0.71, *p*=0.02). No correlation was found with the mean NDVI of the rest of the focus.

The simulation of population synchrony shows that longer transmission ranges (larger *Φ*), not surprisingly, increase mean plague prevalence under these conditions. The effective transmission ranges of plague in this system is currently unknown, but our results suggest that almost regardless of range, a considerable spatial autocorrelation in host densities is needed for the larger epizootics to occur. The mean prevalence (i.e. averaged over each simulation series) peaks for intermediate values of spatial autocorrelation, since plague dies out and has to ‘reinvade’ over longer distances when spatial autocorrelation becomes extremely high, and no alternative reservoirs are present in this simulation model. The maximum magnitude of the epizootics on the other hand continues to increase with increasing spatial synchrony (figure 3). Mean simulated population size is affected neither by *Φ* (transmission range indicator) nor *α* (size of sub-area with highly correlated population dynamics).

## 4. Discussion

The field of disease early warning-systems is currently receiving considerable research activity, but it is considered to be in its infancy, with few published studies using climate factors (as opposed to merely detecting pathogen activity) in a predictive framework (WHO 2004). One problem often encountered when studying zoonotic diseases is that the only comprehensive data available are the human cases (health records, etc.), with only patchy or conjectural information available for the wild host populations. The unique duration and scale of the present data provide a rare opportunity to understand the dynamics of zoonotic diseases in natural systems and their interactions with climatic effects. Overall, we believe that our results are encouraging, suggesting that ecological knowledge and data can be translated into predictive models and systemic understanding of a useful nature.

The effect of *Y. pestis* on gerbil population dynamics detected in these data is likely to be underestimated due to the technical difficulties of diagnosis (see electronic supplementary material), and the sampling frequency, which is on a coarse time-scale compared to the duration of detectable infection in the gerbils. Analyses on finer scales have suggested a negative effect of plague on individual recapture rates (Begon *et al*. 2006) and, though weaker and more variable, a positive effect on colony extinction rates (Davis *et al*. 2007). Thus, despite being relatively asymptomatic to plague (Gage & Kosoy 2004), an effect of plague on gerbils has now been shown on scales from the individual to the whole landscape. However, the effect seems weaker, or at least harder to detect, as spatial scale increases. While our observations are consistent with *Y. pestis* infections contributing to gerbil population declines, the limited predictive value of prevalence, and the highly restricted presence of plague in 1979–1981, when there was a focus-wide decline of gerbils, indicates that plague is far from being the main driving force in the great gerbil density fluctuations. It may well be that infections only appear with detectable prevalence in weakened, already declining, populations.

In arid scrubland with a simple vegetation structure, the NDVI reflects vegetation cover closely (Los *et al*. 2000; Pettorelli *et al*. 2005). The absence of correlation between mean gerbil abundance and mean NDVI is consistent with the assumption (Naumov & Lobachev 1975) that mean great gerbil densities are most strongly related to soil conditions (suitable patches for burrowing). In this light, the fact that we found a higher variability in gerbil densities in less-productive habitats should not be unexpected, as the impact of between-year fluctuations in resources is likely to be greater in less-productive environments when mean population densities are the same.

The high degree of spatial synchrony in gerbil density is striking. While migration may play some role in this, its effect at distances up to 10 times the recorded maximum yearly migration seems unlikely without a time lag. At these scales, climatic fluctuations are much more likely to play a major synchronizing role, generating a ‘Moran effect’ (Moran 1953). This may be mediated through several mechanisms, such as weather conditions affecting food availability, temperature stress, parasite (flea) density and plague transmission (Cavanaugh & Marshall 1972; Stenseth *et al*. 2006). Some avian predators of the steppe (*Buteo rufinus*, *Falco tinnunculus*, *Aquila nipalensis*, *Aquila heliaca* and others; Sanchez-Zapata *et al*. 2003; Sludsky 1978) may range over distances comparable to our study scale, and their abundance would be expected to respond to climate forcing on their prey, linking the synchronizing effects of climate and predation.

As the climatic factors affecting gerbil dynamics (and possibly those of other rodent hosts) in part reflect large-scale fluctuations (Todd & Mackay 2003; Treydte *et al*. 2006), the observed processes seem likely to continue into similar ecosystems in the region. The presence of spatially continuous host populations whose population densities exhibit large-scale autocorrelation seems very likely to be a key factor in allowing large-scale plague epizootics in the region. The co-dependency between NDVI_{April}, gerbil densities and plague, together with predictions from historical climate records (Stenseth *et al*. 2006; K.L. Kausrud *et al*. 2007, unpublished data), suggests that periods of relatively warm and/or moist conditions give rise to periods of high gerbil densities and large epizootics in otherwise dry areas. Global change is predicted to continue towards a warmer climate in arid central Asia, while the heightened intensity of the hydrological cycle increases the variability of monsoon-related climates extending into this region (IPCC 2001). Epizootic magnitude can thus be expected to increase, as the highest prevalence levels seem to occur in the most arid areas and after warm springs. This is supported by the model simulations that suggest greater and more frequent epizootics when spring NDVI increase and population synchrony is high. While extremely warm summers have negative impacts on plague prevalence, this seems likely to be outweighed by the synchronizing effect of extreme climate events and the effect of large rainfalls (and increased glacier runoff) when these occur. On somewhat longer time-scales, the ranges of desert-adapted plague hosts in the area, whose mean densities are relatively insensitive to the degree of aridity, seem likely to expand. When coupled with sporadic bursts of wetter conditions, this is likely to cause plague epizootics to increase in scale.

To summarize, we draw five main conclusions from these analyses.

Great gerbil density fluctuations are highly correlated over large areas, suggesting a climatic forcing as a synchronizing agent. This is probably an important factor causing large-scale plague epizootics in the region.

Great gerbil population densities at large spatial scales can be well predicted 6–12 months ahead when combining spatial environmental effects and intrinsic dynamics. This is important for predicting plague dynamics.

While great gerbil population growth rates exhibit greater variability in areas with low NDVI

_{April}, average population density is not strongly related to average vegetation productivity. This suggests that the gerbils will be capable of maintaining population densities where plague can persist over most of their range even if, as predicted, the climate in central Asia gets increasingly arid.While the presence of plague infection in an area is associated with population decrease over the following months, plague seems unlikely to be the main driving force behind great gerbil density fluctuations.

The magnitude of plague epizootics associated with the great gerbil may be expected to increase under predicted effects of ongoing climate change.

## Acknowledgments

This work was supported by the Norwegian Research Council, the European Union project STEPICA (ICA 2-CT2000-10046) as well as the CEES and the Centre for Biostatistical Modelling in the Medical Sciences (BMMS). Discussions with Tamara Ben Ari, Kung-Sik Chan and Dorothee Ehrich, who also provided valuable help translating between Russian and English, have been greatly appreciated. This project would have been impossible without the cooperation of Vladimir S. Ageyev, Nikolay L. Klassovskiy and Sergey B. Pole as well as hundreds of other Kazakh scientists who collected this vast amount of data over one hundred years. We thank Jan Esper for advice on palaeoclimatology. Finally, we wish to thank the anonymous reviewers who provided valuable feedback on an early version of this paper.

## Footnotes

Electronic supplementary material is available at http://dx.doi.org/10.1098/rspb.2007.0568 or via http://www.journals.royalsoc.ac.uk.

- Received October 22, 2006.
- Accepted April 30, 2007.

- © 2007 The Royal Society