# Scrub typhus in Nan province (Thailand): Seventeen years of data to understand the impact of land cover change

Nolwenn Blache, Karine Chalvet-Monfray, Rawadee Kumlert, Soawapak Hinjoy, Serge Morand

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Nolwenn Blache, Karine Chalvet-Monfray, Rawadee Kumlert, Soawapak Hinjoy, Serge Morand. Scrub typhus in Nan province (Thailand): Seventeen years of data to understand the impact of land cover change. PLoS Neglected Tropical Diseases, 2025, 19 (9), pp.e0013552. 10.1371/journal.pntd.0013552. hal-05271753v2

HAL Id: hal-05271753

https://hal.inrae.fr/hal-05271753v2

Submitted on 29 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH ARTICLE

# Scrub typhus in Nan province (Thailand): Seventeen years of data to understand the impact of land cover change

Nolwenn Blache<sup>1,2*</sup>, Karine Chalvet-Monfray<sup>2</sup>, Rawadee Kumlert<sup>3</sup>, Soawapak Hinjoy<sup>4</sup>, Serge Morand<sup>5,6,7</sup>

1 Université Clermont Auvergne, UMR Territoires, VetAgro Sup, Lempdes, France, 2 Université de Lyon, UMR EPIA, INRAE, VetAgro Sup, Marcy l'Etoile, France, 3 Division of Vector-Borne Diseases, Department of Disease Control, Ministry of Public Health, Nonthaburi, Thailand, 4 Department of Disease Control, Ministry of Public Health, Nonthaburi, Thailand, 5 IRL2021 HealthDEEP, CNRS – Kasetsart University – Mahidol University, Bangkok, Thailand, 6 Faculty of Veterinary Medicine, Kasetsart University, Bangkok, Thailand, 7 Faculty of Tropical Medicine, Mahidol University, Bangkok, Thailand

* nolwenn.blache@vetagro-sup.fr

<figure><img src="image_3.png" /></figure>

## Abstract

### Background

Scrub typhus, caused by *Orientia tsutsugamushi* and transmitted by chigger mites (*Leptotrombidium*), is a major health problem in northern Thailand, particularly in Nan province. Land cover change, by altering the ecosystem, could affect the ecology of the vector and consequently the risk of scrub typhus transmission.

### Methodology/principal findings

This study investigated the impact of land cover changes on scrub typhus transmission in 2.5 km buffer zones around each village of Nan Province between 2003 and 2019. Using the open land cover data of the European Spatial Agency Climate Change Initiative (ESA CCI), we quantified land cover composition and land cover changes and integrated public health data on scrub typhus cases, as well as information on elevation, population, and slope. Generalized Additive Models were applied to assess the effects of land cover changes on annual scrub typhus cases. Scrub typhus cases increased significantly during the study period, peaking in 2012 and 2016, mainly in mountainous areas rather than in the Nan River valley. Land cover associated with cases included shrubland, mosaic land, broadleaf forest, and needle-leaf forest. Cases increased with shrubland and mosaic land, displayed an inverted U-shaped relationship with broadleaf forest, and decreased with needleleaf forest. Key land cover change factors included shrubland transitions, population, and geographic interactions. Reforestation (from shrubland to broadleaf forest) showed an

<figure><img src="image_4.png" /></figure>

OPEN ACCESS

**Citation:** Blache N, Chalvet-Monfray K, Kumlert R, Hinjoy S, Morand S (2025) Scrub typhus in Nan province (Thailand): Seventeen years of data to understand the impact of land cover change. PLoS Negl Trop Dis 19(9): e0013552. https://doi.org/10.1371/journal.pntd.0013552

**Editor:** Husain Poonawala, Tufts Medical Center, UNITED STATES OF AMERICA

**Received:** February 18, 2025

**Accepted:** September 10, 2025

**Published:** September 18, 2025

**Copyright:** © 2025 Blache et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data availability statement:** All relevant data are within the manuscript and its Supporting Information files. Data on human scrub typhus cases from the National Disease Surveillance System are subject to ethical restrictions, as they contain sensitive information and cannot be shared in the supplementary files
or deposited in a public repository. However, these data are appropriately referenced within the article.

**Funding:** This study is supported by "rewild-4Health" project (ANR-24-CE35-5476) funded by the French National Agency for Research (ANR). SM is also supported by the project "Innovative Animal Health" of the Thailand International Cooperation Agency (TICA). NB is supported by VetAgro Sup (doctoral grant, decision by the general management following the scientific council meeting dated 14 June 2022) and Clermont Auvergne Metropole (decision dated 24 June 2022). The funders had no role in study design, either on data collection, analysis, decision to publish, or preparation of the manuscript.

**Competing interests:** The authors have declared that no competing interests exist.

inverted U-shaped relationship with cases, whereas stable broadleaf forest and loss of shrubland to grassland became non-significant. Male population increased cases.

## Conclusions/significance

This study highlights the importance of land cover changes in understanding disease transmission and suggests that landscapes disturbance may create optimal conditions for *O. tsutsugamushi* transmission. This is a novel regional-scale exploration of land cover impacts on scrub typhus in Thailand.

## Author summary

Scrub typhus, caused by the bacterium *Orientia tsutsugamushi*, is a vector-borne disease transmitted to humans exclusively by chigger mites. Despite one million annual cases, it remains a neglected tropical disease due to limited knowledge of vector ecology and transmission risk factors. Thailand is particularly affected by this disease where the pathogen actively circulates mainly in the northern region known as a rural endemic area. Focusing on Nan province, part of the endemic region, we investigated between 2003 and 2019, how the composition of land cover and land cover changes in the surroundings of the villages impact the number of scrub typhus human cases. Using satellite-derived land cover data and public health records, we found that important proportion of shrubland and mosaic land as well as an intermediate proportion of broadleaf forest were at risk for human health. When further investigating the role of land cover change, we discovered that areas with moderate reforestation from shrubland to forest increased scrub typhus cases, whereas stable forest cover from 2003–2019 showed no significant effect. These findings highlight the critical role of land cover and land cover changes in shaping disease risk, providing new insights into the ecological drivers of scrub typhus outbreaks.

## Introduction

Scrub typhus is a vector-borne disease transmitted to humans by chigger mites of the genus *Leptotrombidium*. The pathogen responsible for scrub typhus is a strictly intracellular bacteria called *Orientia tsutsugamushi* from the family of Rickettsiae. It is during the feeding process that the 6-leg larval stage vector also known as "chiggers", transmits the pathogen to humans. The disease manifests itself with flu-like symptoms, a skin eschar and death in 6% of the cases in the absence of antibiotic treatment [1]. A commonly cited estimate records one million new cases per year worldwide and potentially one billion people at threat [2]. The pressure is the highest in the endemic are known as the Tsutsugamushi Triangle encompassing Japan, Southeast Asia and Australia. Thailand is particularly impacted by the disease, with the Ministry of Public Health recording 103,345 cases between 2003 and 2017 [3]
and up to 77% of the villagers estimated to have been exposed to the pathogen in the worst-affected regions [4]. Half of the total scrub typhus cases reported between 2003 and 2018 occurred in the northern region. This area also included the five provinces with the highest annual cases number, making it an significant endemic area of Thailand [3].

