Association between Landscape Factors and Spatial Patterns of Plasmodium knowlesi Infections in Sabah, Malaysia

Forest loss and other environmental changes correlate with increased malaria incidence.

studies have implicated the primarily exophagic Anopheles leucosphyrus group of mosquitoes as the main vector of P. knowlesi and have found relatively high biting rates in farm edges bordering forests and forest areas (11)(12)(13)(14). Human-disturbed environments have been associated with changes in behavior of nonhuman primates and their increased contact with humans (15). Fragmentation of existing habitats can also increase the frequency of disease transmission by creating transition areas with increased spatial overlap among human, mosquito, and wildlife populations or by altering vector ecology (16,17). The effects of these changes at forest edges have been described for malaria and other vector-borne zoonotic diseases (18,19) but not for P. knowlesi. A previous mathematic modeling study highlighted the potential for increased transmission resulting from increased spatial overlap among people, macaques, and mosquitoes at forest edges (20).
Despite these apparent links between land use and P. knowlesi transmission, detailed environmental risk factors for P. knowlesi infections in humans are unknown. Although variability of P. knowlesi risk has been reported at a regional scale, patterns of P. knowlesi transmission have not been described at a subdistrict spatial scale (2,21). Furthermore, studies that formally evaluate associations between P. knowlesi and characteristics of environment and landscape are lacking.
After obtaining approval from the Medical Research and Ethics Committee of the Ministry of Health in Malaysia, we examined the changing incidence of P. knowlesi in Kudat and Kota Marudu districts in northwestern Sabah, Malaysia, on the island of Borneo, an area with relatively high P. knowlesi transmission (22). Our aim was to describe the spatial and temporal patterns of P. knowlesi incidence within these districts and to explore potential associations between village-level incidence and deforestation and other environmental factors. Clarifying these relationships is vital to predicting and responding to future disease outbreaks and understanding the underlying mechanisms of P. knowlesi emergence.

