Dynamics of Plasmodium vivax populations in border areas of the Greater Mekong sub-region during malaria elimination
Malaria Journal volume 19, Article number: 145 (2020)
Countries within the Greater Mekong Sub-region (GMS) of Southeast Asia have committed to eliminating malaria by 2030. Although the malaria situation has greatly improved, malaria transmission remains at international border regions. In some areas, Plasmodium vivax has become the predominant parasite. To gain a better understanding of transmission dynamics, knowledge on the changes of P. vivax populations after the scale-up of control interventions will guide more effective targeted control efforts.
This study investigated genetic diversity and population structures in 206 P. vivax clinical samples collected at two time points in two international border areas: the China-Myanmar border (CMB) (n = 50 in 2004 and n = 52 in 2016) and Thailand-Myanmar border (TMB) (n = 50 in 2012 and n = 54 in 2015). Parasites were genotyped using 10 microsatellite markers.
Despite intensified control efforts, genetic diversity remained high (HE = 0.66–0.86) and was not significantly different among the four populations (P > 0.05). Specifically, HE slightly decreased from 0.76 in 2004 to 0.66 in 2016 at the CMB and increased from 0.80 in 2012 to 0.86 in 2015 at the TMB. The proportions of polyclonal infections varied significantly among the four populations (P < 0.05), and showed substantial decreases from 48.0% in 2004 to 23.7 at the CMB and from 40.0% in 2012 to 30.7% in 2015 at the TMB, with corresponding decreases in the multiplicity of infection. Consistent with the continuous decline of malaria incidence in the GMS over time, there were also increases in multilocus linkage disequilibrium, suggesting more fragmented and increasingly inbred parasite populations. There were considerable genetic differentiation and sub-division among the four tested populations. Temporal genetic differentiation was observed at each site (FST = 0.081 at the CMB and FST = 0.133 at the TMB). Various degrees of clustering were evident between the older parasite samples collected in 2004 at the CMB and the 2016 CMB and 2012 TMB populations, suggesting some of these parasites had shared ancestry. In contrast, the 2015 TMB population was genetically distinctive, which may reflect a process of population replacement. Whereas the effective population size (Ne) at the CMB showed a decrease from 4979 in 2004 to 3052 in 2016 with the infinite allele model, the Ne at the TMB experienced an increase from 6289 to 10,259.
With enhanced control efforts on malaria, P. vivax at the TMB and CMB showed considerable spatial and temporal differentiation, but the presence of large P. vivax reservoirs still sustained genetic diversity and transmission. These findings provide new insights into P. vivax transmission dynamics and population structure in these border areas of the GMS. Coordinated and integrated control efforts on both sides of international borders are essential to reach the goal of regional malaria elimination.
As the global incidence of malaria has been greatly reduced in recent years, Plasmodium vivax has become the main source of malaria infections in many endemic areas where Plasmodium falciparum and P. vivax co-exist . This is also true for most areas of the Greater Mekong Sub-region (GMS) in Southeast Asia , which is actively pursuing regional malaria elimination. This shift in species dominance in the face of intensified control efforts highlights the remarkable adaptive potential and relative resilience of P. vivax to control measures . In recognizing the challenge for eliminating vivax malaria, the GMS countries planned to eliminate falciparum malaria by 2025 and all malaria by 2030 . As the malaria elimination plan progresses in the GMS, malaria displayed increasing heterogeneity in distribution, with transmission being concentrated along international borders . These international borders are porous with intensive human migration, which poses a major threat of parasite introduction [6,7,8]. In addition, border regions have drastically different ecology, vector systems, human populations, and are subject to influences of wars and civil unrest; thus malaria transmission in border regions are often unstable . While malaria surveillance has been strengthened in the GMS countries, it can be complemented by population genetic studies to define parasite diversity, population structure and migration among geographically separated transmission hotspots .
Since a major mechanism of generating genetic diversity in the malaria parasites is meiotic recombination of parasite strains in mosquito vectors, genetic diversity is intrinsically linked to the transmission intensity . Likewise, it is expected that reduction of parasite population in response to intensified interventions will lead to reduced parasite genetic diversity. However, P. vivax seems to defy this principle. In most of its geographic range, P. vivax showed high levels of genetic diversity as revealed by population genetic studies using various genotyping markers [12,13,14,15,16,17,18,19,20,21]. Paradoxically, even in many low transmission settings, P. vivax still maintains high levels of diversity [22,23,24,25,26]. As a result, in areas of P. vivax and P. falciparum co-endemicity, these two parasites tend to exhibit markedly different genetic diversity and population structures, with P. vivax showing more stable transmission patterns [14, 15, 27,28,29]. Newly introduced parasites in areas where malaria had been eliminated may initially display low-level genetic diversity or even clonality [30, 31], but the parasites could rapidly re-establish diverse populations . As a result, the shrinking P. vivax populations in areas with intensified control frequently showed relatively high genetic diversity, but substantial levels of multilocus linkage disequilibrium (LD) and population substructure [15, 17, 23, 33].
The scale-up of malaria control efforts in countries of the GMS has led to substantial changes in malaria epidemiology with a noticeable rising proportion of vivax malaria [5, 34]. In the China-Myanmar border (CMB) areas, P. vivax not only has become the predominant Plasmodium species in recent years , but also has caused several malaria outbreaks [9, 36]. Even in the international border regions, malaria transmission is concentrated in separated hotspots, and among which gene flow is expected to be low . Final attacks to eliminate these hotspots, guided by accurate surveillance, will ensure the success of the regional elimination programme. An improved understanding of the transmission dynamics and the adaptive responses of the parasites to the control measures will also be essential to guide and monitor the progress of the elimination campaign. Recent studies of P. vivax populations from international border regions of Thailand with microsatellite (MS) markers clearly detected parasite population division, consistent with the separation of these parasite populations by the malaria-free central region [13, 38]. Here P. vivax clinical samples collected at different time points of the same study sites were used to determine the effect of intensified control efforts on genetic diversity and population dynamics of P. vivax in two border regions of the GMS.
Written informed consent was obtained from all P. vivax patients or their legal guardians for participants under the age of 18 years. This study was approved by the Institutional Review Boards of Department of Disease Control, Thailand Ministry of Health and the Pennsylvania State University. Use of the de-identified blood samples was approved by the Ethics Review Committee of China Medical University.
Study sites and Plasmodium vivax clinical samples
In order to investigate the change of the P. vivax populations in response to intensified malaria control efforts in the GMS, P. vivax isolates were collected at two time points from two different areas of the CMB and Thailand-Myanmar border (TMB) (Additional file 1: Fig. S1). The CMB area has been traditionally malaria hyper-endemic with both P. falciparum and P. vivax transmission. Although recent malaria control efforts have sharply reduced the P. falciparum incidence , P. vivax incidence remained consistently high and even experienced outbreaks in recent years [9, 36]. At the CMB, 50 samples were collected in 2004 to represent the parasite population prior to the implementation of malaria elimination measures, while 52 samples collected in 2016 represented the parasites after scaling-up of malaria control. These samples were all collected from symptomatic P. vivax patients at the hospital and clinics in the adjacent Laiza Township (Myanmar) and Nabang Township (China). These hospital and clinics are located within 2 km from each other, and are thus considered as a single study site. In Thailand, Tak Province has been among the most malaria-endemic areas. In the past decade, malaria incidence in Tak has experienced more than 10-fold reduction, from 13,706 in 2012 to 1364 in 2016 . In 2012 and 2016, P. vivax accounted for 57.5% and 86.7% of the total cases, respectively (http://22.214.171.124/malariar10/index_newversion.php). Fifty P. vivax clinical samples obtained in 2012 from the Tha Song Yang hospital, Tak Province represented the parasites before the endorsement of malaria elimination plan, and were used to compare with samples collected from the same hospital in 2015 . While the CMB 2004 samples were collected as part of routine surveillance, the rest of the samples were collected from passive case detection implemented in these study sites for the International Center of Excellence for Malaria Research (ICEMR) project. Finger-prick blood samples were collected from symptomatic patients after obtaining informed consent. Malaria diagnosis was based on microscopy of Giemsa-stained thin and thick smears. For confirmed P. vivax cases, ~ 100 μl of finger-prick blood were spotted onto Whatman filter paper, dried and stored in individually zipped plastic bags.
Parasite genomic DNA was isolated from the dried blood spots on filter paper according to the protocol of QIAmp DNA Mini kit (Qiagen, Hilden, Germany). The final purified DNA was eluted into 35 μl of elution buffer and used immediately or stored at − 20 °C until further use.
A volume of 2 μl of purified DNA was used as the template for malaria parasite detection by a genus-specific and species-specific nested PCR targeting the Plasmodium 18S rRNA genes . Ten MS markers (MS1, MS2, MS5, MS6, MS7, MS9, MS10, MS12, MS15, MS20) previously used to differentiate P. vivax populations in Thailand  were used for genotyping these P. vivax samples. A multiplex primary PCR was done using all 10 primer pairs and followed by singleplex secondary PCR with a fluorescently labelled (6-FAM, VIC or NED) forward primer as previously described . PCR products were used for GeneScan fragment analysis on an ABI3730xl capillary electrophoresis platform (Applied Biosystems) using the size standard LIZ500. Genotype calling was facilitated with GeneMapper Version 4.0. The predominant allele and any additional alleles with peak height at least one-third of the height of the predominant allele per locus were scored . Genotyping success was defined as the presence of at least one allele at a given locus in a given sample.
Isolates containing more than one peak for any marker were considered to be multiple clone infections. The multiplicity of infection (MOI) was defined as the maximum number of alleles observed at any of the loci investigated. The mean MOI was calculated from the individual samples for each study site. Alleles were binned using the TANDEM software . Isolates with one allele at all markers, were considered monoclonal infections. An infection was defined as polyclonal if more than one allele was observed at one or more loci. For isolates with more than one allele at any of the loci, the alleles with the highest peak were used to construct the dominant haplotypes as previously described . Input files for the various population genetics software programs were created using CONVERT version 1.31.
Population diversity and differentiation
The indices of genetic diversity within populations, such as the number of polymorphic loci, the number of haplotypes (Nh), the number of alleles (Na), the mean allelic richness, and the expected heterozygosity (HE) were calculated using Arlequin version 3.11  and GenAlEx version 6.5 . In order to assess the genetic differentiation among populations, pairwise comparisons were measured by calculating FST using GenAlEx version 6.5. To estimate the partitioning of genetic variance for different hypothesized population groupings, analysis of molecular variance (AMOVA) was performed using GenAlEx version 6.5.
Linkage disequilibrium (LD)
Multilocus LD in each population was calculated using the program LIAN 3.7  with 50,000 iterations for burn-in and then 100,000 Markov Chain Monte Carlo (MCMC) iterations. Samples with missing data were excluded for this analysis. To avoid detecting false inbreeding resulting from clonal propagation and physical linkage, this analysis was performed for the combined dataset of single and dominant haplotypes and for unique haplotypes only. In addition, LD was analysed using only monoclonal and low-complexity (containing just one multi-allelic locus) samples. Since MS2 and MS5 both localize to chromosome 6 and MS12 and MS15 to chromosome 5, MS2 and MS12 were excluded in separate analysis as they had higher levels of missing data.
Effective population size (Ne)
The effective population size was calculated using the stepwise mutation model (SMM) and infinite alleles model (IAM) as previously described . Mutation rates for P. vivax are lacking and thus the P. falciparum mutation rate of 1.59 × 10−4 (95% confidence interval: 6.98 × 10−5–3.7 × 10−4) was used . In addition, the mutation rate of P. vivax calibrated from an eradicated European strain (5.57 × 10−7), which was ~ 3500× lower than that estimated for P. falciparum, was also used .
A heterozygosity excess test at the population level was used to detect the recent population bottleneck using BOTTLENECK 1.2.02  with the SMM  and the two-phase model (TPM) . The Wilcoxon signed rank test, a robust statistic for testing less than 20 polymorphic loci, was executed in the model in order to ascertain the probability of significant heterozygote excess. Since the method implemented in the bottleneck has low power  unless the decline is greater than 90%, the Garza-Williamson index was computed using Arlequin version 3.11. The Garza-Williamson index is the mean ratio of the number of alleles at a given locus to the range in allele size, i.e., M = (k/r), where k is the number of alleles and r is the allelic range (i.e., the difference in repeat units between the shortest and the longest alleles at a locus) . This measure is based on the assumption that in a bottleneck event, the number of alleles decreases faster than the allelic range because the latter is only reduced if the shortest and/or longest allele is lost, whereas the loss of any allele reduces the former. For the Garza-Williamson index, M < 0.68 indicates a bottleneck, whereas M > 0.80 indicates no reduction of effective population size.
Deviations from Hardy–Weinberg equilibrium could indicate the presence of population structure or inbreeding . To investigate whether haplotypes cluster into distinct genetic populations (K) among the defined geographic areas, a Bayesian analysis of population structure was conducted using STRUCTURE version 2.3.2 . The admixture model was used and the posterior probability of the grouping number (K = 1–10) was estimated by the MCMC method with 10 separate runs to evaluate the consistency of the results. Each run was estimated as 10,000 steps with 100,000-step burn-in. The best-fit number of grouping was evaluated using ΔK in the STRUCTURE HARVESTER version 0.6.93 tool [55, 56]. The CLUMPP version 1.1.2  and DISTRUCT 1.1  were used to display the partitioning clusters.
To visualize genetic relationships among the parasite isolates from the four populations, an individual-based principal coordinate analysis (PCoA) was conducted in GenAlEx version 6.5 using the genetic distances among MS genotypes. In addition, phylogenetic relationships amongst P. vivax isolates were analysed using the Neighbour-Joining method  implemented in MEGA7 . Minimum spanning tree of parasite genotypes constructed by the goeBURST algorithm using the Phyloviz software v1.1 (http://www.phyloviz.net/).
Within-host and population diversity
A total of 206 P. vivax isolates from two areas of the GMS were genotyped at 10 MS markers. Of the 10 MS markers, failure rates for MS7 (16.5%) and MS9 (15.0%) were the top two highest. Complete genotyping data at all 10 loci were obtained for 171 (83.0%) isolates. Since the maximum number of missing MS data for each sample was 4, all samples were included in the analysis. The 10 MS markers all were polymorphic with each one having 2-20 alleles when all parasite populations were considered and allele frequencies varied among the four populations (Additional file 2: Table S1).
Overall, 72 (34.9%) parasite isolates contained polyclonal infections (more than one peak for at least one MS marker), and the proportions of multiclonal infections were significantly different among the four parasite populations (P < 0.05, Pearson Chi square test) (Table 1). At both sites, along with the reduction in malaria incidence, there were substantial decreases in the proportion of multiclonal infections (P < 0.05, Pearson Chi square test). At the CMB, compared with the 48.0% polyclonal infections observed in 2004, this proportion decreased to 30.7% in 2016. Similarly, at the TMB, the percentage of polyclonal infections also experienced a considerable decrease from 40.0% in 2012 to 23.7% in 2016. Accordingly, this was reflected in the reduction of mean MOI from 1.48 in 2004 to 1.33 in 2016 at the CMB, and from 1.50 in 2012 to 1.36 in 2015 in the TMB area (Table 1).
Despite the overall reduction in MOI, the genetic diversity of parasite populations at both sites remained high (Table 1). First, allele richness in these parasite populations remained high (Table 1). At the CMB, allele richness in parasite populations showed a substantial reduction from 12.2 in 2004 to 8.8 in 2016. In contrast, the allele richness in parasite populations at the TMB did not show much change, at 11.2 in 2012 and 12.8 in 2015. Second, haplotype diversity was high; a total of 162 haplotypes were observed and no haplotypes were shared among the four populations. The 2004 CMB and 2015 TMB parasite populations harboured the greatest number of haplotypes (Table 1). A slight decrease in the expected heterozygosity was observed at the CMB from 0.76 in 2004 to 0.66 in 2016, although this did not reach statistical significance (P > 0.05). In comparison, the genetic diversity of parasite populations at the TMB had a slight increase from 0.80 in 2012 to 0.86 in 2015 (Table 1). This high genetic diversity may reflect substantially large parasite populations sustaining continued transmission at these two border areas. With both the SMM and IAM models, the effective population size Ne at the CMB was moderate and showed a ~ twofold decrease over the past decade, whereas Ne in western Thailand remained high and even had a ~ twofold increase in recent years (Table 2). These effective population sizes were much larger when the new P. vivax mutation rate was used (Additional file 3: Table S2).
Multilocus LD and population bottlenecks
The standardized index of association ISA used to measure the degree of inbreeding revealed significant multilocus LD in the P. vivax populations from the both CMB (P < 0.0001) and TMB (P < 0.0001) (Table 3). Multilocus LD was also retained when only unique haplotypes or monoclonal and low-complexity samples were considered or only one MS marker per chromosome (excluding MS2 and MS12) was analysed (Additional file 4: Table S3), confirming that LD was not the result of false reconstruction or physical linkage of MS markers. Also noteworthy is that the LD patterns barely changed when unique haplotypes were used for analysis, suggesting that there is no evidence of epidemic-like transmission patterns in any of the sites/time points. Considering temporal changes in multilocus LD, there was an increasing trend of LD in parasite populations at the CMB from 2004 (ISA = 0.0236) to 2016 (ISA = 0.0329) and at the TMB from 2012 (ISA = 0.0366) to 2015 (ISA = 0.0586) (Table 3), which may contribute to the potential effect of population reduction in both sites with the scale-up of malaria control measures.
To detect whether there were significant changes in the parasite population size, BOTTLENECK analysis was performed under SMM and TPM, and the Garza–Williamson index was computed (Table 4). The SMM detected significant deficiency in HE from the mutation-drift equilibrium, indicating events of population size reduction with possible clonal expansion in all four populations. The TPM only detected such events in the TMB populations. Likewise, the mean Garza-Williamson index was less than 0.68 in all the four populations, also suggesting a reduction of the effective population size.
Genetic differentiation and population structure
Genetic differentiation of the P. vivax populations in the GMS were evaluated by estimating the Wright’s fixation index FST. As expected, the contemporary P. vivax populations from the CMB and TMB showed considerable genetic differentiation (FST = 0.169 and 0.237). The CMB P. vivax populations showed moderate differentiation between 2004 and 2016 (FST = 0.081), whereas the TMB parasite populations had substantial differentiation between 2012 and 2015 (FST = 0.133, Table 5). Interestingly, the 2004 CMB parasite population also showed little genetic differentiation from those collected from the 2012 TMB population (FST = 0.064), but significant genetic differentiation from the 2015 TMB population (FST = 0.172). Overall, AMOVA indicated 87% of the variations was attributed to variations within populations, and 13% among populations. Specifically, at the CMB, 93% of the variations was attributed to variations within populations and 7% among populations, whereas at the TMB, 88% of the variation was attributed to variations within populations and 12% among populations.
STRUCTURE analysis showed a clear distribution pattern of parasite haplotypes and demonstrated multiple sub-populations (Fig. 1a). Most parasites from the TMB collected in 2015 were separated from all other populations. Using the delta K method, the most likely number of sub-populations was identified as 2 (Additional file 5: Fig. S2). At K = 2, populations collected from the TMB in 2015 formed a separate cluster from the rest of the populations. At K = 3, in addition to the TMB2015 cluster, the CMB 2016 samples also formed a cluster, with admixture in samples from the CMB 2004 and TMB 2012. At K = 4, admixing between the CMB 2004 and the TMB 2012 parasite populations became more evident (Fig. 1a), which corroborated the weak genetic differentiation between these two populations from FST analysis (FST = 0.064, Table 5). At K = 5, the CMB 2004 and TMB 2012 parasites were substantially admixed (Fig. 1a), suggesting similar ancestry for these two parasite populations.
In line with the STRUCTURE analysis, PCoA confirmed the genetic separation of the four parasite populations. PC1 and PC2 explained 7.37 and 5.16% of the variances, respectively, and grouped the parasites into two main clusters. The three more contemporary samples TMB 2012, TMB 2015, and CMB 2016 formed individual clusters with no overlap, whereas the older CMB 2004 parasite population showed haplotype admixture with the TMB 2012 and CMB 2016 (Fig. 1b). The phylogenetic tree and network analysis further corroborated the findings of substantial regional population structure and local clustering of haplotypes (Fig. 2).
With P. vivax becoming dominant in most parts of the GMS, a better understanding of transmission dynamics and population structure of this parasite is crucial to guide targeted control efforts and to achieve the goal of malaria elimination under the changing settings of malaria epidemiology . This work aimed to study the changes of parasite populations in the GMS during the elimination process, with samples collected in two international border areas before and after the scale-up of malaria control interventions. This study demonstrated that despite intensified control efforts in the GMS, there were only moderate levels of reduction in parasite genetic diversity and the proportion of polyclonal infections. Meanwhile the P. vivax parasite displayed significant LD and substantial genetic differentiation. The possible epidemiological processes driving these patterns of diversity and population structure are discussed below.
The GMS is malaria hypo-endemic with malaria transmission occurring along international borders. In both CMB and TMB, camps for refugees or internally displayed people have been established, where there is constant human population movement. Genetic diversity of P. vivax populations in the two border regions of the GMS was moderately high with the expected heterozygosity ranging from 0.66 to 0.86. Compared to previous studies which used a nearly identical set of MS markers, the HE values in the current study were remarkably close to those from low-to-moderate transmission areas in Asia, such as Vietnam (0.68), Thailand (0.76), central China (0.816), Laos (0.75), and India (0.72), as well as South America, such as Colombia (0.79) and Brazil (0.71) [22, 29, 30, 62, 63]. Despite strengthened malaria control efforts in both regions, HE only decreased slightly at the CMB and remained relatively unchanged in the TMB. This may reflect the differences in the time span of samples collected at the two sites: 12 years at the CMB versus 2–3 years at the TMB. An earlier study conducted using the CMB 2011–2013 samples also revealed stable genetic diversity (HE = 0.82) in these parasite populations . The high genetic diversity in low transmission areas is a common phenomenon that has been reported in numerous P. vivax-endemic areas [22,23,24, 28]. Reasons for this are not completely clear but are considered to be multifactorial . Of direct relevance to the border settings is the mobility of the camp and border populations along the borders, which could introduce new parasites and enhance genetic diversity of the parasite populations. The earlier study conducted at the CMB provided evidence on the existence of extensive gene flow among border communities and across international borders . Unlike that international borders are present as a potential barrier for the spread of P. falciparum strains in the GMS , landscape factors at the CMB do not appear to impede gene flow . Cross-border frequent human migration facilitates the importation of malaria into low-transmission areas and this represents a major challenge to malaria elimination in the GMS .
Regardless of the time span difference, both the CMB and TMB populations showed substantial reduction (17–18%) in the proportion of polyclonal infections and MOI. Yet, the proportion of polyclonal infections was unchanged at the CMB from 29% in 2011–2013  to 30.7% in 2016. Compared to the hypo-endemic settings in central China (5.9%) and Malaysia (17.6–29.2%) and hyper-endemic settings such as Papua New Guinea (63.8–88.3%), Cambodia (89.6%), Indonesia (23–79%), and Ethiopia (35%) [15, 22, 40, 67,68,69], the proportions of polyclonal infections (30.7% and 23.7%) in the CMB and TMB parasite samples from recent years remained moderately high. In addition, parasite populations in these regions remained relatively large, especially at the TMB, although the estimates from the new P. vivax mutation rate seemed not realistic. Multiclonal infections, together with the large population size, will increase the chance of recombination, resulting in the observed high genetic diversity. It is therefore important to further enhance local malaria control efforts to achieve substantial reduction of the parasite populations in these transmission hotspots.
Another noticeable change is the detection of significant multilocus LD in the parasite populations despite high genetic diversity. This again seems to be a common finding in many vivax low-endemicity settings such as Indonesia, Papua New Guinea, South Korea, India, Vietnam, Colombia, and Brazil [15, 32, 62, 64, 70, 71]. Significant LD against a background of high diversity may reflect the existence of multiple spatially clustered infections within a defined population, which might have arisen from rapid reduction in transmission and effective population size as malaria control interventions have intensified in the GMS. The overall trend of increasing LD in both sites and elsewhere in Thailand  is a demonstration of reduced overall parasite genetic complexity during malaria elimination.
Across different border regions of Thailand, a previous MS analysis detected substantial parasite population differentiation . This study also detected considerable divergence of the parasite populations across time and space. Phylogenetic and network analyses showed that most of the parasites from individual populations fell into distinct clades or clusters, indicating that the parasite populations in the GMS have become increasingly differentiated. The high malaria transmission areas along the international border regions of the GMS were largely connected a decade ago . Consistently, the older parasite population from the CMB 2004 showed remarkable degrees of clustering with the CMB 2016 and TMB 2012 samples, which may reflect shared ancestry of these parasites in the recent past. In contrast, the current malaria map illustrated the presence of separated pockets of transmission foci from intensified control efforts of the malaria elimination campaign, which limit gene flow within the GMS . The genetic differentiation between the contemporary parasite populations is an indication of independent evolution of these isolated parasite populations as also documented for the P. falciparum parasites . In addition, heterogeneity within vector populations could affect transmission  and population structure of the parasites . The notable shift in vector populations observed in western Thailand may also be responsible for shaping the parasite population structure . Moreover, the parasite populations studied here have undergone apparent population bottlenecks, suggesting population reduction, or limited vectors or human movement within the defined geographic areas. At the CMB, the epidemic malaria transmission documented in the camps for internally displaced people and surrounding villages in 2013 and the detected reduction of effective population size may result from expansion of certain parasite isolates such as those that are resistant to the frontline treatment [77, 78]. Interestingly, the TMB 2015 samples were quite genetically distinct, and even differed considerably from the 2012 parasite from the same region, which may suggest population replacement. This is plausible as the border parasite populations are subject to extensive parasite introduction , and introduced parasites are capable of re-establishing a more diverse population within a short period . In addition, new malaria interventions such as mass drug administration  may have accelerated the parasite population changes. The establishment of community-based malaria posts equipped with rapid diagnostic tests and antimalarial drug in the eastern Myanmar has led to a sharp decline of Myanmar malaria patients seeking medical service at the Thai side , which may also be related to the parasite population change.
Overall, this investigation demonstrated that the scale-up of malaria control interventions for malaria elimination in the GMS has resulted in both spatial and temporal P. vivax population differentiation. However, the persistently high genetic diversity, moderate levels of polyclonal infections, and considerably large effective population sizes demand concerted efforts of border nations to implement more effective measures targeting this resilient parasite. Especially in western Thailand, both parasite genetic diversity and effective population size have experienced slight increases in recent years. The large degree of population movement, rapid ecological changes and complex vector population dynamics in western Thailand are very conducive to sustained malaria transmission . With a large proportion of infection linked to cross-border travel, it will be important to strengthen malaria control measures on both sides of the border. At these low-endemicity settings, the identification of the source of residual infections and implementation of targeted control may be critical for the final phase of elimination [82, 83]. These findings provide novel insights for malaria surveillance and enhanced control efforts especially targeting the fragmented populations to reach the goal of malaria elimination by 2030.
Availability of data and materials
The datasets supporting the conclusions of this article are available in additional files.
Analysis of molecular variance
- F ST :
Wright’s fixation index
Greater Mekong sub-region
- H E :
Infinite alleles model
Standardized index of association
Multiplicity of infection
The number of alleles
- N e :
The effective population size
Number of haplotypes
Principal coordinate analysis
Stepwise mutation model
Cotter C, Sturrock HJ, Hsiang MS, Liu J, Phillips AA, Hwang J, et al. The changing epidemiology of malaria elimination: new strategies for new challenges. Lancet. 2013;382:900–11.
WHO. World Malaria Report 2019. Geneva, World Health Organization, 2019. https://www.hoint/publications-detail/world-malaria-report-2019.
Howes RE, Battle KE, Mendis KN, Smith DL, Cibulskis RE, Baird JK, et al. Global epidemiology of Plasmodium vivax. Am J Trop Med Hyg. 2016;95:15–34.
WHO. Eliminating malaria in the Greater Mekong Subregion: United to end a deadly disease. Geneva, World Health Organization, 2016. https://www.hoint/malaria/publications/atoz/eliminating-malaria-greater-mekong/en/.
Cui L, Cao Y, Kaewkungwal J, Khamsiriwatchara A, Lawpoolsri S, Soe TN, et al. Malaria elimination in the Greater Mekong Subregion: challenges and prospects. In Towards Malaria Elimination: A Leap Forward. Manguin S, Dev V, Eds. IntechOpen; 2018:179-200.
Parker DM, Carrara VI, Pukrittayakamee S, McGready R, Nosten FH. Malaria ecology along the Thailand-Myanmar border. Malar J. 2015;14:388.
Sriwichai P, Karl S, Samung Y, Kiattibutr K, Sirichaisinthop J, Mueller I, et al. Imported Plasmodium falciparum and locally transmitted Plasmodium vivax: cross-border malaria transmission scenario in northwestern Thailand. Malar J. 2017;16:258.
Zhou G, Sun L, Xia R, Duan Y, Xu J, Yang H, et al. Clinical malaria along the China-Myanmar border, Yunnan Province, China, January 2011-August 2012. Emerg Infect Dis. 2014;20:675–8.
Zhou G, Lo E, Zhong D, Wang X, Wang Y, Malla S, et al. Impact of interventions on malaria in internally displaced persons along the China-Myanmar border: 2011-2014. Malar J. 2016;15:471.
Arnott A, Barry AE, Reeder JC. Understanding the population genetics of Plasmodium vivax is essential for malaria control and elimination. Malar J. 2012;11:14.
Anderson TJ, Haubold B, Williams JT, Estrada-Franco JG, Richardson L, Mollinedo R, et al. Microsatellite markers reveal a spectrum of population structures in the malaria parasite Plasmodium falciparum. Mol Biol Evol. 2000;17:1467–82.
Koepfli C, Rodrigues PT, Antao T, Orjuela-Sanchez P, Van den Eede P, Gamboa D, et al. Plasmodium vivax diversity and population structure across four continents. PLoS Negl Trop Dis. 2015;9:e0003872.
Kittichai V, Koepfli C, Nguitragool W, Sattabongkot J, Cui L. Substantial population structure of Plasmodium vivax in Thailand facilitates identification of the sources of residual transmission. PLoS Negl Trop Dis. 2017;11:e0005930.
Jennison C, Arnott A, Tessier N, Tavul L, Koepfli C, Felger I, et al. Plasmodium vivax populations are more genetically diverse and less structured than sympatric Plasmodium falciparum populations. PLoS Negl Trop Dis. 2015;9:e0003634.
Noviyanti R, Coutrier F, Utami RA, Trimarsanto H, Tirta YK, Trianty L, et al. Contrasting transmission dynamics of co-endemic Plasmodium vivax and P. falciparum: Implications for malaria control and elimination. PLoS Negl Trop Dis. 2015;9:e0003739.
Lo E, Lam N, Hemming-Schroeder E, Nguyen J, Zhou G, Lee MC, et al. Frequent spread of Plasmodium vivax malaria maintains high genetic diversity at the Myanmar-China Border, without distance and landscape barriers. J Infect Dis. 2017;216:1254–63.
Waltmann A, Koepfli C, Tessier N, Karl S, Fola A, Darcy AW, et al. Increasingly inbred and fragmented populations of Plasmodium vivax associated with the eastward decline in malaria transmission across the Southwest Pacific. PLoS Negl Trop Dis. 2018;12:e0006146.
Rodrigues PT, Alves JM, Santamaria AM, Calzada JE, Xayavong M, Parise M, et al. Using mitochondrial genome sequences to track the origin of imported Plasmodium vivax infections diagnosed in the United States. Am J Trop Med Hyg. 2014;90:1102–8.
Hupalo DN, Luo Z, Melnikov A, Sutton PL, Rogov P, Escalante A, et al. Population genomics studies identify signatures of global dispersal and drug resistance in Plasmodium vivax. Nat Genet. 2016;48:953–8.
Pearson RD, Amato R, Auburn S, Miotto O, Almagro-Garcia J, Amaratunga C, et al. Genomic analysis of local variation and recent evolution in Plasmodium vivax. Nat Genet. 2016;48:959–64.
Baniecki ML, Faust AL, Schaffner SF, Park DJ, Galinsky K, Daniels RF, et al. Development of a single nucleotide polymorphism barcode to genotype Plasmodium vivax infections. PLoS Negl Trop Dis. 2015;9:e0003539.
Liu Y, Auburn S, Cao J, Trimarsanto H, Zhou H, Gray KA, et al. Genetic diversity and population structure of Plasmodium vivax in Central China. Malar J. 2014;13:262.
Gunawardena S, Ferreira MU, Kapilananda GM, Wirth DF, Karunaweera ND. The Sri Lankan paradox: high genetic diversity in Plasmodium vivax populations despite decreasing levels of malaria transmission. Parasitology. 2014;141:880–90.
Friedrich LR, Popovici J, Kim S, Dysoley L, Zimmerman PA, Menard D, et al. Complexity of Infection and Genetic Diversity in Cambodian Plasmodium vivax. PLoS Negl Trop Dis. 2016;10:e0004526.
Zhong D, Lo E, Wang X, Yewhalaw D, Zhou G, Atieli HE, et al. Multiplicity and molecular epidemiology of Plasmodium vivax and Plasmodium falciparum infections in East Africa. Malar J. 2018;17:185.
Pacheco MA, Schneider KA, Cespedes N, Herrera S, Arevalo-Herrera M, Escalante AA. Limited differentiation among Plasmodium vivax populations from the northwest and to the south Pacific Coast of Colombia: a malaria corridor? PLoS Negl Trop Dis. 2019;13:e0007310.
Gray KA, Dowd S, Bain L, Bobogare A, Wini L, Shanks GD, et al. Population genetics of Plasmodium falciparum and Plasmodium vivax and asymptomatic malaria in Temotu Province. Solomon Islands. Malar J. 2013;12:429.
Fola AA, Harrison GLA, Hazairin MH, Barnadas C, Hetzel MW, Iga J, et al. Higher complexity of infection and genetic diversity of Plasmodium vivax than Plasmodium falciparum across all malaria transmission zones of Papua New Guinea. Am J Trop Med Hyg. 2017;96:630–41.
Ferreira MU, Karunaweera ND, da Silva-Nunes M, da Silva NS, Wirth DF, Hartl DL. Population structure and transmission dynamics of Plasmodium vivax in rural Amazonia. J Infect Dis. 2007;195:1218–26.
Iwagami M, Fukumoto M, Hwang SY, Kim SH, Kho WG, Kano S. Population structure and transmission dynamics of Plasmodium vivax in the Republic of Korea based on microsatellite DNA analysis. PLoS Negl Trop Dis. 2012;6:e1592.
Putaporntip C, Miao J, Kuamsab N, Sattabongkot J, Sirichaisinthop J, Jongwutiwes S, et al. The Plasmodium vivax merozoite surface protein 3beta sequence reveals contrasting parasite populations in southern and northwestern Thailand. PLoS Negl Trop Dis. 2014;8:e3336.
Kim JY, Goo YK, Zo YG, Ji SY, Trimarsanto H, To S, et al. Further evidence of increasing diversity of Plasmodium vivax in the Republic of Korea in recent years. PLoS ONE. 2016;11:e0151514.
Delgado-Ratto C, Gamboa D, Soto-Calle VE, Van den Eede P, Torres E, Sanchez-Martinez L, et al. Population genetics of Plasmodium vivax in the Peruvian Amazon. PLoS Negl Trop Dis. 2016;10:e0004376.
Baird JK. Asia-Pacific malaria is singular, pervasive, diverse and invisible. Int J Parasitol. 2017;47:371–7.
Li N, Parker DM, Yang Z, Fan Q, Zhou G, Ai G, et al. Risk factors associated with slide positivity among febrile patients in a conflict zone of north-eastern Myanmar along the China-Myanmar border. Malar J. 2013;12:361.
Geng J, Malla P, Zhang J, Xu S, Li C, Zhao Y, et al. Increasing trends of malaria in a border area of the Greater Mekong Subregion. Malar J. 2019;18:309.
Zeng W, Bai Y, Wang M, Wang Z, Deng S, Ruan Y, et al. Significant divergence in sensitivity to antimalarial drugs between neighboring Plasmodium falciparum populations along the eastern border of Myanmar. Antimicrob Agents Chemother. 2017;61:e01689–716.
Congpuong K, Ubalee R. Population genetics of Plasmodium vivax in four high malaria endemic areas in Thailand. Korean J Parasitol. 2017;55:465–72.
Wampfler R, Mwingira F, Javati S, Robinson L, Betuela I, Siba P, et al. Strategies for detection of Plasmodium species gametocytes. PLoS ONE. 2013;8:e76316.
Koepfli C, Timinao L, Antao T, Barry AE, Siba P, Mueller I, et al. A large Plasmodium vivax reservoir and little population structure in the South Pacific. PLoS ONE. 2013;8:e66041.
Anderson TJ, Su XZ, Bockarie M, Lagog M, Day KP. Twelve microsatellite markers for characterization of Plasmodium falciparum from finger-prick blood samples. Parasitology. 1999;119:113–25.
Matschiner M, Salzburger W. TANDEM: integrating automated allele binning into genetics and genomics workflows. Bioinformatics. 2009;25:1982–3.
Excoffier L, Laval G, Schneider S. Arlequin (version 3.0): an integrated software package for population genetics data analysis. Evol Bioinform Online. 2007;1:47-50.
Peakall R, Smouse PE. GenAlEx 6.5: genetic analysis in Excel. Population genetic software for teaching and research–an update. Bioinformatics. 2012;28:2537-9.
Haubold B, Hudson RR. LIAN 3.0: detecting linkage disequilibrium in multilocus data. Linkage Analysis. Bioinformatics. 2000;16:847-8.
Anderson TJ, Su XZ, Roddam A, Day KP. Complex mutations in a high proportion of microsatellite loci from the protozoan parasite Plasmodium falciparum. Mol Ecol. 2000;9:1599–608.
van Dorp L, Gelabert P, Rieux A, de Manuel M, de-Dios T, Gopalakrishnan S, et al. Plasmodium vivax malaria viewed through the lens of an eradicated European strain. Mol Biol Evol. 2020;37:773-85.
Cornuet JM, Luikart G. Description and power analysis of two tests for detecting recent population bottlenecks from allele frequency data. Genetics. 1996;144:2001–14.
Shriver MD, Jin L, Chakraborty R, Boerwinkle E. VNTR allele frequency distributions under the stepwise mutation model: a computer simulation approach. Genetics. 1993;134:983–93.
Di Rienzo A, Peterson AC, Garza JC, Valdes AM, Slatkin M, Freimer NB. Mutational processes of simple-sequence repeat loci in human populations. Proc Natl Acad Sci USA. 1994;91:3166–70.
Williamson-Natesan EG. Comparison of methods for detecting bottlenecks from microsatellite loci. Conserv Genet. 2005;6:551–62.
Garza JC, Williamson EG. Detection of reduction in population size using data from microsatellite loci. Mol Ecol. 2001;10:305–18.
Hartl DL, Clark GC. Principles of Population Genetics. Sunderland: Sinauer Associates; 1997.
Pritchard JK, Stephens M, Donnelly P. Inference of population structure using multilocus genotype data. Genetics. 2000;155:945–59.
Earl DA, VonHoldt BM. STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conservation Genet Resour. 2012;4:359–61.
Evanno G, Regnaut S, Goudet J. Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Mol Ecol. 2005;14:2611–20.
Jakobsson M, Rosenberg NA. CLUMPP: a cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure. Bioinformatics. 2007;23:1801–6.
Rosenberg NA. DISTRUCT: a program for the graphical display of population structure. Mol Ecol Notes. 2004;4(1):137–8.
Saitou N, Nei M. The neighbor-joining method: a new method for reconstructing phylogenetic trees. Mol Biol Evol. 1987;4:406–25.
Kumar S, Stecher G, Tamura K. MEGA7: Molecular Evolutionary Genetics Analysis Version 7.0 for Bigger Datasets. Mol Biol Evol. 2016;33:1870-4.
Barry AE, Waltmann A, Koepfli C, Barnadas C, Mueller I. Uncovering the transmission dynamics of Plasmodium vivax using population genetics. Pathog Glob Health. 2015;109:142–52.
Hong NV, Delgado-Ratto C, Thanh PV, Van den Eede P, Guetens P, Binh NT, et al. Population genetics of Plasmodium vivax in four rural communities in Central Vietnam. PLoS Negl Trop Dis. 2016;10:e0004434.
Imwong M, Nair S, Pukrittayakamee S, Sudimack D, Williams JT, Mayxay M, et al. Contrasting genetic structure in Plasmodium vivax populations from Asia and South America. Int J Parasitol. 2007;37:1013–22.
Fola AA, Nate E, Abby Harrison GL, Barnadas C, Hetzel MW, Iga J, et al. Nationwide genetic surveillance of Plasmodium vivax in Papua New Guinea reveals heterogeneous transmission dynamics and routes of migration amongst subdivided populations. Infect Genet Evol. 2018;58:83–95.
Shetty AC, Jacob CG, Huang F, Li Y, Agrawal S, Saunders DL, et al. Genomic structure and diversity of Plasmodium falciparum in Southeast Asia reveal recent parasite migration patterns. Nat Commun. 2019;10:2665.
Jitthai N. Migration and malaria. Southeast Asian J Trop Med Public Health. 2013;44 Suppl 1:166-200; discussion 306-7.
Abdullah NR, Barber BE, William T, Norahmad NA, Satsu UR, Muniandy PK, et al. Plasmodium vivax population structure and transmission dynamics in Sabah Malaysia. PLoS ONE. 2013;8:e82553.
Orjuela-Sanchez P, Sa JM, Brandi MC, Rodrigues PT, Bastos MS, Amaratunga C, et al. Higher microsatellite diversity in Plasmodium vivax than in sympatric Plasmodium falciparum populations in Pursat. Western Cambodia. Exp Parasitol. 2013;134:318–26.
Getachew S, To S, Trimarsanto H, Thriemer K, Clark TG, Petros B, et al. Variation in complexity of infection and transmission stability between neighbouring populations of Plasmodium vivax in southern Ethiopia. PLoS ONE. 2015;10:e0140780.
Menegon M, Bardaji A, Martinez-Espinosa F, Botto-Menezes C, Ome-Kaius M, Mueller I, et al. Microsatellite genotyping of Plasmodium vivax isolates from pregnant women in four malaria endemic countries. PLoS ONE. 2016;11:e0152447.
Bahk YY, Kim J, Ahn SK, Na BK, Chai JY, Kim TS. Genetic diversity of Plasmodium vivax causing epidemic malaria in the Republic of Korea. Korean J Parasitol. 2018;56:545–52.
Socheat D, Denis MB, Fandeur T, Zhang Z, Yang H, Xu J, et al. Mekong malaria. II. Update of malaria, multi-drug resistance and economic development in the Mekong region of Southeast Asia. Southeast Asian J Trop Med Public Health. 2003;34 Suppl 4:1-102.
Hewitt S, Delacollette C, Chavez I. Malaria situation in the Greater Mekong Subregion. Southeast Asian J Trop Med Public Health. 2013;44 Suppl 1:46-72; discussion 306-7.
Fouet C, Kamdem C, Gamez S, White BJ. Genomic insights into adaptive divergence and speciation among malaria vectors of the Anopheles nili group. Evol Appl. 2017;10:897–906.
Joy DA, Gonzalez-Ceron L, Carlton JM, Gueye A, Fay M, McCutchan TF, et al. Local adaptation and vector-mediated population structure in Plasmodium vivax malaria. Mol Biol Evol. 2008;25:1245–52.
Sriwichai P, Samung Y, Sumruayphol S, Kiattibutr K, Kumpitak C, Payakkapol A, et al. Natural human Plasmodium infections in major Anopheles mosquitoes in western Thailand. Parasit Vectors. 2016;9:17.
Yuan L, Wang Y, Parker DM, Gupta B, Yang Z, Liu H, et al. Therapeutic responses of Plasmodium vivax malaria to chloroquine and primaquine treatment in northeastern Myanmar. Antimicrob Agents Chemother. 2015;59:1230–5.
Xu S, Zeng W, Ngassa Mbenda HG, Liu H, Chen X, Xiang Z, et al. Efficacy of directly-observed chloroquine-primaquine treatment for uncomplicated acute Plasmodium vivax malaria in northeast Myanmar: A prospective open-label efficacy trial. Travel Med Infect Dis. 2019:101499.
von Seidlein L, Peto TJ, Landier J, Nguyen TN, Tripura R, Phommasone K, et al. The impact of targeted malaria elimination with mass drug administrations on falciparum malaria in Southeast Asia: a cluster randomised trial. PLoS Med. 2019;16:e1002745.
Landier J, Parker DM, Thu AM, Lwin KM, Delmas G, Nosten FH. Effect of generalised access to early diagnosis and treatment and targeted mass drug administration on Plasmodium falciparum malaria in Eastern Myanmar: an observational study of a regional elimination programme. Lancet. 2018;391:1916–26.
Kar NP, Kumar A, Singh OP, Carlton JM, Nanda N. A review of malaria transmission dynamics in forest ecosystems. Parasit Vectors. 2014;7:265.
Ferreira MU, Rodrigues PT. Tracking malaria parasites in the eradication era. Trends Parasitol. 2014;30:465–6.
Edwards HM, Sriwichai P, Kirabittir K, Prachumsri J, Chavez IF, Hii J. Transmission risk beyond the village: entomological and human factors contributing to residual malaria transmission in an area approaching malaria elimination on the Thailand-Myanmar border. Malar J. 2019;18:221.
We want to thank the medical staff at the clinics for collecting the parasite samples and the patients for providing blood samples for this study.
This study was supported by grants from the National Institute of Allergy and Infectious Diseases (U19 AI089672), National Institutes of Health, USA.
Ethics approval and consent to participate
Ethical approvals were obtained from the institutional review boards of China Medical University, Thai Ministry of Health and Pennsylvania State University.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Map of GMS showing the two border regions (stars) where Plasmodium vivax clinical samples were collected.
The number of alleles (#) and allelic richness of each microsatellite marker in the four populations.
Effective population size (Ne) of the P. vivax populations estimated using the SMM and IAM models.
Multilocus linkage disequilibrium (ISA) in the P. vivax populations analysed using 8 microsatellite loci.
Estimation of the optimal number of populations (K) using the deltaK/K method.
About this article
Cite this article
Li, Y., Hu, Y., Zhao, Y. et al. Dynamics of Plasmodium vivax populations in border areas of the Greater Mekong sub-region during malaria elimination. Malar J 19, 145 (2020). https://doi.org/10.1186/s12936-020-03221-9
- Plasmodium vivax
- Population genetics
- Greater Mekong sub-region