Genomic epidemiology of Plasmodium knowlesi reveals putative genetic drivers of adaptation in Malaysia

Abstract

Sabah, Malaysia, has amongst the highest burden of human Plasmodium knowlesi infection in the world, associated with increasing encroachment on the parasite's macaque host habitat. However, the genomic make-up of P. knowlesi in Sabah was previously poorly understood. To inform on local patterns of transmission and putative adaptive drivers, we conduct population-level genetic analyses of P. knowlesi human infections using 52 new whole genomes from Sabah, Malaysia, in combination with publicly available data. We identify the emergence of distinct geographical subpopulations within the macaque-associated clusters using identity-by-descent-based connectivity analysis. Secondly, we report on introgression events between the clusters, which may be linked to differentiation of the subpopulations, and that overlap genes critical for survival in human and mosquito hosts. Using village-level locations from P. knowlesi infections, we also identify associations between several introgressed regions and both intact forest perimeter-area ratio and mosquito vector habitat suitability. Our findings provide further evidence of the complex role of changing ecosystems and sympatric macaque hosts in Malaysia driving distinct genetic changes seen in P. knowlesi populations. Future expanded analyses of evolving P. knowlesi genetics and environmental drivers of transmission will be important to guide public health surveillance and control strategies.

Fig 1. Comparable within-isolate genetic complexity across geographic regions.

A: Boxplots depicting the distribution of within-infection diversity (F_WS) across three regions (Sabah, Sarawak, and Peninsular). B: Dot plots of the within-isolate non-reference allele frequencies (NRAF) across the genome for six Malaysian P. knowlesi infections ranging from low diversity (F_WS=0.99) to high diversity (F_WS=0.54) and with varying levels of within-host relatedness.

Fig 2. Geographic overlap and evidence of shared ancestry between the Mf and Mn P. knowlesi clusters.

A: Unrooted neighbour-joining tree based on identity-by-state (IBS) depicting three predominant genomic clusters of P. knowlesi across Malaysia, specifically the Peninsular Malaysia sub-population (Peninsular), and Malaysian-Borneo macaque-associated subpopulations of Macaca fascicularis (Mf) and M. nemestrina (Mn). The new isolates from Sabah are labelled separately (Sabah). B: Map of Malaysia showing the geographic distribution and number of samples, and genomic clusters across Malaysian-Borneo (right) and Peninsular-Malaysia (left). C: Bar plot illustrating the proportionate ancestry to each of 3 (K) subpopulations determined by ADMIXTURE for each sample (bars on x-axis), sectioned by geographic region. The three K populations identified aligned perfectly with the clustering in the NJ tree; K=1 with Mf, K=2 with Mn and K=3 with Peninsular as per the colour-coding. Shapefile made with Natural Earth: https://www.naturalearthdata.com/downloads/50m-cultural-vectors/50m-admin-1-states-provinces/.

Fig 3. Identity-by-descent (IBD)-based analysis reveals greater relatedness amongst Mn and Peninsular than Mf samples.

Each circle reflects an infection, colour-coded by genomic clustering group, and the number of lines between infections reflects relatedness (more lines reflect greater relatedness) at the given connectivity threshold of minimum IBD. Where two circles are not connected by a line, the estimated IBD between those infections was below the given threshold. The three samples from Peninsular Malaysia with >95% IBD represent laboratory-based strains from the 1960s that have been passaged through macaques (SRR2222335, SRR2225467 & SRR3135172).

Fig 4. Infection connectivity is partly driven by geography at the state-level administrative boundary.

A & B: Identity-by-descent (IBD)-based cluster network illustrating the distant relatedness for samples within Mf (A) and Mn (B) clusters collected in two adjacent states; Sabah and Sarawak, at different cut-offs for the proportion of IBD in a paired comparison. C: Map of East Malaysia on the island of Borneo with colours representing the two states being compared in the IBD analysis. Shapefile made with Natural Earth: https://www.naturalearthdata.com/downloads/50m-cultural-vectors/50m-admin-1-states-provinces/.

Fig 5. Genetic differentiation reveals candidate adaptations.

Genome-wide scans of differentiation between Sabah and Sarawak subpopulations within the (A) Mf and (B) Mn clusters using the between-population fixation index (F_ST). Only the Mf cluster shows clear differences in diversity with peaks of differentiation clear at several chromosomes (most notably on chromosomes 8, 11 and 12), whilst the Mn cluster has substantial ‘noise’ across the genome, with high levels of differentiation across the genome.

Fig 6. Determination of introgressed windows between the Mf and Mn P. knowlesi clusters.

A: Dot and contours plot describing potential introgression events within an Mf sample, where x and y axes represent the genetic distance of the sample to the Mf and Mn clusters, respectively. Genetic distance is the proportion of mismatched SNPs per sliding window (10kb) when comparing the called allele in the sample to the major allele for a cluster at each position. The contours represent the density of genetic distances for the three clusters. B: Dot and contours plot of the same sample above, subset to those windows deemed to be introgressed from the Mn cluster. Possible introgression events are sliding windows that fall outside the major contours of the samples own cluster and within the major contours of another, representing greater similarity in genetic distance to the other cluster. C: Unrooted neighbour-joining tree based on identity-by-state (IBS) of window 1504 on chromosome 08 (950000-959999) and overlapping the PKA1H_080026000 gene (encodes the oocyst capsule protein). The Mf samples/branches clustering within the Mn branches (depicted by asterisk) provides further evidence that introgression of this window has occurred in these samples. D: SNP barcode plot of window 1504 on chromosome 08 (950000-959999) showing greater genetic similarity between several Mf samples (depicted by asterisk) and the Mn cluster, where the colours reflect those in the legend on panel C, and the alpha represents the allele.