Human cases of scrub typhus have been correlated with vector abundance in several studies [5]. Species richness of chigger mites and their infection rates with *Orientia tsutsugamushi* have also been identified as risk factors [6,7]. Chiggers can become infected by feeding on an infected host [8]. However, laboratory studies have suggested that they are primarily infected by transstadial and transovarial transmission, allowing the infection to be maintained over multiple generations [9]. *Leptotrombidium* mites spent 99% of their time living free in the environment and 1% feeding on vertebrate hosts. When free-living, chiggers require stable and suitable temperature and humidity conditions for their survival. They also need to encounter a host on which they can feed to complete their cycle. Early local studies in several regions of Southeast Asia identified chigger habitats as forested and shrubland areas, fallow land, neglected gardens or plantations, flooded meadows, hedgerows, marginal habitats and ecotones [10–12]. Secondary vegetation growth resulting from land cover change and favouring chiggers has also been identified as a risk factor in local studies [13]. Vectors may then benefit from the forest and associated tree cover because of the microclimate they provide and the small mammals and birds they support, making them a risky habitat if humans are exposed [14]. Few studies have investigated the impact of land cover and land cover change on scrub typhus at the regional scales, such as in China [15,16], Taiwan [17], South Korea [18] or Thailand [3]. However, to our knowledge, no studies have been conducted at local or village level which could be relevant given the low mobility of chiggers and the high heterogeneity of human cases within regions.

This study aimed to quantify the role of land cover and land cover change on scrub typhus cases in an endemic area of North Thailand. Instead of aggregating the scrub typhus cases at the subdistrict level, we used a 2500m radius buffer zone around each village of Nan Province that could encompass areas where villagers are mainly exposed during their agricultural activities. Based on previous literature, we hypothesized that 1) forest-related land covers have a higher burden of scrub typhus; 2) fast changing land covers such as secondary growth vegetation or fragmented landscapes increase the burden of scrub typhus compared with more stable land covers. Land cover and land cover transitions between 2003 and 2019 were quantified within these buffer zones. We combined passive public health surveillance data of scrub typhus cases with land cover, land cover transitions, and topographic data in two spatiotemporal models: one assessing the impact of land cover considering years and the other examining the effects of land cover transitions.

## Methods

### Study area

Nan Province (19.15° N, 100.83° E) is located in the northern part of Thailand (Fig 1A), bordered by Phayao and Phrae Provinces to the West, Uttaradit Province to the South and Lao PDR to the Northeast. Nan Province covers approximately 12,100 km², with 15 districts, 99 sub-districts, and 905 villages. The population reached 478,227 inhabitants in 2019 [19] and most of the employed people work in agriculture [20,21]. The Nan River flows through the province from north to south, creating a flat agrarian valley surrounded by mountainous and foothill areas (Fig 1B). The average temperature is 24.8°C with a minimum temperature of 16.4°C (January) and a maximum temperature of 33.8°C (recorded in April) for the period 1991–2021. The average annual rainfall for the same 30-year period was 1,490 mm, with most of the falling during the wet season from May to September. Nan Province is recognized as being severely affected by scrub typhus. In 2019, Nan Province had the second highest reported cases per 100,000 inhabitants of Thailand, reaching 114.71 [22].

### Data collection

We retrieved the number of scrub typhus human cases between 2003 and 2019 from the National Disease Surveillance System (Division of Epidemiology, Department of Disease Control, Ministry of Public Health of Thailand) [23]. The number of human scrub typhus cases was aggregated annually for each village in the province. The cumulated cases per village
<figure><img src="image_5.png" /><figcaption>**Fig 1.** Maps of (A) Thailand with Nan province highlighted in orange and (B) Nan province elevation and villages (n=905). The red dots represent the villages positions. The base layer of the map is available at https://data.humdata.org/dataset/cod-ab-tha.</figcaption></figure>

https://doi.org/10.1371/journal.pntd.0013552.g001

between 2003 and 2019 was also calculated. Human cases included both confirmed and probable cases as defined by the ICD-10: A75.3 and explained in Wangrangsimakul et al. [3]. The location (longitude, latitude) of each village of the province as well as the population per village were retrieved from national census data. We estimated the elevation in meters of each village and the mean slope in degrees for each village buffer from the Environment Operations Center (www.gms-eoc.org) based on Version 4.1 of NASA's Shuttle Radar Topographic Mission (SRTM) elevation dataset. We used the land use land cover time series maps provided by the European Space Agency Climate Change Initiative (ESA CCI) to extract land cover data [24]. These global land cover maps were built at a spatial resolution of 300m on an annual basis between 1992–2020. We downloaded the 2003–2019 raster maps to match the temporal range of scrub typhus cases. The maps initially encompass 37 global or regional land cover classes numbered from 0 to 220 shown in Supporting Information S1 Fig. We selected and aggregated the following nine land cover classes: Forest broadleaf (corresponding to the n°50, 60, 61 and 62 in the ESA CCI map products classification); Forest needleleaf (n°70, 71, 72, 80, 81 and 82); Cropland rainfed (n°10); Cropland irrigated (n°20); Shrubland (n°120, 121, 122, 152 and 180); Grassland (n°153 and 120); Mosaic (n°30, 40, 100, 110, 150 and 151); Urban areas (n°190) and Water bodies (n°210). For each village and for each year, the area of the nine new land cover classes was extracted in the 2500 m buffer using R software [25] and the sf package v1.0-15 [26] and terra packages v1.7-55 [27]. By comparing raster maps from 2003 and 2019, we characterized 25 different land cover transition patterns that occurred between 2003 and 2019. These new variables allowed us to separate land cover that had remained stable during this period from newly formed land cover. We estimated the area of each land cover transition pattern. All variables are shown in Table 1.

## Statistical analysis

