Introduction

Malaria continues to inflict a great health and socio-economic burden on humanity, with an estimated 3.2 billion people at risk of being infected1. In 2018, globally there were 228 million cases and 405,000 deaths, around 67% (272,000) of deaths were in children aged under 5 years2. However, in 2018, there were 23 million fewer cases as compared to 20102. In 2016, malaria remained endemic in 91 countries and territories as compared to 108 in 20003. The World Health Organization (WHO) African Region accounts for around 90% of malaria cases globally, followed by the South-East Asian Region (SEAR) (5%) and the Eastern Mediterranean Region (2%)4. Some of the factors that have led to the observed reductions in malaria incidence since 2000, are intensification of malaria control interventions supported by unprecedented financial support, socio-economic improvement in endemic countries and increasing urbanization5,6,7,8. In 2018, total investment for malaria control and elimination was US$ 2.7 billion2.

WHO developed the Global Technical Strategy for Malaria 2016–2030 (GTS)5 with an aim to fast track progress towards malaria elimination. This strategy is complemented by the Roll Back Malaria advocacy plan, Action and Investment to Defeat Malaria 2016–2030 (AIM)9. GTS and AIM set a global goal to eliminate malaria in at least 21 countries by 2020, known as E-2020 countries and 35 countries by 20303,9.

Malaria is reported from seven districts of Bhutan along the southern border with India. These districts are Chukha, Dagana, Pemagatshel, Samdrup Jongkhar, Samtse, Sarpang and Zhemgang (Fig. 1). Malaria control activities in Bhutan are based on: (1) Early diagnosis and prompt treatment with artemisinin-based combination therapy (ACT), (2) Protection of at-risk populations with long-lasting insecticide nets (LLINs) and indoor residual spraying (IRS), and (3) integrated vector management (IVM). With the dwindling of malaria cases, Bhutan was aiming to eliminate malaria by 201810. However, there are still malaria cases and an updated elimination goal date has been set to 2020.

Figure 1
figure 1

Map of Bhutan with malaria transmitting districts.

Full size image

Small-scale geographic variation in transmission becomes increasingly important in the elimination phase as malaria cases are confined to hard-to-reach populations, often in border areas, that could act as source of reintroduction into areas which have achieved elimination11,12. Therefore, the aims of this study were to quantify the spatial and temporal patterns of malaria and associations between malaria risk and climatic factors. The information from this study can be used to intensify control measures to achieve malaria elimination goals in Bhutan.

Results

Descriptive analysis

A total of 2,062 P. falciparum and 2,284 P. vivax cases were recorded during the study period, of which 328 (8.2%) were mixed infection with both species. There were more cases amongst males older than 5 years than females and male children <5 years. In 2006, there were 1,721 cases and numbers subsequently reduced over the study period, with fewer than 20 cases in 2013 and 2014 respectively. There were no cases in the <5 years age group in 2013 and 2014 (Table 1). Maps of crude standardized morbidity ratios (SMR) of P. falciparum and P. vivax revealed higher risk in Samdrup Jongkhar and Sarpang districts, particularly in sub-districts adjacent to the international borders (Fig. 2).

Table 1 Malaria incidence stratified by sex and age group during the study period (2006–2014).
Full size table
Figure 2
figure 2

Raw standardised morbidity ratios of (a) Plasmodium falciparum and (b) Plasmodium vivax by sub-districts in Bhutan, 2006–2014.

Full size image

Time-series decompositions

For both P. falciparum and P. vivax, time-series decompositions of raw data showed a clear seasonal pattern. Large peaks occurred for each parasite species in May of each year, with incidence dropping off during the remainder of the rainy season. The inter-annual pattern showed a large peak in 2006 and a small peak in 2010 with lower incidence in the intervening and subsequent years for both species of malaria (Fig. 3).

Figure 3
figure 3

Decomposed Plasmodium falciparum and Plasmodium vivax of Bhutan, 2006–2014.

Full size image

Spatial autocorrelation analysis

Hot spots of both P. falciparum and P. vivax were reported in central and eastern parts of the Bhutan-India border region, while high-high clustering was located in the central part of the border region (Figs. 4 and 5). Eight sub-districts in Sarpang were in high-high clusters for P. falciparum and P. vivax. Fifteen sub-districts were in low-low clusters in three districts of Chukha, Dagana and Pemagatshel for both P. falciparum and P. vivax (Supplementary Fig. 1).

Figure 4
figure 4

Time series hot spot analysis of Plasmodium falciparum of Bhutan, 2006–2014.

Full size image
Figure 5
figure 5

Time series hot spot analysis of Plasmodium vivax of Bhutan, 2006–2014.

Full size image

Spatio-temporal model

