Identification of Chinese plague foci from long-term epidemiological data

Clustering.

Here, the entire dataset in its gridded form is used to identify large-scale patterns of plague occurrence dynamics. Our method (SI Methods contains descriptions) reveals spatial coherency of local TS (Fig. 1A). The distribution of these groups is mostly independent from the number of groups, N. Precisely, for n ≥ 3, groups emerge that gather a large number of grid cells: SW, the western part of SE and NW for group 1 TS, NC for group 2, and the eastern part of NE plus SE for group 3 TS. Groups 1 and 3 thus cover the four suspected nonepidemic plague territories revealed by Fig. S2, whereas group 2 TS are associated with the 1894–1952 epidemic period. The robustness of the aggregation is demonstrated by the comparison between grouping results for n = 3 and n = 6. Adding more groups only removes a very limited number of outliers from each group (not shown for clarity). Other parameters than N could also influence our results. To check their robustness, we carried out a series of sensitivity tests. A small number of TS (<10%) may belong to one group or the other depending on changes in the parameters (e.g., N), which does not invalidate the main results above or the interpretations made in the discussion.

TS for each of the three groups are presented in Fig. 1B and reflect the importance of time patterns in the grouping. Main characteristics of each group TS are as follows: reports distributed throughout the whole time period with a long period of above-normal (epidemic) conditions for group 1, a vast majority of plague reports between 1937 and 1952 for group 2 and finally, two to three bursts in 1917–1918, 1928, and 1932 with recurrent low incidence reports after 1900 (no cases before) for group 3. Noticeably, the two coherent and geographically distinct sets that form group 3 remain aggregated for up to n = 15. Also, rerunning the clustering for that subset does not indicate major dissimilarities within it (this is also true for groups 1 and 2, not shown for clarity). Connection between southeast and northeast territories is an important result of this method that we downplay in Discussion.

Ecological Niche Modeling (ENM).

ENM allows the characterization of environmental differences between plague territories. To avoid including accidental/nonsignificant occurrences and imported cases, we use plague occurrences having more than one plague occurrence over the nonepidemic periods (henceforth labeled “valid”) as model inputs. Here we present results obtained using the Maxent software but we also used GARP. Maxent is known for overfitting the data and GARP for overpredicting the extension of suitable areas (21, 22) (Fig. S3 provides a visual comparison of the outputs of two runs with identical inputs). However, the fact that both models predict essentially similar suitable territories (not shown for clarity) when run over regional input subsets (black dots in Fig. 2) gives us confidence in the robustness of the conclusions.

We first investigate environmental differences between north and south territories by dividing the valid set of plague occurrence points into northern and southern subsets (Fig. 2, Middle black dots). Two well-circumscribed niches are predicted by ENM that each roughly coincide with their respective occurrence points. Most noticeably, none of the areas predicted by any individual subset encompasses another subset’s realm. They are well separated from each other by the central China area that our ENM model successfully predicts as plague-free. The area predicted by the southern subset encompasses northern Viet Nam where plague is known to have occurred (23, 24) and parts of Nepal and northeastern India (which is not inconsistent with ref. 25). Other areas are predicted as suitable for plague but cannot be completely confirmed due to the scarcity of data (e.g., Mongolia or Russian borders). In Discussion we elaborate on why a territory could be suited for plague but have no known expression in humans. For several retained environmental variables, little or no overlap is observed between northern and southern niches. Plague is predicted in very dry to arid areas in the north. The southern niche is generally wetter (annual precipitation ranging from ∼750–1,700 mm.y⁻¹ in the south versus 0–700 mm.y⁻¹ in the north) and warmer (8–27 °C versus −5–8 °C). This niche has significantly fewer marked temperature seasonal fluctuations as revealed by the annual temperature range (20–25 °C with outliers around 30 °C compared with 35–50 °C in the north). Interestingly, the total predicted niche (north + south) covers roughly the entire available value range for most of the input chosen WorldClim (26) variables. One exception is that endemic plague is not predicted for the coldest temperatures (i.e., areas with an average annual temperature below −5 °C that also have the highest annual temperature range).

