Landscape drives zoonotic malaria prevalence in non-human primates

  1. Emilia Johnson  Is a corresponding author
  2. Reuben Sunil Kumar Sharma
  3. Pablo Ruiz Cuenca
  4. Isabel Byrne
  5. Milena Salgado-Lynn
  6. Zarith Suraya Shahar
  7. Lee Col Lin
  8. Norhadila Zulkifli
  9. Nor Dilaila Mohd Saidi
  10. Chris Drakeley
  11. Jason Matthiopoulos
  12. Luca Nelli
  13. Kimberly Fornace
  1. School of Biodiversity, One Health and Veterinary Medicine, University of Glasgow, United Kingdom
  2. Department of Disease Control, London School of Hygiene & Tropical Medicine, United Kingdom
  3. Centre on Climate Change and Planetary Health, London School of Hygiene & Tropical Medicine, United Kingdom
  4. Faculty of Veterinary Medicine, Universiti Putra Malaysia, Malaysia
  5. Lancaster University, Bailrigg, United Kingdom
  6. Liverpool School of Tropical Medicine, Pembroke Place Liverpool, United Kingdom
  7. School of Biosciences, Cardiff University, United Kingdom
  8. Wildlife Health, Genetic and Forensic Laboratory, Sabah Wildlife Department, Wisma Muis, Malaysia
  9. Danau Girang Field Centre, Sabah Wildlife Department, Malaysia
  10. Department of Infection Biology, London School of Hygiene & Tropical Medicine, United Kingdom
  11. Saw Swee Hock School of Public Health, National University of Singapore, Singapore

eLife assessment

This useful study presents findings regarding the impact of forest cover and fragmentation on the prevalence of malaria in non-human primates. The evidence supporting the claims of the authors is solid.

https://doi.org/10.7554/eLife.88616.4.sa0

Abstract

Zoonotic disease dynamics in wildlife hosts are rarely quantified at macroecological scales due to the lack of systematic surveys. Non-human primates (NHPs) host Plasmodium knowlesi, a zoonotic malaria of public health concern and the main barrier to malaria elimination in Southeast Asia. Understanding of regional P. knowlesi infection dynamics in wildlife is limited. Here, we systematically assemble reports of NHP P. knowlesi and investigate geographic determinants of prevalence in reservoir species. Meta-analysis of 6322 NHPs from 148 sites reveals that prevalence is heterogeneous across Southeast Asia, with low overall prevalence and high estimates for Malaysian Borneo. We find that regions exhibiting higher prevalence in NHPs overlap with human infection hotspots. In wildlife and humans, parasite transmission is linked to land conversion and fragmentation. By assembling remote sensing data and fitting statistical models to prevalence at multiple spatial scales, we identify novel relationships between P. knowlesi in NHPs and forest fragmentation. This suggests that higher prevalence may be contingent on habitat complexity, which would begin to explain observed geographic variation in parasite burden. These findings address critical gaps in understanding regional P. knowlesi epidemiology and indicate that prevalence in simian reservoirs may be a key spatial driver of human spillover risk.

eLife digest

Zoonotic diseases are infectious diseases that are transmitted from animals to humans. For example, the malaria-causing parasite Plasmodium knowlesi can be transmitted from monkeys to humans through mosquitos that have previously fed on infected monkeys. In Malaysia, progress towards eliminating malaria is being undermined by the rise of human incidences of ‘monkey malaria’, which has been declared a public health threat by The World Health Organisation.

In humans, cases of monkey malaria are higher in areas of recent deforestation. Changes in habitat may affect how monkeys, insects and humans interact, making it easier for diseases like malaria to pass between them. Deforestation could also change the behaviour of wildlife, which could lead to an increase in infection rates. For example, reduced living space increases contact between monkeys, or it may prevent behaviours that help animals to avoid parasites.

Johnson et al. wanted to investigate how the prevalence of malaria in monkeys varies across Southeast Asia to see whether an increase of Plasmodium knowlesi in primates is linked to changes in the landscape. They merged the results of 23 existing studies, including data from 148 sites and 6322 monkeys to see how environmental factors like deforestation influenced the amount of disease in different places.

Many previous studies have assumed that disease prevalence is high across all macaques, monkey species that are considered pests, and in all places. But Johnson et al. found that disease rates vary widely across different regions. Overall disease rates in monkeys are lower than expected (only 12%), but in regions with less forest or more ‘fragmented’ forest areas, malaria rates are higher. Areas with a high disease rate in monkeys tend to further coincide with infection hotspots for humans. This suggests that deforestation may be driving malaria infection in monkeys, which could be part of the reason for increased human infection rates.

Johnsons et al.’s study has provided an important step towards better understanding the link between deforestation and the levels of monkey malaria in humans living nearby. Their study provides important insights into how we might find ways of managing the landscape better to reduce health risks from wildlife infection.

Introduction

Zoonotic infectious diseases arise from the spillover of pathogens into human populations, typically from a reservoir in wildlife hosts. Anthropogenic land use and land cover change have now been widely linked to infectious disease outbreaks (Brock et al., 2019; Davidson et al., 2019a; Loh et al., 2016). Such practices, including deforestation, logging, clearing for cash-crop plantations or conversion of intact forest into arable land, are accelerating across tropical forests of Southeast Asia (Fornace et al., 2021; Imai et al., 2018; Fornace et al., 2021; Imai et al., 2018). Mechanisms that underly the association between habitat disturbance and spillover risk from wildlife hosts are complex and occur over multiple spatial scales (Brock et al., 2019). In Brazil, re-emergence of Yellow Fever Virus in both NHPs and humans has been linked to areas with highly fragmented forest (Ilacqua et al., 2021). In part, an increase in ‘edge’ habitat in fragmented or mosaic landscapes can facilitate spatial overlap and altered contact patterns between wildlife, vectors, and humans (Lehman et al., 2006). Such ecological interfaces are also thought to contribute to parasite spillover in other vector-borne diseases including Zika (Li et al., 2021a), Babesiosis and Lyme disease (Simon et al., 2014), Trypanosoma cruzi (Vaz et al., 2007) and zoonotic malaria (Brock et al., 2019; Grigg et al., 2017). At the same time, habitat fragmentation can have detrimental impact on wildlife population viability, with reduced host species occupancy and reduced disease burden in highly disturbed habitats (Hanski and Ovaskainen, 2000). Disentangling this interplay is essential to inform ecological strategies for surveillance and mitigation of diseases in regions undergoing landscape change (Fornace et al., 2021).

Zoonotic P. knowlesi is a public health threat of increasing importance across Southeast Asia, following the identification of a prominent infection foci in Borneo in 2004 (Singh et al., 2004). P. knowlesi is a zoonosis, with a sylvatic cycle circulating in non-human primates (NHPs). Human cases currently occur only from spillover events (Ruiz Cuenca et al., 2022; Fornace et al., 2022; Fornace et al., 2023; Lee et al., 2011). Human transmission requires bites from infective mosquitos, primarily anopheline mosquitos of the Leucosphyrus Complex (Anopheles balabacensis, An. latens, An. introlactus) and Dirus Complex (An. dirus, An. cracens) (Moyes et al., 2016; Vythilingam et al., 2006; Wong et al., 2015). Natural hosts for P. knowlesi are typically Long-tailed macaques (Macaca fascicularis) and Southern Pig-tailed macaques (M. nemestrina) (Moyes et al., 2016), both occurring widely across Southeast Asia. Currently, distribution of P. knowlesi cases is thought to be restricted to the predicted ranges of known vector and host species (Davidson et al., 2019b), although recent studies have also identified other NHPs found to be harbouring P. knowlesi. This includes Stump-tailed macaques (M. arctoides), which are now considered to be another natural reservoir (Fungfuang et al., 2020).

Progress towards malaria elimination in Malaysia has been stymied by a recent rise in human incidence of P. knowlesi malaria. Even after accounting for increases in surveillance and diagnostic improvements it is now recognised as the most common cause of clinical malaria in Malaysia (Cooper et al., 2020). Indeed, Malaysia was the first country not to qualify for malaria elimination due to ongoing presence of zoonotic malaria and the WHO updated the guidelines to reflect zoonotic malaria as a public health threat (World Health Organization, 2021). Emergence of Plasmodium knowlesi infections has been linked to changes in land cover and land use (Fornace et al., 2021). While sporadic cases have been reported across Southeast Asia, including in Indonesia (Setiadi et al., 2016), the Philippines (Fornace et al., 2018), Vietnam (Maeno et al., 2015), Brunei (Koh et al., 2019), and Myanmar (Ghinai et al., 2017), the majority of P. knowlesi cases are found in East Malaysia (Borneo) with hotspots in the states of Sabah and Sarawak (Jeyaprakasam et al., 2020), areas that have seen extensive deforestation and landscape modification. In Sabah, human prevalence of P. knowlesi infection has recently been shown to be specifically associated with recent loss of intact forest, agricultural activities, and fragmentation across multiple localised spatial scales (Brock et al., 2019; Fornace et al., 2019b; Fornace et al., 2016).

Prevalence of the pathogen in reservoir hosts is one of three crucial factors determining the force of infection in zoonotic spillover events (Murray and Daszak, 2013). Despite this, very little is known of the impact of rapid landscape change on the distribution of P. knowlesi in NHPs. Literature on the impacts of fragmentation on primates tends to focus on primate density and abundance (Link et al., 2010; Zunino et al., 2007). What is known is that effects of land cover changes on primate-pathogen dynamics are highly variable and context-specific. Although the vector species responsible for sylvatic transmission remain unknown, the Anopheles leucospryphus group, the only vector group implicated in P. knowlesi transmission, is widely associated with secondary, disturbed forest (Brant, 2011; Hawkes et al., 2019; Wong et al., 2015). Macaques have been known to preferentially rely on fringe habitat, a behaviour that may be exaggerated in response to habitat fragmentation and facilitate exposure to vectors (Lehman et al., 2006; Stark et al., 2019). Changes to land composition can also create the biosocial conditions for higher rates of parasitism in primates. Under conditions of limited resources and reduction in viable habitat, conspecific primate density may increase as troops compete for available space. In turn, this can favour transmission via intra-species contact or allow the exchange of pathogens between troops dwelling in interior forest versus edge habitat (Faust et al., 2018; Stark et al., 2019). Habitat use may also become more intensive, preventing parasite avoidance behaviours (Nunn and Dokey, 2006). Land cover change is also known to favour more adaptable, synanthropic species such as M. fascicularis (McFarlane et al., 2012). Considering the spillover risk posed by wildlife reservoirs of P. knowlesi, clarifying any relationships between environmental factors and parasitaemia in key host species may contribute to a more comprehensive understanding of P. knowlesi transmission patterns.

