Within-host and Population Diversity
A total of 206 P. vivax isolates from two areas of the GMS were genotyped at 10 MS markers. Complete genotyping data for all 10 loci were obtained for 166 (80.5%) isolates. Samples with fewer than five MS markers genotyped were removed from analysis. The predominant allele of each locus in each sample was used for population genetic analysis. The 10 MS markers all were highly polymorphic with each one having 2 – 20 alleles when all parasite populations were considered (S1 Table). Allele frequencies varied among the four populations. In the CMB samples, MS2 was the most diverse in 2004 with 20 alleles, whereas MS7 was the most diverse in 2016 with 16 alleles, followed by MS9 with 15 alleles. In the 2012 samples in western Thailand, MS20 was the most diverse with 18 alleles, followed by MS2 with 16 alleles. The 2015 samples from three provinces of Thailand also had MS2 as the most diverse marker [13]. Therefore, MS2 appeared to be among the most diverse markers in malaria-endemic areas in the GMS. In comparison, the least diverse MS markers differed among the populations: MS6 in the 2004 CMB population with 6 alleles, MS20 in the 2016 CMB population with two alleles; MS6 and MS12 equally least diverse in 2012 western Thailand with 7 alleles; and MS1 in 2015 western Thailand with 6 alleles (S1 Table).
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, in western Thailand, 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 in CMB, and from 1.50 in 2012 to 1.36 in 2015 in western Thailand (Table 1). The genetic diversity was correspondingly reflected in allele richness in these parasite populations (Table 1). At the CMB, comparison of allele richness in parasite populations in 2004 and 2016 showed a substantial reduction from 12.2 to 8.8. In contrast, the allele richness in parasite populations in western Thailand did not show much changes between 2012 and 2015.
Despite the overall reduction in MOI, the genetic diversity of parasite populations at both sites remained high (Table 1). Haplotype diversity was high; a total of 162 haplotypes were observed and no haplotypes were shared among the four populations. The 2015 Thailand parasite population harbored 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 population in western Thailand had a slight increase from 2012 to 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 ~2-fold decrease over the past decade, whereas Ne in western Thailand remained high and even had a ~2-fold increase in recent years (Table 2).
Multilocus LD and Population Bottlenecks
Though haplotype sharing among the four parasite populations was not detected, the 2004 CMB and 2015 Thailand populations had low degrees of haplotype sharing within the populations. In the 2004 CMB population, two parasite isolates carried the same haplotypes 13 and 22, respectively (S2 Table). In the 2015 Thailand parasite population, identical haplotypes 25, 79, 92, and 109 were shared by two, two, three, and four parasites, respectively [13]. 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 western Thailand (P < 0.0001) (Table 3). Multilocus LD was also retained when only unique haplotypes or only one MS marker per chromosome (excluding MS2 and MS12) was used, suggesting the presence of significant inbreeding in all these populations. 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 in Tak Province of Thailand 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, we performed BOTTLENECK analysis under SMM and computed the Garza-Williamson index statistic (Table 4). The model detected significant deficiency in HE from the mutation-drift equilibrium, indicating events of population size reduction with possible clonal expansion in all four populations. Likewise, the mean G-W 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 western Thailand 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 western Thailand 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 Thailand population (FST = 0.064), but significant genetic differentiation from the 2015 Thailand population (FST = 0.172). The Mantel test on correlation between pairwise genetic and geographical distances indicated a weak, insignificant positive association between genetic and geographical distances (Mantel rank test, r2 = 0.0297, p = 0.44). Analysis of molecular variation (AMOVA) indicated 87% of the variation was attributed to variations within populations, and 13% among populations.
STRUCTURE analysis showed a clear distribution pattern of parasite haplotypes and demonstrated multiple sub-populations (Fig 1A). Most parasites from Thailand 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 (S2 Fig). At K = 3, populations collected in Thailand in 2015 formed a separate cluster, as well as the CMB samples collected in 2016, with admixture in samples from the CMB in 2004 and Thailand 2012. At K = 4, admixing between the CMB 2004 and Thailand 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, parasites in 2004 from the CMB and 2012 from Thailand were substantially admixed (Fig 1A), suggesting similar ancestry for parasites from the CMB and 2012 Thailand.
In line with the STRUCTURE analysis, PCA 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 Thailand 2015, Thailand 2012 and CMB 2016 formed individual clusters with no overlap, whereas the older parasite population CMB 2004 showed haplotype admixture with Thailand 2012 and CMB 2016 (Fig 1B). Such a finding was further illustrated with phylogenetic and network analyses.
The phylogenetic tree also supports the finding of substantial regional population structure and local clustering of haplotypes (Fig 2A). Clades 1 and 2 contained parasites predominantly from the CMB 2004 and 2016, which more or less reflected the lower right coordinate of the PCA analysis to show extensive clustering of parasites collected most recently from Myanmar with earlier parasite samples from the rest of the GMS (Fig 1B). Clade 3 and 4 parasites were predominantly from Thailand collected in 2012 and 2015, respectively. Similarly, network analysis roughly divided parasites into four clusters, with each cluster being dominated by parasites from a single population (Fig 2B). Cluster 4 contained mostly parasites from western Thailand 2015. Cluster 3 contained mostly parasites from western Thailand collected in 2012, while it also contained a few haplotypes from CMB border in 2004. Clusters 1 and 2 were dominated by haplotypes from CMB 2016 and 2004, respectively. Parasites collected most recently in 2016 from CMB were mainly found in clusters 1 and 2, where they showed close relationships with haplotypes from the earlier parasite population collected in 2004.