DOI: https://doi.org/10.21203/rs.3.rs1218122/v1
Background: Global warming has compromised human health by increasing the incidence of infectious diseases. Scientific evidence is required to expand the knowledge of the association between meteorological factors and the incidence of infectious diseases. Our study focused on meteorological factors such as ambient temperature (AT), humidity, and scrub typhus incidence.
Objective: We aimed to investigate the longterm effects of AT and relative humidity (RH) elevation on scrub typhus incidence in South Korea.
Method: Meteorological variables were combined with scrub typhus cases reported from 2001 to 2019 in South Korea. A generalized additive model was used to explore the nonlinearity of the lagged association between meteorological variables and weekly scrub typhus incidence. To explore the longterm association between meteorological factors and scrub typhus incidence, the difference between annual mean AT or the annual number of heatwave days during 2001 to 2019 and those of the period 1971–2000 were linearly regressed on annual cumulative scrub typhus incidence.
Results: Association between weekly scrub typhus incidence and mean temperature or relative humidity with 15 weeks lags showed V or U shaped relationship. Above the threshold temperature (14.9°C to 17.0°C), scrub typhus incidence increased by 1.85% (95% CI: 1.5, 2.2) per 0.1°C elevation in mean temperature with 15 weeks lags. Scrub typhus incidence increased by 3.7% (95% CI: 2.7%, 4.8%) and 2.5% (95% CI: 1.6%, 3.3%) per 0.1°C increase in annual mean temperature and per one day increase in heat wave days.
Conclusion: Mean AT elevation and RH increase in summer were associated with an increased incidence of scrub typhus in the fall. Increases in annual mean AT and annual number of heatwave days were associated with an increase in the incidence rate during 2001–2019.
Global mean ambient temperature (AT) has increased by approximately 1°C as compared with the preindustrial era. The Intergovernmental Panel on Climate Change (IPCC) has warned of more frequent flooding, reductions in crop production, and disruptions of ecosystem functions ^{1}. Climate change alters environmental conditions and, consequently, increases the incidence of infectious diseases, including vectorborne diseases ^{2,3}. Increased rainfall may induce an increase in vector density, and high temperatures may lead to an extended incubation period for the vectors. AT and relative humidity (RH) ^{4} may be related to determining vector survival, such as affecting the distribution and life cycle of vectors ^{5–7}. Ticks carry Rickettsia as an endosymbiont and transmit these pathogens to humans through human–vector interactions.
Mites are the most common vectors of infectious diseases worldwide ^{8,9}. Among these, scrub typhus, denoted as “tsutsugamushi disease” in Japan, has become a critical concern in endemic regions as well as in other areas as traffic to and from the endemic regions increases in numbers. ^{10,11}.
AT and RH are known to be key factors in the development of Leptotrombidium, the vector mites of scrub typhus ^{12}. The response of chiggers (“mites”) to climate change is associated with scrub typhus incidence. In a previous study, climate factors such as AT, RH, mean wind speed, precipitation, and sunshine duration ahead of the outbreak season were associated with scrub typhus incidence during the peak periods ^{13,14}. The geographical distribution of scrub typhus vectors across the Republic of Korea (ROK) ^{15}, which has been studied since 1980, has been dominated by two species, Leptotrombidium pallidum in the northeastern provinces and L. scutellare in the southwestern provinces ^{16}. Japanese studies found that L. scutellare grew from mites to adults between April and August, whereas L. pallidum emerged from June and peaked in late July, August, and early September ^{17}.
Scrub typhus usually causes symptoms such as fever, malaise, skin rash (e.g., eschar), and conjunctival injection. However, severe illnesses such as secondary acute meningoencephalitis with neurological sequelae or even fatalities can occur in rare cases ^{18–21}. To date, the geographical spread of the disease and its vectors has expanded ^{10}, and climate change and global warming are suspected to be the cause of the spread of scrub typhus ^{2,3,22}.
Considering the geographical distribution of vectors in the country, we assumed that the incidence of scrub typhus in the ROK might be influenced by climatic conditions. Hence, we aimed to explore the association between climatologically preferable conditions and scrub typhus incidence from 2001 to 2019 in the ROK, with a focus on AT and RH.
1. Seasonal variation of climate factors and scrub typhus
Weekly scrub typhus incidence and meteorological variables showed seasonal variations from January 1, 2001 to December 31, 2019 (Figure 1). For the first few years from the beginning of the online system activation, low numbers of scrub typhus cases have been reported in literature. Nonetheless, the annual peak of scrub typhus incidence remained constant, and the incidence tended to increase up to 2019, with a few exceptions in some years (2016–2017). Between epidemiological weeks 43 (midOctober) and 47 (midNovember), weekly cases reached the 50th percentile or above in most regions, independent of latitude. Meteorological variables such as the minimum, mean, and maximum temperatures, mRH, and weekly cumulative precipitation showed seasonal trends, with maximum observation values in the summer (June, July, and August) and minimum values observed in the winter (December, January, and February). The magnitudes of the seasonal variation for weekly cumulative sunshine duration, DTR, and MWS were comparatively small (Table S1, Figure 1, Figure S1).
Descriptive statistics
Weekly cumulative scrub typhus cases peaked in the epidemiologic weeks between 39 and 52/53 week (mean: 27 cases, IQR: [2, 32]). In addition, the weekly case reporting was less than one, if averaged. The medians of the minimum, mean, and maximum temperatures (°C) across the country were 20.8°C (IQR: [18.5, 23.0]), 24.0°C (IQR: [22.2, 26.1]), and 28.5°C (IQR: [26.7, 30.5]) in the summer (from June to October, epidemiologic weeks 22 to 35). In winter, these values were 2.9°C (IQR: [5.6, 0.1]), 1.5°C (IQR: [1.2, 4.2]), and 6.4°C (IQR: [3.7, 9.3]), respectively (from December to February, epidemiologic weeks 1–8 and 48–52). The DTR during winter was slightly wider than that during summer (median: 9.3°C vs. 7.8°C). The mean RH during the summer reached 98.6% (4th quartile), whereas it decreased to 22.8% (1st quartile) during winter (Table 1). Likewise, the weekly cumulative precipitation, especially during the monsoon season in the summer, rose to 695.5 mm (4th quartile). The weekly cumulative precipitation in winter was approximately oneseventh of the value in summer (mean: 7.0 mm vs. 47.7 mm). Weekly cumulative sunshine duration and MWS had the least variance compared with other meteorological variables (mean ± SD in summer vs. winter: 35.5 ± 19.1 vs. 35.4 ± 14.6 for weekly cumulative sunshine duration; 2.1 ± 0.7 vs. 2.3 ± 1.0 for MWS).
Province 
Mean temperature (°C) 
Mean relative humidity (%) 
Scrub typhus cases (cases per week) 