Earth Observation (EO) data provides novel opportunities to investigate epidemiological patterns of diseases which are linked to environmental drivers (Kalluri et al., 2007). In relation to P. knowlesi, utility of fine-scale remote-sensing data has been demonstrated: examples include satellite-derived data used to examine household-level exposure risk in relation to proximate land configuration (Fornace et al., 2019b), UAV-imagery used to link real-time deforestation to macaque host behavioural change (Stark et al., 2019), and remote-sensing data used to interrogate risk factors for vector breeding sites (Byrne et al., 2021). Although macroecological studies that utilise geospatial data are often confounded by issues of matching temporal and spatial scales, as well as by the quality and accuracy of available georeferencing, measures can be taken to account for this when examining the role of environmental factors in modulating disease outcomes. Furthermore, ecological processes occur and interact over a range of distances, or ‘spatial scales’ (Brock et al., 2019; Fornace et al., 2016; Loh et al., 2016). This applies to determinants of vector-borne disease ecology, from larval breeding microclimate to wildlife host foraging behaviour. As multiple influential variables are rarely captured by a single scale (Cohen et al., 2016), data-driven methods can be applied to examine risk factors over multiple scales and identify covariates at their most influential extent (Byrne et al., 2021).

We hypothesise that prevalence of P. knowlesi in primate host species is spatially heterogeneous and that higher prevalence is partially driven by forest loss and fragmentation, contributing to the strong associations described between land use, land cover and human P. knowlesi risk. This study is the first to systematically assess P. knowlesi prevalence in NHPs at a regional scale, and across a wide range of habitats. In conceptual frameworks and transmission models, it is often assumed that P. knowlesi infections in NHPs are chronic (low level, persistent infection) and ubiquitous (uniformly distributed across populations; Brock et al., 2016; Jeyaprakasam et al., 2020). No studies have systematically assessed the extent and quality of all available data on P. knowlesi in NHPs. Independent studies investigating P. knowlesi in primates are typically constrained by small sample sizes and confined geographic areas, limiting inference that can be made about relationships between infection dynamics and landscape characteristics. Systematic tools developed for epidemiological studies of disease prevalence in human populations are rarely applied to the study of wildlife disease prevalence; however, such tools can be used to capture the scale and contrast required in macroecological studies to quantify disease burdens regionally. Furthermore, while recent research has shown the impact of deforestation on the distribution of macaques in the context of P. knowlesi (Moyes et al., 2016; Stark et al., 2019), associations between landscape and variation in the prevalence of simian Plasmodium spp. in primates have not been explored. We aimed to (1) assemble a georeferenced dataset of P. knowlesi in NHPs; (2) evaluate variation in NHP P. knowlesi prevalence by geographic region; and (3) assess environmental and spatial risk factors for P. knowlesi prevalence in NHPs across Southeast Asia.

Results

A systematic literature review was conducted in Medline, Embase, and Web of Science to identify articles reporting prevalence of naturally acquired Plasmodium knowlesi in NHPs. Twenty-three research articles were identified (Akter et al., 2015; Amir et al., 2020; Chang et al., 2011; Fungfuang et al., 2020; Gamalo et al., 2019; Ho et al., 2010; Zamzuri et al., 2020; Jeslyn et al., 2011; Kaewchot et al., 2022; Lee et al., 2011; Li et al., 2021b; Muehlenbein et al., 2015; Nada-Raja et al., 2022; Putaporntip et al., 2010; Saleh Huddin et al., 2019; Salwati et al., 2017; Seethamchai et al., 2008; Shahar, 2019; Gilhooly et al., 2016; personal correspondence, 2013 and 2015 ; Vythilingam et al., 2008; Yusuf et al., 2022; Zhang et al., 2016), containing 148 unique primate survey records to form the dataset for analyses (see Appendix 2 for details of JBI Critical Assessment) (Munn et al., 2015). Year of sampling ranges from 2004 to 2019. No primatological studies were identified from Vietnam, Brunei, or Timor-Leste. Full characteristics of the articles and individual study methodologies are reported in Appendix 1—table 2. Spatial resolution of the survey sites varied from GPS point coordinates to country-level administrative boundaries (Appendix 5—table 1). Geographic distribution of sampling is illustrated in Figure 1.

Sampling sites and primate species sampled across Southeast Asia.

‘Other’ includes Trachypithecus obscurus and undefined species from the genus Presbytis. Total surveys (n) = 148.

Overall, records report on a total of 6322 primates, with the largest proportion sampled from Peninsular Malaysia (48.5%, n=3069/6322). Primate surveys were primarily conducted on Long-tailed macaques (Macaca fascicularis) (90.5%, n=5720/6322) followed by Pig-tailed macaques (M. nemestrina; n=532/6322; Amir et al., 2020; Lee et al., 2011; Muehlenbein et al., 2015; Putaporntip et al., 2010; Appendix 1—table 3). Reported prevalence of Plasmodium knowlesi in NHPs ranged from 0% to 100%. Only 87 of the surveys (58.8%, n=87/148) reported a positive diagnosis, with the remaining 61 sites finding no molecular evidence of P. knowlesi infection (41.2%) in any primates tested. A full breakdown of P. knowlesi infection rates according to reported primate characteristics can be found in Appendix 1.

Meta-analysis of P. knowlesi prevalence

To quantify regional heterogeneity in simian cases of P. knowlesi, a one-stage meta-analysis of prevalence (number positive out of the number sampled) was conducted on primate malaria survey data. Overall pooled estimate for P. knowlesi prevalence was 11.99% (CI95% 9.35–15.26). Overall heterogeneity was assessed using the I2 statistic. Substantial between-study heterogeneity (I2 ≥75%) was found across all prevalence records (I2=80.5%; CI95% 77.3–83.1). In the sub-group analysis by region, pooled prevalence estimates are consistently low for Thailand (2.0%, CI95% 1.1–3.5%), moderate in Peninsular Malaysia (14.3%, CI95% 11.1–18.2) and elevated in Singapore (23.3%, CI95% 11.0–42.8) and Malaysian Borneo (41.1%, CI95% 20.8–64.9) (Figure 2). Sub-group heterogeneity was assessed using prediction intervals, derived from τ 2 statistic used to describe between-study variability. Prediction intervals indicate high heterogeneity of estimates within regions, consistent with expectations of high variability of prevalence across individual study sites. Detailed forest plots for individual prevalence estimates can be found in Appendix 3—figure 2.

Random-effects meta-analysis of P. knowlesi prevalence across Southeast Asia.

(A) Forest plot of pooled estimates for P. knowlesi prevalence (%) in all non-human primates tested (n=6322) across Southeast Asia, disaggregated by species and sampling site (k=148). Random-effects meta-analysis sub-grouped by region, with 95% confidence intervals and prediction intervals. (B) Map of regional prevalence estimates for P. knowlesi prevalence in NHP in Southeast Asia from meta-analysis. Point colour denotes pooled estimate (%). Size denotes total primates tested per region (n). Shading indicates data availability.

Risk factor analysis

Covariate data and P. knowlesi prevalence data were used to fit additional models to explore the relationships between localised landscape configuration and NHP malaria prevalence. Environmental covariates were extracted from satellite-derived remote sensing datasets (Table 1) at either true sampling sites (GPS coordinates) or 10 random pseudo-sampling sites to account for geographic uncertainty in prevalence data. Host species was grouped as ‘Macaca fascicularis’ or ‘Other’ due to sample counts of <10 for certain primate species. Only 57.4% (n=85/148 records) of data included year of sampling, deemed to be insufficient to assess temporal patterns in prevalence. Tree canopy cover ranged from negligible to near total cover (100%) within buffer radii (Appendix 4—table 2). Details of covariate data processing is illustrated in Appendix 4.

Table 1
Spatial and temporal resolution (res.) and sources for environmental covariates.

Summary metrics extracted within 5, 10 and 20 km circular buffers.

CovariateSpatial res.Temporal res.Source
Human density (p/km2)1 km2012WorldPop, 2018
Elevation (m)1 km2003SRTM 90 m Digital Elevation v4.1 Jarvis et al., 2008
Tree cover (1/0)*30 mAnnualHansen’s Global Forest Watch Hansen et al., 2013
  1. *

    Derivatives: Proportion canopy cover (%), Perimeter: area ratio (PARA >0)

Following a two-stage approach for selection of explanatory variables, tree cover and fragmentation (measured by perimeter: area ratio, PARA) were retained at 5 km as linear terms, human population density was retained at both 5 km and 20 km and primate species was retained as a categorical variable. Spearman’s rank tests for residual correlation between final variables at selected scales indicates a strong negative correlation between tree cover and fragmentation index (PARA; ρ = –0.75; Appendix 6—figure 2).

Adjusting for all other covariates in the model, we identified strong evidence of an effect between increasing tree canopy cover and higher prevalence of P. knowlesi in NHPs within a 5 km radius (aOR = 1.38, CI95% 1.19–1.60; p<0.0001). Evidence was also found for an association between likelihood of P. knowlesi and higher degrees of habitat fragmentation (PARA) within 5 km (aOR = 1.17, CI95% 1.02–1.34, p<0.0281). Evidence suggests that human population density within a 5 km radius is associated with risk of P. knowlesi in NHP (aOR = 1.36, CI95% 1.16–1.58, p=0.0001) whilst human density within 20 km has an inverse effect on likelihood of P. knowlesi (aOR = 0.56, CI95% 0.46–0.67, p<0.0001). M. fascicularis is also associated with higher prevalence relative to all other non-human primate species (aOR = 2.50, CI95% 1.31–4.85; p=0.0051). Additional complexity did not improve optimal model fit and effect modification was not pursued. In sensitivity analyses removing data points with excessive spatial uncertainty or restricting data points only to areas with high probability of macaque occurrence, evidence was consistently found that tree canopy cover (5 km) and host species exhibit a strong positive association with prevalence of P. knowlesi in NHP (Appendix 6). Final adjusted OR for the full multivariable model can be visualised in Figure 3.

