Plasmodium vivax Population Structure and Transmission Dynamics in Sabah Malaysia

Despite significant progress in the control of malaria in Malaysia, the complex transmission dynamics of P. vivax continue to challenge national efforts to achieve elimination. To assess the impact of ongoing interventions on P. vivax transmission dynamics in Sabah, we genotyped 9 short tandem repeat markers in a total of 97 isolates (8 recurrences) from across Sabah, with a focus on two districts, Kota Marudu (KM, n = 24) and Kota Kinabalu (KK, n = 21), over a 2 year period. STRUCTURE analysis on the Sabah-wide dataset demonstrated multiple sub-populations. Significant differentiation (F ST  = 0.243) was observed between KM and KK, located just 130 Km apart. Consistent with low endemic transmission, infection complexity was modest in both KM (mean MOI  = 1.38) and KK (mean MOI  = 1.19). However, population diversity remained moderate (H E  = 0.583 in KM and H E  = 0.667 in KK). Temporal trends revealed clonal expansions reflecting epidemic transmission dynamics. The haplotypes of these isolates declined in frequency over time, but persisted at low frequency throughout the study duration. A diverse array of low frequency isolates were detected in both KM and KK, some likely reflecting remnants of previous expansions. In accordance with clonal expansions, high levels of Linkage Disequilibrium (I A S >0.5 [P<0.0001] in KK and KM) declined sharply when identical haplotypes were represented once (I A S  = 0.07 [P = 0.0076] in KM, and I A S = -0.003 [P = 0.606] in KK). All 8 recurrences, likely to be relapses, were homologous to the prior infection. These recurrences may promote the persistence of parasite lineages, sustaining local diversity. In summary, Sabah's shrinking P. vivax population appears to have rendered this low endemic setting vulnerable to epidemic expansions. Migration may play an important role in the introduction of new parasite strains leading to epidemic expansions, with important implications for malaria elimination.


Background
With an estimated 2.85 billion people living at risk of infection [1], Plasmodium vivax is the most widely distributed of the Plasmodium species infecting humans. Accumulating reports of resistance to chloroquine and associated severe and fatal disease complications [2] demonstrate the importance of addressing the particular challenges faced in controlling and eliminating this parasite. P. vivax differs from P. falciparum in having a dormant hypnozoite stage, greater asymptomatic asexual carriage and early gametocyte production, these factors enhancing the parasite's potential for local transmission. For this reason, the elimination of P. vivax faces far greater obstacles than that of P. falciparum [3,4]. With concerted efforts to eliminate P. falciparum and resultant decline in prevalence, P. vivax is becoming the dominant species of malaria in many coendemic regions.
Malaysia has demonstrated its commitment to the control and elimination of malaria, with national programs dating back to the 1960s. There has been a substantial decline in annual cases from an estimated 181,495 in 1967 to 55,000 in 1990 and 6,426 in 2010 [5]. Malaria control efforts have been particularly effective in Peninsular Malaysia, whilst the remote regions of Sabah and Sarawak, on the island of Borneo, have proven more challenging owing to constraints in accessibility to effective diagnostics and treatment. Malaysia is now committed to eliminating malaria by 2015 in the Peninsular region and 2020 on the island of Borneo [5].
Four species of malaria (P. falciparum, P. vivax, P. malariae and P. knowlesi) are endemic in Malaysia, with P. vivax contributing to more than 50% of the malaria cases nationwide in 2010 [5]. However, exact proportions remain uncertain owing to discrepancies between microscopy and PCR-based diagnoses [6,7]. Recent data indicates that the region most widely affected by malaria is Sabah, with an estimated 24.5% of the population living at risk of infection [8]. In 2011, P. vivax accounted for 31% of all malaria notifications [9]. Sabah is located in the north of Borneo, sharing borders with the Malaysian state Sarawak in the southwest, Kalimantan (Indonesia) in the south-east, and the islands of the Philippines to the north. P. vivax transmission persists in all of these border regions [5]. In this context, imported malaria is a major threat to elimination efforts, potentially undermining local intervention strategies, enhancing the parasite population diversity and adaptation potential, and increasing the risk of drug resistance spread and outbreaks in host populations with insufficient immunity.
An assessment of P. falciparum diversity in Sabah undertaken in 2005 demonstrated a fragmented population structure associated with the declining parasite endemicity [10]. However, the impact on the transmission dynamics and subsequent structure and diversity of the local P. vivax population is unknown.
The Asia Pacific Malaria Elimination Network (APMEN) was established in 2009 to address the specific challenges of malaria elimination in the region, including a high prevalence of P. vivax infections and the burden of imported malaria [11]. The APMEN Vivax Working Group (VxWG) is a body of 14 country partners that have acknowledged the utility of genotyping to inform on the parasite's complex transmission dynamics and patterns of spread within and across borders, developing a consensus with which to compare and contrast parasite population from a variety of locations. The protocol proposes nine short tandem repeat markers including eight putatively neutral microsatellites (MS1, MS5, MS8, MS10, MS12, MS16, MS20, Pv3.27) and one variable surface antigen marker (msp1F3) [12,13,14]. In order to assess the impact of current malaria intervention strategies on the local P. vivax transmission dynamics, P. vivax isolates from individuals residing in several regions of Sabah were genotyped using the consensus marker set. The resultant data were used to gauge the parasite diversity, population structure and transmission dynamics in the region.