To confirm the regional coherence found with the cluster analysis (Clustering) and the visual pooling presented in Preamble on the Data, we now subdivide the set of valid plague occurrence points into the four patches that geographically stand out in Fig. S2 B and D. ENM-predicted plague niches for each regional subset are represented in Fig. 2. The SW niche predicts an area corresponding to the Yunnan province. The SE niche covers Fujian, Guangdong, Guangxi, and Guizhou, the latter being outside the epidemiological plague area. The NE niche almost exactly matches the region delimited by its associated plague occurrence points. The NW niche extends over a much broader territory that covers most of northern China. The predicted area spans the NE territory and successfully includes the NC area where plague is observed only during the epidemic period and was thus excluded from the ENM input subsets (but see finer scale investigation in Fig. S4).

Distribution ranges for the seven input environmental variables across subsets (Fig. S5) show that differences are comparatively smaller within the northern/southern niches than between them. This is not unexpected because most relevant ENM variables vary primarily with latitude. We also find well-defined and separated ranges for climate variables associated with each subset. In particular, the ranges for the NW niche are surprisingly restricted in regard to its geographical coverage. This niche corresponds to the coldest and driest plague territory (Fig. S6). Conversely, the warmest and wettest territory is the SE, which is noticeably different from NE with respect to average temperature and precipitation.

North/South Distinction.

The most evident characteristic of Chinese plague is the existence of at least two distinct plague territories; one in the north and one in the south. Support for this can first be found in the ENM analysis that predicts two well-separated niches corresponding to northern and southern subsets and distinct ranges for their environmental variables. This is also corroborated by the fact that our data contain no plague reports in the in-between area (from about 30–35° N with small variations depending on the longitude) although it is densely populated. Field observations also confirmed that no suitable territory exists in central China (20). This indicates that there is no evidence of a plague pathway joining northern and southern plague territories. In addition, ENM prediction results show that none of the four subsets has environmental conditions that project into central China (Fig. 2 A, C, D, and F). Biomes and WorldClim data are consistent with the existence of two environmentally distinct territories: an herbaceous cover with increasing elevation toward the west in the north versus evergreen tree or shrub cover in the south (Fig. 3). Therefore, a conclusion from our findings is that long-term coexistence of all three plague components (i.e., pathogens, rodents, and fleas) is possible under distinct environmental conditions.

Northeast/Southeast Puzzle.

As group 3 (green dots in Fig. 1A) encompasses both northeast and southeast territories, the clustering analysis may seem more ambiguous about a north/south disconnect. However, the procedure through which TS have been aggregated needs to be considered. In fact, visual inspection of the local TS (Fig. 1B) indicates that TS composing group 3 share the following characteristic: a pronounced reported incidence of plague in the short period 1937–1952 not seen in the other groups (fraction of gridded counts of infected villages during 1937–1952 over the entire period are 0.65 for this group versus 0.25 for group 1 and 0.23 for group 2). The grouping of northeast/southeast territories together in group 3 may thus occur for circumstantial reasons that should not be overinterpreted. These areas simultaneously underwent Japanese invasion during 1937–1945. In fact, according to Liu (20), plague was repeatedly spread by the Japanese army in the southeast during 1940–1948. The complete absence of overlap between SE and NE ENM niches also provides evidence against a plague focus spanning the entire group 3; they do not share similar environmental traits (Figs. S5 and S6). In fact, ENM analysis further suggests that the southern subset of group 3 may not even be part of a plague enzootic territory. Indeed, the corresponding infected gridded incidences only occurring during the war period sit outside the niche predicted by southern and SE occurrence points (green dots in Fig. 2 E, Inset and F). This is generally consistent with the reconstruction of plague spread proposed by Benedict (4) that involves contact diffusion (i.e., spread outward from a central focus) from Fujian to Zhejiang during the epidemic period of 1884–1949. This is also consistent with the idea that, during epidemics, plague spills over from favorable (i.e., where plague is enzootic or at least can maintain itself over relatively long periods of times without being regularly imported) to unfavorable environments such as the area where the southern green dots are located (Fig. 2E, Inset). The southern points that make up group 3 do not fulfill the conditions to qualify as an enzootic plague territory and do not refute the existence of a north/south plague disconnect.

Southern Plague Region.