Table 2 describes the three spatio-temporal models: Model I, containing the unstructured random effects and Model III containing both the unstructured and structured random effects, had similar fit (indicated by a difference of deviation information criterion [DIC] of <3), and were better fitting than Model II containing the structured random effect only, for both P. falciparum and P. vivax. Here we present the results for Model I which is the simpler model. For P. falciparum, there was an estimated increase of 0.7% (95% credible intervals [CrI] 0.3%, 0.9%) for a one mm increase in rainfall. P. falciparum cases decreased by 2% (95% CrI 0.4%, 4.0%) for each 1 °C increase of maximum temperature, but this was not statistically significant. Similarly, a one mm increase in rainfall and 1 °C increase of maximum temperature was associated with a 0.2% increase (95% CrI −0.1%, 0.6%) and 0.1% increase (95% CrI −0.3%, 1.0%), respectively, in P. vivax incidence, but were statistically not significant. Sex and age were not statistically significant predictors of either species of malaria. The variance of the random effects for P. falciparum and P. vivax in this model were 0.16 (95% CrI 0.10, 0.24) and 0.20 (95%CrI 0.13, 0.29), respectively (Table 2). The means of the spatially unstructured random effects when mapped showed evidence of spatial clustering after accounting for the model covariates, despite these random effects not being spatially structured in the model (Fig. 6). Some sub-districts had a high probability of being above (or below) the overall mean residual risk for the study period (Supplementary Fig. 2). In general, these sub-districts matched the high-high and low-low clusters identified in the exploratory Local Indicators of Spatial Association (LISA) analysis.

Table 2 Regression coefficients and 95% CrI from Bayesian spatial and non-spatial models of Plasmodium falciparum and P. vivax cases reported by month and sub-districts, Bhutan, 2006–2014.
Full size table
Figure 6
figure 6

Spatial distribution of the posterior means of unstructured random effects for (a) Plasmodium falciparum and (b) Plasmodium vivax in Model I.

Full size image

Discussion

This study aimed to describe the spatial and temporal epidemiology of malaria in the pre-elimination setting of Bhutan, which is aiming to eliminate malaria by 202013. Time series decomposition showed seasonal peaks for both P. falciparum and P. vivax. Rainfall was associated with an increased risk of P. falciparum. Malaria clusters and hot spots of both species were located in Sarpang district.

Malaria transmission in Bhutan seems to have been interrupted and the risk of malaria is equal in both genders and across age groups (typical of low-transmission settings where there is little acquired immunity14). This finding is in line with an earlier cross sectional survey which showed zero prevalence15. Ongoing cases might relate to small localised clusters of transmission created by movement of infected individuals across the international border. Transmission interruption could be attributed to a number of factors including good coverage and use of LLINs15,16, prompt diagnosis and treatment and enhanced surveillance17.

Malaria transmission hot spots for both species continued to occur in the same region over 10 years18. Some of the reasons outlined for these hot spots were their proximity to a very porous border16,19 with the Indian state of Assam, which reported some of the highest malaria incidence rates in India20,21,22,23,24. However, there has been a significant reduction in the number of cases in Assam since 2016 as India has accelerated towards elimination by 203025.

The seasonal peak of malaria corresponds well with the monsoon in Bhutan, which begins in June and lasts until September. Temperature plays a crucial role in the transmission cycle of the malaria parasite and mosquito survival26,27. Studies found that at a temperature of 22 °C, the life cycle of malaria parasite development in mosquito vector is completed in less than 3 weeks28. The biting rate and gonotrophic processes are also temperature dependent29,30. Other studies have reported rainfall as an important driver of malaria transmission31,32. However, the only climate association found in the current study was between rainfall and P. falciparum risk, perhaps reflecting the disruption to pre-intervention transmission dynamics caused by high LLIN and insecticide coverage in the study area.

Plasmodium vivax risk was not associated with any covariates in the model. A plausible reason could be that, unlike P. falciparum, P. vivax infection can hide in the liver, lying dormant (in hypnozoite form) and protected from the external environment until a relapse is activated. P. vivax represents the most frequent cause of malaria outside of Africa33,34. RDTs remain the cornerstone of diagnosis of P. vivax in many countries. However, RDTs are not adequately sensitive to diagnose the hypnozoite in the liver or in pregnant women. This makes it a particularly challenging parasite to identify and eliminate.

P. vivax is treated with cholorquine (three doses) and radical cure is done using primaquine over two-week period34. A recent study showed that choloroquine and primaquine are still effective against P. vivax in Bhutan34. Therefore, ascertaining adherence to the radical cure is crucial as Bhutan attempts to eliminate malaria by 2020. Not adhering to the treatment continues to be an important factor that leads to continued P. vivax transmission in other parts of the world35,36. It is recommended to undertake operational research to understand barriers and enablers to adherence to two weeks of radical cure of P. vivax in Bhutan to develop strategies to improve adherence.