Multivariable regression results.

Spatial scale denoted in square bracket. Canopy cover = %. Adjusted odds ratios (OR, dots) and 95% confidence intervals (CI95%, whiskers) for factors associated with P. knowlesi in NHPs at significant spatial scales. N=1354, accounting for replicate pseudo-sampling.

Discussion

Land use and land cover change is widely linked to spillover of zoonotic pathogens from sylvatic reservoirs into human populations, and pathogen prevalence in wildlife host species is key in driving the force of infection in spillover events. Our initial analyses found that for Plasmodium knowlesi, there is substantial spatial heterogeneity and prevalence in non-human primates varies markedly between regions of Southeast Asia (Zhang et al., 2016). Consistent with our hypothesis that parasite density in primate hosts would be higher in areas experiencing habitat disturbance, we identified strong links between P. knowlesi in NHPs and measures of contemporaneous tree cover and habitat fragmentation. To our knowledge, this is the first systematic study to find evidence of landscape influencing the distribution of P. knowlesi prevalence in NHPs. Results offer evidence that P. knowlesi infection rates in NHPs are linked to changes in landscape across broad spatial scales, and that prevalence of P. knowlesi in reservoir species may be driving spillover risk across Southeast Asia. These findings could provide insight to improving surveillance of P. knowlesi and to the development of ecologically targeted interventions.

While previous studies have estimated that P. knowlesi infection would be chronic in all macaques, or as high as 50–90% for modelling P. knowlesi transmission in Malaysia (Brock et al., 2016), this data strongly suggests that this is not the case. Overall prevalence of P. knowlesi infection in all NHPs is markedly lower than usual estimates, emphasising the importance of accounting for absence data in estimations of prevalence. Considerable heterogeneity was identified between and within regional estimates for P. knowlesi across Southeast Asia, which likely reflects genuine differences according to distinct climates and habitats (Shearer et al., 2016). Malaysian Borneo was found to have an estimated prevalence over five-fold higher than West Malaysia. Crucially, such extreme prevalence estimates for NHPs in Borneo align with the known hotspot for human incidence of P. knowlesi (Cooper et al., 2020). By comparison, for Peninsular Malaysia, estimated prevalence is far lower than anticipated. Cases of human P. knowlesi do occur in West Malaysia, although transmission has been found to exhibit spatial clustering (Phang et al., 2020) which may correspond to pockets of high risk within the wider context of low prevalence of P. knowlesi in macaque populations. Regional trends in P. knowlesi also mask differences in infection rates between sample locations, driven by more localised factors. Multiple studies reported finding P. knowlesi infections in wild macaques to be low or absent in peri-domestic or urbanised areas, attributed to the absence of vector species typically found in forest fringes (Brant et al., 2016; Chua et al., 2019; Manin et al., 2016). This pattern is seen in reports from Peninsular Malaysia (Saleh Huddin et al., 2019; Vythilingam et al., 2008), Singapore (Jeslyn et al., 2011; Li et al., 2021b), and Thailand (Fungfuang et al., 2020; Putaporntip et al., 2010). The high heterogeneity of reports here suggests that the picture is even more complex. P. knowlesi infections may even vary between troops within a single study site, as was seen in the Philippines (Gamalo et al., 2019). Fine-scale interactions are unlikely to be captured by the scale of this study.

Ecological processes determining P. knowlesi infection are influenced by dynamic variables over multiple spatial scales (Cohen et al., 2016). We utilised a data-driven methodology to select variables at distances that capture maximum impact on P. knowlesi prevalence (Byrne et al., 2021; Fornace et al., 2019b), with tree cover and fragmentation influential at localised scales and human population density also exerting influence within wider radii. Contrary to previous studies on risk factors for human incidence of P. knowlesi (Fornace et al., 2019b; Fornace et al., 2016), elevation was not found to be associated with P. knowlesi in NHPs at any scale. Vector and host species composition vary substantially across tropical ecotones, and it is likely that the study extent encompasses a range of putative vectors across different landscapes, such as those of the Minimus Complex in northern regions (Parker et al., 2015) or the recently incriminated An.-collessi and An.-roperi from the Umbrosus Group (De Ang et al., 2021). Given that the vector species driving sylvatic transmission remain elusive, it is conceivable that the elevation range covers multiple vector and host species niches and explains the lack of observed relationship between elevation and P. knowlesi in NHPs. Human population density was found to be significant at multiple distances, with contrasting effects on parasite prevalence in NHP. Previous studies have found a negative association between human density and vector density and biting rates in forested landscapes (Fornace et al., 2019a). Across wide spatial scales, increased vector density in less populated, more forested areas could generate higher parasite prevalence in NHPs. At the same time Long-tailed macaques, a species shown here to have higher prevalence rates, are notorious as nuisance animals and many of the available samples were collected opportunistically in urban areas, which might underly the observed positive association between localised high human density and higher prevalence in NHP. Whilst more data would be needed to understand this interaction, this further demonstrates the importance of using approaches to identify disease dynamics across multiple spatial scales (Brock et al., 2019).

A key finding is the link between high prevalence of P. knowlesi in primate host species with high degrees of habitat fragmentation. Habitat fragmentation is a key aspect of landscape modification, where large contiguous areas of habitat (for example, forests) are broken into a mosaic of smaller patches. This disturbs the ecological structure by increasing the density of fringes or ‘edges’, dynamic habitat often at the boundaries between natural ecosystems and human-modified landscapes (Borremans et al., 2019). Other studies have linked habitat fragmentation to increased generalist parasite density in primates. In Uganda, a higher prevalence and infection risk of protozoal parasites was observed in wild populations of red colobus primates (Procolobus rufomitratus) inhabiting fragmented forests compared to those in undisturbed habitat (Gillespie and Chapman, 2008). For P. knowlesi, creation of edge habitat is thought to favour vectors of the Leucosphyrus Complex (Davidson et al., 2019a; Hawkes et al., 2019). Anopheles spp. presence can be predicted by indices of fragmentation in Sabah, Borneo, with land cover changes creating more suitable micro-climate for larval habitats (Byrne et al., 2021), and an increased abundance of An. balabacensis found in forest fringes (Hawkes et al., 2019; Wong et al., 2015). Increasing landscape complexity results in increased density of edge habitat, with conceivably higher density of vectors in forest fringes. Therefore, preferential use of fringe habitat and high exposure to vectors in forest fringes may contribute to higher conspecific transmission of P. knowlesi between primates in increasingly fragmented habitats. This finding also lends clarity to landscape fragmentation as a risk factor for human exposure to P. knowlesi in Malaysian Borneo (Brock et al., 2019; Fornace et al., 2019b), with changes in relative host density, vector density and wildlife parasite prevalence in nascent forest fringes potentially enhancing the spillover of this disease system into human populations in fragmented habitats.

Conversely, we saw a strong association between high parasite prevalence and high tree canopy coverage. Given that a strong inverse relationship with fragmentation was observed, with high tree density correlating to low fragmentation indices and vice versa, this speaks to a trade-off between dense canopy cover and high habitat complexity and suggests an ‘ideal’ amount of habitat fragmentation that facilitates prevalence in primate hosts. For animals with larger home ranges, individual-based disease models combined with movement ecology approaches have shown that the most highly fragmented areas are less favourable for maintaining parasite transmission (White et al., 2018). In Sabah, individual macaques were shown to increase ranging behaviour in response to deforestation (Stark et al., 2019). Forest edge density also peaks at intermediate levels of land conversion (Borremans et al., 2019). With smaller habitat patches in maximally fragmented landscapes potentially insufficient to support macaque troops, this interplay between disease ecology and metapopulation theory may explain why both tree density and habitat fragmentation appear to pose a greater risk for simian P. knowlesi. Likewise, this may relate to the finding that in Borneo, larger forest patches (lower fragmentation indices) were associated with P. knowlesi spillover in Borneo (Fornace et al., 2019b). Overall, this finding offers an insight to mechanisms that underpin the increased force of infection of P. knowlesi that is associated with landscape change.

There are limitations to consider in the available data and interpretation of these findings. ‘Small-study effects’ were observed in the dataset, suggestive of a bias toward positive effect estimates (Stewart et al., 2012). This may be a result of data disaggregation and small studies creating artefactually higher estimates or may reflect true bias in data collection toward areas known to be endemic for P. knowlesi and convenience sampling of macaques. Assumptions have also been made that sample site equates to habitat, which may not reflect actual habitat use, and even accurate georeferenced data points are unlikely to entirely reflect surrounding habitat within the macaque home range. Variability in study designs and data reporting also impacted geospatial accuracy. Steps were taken to account for spatial bias by extracting covariates at randomly generated pseudo-sampling points. Whilst uncertainty cannot be eliminated, we demonstrate a robust methodology to accommodate for geographical uncertainty in ecological studies. Future investigations should prioritise systematic, georeferenced sampling across a range of landscape scenarios.

Results show important regional ecological trends, but broad geographic patterns may not be generalisable at individual levels, or to all putative host species in all geographic contexts (Zhang et al., 2016). Follow-up studies should be conducted at higher spatial and temporal resolution to characterise the effect of local landscape configuration on wildlife P. knowlesi prevalence. Effects of fragmentation are likely to be dependent on land conversion type, species composition and surrounding matrix habitat (Fornace et al., 2019b). Use of perimeter: area ratio (PARA) as a fragmentation index was justified given high canopy coverage in study sites (Wang et al., 2014), although Edge Density (ED) or normalised Landscape Shape Index (nLSI) might be more appropriate in future analyses to account for variation in forest abundance. Specific land configurations have previously been linked to P. knowlesi exposure in Borneo (Fornace et al., 2019b), notably in areas where palm oil plantation is a dominant industry. Given this, broad forest classifications used here may mask important differences in P. knowlesi prevalence between land classes. As it was not possible to include contemporary land cover classifications in this analysis, future studies would also benefit from looking at specific habitat type (e.g., primary forest, agroforest, plantation).

Concluding remarks