Study Site and Population
This study was conducted in the districts of Kudat and Kota Marudu in northwestern Sabah, Malaysia (07.38°-06.19°N, 116.62°-117.46°E), an area of 3,204 km 2 with a population of ≈120,000 persons predominantly of Rungus and Dusun ethnicities (23) (Figure 1). The climate is tropical, with no dry season and increased rainfall during November-March; the area has both coastal and inland regions and elevations ranging from sea level to 1,000 m above sea level. Substantial environmental change is ongoing in the region because of conversion of land for oil palm plantations and other agricultural activities (24 and treatment free of charge. All clinics refer patients to the central district hospital, where hospitalization is mandatory for malaria patients until a negative blood smear for malaria parasites has been obtained.

Geolocation of Patients and Villages
We conducted a retrospective review of malaria patients reported by Kudat and Kota Marudu district hospitals during 2008-2012. Data on diagnosis, admission date, demographics, and address of all malaria patients were obtained from hospital laboratory microscopy records. Villages and populations were identified from the 2010 Population and Housing Census in Malaysia (23), the most recent census conducted in this area, and the global positioning system coordinates of village centroids were recorded as part of a larger interdisciplinary study (MONKEYBAR ESEI project; http://malaria.lshtm.ac.uk/MONKEYBAR). Village populations were updated by using published population growth rates for Malaysia (25). All locations were confirmed by using imagery available through Google Earth (https://www.google.com/earth/) or other freely available satellite data. Patient addresses were matched to the census data or to the nearest reported village from the census data. Administrative boundaries were used to define the extent of urban areas, within which village data were combined.
Calculating the Proportion of Malaria Patients with P. knowlesi P. knowlesi is microscopically similar to P. malariae but can also be misdiagnosed as P. falciparum or other human malaria species. In hospital microscopy records, the species of malaria is recorded, as determined by morphology, but no separate listing for P. knowlesi exists. Consequently, uncertainty in diagnosis resulted in some infections being recorded as P. malariae or P. malariae/knowlesi. To estimate the true proportion of P. knowlesi, we calculated the sensitivity and specificity of microscopy diagnosis of P. malariae for a subset of 539 malaria patients for whom both microscopy and molecularly confirmed results were available, including all patients from Kudat and Kota Marudu hospitals who were referred to a tertiary care hospital during this period and all patients recruited at these hospitals during 2013 and 2014 (22,26). The proportion of malaria cases reported as P. knowlesi per village per year was adjusted for this sensitivity and specificity by using a Bayesian estimation of true incidence from apparent incidence obtained by testing individual samples (27). This estimation is calculated as follows: Models were fitted in R (http://www.R-project.org) through the prevalence and rjags packages interfacing with JAGS version 3.4.0 (http://mcmc-jags.sourceforge.net/) by using 2 chains containing 1,000 burn-in samples and 5,000 retained samples. Sensitivity and specificity parameters were determined by using β-PERT distributions of the minimum, maximum, and most likely values.

Environmental Data
Topography and land use data were extracted from various datasets derived from satellite-based remote-sensing data. These data were evaluated for buffer areas (i.e., areas within a certain radius of a village center) with a radius of 1, 2 and 5 km from the center of each village; these distances were chosen to explore a range of spatial scales at which environment may be relevant on the basis of the typical distribution of households, farming land, and local human and animal movements. Elevation data with a spatial resolution of 30 m were obtained from the ASTER Digital Global Elevation Model (28). The average annual normalized difference vegetation index (NDVI), which quantifies the greenness of vegetation, was calculated from Moderate Resolution Imaging Spectroradiometer 16-day composites at 250-m resolution (29). The NDVI is influenced by climatic factors (e.g., rainfall and temperature) and has been used extensively to predict malaria incidence and develop early warning systems in other contexts (30,31).
Tree cover data, derived from classified Landsat imagery at 30-m resolution, were obtained from Hansen et al. (32). Annual forest cover maps for the study districts were produced; forest was defined as >50% tree crown cover density. Although this land classification represents forested areas, it cannot distinguish types of forest or agroforestry such as rubber or oil palm. The proportion of forest coverage, proportion of forest loss during the year for which incidence was estimated, and proportion of cumulative forest loss for the previous 5 years (i.e., total forest loss for 2006-2010 was evaluated by using 2010 incidence) were calculated for each buffer radius for each village and time point. We used the Landscape Ecology Statistics plugin for Quantum GIS (33) to evaluate the effect of forest configuration as the number of forest patches per radius, a standard metric representing landscape fragmentation. Distributions of these variables were examined, and quartiles were used to categorize variables.
Because this analysis relied on passive reporting of malaria, we included travel time to the nearest clinic where patients would seek treatment for a febrile illness as a measure of access to care. Travel time to the clinic from each village was estimated by using travel times reported in community interviews and by patients recruited as part of a population-based case-control study (26).
Annual P. knowlesi incidence for each village was mapped and smoothed incidence maps produced to visualize the data by using a kernel density estimation method, a standard method for interpolating point location data. Because the data were skewed relative to Poisson distribution, potential associations between environmental factors and reported P. knowlesi patients at the village level were assessed by using general linearized mixed models with a negative binomial distribution and an offset for population in R (34). To account for correlation between repeat measurements for the same village, we included the village variable as a random effect. Bivariable analysis was conducted for each covariate; variables for which p<0.2 were included in multivariable models. We used likelihood ratio tests to assess the significance of single variables and the Akaike Information Criterion (http://www. modelselection.org/aic/) for the final model selection. For correlated variables (e.g., mean elevation at different buffer radiuses), a single variable was selected for inclusion on the basis of marginal Akaike Information Criterion values. Potential bias from residual spatial autocorrelation in the model was explored for island and mainland areas through Moran's I. On the basis of this statistic, the negative binomial model was fit with a spatial correlation component estimated by using a distance-based Matern correlation function.
Most villages (60%, 245/405) reported at least 1 malaria patient during this period. The number of malaria patients reported varied annually, with marked seasonal variations in numbers of patients and amount of rainfall. Almost half (47%, 878/1,847) of reported malaria patients had suspected P. knowlesi or P. malariae infections diagnosed by microscopy. P. falciparum and P. vivax malaria were diagnosed in another 27% (512/1,847) and 25% (457/1,847) of patients, respectively.

Association with Environmental Variables
On the basis of estimates from remote sensing data (32), forest cover declined by 4.8% in Kudat and Kota Marudu during 2008-2012 ( Figure 2, panel B). Although this loss was often highly localized, large tracts of forest were cleared in the interior of Pulau Banggi. Overall, substantial environmental change was observed in the districts; 39% (157/405) of villages lost >10% of forest cover within a 1-km radius during this 5-year period; 44% (179/405) lost >10% within a 2-km radius; and 51% (206/405) lost >10% within a 5-km radius.
Lower elevations, >65% of forest cover (within a 2-km radius), and higher historical forest loss were associated with greater incidence of P. knowlesi infections. Travel times to clinics, average annual NDVI, and number of forest patches were not significant at any radius in the bivariate analysis.

Spatial Patterns
Maps of interpolated P. knowlesi incidence show distinctive spatial patterns that vary from year to year, with some areas of high incidence persisting over time (Figure 2, panel A). To assess the degree of residual spatial correlation after fitting the negative binomial regression, the residuals of the best-fitting model were mapped and Moran's I was calculated (mainland Borneo, Moran's I 0.07, p<0.0001; islands off Borneo's coast, Moran's I 0.07, p = 0.07). Because Moran's I was significant, the final model was adjusted for spatial autocorrelation (Tables 2, 3).

Discussion
This study aimed to describe spatial and temporal patterns of P. knowlesi infection in northern Sabah, Malaysia, and to evaluate potential associations between village-level P. knowlesi incidence and key environmental factors. Although land use changes can affect emergence of infectious diseases, such associations have not been previously evaluated for P. knowlesi. We found links among deforestation, environmental characteristics, and reported incidence of P. knowlesi malaria in this region. Village-level malaria data show marked spatial and temporal heterogeneity in P. knowlesi incidence in Kudat and Kota Marudu districts. After adjustment for the sensitivity and specificity of microscopy, P. knowlesi was the most common cause of human malaria, a finding consistent with other studies for this geographic area (4,5,22). Quantification of annual forest loss indicates that substantial environmental changes occurred during 2008-2012, with many villages losing substantial proportions of surrounding forest cover.
The proportion of forest surrounding a village was associated with incidence of P. knowlesi infections in the final model (Table 3); this association potentially reflects the role of forest environments as habitats of macaques and mosquito vectors. Studies report higher vectorial capacity and sporozoite rates for P. knowlesi in forest environments than in agricultural and settled lands (12,35). Long-tailed macaques have also been reported in various environments, including degraded secondary forest areas (36). The association of increased P. knowlesi incidence with both forest and forest loss likely confirms findings that transmission is occurring in forested areas undergoing substantial change, as with previously described frontier malaria (37).
Higher incidence of P. knowlesi was associated with higher proportions of forest loss surrounding villages during the 5-year period before year of reported incidence. This association could result from changes in macaque or mosquito habitats and from increased levels of human activity. Increased density of long-tailed macaques has been reported as a response to deforestation; loss of previous habitats can result in crowding within forest patches, with potential implications for disease transmission (36). Land use changes have also been shown to affect abundance and community composition of potential vectors (17,38). Deforestation and associated agricultural development are also associated with changes in human risk because of altered distribution and behavior of humans; these changes result from more employment opportunities and shifts in human movement patterns because of forest clearing and agricultural activities (19,37). Although clearing of forests may initially deplete vector populations and thereby reduce malaria transmission, this reduction may be followed by colonization of cleared areas by more efficient vector species and subsequent increases in transmission (9). Historical forest loss was more significantly associated with P. knowlesi incidence than forest loss occurring during the same year of reported incidence (Table 2), suggesting that increased transmission is related to long-term changes in vector, host, or human populations involved. Additional longitudinal studies are required to investigate this hypothesis. The effect of habitat may also be reflected in associations with elevation. Elevation was negatively correlated with P. knowlesi incidence, although the range of elevations within the study site was limited; most villages were at elevations <100 m. Both macaques and vectors are more frequently described in low elevation areas but have been reported at higher elevations (36,39).
Neither forest configuration nor vegetation indices appeared to be strongly associated with P. knowlesi incidence. Because NDVI is a measure of vegetation greenness, this measure may not differentiate between forest and other types of agriculture (e.g., oil palm) because NDVI tends to become saturated in tropical environments. Studies of other zoonotic diseases have found that the fragmentation and configuration of other types of land cover, in addition to forest, influence disease transmission (18,19). Other fragmentation metrics could be included to evaluate the effects of forest patch size and shape on P. knowlesi transmission.
The buffer sizes that we evaluated represent a range of potential scales at which variables (e.g., mean elevation) within different radiuses may affect P. knowlesi transmission. Human infections result from many factors interacting across different spatial scales, and the strength of association of these factors likely varies by distance. Our study covers an extensive area in which behaviors, distribution of villages, and ecology vary. Villages are typically within small spatial areas, but farming practices range from small-scale swidden farming to large-scale plantations and affect the scale of human interactions with the environment. Future analyses would benefit from including more detailed spatial data on household locations and human movement patterns.
The main limitation of this study is the reliance on records of malaria patients who sought care at hospitals; these patients may represent only a portion of malaria patients in the community. Although malaria is a notifiable disease and distance to a hospital was included as a measure of access to care, asymptomatic malaria cases and symptomatic cases that resolved without treatment are unaccounted for. Previous studies of other malaria species in similar transmission settings have described a large proportion of asymptomatic carriage within communities; however, asymptomatic carriage has not been evaluated for P. knowlesi (40). A cross-sectional survey of 2,019 persons in central Vietnam identified 3 persons whose samples were positive for P. knowlesi, yet all 3 were asymptomatic at the time of the survey and for the subsequent 6   This study was also limited by the environmental data used. Data on forest cover and forest loss were aggregated by year, and finer-scale temporal associations between land cover (i.e., forest and other land types) and incidence could not be explored. This dataset was limited by defining forest as canopy cover, a definition that does not enable differentiation between types of forest or crops or between patches of different types of land. In addition, analysis was limited by the spatial resolution of satellitebased remote-sensing data, the use of a centroid point to represent village location, and the use of circular buffers rather than buffers of actual village shape. Exploratory spatial analysis suggested spatial heterogeneity in villagelevel data, and these spatial effects were included in the model. Additional work with more spatially specific outcomes and environmental data are needed to investigate these spatial patterns in more detail.
Despite inherent limitations in the outcome and covariate data used in this study, results strongly suggest a link between environmental change and reported incidence of emerging P. knowlesi in northern Sabah. Spatial analysis of environmental factors affecting disease emergence can be used to target surveillance and public health activities to areas expected to have increasing disease risk. Although additional population-based studies are needed to define environmental risk factors, this study indicates that deforestation is associated with human cases of P. knowlesi within northwestern Sabah, Malaysia.