The main strength of this study was the fine resolution of spatial analysis at the sub-district level over a long time series (108 months). However, there are some limitations that are worth noting. One of the major limitations in using surveillance data is that the completeness and representativeness of such data cannot be ascertained. Secondly, populations of sub-districts were projected and may have led to over or under estimation. Thirdly, there was no reconciliation to accommodate the different level of the climate variables (district) to the disease data (sub-district), and the climate conditions were assumed to be homogeneous within a district. Fourthly, all malaria cases were confirmed using microscopy. Even though conventional light microscopy is considered a “gold standard” for malaria parasite identification and confirmation, light microscopy has a detection limit of 5–50 parasites/μL37,38. Therefore, patients with low parasite intensity (<5 parasites/μL) could have been missed. Fifthly, P. vivax cases could have been inflated due to relapse because this was not accounted for when diagnosing P. vivax. Lastly, unmeasured risk modifiers, such as socio-economic development, living standards, treatment, localised behavioural patterns, population mobility, and bed net use and residual indoor insecticide coverage were unaccounted for in this study.

Conclusion

Plasmodium falciparum transmission in Bhutan was associated with rainfall. Hot spots of both P. falciparum and P. vivax were isolated in sub-districts of Sarpang district. A high residual risk area of malaria transmission was identified in the same sub-districts. Targeted distribution of resources, including intensified interventions in this part of the sub-district will be required for local malaria elimination. Cross-border surveillance needs to be strengthened because of the risk of cross-border malaria to elimination and maintenance of elimination. Maps could be used to target surveillance.

Methods

Study site and data

Bhutan is located in the Eastern Himalayas, bordering China in the north and India in the east, south and west. The country is divided administratively into 20 districts and 205 sub-districts. In 2017, the total population of Bhutan was 681,72039. Historically, malaria transmission occurred in 82 sub-districts of seven districts with 209,090 people living in these districts (Fig. 1).

In this study, malaria data of the seven districts from 2006 to 2014 were obtained from the national malaria surveillance system, hosted by the Bhutan Vector-borne Disease Control Program (VDCP). All other districts were assumed to be non-endemic and were excluded from the analysis. The dataset contained laboratory-confirmed malaria cases, defined as clinically diagnosed cases with either malaria parasite species, confirmed by microscopy. The infections were categorised by species: P. falciparum and P. vivax and stratified by gender and two age groups (<5 and ≥5 years) at the sub-district level.