Descriptive analyses were conducted using R software [25] to explore the land cover and spatial distribution and evolution of scrub typhus. A matrix of change was used to examine the land cover change between 2003 and 2019. We identified the villages that overcame land cover changes within the 2500 m buffer zone and the villages that did not. A Mann-Whitney-Wilcoxon non-parametric test was performed to compare the difference in cumulative scrub typhus cumulated cases between these two groups.
<table><caption>Table 1. Description of variables and their use in the models. The area is defined as the 2500 m buffer zone.</caption><thead><tr><th>Variable category</th><th>Variable</th><th>Description</th><th>Used in models</th></tr></thead><tbody><tr><td rowspan="2">Dependent variables</td><td>Scrub_cases</td><td>Number of scrub typhus cases per year and per village</td><td>GAM* land cover (1)**</td></tr><tr><td>Scrub_tot</td><td>Number of cumulated scrub typhus cases per village</td><td>GAM land cover transition (2)***</td></tr><tr><td rowspan="3">Spatio-temporal variables</td><td>Year</td><td>Year of the variable record</td><td>GAM land cover (1)</td></tr><tr><td>Lon</td><td>Longitude of the village</td><td>GAM land cover (1) and GAM land cover transition (2)</td></tr><tr><td>Lat</td><td>Latitude of the village</td><td>GAM land cover (1) and GAM land cover transition (2)</td></tr><tr><td rowspan="2">Topographic variables</td><td>Elevation</td><td>Elevation of the village in m</td><td>GAM land cover (1)</td></tr><tr><td>Slope_mean</td><td>Mean slope of the buffer zone of each village</td><td>GAM land cover (1)</td></tr><tr><td>Demographic variables</td><td>N_men</td><td>Number of men per village</td><td>GAM land cover (1) and GAM land cover transition (2)</td></tr><tr><td rowspan="9">Land cover variables</td><td>Rainfed Cropland cover</td><td>Yearly area of rainfed cropland cover</td><td rowspan="9">GAM land cover (1)</td></tr><tr><td>Grassland cover</td><td>Yearly area of grassland cover</td></tr><tr><td>Mosaic cover</td><td>Yearly area of mosaic cover</td></tr><tr><td>Irrigated cropland cover</td><td>Yearly area of irrigated cropland cover</td></tr><tr><td>Broadleaf forest cover</td><td>Yearly area of broadleaf forest cover</td></tr><tr><td>Urban cover</td><td>Yearly area of urban cover</td></tr><tr><td>Shrubland cover</td><td>Yearly area of shrubland cover</td></tr><tr><td>Needleleaf forest cover</td><td>Yearly area of needleleaf forest cover</td></tr><tr><td>Water cover</td><td>Yearly area of water cover</td></tr><tr><td rowspan="20">Land cover transition variables</td><td>Irrigated cropland to Mosaic</td><td>Area of irrigated cropland that became mosaic</td><td rowspan="20">GAM land cover transition (2)</td></tr><tr><td>Stable rainfed cropland</td><td>Area of rainfed cropland that has been stable</td></tr><tr><td>Rainfed cropland to urban areas</td><td>Area of rainfed cropland that became urban areas</td></tr><tr><td>Rainfed cropland to shrubland</td><td>Area of rainfed cropland that became shrubland</td></tr><tr><td>Rainfed cropland to mosaic</td><td>Area of rainfed cropland that became mosaic</td></tr><tr><td>Rainfed cropland to broadleaf forest</td><td>Area of rainfed cropland that became broadleaf forest</td></tr><tr><td>Grassland to shrubland</td><td>Area of grassland that became shrubland</td></tr><tr><td>Grassland to mosaic</td><td>Area of grassland that became mosaic</td></tr><tr><td>Grassland to broadleaf forest</td><td>Area of grassland that became broadleaf forest</td></tr><tr><td>Stable shrubland</td><td>Area of shrubland that has been stable</td></tr><tr><td>Shrubland to irrigated cropland</td><td>Area of shrubland that became irrigated cropland</td></tr><tr><td>Shrubland to rainfed cropland</td><td>Area of shrubland that became rainfed cropland</td></tr><tr><td>Shrubland to grassland</td><td>Area of shrubland that became grassland</td></tr><tr><td>Shrubland to mosaic</td><td>Area of shrubland that became mosaic</td></tr><tr><td>Shrubland to broadleaf forest</td><td>Area of shrubland that became broadleaf forest</td></tr><tr><td>Shrubland to needleleaf forest</td><td>Area of shrubland that became needleleaf forest</td></tr><tr><td>Stable mosaic</td><td>Area of mosaic that has been stable</td></tr><tr><td>Mosaic to shrubland</td><td>Area of mosaic that became shrubland</td></tr><tr><td>Mosaic to broadleaf forest</td><td>Area of mosaic that became broadleaf forest</td></tr><tr><td>Mosaic to needleleaf forest</td><td>Area of mosaic that became needleleaf forest</td></tr><tr><td>Stable broadleaf forest</td><td>Area of broadleaf forest that has been stable</td></tr><tr><td>Broadleaf forest to shrubland</td><td>Area of broadleaf forest that became shrubland</td></tr><tr><td>Broadleaf forest to mosaic</td><td>Area of broadleaf forest that became mosaic</td></tr><tr><td>Stable needleleaf forest</td><td>Area of needleleaf forest that has been stable</td></tr><tr><td>Needleleaf forest to mosaic</td><td>Area of needleleaf forest that became mosaic</td></tr></tbody></table>

The study and transition periods are between 2003 and 2019.

*GAM = Generalized Additive Model.

** Referring to model number (1) described in the statistical analysis section.

***Referring to model number (2) described in the statistical analysis section.

https://doi.org/10.1371/journal.pntd.0013552.t001
We assessed the spatial autocorrelation of scrub typhus cases using semi-variograms. We performed spatial interpolation of cases using Gaussian process regression (kriging) with the gstat package v 2.1-1 [28] and sp package v1.6-0 [29].

Based on these preliminary analyses, we performed two General Additive Models (GAM) using the mgcv package v1.9-1 [30]. The first considered the spatiotemporal land cover dynamics to explain scrub cases per village and per year. This model included the nine land cover area classes, the number of males per village and topographic variables (Table 1). The structure of the tested GAM is given by the following formula:

$$scrub_{cases} = \beta_0 + f_{year} (Year) + f_{n\_men} (n_{men}) + f_{LonLat} (lon, lat) + f_1 (x_1) + \dots + f_{11} (x_{11}) \quad (1)$$

Where β0 is the intercept, f the unrestricted splines and x1 to x11 were the nine land cover plus the two topographic variables to select, described in Table 1.

We used the number of males per village in the model because scrub typhus affects mainly males than females in Thailand for the period 2001–2018 [22,31]. The sex ratios between our dataset and Thailand are comparable for this period of time (census of 2010) with respectively 98.9 males to 100 females compared to 96.2 males to 100 females [32]. Additional models with the whole population are available in Supporting information S2 Fig. We set a negative binomial function link based on the assessment of the distribution of scrub typhus cases using fitdistrplus package v1.1-11 [33]. All land cover and topographic variables were scaled. We adopted a parsimonious approach to select the final model, based on the AIC indices and we checked model for validity. The basic model considered the relationship between scrub typhus, village population (i.e., number of males) to control for the demographic effect, year to control for the temporal autocorrelation as well as special meteorological events, and longitude and latitude to control for spatial autocorrelation revealed in the previous steps of the analysis. As latitude and longitude interaction within a GAM may introduce spatial confounding, we implemented the spatial+ method described in [34] to account for this bias. Models presented in this article are issued from this method. We tested the explanatory factors one by one, selecting the factors that caused the AIC index to drop the most (See Supporting information S3 Fig. for the results of this selection). Concurvity between variables was then checked, variables were excluded when superior to a threshold of 0.8. We used mgcv package v1.9-1 [30] and DHARMA package v0.4.6 [35] to test the validity of the models and the concurvity. Finally, we checked the quality of the spatiotemporal autocorrelation using ACF plot, the semi-variogram and the maps of residuals.

The second GAM model aimed to explain the cumulative cases of scrub typhus between 2003 and 2019 by the global land cover change in the 2500 m buffer zone that occurred during the same period. The explicative variables included the 25 types of land cover change area, the male population and the interaction between longitude and latitude to account for spatial autocorrelation. The structure of this second GAM was given as follow:

$$scrub_{tot} = \beta_0 + f_{n_{men}} (n_{men}) + f_{LonLat} (lon, lat) + f_1 (x_1) + \dots + f_{25} (x_{25}) \quad (2)$$

Where β0 is the intercept, f the unrestricted splines and x1 to x25 were the 25 land cover transitions and topographic variables to select, described in Table 1.

The family was defined as negative binomial, and we also used a parsimonious approach to select the best model based on the AIC indices. Concurvity was also checked following the same method. Once the GAM model was defined, we replaced the unidentified splines with a polynomial function ranging from 0 to 4 degrees to assess the mathematical equation of the effects.

# Results

## Burden of scrub typhus in Nan province

Between 2003 and 2019, a total of 8,045 cases of scrub typhus were referenced by the passive health surveillance system. The burden of scrub typhus was low in 2003 with 119 cases and increased until 2019 (499 cases), reaching two peaks in
<figure><img src="image_6.png" /></figure>

2012 (with a number of scrub typhus cases of 906) and 2016 (with a number of scrub typhus cases of 721) (Fig 2). The number of cases per village and per year followed a negative binomial distribution and were spatially heterogeneously distributed. As shown in Fig 3A and 3B, the burden of scrub typhus infection between 2003 and 2019 was rather low in the area of Nan River valley, with most cases located in the mountainous and foothill areas of the Northwest, North and East of the valley.

## Land cover within the 2,500 m buffer zone