Strong links have been identified between land use and land cover change and ecosystem perturbation that favours the transmission of vector-borne diseases (Loh et al., 2016). Prevalence of P. knowlesi in macaques is likely to be a crucial determinant of human infection risk, and more representative estimates of P. knowlesi prevalence derived here can better inform regional transmission risk models. This study also characterises landscape risk factors for heightened prevalence of P. knowlesi in NHPs. Findings provide evidence that P. knowlesi in primate hosts is partly driven by landscape modification across Southeast Asia. While the full complexity is not captured by the covariates used, it is clear that P. knowlesi infection in NHPs is not restricted to densely forested areas. This study also demonstrates the utility of systematic meta-analysis tools and remote-sensing datasets in the investigation of macroecological disease trends, in conjunction with methods to standardise a spatially heterogeneous dataset and data-driven selection of spatial scales. Gaps identified in data reporting should inform more systematic and localised primatological surveys to disentangle precise mechanisms. Notwithstanding limitations, this study highlights the marked spatial heterogeneity and role of landscape complexity in driving P. knowlesi infection rates in NHPs. Given the clear intersection between human epidemiology and wildlife ecology, it is essential that infection dynamics within wildlife reservoirs are considered in future public health interventions.

Methods

Study site

This study focused on the simian malaria Plasmodium knowlesi across Southeast Asia, within 28°30'00.0"N, 92°12'00.0"E and 11°00'00.0"S, 141°00'00.0"E. Climate mainly corresponds to the equatorial tropical zone, with high temperatures and high humidity.

Data assembly

A systematic literature review was conducted under the CoCoPop framework (Condition, Context, Population) (Ruiz Cuenca et al., 2022; Munn et al., 2015). All studies identified in the literature review were screened for data on NHPs with a confirmed P. knowlesi diagnosis or absence data (zero counts of P. knowlesi with appropriate diagnostic methods). Exclusion criteria included (a) studies exclusively relying on microscopy (Antinori et al., 2013) (b) laboratory, animal model or experimental infection studies (c) data from outside of Southeast Asia. No limit was set on the temporal range for primate survey records. Duplicate records reporting results from the same surveys were removed, with one record per survey retained. Critical appraisal of the studies was conducted using the Joanna Briggs Institute (JBI) checklist for prevalence studies (Munn et al., 2015; see Appendix 1 for details and criteria). A flowchart of the selection process is illustrated in Appendix 1—figure 3, with a full list of articles included provided in Appendix 1—table 2.

Primary outcome was defined as P. knowlesi prevalence (p, proportion positive for P. knowlesi infection from n sampled NHPs). For each independent primate study, the following variables were extracted: year of data collection, primate species sampled, primate status (wild/captive), diagnostic test (PCR/sequencing) and target gene(s), sampling method (routine/purposive), number of P. knowlesi positive samples, number of Plasmodium spp. positive samples, total number of primates tested and geographical information.

In most studies identified, study site was only geolocated to a geographic area or descriptive location. Geolocation was assigned at the lowest available level of administrative polygon (i.e. district/state/country) by cross-referencing reported sampling location with GADM (v3.6) administrative boundaries. If specific location was given, GPS coordinates were assigned via Google Maps. For data visualisation, point coordinates were plotted in QGIS (3.10.14) and R (4.1.0) software.

Meta-analysis of P. knowlesi prevalence

Meta-analysis was conducted using methods that are standard in the analysis of human disease prevalence for individual participant datasets (IDP) (Liberati et al., 2009; Stewart et al., 2012). Data were disaggregated by geographic location (site) and primate species, to illustrate variance in prevalence by survey unit (Stewart et al., 2012). One-stage meta-analysis is considered appropriate for studies where the outcome may be infrequent, so data was included in a single model under the ‘DerSimonian and Laird’ variance estimator (Munn et al., 2015). Sensitivity analyses were conducted to compare methods for the back-transformation of prevalence estimates. For studies where prevalence estimates tend towards 0% or 100%, variance tends towards 0. To stabilise the variance and enable back-transformation of zero prevalence records, logit method was selected for the transformation of prevalence, with the inverse variance method used for individual study weights (see Appendix 3 for details).

Overall heterogeneity of prevalence records was assessed using the I2 statistic (von Hippel, 2015), a relative estimate of true between-study variance. Sub-group analysis was conducted according to geographic region, with the heterogeneity of reported prevalence within regional sub-groups assessed using prediction intervals derived from the τ 2 statistic. Small-study effects, including selection and publication biases, were assessed by examining funnel plots and imputing ‘missing’ estimates using the trim-and-fill method (Lin and Chu, 2018). Full rationale and details of small-study effect assessments can be found in Appendix 3.

Remote sensing data

Satellite-derived remote sensing datasets were used to assemble local environmental and anthropogenic covariates. Gridded UN-adjusted human population estimates were assembled at 1 km resolution from WorldPop, 2018. Elevation data was obtained from NASA SRTM 90 m Digital Elevation Database v4.1 (CGIAR-CSI) (Jarvis et al., 2008) with a spatial resolution of 1 km. Contemporaneous tree cover was derived from Hanson’s Global Forest Watch (30 m) (Hansen et al., 2013), extracted for every year between 2006 and 2020.Tree cover was classified as ≥50% crown density, and then matched to primate data by sample site geolocation and by year of sample collection to account for rapid forest loss (Appendix 4—figure 1). Where a broad timeframe of sampling was provided (≥3 years), median year was used. Full details for variable selection and processing can be found in Appendix 4.

Perimeter: area ratio (PARA, ratio of patch perimeter length to patch surface area) of given land class is a key metric for habitat conversion, where a higher PARA provides a measure of boundary complexity and indicates a more fragmented landscape (McGarigal et al., 2021). Mean PARA was extracted from canopy cover within circular buffers. Habitat fragmentation has been shown to correlate with disease transmission parameters (Borremans et al., 2019; Faust et al., 2018), but definitions often lack precision and can be considered with respect to ‘separation effects’ (division and isolation of patches) and ‘geometric effects’ (changes to ratios of perimeter and core habitat; Wilkinson et al., 2018). PARA provides a measure of edge density within the buffer area (PARA >0) and has been shown to provide a good index of fragmentation and good discrimination of spatial aggregation across areas where habitat abundance (tree canopy cover) is high (Wang et al., 2014; Appendix 4—table 2, Appendix 4—figure 4).

Covariate assembly

For studies with exact GPS coordinates, precise environmental data at a single site could be obtained. For surveys published without GPS coordinates, there is considerable geographic uncertainty in the exact sampling location (Appendix 5). Uncertainty in the spatial and environmental determinants of prevalence generates a sampling bias, with the precision of covariates correlated to certain studies. Use of a single centroid proxy site is standard procedure, but often generates erroneous estimates in large or heterogenous sampling units (Cheng et al., 2021). Alternative strategies were employed to account for and mitigate the effect of spatial uncertainty and spatial bias. Each prevalence observation was replicated and assigned a random sample of environmental realisations. 10 random sampling points were generated within the sampling area provided by the study, and covariates were extracted at each proxy sampling site (Appendix 5—figure 1). Selection of random points was validated by visual inspection of the stability of model coefficients with the inclusion of an increasing number of points. Number of points was selected conservatively at the point where coefficients stabilised (n=10).

For every georeferenced sampling point, mean values for all selected covariates were extracted within buffer radii at 5 km, 10 km, and 20 km (Appendix 4). Buffer area sizes were selected to investigate multiple spatial scales over which associations between risk factors and P. knowlesi prevalence might occur. A minimum radius of 5 km was chosen to approximate the maximum ranging distance for M. fascicularis (Waxman et al., 2014), with wider radii (10–20 km) included to account for the geographic uncertainties in areal data. Flowchart of data processing chain can be found in Appendix 4—figure 2.

Analysis of environmental risk factors

Generalised linear mixed-effect regression models (GLMM) were fitted to NHP prevalence data using a binomial distribution with a logit link. To account for within-study correlation in reported average prevalence, a unique identifier combining author and study was included as a random intercept in all models. Artificial inflation of sample size in the replicated data (10 pseudo-sampling sites for data geolocated to administrative areas) was accommodated by reducing individual observation weights to 1/10th within the model.

Each covariate at each spatial scale was assessed for inclusion in the multivariable model based on bivariable analysis and a criterion of p>0.2 under likelihood ratio tests (LRT; Appendix 6—table 1). A quadratic term for the fragmentation index ‘PARA’ was included to account for possible nonlinearity. Multicollinearity among independent predictors at multiple scales was examined via variance inflation factors (VIF). The VIF of each predictor variable was examined following a stepwise procedure, starting with a saturated model and sequentially excluding the variable with the highest VIF score from the model. Stepwise selection continued in this manner until the entire subset of explanatory variables in the global model satisfied a moderately conservative threshold of VIF ≤6 (Rogerson, 2001). Qualifying variables obtained were then assessed for model inclusion using a backward stepwise strategy, removing variables with the highest p value (LRT) until a pre-defined threshold of α<0.05. Spearman’s rank tests were conducted on the selected variables to observe residual correlation, plotted as a correlation matrix (Appendix 6—figure 2).

Fully adjusted odds ratios (OR) for associations between environmental covariates and P. knowlesi prevalence were derived from the final multivariable GLMM with p values derived from LRT. Spatial sensitivity analyses were conducted by excluding data points from administrative boundaries outside a reasonable size or above a reasonable threshold of environmental certainty, according to the standard deviation (SD) of the covariate values within each set of 10 environmental realisations (Appendix 5—figure 2, Appendix 6—tables 4–6). Ecological sensitivity analyses were conducted by removing data points that fall outside areas with high predicted probability of occurrence for Macaca fascicularis, Macaca nemestrina, and Macaca leonina and running regression analyses on the constrained dataset (Moyes et al., 2016; Appendix 6—figures 35, Appendix 6—tables 7–9).

Appendix 1

Data assembly

Prior to conducting the study, a review of current literature was constructed to find articles related to ‘Plasmodium knowlesi’ or to both ‘malaria’ and ‘primate’, including synonyms and sub-headings. The search was elaborated to specify environmental factors (Appendix 1—figure 1). The following databases were searched:

  • Medline

  • Embase

  • Web of Science

Provisional data were extracted using a standardised form (Appendix 1—figure 2) using standardised definitions (Appendix 1—table 1), from which an initial set of studies were identified for this investigation. Duplicate records (confirmed/suspected to be the same specimens) were removed, with one record retained (Appendix 1—figure 3)