Min 
Mean ± SD 
Max 
Min 
Mean ± SD 
Max 
Min 
Mean ± SD 
Max 

Seoul 
10.1 
13.0 ± 10.3 
32.3 
25.6 
60.2 ± 11.0 
95.6 
0 
4.0 ± 10.4 
79 

Busan 
1.8 
15.1 ± 7.9 
31.7 
22.8 
62.3 ± 15.1 
98.6 
0 
9.7 ± 27.5 
200 

Daegu 
4.8 
14.6 ± 9.3 
31.7 
26.1 
58.3 ± 13.1 
97.1 
0 
3.8 ± 12.4 
161 

Incheon 
9.7 
12.5 ± 9.6 
29.5 
36.5 
69.5 ± 11.7 
98.5 
0 
1.3 ± 3.5 
35 

Gwangju 
4.1 
14.4 ± 9.2 
30.8 
33.1 
67.8 ± 10.5 
98.0 
0 
4.2 ± 12.0 
117 

Daejeon 
7.2 
13.3 ± 9.8 
31.8 
34.1 
67.4 ± 11.0 
96.4 
0 
4.7 ± 14.2 
118 

Ulsan 
3.4 
14.6 ± 8.4 
31.1 
23.6 
63.5 ± 14.9 
97.9 
0 
5.5 ± 17.2 
164 

Kyeonggido (KK) 
10.4 
11.8 ± 10.4 
30.7 
30.9 
67.5 ± 10.6 
96.6 
0 
10.6 ± 29.0 
223 

Kangwondo (KW) 
9.4 
11.3 ± 9.7 
28.9 
30.0 
65.9 ± 11.8 
94.4 
0 
3.2 ± 12.3 
168 

Chungcheongbukdo (CB) 
8.8 
11.9 ± 10.1 
29.9 
29.7 
65.9 ± 10.4 
90.8 
0 
4.5 ± 13.6 
104 