Nan Province has overcome several land cover changes during the period 2003–2019 (see Supporting information S4 Fig). The province is a rural mountainous area mainly composed of mosaic land (median = 4,203 km²), broadleaf forest (median = 3,249 km²), cropland rainfed (median = 2,394 km²) and shrubland (median = 1,998 km²) between 2003 and 2019. Mosaic land was also predominant within the 2,500 m buffer zones followed by cropland rainfed, broadleaf forest and shrubland.

The Fig 4 depicts the land cover changes. The area of rainfed cropland increased from 2003 to 2008 and reached a plateau. Irrigated cropland decreased overall, with an acceleration from 2010. Grassland increased overall and reached a temporary stable area between 2009 and 2017. Broadleaf and needleleaf forests started to increase in 2009 and 2010, respectively. Mosaic land increased until 2015 and started to decrease slightly thereafter. Finally, shrubland decreased sharply until 2019. Water body areas did not change. Finally, urban areas increased and reached a stable area in 2015.

Spatiotemporal effect of land cover on scrub typhus: GAM model results.

<figure><img src="image_7.png" /><figcaption>Fig 2. Description of scrub typhus human cases: trend of scrub typhus cases per year.</figcaption></figure>

https://doi.org/10.1371/journal.pntd.0013552.g002
A

<figure><img src="image_8.png" /></figure>

B

Mean Scrub
typhus cases
-1 to 0
0 to 1
1 to 2
2 to 3
3 to 4
4 to 5
5 to 6
6 to 7
7 to 8
8 to 9
9 to 10
10 to 11
11 to 12
12 to 13

0 100 150 200 250 300 350 400 km

<caption>Fig 3. Maps of Nan Province: (A) spatial distribution of cumulative scrub typhus cases between 2003 and 2019 and (B) interpolation of mean scrub typhus cases at village level between 2003 and 2019 by kriging using a semi-variogram based on the village position. The base layer map for administrative boundaries is available at https://data.humdata.org/dataset/cod-ab-tha.</caption>

https://doi.org/10.1371/journal.pntd.0013552.g003

A first GAM model (1) was developed considering the yearly scrub typhus cases per village according to the nine yearly land cover surface areas within the 2,500 m buffer zones. The best model was selected according to the Akaike Information Criteria AIC. Intermediate models and their respective AICs are presented in the process of selection described in the Supporting information S3 Fig. The best model explained 64.3% of the deviance (pseudo-R² = 0.54). Rainfed cropland and elevation were excluded from the model due to concurvity problem. The land covers that were significantly associated with the number of scrub typhus cases per village were mosaic land, broadleaf forest (deciduous forest), needleleaf forest and shrubland (Fig 5 and Table 2). The number of scrub typhus cases increased with shrubland and mosaic land areas in the village buffer zones. An inverted U-shape relationship was observed between scrub typhus cases and broadleaf forest area and a negative relationship was observed with needleleaf forest. Mapping of semi-variogram and ACF of the residuals (see Supporting information S5 Fig) showed that the GAM model controlled well for spatial and temporal autocorrelation of scrub typhus cases. The diagnostic of the GAM model (1) is given in Supporting information S6 Fig.

## Land cover change in Nan province between 2003 and 2019

The land cover change matrix (Fig 6) showed that Nan Province was in a process of forest encroachment between 2003 and 2019. By 2019, shrubland cover had decreased by 25% of its initial in 2003. This area had been colonized by mosaic land first, then by rainfed cropland followed by broadleaf forest.

The study of the land cover transition model helps to understand the contribution of land cover to explain scrub typhus cases. Here, we considered only the land cover transition of interest for scrub typhus, identified in the first GAM model (1) as: rainfed crop-land (at first selected but then excluded because of concurvity problem), mosaic land, broadleaf forest (deciduous forest), needleleaf forest and shrubland. Among the 905 villages, 580 of them experienced significant land cover changes between 2003 and 2019.

## Impact of specific land cover transition on scrub typhus

Villages were categorized as whether they had experienced or not a land cover transition for the five relevant land cover classes. The non-parametric Mann-Whitney-Wilcoxon test showed that villages that overcame relevant land cover changes had more scrub typhus cases than stable villages (median of 2 compared to 1, p-value <0.05).
<figure><img src="image_9.png" /><figcaption>**Fig 4.** Change in land cover surface area between 2003 and 2019 per land cover class. The y-axes represent the number of cells, which is a square area of 300 meters by 300 meters.</figcaption></figure>

https://doi.org/10.1371/journal.pntd.0013552.g004

## Association between cumulative scrub typhus cases and detailed trajectory of land cover change: GAM model results

A second GAM model (2) was developed considering the cumulative scrub typhus cases between 2003 and 2019 per village according to the 25 land cover transitions previously identified in the buffer zones between 2003 and 2019. The best model selected according to AIC is shown in the Table 3 and Fig 7, using a negative binomial linked function (theta = 2.74), with an explained deviance of 81.6% (pseudo-R² of 0.71). The selected variables were stable broadleaf forest, shrubland to grassland change, shrubland to broadleaf forest change, as well as the number of males by village and the longitude and latitude interaction term. Stable rainfed cropland and the transition from rainfed cropland to urban areas were at first selected but then excluded because of concurvity problem. The model (2) diagnostic is given in the Supporting information S7 Fig. The Table 4 shows the fitted polynomial equations for each variable. The GAM output with polynomial functions and its diagnostic are available in the Supporting Information S8 Fig.

As shown in Fig 7, stable broadleaf forest area was not a significant risk factor. However, the reforestation process behind the transition from shrubland to broadleaf forest showed an inverted U-shape relationship with scrub typhus human cases. The relationship shows an initial increase at low levels of reforested area, followed by a decrease at higher reforested surface values, likely driven by a single data point with substantial reforestation and greater uncertainty. A polynomial function of degree 3 was found to give the best fit. Loss of shrubland to grassland was not significant (coefficient = -0.018, p-value = 0.28197 in the fitted GAM-polynomial model). Male population was positively correlated with the cumulative scrub typhus cases (polynomial function of degree 2).
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /><figcaption>**Fig 5.** General Additive Modeling (GAM) (1) results of the selected land covers explaining the number of scrub typhus cases in Nan Province between 2003–2019 by villages and year, using a binomial negative link function (theta=1.075). The smoothed variables selected in the best GAM were: (A) year, (B) number of males per village, (C) shrubland cover, (D) mosaic land cover, (E) broadleaf forest cover, (F) needleleaf forest cover and (G) geographical distribution of villages (longitude and latitude).</figcaption></figure>

https://doi.org/10.1371/journal.pntd.0013552.g005

Discussion

Summary of main results

We investigated several environmental factors associated with scrub typhus cases at the village level. First, we confirmed the effect of topography, with relatively few cases in the Nan River valley and most cases in villages localized at higher altitude in mountainous and foothill areas (although slope was not selected in the best model). Our results support our first hypothesis that forest cover was related to scrub typhus cases, but with slight differences. Shrubland and mosaic land were the most important land cover types associated with scrub typhus. Interestingly, needleleaf forest appeared to reduce the risk. Broadleaf forest cover showed an inverted U-shape relationship with scrub typhus cases, suggesting that a very small surface or a large area of forest cover surrounding a village was associated with a lower risk of scrub
<table><caption>Table 2. Results of general additive modelling (GAM) (1) explaining the number of cases of scrub typhus per village in Nan Province using a negative binomial link (theta=1.075), with approximate significance of smooth terms.</caption><thead><tr><th>Explanatory variables</th><th>Estimated degrees of freedom</th><th>Chi square</th><th>P-value</th></tr></thead><tbody><tr><td>S(Year)</td><td>8.557</td><td>708.32</td><td>&lt;2e-16 ***</td></tr><tr><td>S(Mosaic)</td><td>3.016</td><td>22.50</td><td>1.61e-04***</td></tr><tr><td>S(ForestB)</td><td>4.310</td><td>34.59</td><td>3.04e-06***</td></tr><tr><td>S(Shrubland)</td><td>7.362</td><td>60.04</td><td>&lt;2e-16***</td></tr><tr><td>S(ForestN)</td><td>4.762</td><td>49.38</td><td>&lt;2e-16***</td></tr><tr><td>S(Lon, Lat)</td><td>28.177</td><td>2515.86</td><td>&lt;2e-16***</td></tr><tr><td>S(N_men)</td><td>6.671</td><td>5963.26</td><td>&lt;2e-16***</td></tr></tbody></table>

