Dihydrofolate-Reductase Mutations in Plasmodium knowlesi Appear Unrelated to Selective Drug Pressure from Putative Human-To-Human Transmission in Sabah, Malaysia

Background Malaria caused by zoonotic Plasmodium knowlesi is an emerging threat in Eastern Malaysia. Despite demonstrated vector competency, it is unknown whether human-to-human (H-H) transmission is occurring naturally. We sought evidence of drug selection pressure from the antimalarial sulfadoxine-pyrimethamine (SP) as a potential marker of H-H transmission. Methods The P. knowlesi dihdyrofolate-reductase (pkdhfr) gene was sequenced from 449 P. knowlesi malaria cases from Sabah (Malaysian Borneo) and genotypes evaluated for association with clinical and epidemiological factors. Homology modelling using the pvdhfr template was used to assess the effect of pkdhfr mutations on the pyrimethamine binding pocket. Results Fourteen non-synonymous mutations were detected, with the most common being at codon T91P (10.2%) and R34L (10.0%), resulting in 21 different genotypes, including the wild-type, 14 single mutants, and six double mutants. One third of the P. knowlesi infections were with pkdhfr mutants; 145 (32%) patients had single mutants and 14 (3%) had double-mutants. In contrast, among the 47 P. falciparum isolates sequenced, three pfdhfr genotypes were found, with the double mutant 108N+59R being fixed and the triple mutants 108N+59R+51I and 108N+59R+164L occurring with frequencies of 4% and 8%, respectively. Two non-random spatio-temporal clusters were identified with pkdhfr genotypes. There was no association between pkdhfr mutations and hyperparasitaemia or malaria severity, both hypothesized to be indicators of H-H transmission. The orthologous loci associated with resistance in P. falciparum were not mutated in pkdhfr. Subsequent homology modelling of pkdhfr revealed gene loci 13, 53, 120, and 173 as being critical for pyrimethamine binding, however, there were no mutations at these sites among the 449 P. knowlesi isolates. Conclusion Although moderate diversity was observed in pkdhfr in Sabah, there was no evidence this reflected selective antifolate drug pressure in humans.