Population estimates used in this study were from publications from the National Statistical Bureau and the Office of the Census Commissioner of Bhutan40,41. An electronic map of district boundaries in shapefile format was obtained from Global Administrative Areas database (http://www.gadm.org/country).

Exploration of seasonal patterns and inter-annual patterns

The average monthly number of malaria cases was calculated from the full time series (January 2006-December 2014). The time series of malaria incidence was decomposed using seasonal-trend decomposition based on locally (STL) weighted regression to show: the seasonal pattern, inter-annual patterns and the residual variability. The STL model was structured as follows:

$${Y}_{t}={S}_{t}+{T}_{t}+{R}_{t}$$

where Yt represents numbers of local malaria cases with logarithmic transformation, St is the additive seasonal component, Tt is the trend, Rt is the “remainder component” and t is time in months42,43,44,45.

Spatial autocorrelation analysis

At a global (study area) scale, Moran’s I statistic was used to explore spatial autocorrelation and its strength and to test the assumption of spatial independence. Anselin Local Moran’s I statistic and Getis-Ord Gi* statistic were used to undertake local (sub-districts) level clustering and hot spots analysis46. Clustering is the occurrence of unusual aggregation of events in a sub-district. Hot spot is a form of clustering where the target sub-district reports higher rates than the study area average.

Anselin Local Moran’s I (LISA) was used to identify high-risk clusters, low-risk clusters and outliers (high-low and low-high). A positive Local Moran’s I (significant cluster) value shows that the target sub-district with high-risk cluster is surrounded by sub-districts with high-risk cluster or the target sub-district with low-risk cluster is surrounded by sub-districts with low-risk cluster. In a negative Local Moran’s I (outliers) value, the target sub-district with high-risk cluster is surrounded by sub-districts with low-risk or low-risk cluster is surrounded by high-risk cluster sub-districts46. Clusters are determined by comparing Moran’s I values of the target sub-districts and its neighbouring sub-districts to Moran’s I values of all sub-districts in the study area46. In significant high-risk cluster (p  ≤  0.05), the observed Moran’s I value is larger than the expected Moran’s I value. Where as for low-risk cluster, an opposite relationship is observed (observed Moran’s I is smaller than the expected Moran’s I value). In the case of spatial outliers (not significant) Moran’s I value remains in a neutral class46. These classifications of clusters were represented with sub-district boundaries. The Getis-Ord Gi* and Anselin Local Moran’s I cluster analytical methods were used to generate clusters and this served as a sensitivity analysis.

Crude standardized morbidity ratios

Crude standardized morbidity ratios (SMRs) analyses of both species were undertaken to describe the malaria incidence by sub-districts across the study period (9 years). SMR was calculated from:

$${Y}_{i}=\frac{{O}_{i}}{{E}_{i}}$$

Yi-is the overall SMR in sub-district i, Oi - the total number of reported malaria cases in the sub-district i and Ei- expected number of malaria cases in the sub-district i. The Ei was derived by multiplying the average population for sub-district i with the national incidence of malaria12.

Independent variable selection

Best fit climatic covariates were selected using a Poisson regression and Akaike’s information criterion (AIC). Climatic variables (maximum and minimum temperature and rainfall) without a lag, and with one and two-month lag times, plus altitude, were fitted into univariate Poisson regression models. Maximum temperature, unlagged rainfall and altitude showed the best fit with the lowest AIC values. However, maximum and minimum temperature were found to be highly co-linear when tested using variance inflation factors (VIF). Therefore, minimum temperature was dropped from the final model, which included variables rainfall and maximum temperature.

Spatio-temporal model

The number of zero counts for P. falciparum and P. vivax was 34,142 (96.4%) and 34,340 (96.9%). Therefore, we ran an analysis to determine the best model. Zero-inflated Poisson (ZIP) regression was selected over the standard Poisson regression because ZIP had a better fit with lower AIC and BIC as compared to Poisson regression and a Vuong test showed the two models were statically different (Supplementary Tables 1–4). Bayesian statistical software WinBUGS version 1.4 (Medical Research Council, Cambridge, UK and Imperial College London, UK) was used to run ZIP regression models for P. falciparum and P. vivax. Three different models were tested for each species; first model included only climatic variables (rainfall and maximum temperature) as explanatory covariates. Second model contained spatially structured random effects. Third model contained both climatic covariates and spatially structured random effects. The model with the lowest DIC was selected as the final explanatory model for each species.

The most comprehensive model, which had an outcome as the observed counts of malaria, Y, for ith sub-district (i = 1…82) in the jth month (January 2006-December 2014), age group k, and sex group l was structured as follows:

$$P({Y}_{ijkl}={y}_{ijkl})=\{\begin{array}{c}\omega +1(1-\omega ){e}^{-\mu },\,{y}_{ijkl}=0\\ (1-\omega ){e}^{-\mu }\,{\mu }_{ijkl}^{{y}_{ijkl}}/{y}_{ijkl},\,{y}_{ijkl} > 0;\end{array}$$

Yijkl   ̴Poisson (μijkl)

log (μijkl) = log(Eijkl) + θijkl

θijkl = α + β1 × Agek + β2 × Sexl + β3 × rainfallij + β4 × Tempmaxij + ui + si

where E is the expected number of cases (acting as an offset to control for population size) and θ is the mean log relative risk (RR); α is the intercept, and β1, β2, β3 and β4 the coefficients for age (<5 years reference), sex (male reference), rainfall and maximum temperature, respectively; ui is the unstructured random effect (assumed to have a mean of zero and variance σu2) and si is the spatially structured random effect (assumed to have a mean of zero and variance σs2).

Spatially structured random effect was modelled with a conditional autoregressive (CAR) prior structure. An adjacency weights matrix was used to calculate the spatial relationships between the sub-districts, a weight of “1” was assigned when two sub-districts shared a border and “0” if they did not. For the intercept, a flat prior distribution was specified, whereas a normal prior distribution was specified for the coefficients. The priors for the precision of unstructured and spatially structured random effects were specified using non-informative gamma distributions with shape and scale parameters equal to 0.01. Models were also developed without the structured and unstructured random effects to assess whether the inclusion of these components improved model fit.

An initial burn-in of 10,000 iterations was run and these iterations were discarded. Subsequent blocks of 20,000 iterations were run and examined for convergence. Convergence was assessed by visual inspection of posterior density, history plots and Gelman-Rubin statistics, and occurred at approximately 100,000 iterations for each model. Hundred thousand values from the posterior distributions of each model parameter were stored and summarised for the analysis (posterior mean and 95%CrI).

In all analyses, an α-level of 0.05 was adopted to indicate statistical significance (as indicated by 95% CrI for relative risks [RR] that excluded 1). ArcMap 10.5 software (ESRI, Redlands, CA) was used to generate maps of the posterior means of the unstructured and structured random effects and the spatiotemporal random effects obtained from the three models.