For the best selected model, the deviance explained = 64.3%, R² = 0.54, restricted maximum likelihood (REML) = 8771.8, AIC = 17366.

https://doi.org/10.1371/journal.pntd.0013552.t002

<figure><img src="image_18.png" /><figcaption>Fig 6. Land cover change matrix of Nan Province between 2003 and 2019. "Wat"=Water area, "Urb"=Urban area, "Cropl"=Irrigated cropland, "CropR"=Rainfed cropland, "Grass"=Grassland, "Shrub"=Shrubland, "Mos"=Mosaic, "ForB"=Broadleaf forest, "ForN"=Needleleaf Forest. The coloured bar represents the proportion of land cover of 2019 compared to the 2003 period. For example, Water surface area had not changed and still represents 100% of its surface area of 2003. A 100% of the urbanised areas in 2003 stayed urbanised in 2019, however urban areas also expanded on the cropland irrigated field.</figcaption></figure>

https://doi.org/10.1371/journal.pntd.0013552.g006

typhus transmission. By studying land cover dynamics and transitions, we showed that broadleaf forest cover, which had remained stable over the period, was no longer an important factor in scrub typhus transmission. Conversely, the reforestation process identified in the transition from shrubland to forest appeared significantly at risk for scrub typhus.
<table><caption>Table 3. Results of general additive modelling (GAM) (2) of land cover transition explaining the number of cumulative scrub typhus cases per village in Nan Province using a negative binomial link (theta=2.74), with approximate significance of smooth terms.</caption><thead><tr><th>Explanatory variables</th><th>Estimated degrees of freedom</th><th>Chi square</th><th>P-value</th></tr></thead><tbody><tr><td>S(N_men)</td><td>2.501</td><td>189.027</td><td>&lt;2e-16***</td></tr><tr><td>S(ForestB_stable)</td><td>2.487</td><td>4.251</td><td>0.236</td></tr><tr><td>S(Shrubland_Grassland)</td><td>1.002</td><td>1.619</td><td>0.204</td></tr><tr><td>S(Shrubland_ForestB)</td><td>2.620</td><td>15.287</td><td>0.002**</td></tr><tr><td>S(Lon,Lat)</td><td>26.997</td><td>1589.591</td><td>&lt;2e-16***</td></tr></tbody></table>

For the best selected model, the deviance explained = 81.6%, R² = 0.71, restricted maximum likelihood (REML) = 1408.9, AIC = 2738.

https://doi.org/10.1371/journal.pntd.0013552.t003

<figure><img src="image_19.png" /><figcaption>Fig 7. Results of General Additive Modeling (GAM) (2) of selected land cover transition explaining the total number of scrub typhus cases in Nan Province between 2003–2018 at the village level, using a binomial negative link function (theta=2.74). The smoothed variables selected in the best GAM were: (A) transition from shrubland to grassland, (B) transition from shrubland to broadleaf forest, (C) stable broadleaf forest, (D) the number of males per village and (E) geographical distribution of villages (longitude and latitude).</figcaption></figure>

https://doi.org/10.1371/journal.pntd.0013552.g007

Villages that experienced land cover changes during the period 2003–2019 also recorded significantly higher scrub typhus cases than stable villages. However, these stable villages were not evenly distributed in Nan Province and were located in agrarian valley rather than in foothill and mountainous areas. These results support our second hypothesis that unstable land covers such as secondary growth vegetation or fragmented landscape may increase the burden of scrub typhus
<table><caption>Table 4. Results of general additive modelling (GAM) of land cover transition explaining the number of cumulative cases of scrub typhus per village in Nan Province using a negative binomial link (theta=2.712), with polynomial adapted terms.</caption><thead><tr><th>Explanatory variables</th><th>Estimate</th><th>Std Error</th><th>z-value</th><th>P-value</th></tr></thead><tbody><tr><td>(Intercept)</td><td>-0.062</td><td>0.153</td><td>-0.402</td><td>0.6877</td></tr><tr><td>poly(N_men, 2, raw=T)1</td><td>0.005</td><td>6.578e-04</td><td>7.434</td><td>1.05e-13***</td></tr><tr><td>poly(N_men, 2, raw=T)2</td><td>-1.766e-06</td><td>6.840e-07</td><td>-2.581</td><td>0.0098**</td></tr><tr><td>poly(ForestB_stable, 3, raw=T)1</td><td>0.012</td><td>7.061e-03</td><td>1.765</td><td>0.07756</td></tr><tr><td>poly(ForestB_stable, 3, raw=T)2</td><td>-1.732e-04</td><td>1.283e-04</td><td>-1.350</td><td>0.17687</td></tr><tr><td>poly(ForestB_stable, 3, raw=T)3</td><td>5.054e-07</td><td>6.203e-07</td><td>0.815</td><td>0.41522</td></tr><tr><td>Shrubland_Grassland</td><td>-0.018</td><td>0.016</td><td>-1.076</td><td>0.28197</td></tr><tr><td>poly(Shrubland_ForestB, 3, raw=T)1</td><td>0.080</td><td>0.032</td><td>2.460</td><td>0.01391*</td></tr><tr><td>poly(Shrubland_ForestB, 3, raw=T)2</td><td>-2.753e-03</td><td>2.395e-03</td><td>-1.150</td><td>0.25033</td></tr><tr><td>poly(Shrubland_ForestB, 3, raw=T)3</td><td>1.859e-05</td><td>3.276e-05</td><td>0.567</td><td>0.57038</td></tr></tbody></table>

For the best selected model, the deviance explained = 81.5%, R² = 0.662, restricted maximum likelihood (REML) = 1466.1, AIC = 2736.766.

https://doi.org/10.1371/journal.pntd.0013552.t004

compared to stable land covers. The spatial interaction term in the models could have introduced spatial confounding [34,36]. We checked for such bias and found that the corrected models presented in the article were highly similar to the initial models (see Supporting Information S9 Fig) which strengthen the explanatory power of the land cover and land cover transition variables.

The scrub typhus vector mite, may not tolerate all types of land cover change. Gain in forest cover, as illustrated by the transition from shrubland to forest, appears to increase the risk of scrub typhus transmission when the surface area is intermediate (i.e., inverted U-shape relationship). Shrubland to grassland change studied in the transition model (2) suggests that less bushy coverage might reduce the transmission of the disease (yet not significant). This may indicate that the vector of scrub typhus requires some tree cover to thrive and may not tolerate too much habitat degradation. The ecological states and the functional processes of forest ecosystems are important for scrub typhus dynamics. Further research on vector ecology is needed to understand how.

