Climate drives the spatiotemporal dynamics of scrub typhus in China

Abstract Scrub typhus is a climate‐sensitive and life‐threatening vector‐borne disease that poses a growing public health threat. Although the climate‐epidemic associations of many vector‐borne diseases have been studied for decades, the impacts of climate on scrub typhus remain poorly understood, especially in the context of global warming. Here we incorporate Chinese national surveillance data on scrub typhus from 2010 to 2019 into a climate‐driven generalized additive mixed model to explain the spatiotemporal dynamics of this disease and predict how it may be affected by climate change under various representative concentration pathways (RCPs) for three future time periods (the 2030s, 2050s, and 2080s). Our results demonstrate that temperature, precipitation, and relative humidity play key roles in driving the seasonal epidemic of scrub typhus in mainland China with a 2‐month lag. Our findings show that the change of projected spatiotemporal dynamics of scrub typhus will be heterogeneous and will depend on specific combinations of regional climate conditions in future climate scenarios. Our results contribute to a better understanding of spatiotemporal dynamics of scrub typhus, which can help public health authorities refine their prevention and control measures to reduce the risks resulting from climate change.


| INTRODUC TI ON
The potential impacts of climate change on public health have become one of the most challenging environmental issues (Patz et al., 2005;Rocklöv & Dubrow, 2020;Rogers & Packer, 1993). In particular, the effects of climate change on spatiotemporal dynamic patterns of vector-borne infectious diseases have garnered substantial research interest. Given that the viral and bacterial infectious agents and their associated vector organisms (e.g., mosquitoes, ticks, and mites) are largely devoid of thermostatic mechanisms, many aspects of their biology (i.e., reproduction and survival rates) are highly susceptible to changing climate (Caldwell et al., 2021;Franklinos et al., 2019;Iwamura et al., 2020). This fact, combined with the constant threat that vector-borne diseases pose to public health and the economy in many regions (Kilpatrick & Randolph, 2012), renders it crucial to quantify the impacts of climate change on such diseases.
Scrub typhus is a life-threatening vector-borne disease caused by the bacterium Orientia tsutsugamushi, which is occasionally transmitted to humans during feeding of larval mites (Lim et al., 2020). The disease was traditionally thought to be restricted to a well-defined geographic region called the "Asia-Pacific tsutsugamushi triangle," extending from far eastern Russia to northern Japan, northern Australia, the islands of the southwestern Pacific, and Afghanistan in the west (Kelly et al., 2009). It is estimated that scrub typhus potentially threatens more than 1 billion persons and causes 1 million clinical cases per year globally (Weitzel et al., 2016). The broad range of its clinical manifestations includes eschar, headache, fever, chills, rashes, enlarged lymph nodes, and mental changes (Cho et al., 2007;Paris & Day, 2020). Human infection may cause severe complications, including meningitis, meningoencephalitis, acute respiratory distress syndrome, myocarditis, multiple-organ failure, and bleeding diatheses, all of which can be fatal if not appropriately treated (Dittrich et al., 2015;Olsen et al., 2015). In recent decades, under the situation where no effective vaccines are available against the disease, scrub typhus has been jeopardizing public health with the growing number of cases of human infection and the geographical range of infected areas (Walker, 2016;Zheng et al., 2018).
Previous studies have documented that the spatiotemporal dynamics of scrub typhus are sensitive to climate factors Mayxay et al., 2013). For instance, transmission in southern China and northern Japan is evidently seasonal, occurring primarily in summer (Liu et al., 2006;Seto et al., 2017). On the other hand, in northern China and Korea, transmission is more prone to occur during autumn and winter (Kim & Jang, 2010;Liu et al., 2006). In addition, several attempts have been made to explore the climate-epidemic association by using statistical models (Kim & Kim, 2018;Seto et al., 2017;Wei et al., 2017;Zheng et al., 2018). For example, Seto et al. (2017) adopted a negative binomial regression to examine the effects of climate factors on scrub typhus in the Yamagata Prefecture of Japan. Kim and Kim (2018) used a hierarchical Bayesian Poisson model to reveal a positive correlation between scrub typhus cases and precipitation in summer in South Korea. Wei et al. (2017) employed two quantitative models to identify the key roles that climate factors (i.e., precipitation and relative humidity) played in driving the seasonal activity of scrub typhus in Guangzhou. Zheng et al. (2018) analyzed the spatial heterogeneity of scrub typhus in southern China from 2007 to 2017 through a boosted regression tree modeling procedure with several spatial correlates and concluded that the transmission of the disease is highly dependent on temperature and relative humidity.
Although the scientific community has studied the climateepidemic association for decades at different temporal and spatial scales (Kim & Kim, 2018;Kim & Jang, 2010;Seto et al., 2017;Wei et al., 2017;Zheng et al., 2018) (confirmed, clinically diagnosed, suspected). Individual-level data on human cases had been de-identified to protect patient confidentiality. Patients with the suspected classification were excluded from this study due to the uncertainty. All clinically diagnosed and laboratory-confirmed human cases were included in our analysis, with a total number of 166,839 from 2010 to 2019. It was determined by the National Health Commission of China that the collection of data on scrub typhus cases was conducted as part of ongoing public health surveillance of infectious disease and was, therefore, exempt from institutional review board assessment.