Appendix 1—figure 1
Search strategy for background research.
Appendix 1—figure 2
WHO Report primate data extraction form.
Appendix 1—figure 3
Flow chart illustrating study selection process.
Appendix 1—table 1
Standardised definitions for qualitative primate characteristics.
VariableCategoryDefinition
SamplingRoutineAnimals collected for surveillance purposes or extracted from human–conflict zones; data collected opportunistically Shearer et al., 2016
StudyAnimals captured and sampled specifically for a study of Plasmodium spp and/or P. knowlesi
StatusCaptiveAnimal resident in sanctuary or conservation park
WildFree-living animal, not registered/resident in any sanctuary
SanctuaryA wildlife sanctuary/rehabilitation centre housing key primate species Gamalo et al., 2019
AreaForestAreas that are uninhabited with extensive tree cover
Peri–domesticAs defined by the author. Example definitions as follows:
AgriculturalAnimal located in agricultural areas, predominantly monoculture (e.g. orchard, plantation) Shahar, 2019
UrbanAs defined by the author. Generally, areas with high human population density Li et al., 2021b
Appendix 1—table 2
Characteristics of the included studies.
AuthorYear(s)Country/regionN*SampleDiagnosticTarget gene(s)Primer
Lee et al., 20112004–2008Malaysia/Borneo108StudyNested PCRSSU-rRNA/csp/mtDNAKn1f/Kn3r
Seethamchai et al., 20082006Thailand99StudySequencingA-type-SSU-rRNA/cytb
Vythilingam et al., 20082007Malaysia/Peninsular145StudyPCR/SequencingSSU-rRNA/cspPmk8/Pmkr9
Zhang et al., 20162007Singapore40StudyPCR
2007–2010Indonesia/Sumatra70StudyPCR
2011Cambodia54StudyPCR
2012Philippines68StudyPCR
2015Laos44StudyPCR/SequencingSSU-rRNAPK18SF/PK18SRc
Jeslyn et al., 20112008Singapore13RoutinePCR/SequencingSSU-rRNA /cspPmk8/Pmkr9
Ho et al., 20102008Malaysia/Peninsular107RoutineNested PCRSSU-rRNAPmk8/Pmkr9
Li et al., 2021b2008–2017Singapore1039RoutineNested PCRSSU-rRNAPmk8/Pmkr9
Putaporntip et al., 20102009Thailand655StudySequencingcytb
Chang et al., 20112010Myanmar45StudyPCRSSU-rRNA
Muehlenbein et al., 20152010Malaysia/Borneo41StudyPCRmtDNA/AMA-1/MSP-1
Shahar, 2019 2010–2017Malaysia/Peninsular1587RoutineNested PCRSSU-rRNA
§2013Malaysia/Peninsular15StudyPCR
§2013–2016Malaysia/Borneo25StudyNested PCRcytbPKCBF/PKCBR
Saleh Huddin et al., 20192014Malaysia/Peninsular415StudyPCR/SequencingSSU-rRNAPmk8/Pmkr9
Akter et al., 20152015Malaysia/Peninsular70RoutinePCR/SequencingA-type-SSU-rRNAPmk8/Pmkr9
Amir et al., 20202016Malaysia/Peninsular103RoutineNested PCRSSU-rRNAPkF1140/PkR1550
Gamalo et al., 20192017Philippines95StudyNested PCRSSU-rRNAKn1f/Kn3r
Fungfuang et al., 20202017–2019Thailand93StudyNested PCRSSU-rRNAKn1f/Kn3r
Nada-Raja et al., 20222018Malaysia/Borneo73StudyNested PCRSSU-rRNA/csp/mtDNAKn1f/Kn3r
Yusuf et al., 20222016–2019Malaysia419StudyNested PCRSSU-rRNAKn1f/Kn3r
Zamzuri et al., 20202018Malaysia/Peninsular212RoutinePCR
Kaewchot et al., 20222019Thailand649StudyNested PCRSSU-rRNAPmk8/Pmkr9
Salwati et al., 20172015Indonesia/Sumatra38StudyPCR/Sequencing
  1. *

    N=number of primates sampled.

  2. Animal trapped either on routine or study basis.

  3. Unpublished, personal correspondence (p/c).

  4. §

    Danau Girang Field Centre, p/c from Dr Salgado Lynn.

Of the 87 records reporting presence of P. knowlesi, only 22 records (containing 248 P. knowlesi positive macaques) report whether P. knowlesi infection was a mono-infection or mixed infection with other simian Plasmodium spp. With a low proportion of data represented, this was deemed insufficient to conduct any further investigations.

Macaca fascicularis is the predominant species tested. However, reports also include M. nemestrina (6.1%, n=9/148; 527 macaques) (Amir et al., 2020; Lee et al., 2011; Putaporntip et al., 2010; Muehlenbein et al., 2015), M. arctoides (1.4%, n=2/148; 36 macaques) (Fungfuang et al., 2020; Putaporntip et al., 2010), M.-leonina (n=1/148; 25 macaques) (Fungfuang et al., 2020), Trachypithecus obscurus (Dusky leaf monkey) (n=1/148; 7 tested) and unspecified species from the Presbytis genus (n=1/148; 2 tested) (Appendix 1—table 3). One study additionally sampled 1 Presbytis melalophos (Black-crested Sumatran langur) (Vythilingam et al., 2008), but species-specific P. knowlesi was not reported (Appendix 1—table 4).

Appendix 1—table 3
Published studies of P. knowlesi infections (n) in non-human primate species collected (N) in Southeast Asia, grouped by region and author.
RegionSpeciesTotalRef
M. fascicularisM. nemestrinaM. arctoidesOther
Peninsular25/107473/3069Ho et al., 2010
Malaysia48/415Saleh Huddin et al., 2019
21/70Akter et al., 2015
11/980/5Amir et al., 2020
0/15
10/145Vythilingam et al., 2008
215/1587Shahar, 2019
66/415Yusuf et al., 2022
74/2073/5Zamzuri et al., 2020
Borneo4/262/15119/251Muehlenbein et al., 2015
71/8213/26Lee et al., 2011
18/25
7/452/28Nada-Raja et al., 2022
2/4Yusuf et al., 2022
Sumatra0/706/108Zhang et al., 2016
5/321/40/2Salwati et al., 2017
Thailand1/1955/4490/41/7 8/1496Putaporntip et al., 2010
0/21Seethamchai et al., 2008
0/4Fungfuang et al., 2020
0/320/25*1/32Fungfuang et al., 2020
0/78Seethamchai et al., 2008
0/649Kaewchot et al., 2022
Philippines18/9518/163Gamalo et al., 2019
0/68Zhang et al., 2016
Singapore3/13148/1092Jeslyn et al., 2011
145/1039Li et al., 2021b
0/40Zhang et al., 2016
Laos1/441/44Zhang et al., 2016
Cambodia0/540/54Zhang et al., 2016
Myanmar0/450/45Chang et al., 2011
Total743/572026/5571/361/9773/6322
  1. *

    Macaca-leonina (Northern Pig-tailed macaque, recently classified as separate species).

  2. Presbytis spp.

  3. Trachypithecus obscurus (Dusky leaf monkey).

Appendix 1—table 4
Characteristics of primates tested and number/percentage of confirmed P. knowlesi infections (Pk+).
N(%)*Pk+Pk+ (%)CI95%
SpeciesM. fascicularis5720(90.5%)74513.0%(12.2–13.9)
M. nemestrina532(8.4%)264.9%(3.4–7.1)
M. leonina25(0.4%)00.0%(0.0–13.3)
M. arctoides36(0.6%)12.8%(0.5–14.2)
T. obscurus7(0.1%)114.3%(2.6–51.3)
Presbytis spp.2(0.03%)00.0%(0.0–65.8)
AreaForest1740(27.5%)25314.5%(13.0–16.3)
Agriculture491(7.8%)7214.7%(11.8–18.1)
Peri-domestic2192(34.7%)34115.6%(14.1–17.1)
Urban1143(18.1%)564.9%(3.8–6.3)
Sanctuary109(1.7%)54.6%(2.0–10.3)
Unspecified647(10.2%)467.1%(5.4–9.4)
StatusWild6183(97.8%)76812.4%(11.6–13.3)
Captive139(2.2 %)53.6%(1.5–8.1)
RegionPen. Malaysia3069(48.5%)47315.4%(14.2–16.7)
Borneo251(4.0%)11947.4%(41.3–53.6)
Sumatra108(1.7%)65.5%(2.6–11.6)
Thailand1496(23.7%)80.5%(0.3–1.1)
Philippines163(2.6 %)1811.0%(7.1–16.8)
Singapore1092(17.3%)14813.6%(11.7–15.7)
Cambodia54(0.9%)00.0%(0.0–6.6)
Laos44(0.7%)12.3%(0.4–11.8)
Myanmar45(0.7%)00.0%(0.0–7.9)
Total6322(100%)77312.2%(11.4–13.1)
  1. *

    Percentage of total number of primates tested (column %).

  2. Proportion of N positive for P. knowlesi (row %).

  3. 95% confidence interval (CI95%) calculated in R using count and sample size (binomial distribution).

Appendix 2

Quality appraisal

Quality was assessed using the Joanna Briggs Institute (JBI) Critical Appraisal tool for prevalence studies (Munn et al., 2015). Studies were assessed on nine standardised criteria used to inform inclusion in the meta-analysis. Full criteria and examples of scoring are given in Appendix 2—tables 1 and 2. Studies assessed to be of lower quality were those that omitted key information about the sampling method . Two studies were considered to be of higher quality owing to robust sampling and completeness of evidence (Lee et al., 2011; Saleh Huddin et al., 2019). Given the objective to assess variation in reported prevalence, and in considering the limited usefulness of criteria designed for human participants, the reliable diagnostic methods identified in all studies and the appreciable limitations in surveying wild animals, data from all studies (n=148 estimates) were included for further analyses.

Appendix 2—table 1
JBI criteria for assessing bias in meta-analyses of prevalence studies.
CriteriaYesNoUnclearN/A
Q1Was the sample frame appropriate to address the target population?
Q2Were study participants sampled in an appropriate way?
Q3Was the sample size adequate?
Q4Were the study subjects and the setting described in detail?
Q5Was the data analysis conducted with sufficient coverage of the identified sample?
Q6Were valid methods used for the identification of the condition?
Q7Was the condition measured in a standard, reliable way for all participants?
Q8Was there appropriate statistical analysis?
Q9Was the response rate adequate, and if not, was the low response rate managed appropriately?
Appendix 2—table 2
Example rationale for quality appraisal.
Sample questionExampleAssessment
Q1Was the sample frame appropriate to address the targetWild animalYes
population?Captive animalNo
Not specifiedUncertain
Q2Were study participants sampled in an appropriate way?Trapped for studyYes
Routine collectionNo
Not specifiedUncertain
Q9Was the response rate adequate?Primate dataN/A