The importance of forest cover and vegetation highlighted in this study had been previously observed in several studies in China and Taiwan [15–17,37–39]. Our results are also in accordance with Wardrop et al. [17] and Li et al. [15] who identified transitional land covers and in particular mosaic land and shrubland as significant variables explaining scrub typhus cases. The impact of land cover dynamics on scrub typhus had been poorly studied [39], but the study of Min et al. [18] suggested that scrub typhus thrives with secondary growth vegetation. Min et al. [18] observed this phenomenon in a deforestation context in South Korea based on Global Forest Watch data. Deforestation was not significant in our modelling results although we observed that stable forest was less at risk than newly grown forest from 2003 to 2019. Wangrangsimakul et al. [3] investigated scrub typhus incidence in Chiang Rai Province in Northern Thailand, near Nan Province. Consistent with our study, they also observed a positive association between elevation and scrub typhus. These authors also used land cover data from the European Spatial Agency, though at 1 km resolution. However, in contrast to our results, they observed that forest cover and mosaic land decreased scrub typhus incidence while habitat complexity increased it. These differences could be explained by the differences in scales (subdistrict level versus villages level) and the presence of different variables (such as landscape complexity).

To explain the importance of secondary growth vegetation, further studies should investigate the vectors (chiggers), hosts and pathogen in relation to the environment. Kuo et al. [13] compared rodent infestation by chiggers and ticks in fallow (abandoned fields) and ploughed fields in Taiwan. The burden of chiggers was three times higher in fallows where
<figure><img src="image_20.png" /></figure>

secondary vegetation provided suitable microhabitats. An early study revealed that in Malaysia in 1974, the distribution of the vector mirrored the distribution of its main host, *Rattus tanemuzi* [40]. This rodent species is generalist and synanthropic and should be favoured in disturbed areas. The importance of this rodent host for the chiggers and the pathogen has been confirmed in more recent studies in Southeast Asia [41] and more specifically in northern Thailand [42]. Network analysis showed that *Rattus tanezumi* with the other synanthropic rodent *Bandicota indica* as well as the chigger species *Leptotrombidium deliense* and *Walchia kritochaeta* emerged as central nodes in rodents-chiggers networks [42]. *Orientia tsutsugamushi* positive chiggers and rodents were more abundant in lowland than in forested areas. However, another study in Thailand found that chigger species richness was higher in forested areas rather than in human-disturbed habitats [6]. A similar study was also conducted in Thailand on the association between human land use and occurrence of *Orientia tsutsugamushi* in rodents [43]. Rodent infected with *Orientia tsutsugamushi* were also more likely to be found in forested habitats. The authors suggested that rodents were likely infected in habitats such as houses, fallow land or rice fields when they were close to forested areas. The results on chiggers and rodents led the authors to suggest that a minimum threshold of biodiversity was required for vectors, hosts and the pathogen. However, these authors also showed that scrub typhus incidence was negatively associated with host-parasite network connectance, suggesting that high complexity of interactions might reduce human exposure, or by giving non-vector species a better chance to dominate. This may be reflected in the inverted U-shaped relationship between broadleaf forest cover and scrub typhus cases. An alternative explanation is that disturbed areas such rainfed cropland, may not provide the micro-climate requirements for the vectors [14]. Another local study in Malaysia highlighted the importance of habitat complexity and ecotone as a potential risk factor [44]. In this study, the authors found that chigger species richness on rodents was highest at forest edges compared to four other land cover classes. However, none of the rodents were infected with *Orientia tsutsugamushi* in this habitat, while the highest prevalence was observed in rodents trapped in oil palm plantation. To further support the importance of forest edge and transitional land covers, a previous empirical study in a temperate area (Nebraska) indicated that chiggers were found all along forest edges. But the vectors particularly thrived in short to tall-grass transition zones rather than in the understorey with important tree canopy [45]. Understanding scrub typhus ecology in term of habitat remains very challenging. Further empirical data on vector distribution, combined with detailed habitat description, are needed.

## Limitations

Our studies had several limitations. First, the type of satellite images available showed mixed results. The literature on land use change has shown that Nan Province suffers from deforestation, mainly due to the expansion of maize and commercial tree plantations [46]. These results were not observed using ESA CCI land cover data but using a combination of different satellite images with much higher resolution. One explanation is that what is recognised as mosaic land and shrubland cover in ESA CCI land covers was often considered forest cover in some studies [20,46]. Another explanation is that the rubber plantation cover, which doubled in area between 1990 and 2019 [46], was recognised as forest or mosaic land and not as agricultural land. This leads to the second limitation, commercial plantations such as rubber, teak or orchard are not classified, making it difficult to validate or invalidate the observation of a positive association between scrub typhus cases and commercial plantations [47]. The same problem applies to some crops such as maize, which has increased dramatically between 2003 and 2019. A third limitation is associated with agricultural practices, which make land covers unstable, such as slash-and-burn cultivation, crop rotation with fallow that reduces over the years [48]. Therefore, it remains difficult to study land cover changes in an area with such low resolution, either in area or classification. A fourth limitation concerns the impact of abiotic factors. Abiotic factors such as rainfall and temperature have a significant impact on scrub typhus transmission [49]. Our results showed (Fig 2) that scrub typhus cases increased significantly in 2010 and 2016 during ENSO events [50]. The year variable in GAM model (1) accounted for such climatic variation during the period 2003–2019. However, no data of good quality was available at such a fine resolution. Nan province possesses only
two meteorological stations and the prediction of temperature and humidity at the village level did not enabled us to use this data with good satisfactory along the land cover variables. A fifth limitation is the lack of information on the behaviour of populations and their exposure to risky habitats around their villages. Socio-economic and vocation data that could explain such information were not available at the village level which prevented us to use it in the model. Finally, diagnostic limitations and passive surveillance likely led to underreporting. Most hospitals lack confirmatory tests, in reliance on clinical criteria with moderate accuracy. Incomplete reporting and variable care-seeking behaviour further underestimate the true burden.

Choosing the immediate vicinity of the villages as the area of interest allowed us to link local ecological studies on vectors, hosts, and pathogens with regional epidemiological analysis of scrub typhus burden in humans. As sampling the vector when it is free-living in the environment is a high effort with low efficiency, this study represents a further step towards a comprehensive understanding of scrub typhus risk habitats and transmission ecology. In particular, it contributes to the knowledge of the scrub typhus situation in northern Thailand, an endemic area. Understanding the relationship between land cover and human cases offers new opportunities for developing strategies to reduce the human burden of disease. To summarize, our study highlighted tree, shrub, second growth vegetation and unstable land cover as risky habitats where needleleaf forests and broadleaf (deciduous) forests, if stable, may have a protective effect. A higher resolution of the land cover map could refine the conclusions, especially regarding agricultural land and process of forest fragmentation we could not assess here.

## Supporting information

**S1 Fig. Land use–land cover classification for Nan province.** Extracted form Land cover CCI product user guide version 2.0.

(DOCX)

**S2 Fig. Charts of general additive model (GAM (1)) explaining the number of scrub typhus cases according to land cover surface area, elevation, slope and the number of people (male and female) per villages and GAM (2) explaining the cumulated number of scrub typhus cases according to land cover change surface area, elevation and the number of people (male and female) per villages.**

(PDF)

**S3 Fig. Selection of Gam model based on AIC.**

(PDF)

**S4 Fig. Maps of Nan Province:** (A) land cover in 2003, (B) land cover in 2019, and (C) land cover change between 2003 and 2019 (*LC=land cover). The base layer of the map of administrative boundaries is available at https://data.humdata.org/dataset/cod-ab-tha.

(PDF)

**S5 Fig.** (A) Map of the GAM residuals, (B) ACF plot verifying temporal autocorrelation of GAM land cover residuals, (C) semivariogram of GAM land cover residuals verifying spatial autocorrelation. The base layer of the map is available at https://data.humdata.org/dataset/cod-ab-tha.