Chungcheongnamdo (CN) 
7.4 
12.4 ± 9.7 
30.0 
38.9 
71.3 ± 8.6 
93.8 
0 
14.2 ± 39.6 
281 

Jeollabukdo (JB) 
6.6 
12.8 ± 9.6 
29.6 
49.1 
71.6 ± 8.1 
93.3 
0 
14.9 ± 40.1 
287 

Jeollanamdo (JN) 
3.5 
13.7 ± 8.6 
29.9 
44.3 
72.6 ± 9.3 
93.4 
0 
16.4 ± 42.0 
396 

Gyeongsangbukdo (GB) 
6.0 
12.6 ± 9.2 
29.6 
30.5 
65.6 ± 12.1 
90.6 
0 
8.6 ± 24.2 
177 

Gyeongsangnamdo (GN) 
3.3 
14.0 ± 8.8 
31.0 
35.8 
65.7 ± 11.8 
96.6 
0 
19.3 ± 53.6 
562 

Jejudo (JJ) 
1.8 
16.2 ± 7.3 
30.0 
44.0 
71.6 ± 10.2 
96.1 
0 
1.1 ± 3.0 
25 
2. Association between climate factors and scrub typhus incidence (threshold analysis)
Weekly scrub typhus incidence was correlated with weekly mean temperature, mRH, and weekly cumulative precipitation with lags of 10–15 weeks, and with Pearson’s correlation coefficients ranging from 0.22 to 0.44, 0.22, 0.28, and 0.08 to 0.28, respectively. DTR, weekly cumulative sunshine duration, and MWS were correlated with weekly scrub typhus incidence, with a lag of approximately 20 weeks. The linear association between scrub typhus incidence for each province and meteorological variables above or below the threshold point for all provinces from GAM are shown in Figure S2.
The GAM for the association between mean AT or mean RH and weekly scrub typhus incidence per 100,000 persons in Busan is shown in Figure 2. Association plots between climate factors (i.e., mean AT and mean RH) and weekly scrub typhus incidence (per 100,000 persons) with a simple lag of 15 weeks showed a nonlinear relationship. In fact, mean AT and mean RH had a delayed impact (15week lag) on the increase in scrub typhus incidence above a threshold point (Figure S6, S7). From Seoul to Jejudo, the thresholds of mean ATs ranged from 14.9°C (Jejudo) to 17.4°C (Incheon). The threshold points for mean RH for each province ranged from 57.3% (GB) to 67.7% (JN); however, in 9 out of 16 provinces, we observed a linear association between the mean RH and scrub typhus incidence (Table 2).
Mean temperature (°C) 
Mean relative humidity (%) 


Threshold 
Beta coefficient (95% CI) 
pvalue 
Threshold 
Beta coefficient (95% CI) 
pvalue 