Plague habitats in the south are quite similar (Fig. 3): tropic/subtropic forest zone with homogeneous annual temperature range and average and, although lower in Yunnan, comparable annual precipitation (Fig. S5 A, D, and E and also ref. 20). Despite this relative coherency, SW is characterized by higher elevations probably responsible for observed differences in mean diurnal range and precipitation seasonality (Fig. S5 B and G).

Our epidemiological data suggest that at least two plague regions exist in southern China with some degree of independence between them and most of the Guanxi province (plague-free during the entire study period) as a physical separation. The first region corresponds to Yunnan (SW) and the second one to Guangdong and Fujian (S and SE). Support for the existence of two separate regions can be found in the relative number of plague reports as a function of time: the first plague record outside Yunnan is in Beihai (Juangxi province) in 1867 (i.e., plague remained within the boundaries of Yunnan for at least a hundred years after the first confirmed record in its borders) and the first occurrence in Fujian is in 1884 at a time of mild epidemic conditions in Yunnan. On the other hand, no plague cases were reported in Yunnan from 1956 until 1964, whereas concomitant sporadic occurrences were reported in SE. ENM analysis would also be consistent with the existence of two subregions in the south. The SW predicted niche does not include eastern territories, whereas the SE subset of occurrence points predict a wider niche that spans into southeast Asia and some areas in Yunnan, albeit with a lower suitability index. Interestingly, field studies conclude that Rattus flavipectus is the unique dominant plague host for the whole south (Fig. 3) but with variation: Apodemus chevrieri as a secondary host in some areas in the south and SE and Eothenomys miletus in the highland of western Yunnan (20). It is difficult, at this resolution to strongly argue in favor of either a unique or two plague foci in the south.

Physiographic macroregions (28) are self-contained regions with major geomorphologic features (constraining travel and exchanges, thus coherent in terms of economic resources and development) and drainage basins. The totality of what we have called southern China is separated into four macroregions: Yungui, roughly corresponding to the Yunnan province; Lignan, to Guangxi and Guandong; southeast coast, roughly to Fujian and Zhejiang; and part of lower Yangzi, roughly to Jiangsu. According to Benedict’s work (4), plague first slowly spread within Yungui. During the second part of the 19th century, coincident with the Muslim rebellion (1855–1873), it spread outside of its source province and spread to Lignan, the southeast coast and lower Yangzi. The reconstruction work by Benedict shows how plague spread along trade routes and within the boundaries of the regional trading system, explaining the link between plague dynamics and Chinese macroregions. Cities and ports served as distribution centers where infected rodents or fleas were passively transported in merchandise (4). According to Benedict, a secondary spread then would have occurred from newly infected cities inland. This is consistent with a spread from the west to the east we have described and a subsequent settlement of plague within the contours of a favorable plague niche found in our ENM analysis (Fig. 2). This is also consistent with the similarity observed for the south plague system (similar dominant hosts and fleas, interannual plague dynamics, and environmental characteristics).

Northern Plague Territories.

In contrast to the south, northern plague territories exhibit important heterogeneities and a wider variety of biomes (Fig. 3). Fig. S2 reveals that there could be three distinct patches of plague occurrence in the north: one in the east (NE as captured by group 3 in the cluster analysis; Fig. 1A) and two other, northwest (NW) and north central (NC, essentially group 2) that are somewhat geographically isolated from each other. Only two (NE and NW) would qualify as endemic Fig. S2D) because NC has no reporting outside the epidemic period.

NE local TS mostly belong to cluster group 3 because they have a concentration of reports in the 1937–1952 period (although less so than for the southern part of this group, as we have verified). We suspect this concentration should be even less marked had the 1910–1911 and 1920–1921 epidemics been fully reported (10, 29, 30) (partial control of the Russians and Japanese over the region may be responsible for underreporting or loss of information). Historical reports suggest an endemic focus corresponding to NE with no or little expression in indigenous humans; pneumonic outbreaks later occurred in relation to an increased hunting of marmots (31). Infection apparently followed rail routes throughout Heilongjiang and northeastern Inner Mongolia. There are several clues suggesting that the NE plague territory is isolated from the rest of the northern plague region. Most importantly, the ENM analysis using the NE subset predicts a very restricted territory corresponding to a cultivated wet/subwet steppe biome (Fig. 3 B and D) that differs from NC and NW in this regard.