(PNG)

**S6 Fig. Land cover GAM diagnostic.**

(PDF)

**S7 Fig. Land cover transition GAM diagnostic.**

(PDF)
**S8 Fig.** Charts of general additive model (GAM) explaining the number of cumulated cases of scrub typhus per villages in Nan province using a negative binomial link with polynomial adapted terms and its diagnostic. Adapted GAM-polynomial model.

(PDF)

**S9 Fig.** Land cover and land cover transition GAMs and their diagnostic without spatial confounding correction (GAM (1): pages 1&2, GAM (2): pages 3&4).

(PDF)

**Acknowledgments**

We would like to thank Areeya Kriengudom for her significant help in getting the most updated version of Nan province data. We also thank Dr. Theirno Madiou Bah for sharing his precious knowledge on GAM and model selection methodology. We also thank Dr. Rodolphe Loubaton for his important help on the understanding and implementation of the spatial+ method.

**Author contributions**

**Conceptualization:** Nolwenn Blache, Karine Chalvet-Monfray, Serge Morand.

**Data curation:** Nolwenn Blache.

**Formal analysis:** Nolwenn Blache, Karine Chalvet-Monfray.

**Investigation:** Nolwenn Blache, Karine Chalvet-Monfray, Rawadee Kumlert, Soawapak Hinjoy, Serge Morand.

**Methodology:** Nolwenn Blache, Karine Chalvet-Monfray, Serge Morand.

**Resources:** Rawadee Kumlert, Soawapak Hinjoy, Serge Morand.

**Software:** Nolwenn Blache, Karine Chalvet-Monfray, Serge Morand.

**Supervision:** Karine Chalvet-Monfray, Serge Morand.

**Validation:** Nolwenn Blache, Karine Chalvet-Monfray, Rawadee Kumlert, Soawapak Hinjoy, Serge Morand.

**Visualization:** Nolwenn Blache.

**Writing – original draft:** Nolwenn Blache.

**Writing – review & editing:** Nolwenn Blache, Karine Chalvet-Monfray, Rawadee Kumlert, Soawapak Hinjoy, Serge Morand.

**References**

1. Taylor AJ, Paris DH, Newton PN. A Systematic Review of Mortality from Untreated Scrub Typhus (Orientia tsutsugamushi). PLoS Negl Trop Dis. 2015;9(8):e0003971. https://doi.org/10.1371/journal.pntd.0003971 PMID: 26274584

2. Watt G, Parola P. Scrub typhus and tropical rickettsioses. Curr Opin Infect Dis. 2003;16(5):429–36. https://doi.org/10.1097/00001432-200310000-00009 PMID: 14501995

3. Wangrangsimakul T, Elliott I, Nedsuwan S, Kumlert R, Hinjoy S, Chaisiri K, et al. The estimated burden of scrub typhus in Thailand from national surveillance data (2003-2018). PLoS Negl Trop Dis. 2020;14(4):e0008233. https://doi.org/10.1371/journal.pntd.0008233 PMID: 32287307

4. Suputtamongkol Y, Suttinont C, Niwatayakul K, Hoontrakul S, Limpaiboon R, Chierakul W, et al. Epidemiology and clinical aspects of rickettsioses in Thailand. Ann N Y Acad Sci. 2009;1166:172–9. https://doi.org/10.1111/j.1749-6632.2009.04514.x PMID: 19538278

5. Elliott I, Pearson I, Dahal P, Thomas NV, Roberts T, Newton PN. Scrub typhus ecology: a systematic review of Orientia in vectors and hosts. Parasit Vectors. 2019;12(1):513. https://doi.org/10.1186/s13071-019-3751-x PMID: 31685019

6. Chaisiri K, Gill AC, Stekolnikov AA, Hinjoy S, McGarry JW, Darby AC, et al. Ecological and microbiological diversity of chigger mites, including vectors of scrub typhus, on small mammals across stratified habitats in Thailand. Anim Microbiome. 2019;1(1):18. https://doi.org/10.1186/s42523-019-0019-x PMID: 33499969
7. Wei C-Y, Wang J-K, Shih H-C, Wang H-C, Kuo C-C. Invasive plants facilitated by socioeconomic change harbor vectors of scrub typhus and spotted fever. PLoS Negl Trop Dis. 2020;14(1):e0007519. https://doi.org/10.1371/journal.pntd.0007519 PMID: 31961870

8. Chaisiri K, Linsuwanon P, Makepeace BL. The chigger microbiome: big questions in a tiny world. Trends Parasitol. 2023;39(8):696-707. https://doi.org/10.1016/j.pt.2023.05.002 PMID: 37270375

9. Traub R, Wisseman CL Jr, Jones MR, O'Keefe JJ. The acquisition of Rickettsia tsutsugamushi by chiggers (trombiculid mites) during the feeding process. Ann N Y Acad Sci. 1975;266:91-114. https://doi.org/10.1111/j.1749-6632.1975.tb35091.x PMID: 829479

10. Audy JR, Harrison JL. A review of investigations on mite typhus in Burma and Malaya, 1945-1950. Trans R Soc Trop Med Hyg. 1951;44(4):371-404. https://doi.org/10.1016/s0035-9203(51)80018-x PMID: 14817817

11. Audy JR. A summary topographical account of scrub typhus, 1908-1946. Kuala Lumpur: Printed at the Govt. Press. 1949.

12. Dohany AL, Shirai A, Lim BL, Huxsoll DL. Variation in populations of chigger vectors of scrub typhus in developing oil palm areas of different ages. Jpn J Med Sci Biol. 1980;33(5):263-70. https://doi.org/10.7883/yoken1952.33.263 PMID: 7300038

13. Kuo C-C, Huang J-L, Shu P-Y, Lee P-L, Kelt DA, Wang H-C. Cascading effect of economic globalization on human risks of scrub typhus and tick-borne rickettsial diseases. Ecol Appl. 2012;22(6):1803-16. https://doi.org/10.1890/12-0031.1 PMID: 23092017

14. Misumi H, Takahashi M, Urakami H, Matsumoto I. Distributions of infective spots composed of unfed larvae infected with Orientia tsutsugamushi in Leptotrombidium mites and their annual fluctuations on the soil surface in an endemic area of tsutsugamushi disease (Acari: Trombiculidae). Med Entomol Zool. 2002;53(4):227-47. https://doi.org/10.7601/mez.53.227

15. Li X, Wei X, Yin W, Soares Magalhaes RJ, Xu Y, Wen L, et al. Using ecological niche modeling to predict the potential distribution of scrub typhus in Fujian Province, China. Parasit Vectors. 2023;16(1):44. https://doi.org/10.1186/s13071-023-05668-6 PMID: 36721181

16. Wei X, He J, Yin W, Soares Magalhaes RJ, Wang Y, Xu Y, et al. Spatiotemporal dynamics and environmental determinants of scrub typhus in Anhui Province, China, 2010-2020. Sci Rep. 2023;13(1):2131. https://doi.org/10.1038/s41598-023-29373-7 PMID: 36747027

17. Wardrop NA, Kuo C-C, Wang H-C, Clements ACA, Lee P-F, Atkinson PM. Bayesian spatial modelling and the significance of agricultural land use to scrub typhus infection in Taiwan. Geospat Health. 2013;8(1):229-39. https://doi.org/10.4081/gh.2013.69 PMID: 24258898

18. Min K-D, Lee J-Y, So Y, Cho S-I. Deforestation Increases the Risk of Scrub Typhus in Korea. Int J Environ Res Public Health. 2019;16(9):1518. https://doi.org/10.3390/ijerph16091518 PMID: 31035715