Appendix 3

Meta analysis

Sub-group analysis by region was conducted under a random-effects model. Pooled estimates are then back-transformed for interpretation. The Freeman-Tukey double arcsine method is recommended in the transformation of prevalence (Munn et al., 2015). However, recent studies have found that back-transformation of the Freeman-Turkey method can generate misleading results, owing to the requirement for a global sample size for inversion (Schwarzer et al., 2019). A sensitivity analysis conducted using the logit transformation and untransformed proportions (Appendix 3—table 1) revealed a deficit in the back-transformation of the pooled prevalence estimate for Thailand under the Freeman-Tukey transformation, generating a null point estimate. To avoid this error, meta-analysis was conducted under the logit transformation with the inverse variance estimator to account for individual study weighting.

Appendix 3—table 1
Sensitivity analysis for transformation of P. knowlesi prevalence estimate under random-effects model, shown overall and for Thailand subgroup analysis.
MethodOverall (k=148)Subgroup (Thailand, k=21)
P*CI95%PCI95%
Freeman-Turkey double arcsine0.0943(0.0641–0.1284)0.0000(0.0000–0.0000)
Logit0.1199(0.0935–0.1526)0.0199(0.0113–0.0346)
Untransformed0.1415(0.1101–0.1730)0.0022(0.0000–0.0059)
  1. *

    Estimated proportion

Small study effects, including selection and publication bias, were assessed by examining funnel plots and imputing ‘missing’ estimates using the trim-and-fill method (Lin and Chu, 2018). Funnel plots were generated using both SE and study size as metrics of variance, as study size has been shown to be more accurate for meta-analyses of proportions where raw estimates tend towards 0 or 1 (Munn et al., 2015). Funnel plots for the disaggregated dataset are shown using SE as the variance estimate in Appendix 3—figure 1. Asymmetry is highlighted by the trim-and-fill interpolation method, which provides an estimate of missing data and added an additional 54 imputed points to the plot.

Appendix 3—figure 1
Assessment of small study effects in meta-analysis.

(A) Funnel plot of transformed prevalence (%) against standard error (SE) for study sites (B) Funnel plot with imputed data to illustrate asymmetry using trim-and-fill method.

Meta-analysis was conducted with data disaggregated by survey location and primate species (k=148). Forest plot of individual study prevalence, presented with pooled regional prevalence estimates and relative sampling effort for the disaggregated data can be visualised in Appendix 3—figure 2.

Appendix 3—figure 2
Forest plot of P. knowlesi prevalence (%) in all species of NHPs in Southeast Asia, disaggregated by species and sampling site, including 95% confidence intervals and individual study weighting.

Random-effects analysis, sub-grouped by region. N=148. Zarith et al. refers to personal correspondence derived from the reference Shahar, 2019.

Appendix 4

Remote sensing data and covariate assembly

Environmental covariates were extracted from satellite-derived remote-sensing datasets. Elevation can be used as a proxy for vector range, with malaria transmission patterns often correspond to altitudinal ranges of vector species (Hawkes et al., 2019). Estimates of human population density provide a measure of urbanisation, used to examine risks related to human settlement proximity. Metrics of deforestation including canopy cover (cumulative loss) and degree of fragmentation, which are key determinants of macaque habitat selection (McGarigal et al., 2021) and of mosquito vector breeding sites and were derived from Hansen’s Global Forest Watch (Hansen et al., 2013) Land use classification maps derived from the Intact Forest Landscape (IFL) (Potapov et al., 2008) and Copernicus Global Land Service (100 m) (Buchhorn et al., 2020) were also explored to provide more detailed information on specific landscape composition. However, as 77.0% of NHP records were collected before the earliest land classification in 2015 (N=114/148), datasets were of limited utility and not pursued further.

Appendix 4—figure 1
Recent forest loss in Peninsular Malaysia (first row) and Malaysian Borneo (secoond row), shown for the years (A) 2006 (B) 2012 and (C) 2019.

Gridded UN-adjusted human population estimates were assembled at 1 km resolution from WorldPop, 2018 for multiple timepoints as a measure of urbanisation, a proxy for risks related to human settlement proximity. As minimal variation between timepoints was observed, only 2012 (median year of primate data collection) was retained.

Appendix 4—table 1
Environmental covariates assembled for regression analysis.

Summary values extracted for each covariate within 5, 10, and 20 km circular buffers during processing.

CovariateDescriptionMetricResolutionProcessingSource
SpatialTemporal
PopulationUN-adjusted gridded posterior population model estimates at 30 arc-seconds resolutionPerson count/1 km21 km2000 2012 2019Population density reclassified as high/low (≤300 persons/km2) in QGISWorldPop, 2018 Downloaded as tiff files per country in AOI for years 2000/2012/2019
ElevationMean height above sea levelm1 km2003Mean and SD of continuous elevation per radii. Mean-centred and scaled. Categorised into discrete classifications: low (≤200 m), moderate (200–500 m) or high elevation (>500 m)NASA SRTM 90 m Digital Elevation Database v4.1 (CGIAR-CSI) Jarvis et al., 2008. Downloaded as a tiff file at 1 km resampled resolution
ForestPercentage canopy cover per grid cell. Derived from tree cover (vegetation >5 m) and loss (forested to non-forested)0–130 mAnnual 2006–2020Tree cover classified as ≥50% crown density per raster cell, generating binary raster (1=forest, 0=non-forest). Annual cover calculated by subtracting cumulative loss per year 2006–2019. Data records matched to reclassified tile by geolocation and year. Posterior proportions categorised as high (>50%) medium (20–50%) or low (≤20%)Hansen’s Global Forest Watch, 30 m resolution Landsat imagery Hansen et al., 2013. Tiles downloaded as tiff files for each year 2006–2019 to cover AOI
Fragmentation (perimeter: area ratio, PARA)Perimeter length (m) to patch area (m2) ratio for contiguous forest cover McGarigal et al., 2021 within bufferPARA >030 mAnnual 2006–2020Extracted from annual reclassified tree cover rasters within 5, 10 and 20 km circular buffers Output categorised into quartilesHansen’s Global Forest Watch Hansen et al., 2013 (as above)
Appendix 4—figure 2
Flowchart of data processing.

Details of pseudo-sampling and environmental covariate extraction at multiple spatial scales to create final weighted dataset (N=1354).

Appendix 4—figure 3
Example covariate resolutions in Peninsular Malaysia.

(A) Data point and (B) 20 km buffer over population density layer, 1 km resolution. (C) Data point and (D) 20 km buffer over SRTM elevation layer, 1 km resolution.

Elevation, population density and forest cover all varied markedly across surveyed sites. Forest cover ranges from negligible to near total cover within 5 km, and up to 99.96% and 99.64% within 10 km and 20 km respectively (Appendix 4—table 2, Appendix 4—figure 4). Within a 20 km buffer, 46.1% of sites have dense forest cover ≥50% (n=683/1480) and 83.85% have moderate or high forest cover (≥20%), with similar distributions over 5 km and 10 km. Example buffers over forest cover data can be visualised in Appendix 4—figure 5.

Appendix 4—table 2
Summary of forest cover data (N=1480).
MeanSDRange
Forest cover (5 km)50.20%±29.29%0.00–100.00%
Forest cover (10 km)49.68%±27.30%0.00–99.96%
Forest cover (20 km)48.29%±25.35%0.00–99.64%
Total1480 (100%)
Appendix 4—figure 4
Boxplots showing distribution and interquartile range (IQR) of proportional forest cover (0–1) for sampling sites within 5, 10 and 20 km circular buffers across all sites (N=1354).
Appendix 4—figure 5
Examples of buffer zones around macaque sample sites.

Shown over forest cover for 2019 (Hansen et al., 2013).

Appendix 5

Spatial uncertainty

Available spatial resolution of the survey sites varied. 14 records (9.5%, n=14/148) could be geolocated to a point, using geographic coordinates provided or inferred. The remaining 134 were geolocated to the lowest administrative polygon according to GADM boundary definitions (Appendix 5—table 1).

Appendix 5—table 1
Geo-positioning of available primate survey data.
Resolution/GADM*Records/nPrimates/NMin. (km2)Max. (km2)
PolygonCountry/GID06 (4.9%)853 (17.3%)70077,650
State/GID140 (22.0%)2699 (32.2%)13087,860
District/GID288 (61.8%)2433 (43.6%)27015,890
PointGCS 14 (11.4%)337 (6.8%)
Total148 (100%)4931 (100%)
  1. *

    Administrative boundaries, as classified by GADM (v3.6).

  2. Minimum and maximum size (km2) of polygons containing P. knowlesi data at each admin level.

  3. Geographic Coordinate System.

Crude sensitivity analyses were initially conducted to evaluate use of centroids vs random points to approximate macaque survey site. Using GADM classifications, the largest polygon containing NHP data was identified at each administrative level. 10 points were randomly generated within each polygon in QGIS, with buffers at 5/10/20 km. Proportion of forest cover per buffer was extracted, categorised and compared to the forest cover for the centroid (Appendix 5—figure 1).

Appendix 5—figure 1
Sensitivity analysis comparing centroid forest cover to 10 randomly generated points, shown per radius for the largest polygon at each GADM level.

Results indicate that at district level (GID2), minimal change in forest cover between the points was observable (5 km: 0.68–0.96). However, for both the state of Southern Sumatra (GID1 5 km: 0.26–0.97) and for Myanmar (GID0 5 km: 0.16–0.99) variation was observed, with several points classifying as ‘moderate’ (20–50% cover) rather than ‘high’ (>50%) as suggested by the centroid. Overall, results show a disparity between covariates obtained at a central point compared with random points within larger administrative polygons, indicating inadequate sensitivity of centroids as a proxy for local landscape variables where there is spatial uncertainty (Cheng et al., 2021).

Appendix 5—figure 2
Boxplots of standard deviation in repeat sampling of covariates at multiple buffer and boundary sizes.