Ethics
All samples were collected with written informed consent from the patient, parent or legal guardian (individuals ,18 years of age

Study Sites and sample collection
Sabah is divided into 5 administrative divisions (Kudat, Interior, Sandakan, Tawau and West Coast), with heterogeneity observed in the estimated prevalence of endemic P. vivax malaria [15]. Kudat Division spans an area of 4,623 square kilometres split between three districts, Kudat, Kota Marudu and Pitas. The West Coast Division covers an area of 7,588 square kilometres, divided into the districts of Ranau, Kota Belud, Tuaran, Penampang, Putatan, Papar, and the state capital Kota Kinabalu. The current study focused on Kota Kinabalu and Kota Marudu districts, where established collaborations with local health researchers at major hospitals and health centers enabled representative sampling from patients in the local districts. These sites, located approximately 130 Km apart, represent low (Kota Marudu; incidence = 1 case per 1000 population in 2010) and very low (Kota Kinabalu; incidence = 0.2 cases per 1000 population in 2010) P. vivax endemic settings ( Figure 1). Kota Marudu district spans a region of approximately 1,917 square kilometres, with an estimated population size of 80,900 in 2010 [16]. Kota Kinabalu district is home to the capital of Sabah, one of the main industrial and commercial centres in East Malaysia. The Malaysian Census 2010 Report estimated a population of 452,058 inhabitants in Kota Kinabalu [16]. Malaria is not markedly seasonal in these regions, which experience a typical equatorial climate with relatively constant temperature and high humidity year round.
With the exception of Tawau, P. vivax-infected individuals were recruited by passive case detection of patients presenting at Queen Elizabeth Hospital (QEH) [17], Kota Kinabalu, Kudat District Hospital, or Kota Marudu District Hospital between September 2010 and February 2013. Active case detection was undertaken in villages, palm plantation and logging camps in Tawau between June 2008 and November 2009. Capillary or venous blood samples were collected from patients with a positive rapid diagnostic test (RDT) (Paramax-3, Zephyr Biomedicals, India) during active case detection, or Plasmodium spp. positive thick blood smears on microscopy during passive case detection. Capillary sampling involved spotting ,200 ml blood onto an FTA card (Whatman), and storage as per the manufacturer's instructions. Venous blood sampling involved withdrawal of 2-5 ml blood into an EDTA-containing tube followed by transfer of a ,200 ml aliquot to an Eppendorf tube, and storage at 220uC or 280uC until DNA extraction.
Clinical management was based on microscopic diagnosis, in accordance with the national treatment policy for uncomplicated malaria: chloroquine (3 doses, 25 mg/kg body weight) and a 14 day course of unsupervised primaquine (0.5 mg/kg body weight daily) for patients with P. vivax monoinfection. Although patients were not prospectively followed after treatment, patients representing to the same clinic within the study period could be matched with their preceding history of malaria.

DNA extraction and PCR-based species identification
DNA extraction was undertaken on 200 ml samples using Qiagen's QIAamp DNA Blood mini kit following the manufacturer's instructions, with an elution volume of 200 ml. Plasmodium species was confirmed by PCR. Detection for P. vivax, P. falciparum, P. malaria and P. ovale parasites was undertaken using a version of Padley et al. [18] with the modification that each species was diagnosed in a separate (non-multiplex) assay. Detection for P. knowlesi parasites was undertaken using the method of Imwong et al. [19].

Genotyping
The MS1, MS5, MS8, MS10, MS12 and MS20 assays involved a single round of PCR using a slightly modified version of Karunaweera et al. [20]: 6.4 pmol of each primer, 0.1 mM each dNTP (Bioline), 16PCR buffer, 2.5 mM MgCl 2 , 1 Unit Taq DNA polymerase (Qiagen) and 2 ml genomic DNA in a final reaction volume of 20 ml. Thermocycling involved an initial denaturation at 94uC for 2 min, followed by 40 cycles of denaturation at 94uC for 30 sec, annealing at 58uC for 40 sec, elongation at 72uC for 30 sec, followed by a final elongation at 72uC for 5 min.
The MS16, Pv3.27 and msp1F3 assays involved nested and semi-nested PCRs as described by Koepfli et al. [14]. PCRs were performed in 20 ml reactions containing 0.25 mM of each primer, 2 ml of 10 x buffer, 0.2 mM each dNTP (Bioline), 2 mM MgCl 2 , 1.5 Unit Taq DNA polymerase (Qiagen) and 1 ml genomic DNA. Thermocycling conditions included an initial denaturation at 95uC for 5 min, followed by 30 cycles of denaturation at 95uC for 1 min, annealing (56uC for MS16; 58uC for Pv3.27; 59uC for msp1F3) for 1 min, elongation at 72uC for 1 min, followed by a final elongation at 72uC for 5 min. PCR products were diluted as necessary, and 1 ml of diluted product was used as template for the nested (MS16, msp1F3) or semi-nested (Pv3.27) assay, in a 20 ml volume with reagent concentrations as for the primary reaction.
Thermocycling conditions for the nest reaction included an initial denaturation at 95uC for 5 min, followed by 25 cycles of denaturation at 95uC for 1 min, annealing (57uC for MS16; 58uC for Pv3.27; 60uC for msp1F3) for 1 min, elongation at 72uC for 1 min, followed by a final elongation at 72uC for 5 min. The final PCR products were diluted as necessary and sized by denaturing capillary electrophoresis on an ABI 3130xl Genetic Analyzer (Applied Biosystems, Foster City, CA) with GeneScan LIZ-600 (Applied Biosystems) internal size standards. Genotype calling was facilitated with GeneMapper version 4.0. software (Applied Biosystems). A minimum arbitrary fluorescence intensity of 100 fluorescence units was set to exclude potential false peaks arising from background 'noise'. All electropherogram traces were additionally inspected manually. For each isolate, at each locus, the predominant allele (maximum height peak), and any additional alleles with peak height $33% the height of the predominant allele were scored. The implementation of a peak intensity threshold facilitates comparability between samples and populations by reducing potential bias in polyclonality estimates arising from differences in sample DNA quantity [21]. Genotyping success was defined as the presence of at least one allele at a given locus in a given sample.

Population Genetic Analysis
An infection was defined as polyclonal if two or more alleles were observed at one or more of the 9 loci investigated. The Multiplicity of Infection (MOI) for a given sample was defined as the maximum number of alleles observed at any of the 9 loci investigated. The average MOI was calculated by dividing the sum of MOIs across all samples by the number of samples.
'Infection haplotypes' were reconstructed from the predominant allele at each locus, as in other studies [21,22,23]. Population-level genetic diversity was characterised using a measure of the expected heterozygosity (H E ). Only the predominant allele at each locus in each isolate was used to ensure an unbiased estimate of the minor allele frequency [21]. H E provides a measure of the probability that two unrelated parasites will exhibit different genotypes at a given locus. H E was calculated for each locus in each population using the following formula: where n is the number of isolates analyzed and pi is the frequency of the ith allele in the population. The inclusion of 'n/(n21)' to adjust for sample size facilitated cross-population comparability. H E was averaged across the 9 loci in each population to provide a measure of population diversity.
Using infection haplotypes, multi-locus linkage disequilibrium (LD) was assessed by the standardised index of association (I A S ) using the web-based LIAN 3.5 software [24]. Testing the null hypothesis of linkage equilibrium, the significance of the I A S estimates was assessed using 10,000 random permutations of the data. Only isolates with no missing data at any of the loci investigated could be included in this analysis.
The differentiation between pairs of populations was measured using the pairwise F ST metric in the FSTAT software version 2.9.3, with Bonferroni correction for multiple tests, and a significance level of 0.05. Pairwise tests were undertaken for 1) Kota Kinabalu district versus Kota Marudu district, and 2) Kudat Division versus West Coast Division. Only the predominant allele at each locus in each isolate was used.
STRUCTURE software version 2.3.3 [25] was used to determine the most likely number of populations (K) in the total sample, and derive the probability of ancestry of each isolate to each of the K populations. The software required infection haplotypes and tolerated missing data. For optimal performance, samples missing data at 3 or more loci were excluded. Model parameters were admixture with correlated allele frequencies.
Twenty replicates, with 100,000 burnins and 100,000 iterations were run for each of K from 1-10. The most probable K was derived by calculating DK as described elsewhere [26] for each of K = 2-9. STRUCTURE results were displayed in barplots using distruct software version 1.1 [27].
The genetic relatedness between sample pairs was assessed by a simple measure of the proportion of alleles shared between haplotype pairs (ps), and using (1-ps) as a measure of genetic distance [28]. Samples with missing data at one or more loci were excluded from analysis. An unrooted neighbour-joining tree [29] was generated from the distance matrix using the APE (Analysis of Phylogenetics and Evolution) package in R [30,31].
With the exception of the neighbour-joining tree, recurrent infections (potentially non-independent) were excluded from the population genetic analyses.
Analysis of LD, population differentiation, structure, and genetic relatedness, were all undertaken with and without the msp1F3 locus. Exclusion of the locus had negligible impact on the results (same as test results excluding other, putatively neutral, loci) and, thus, the results are presented with inclusion of the msp1F3 data.
Temporal analysis of haplotype persistence was conducted on samples collected over a 29 months period from September 2011 to January 2013, in which all non-recurrent (i.e. independent isolates) were included with the exception of those from Tawau, which were detected by active surveillance over a shorter period of time.