| Terrain data
Previous studies have suggested the apparent influence of elevation on the distribution of small mammals and their ectoparasitic chigger mites. For example, there are many more species of small mammals and mites found in Mountainous land than in cultivated flatland in Yunnan, China (Peng et al., 2018). Additionally, Zheng et al. (2018) reveal that elevation is one of the associate factors for scrub typhus in southern China. Thus, we adopted elevation as a covariate into our model. The gridded terrain data were downloaded from the Data Center for Resources and Environmental Sciences, Chinese Academy of Science (RESDC) (http://www.resdc.cn), with a spatial resolution of 1 km × 1 km.

| Historical climate data
Several lines of evidence suggest that the abundance and activity of chigger mites are strongly associated with climate conditions (Gubler et al., 2001;Jenkins, 1948;Sasa, 1961), which determine the chance of acquiring scrub typhus. For example, a positive association was observed between the rate of larval movement and temperature (Jenkins, 1948). Another example is the moderate amount of precipitation favouring humidity in habitats and food availability of chiggers, as opposed to the possible damage done to chigger mite eggs or scrub typhus nests by heavy precipitation to some extent (Gubler et al., 2001). In addition, various studies have explored the relationship between climate factors and scrub typhus, illustrating that the spatial distribution of the disease is affected by climate (Zheng et al., 2018). For instance, a study conducted in Guangzhou of China stated that every one-degree and one-millimetre increase in temperature and precipitation corresponded to an increase of 15% and 0.1%, respectively, in the number of monthly scrub typhus cases (Tiegang Li et al., 2014). To explain the spatial-temporal patterns of scrub typhus in China, mean temperature (°C), relative humidity (%), and precipitation (mm) were selected as covariates for Afterward, the projected daily climate data were aggregated to the month-county level, regarded as the future climate variables to predict climate-related scrub typhus situation in the 2030s, 2050s, and 2080s under various RCPs.

| Model specification
The association between local climate conditions and scrub typhus cases was estimated by a general additive model with a quasi- The term denotes the intercept; the dependent variable Y i,t was the number of scrub typhus cases at county i in time t (monthly/2010-2019); l was lag months; Elevation i represents the average altitude of county i; ε was defined as the error term; X i,t−l indicate climatic predictors (precipitation, mean temperature, and relative humidity) for month t − l; f() was smooth function defined by thinplate splines in this study. Considering the amount of time required for the mite lifecycle and the incubation period of this disease, lags of 0-3 months (l = {0,1,2,3}) were used to examine the delayed effects of climate factors on scrub typhus cases. We fitted the association between climate variables and scrub typhus cases under 0-3 month time lags. More specifically, we used the climate variables in the preceding 0-3 months to predict the number of cases in the current month. A term of fixed effect (Area i ) was included to account for the effects of unknown or unavailable variables in the model. We tested with both province and city fixed effects in our study to explain differences in the baseline level among areas. In summary, a total of eight model formulae were fitted in our study ( Table 1). The mgcv package of R was implemented for general additive model analyses. We used values of GCV, R, and deviance explanation as model evaluation criteria. A model with lower GCV, higher R, and deviance explanation was determined to be more predictive and hence selected.
The following steps of data manipulation and analyses were employed to fit the model. (a) Monthly aggregated scrub typhus cases were calculated from 2010 to 2019 at county level; (b) terrain and climate data were compiled from 2010 to 2019 as spatial covariates, which were extracted for each county using ArcGIS 10.6 software (ESRI Inc.); (c) the general additive model models described above were fitted based on assembled data to explore the relationship between climate covariates and scrub typhus cases, and the best-fitted models were selected judging by their performance; (d) projected climate data were used under three RCPs from four GCM models to predict scrub typhus cases in mainland China in the 2030s, 2050s, and 2080s based on the best-fitted model; (e) an ensemble approach was taken to average outputs from multi-GCMs under future climate conditions.

| Seasonality of scrub typhus
To quantify the characteristics of scrub typhus seasonality, we calculated the duration time and peak month of scrub typhus case numbers. The duration time was defined as follows: the cumulative case numbers captured in the interval between the start and end accounted for >95% of the total cases throughout the epidemic. The peak month was defined as the month with the highest number of scrub typhus cases.

| Historical spatiotemporal pattern of scrub typhus
Our analysis comprised all clinically diagnosed and laboratory-  Table S1.

| Climate driving factors on the spatiotemporal dynamics of scrub typhus
In this study, we build eight climate-driven generalized additive mixed models to simulate the spatiotemporal dynamics of TA B L E 1 Formulae and performance of climate-driven generalized additive mixed models. Based on the fixed effects of province (models 1-4) and city (models 5-8), the performances of models with 0-month (models 1 and 5), 1-month (models 2 and 6), 2-month (models 3 and 7), and 3-month (models 4 and 8) lags of climate conditions were quantified by values of the coefficient of determination (R), generalized cross validation (GCV) criterion and the proportion of deviation explained scrub typhus (see Section 2). The formulae and performance of these models are shown in Table 1. According to the results from Table 1, Model 3 and Model 7 had the best performances at their respective scales, with relatively higher deviance explained values and lower generalized cross validation (GCV) values, indicating that local climate conditions with a 2-month lag provided the best fits compared with 0, 1, or 3-month lags. The performance of models with city fixed effects (models 5-8) was generally better than those with provincial fixed effects (models 1-4) under the same climate conditions. Overall, Model 7 was the optimal choice for simulating the spatiotemporal dynamics of scrub typhus, with values for R, GCV, and deviance of 0.802%, 1.149%, and 64.3%, respectively.
Significant associations were found between climate factors and the spatiotemporal dynamics of scrub typhus (Table S2). The climate-disease nonlinear links are depicted in Figure 2 using the optimal model (model 7). There was an inverted U-shaped relationship between precipitation and the number of scrub typhus cases (F = 90.55, p < .001) with a peak at nearly 100 mm/month and a negative correlation for precipitation above 100 mm/month.
The relationship between mean temperature and the number of scrub typhus cases was complex (i = 2174.21, p < .001). For example, there was a significant negative correlation between mean temperature and the number of scrub typhus cases when the mean temperature was between −8 to 0°C, and a generally positive association when it was above 0°C. However, a slight decrease appears when mean temperatures reach 25°C. There was also a nonlinear but generally increasing association between relative humidity (F = 890.03, p < .001) and the 2-month lagged number of scrub typhus cases. Additionally, the potential nonlinear effects of ambient climate factors on the numbers of scrub typhus cases across lags of 0-3 months at two different scales are presented in Figure S3.

| Predicting scrub typhus dynamics under future climate changes
We predicted the spatial-temporal patterns of scrub typhus in the 2030s, 2050s, and 2080s under various RCPs using the best-fitting model (model 7) based on projected climate data derived from four GCMs (see Section 2). Figure 3 depicts the predicted total numbers of scrub typhus cases in mainland China under different RCPs, indicating that the projected number of scrub typhus cases is likely to peak in the 2080s under RCP4.5, increasing at 39% from baseline 2010s (Table S3). Under scenarios of RCP6.0 and RCP8.5, the predicted numbers of cases will peak in the 2050s, with an increase of 37% and 35% compared with the 2010s, respectively (Table S3).
Overall, we predicted minimal changes in the geographical patterns of predicted scrub typhus cases in mainland China, but signif- (3) Shandong, Anhui, and Jiangsu provinces, accounting for about 44%, 25%, and 15% of total projected cases, respectively (Table S4).
Changes in cases number differed among regions (Figure 4d-f 2080s under RCP8.5 ( Figure S5). In addition, we also quantified the model uncertainty in spatial predictions for scrub typhus based on the coefficient of variation values calculated for each county across the model ensemble ( Figure S6). The uncertainty analysis reveals low prediction uncertainty in most epidemic areas.
The estimated seasonality profiles for each province in the future under RCP6.0 present a generally extended epidemic duration compared with the 2010s. The epidemic season in northern China was predicted to extend in the next few decades compared to that in the 2010s. For example, the predicted seasonality in Shandong in the 2080s lasts from July to December, 1 month longer compared with its historical epidemic duration (from July to November). In addition, the scrub typhus cases tend to peak earlier in southern China.
For instance, the month when the case numbers peak in Guangxi and Guangdong changed from August in the 2010s to July in the 2080s.

| DISCUSS ION
Our study builds on previous research that has linked spatiotemporal dynamics of scrub typhus to climate variables (Kim & Jang, 2010;Seto et al., 2017) (Marks et al., 1996). More specifically, when the temperature drops too low or rises too high, it will be difficult for chiggers to feed with reduced opportunity for attaching to humans (Elliott et al., 2019;Moniuszko & Mąkol, 2016), while a temperate climate may trigger the development of eggs into larvae and thus increases exposure risk (Jameson Jr, 1967). Warmth could also enhance farming and recreational activities that contribute to high-risk exposure (Yao et al., 2019). Although the inverse U-shaped relationship between precipitation and the disease risk differs from previous studies (Chen et al., 2016;Yang et al., 2014), it may be explained by two opposing effects: precipitation increases the population of the larvae while reducing human outdoor activities. Humid environments are beneficial to the reproduction of most chigger mites, which may result in a positive association between relative humidity and the spread of scrub typhus (Lv et al., 2020). It is also important to note that the 2-month lagged effect may be due to the sum of the duration of the larvae development and the incubation period of the pathogen, which is nearly 2 months (Newton & Day, 2020;Shatrov & Kudryashova, 2006).
Previous studies have suggested that the future spatiotemporal distribution of vector-borne disease can be predicted using either a mechanistic model or a statistical model (Altizer et al., 2013;Gething et al., 2010;Messina et al., 2019;Morin & Comrie, 2013;Ryan et al., 2021). Compared with mechanistic modeling approaches that are highly dependent on the well-identified and temperature-based biological processes related to vector-borne disease spread (Messina et al., 2015), which ignore the effects of precipitation and relative humidity, statistical modeling approaches can better fit the empirical relationship observed from the more comprehensive historical data to predict the future risk. There are several studies about predicting the potential F I G U R E 3 The multi-GCM ensemble mean of the predicted number of national scrub typhus cases under different climate change scenarios (RCP4.5, RCP6.0, and RCP8.5) in the 2030s, 2050s, and 2080s. Rhombic points show the mean estimates of national scrub typhus cases and error bars are defined as the range.
risk of scrub typhus across different geographic locations Yao et al., 2019;Zheng et al., 2018). However, such predictions derived from yearly mean average values of covariates that only represented a long-term average distribution of the disease risk and could not infer seasonal patterns of scrub typhus. collection, data analysis, data interpretation, in preparing the paper or in our decision to publish. This work is funded by the National Natural Science Foundation of China (Grant No. 42201497). For the purpose of Open Access, the author has applied a CC BY public copyright license to any Author Accepted Manuscript version arising from this submission.

CO N FLI C T O F I NTE R E S T
The authors declare no competing interests.

DATA AV I A L A B I L I T Y S TAT E M E N T
The derived data that support the results and figures in this studyare openly available in Dryad (https://doi.org/10.5061/dryad.x69p8 czmz). F I G U R E 4 Spatiotemporal pattern of the predicted scrub typhus cases and its change over time under RCP6.0 scenario. The left panel depicts the spatial and temporal pattern of the four-GCM ensemble mean of the predicted scrub typhus cases for the 2030s (a), 2050s (b), and 2080s (c). Radar-Bar map denotes the seasonality characteristics of scrub typhus cases by province. The circumference is divided into 15 provinces in a clockwise direction, and the radius from inside to outside represents a particular month from January to December. Bar area represents the epidemic duration of scrub typhus, the black dot means the peak month, which defines as the month with the highest number of scrub typhus cases. The right panel shows changes in projected scrub typhus cases from the 2010s to 2030s (d), from the 2030s to 2050s (e), from the 2050s to 2080s (f), respectively. Colors determine different change types in cases number. Blue represents a decline, meaning that the number of cases decreased by more than 20. Red represents ascend, meaning that the number of cases increased by more than 20. Yellow represents stability, meaning that the number of cases fluctuated by no more than 20. Sunburst chart shows the proportion of provinces for each changes type. Map lines delineate study areas and do not necessarily depict accepted national boundaries.