Seoul 
17.0 
0.117 (0.085, 0.150) 
< 0.001 
NA 
0.056 (0.043, 0.069) 
< 0.001 
Busan 
16.7 
0.204 (0.154, 0.255) 
< 0.001 
NA 
0.058 (0.045, 0.070) 
< 0.001 
Daegu 
16.1 
0.309 (0.139, 0.480) 
< 0.001 
NA 
0.064 (0.046, 0.082) 
< 0.001 
Incheon 
17.4 
0.103 (0.061, 0.144) 
< 0.001 
69.6 
0.066 (0.048, 0.084) 
< 0.001 
Gwangju 
16.2 
0.225 (0.124, 0.326) 
< 0.001 
NA 
0.054 (0.032, 0.077) 
0.002 
Daejeon 
15.2 
0.254 (0.136, 0.372) 
< 0.001 
NA 
0.050 (0.028, 0.072) 
< 0.001 
Ulsan 
14.8 
0.283 (0.213, 0.353 
< 0.001 
NA 
0.041 (0.024, 0.058) 
0.539 
GG 
15.7 
0.151 (0.117, 0.185) 
< 0.001 
55.9 
0.042 (0.027, 0.057) 
< 0.001 
GW 
16.1 
0.080 (0.054, 0.107) 
< 0.001 
56.8 
0.015 (0.001, 0.030) 
0.070 
CB 
15.0 
0.279 (0.086, 0.472) 
< 0.001 
NA 
0.035 (0.007, 0.077) 
0.102 
CN 
15.5 
0.242 (0.171, 0.313) 
< 0.001 
NA 
0.032 (0.007, 0.058) 
0.013 
JB 
16.1 
0.240 (0.172, 0.308) 
< 0.001 
NA 
0.016 (0.010, 0.043) 
< 0.225 
JN 
16.3 
0.105 (0.081, 0.129) 
< 0.001 
67.7 
0.045 (0.029, 0.061) 
< 0.001 
GB 
15.0 
0.229 (0.167, 0.292) 
< 0.001 
57.3 
0.043 (0.024, 0.062) 
0.002 
GN 
16.2 
0.172 (0.136, 0.208) 
< 0.001 
57.5 
0.042 (0.026, 0.058) 
< 0.001 
JJ 
14.9 
0.196 (0.136, 0.256) 
< 0.001 
57.8 
0.001 (0.021, 0.023) 
0.918 
Above a threshold point, weekly scrub typhus incidence (per 100,000 persons) constantly increased as the mean AT or mean RH increased. This tendency was consistently observed in all the provinces (Figure S6 and S7). GLM adopted for the linear section above a threshold mean AT resulted in a statistically significant association (beta, the slope of a linear section) for all provinces (Table 2). The beta coefficient for a 1°C increase in mean AT ranged from 0.103 (Incheon) to 0.309 (Daegu). Similarly, the beta coefficient of a 1% increase in mean RH ranged from 0.001 (Jeju, 95% CI: 0.021, 0.023) to 0.066 (Incheon, 95% CI: 0.048, 0.084) (Table 2). When pooled by metaanalysis, percent change in weekly scrub typhus incidence per 100,000 persons with a simple lag of 15 weeks per 0.1°C increase in mean AT was 1.85% (95% CI: 1.50, 2.20) across the country. Correspondingly, the percent change in weekly scrub typhus incidence per 100,000 persons with a simple lag of 15 weeks per 1% increase in mean RH in the country was 4.33% (95% CI: 3.41, 5.26) by pooling provincial data in the same manner (Figure 3A, 3C).
For sensitivity analysis, to smoothen the association plots between mean AT or mean RH and scrub typhus incidence, we adjusted for seasons by applying harmonics. To view the robustness of the applied GAM formula, we drew different GAM plots with or without different harmonics (2, 4, 6, 8, and 10). Adopting harmonics showed a smoother association. When the harmonics were given a value of six, the smoothed GAM plots remained unchanged, which supported the robustness of the model (Figure S8). After the simple lag was confined to 15 weeks, GAM and GLM were repeatedly performed by replacing the df values for a natural cubic spline with different values (Figure S9). Different df values in the natural cubic spline estimated a nearly unchanged beta coefficient. To estimate an increase in scrub typhus incidence following a unit increase in mean AT or mean RH, the GLM adopted covariates from several meteorological variables. Model I included two meteorological variables (mean AT or mean RH reciprocally and weekly cumulative precipitation). Model II was additionally adjusted for weekly sunshine duration and mean wind speed. The two models resulted in similar outcomes in the estimation of scrub typhus incidence (Figure 3B, 3D).
3. Longterm association between AT change and scrub typhus incidence
Compared with the previous 30 years, the annual mean AT difference tended to increase between 2001 and 2019. The number of heatwave days per year increased in a similar manner. When visualized in fiveyear terms, the longterm mean AT difference showed an increasing trend. In Figure 5, we plotted the longterm trends of meteorological variables (e.g., annual mean AT, annual mean RH), differences in annual mean AT or yearly total number of heatwave days, and cumulative scrub typhus incidence per year. Through LOESS smoothing of the trends, the plots showed that weekly sunshine duration and the number of heatwave days had gradually increased. Cumulative scrub typhus incidence, annual mean AT, and the different values of mean AT or the number of heatwave days slightly increased (Figure 4 and 5). The fitted linear regression model estimated an annual mean AT elevation of 1.17°C between 2001 and 2019 (95% CI: 1.02°C, 1.33°C). Likewise, the number of heat wave days per year compared with the past 30 years increased by approximately 2.82 days (95% CI: 2.04, 3.60). Simple linear regression analysis revealed that “year” had a stronger association than the other two variables (difference in annual mean AT or the number of heatwave days per year) with total outbreak cases. In a simple linear regression, a 1°C elevation of annual mean AT compared with the previous 30 years was associated with an additional 453.7 total outbreak cases (95% CI: 400.6, 506.8). In addition, a oneday increase in the number of heatwave days per year over the first 19 years (2001–2019) was associated with 429.5 more total outbreak cases (95% CI: 378.7, 480.4). In GLME, an annual mean AT elevation of 1°C explained a 36.9% increase in annual cumulative scrub typhus incidence per 100,000 persons in a country. In addition, one more heatwave day in a year was associated with a 2.5% increase in annual cumulative scrub typhus incidence per 100,000 persons at the country level.
Weekly scrub typhus incidence was associated with mean AT and mean RH with a simple lag of 15 weeks in the ROK. This association was nonlinear, namely Ushape or Vshape; thus, we approximated a threshold point above or below which scrub typhus incidence was linearly associated with mean AT or mean RH. The difference in the annual mean AT and the annual number of heatwave days in the period 1971–2000 increased over time from 2001 to 2019, and the annual cumulative scrub typhus incidence rate (per 100,000) was positively associated with the difference in annual mean AT or annual number of heatwave days in the past 30 years (1971–2000).
The Infectious Disease Control and Prevention Act in Korea has classified scrub typhus as a class 2 infectious disease that requires mandatory reporting since 1994. All confirmed or clinically suspected scrub typhus cases are reported to local health bureaus and then to the Korea Centre for Disease Control and Prevention ^{23}. This means that the total number of reported cases to the health agency matched our target population. The ROK allows a single national health insurance service to be operated by a public agency under the Ministry of Health and Welfare (MOHW). Health claims data might provide corresponding information on scrub typhus incidence, considering that the Korea Standard Code of Disease 10 (KCD 10) was used when the claim data were initially collected. Compared with the national mandatory reporting system of infectious diseases, health claims data relies on health professionals’ inputs, mostly based on clinical diagnosis. Therefore, serologically unidentified clinical cases could be included in the health claims data. In comparison, surveillance data comprised of clinically confirmed cases, which did not allow for the overestimation of the real numbers. Hence, for clarity and completeness of the data, we used surveillance data for scrub typhus from the National Infectious Disease Portal.
The Responsible Division for Zoonotic Infectious Diseases, previously denoted as the Division of Vectors and Parasitic Diseases of the Korea National Institute of Health (KNIH), has annually conducted vector surveillance for scrub typhus in spring and autumn. This was conducted at evenly distributed locations by reporting the ratio of infected chiggers ^{24–26}. At least 32 species of trombiculid mites transfer Orientia tsutsugamushi to hosts, of which nine representative Leptotrombidium species distributed across the AsiaPacific region, are mainly reported in the ROK and Japan ^{27,28}. The detection of infected chiggers and case reports in Seoul indicate that scrub typhus has already been urbanized across the country. Seven Trombiculidae species are domestically wellknown vectors, with the dominant species being L. pallidum in the northeastern provinces and L. scutellare in the southwestern provinces. The reason for why the nonlinear curves obtained from GAM for mean AT or mean RH differed for each province (Figure S5 and S6) might be due to the proportion of dominant vector species and the difference in dominant vectors’ favorable climate conditions.
Consistent with previous studies on the incidence of mosquitoborne diseases ^{29–31}, scrub typhus incidence has a nonlinear association with mean AT or mean RH in the ROK ^{32}. For scrub typhus, ranges of mean AT or mean RH, in which weekly scrub typhus incidence increased accordingly, were above specific threshold points ^{30}. Modelling approaches such as transmission rate calculations using temperaturedependent entomological parameters and mosquito biting rates proved that there was an upper limit of optimum temperatures ^{33}. Furthermore, other infectious disease incidences (e.g., herpes zoster) or elderly mortality had a similar correlation pattern. The incidence of infectious disease or chronic disease incidence increased linearly ^{34,35}.
AT affects the development of Leptotrombidium deliense (a vector species of scrub typhus) during its life cycle. Previous studies have found critical temperatures for the development and reproduction of L. delicense, as well as optimum AT ranges for larval activity and disease transmission. Under suitable AT conditions, egg hatching and chigger development speed up as AT increases ^{12,36−38}.
Warm humid summers enhance the transmission of diseases such as Lyme borreliosis or tickborne encephalitis via an increased number of vectors and their exposure to humans ^{39–41}. High summer temperatures are likely to restrict larval activity in midsummer and promote nymphal activity in late autumn or early spring of the following year ^{42,43}. In addition, low humidity depresses nymphal and adult tick activity ^{42,44,45}. Our study, which is associated with summer meteorological factors (i.e., mean AT and mean RH) with the later incidence of scrub typhus in autumn, provides supporting evidence for previous study findings or suggestions accordingly.
One of our study’s main outcomes was the association between meteorological factors (mean AT and mean RH) and scrub typhus incidence. Previous studies assumed a linear relationship between scrub typhus incidence and meteorological factors to estimate the increase in scrub typhus incidence per unit change in each meteorological variable ^{46}. To the best of our knowledge, this study is the first to approximate the threshold points for U or Vshaped nonlinear associations. Furthermore, the annual cumulative incidence of scrub typhus showed a linear relationship with the difference between the annual mean AT or the annual number of heatwave days and those of the period 1971–2000.
Our study results are in line with previous experiences regarding the impact of climate change, especially for global warming ^{2,47−51}. The countrywide annual mean AT gradually increased over the past 19 years, from 2001 to 2019. Extreme hot weather, depicted by the number of heatwave days, increased in parallel. We intended to prove that warmer temperatures or extreme hot summer weather influenced the outbreak curve in the following autumn season, which was quantified by the cumulative outbreak cases. Milder and shorter winters are known to induce longer tickactivity seasons ^{52}. In contrast, an extremely cold winter in Germany in 2012 was followed by a significantly low rate of hostseeking activity in Ixodes ricinus nymphs ^{53}. Retrospective data from Europe proves that the development and activity of ticks in summer are influenced by high temperatures, and global warming is presumably related ^{40}. Similarly, extreme weather conditions (i.e., extreme cold or hot temperature events) are related to the increase in noncommunicable diseases ^{54} and to an increase in the incidence of chronic infectious diseases ^{55}.
To provide an early alarm to the public as well as the public health sector, predictive models and projectile tools such as the Early Warning and Response System ^{56–60} should be created in the near future ^{56–60}. In addition, we may need to thoroughly acknowledge the interactions between environments, pathogens, vectors, and hosts in order to prepare preemptively. A strength of our study is that it provides useful data for future research.
The goal of this study was to relate components of an ecosystem, which consists of environmental factors, pathogens, vectors, and hosts. To do so, we repeatedly ran GAMs by applying several different models, including simple and moving average lags in weeks. We might be able to better demonstrate the real situation if the cumulative effects of climate factors preceding outbreaks were combined accordingly. However, the increased infection risk of scrub typhus due to changes in climate factors can only be quantified by drawing on the association between them.
We approximated the increased incidence of scrub typhus in the ROK in recent years (2001–2019) when mean AT, mean RH, annual mean AT, and annual number of heatwave days increased. Compared with the previous 30 years, the annual mean AT increased by just 1°C from 2001 to 2019. The number of heatwave days also increased by almost three days during the same period. Global warming has influenced the population dynamics of arthropod vectors and the transmission dynamics of the diseases they carry ^{61}. Our results support this idea by presenting an increased incidence rate. Longterm trends consistently showed a positive association between mean AT or heatwave days and scrub typhus incidence. Future climate change is likely to extend the transmission period of vectorborne diseases and expand the geographical range of vector habitats, which may cause new emergence or reemergence of the diseases ^{62}. Further studies to estimate the impact of climate variability on the incidence of vectorborne diseases, such as scrub typhus, are needed, as some of these diseases might threaten vulnerable populations at risk or continue to influence human lives ^{63}.
Study population
The number of weekly scrub typhus cases was obtained from the National Infectious Disease Portal, provided by the Korea Disease Control and Prevention Agency (KDCA), and the period of interest ranged from 2001 to 2019 since data acquisition from this platform was only available from 2001. Surveillance system of scrub typhus was based on mandatory reports of all suspected and confirmed cases districtwise. We obtained weekly reported suspected and confirmed cases in seven metropolitan cities and nine provinces. Suspected cases were defined by clinical manifestations of scrub typhus, such as acute febrile illness, lymphadenopathy and/or skin eschar, and epidemiological association. Confirmed cases were defined based on clinical symptoms and laboratory confirmation. The year 2020 was excluded to avoid the influence of the coronavirus disease in 2019. Our study only used publicly available data; therefore, the requirement for informed consent was waived off (IRB No. 20031311110).
Meteorological data acquisition
We obtained meteorological data from the Weather Data Open Portal provided by the Korea Meteorological Administration. Meteorological data such as AT, RH, precipitation, wind speed, and sunshine duration on a minutetominute basis were collected from a total of 102 monitoring stations across the country. We computed the daily average values of these meteorological variables by city and province.
We obtained information on the daily mean, maximum, and minimum AT, mean RH, mean wind speed, sunshine duration, and precipitation. The diurnal temperature range (DTR) was defined as the difference between the daily maximum and minimum temperatures. We then computed the weekly mean AT (mean, minimum, and maximum), mean RH, DTR, mean wind speed, and weekly cumulative sunshine duration and weekly cumulative precipitation.
We also computed the difference between the annual mean AT each year (2001–2019) and the mean AT for the period 1971–2000. Likewise, we computed the difference between the number of heatwave days each year (2001–2019) and the mean number of heatwave days during the years 1971–2000. The number of heatwave days per year was defined by the number of days in which the daily maximum AT was ≥33°C.
Statistical analysis
1. Descriptive analysis
We adopted locally estimated scatterplot smoothing (LOESS) curve fitting to determine whether scrub typhus incidence showed seasonality ^{64,65}. Time series plots of climate factors and weekly incidence of scrub typhus from January 2001 to December 2019 are shown in Figure 1. We then conducted a bivariate analysis between the meteorological variables and weekly scrub typhus incidence. As there were lags of up to 7 months between the annual peaks of the mean AT and the peaks of scrub typhus incidence, we decided to use lags from 0 to 28 weeks (Figure S1).
The correlation between each meteorological variable (mean, maximum, and minimum AT, mean RH, mean wind speed, weekly cumulative sunshine duration, weekly cumulative precipitation, and DTR) and weekly scrub typhus incidence was explored by applying a generalized additive model (GAM). Minimum, mean, and maximum AT, mean RH, and weekly cumulative precipitation were correlated with weekly scrub typhus incidence. Since the minimum, mean, and maximum AT were correlated, we only selected mean RH and weekly cumulative precipitation as covariates for GAM.
We constructed GAM for the association between AT and weekly scrub typhus incidence by adjusting for mean RH and precipitation, using a quasiPoisson function, and simple lags and moving average lags from 0 to 28 weeks. Meteorological variables and weekly scrub typhus incidence showed seasonality between 2001 and 2019 (Figure 1 and S1). According to previous literature, Fourier terms up to six harmonics per year were included to reflect seasonality in the GAM ^{66,67}. The GAM formula used is as follows:
Ln [E(Y)] = β_{0} + β_{1}X_{t} + Σ S_{i} (X_{j}) + sin (2πt_{j} /N) + cos (2πt_{j}/N) + year,
where E(Y) denotes expected scrub typhus incidence; β_{1} is the coefficient (slope) of a meteorological variable for week t; S_{i}(X_{j}) is a smoothing function for an individual covariate i; t_{j} indicates the week of the year j (j = 1, 2, ... ,52); N is a harmonic for the year (N = 1, 2, ..., 6); and year depicts the year when the scrub typhus cases were reported on an ordinal scale.
Associations between AT and weekly scrub typhus incidence were explored for the seven cities and nine provinces with different lags, from 10 to 15 weeks (Figure S5). Association plots obtained by applying GAM showed that there was a nonlinear relationship between a meteorological factor 10–15 weeks ahead and weekly scrub typhus incidence (Figure S5). Compared with other lags, such as simple lags of 10–14 weeks or moving average lags of 10–15 weeks, GAM plots with a simple lag of 15 weeks consistently showed V or U shape in most of the cities or provinces (Figure S5).
2. The association between weekly AT and scrub typhus incidence (threshold analysis)
To estimate the association between the increase of scrub typhus incidence and AT or RH, which showed “V” or “U” shape, we attempted to find a threshold where the slope of the graph changes from negative to positive (Figure S6 and S7). The R HEAT package uses a generalized linear model to approximate a TP if the beta coefficient (slope) is converted from negative (positive) to positive (negative) ^{68}. In particular, the threshold function in the HEAT package uses piecewise linear regression to approximate a threshold. For instance, β_{1} was calculated below a threshold, where weekly scrub typhus incidence decreased as mean AT increased and β_{2} was calculated above a threshold, where scrub typhus incidence increased as mean AT increased. After visual approximation of a threshold for mean AT and mean RH from GAM, we applied the threshold () function to quantify the percent changes in scrub typhus incidence following a unit change of each meteorological variable above or below the threshold. The GLM formula used is as follows:
Ln [E(Y)] = β_{0} + β_{1}X_{t} + Σ S_{i} (X_{j}) + sin (2πt_{j} /N) + cos (2πt_{j}/N) + year, if X_{t} <TP.
Ln [E(Y)] = β_{0} + β_{2}X_{t} + Σ S_{i} (X_{j}) + sin (2πt_{j} /N) + cos (2πt_{j}/N) + year, if X_{t} ≥TP,
where E(Y) denotes the expected scrub typhus incidence; β_{1} is the coefficient (slope) of a meteorological variable for week t below the threshold point (TP) of X_{t} and β_{2} is the coefficient (slope) of a meteorological variable for week t above the threshold point (TP) of X_{t}; S_{i}(X_{j}) is a smoothing function for an individual covariate i; t_{j} indicates the week of the year j (j = 1, 2, ... ,52); N is a harmonic for the year (N = 1, 2, ..., 6); and year depicts the year when scrub typhus cases were reported.
Since data distributions for meteorological variables and scrub typhus incidence were skewed, we adopted quasiPoisson regression in the GAM and GLM models. The smoothing parameter was determined by natural cubic spline by selecting an optimal degree of freedom ^{60}, where the unbiased risk estimator (UBRE) ^{69,70} in the GAM and Akaike information criterion (AIC) ^{71,72} values in the GLM for each province approached their minimum. Considering all provinces’ UBRE values following increasing degrees of freedom (df), our GAM model used 5 as its df.
For the sensitivity analysis, we used different dfs for each city and province. We also compared GAM with and without harmonics. Harmonics of 2, 4, 6, 8, and 10 were used for further comparison. Furthermore, we compared the GLM by additionally adjusting for mean wind speed, weekly precipitation, and weekly sunshine duration.
3. Longterm association between meteorological factors and scrub typhus incidence
We investigated the longterm association between meteorological factors and scrub typhus incidence using annual cumulative scrub typhus cases. The longterm mean AT difference was defined as the difference between the mean annual AT each year (2001–2019) and the mean AT of the previous 30 years (1971–2000). The difference in the number of heat wave days was defined as the difference between the number of heat wave days each year (2001–2019) and the mean number of heat wave days during the years 1971–2000.
We plotted the longterm trends of meteorological variables from 2001 to 2019 using the LOESS regression. The difference in annual mean AT or annual number of heatwave days from the mean AT or mean number of heatwave days between 1971 and 2000 was grouped by fiveyear terms to detect longterm trends (Figure 4A, 4B). A simple linear regression model overlaid the point plots to determine if a linear correlation existed between the variables. The fitted model slope, pvalue, and coefficient of determination (R^{2}) were calculated accordingly. Although the time variable year best explained increases in cumulative scrub typhus cases in autumn outbreaks in the model selection process, we intended to correlate cumulative scrub typhus cases in autumn outbreaks with differences in the annual mean AT or the number of heatwave days from those of the years from 1971 to 2000. Therefore, point plots, linear regression lines, 95% confidence intervals, slopes, pvalues, and R^{2} values were obtained in the same manner (Figure 4).
To analyze the association between the difference in annual mean AT of each year and mean AT of the past 30 years (1971–2000), we used a generalized linear mixedeffects model (GLME). Covariates were selected from the directed acyclic graph (DAG) (www.dagitty.net) (Figure S10). We considered changes in population and land use or greenness over the last 19 years; however, we did not include them as covariates since population changes were not directly associated with AT, and land use was a causal mediator rather than a confounder (Figure S10). Instead, we used the incidence rate (annual incidence of scrub typhus per 100,000 population) to consider population factors in the regression model. We put the province as a random intercept in the glmer() to diminish spatial interference in the estimation. The analysis was repeated after replacing the difference in annual mean AT by the difference in the number of heatwave days.
All statistical analyses were conducted using the statistical software packages SAS (version 9.4) and RStudio software (ver. 1.2.5042).
Ethical statement
All the methods used in this study followed relevant guidelines and regulations. Our study only used publicly available data; therefore, the requirement for informed consent was waived off.
Data availability
The dataset generated and/or analyzed during the current study are available from the corresponding author on reasonable request.
Author contributions
D.S. and Y.C. contributed equally as first authors. D.S. and Y.C. collected and analyzed data, and wrote the first draft. All authors contributed to the conceptualization of study design, writing and approval of the final version of manuscript.
Additional Information
The authors declare no competing interests.
Funding
This research did not receive any specific grant from funding agencies in the public, commercial, or notforprofit sectors.