In contrast, the NW/NC plague region is more heterogeneous. Although plague reports are essentially contiguous, our analyses suggest it is made of at least two distinct plague territories. Village reports for NC (corresponding to the Inner Mongolia, Shanxi, Shaanxi, and small parts of Gansu provinces) mostly belong to group 2, whereas those in the far west mostly belong to group 1 (Qinghai province). NW group 1 TS have plague reports over the entire study period, including at the very beginning of the time series, an indication of the enzootic nature of this plague territory. On the other hand, group 2 reports occur during the epidemic period with two major peaks in 1917–1918 and 1931, indicative of epidemic conditions (>180 villages for each peak). However, low levels of plague occurrences are also found throughout the 1900–1950 period (first and last reported cases in the region); in 4–10 y per decade, 1–17 villages reported plague. This is compatible with the existence of an enzootic area (ancestral as it has been suggested elsewhere (10) or more recently established due to successful and sustainable importation of plague) whose precise location can be inferred from the positions of villages infected during these low plague level years (Fig. S4). A complementary ENM analysis with NC group 2 locations as input predicts a niche that is roughly consistent with this set of points. The Mongolian gerbil (Meriones unguiculatus) and the Daurian ground squirrel (Spermophilus dauricus) could be the dominant plague hosts in most of the NC plague territory (Fig. 3). However, because NC occurrences are reported during the epidemic period, the geographical range of NC points is undoubtedly wider than that of the NC actual endemic core. In terms of biomes, NC occurrences occupy a region straddling several environments (desert, arid to dry steppes; herbaceous cover to bare areas). On the other hand, NW points that belong to group 1 are associated with an alpine-type biome with marmots (Marmota himalayana and M. baibacina) as main hosts.

All these elements (behavior of the time series and hence grouping; ENM niches, hosts, and biomes) suggest the existence of three plague foci in the north with an important degree of independence between them. Interestingly, there are three physiographic macroregions in the north: NE corresponds well to the northeast macroregion; both NW and parts of NC are included in the northwestern macroregion; to the contrary, the north macroregion as defined in ref. 28 is plague-free in our dataset, except for a few single occurrences that we deem insufficient to support the existence of an endemic plague focus there.

Note however that according to different sources (9, 29–31), plague was endemically present in northern territories and went unnoticed before the third pandemic when it spread as a result of population increase and transit system development (29, 30). Similarly, it is probable that plague presence in the NW was underreported because this territory is scarcely populated. Overall, we conclude the possibility of at least three independent foci in the north.

Classification Method.

The 744 time series, corresponding to the aggregated yearly number of villages reporting plague in a 0.5 × 0.5° grid, are, as most often the case with epidemiological time series, nonstationary (i.e., their mean, variance, or frequency components change in time). We thus apply wavelet analysis where the original time series are decomposed into oscillatory components. In contrast to Fourier analysis, this time–frequency analysis captures the amplitude of frequency components as a function of time (32, 33). For further details see SI Methods.

Ecological Niche Modeling.

We characterized spatial patterns of the landscape that are necessary for plague maintenance or transmission. Our approach is based on ecological niches defined as the set of environmental conditions under which a species is able to maintain populations without immigration (34, 35) (or, in this case, the conditions that prevail in plague endemic territories). Known occurrences of species (here plague occurrence points) can be related to raster coverage summarizing environmental variation (36). We used available bioclimatic layers at resolution ∼10–14 km from the WorldClim (11) data archive (1950–2000). For further details see SI Methods.

To avoid fitting models in too many environmental dimensions, we chose the 7 least correlated (P < 0.6) of the 19 available variables (Fig. S5). These are average and range of annual temperature, average diurnal temperature range, isothermality, annual precipitation, precipitation of the driest month, and precipitation seasonality. We used two algorithms for ENM development based on presence-only data (no real absence points): the genetic algorithm for rule-set prediction (GARP) using the software DesktopGarp (38, 39) and via maximum entropy optimization approaches using the program Maxent (40–42). For further details see SI Methods.