19. Registration Office Department of the Interior, Ministry of the Interior. รายงานสถิติจำนวนประชากรและบ้านประจำปี Statistics, population and house statistics for the year 2019. 2019.

20. Trisurat Y, Shirakawa H, Johnston JM. Land-Use/Land-Cover Change from Socio-Economic Drivers and Their Impact on Biodiversity in Nan Province, Thailand. Sustainability. 2019;11(3):649. https://doi.org/10.3390/su11030649 PMID: 33354352

21. Nan Provincial Statistical Office, National Statistical Office. รายงานสถิติจังหวัดน่าน Nan Provincial Statistical Report 2018. n.d.

22. Bureau of Epidemiology. Annual Epidemiology Surveillance Report. Thailand: Ministry of Public Health; 2001.

23. สำนักระดับวิทยา :: Bureau of Epidemiology, Thailand. n.d. [cited January 16, 2025]. Available from: https://apps-doe.moph.go.th/boe/diseases.php?ds_key=dGJz&dsid=MTE=&ds=U2NydWlgdHlwaHVz

24. ESA. Land Cover CCI Product User Guide Version 2 2017.

25. R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing. 2022.

26. Pebesma E, Bivand R. Spatial Data Science: With Applications in R. CRC Press. 2023.

27. Hijmans RJ, Bivand R, Cordano E, Dyba K, Pebesma E, Sumner MD. terra: Spatial Data Analysis. 2024.

28. Pebesma EJ. Multivariable geostatistics in S: the gstat package. Computers & Geosciences. 2004;30(7):683-91. https://doi.org/10.1016/j.cageo.2004.03.012

29. Pebesma EJ, Bivand RS. Classes and methods for spatial data in R. R News. 2005;5:9-13.

30. Wood SN. Fast Stable Restricted Maximum Likelihood and Marginal Likelihood Estimation of Semiparametric Generalized Linear Models. Journal of the Royal Statistical Society Series B: Statistical Methodology. 2010;73(1):3-36. https://doi.org/10.1111/j.1467-9868.2010.00749.x

31. Gonwong S, Mason CJ, Chuenchitra T, Khanijou P, Islam D, Ruamsap N, et al. Nationwide Seroprevalence of Scrub Typhus, Typhus, and Spotted Fever in Young Thai Men. Am J Trop Med Hyg 2022;106:1363-9. https://doi.org/10.4269/ajtmh.20-1512

32. National Statistical Office of Thailand. Executive summary. 2010.

33. Delignette-Muller ML, Dutang C. fitdistrplus: AnRPackage for Fitting Distributions. J Stat Soft. 2015;64(4). https://doi.org/10.18637/jss.v064.i04

34. Dupont E, Augustin NH. Spatial Confounding and Spatial+ for Nonlinear Covariate Effects. JABES. 2023;29(3):455-70. https://doi.org/10.1007/s13253-023-00586-7

35. Hartig F. DHARMA: Residual Diagnostics for Hierarchical (Multi-Level/ Mixed) Regression Models. 2022.

36. Clayton DG, Bernardinelli L, Montomoli C. Spatial correlation in ecological analysis. Int J Epidemiol. 1993;22(6):1193-202. https://doi.org/10.1093/ije/22.6.1193 PMID: 8144305

37. Huang X, Xie B, Long J, Chen H, Zhang H, Fan L, et al. Prediction of risk factors for scrub typhus from 2006 to 2019 based on random forest model in Guangzhou, China. Trop Med Int Health. 2023;28(7):551-61. https://doi.org/10.1111/tmi.13896 PMID: 37230481
38. Liu L, Xiao Y, Wei X, Li X, Duan C, Jia X, et al. Spatiotemporal epidemiology and risk factors of scrub typhus in Hainan Province, China, 2011-2020. One Health. 2023;17:100645. https://doi.org/10.1016/j.onehlt.2023.100645 PMID: 38024283

39. Blache N, Chalvet-Monfray K, Déprés C, Morand S. A scoping review of the impacts of forest dynamics on acari-borne diseases: beyond forest fragmentation 2024. https://doi.org/10.21203/rs.3.rs-4185885/v1

40. Traub R, Wisseman CL Jr. The ecology of chigger-borne rickettsiosis (scrub typhus). J Med Entomol. 1974;11(3):237-303. https://doi.org/10.1093/jmedent/11.3.237 PMID: 4212400

41. Morand S, Bordes F, Blasdell K, Pilosof S, Cornu J, Chaisiri K, et al. Assessing the distribution of disease-bearing rodents in human-modified tropical landscapes. Journal of Applied Ecology. 2015;52(3):784-94. https://doi.org/10.1111/1365-2664.12414

42. Elliott I, Thangnimitchok N, Chaisiri K, Wangrangsimakul T, Jaiboon P, Day NPJ, et al. Orientia tsutsugamushi dynamics in vectors and hosts: ecology and risk factors for foci of scrub typhus transmission in northern Thailand. Parasit Vectors. 2021;14(1):540. https://doi.org/10.1186/s13071-021-05042-4 PMID: 34663445

43. Chaisiri K, Cosson J-F, Morand S. Infection of Rodents by Orientia tsutsugamushi, the Agent of Scrub Typhus in Relation to Land Use in Thailand. Trop Med Infect Dis. 2017;2(4):53. https://doi.org/10.3390/tropicalmed2040053 PMID: 30270910

44. Alkathiry H, Al-Rofaai A, Ya'cob Z, Cutmore TS, Mohd-Azami SNI, Husin NA, et al. Habitat and Season Drive Chigger Mite Diversity and Abundance on Small Mammals in Peninsular Malaysia. Pathogens. 2022;11(10):1087. https://doi.org/10.3390/pathogens11101087 PMID: 36297144

45. Clopton RE, Gold RE. Distribution and seasonal and diurnal activity patterns of Eutrombicula alfreddugesi (Acari: Trombiculidae) in a forest edge ecosystem. J Med Entomol. 1993;30(1):47-53. https://doi.org/10.1093/jmedent/30.1.47 PMID: 8433345

46. Kruasilp J, Pattanakiat S, Phutthai T, Vardhanabindu P, Nakmuenwai P. Evaluation of Land Use Land Cover Changes in Nan Province, Thailand, Using Multi-Sensor Satellite Data and Google Earth Engine. Environ Nat Resour J. 2023;21(2):1-12. https://doi.org/10.32526/ennrj/21/202200200

47. Shah HA, Huxley P, Elmes J, Murray KA. Agricultural land-uses consistently exacerbate infectious disease risks in Southeast Asia. Nat Commun. 2019;10(1):4299. https://doi.org/10.1038/s41467-019-12333-z PMID: 31541099

48. Turkelboom F, Poesen J, Trébuil G. The multiple land degradation effects caused by land-use intensification in tropical steeplands: A catchment study from northern Thailand. CATENA. 2008;75(1):102-16. https://doi.org/10.1016/j.catena.2008.04.012

49. Wang Q, Ma T, Ding F-Y, Lim A, Takaya S, Saraswati K, et al. A systematic review of environmental covariates and methods for spatial or temporal scrub typhus distribution prediction. Environ Res. 2024;263(Pt 1):120067. https://doi.org/10.1016/j.envres.2024.120067 PMID: 39341542

50. Chueasa B, Humphries UW, Waqas M. Influence of El Niño southern oscillation on precipitation variability in Northeast Thailand. MethodsX. 2024;13:102954. https://doi.org/10.1016/j.mex.2024.102954 PMID: 39315397