The estimated burden of scrub typhus in Thailand from national surveillance data (2003-2018)

Background Scrub typhus is a major cause of acute febrile illness in the tropics and is endemic over large areas of the Asia Pacific region. The national and global burden of scrub typhus remains unclear due to limited data and difficulties surrounding diagnosis. Methodology/Principal findings Scrub typhus reporting data from 2003–2018 were collected from the Thai national disease surveillance system. Additional information including the district, sub-district and village of residence, population, geographical, meteorological and satellite imagery data were also collected for Chiangrai, the province with the highest number of reported cases from 2003–2018. From 2003–2018, 103,345 cases of scrub typhus were reported with the number of reported cases increasing substantially over the observed period. There were more men than women, with agricultural workers the main occupational group affected. The majority of cases occurred in the 15–64 year old age group (72,144/99,543, 72%). Disease burden was greatest in the northern region, accounting for 53% of the total reported cases per year (mean). In the northern region, five provinces–Chiangrai, Chiangmai, Tak, Nan and Mae Hong Son–accounted for 84% (46,927/55,872) of the total cases from the northern region or 45% (46,927/103,345) of cases nationally. The majority of cases occurred from June to November but seasonality was less marked in the southern region. In Chiangrai province, elevation, rainfall, temperature, population size, habitat complexity and diversity of land cover contributed to scrub typhus incidence. Interpretation The burden of scrub typhus in Thailand is high with disease incidence rising significantly over the last two decades. However, disease burden is not uniform with northern provinces particularly affected. Agricultural activity along with geographical, meteorological and land cover factors are likely to contribute to disease incidence. Our report, along with existing epidemiological data, suggests that scrub typhus is the most clinically important rickettsial disease globally.