Samples and Genotyping
A total of 97 P. vivax samples collected between September 2010 and January 2013 were included in the study (summary in Table 1), largely representing individuals residing in Kudat Division (n = 40, 2 in recurrent infections) or West Coast Division (n = 40, 6 in recurrent infections). Additional P. vivax samples collected from the Tawau (n = 12), Sandakan (n = 3) and Interior (n = 2) divisions were too small to gauge local diversity patterns and were thus included only in the Sabah-wide measures. Overall 98% (93/97) of the isolates exhibited successful genotype calls at 7-9

Within-host and Population Diversity
Across Sabah, the prevalence of polyclonal infections was 25.8% (23/89; Table 2) and a large proportion of these infections (60.9% [14/23]) exhibited multiple alleles at just one of the 9 loci investigated, indicative of highly related clones ( Figure S1). The prevalence of polyclonal infections in each of the divisions and districts remained below 30% (Table 2). On the basis of the apparent diversity of each of the 9 markers, the majority of polyclonal infections (.90%) would have been detected with the 6 most diverse markers, and 100% of polyclonal infections would have been detected with the 8 most diverse markers ( Figure S2). In accordance with the modest prevalence of polyclonal infections, infection complexity, as measured by the MOI, was modest in the Sabah-wide sample set (mean MOI = 1.30), as well as at the division (mean MOI = 1.37 in Kudat and 1.21 in the West Coast Division), and district levels (mean MOI = 1.38 in Kota Marudu and 1.19 in Kota Kinabalu) (

Population Structure and Differentiation
Pair-wise measures of the fixation index revealed significant differentiation between Kudat and the West Coast Division (F ST = 0.125), this differentiation being more apparent at the district level between Kota Marudu and Kota Kinabalu (F ST = 0.243) ( Table 3).
STRUCTURE analysis on the Sabah-wide dataset demonstrated multiple sub-populations. The most likely number of subpopulations was estimated to be four ( Figure S3). None of these sub-populations were observed exclusively in a given district, but the representation of three of the sub-populations (K1, K2 and K3) appeared to be moderately disproportionate between certain districts ( Figure 2 A neighbour-joining tree based on a pairwise distance matrix at 6 loci (MS10, MS8 and MS20 excluded) was generated from all samples with no missing data at these loci (n = 92) (Figure 3). The subset of 6 loci was selected to optimize sample size and geographic spread as well as marker number as the distance matrix analysis does not permit missing data. The impact of different marker subsets was investigated, and no notable alterations to the general clustering pattern of the tree were observed (data not presented). The neighbour-joining tree enabled further visual assessment of the population structure and genetic relatedness amongst the Sabahan samples. A large cluster of isolates, predominantly from Kota Kinabalu, which appeared to reflect the K3 sub-population, demonstrated identical or nearidentical haplotypes. The K1 sub-population appeared to be split into one large cluster of identical isolates, almost exclusively from Kota Marudu, and a second, cluster of identical or near-identical isolates derived from a range of districts spanning the Kudat and West Coast Divisions. The two clusters were not highly related to one another but possibly diverged from a recent common ancestor, with a recent recombination event breaking down the proportion of shared alleles. Aside from the large clonal expansions of K1 and K3, throughout the sample set, smaller groups of isolates with identical or near-identical haplotypes were observed, interspersed with a range of less closely related isolates reflecting the K4 sub-population.

Linkage Disequilibrium
When all samples were included in the analysis, high levels of linkage disequilibrium (LD) were observed Sabah-wide and at the division and district levels (I A S range 0.320 to 0.690; all P,0.0001) ( Table 4). However, when each unique haplotype was represented only once, marked reductions in LD were observed in all of the datasets (I A S range -0.003 to 0.109), with loss of significance in the West Coast division (P = 0.283) and Kota Kinabalu (P = 0.606).

Temporal Haplotype Dynamics
A total of 34 distinct haplotypes were observed in 72 samples collected between September 2011 and January 2013. A high frequency of these haplotypes (29.4% [10/34]) were observed more than once in independent samples, the frequency of multiple  incidences ranging from 2-13 ( Figure 4). The most frequently observed haplotypes were haplotype 15 (13 incidences The other multiple incidence haplotypes observed in the duration of the study generally presented at more distal time points, all observed at frequencies of 3 or less within a given 3 month period. Within the duration of the study period assessed (September 2011 to January 2013), the longest detected persistence of the same haplotype (haplotype 9) was 796 days. As observed with several other haplotypes, haplotype 9 was detected in a range of geographic locations, appearing first in Tuaran (West Coast Division), followed by Kota Kinabalu, Kota Marudu, and then Pitas.

Recurrent Infections
During the study period, 7 patients presented with recurrent P. vivax infections (one patient twice, all others once) between 28 and 104 days following primary therapy (Table 5). P. vivax parasitaemia was not detected by microscopy at 3 of the 8 'primary' presentations and thus primaquine (PQ) radical therapy was not prescribed for these cases. All sample pairs demonstrated homologous genotype calls at all of the successfully genotyped loci. The recurrences were sourced from individuals residing in a range of districts, and exhibited a range of haplotypes, including three haplotypes (9, 15 and 22) that had been observed in more than one independent (i.e. non-recurrent) infection (Figure 4).

Discussion
Effective reduction in local malaria transmission can lead to complex parasite transmission dynamics, with small populations being more effected by the impact of migration and random genetic drift [32]. The consequences of these processes, including the risks of epidemic expansions, and enhanced vulnerability to the emergence of drug resistant strains have clear implications for malaria elimination [33]. We investigated the population structure and transmission dynamics in P. vivax isolates sourced largely from the Kudat and West Coast divisions of Sabah. These two divisions account for approximately 36% of Sabah's population. As malaria elimination approaches, our results demonstrate that the low level of local transmission is associated with frequent clonal expansions, possibly resulting from the presentation of new parasite strains introduced by migration or generated locally by mutation or recombination. Although the incidence of these strains appeared to decline over time, several strains were detected in the population more than a year after the initial expansion, possibly reflecting long-term persistence. The available evidence indicated that, amongst other factors, homologous infections which we hypothesise to be relapses might facilitate the persistence of these strains in the local population.
Inevitable sample size constraints in pre-elimination settings present significant analytical challenges. In 2011, only 628 notifications of P. vivax were recorded across Sabah, with 151 (24%) of these reported from the Kudat and West Coast Divisions [9]. Nonetheless, for the purpose of population genetic characterisation, we were able to capture sufficient sample numbers from two districts, Kota Marudu and Kota Kinabalu, together representing 10% of the reported vivax cases in Sabah in 2011 [9]. Additional samples captured from other districts enabled a rough assessment of the relatedness between P. vivax isolates across the state.
Our analyses suggested that for most intents, Sabah should not be considered a single population. Rather, as seen in P. falciparum, the P. vivax population demonstrated multiple sub-populations [10]. STRUCTURE analysis revealed evidence of four subpopulations amongst the sites investigated. Although none of these sub-populations were exclusive to a given district, Kota Marudu and Kota Kinabalu exhibited moderately different STRUCTURE profiles, largely reflecting different representation of two subpopulations (K1 and K3). The pair-wise fixation index between Kota Marudu and Kota Kinabalu, located just 130 Km apart, demonstrated significant differentiation (F ST = 0.243), approach-   ing levels observed in international comparisons [12]. Although these samples represent a large proportion of the P. vivax isolates presenting in the Kudat and West Coast Divisions (based on the 2011 notifications described above), we should nonetheless acknowledge that sample size may influence the detection of population structure. Despite differences in P. vivax incidence, similar patterns of diversity were observed in Kota Marudu and Kota Kinabalu, suggesting that a common mechanism other than local transmission had a strong impact in both settings. In both districts, and across Sabah as a whole, rates of polyclonal infection were modest (less than 30%), and a large proportion of these polyclonal infections only displayed multiple alleles at one of the nine loci investigated (,60%). By genotyping the 6 most diverse (highest H E ) markers, we demonstrated that more than 90% of the polyclonal infections identified with the full set of 9 markers would have been detected. Therefore, we anticipate that genotyping at additional loci would not lead to a major increase in detection of polyclonal infections. Moreover, if extra markers are required to identify a polyclonal infection, by definition, it remains that the clones within that infection are highly related (i.e. overall infection diversity is modest). In accordance with the modest rates of polyclonal infection, the mean MOI was modest in all sites tested (less than 1.5 clones per infection). These measures of within-host diversity are generally lower than those reported for other tropical endemic sites in the Asia-Pacific region [12,20,23,34,35], and largely indicative of low levels of local transmission, as has been observed in other locations after aggressive malaria control interventions [36].
Despite the modest rates of polyclonal infection, population level diversity remained moderate in both Kota Marudu and Kota Kinabalu, which may reflect frequent introductions of new parasite strains. With small populations generally being more affected by the impact of migration [32], this may present an important reservoir for the introduction of new parasite strains, which may then result in small outbreaks. Neighbour-joining analysis illustrated large clusters of identical or nearly identical isolates in each of Kota Marudu and Kota Kinabalu. When these isolates were represented just once in linkage disequilibrium (LD) analysis, the high levels of LD observed in the presence of all isolates (I A S .0.5) declined sharply (I A S #0.07). The short temporal span within which peak frequencies of identical strains were observed, and the moderate geographic confinement of these incidences, further supports outbreak dynamics. Sabah's rapidly shrinking P. vivax population may thus have significant adverse public health consequences, including a greater impact of migration and increased risk of epidemics. In this low endemic Haplotypes with a minimum of 8 identical alleles at the 9 loci investigated and belonging to the same STRUCTURE sub-population (at K = 4) were grouped together as a single haplotype [22]. Singleton haplotypes were grouped together for simple visual representation in the plot (black bars). doi:10.1371/journal.pone.0082553.g004 setting, low levels of acquired immunity in the host population may have compounded the risk of epidemics. Indeed, in contrast to high transmission settings such as the island of New Guinea, where the majority of patients presenting with symptomatic P. vivax infection are young children [37], in Sabah, the majority of symptomatic patients were adults, with a median age of 24 years (Table 1). Another public health risk is presented by the persistence of specific strains in the population, with concurrent risks of parasite spread and maintenance of population diversity. In a slightly higher transmission setting in Brazil, only short short-term persistence of parasite lineages (up to 5 months) was observed [22]. In contrast, in Sabah, persistence of certain haplotypes was observed for up to two years, and this may be an underestimation limited by the duration of the study. As suggested by Ferreira and colleagues, variation in the persistence of parasite lineages might reflect a number of factors including differences in parasite mutation, recombination, and migration rates, population size, and acquired variant-specific immunity in the host population [22]. Homologous relapses might support the persistence of certain parasite lineages. Eight pairs of primary and recurrent infections were captured during the study, all representing as homologous haplotypes. The majority of these cases are likely to be relapses. Firstly, the timing of the recurrences (28 to 104 days after primary infection), were in accordance with the expected timing of relapses in tropical parasite strains [38]. Secondly, low local transmission levels in the regions concerned ensured a low probability of reinfection, although the potential impact of transmission hot-spots cannot be excluded [39]. Lastly, in at least 3 cases (where the primary infection was misdiagnosed as a species other than P. vivax [7]), adequate anti-relapse therapy was not administered thus increasing the risk of relapse which occurs naturally in this equatorial region in almost 80% of cases [38].
We observed incidences of the same haplotype identified in different districts and divisions at different time points, compatible with a process mediated by relapsing infections. This offers some promise for the ability to trace at least the short-term movements of particular parasite strains, with potential utility for transmission mapping efforts. However, better approaches for identifying recent mutations and recombinants are required. High-throughput genome sequencing efforts should facilitate this objective [40,41].
In the current study, interspersed amongst the clusters of identical haplotypes, we observed a number of singletons, many exhibiting very different haplotypes from one another. Some of these isolates might reflect remnants of previous expansions following the introduction of new strains which have persisted in the population. The diversity of these isolates largely explains the population diversity observed in Kota Marudu (H E = 0.583) and Kota Kinabalu (H E = 0.667). It follows that, despite lower P. vivax incidence, as a major industrial, administrative and commercial center, the higher level of P. vivax diversity observed in Kota Kinabalu might reflect more frequent immigration of novel genotypes. Further transmission mapping efforts over longer periods of time should enable greater insight into the role of migration in maintaining local P. vivax diversity.
In conclusion, our study highlights the potential utility of genotyping as a tool to inform P. vivax transmission dynamics, population diversity and structure. A better understanding of these dynamics in pre-elimination settings will inform on the underlying processes that sustain the dynamic transmission of parasite populations, thus assisting malaria control programs to prioritise activities and mobilise appropriate resources.