Introduction
Zoonotic transmission of P. knowlesi from its simian hosts to humans has likely been occurring since the human population first expanded in Southeast Asia around 40,000 years ago [1,2]. Despite microscopy misdiagnosis of P. knowlesi with other human Plasmodium species [3,4], accurate molecular detection has validated recent reports of an increasing incidence of P. knowlesi malaria cases in Eastern Malaysia [5][6][7], and the continued presence of low P. knowlesi endemicity in other countries in Southeast Asia [8][9][10][11][12][13][14]. The large geographical distribution of P. knowlesi is constrained only by the range of its natural macque hosts and Leucosphyrus group vector [15].
Human-to-human (H-H) transmission may also have been present since humans were first infected, with studies conducted in the 1960s demonstrating that Anopheles balabacensis (the most common malaria vector in Sabah, Malaysia [16]) can transmit P. knowlesi competently between human hosts [17]. Despite this, evidence to date suggests that P. knowlesi transmission remains primarily zoonotic. P. knowlesi circumsporozoite protein (CSP) gene and mitochondrial genome sequence data from studies in the Malaysian state of Sarawak [1], Peninsular Malaysia [18], Singapore [19] and Thailand [10] have found no evidence of divergence between the infections in human versus macaque hosts. In addition, ongoing P. knowlesi transmission has not been reported in areas where the long-tailed or pig tailed macaque hosts [20] are not present, consistent with recent modelling indicating a very low likelihood of sustained H-H transmission together with a low human reproductive rate (R OH = 1.04) [21]. However, with an increasing incidence of knowlesi malaria in Malaysia [6], along with evidence of peri-domestic transmission [22], it is possible that H-H transmission is occurring to at least some degree. Molecular markers, such as monitoring for effects of selective antimalarial drug pressure on parasite populations, may provide a novel indicator, and would allow evaluation of spatio-temporal trends as seen in mapping of antimalarial drug-resistant genotypes in other areas [23,24]. Resistance to antifolate medications such as pyrimethamine arises readily in P. falciparum [25,26] and P. vivax [27][28][29] and has been well documented throughout Southeast Asia, including Malaysia [30,31]. In Malaysia, proguanil was deployed from the late 1940s, and resistance documented almost immediately. Pyrimethamine resistance was documented from the mid-1970s onwards [32], although sulfadoxine-pyrimethamine (SP) continued to be used as a first-line medication for the treatment and prophylaxis of P. falciparum for over 30 years until artemisinin-combination therapy (ACT) was introduced after 2009. The use of single dose pre-referral SP for falciparum malaria continues in remote areas of Sabah and Sarawak in Eastern Malaysia. As P. knowlesi is commonly misdiagnosed as P. falciparum due to morphological similarities in the early ring stages [3,7], and as co-infections of P. knowlesi with P. falciparum or P. vivax do occur [33], it is highly likely P. knowlesi parasite populations in humans would have been exposed to SP in Malaysia. Additional antifolate drug exposure from medications such as trimethoprim, used commonly for non-malarial infections, has demonstrated pyrimethamine cross-resistance in P. falciparum [34].
Pyrimethamine acts on Plasmodium parasites by inhibiting the dihydrofolate-reductase (dhfr) enzyme involved in the folate biosynthesis pathway. Non-synonymous point mutations in the dhfr gene affecting drug-enzyme binding confer resistance to treatment [25]. In P. falciparum dhfr (pfdhfr) this is manifest firstly with the acquisition of the S108N mutation, followed by increasing resistance associated with the addition of mutations located at residues 51, 59 and 164 [35]. For P. vivax dhfr (pvdhfr) orthologous mutations to pfdhfr are located at codons 50, 58, 117, and 173 respectively [36] although these may not confer the same degree of functional resistance to pyrimethamine [37]. Sequencing of pkdhfr to date has only been reported from clinical isolates in the Andaman and Nicobar Islands in India, which have a very low P. knowlesi endemicity predominantly consisting of mixed infections with either P. falciparum or P. vivax [38]. Analysis revealed discordant mutations between Plasmodium species with only wild type pkdhfr reported together with multiple mutations of pfdhfr or pvdhfr, consistent with zoonotic-only transmission of P. knowlesi [38].
To assess selection pressure provided by antifolates as a possible surrogate marker of naturally occurring H-H transmission we undertook P. knowlesi dhfr sequencing in 449 human malaria cases in Sabah, Malaysia. As only human hosts would be expected to have SP drug exposure, the presence of similar functional pkdhfr mutation genotypes in spatio-temporal clusters would support the likelihood of H-H transmission. Sustained transmission of P. knowlesi genotypes with dhfr mutations would suggest that competent transmission from humans back to the monkey host reservoir is occurring naturally, and hence increase the possibility of H-H transmission. In addition it was hypothesised that clinical correlates of human adapted P. knowlesi parasites would be associated with pkdhfr mutations. Early neurosyphillis studies demonstrated repeated blood passage of P. knowlesi through humans resulted in higher parasite counts and risk of severe disease [39,40]. This is consistent with P. knowlesi adapting to invade a wider age range of human RBCs over time in an in vitro culture system [41], with hyperparasitaemia a known independent predictor of severe knowlesi malaria [42]. P. knowlesi invasion gene variants Pknbpxa and Pknbpxb have also been associated with increased risk of hyperparasitaemia and manifestations of severe disease [43]. A study also described an admixture of two distinct parasite sub-populations arising from the separate long-tailed and pigtailed macaque hosts, including in samples from Sabah, with increased hybridisation between these two sub-populations potentially associated with parasite adaptation for humans [44]. Evaluating H-H transmission remains problematic but is important in the context of malaria elimination goals in South-East Asia, particularly due to the inability to manage the simian reservoir of P. knowlesi with conventional public health measures.