Standard deviation of environmental covariates across 10 sampling site realisations within 5/10/20km buffers, grouped by administrative boundary size or GPS coordinates. Shown for (A) canopy cover (%) (B) forest fragmentation (P: A ratio) and (C) human population density (p/km2).

Appendix 6

Regression analysis

Appendix 6—table 1
Bivariable analysis for P. knowlesi in NHP against all covariates at all spatial scales (N=1354).
Bivariable analysis
VariableCrude ORCI95%p value
Elevation (m) *
≤5 km1.18(1.07–1.28)0.000562
≤10 km1.20(1.09–1.31)0.0001246
≤20 km1.22(1.11–1.33)7.23E-05
Human density (p/km2) *
≤5 km0.84(0.77–0.92)4.72E-05
≤10 km0.75(0.68–0.82)2.70E-12
≤20 km0.71(0.63–0.79)1.37E-12
Forest cover (%) *
≤5 km1.34(1.21–1.49)1.86E-08
≤10 km1.41(1.26–1.57)6.51E-10
≤20 km1.47(1.30–1.67)8.66E-10
Fragmentation (PARA) *
≤5 km0.85(0.76–0.95)0.003944
≤10 km0.69(0.60–0.80)1.80E-07
≤20 km0.67(0.57–0.79)5.14E-07
PARA2 *Quadratic term
≤5 km0.69(0.60–0.80)0.102.06E-06
≤10 km0.64(0.55–0.74)0.082.65E-08
≤20 km0.67(0.57–0.78)0.032.78E-06
Host species
OtherRef
M. fascicularis2.37(1.25–4.60)0.007971
  1. *

    Continuous variable, mean-centred and scaled.

  2. p value derived from Likelihood ratio test (LRT).

Appendix 6—figure 1
Spearman’s correlation matrix for all candidate covariates at all spatial scales (n=1354).
Appendix 6—figure 2
Spearman’s correlation matrix for covariates at selected spatial scales for final model inclusion (n=1354).

Percentage forest cover (5 km) and forest fragmentation (PARA, 5 km) show strong negative correlation (ρ=–0.75).

Appendix 6—table 2
Multivariable binomial regression analysis of P. knowlesi prevalence in NHP with environmental covariates at influential spatial scales, full dataset (N=1354).

AIC = 1229.8.

Multivariable analysis
VariableRadiusaOR *CI95%p value
Human density (p/km2)
≤5km1.36(1.16–1.58)1.082E-04
≤20 km0.56(0.46–0.67)1.311E-10
Forest cover (%)
≤5 km1.38(1.19–1.60)2.046E-05
Fragmentation (PARA)
≤5 km1.17(1.02–1.34)0.0281
Host species
OtherRef
M. fascicularis2.50(1.31–4.85)0.005121
  1. *

    Odds Ratios adjusted for all other variables in the table (aOR). Radius calculated as distance from sample point.

  2. p value derived from Likelihood ratio test (LRT).

  3. Continuous variable, mean-centred and scaled. OR shown per 1 SD increase.

For observations with high geographic uncertainty, random pseudo-sampling of 10 sites (as described in the manuscript and Appendix 4—figure 2) was used to avoid overgeneralisations and biases typical when centroids are used as proxy sampling sites35, with data weighted accordingly. However, this has the potential to introduce extreme or unrealistic values by generating points in landscapes that are outside reasonable estimations of primate study site. Given this, further sensitivity analyses were conducted to validate the results against geographic precision of observations. Data were first truncated to include only data geolocated to administrative boundaries for relatively small area size (see Appendix 6—table 3) and exclude highly variable data from country-level boundaries (GID0). Singapore was retained as small administrative unit. GLMM regression models were fit to the truncated dataset.

Appendix 6—table 3
Admin boundary sensitivity analysis.

Binomial regression analysis of P. knowlesi prevalence in NHP for datapoints assigned to GPS or small sized administrative boundaries (excluding country data) (N=1324).

AIC = 1221.9Multivariable analysis
aORCI 95%p value (Wald test)*
Human density [5 km]1.36(1.16–1.58)***
Human density [20 km]0.56(0.46–0.67)***
Forest cover (%) [5 km]1.38(1.19–1.60)***
Fragmentation (PARA) [5 km]1.18(1.02–1.34)*
Host group Other
M. fascicularis
REF
2.51
(1.–.31–4.85)**
  1. *

    Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1.

Administrative boundaries are arbitrary categories and vary considerably in size and landscape consistency. To better evaluate environmental uncertainty associated with each observation, standard deviation (SD) of the covariate values within each set of 10 environmental realisations was calculated (resulting in a single standard deviation value for each covariate at each scale for every prevalence data point). Studies for which the uncertainty (SD) of covariates exceeded half of the maximum standard deviation were censored to avoid spurious associations derived from unreliable/extreme values for both forest cover (5 km) and fragmentation (5 km). Regression models were fit to the winsorized dataset and compared to results from the full dataset to ensure that associations are robust.

Appendix 6—table 4
Distribution of standard deviations across 10 environmental covariates per prevalence data point for landscape variables at all spatial scales (N=1354).
CovariateMeanRangeMedianIQR
Canopy [5 km]*0.15880.0000–0.42370.15340.1039–0.2277
Canopy [10 km]0.13190.0000–0.41240.11770.0827–0.1891
Canopy [20 km]0.10510.0000–0.39990.09210.0541–0.1635
Fragmentation [5 km]*0.00830.0000–0.03750.00630.0041–0.0094
Fragmentation [10 km]0.00610.0000–0.04550.00430.0025–0.0071
Fragmentation [20 km]0.00410.0000–0.03160.00270.0017–0.0047
  1. *

    Spatial scales selected in final variables.

Appendix 6—table 5
Tree canopy cover sensitivity analysis.

Binomial regression of P. knowlesi prevalence in NHP for datapoints, with data where SD < ½ the maximum for tree canopy within 5 km (N=814).

AIC = 771.1Multivariable analysis
aORCI 95%p value (Wald test)*
Human density [5 km]0.90(0.67–1.20)-
Human density [20 km]0.72(0.50–1.01).
Forest cover (%) [5 km]1.70(1.30–2.24)***
Fragmentation (PARA) [5 km]1.38(1.01–1.88)*
Host group Other
M. fascicularis
REF
2.63
(1.35–5.21)**
  1. *

    Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘-’ 1.

Appendix 6—table 6
Landscape fragmentation sensitivity analysis.

Binomial regression of P. knowlesi prevalence in NHP for datapoints, with datapoints where SD < ½ the maximum for fragmentation within 5 km (N=1134).

AIC = 982.9Multivariable analysis
aORCI 95%p value (Wald test)*
Human density [5 km]0.91(0.69–1.18)-
Human density [20 km]0.69(0.52–0.93)*
Forest cover (%) [5 km]1.31(1.08–1.60)**
Fragmentation (PARA) [5 km]1.18(0.90–1.54)-
Host group Other
M. fascicularis
REF
2.52
(1.33–4.87)**
  1. *

    Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1.

Sensitivity analysis has shown that the trends are robust when the data is constrained according to small administrative boundaries or by measures of spatial uncertainty in the environmental variables. However, given that a proportion of points are randomly generated, questions remain about how suitable the resulting sites are for macaque species and consequently whether the associations observed are a realistic indication of ecological trends. To address this, points were subset according to macaque species habitat suitability maps, derived from Moyes et al., 2016.

Predicted occurrence maps were combined for three species Macaca fascicularis, Macaca nemestrina and Macaca leonina to create a joint macaque extent for Southeast Asia. Binary maps of predicted habitat extent were then generated using thresholds of moderate (predictions of 0.5 and above), high (>0.75) and very high predicted probability of occurrence (>0.9) (Moyes et al., 2016; Appendix 6—figures 35). Datapoints from the main analysis were then overlayed with each map, and any points (±5 km buffer) that occurred outside of predicted habitat extent for macaque species were removed. Regression models were then fit to the reduced datasets to assess whether associations observed are plausible according to macaque ecology (Appendix 6—tables 7–9).

Appendix 6—figure 3
Distribution and habitat range of dominant macaque species (M. fascicularis, M. nemestrina, M. leonina) according to predicted probability of occurrence ≥0.5 (on a scale of 0–1.0) per 5x5 km pixel.
Appendix 6—table 7
Macaque habitat suitability sensitivity analysis.

Binomial regression of P. knowlesi prevalence in NHP for datapoints, including only datapoints with 5 km buffers that intersect with areas with ≥0.5 probability of predicted macaque occurrence (N=1331).

AIC = 1197.2Multivariable analysis
aORCI 95%p value (Wald test)*
Human density [5 km]1.32(1.13–1.54)***
Human density [20 km]0.55(0.45–0.66)***
Forest cover (%) [5 km]1.30(1.12–1.52)***
Fragmentation (PARA) [5 km]1.12(0.97–1.29)-
Host group Other
M. fascicularis
REF
2.48
(1.31–4.82)**
  1. *

    Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1.

Appendix 6—figure 4
Distribution and habitat range of dominant macaque species (M. fascicularis, M. nemestrina, M. leonina) according to predicted probability of occurrence ≥0.75 (on a scale of 0–1.0) per 5x5 km pixel.
Appendix 6—table 8
Macaque habitat suitability sensitivity analysis.

Binomial regression of P. knowlesi prevalence in NHP for datapoints, including only datapoints with 5 km buffers that intersect with areas with ≥0.75 probability of predicted macaque occurrence (N=1177).

AIC = 1115.5Multivariable analysis
aORCI 95%p value (Wald test)*
Human density [5 km]1.34(1.14–1.58)***
Human density [20 km]0.57(0.47–0.69)***
Forest cover (%) [5 km]1.23(1.04–1.47)*
Fragmentation (PARA) [5 km]1.04(0.86–1.24)-
Host group Other
M. fascicularis
REF
2.69
(1.38–5.38)**
  1. *

    Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1.

Appendix 6—figure 5
Predicted distribution and habitat range of all macaque species (M. fascicularis, M. nemestrina, M. leonina) according to predicted probability of occurrence ≥0.9 (on a scale of 0–1.0) per 5x5 km pixel.
Appendix 6—table 9
Macaque habitat suitability sensitivity analysis.