Introduction Scrub typhus in Asia, Australia and the Islands of the Indian and Pacific Oceans is a major cause of acute non-malarial febrile illness and a potentially severe but treatable zoonotic disease caused by the obligate intracellular Gram negative bacterium Orientia tsutsugamushi [1][2][3][4]. The disease is prevalent across an area of at least 13,000,000 km 2 , encompassing large tracts of the Asia Pacific region, with emerging reports from the Middle East, South America, Africa and Europe suggesting the disease is more widespread than previously thought [5][6][7][8][9]. Median mortality rates have been estimated at 6.0% and 1.4% for untreated and treated disease, respectively, with a wide range of 0-70% depending on the study and the region [10,11]. Transmission occurs through the bite of an infected larval stage (chigger) of the trombiculid mites [12]. Disease incidence may be related to seasonal variation in temperature, rainfall and humidity; socioeconomic factors, terrain, agricultural or other activities exposing humans to suitable habitat or transitional land cover type (mixed cropland and wild vegetation) favouring increased rodent and vector chigger density [13][14][15][16]. Although the disease is more prevalent in rural regions, cases from urban centres continue to be reported [17].
The incubation period ranges from 6 to 21 days and clinically, scrub typhus is characterised by fever, eschar, maculo-papular rash, headache, cough, nausea, vomiting, myalgia and lymphadenopathy [18]. The eschar at the site of the mite bite is the most characteristic feature but is not always present, being dependent on the degree of pre-existing immunity to the infecting Orientia tsutsugamushi strain [19,20]. Severe disease may present with complications affecting multiple organ systems [18]. The standard antibiotic treatment is doxycycline for 7 days with alternative options including tetracycline, chloramphenicol, azithromycin and rifampicin [21,22]. In the majority of patients, fever will clear within 48-72 hours of treatment initiation. Reports of scrub typhus infections poorly responsive to treatment in northern Thailand and the lack of effective preventative measures, mainly vaccines with long-term multistrain protection, remains a concern [23,24].
Efforts to estimate the global burden of scrub typhus are hindered by limited data and difficulties surrounding diagnosis [11]. The frequently quoted estimate of one million cases annually lacks supportive evidence [25]. In Thailand, the first human case of scrub typhus was reported in 1952 although captured records from World War II suggest that cases among Japanese troops occurred from 1943 to 1944 [26,27]. Subsequently, the disease has been shown to be highly prevalent and a major cause of acute febrile illness in the country [4,[28][29][30]. It is notifiable to the National Disease Surveillance System (R506), Bureau of Epidemiology (BoE), Department of Disease Control (DDC), Ministry of Public Health (MoPH), with government and private healthcare facilities in each province reporting to the public health office [31]. There are 76 provinces in Thailand (plus Bangkok, a special administrative area), grouped into 4 administrative regions (North, Northeast, Central and South), with each province being further divided into districts, sub-districts and villages. Within each province, the population is served by a central provincial hospital, district hospitals and primary care or health promoting units located in every sub-district. There has been an increase in notified cases reported from other countries with established passive notification systems such as South Korea, China and Bhutan, while results from Japan and more recently Taiwan suggest the disease incidence has stabilised [15,[32][33][34][35][36][37].
To date, limited data have been published on the national burden of scrub typhus in Thailand. In this study, we describe the burden of disease using national surveillance data from 2003-2018 and investigate the effect of geography, population, rainfall, temperature, habitat complexity and diversity of land use/land cover on disease incidence.

Ethics statement
Publically available disease reporting data were obtained from the National Disease Surveillance System (R506), Bureau of Epidemiology (BoE), Department of Disease Control (DDC), Ministry of Public Health (MoPH). All data analysed were anonymised.

Data collection
Detailed data for scrub typhus from 2003-2018 were obtained from the National Disease Surveillance System (R506) [31]. Brief annual summary reports from 1980-2002 were also available and collected. Probable and confirmed cases were reported according to the case definition, based on ICD-10: A75.3, and outlined in Table 1 below.
Cases were reported from governmental healthcare facilities including provincial hospitals, district hospitals and primary care units. Reporting from private healthcare facilities also occurred, albeit to a lesser extent. Reports for each year from 2003-2018 were obtained and collated into a single dataset. Information on the number of cases per province each month, population, annual incidence rate per 100,000 population (AIR), age groups, sex and occupation were collected. Additionally, less detailed information from annual summary reports from 1980-2002 were extracted and data collated into the main dataset. Regional administrative boundaries were obtained from the National Statistical Office, Ministry of Information and Communication Technology. As a case study, we obtained detailed reporting data for Chiangrai province from the Ministry of Public Health, the province with the highest number of reported cases from 2003-2018. Data included patient demographics, location and type of healthcare facilities, and residential data to the sub-district and village level. Administrative boundaries of the province including districts and sub-districts as well as population size were obtained from the Ministry of Natural Resources and Environment, Thailand (2012).
Average monthly temperature (˚C) and total monthly rainfall (mm) data for Chiangrai at a central representative weather station were obtained from the Thai Meteorological Department, Ministry of Information and Communication Technology. Land use information was obtained from the GlobCover 2009 with a resolution of 1 km (European Space Agency, http:// due.esrin.esa.int/page_globcover.php, classification given in S1 Document and S1 Fig). Data were obtained at different administrative levels for Chiangrai province: scrub typhus incidence at village level, population size at sub-district level, rainfall at provincial level and land use at 1 km resolution. We integrated all data at the sub-district level, aggregating human cases obtained at village level and cropping raster data (land use) using sub-district boundaries.

Statistical analyses
Proportions, percentages and averages (median and interquartile range [IQR] or mean and standard deviation [SD] as appropriate) were calculated controlling for any missing data. Seasonality was assessed by calculating proportions of cases (and 95% confidence intervals) reported during discrete time-periods and assessing for overlap as well as performing twosample test of proportions. These descriptive analyses were performed using STATA 15 software (College Station, Texas, USA). Spatio-temporal analyses were conducted using R software (R Core Team, 2018) [38]: i. Spatial autocorrelation of scrub typhus incidence and spatial interpolation using Gaussian process regression (kriging). ii. Temporal incidence and temporal autocorrelation of scrub typhus cases.
iii. Spatiotemporal distribution of incidence using autocorrelation method, wavelet analysis and correlation between scrub typhus cases, rainfall and temperature.
iv. Analysed the links between scrub typhus cases and the agro-environmental data set using general linear modelling (GLM).
v. Performed general additive modelling (GAM), taking into account the spatiotemporal dynamics obtained from I, II and III and starting with the initial model that retained the potential explanatory factors of IV.
Further details on spatio-temporal analyses are outlined in S2 Document. Maps were drawn using mapchart.net ©, tmap and tmaptools [39,40]. The months with the highest average number of cases (mean, SD) were October, July and September at 809 (361), 772 (315), and 757 (342) cases, respectively. The majority occurred between June and November, corresponding to the rainy and early cool or dry seasons in Thailand. The proportion of cases (95% CI) from June to November and December to May were calculated: 0.678 (0.675-0.681) and 0.322 (0.319-0.325), p<0.001. Detailed occupational data to the regional and provincial level were available from 2006-2018. In the North, agricultural workers made up the largest percentage of cases (median, IQR) at 40.35% (39.72-42.29), followed by labourers at 19.36% (18.79-20.94) and students at 16.45% (15.10-19.10). In the Northeast, there was a larger percentage of agricultural workers at 57.78% (52.82-60.93), followed by labourers and students at 11.31% (10.36-12.57) and 10.47% (9.56-13.79), respectively. A similar pattern was seen in the South, albeit to a lesser extent, with agricultural workers at 33.33% (31.67-33.70), labourers at 22.72% (20.74-25.37) and students at 15.71% (13.68-17.65). A difference was seen in the central region where
Regional differences were also observed in the seasonality of cases reported from 2003-2018 as shown in Fig 4. In all four regions, average monthly reported cases were at their lowest from February to April before rising in May. In the North and Northeast, cases peaked in July and October, respectively, before falling in November and December, reaching a trough from February to April. In the southern and central regions, the decline in cases at the end of the year was interrupted by another peak in January prior to falling to a nadir in April. Seasonality of scrub typhus cases was less marked in the southern region when compared to other regions. At the provincial level, the greatest burden of scrub typhus was in the northern provinces ( Fig 5).

Single province case study-Chiangrai province
Chiangrai province was chosen as it represented the province with the highest average number of scrub typhus cases per year along with the 4 th highest disease incidence rate per 100,000 population. The province is the northernmost province in Thailand and is divided into 18 districts, 124 sub-districts and 1,816 villages. I. Spatial autocorrelation, semi-variogram and kriging interpolation of scrub typhus cases. A weak spatial autocorrelation was observed when pooling scrub typhus cases for 2003-2018 for Chiangrai province but was no longer significant after 30 km (S3 Fig). The spatial distribution of the scrub typhus cases among sub-districts was analysed using semi-variogram analysis and depicted in Fig 6B. The kriging interpolation, using the results of the semi-

PLOS NEGLECTED TROPICAL DISEASES
variogram analysis, is represented in Fig 6C. Spatial interpolation of scrub typhus cases showed hotspots of cases in high elevation area of the province (Fig 6A and 6C).
II. Time series analysis of scrub typhus cases, rainfall and temperature. There was an increasing trend in reported scrub typhus cases per month for Chiangrai from 2003 to 2018

PLOS NEGLECTED TROPICAL DISEASES
The burden of scrub typhus in Thailand with a marked increase after 2010 (Fig 7). There was a strong seasonal pattern with large numbers of scrub typhus cases during the rainy season. The wavelet analysis (S4 Fig) confirmed a seasonal pattern starting in 2008. A strong seasonal pattern for both rainfall and temperature was also observed (Fig 7). A slight increase in average temperature from an average of 24.5˚C in 2003 to 25.5˚C in 2018 was observed, coinciding with the increase in reported cases. Wavelet power spectrum (S4 Fig) confirmed the above time series analyses by revealing significant 12-month periodicity over the entire study period.

IV. Association between scrub typhus cases and explanatory variables using General Linear Modelling (GLM).
A GLM analysis was performed using the following explanatory variables: habitat complexity, habitat fragmentation, forest cover (%), forest open cover (%), grassland open cover (%), mosaic habitat cover (%), rain-fed cropland cover (%), flooded-irrigated land cover (%), population size, rainfall and temperature (S2 Document). The top best model selected showed that scrub typhus cases per sub-district per month were significantly associated with all variables included in the initial model supported by an Akaike weights (w r ) value of 0.36 (S2 Document, Table 2).
Multi-collinearity assessed using Variance Inflation Factor (VIF) showed that most VIF values were inferior to 10, suggesting lack of collinearity, except for forest open cover with a VIF value of 22.55, suggesting high collinearity with other land use characteristics. The percentage of deviance explained by this model was estimated by maximum likelihood R 2 to 0.25 indicating a good prediction of this best model.
A positive association was found between population size and scrub typhus cases. The number of cases were also positively correlated with temperature and rainfall. Higher numbers of scrub typhus cases were found associated with habitat complexity characterised by large cover of open forested habitat. Lower numbers of scrub typhus cases were associated with fragmented habitat characterised by mosaic habitat or open grassland habitat. Negative associations were also observed between scrub typhus cases and forested, rain-fed and floodedirrigated lands. The above relationships are depicted in S6 Fig. V. Association between scrub typhus cases and explanatory variables using General Additive Modelling (GAM). An initial GAM model was developed that took into account the spatiotemporal dynamics of scrub typhus cases, the temporal dynamics of rainfall and temperature, and the potential explanatory variables selected by the GLM, with the exception of open forested habitats (due to high VIF value), using a negative binomial function (with theta Table 2 estimated as above). The best GAM model selected using Akaike Information Criteria (AIC) is shown in Table 3 and Fig 8, which did not include the variables habitat fragmentation and grassland open cover (%). The model explained 50.4% of the deviance (R 2 = 0.37) with a better performance than the above best GLM. The best GAM model confirmed the different results obtained above with a significant influence of geography (matrix of geographic coordinates of sub-districts), rainfall, temperature, population, habitat complexity, mosaic habitat, rain-fed and flooded-irrigated lands ( Table 3, Fig 8).  Table 3. Results of general additive modelling (GAM) explaining the number of cases of scrub typhus per sub-district in Chiangrai province using a negative binomial link (theta = 2.12), with approximate significance of smooth terms. For the best selected model, the deviance explained = 50.4%, R 2 = 0.37, maximum likelihood = 8868.2, AIC = 17689.8 (see Fig 8).

Discussion
The results of this study reveal a high burden of scrub typhus in Thailand that has increased significantly over the last two decades. Although cases have been reported from all 77 provinces from 2003-2018, the burden of disease is not uniform with a handful of northern provinces particularly affected. We have also shown that within Chiangrai province, scrub typhus is highly prevalent in the rural and mountainous areas to the western, northwestern and, to a lesser extent, eastern parts of the province.
Our results from Thailand are not unexpected when compared to other countries with reported surveillance data. There have been increases in scrub typhus cases reported in South Korea  [32][33][34][35][36]41]. Similar to Thailand, there were more cases in men than women in Japan and Taiwan but the opposite was seen in China and South Korea. In Thailand, the disease burden was highest in adults of working age while the proportion of children affected continues to fall, reflecting the accelerated increase in disease burden in adults compared to children. Scrub typhus was more common in older age groups in the other five countries, particularly Japan and South Korea with its ageing farming communities [42]. Unsurprisingly, high disease incidence among farmers in China, eastern Taiwan and South Korea were also described and in conjunction with our data, highlights the ongoing risk of disease in agricultural workers [32][33][34]. Increasing urbanisation of scrub typhus has recently been reported from Seoul, South Korea, which may either reflect increased human exposure to rural habitats as urban centres expand or the increase in development of parkland and gardens within municipal boundaries [17].
Seasonality of scrub typhus was observed in China, Taiwan, South Korea and Japan. In Japan, the disease incidence peaked in November and December with a smaller peak in May and June [36]. South Korea reported a disease peak from October to December with lower latitudes associated with a later peak in incidence [15,32]. In Taiwan, the disease incidence was biphasic with the main peak in June and July and a second smaller peak in September and October [41]. In mainland China, intra-country variation in seasonality was observed with cases in the northern provinces occurring mainly in autumn and early winter while in the southern provinces, the scrub typhus "season" began in spring and extended to the end of the year [33]. The disease pattern in southern China and Taiwan most closely resemble the results from northern and northeastern Thailand which is not unexpected given their similar latitudes. However, our study clearly demonstrates the change in disease seasonality in Thailand as latitude decreases. This effect was more pronounced than the study from South Korea and may reflect Thailand's proximity to the equator along with its greater latitude range [15]. In southern Thailand, there are two main seasons (wet and dry) and the temperature range is narrower than northern Thailand. Rainfall tend to be higher in the South, with the southwest monsoon affecting the west coast of the peninsula earlier in the season (April to October) while the east coast tends to have higher rainfall later (September to December). In higher latitudes, winters are frequently cold and summers hot, the temperature range affecting the population of vector chiggers [43]. This is demonstrated in Japan, where monthly incidence is more evenly distributed moving southwards [44]. In tropical regions, temperature varies less and here rainfall appears to influence chigger abundance and thus human risk [45,46].
In Chiangrai province, modelling revealed the partial contribution of geography (elevated or mountainous areas), rainfall, temperature, population size, habitat complexity and diversity of land cover to the number of reported scrub typhus cases per sub-district. The lag time of 2 months for temperature and 1 month for rainfall in the cross correlation analysis could be explained by agricultural activity, such as the clearing of scrub or new vegetation prior to planting in June/July, or chigger activity during the rainy season. Habitats characterised by growth of new or transitional vegetation and associated with the presence of infected rodents and chigger mites are important ecological factors [27,47]. These habitats are diverse and include virgin forests, gardens, fringe habitats or scrubland, water-meadows, beachheads, rice paddies, bamboo patches and oil palm or rubber plantations, many of which are present in Chiangrai. Although there was a significant effect of population on disease incidence, the proportion of the population exposed to the habitats described is likely to be more predictive but will be harder to obtain. The association of elevation and scrub typhus incidence in our study reflect findings from southern China [16]. This may reflect the reduced accessibility and greater habitat complexity of these regions. In addition, an increase in average monthly temperature from 24.5˚C in 2003 to 25.5˚C in 2018 correlated with an increase in scrub typhus incidence. This finding is also similar to reports from southern China and together, are suggestive of the influence of climate change on disease burden [16,48,49].
Diagnostic capabilities for scrub typhus at district hospitals are limited and most provincial hospitals are not able to perform confirmatory assays in-house, requiring the aid of reference laboratories for confirmation of diagnosis [50]. Antibody detection-based rapid diagnostic tests (RDTs) for scrub typhus have been introduced (e.g. since 2008 in Chiangrai province) but are sub-optimal [51]. The majority of scrub typhus cases reported in Thailand are probable cases based on the clinical criteria of fever and eschar (with or without positive RDT results). This combination of clinical findings has a sensitivity and specificity of 47% and 81% in adults and 60% and 92% in children for scrub typhus in Chiangrai [52,53]. Despite this weakness, a scrub typhus surveillance system based on a similar clinical criteria for probable cases can play a valuable role in allowing the disease burden to be extrapolated, particularly in regions where scrub typhus is the dominant eschar-associated disease (alternatives include spotted fever rickettsial infections or necrotic spider bite wounds).
Other limitations include the passive nature of the surveillance system and the effect of healthcare-seeking behaviour on the data used in this study. Passive surveillance systems are reliant on facilities and healthcare and public health staff for data completeness and quality. Although reporting data from private healthcare facilities are included in the analysis, it is unclear whether cases are consistently reported from the private sector. These limitations suggest that the burden of scrub typhus reported in this study is likely to be a gross underestimate of the true disease burden. Additionally, we were unable to measure the influence of other factors including the increase in disease awareness, availability of diagnostic tests, agricultural activity, human behaviour, and chigger abundance and activity in influencing the impressive rise in reported cases since 1980. Finally, although the data on patients' villages of residence were collected, accurate geolocation data proved particularly difficult to obtain which limited the utility of this high-resolution data. However, the data will allow public awareness and engagement programmes to target villages and sub-districts with the highest burden of scrub typhus. Paradoxically, despite the high incidence of scrub typhus in these areas, awareness among community health workers remains low and is almost completely absent in villagers.
In this report, we have estimated burden of scrub typhus in Thailand from national surveillance data for 2003-2018. Considering that currently available epidemiological data are limited, these results should contribute significantly towards our understanding of the burden of disease both locally and globally. Additional work on the epidemiology of scrub typhus and the development of accurate and cost-effective diagnostic tests should be prioritised. The single province case study has shown how geography, rainfall, temperature and landscape complexity may partly contribute to disease incidence. Further data collection and analyses to investigate the impact of geographical and meteorological factors on disease incidence for the remaining provinces are ongoing. The burden of scrub typhus from Thailand and other countries within and beyond the traditional endemic region suggest it is the most clinically important rickettsial disease globally.
Supporting information S1 Document. Land use-land cover classification for Chiangrai province. GlobCover 2009 with a resolution of 1 km was used (http://due.esrin.esa.int/page_globcover.php). for their help with data collection and permission to utilise the data. We thank staff at the Ministry of Natural Resources and Environment and Thai Meteorological Department, Ministry of Information and Communication Technology for their help with data collection. Additionally, we thank Nidanuch Tasak and Piangnet Jaiboon for their role in obtaining data for Chiangrai province.