Study sites and patient enrolment
The study was conducted at Queen Elizabeth Hospital (QEH), Kudat District Hospital (KDH) and Kota Marudu District Hospital (KMDH). QEH is an adult tertiary-referral hospital located in Sabah's capital city Kota Kinabalu, and serves as a referral hospital for the West Coast and Kudat Division of Sabah, comprising 6 district hospitals (including KDH and KMDH) with a catchment population of 1.14 million. Patients with severe malaria are generally referred from district hospitals to QEH. KDH and KMDH are located in adjacent districts in Kudat Division, northeast Sabah, with a combined catchment area of 3200 square kilometres and a population of 150,000. Details of the three study sites and referral practices have been described elsewhere [42,45].
Patients at QEH were enrolled from September 2010 -June 2014 alongside a prospective study involving consecutive non-pregnant patients 12 years old who were admitted with PCR-confirmed malaria, were within 18 hours of commencing malaria treatment, had no major co-morbidities or concurrent illness, and had not been previously enrolled [46]. For the current study, patients with a PCR-confirmed P. knowlesi monoinfection were included, in addition to 48 randomly selected patients with a P. falciparum monoinfection. At KDH and KMDH non-pregnant patients 1 year old admitted with PCR-confirmed knowlesi malaria were enrolled from December 2012 -June 2014 as part of separate prospective studies described in detail elsewhere [45,47].

Study procedures
Baseline epidemiological and clinical information were recorded on standardised forms. Severe knowlesi malaria was defined according to modified WHO 2010 criteria for severe falciparum malaria, as previously described [46]. Pre-treatment blood slides were obtained on enrolment and read by an experienced research microscopist, with parasite counts per microlitre reported [48].
PCR amplification and sequencing of pkdhfr and pfdhfr DNA was extracted from 100 μL of packed red blood cell pellets using the QIAamp1 DNA Blood Kit (Qiagen, Australia) according to the manufacturer's instructions. Plasmodium species were confirmed as described previously [49,50]. A nested PCR approach using the same assay and cycling conditions for both primary and nested PCR was used to amplify pfdhfr and pkdhfr (i.e., 95°C for 3 min, followed by 25 cycles of 95°C for 30 sec, 52°C for 90 sec, and 72°C for 90 sec). The primary PCR was performed in a 50 μL reaction buffer volume, the nested PCR in a 100 μ L reaction buffer volume containing: 200 μM each deoxyribonucleotide triphosphate (dNTP; Bioline, Australia), 3  ]. The primary PCR was run using 2.5 μL DNA; 5 μL primary PCR was used in the nested PCR reaction. DNA sequencing of the nested PCR amplicons was performed by Macrogen (Seoul, Korea) [S1 File]. Sequences were analysed using Chromas Pro software (Technolysium Ltd., Australia) and BioEdit Sequence Alignment Editor [51] against the reference P. knowlesi genome [52]. Polyclonal infections, as represented by double peaks in the electropherogram, were observed in 7 infections (1x 34R/L 1x 52L/V, 1x 91P/A, 2x 91P/T, and 2x 149A/V) and categorized as mutant alleles [35].

DNA sequence diversity and neutrality tests
DnaSP (version 5.10.01) was used to calculate the total number of polymorphic sites, singletons and haplotypes in the pkdhfr and pfdhfr fragments. Tajima's neutrality test (D test) was also implemented with DnaSP software [53]. The nucleotide diversity at synonymous (dS), non-synonymous (dN), and all nucleotide sites (π) was calculated using the Nei-Gojobori method [54] with Jukes-Cantor correction [55] as implemented in the MEGA software (version 6.06). The standard error was determined by 1,000 bootstrap replications, and the rates of synonymous versus non-synonymous substitutions were compared using the Z-test of selection (MEGA 6.06).

3D pkdhfr protein modelling
Pkdhfr amino acid sequence was subjected to Basic Local Alignment Search (BLAST) [56], and further alignment was performed using ClustalW [57]. The sequence with maximum identity of 82.5%, pvdhfr, (PDB ID: 2BL9) was used as a template for homology modeling using SWISS MODEL server in http://swissmodel.expasy.org [58][59][60]. The constructed model was validated by PROCHECK [61]. The effects of both non-synonymous orthologous and nearby pkdhfr point mutations found in the current study on the pyrimethamine drug binding pocket were then evaluated by Autodock Vina [62].

Mapping of geographic distribution of pkdhfr genotypes
All patients with knowlesi malaria enrolled at district sites had their household and central village locations determined using a hand-held global positioning system (GPS; Garmin TM ; model 62sc) by local research fieldworkers. Patients enrolled at QEH had the central point of their village located and GPS coordinates recorded using Google Earth (version 7.1.2; Landsat images, July 2013), with assistance and corroboration via government public health and 2010 census maps, local fieldworkers, and input from patients where necessary. Final mapping of pkdhfr genotypes at village locations was performed using Quantum GIS software (version 2.4.0, www.qgis.org).