Binomial regression of P. knowlesi prevalence in NHP for datapoints, including only datapoints with 5 km buffers that intersect with areas with ≥0.9 probability of predicted macaque occurrence (N=567).

AIC = 685.2Multivariable analysis
aORCI 95%p value (Wald test)*
Human density [5 km]1.86(1.49–2.32)***
Human density [20 km]0.36(0.26–0.49)***
Forest cover (%) [5 km]1.47(1.14–1.90)**
Fragmentation (PARA) [5 km]1.35(1.02–1.77)*
Host group Other
M. fascicularis
REF
3.13
(1.50–6.75)**
  1. *

    Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1.

Data availability

All data of P. knowlesi prevalence in non-human primates were sourced from systematic literature review of publicly available research articles. Published primate datasets are available in Ruiz Cuenca et al., 2022 and Shahar, 2019. Full reference list for all publicly available primate infection data including details of individual articles and data types are available in Appendix 1. Earth Observation (EO) data were derived from open-source environmental datasets. References, links and details of scale and resolution and manipulation are provided in manuscript (Table 1 in Results and details in Methods) with further information provided in the Appendix 4. Example code for manipulation of environmental data and data analysis are available on an associated GitHub repository (copy archived at Johnson, 2024).

References

  1. Book
    1. Brant HL
    (2011)
    Changes in Abundance, Diversity and Community Composition of Mosquitoes Based on Different Land Use in Sabah, Malaysia
    Imperial College London.
    1. Gilhooly L
    2. McIntyre A
    3. Grigg M
    4. Fornace K
    5. Cox J
    6. Drakeley C
    7. Lynn MS
    8. Stark D
    (2016)
    The Monkeybar Project: population density of long-tailed macaques (Macaca fascicularis) in two different forest types in Kudat District, Sabah, Malaysia
    American Journal of Physical Anthropology 159:155.
  2. Conference
    1. Ho GC
    2. Lee CL
    3. Sharma RSK
    (2010)
    Association of Institutions for Tropical Veterinary Medicine (AITVM)
    14th International Conference of the Association of Institutions for Tropical Veterinary Medicine.
    1. Loh EH
    2. Murray KA
    3. Nava A
    4. Daszak P
    (2016)
    Evaluating the links between biodiversity, land-use change, and infectious disease emergence
    Tropical Conservation pp. 79–88.
  3. Thesis
    1. Shahar ZS
    (2019)
    Molecular prevalence, spatial distribution and epidemiological risk factors of Plasmodium knowlesi infection among Macaca fascicularis raffles in Peninsular Malaysia
    Universiti Putra Malaysia.

Article and author information

Author details

  1. Emilia Johnson

    1. School of Biodiversity, One Health and Veterinary Medicine, University of Glasgow, Glasgow, United Kingdom
    2. Department of Disease Control, London School of Hygiene & Tropical Medicine, London, United Kingdom
    3. Centre on Climate Change and Planetary Health, London School of Hygiene & Tropical Medicine, London, United Kingdom
    Contribution
    Conceptualization, Data curation, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing - original draft, Writing – review and editing
    For correspondence
    emilia.johnson@glasgow.ac.uk
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-5460-1715
  2. Reuben Sunil Kumar Sharma

    Faculty of Veterinary Medicine, Universiti Putra Malaysia, Selangor, Malaysia
    Contribution
    Conceptualization, Resources, Data curation, Supervision, Writing – review and editing
    Competing interests
    No competing interests declared
  3. Pablo Ruiz Cuenca

    1. Department of Disease Control, London School of Hygiene & Tropical Medicine, London, United Kingdom
    2. Lancaster University, Bailrigg, Lancaster, United Kingdom
    3. Liverpool School of Tropical Medicine, Pembroke Place Liverpool, Liverpool, United Kingdom
    Contribution
    Data curation, Supervision, Validation, Investigation, Methodology, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-2180-9509
  4. Isabel Byrne

    Department of Disease Control, London School of Hygiene & Tropical Medicine, London, United Kingdom
    Contribution
    Formal analysis, Methodology, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-7800-3733
  5. Milena Salgado-Lynn

    1. School of Biosciences, Cardiff University, Cardiff, United Kingdom
    2. Wildlife Health, Genetic and Forensic Laboratory, Sabah Wildlife Department, Wisma Muis, Kota Kinabalu, Malaysia
    3. Danau Girang Field Centre, Sabah Wildlife Department, Kinabalu Sabah, Malaysia
    Contribution
    Conceptualization, Resources, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-1769-6465
  6. Zarith Suraya Shahar

    Faculty of Veterinary Medicine, Universiti Putra Malaysia, Selangor, Malaysia
    Contribution
    Resources, Data curation
    Competing interests
    No competing interests declared
  7. Lee Col Lin

    Faculty of Veterinary Medicine, Universiti Putra Malaysia, Selangor, Malaysia
    Contribution
    Resources, Data curation
    Competing interests
    No competing interests declared
  8. Norhadila Zulkifli

    Faculty of Veterinary Medicine, Universiti Putra Malaysia, Selangor, Malaysia
    Contribution
    Resources, Data curation
    Competing interests
    No competing interests declared
  9. Nor Dilaila Mohd Saidi

    Faculty of Veterinary Medicine, Universiti Putra Malaysia, Selangor, Malaysia
    Contribution
    Resources, Data curation
    Competing interests
    No competing interests declared
  10. Chris Drakeley

    Department of Infection Biology, London School of Hygiene & Tropical Medicine, London, United Kingdom
    Contribution
    Conceptualization, Funding acquisition, Investigation, Methodology, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-4863-075X
  11. Jason Matthiopoulos

    School of Biodiversity, One Health and Veterinary Medicine, University of Glasgow, Glasgow, United Kingdom
    Contribution
    Conceptualization, Formal analysis, Validation, Investigation, Methodology, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-3639-8172
  12. Luca Nelli

    School of Biodiversity, One Health and Veterinary Medicine, University of Glasgow, Glasgow, United Kingdom
    Contribution
    Formal analysis, Validation, Investigation, Visualization, Methodology
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-6091-4072
  13. Kimberly Fornace

    1. School of Biodiversity, One Health and Veterinary Medicine, University of Glasgow, Glasgow, United Kingdom
    2. Department of Disease Control, London School of Hygiene & Tropical Medicine, London, United Kingdom
    3. Centre on Climate Change and Planetary Health, London School of Hygiene & Tropical Medicine, London, United Kingdom
    4. Saw Swee Hock School of Public Health, National University of Singapore, Singapore, Singapore
    Contribution
    Conceptualization, Resources, Formal analysis, Supervision, Funding acquisition, Validation, Investigation, Methodology, Project administration, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-5484-241X

Funding

Wellcome Trust (221963/Z/20/Z)

  • Kimberly Fornace

Royal Society (221963/Z/20/Z)

  • Kimberly Fornace

Ministry of Higher Education Malaysia (LRGS/1/2018/UM/01/1)

  • Reuben Sunil Kumar Sharma

World Health Organization

  • Kimberly Fornace
  • Chris Drakeley

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication. For the purpose of Open Access, the authors have applied a CC BY public copyright license to any Author Accepted Manuscript version arising from this submission.

Acknowledgements

Research was supported by the Sir Henry Dale Fellowship, jointly funded by the Wellcome Trust and the Royal Society (Grant Number 221963/Z/20/Z). Data collected in Peninsular Malaysia by the Faculty of Veterinary Medicine, Universiti Putra Malaysia, was funded by the Ministry of Higher Education Malaysia (Grant LRGS/1/2018/UM/01/1/). Additional funding was supported by the World Health Organization.

Version history

  1. Preprint posted:
  2. Sent for peer review:
  3. Reviewed Preprint version 1:
  4. Reviewed Preprint version 2:
  5. Reviewed Preprint version 3:
  6. Version of Record published:

Cite all versions

You can cite all versions using the DOI https://doi.org/10.7554/eLife.88616. This DOI represents all versions, and will always resolve to the latest one.

Copyright

© 2023, Johnson et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 1,107
    views
  • 138
    downloads
  • 2
    citations

Views, downloads and citations are aggregated across all versions of this paper published by eLife.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Emilia Johnson
  2. Reuben Sunil Kumar Sharma
  3. Pablo Ruiz Cuenca
  4. Isabel Byrne
  5. Milena Salgado-Lynn
  6. Zarith Suraya Shahar
  7. Lee Col Lin
  8. Norhadila Zulkifli
  9. Nor Dilaila Mohd Saidi
  10. Chris Drakeley
  11. Jason Matthiopoulos
  12. Luca Nelli
  13. Kimberly Fornace
(2024)
Landscape drives zoonotic malaria prevalence in non-human primates
eLife 12:RP88616.
https://doi.org/10.7554/eLife.88616.4

Share this article

https://doi.org/10.7554/eLife.88616

Further reading

    1. Ecology
    Mercury Shitindo
    Insight

    Tracking wild pigs with GPS devices reveals how their social interactions could influence the spread of disease, offering new strategies for protecting agriculture, wildlife, and human health.

    1. Ecology
    2. Neuroscience
    Ralph E Peterson, Aman Choudhri ... Dan H Sanes
    Research Article

    In nature, animal vocalizations can provide crucial information about identity, including kinship and hierarchy. However, lab-based vocal behavior is typically studied during brief interactions between animals with no prior social relationship, and under environmental conditions with limited ethological relevance. Here, we address this gap by establishing long-term acoustic recordings from Mongolian gerbil families, a core social group that uses an array of sonic and ultrasonic vocalizations. Three separate gerbil families were transferred to an enlarged environment and continuous 20-day audio recordings were obtained. Using a variational autoencoder (VAE) to quantify 583,237 vocalizations, we show that gerbils exhibit a more elaborate vocal repertoire than has been previously reported and that vocal repertoire usage differs significantly by family. By performing gaussian mixture model clustering on the VAE latent space, we show that families preferentially use characteristic sets of vocal clusters and that these usage preferences remain stable over weeks. Furthermore, gerbils displayed family-specific transitions between vocal clusters. Since gerbils live naturally as extended families in complex underground burrows that are adjacent to other families, these results suggest the presence of a vocal dialect which could be exploited by animals to represent kinship. These findings position the Mongolian gerbil as a compelling animal model to study the neural basis of vocal communication and demonstrates the potential for using unsupervised machine learning with uninterrupted acoustic recordings to gain insights into naturalistic animal behavior.