Spatial and temporal clustering
Analysis of spatial, temporal and combined space-time clustering of pkdhfr genotypes was performed with SaTScan TM software (version 9.2, http://satscan.org/), as previously described [60,63]. Data was inputted in a 0/1 event Bernoulli model with a case file including the number of patients with the specific pkdhfr genotype of interest at each village location and date, and a control file with the number, location and date of all other pkdhfr genotypes (the background distribution population). This software then uses the Kulldorf spatial [64], temporal or retrospective space-time scan statistic [65] to assess whether the specified pkdhfr genotypes are randomly distributed in space and/or time. Both large and small clusters are detected, with up to 50% of the study population allowed in a potential cluster within a maximum 30km radius. The time period for the temporal and space-time analysis was defined as 30 days, based on the average time of parasite maturation in the An. leucosphyrus vector (10 days) [66], the human pre-patent period (9-12 days) [17], and median fever duration before presentation (5 days) [46].

Statistical analysis
Data were analysed using Stata TM (version 12) software. Continuous variables were compared using Student t test or Mann-Whitney test depending on distribution, and proportions were compared using Chi squared or Fisher's exact test. The Cuzick's test for trend was used to assess the change in proportion of pkdhfr mutations over time, with enrolments at QEH (September 2010 -June 2014) divided into 5 x 9-month time-periods and enrolments at the district hospitals (December 2012 -June 2014) divided into 3 x 6-month time-periods.

Study population
From September 2010 -June 2014, 504 patients with PCR-confirmed P. knowlesi monoinfection were enrolled at the 3 study sites, with pkdhfr sequence data available from 449 (89%) patients. Of those with pkdhfr sequence data, 263 (59%) patients were enrolled at QEH and 186 (41%) were enrolled at KDH/KMDH. Baseline demographics are summarised in Table 1.
Overall, 78% of patients were male. Median age was 44 years at QEH adult hospital (range 14-94 years) and 30 years at KDH/KMDH (range 1-78 years). At KDH/KMDH, 40 (18%) children <12 years old were enrolled, with sequencing of the pkdhfr gene performed on 29 of these patients. From all study sites there were 190 (42%) farmers or plantation workers. Nearly all patients (96%) reported some forest or plantation exposure during the preceding 6 weeks, and 55% reported having seen a monkey during this time. Thirteen patients enrolled into the study at KDH/KMDH were subsequently transferred to QEH because of severe malaria. Sequencing of the pkdhfr gene was performed on all of these patients.

DNA sequence diversity and tests of neutrality
After excluding sequences with missing data, a total of 446 pkdhfr sequences and 45 pfdhfr sequences were available for analysis of sequence diversity. High levels of diversity were observed across the pkdhfr fragment, with 113 sites displaying evidence of single nucleotide variation within Sabah ( Table 4). A total of 188 haplotypes were observed. The overall haplotype diversity (Hd) and nucleotide diversity (π) was estimated to be 0.983 and 0.006. The pairwise nucleotide diversity was significantly higher at synonymous (dS = 0.018) versus nonsynonymous (dN = 0.003) sites (p = 0.007). The Tajima's D value for the pkdhfr region was -2.24 (p<0.001), reflecting an excess of low frequency polymorphism relative to expectation. In marked contrast to pkdhfr, only 2 pfdhfr sites displayed evidence of polymorphism within Sabah (108N and 59R were fixed in the population), both representing non-synonymous variants (dN = 0.001 and dS = 0). A total of 3 haplotypes were observed, with accordingly low haplotype diversity (Hd = 0.244). Tajima's D value was -0.820 but the result was not significant, likely reflecting the limited informative variation.  1 Number of polymorphic sites. 2 Number of haplotypes (standard deviation). 3 Haplotype diversity (standard error). 4 Nucleotide diversity across synonymous and non-synonymous sites (standard error). 5 Nucleotide diversity at non-synonymous sites (standard error). 6 Nucleotide diversity at synonymous sites (standard error). There was also no increase over time in the proportion of patients with the 2 most common mutations, R34L and T91P.

Spatial and temporal clustering analysis
Of patients with pkdhfr mutations, a non-randomly distributed spatial cluster was detected among a subset of those with the R34L mutation (n = 9 [25%]; RR 8.94; p = 0.00051) (Fig 1  and Table 5). In addition, all patients mapped with an E119V mutation were found in a single cluster (n = 5 [100%]; RR>100; p = 0.0017) geographically overlapping with that of the R34L cluster (Fig 2). None of the other pkdhfr mutations were detected in statistically significant spatial clusters, including the most common mutation, T91P (most likely cluster; RR 4.23; p = 0.296), or the pkdhfr wild-type as an independent group (most likely cluster; RR 1.51; p = 0.058). Cases with a pkdhfr E119V mutation were the only group which demonstrated statistically significant clustering on temporal analysis (p = 0.009), with no significant temporal association seen for the two most common pkdhfr mutations T91P and R34L, or the wild-type. Clusters in both the pkdhfr R34L and E119V cases also appeared non-randomly distributed in the combined space-time analysis, with p values of 0.0075 and <0.0001, respectively.

Epidemiology of pkdhfr E119V mutation cluster
Five patients with the E119V mutation were located within a 29.7 km radius, although all were from different villages. Median age was 49 years (IQR 25-29) and 3 were female ( Table 6). All five patients presented within a 37-day period from January to February 2014, with two females presenting on the same day. Two patients had severe malaria, while the median parasite count for all patients was 5,034/μL (IQR 774-5384). All patients reported recent forest and/or plantation exposure.

Epidemiology of patients within the pkdhfr R34L mutation cluster
Nine patients with the R34L mutation were located within a 27.4 km radius. Median age was 55 years (IQR 51-68 years), five were male, and all were of either the closely related Dusun or Kadazan ethnicity. One patient had severe disease, with parasitaemia of 693,061/μL. Five reported a previous history of malaria. None were found to be family members, and it is not known if any of these patients had been exposed to one another. Two patients presented to the hospital 34 days apart between June and July 2011, while another three presented within a 44-day period between March and April 2012. All but one reported recent plantation and/or forest exposure.

Effect of pkdhfr mutations on pyrimethamine binding
The 3D structures of 19 pkdhfr amino acid sequences derived from patients with a P. knowlesi monoinfection were constructed using the SWISS MODEL server. Autodock vina was used to assess the effect of mutations on pyrimethamine binding. This encompassed both pkdhfr wildtype and all single monomorphic mutations (with the exception of D157Y), and also all genotypes with multiple mutations (with the exception of R34L-D157N). The results of molecular docking analysis indicated that five residues are in direct contact with the inhibitor pyrimethamine: I13, D53, S120 and I173. None of the mutations detected in this study were involved in inhibitor interactions and therefore, had no modelled effect on pyrimethamine binding (Fig 3).

Discussion
This study demonstrated a high rate of pkdhfr mutations in human knowlesi malaria, including a number of common mutations such as located at R34L, which along with E119V were associated with clusters of cases at adjacent villages. Clinical correlates hypothesized to be more likely with P. knowlesi H-H transmitted infections such as a higher likelihood of hyperparasitaemia Dihdyrofolate-Reductase Mutations in Plasmodium knowlesi and severe disease were not related to pkdhfr mutations [39,40]. The non-synonomous pkdhfr mutations described in this study do not provide evidence of H-H transmission, with homology modelling demonstrating that they are very unlikely to affect pyrimethamine binding and hence, are not likely to be a result of selective antifolate drug pressure. Antifolate resistance can be selected readily in experimental malarias including simian P. knowlesi [67]. The utility of antifolate-resistance associated mutations as a marker of H-H transmission in this study requires both sufficient SP drug exposure in humans to select for resistant P. knowlesi parasite populations, and sustained transmission from humans to other human or monkey hosts. In the past P. knowlesi was often misdiagnosed as P. malariae [3,7], and therefore would have been nominally treated with chloroquine rather than SP according to longstanding Malaysian Ministry of Health guidelines. However, before the current pre-elimination malaria setting in Sabah where hospital admission is mandatory for all malaria cases, fever was frequently treated empirically with single dose SP by public health officials in the community for presumed non-severe malaria. The accumulation of dhfr mutations associated with SP-resistance in P. falciparum over 30 years of SP usage is well documented in Sabah [31,68,69]. With P. knowlesi endemic during the same period and commonly misdiagnosed as P. falciparum [3], the parasites infecting humans would have been exposed to some degree of within-host SP selection pressure [5,6]. For this to be selected at a population level would have required onward transmission.
One third of P. knowlesi isolates in this study had pkdhfr single point mutations, and 3% had double mutations. There were no triple mutants. The proportions of multiple mutants are less than predicted from a Poisson distribution arguing against stepwise selection. In contrast, in P.falciparum multiple pfdhfr mutations associated with antifolate resistance were found in this study, as previously described in Sabah, indicative of stepwise selection [31,68,69]. In P.falciparum, the primary pfdhfr mutation acquired in this progression is S108N, and in P. vivax, which is more closely related to P.knowlesi it is S117N. However, modelling showed residues involved in pyrimethamine binding in pkdhfr were located at I13, D53, S120, and I173, none of which were found in their mutant form in this study.
In this study, two P. knowlesi genotypes appeared to have non-random spatial distributions. This is consistent with a previous study demonstrating a higher level of genetic relatedness of P. knowlesi populations, as measured by the degree of multi-locus allelic frequencies (F ST indices) from samples collected at different sites within Borneo, with decreasing geographical distance [44]. H-H transmission in these clusters appears unlikely but cannot be definitively excluded. Although the timing of hospital presentation for subsets of individuals in each cluster was consistent with potential exposure to a preceding human case with patent parasitaemia, the presence or duration of subpatent parasitaemia prior to this and hence the infectious period is not fully known. Finer scale clustering with H-H transmission would also be expected, including family or household members from the same village without significant travel history. Nearly all individuals in identified clusters also had forest and/or monkey exposure, and other potentially relevant geo-spatial variables were not controlled for in this analysis, while the relatively small number of samples in this study also limits their interpretation. The use of P. knowlesi spatio-temporal clusters alone to evaluate H-H transmission is constrained by the inability to discern whether a human was bitten by a human-fed mosquito, or if a cluster of humans were bitten by mosquitoes which had fed on the same monkey or monkey group. Assessing the possibility of individual level P. knowlesi H-H transmission using mutations associated with drug selection pressure, or other potentially better suited methods such as microsatellite markers or mitochondrial sequences, is also complicated by drug-exposed parasite populations in humans conceivably being transferred back to monkey hosts [17]. In this context comparative analysis of whether the same P. knowlesi dhfr mutations are present in nearby macaque hosts would also be useful but was not logistically feasible in this study.
The low level of nucleotide diversity and absence of putatively neutral synonymous variation in the pfhdfr region in Sabah is consistent with purifying selection acting upon a gene whose product has an important functional role, allowing little room for variation with the exception of the known pyrimethamine resistance conferring variants observed. In marked contrast, extensive haplotype and nucleotide diversity was observed in the pkdhfr region, with significantly higher diversity in synonymous versus non-synonymous sites. These patterns infer that the variation observed in pkdhfr might reflect neutral, stochastic processes, particularly within the monkey hosts where polyclonal infections are usual and transmission is intense [44]. Tajima's D test also confirmed a significant excess of low frequency variants relative to expectation under neutral evolution consistent with other reports [44,70], potentially reflecting a recent population expansion, but not necessarily excluding purifying selection. A more comprehensive understanding of the genetic diversity across the P. knowlesi genome is needed for effective interpretation of the trends observed in pkdhfr. In addition, further analysis of the patterns of variation in dhfr in a broader range of human and simian Plasmodium species from settings with varying epidemiology and history of pyrimethamine use are required to elucidate the selective forces on dhfr in different Plasmodium spp.

Conclusion
Despite an increasing incidence of P. knowlesi malaria in Sabah, Malaysia, the non-synonomous pkdhr mutations reported in this study do not appear to have been selected by repeated pyrimethamine exposure associated with competent onward transmission from humans, and hence cannot confirm H-H transmission. Despite these findings, low or unsustained H-H transmission cannot be excluded. In addition to establishing the transmission competency of naturally acquired P. knowlesi from human hosts, monitoring for spatio-temporal clusters of identical multi-locus genotypes and more detailed genome-wide studies may provide useful information on the likelihood of H-H transmission occurring.