Ethics statement
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). 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 [35], 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 1,364 in 2016 [13]. In 2012 and 2016, P. vivax accounted for 57.5% and 86.7% of the total cases, respectively (http://203.157.41.215/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 [13]. 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.
Microsatellite genotyping
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 ml of elution buffer and used immediately or stored at -20°C until further use.
A volume of 2 ml 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 [39]. Ten MS markers (MS1, MS2, MS5, MS6, MS7, MS9, MS10, MS12, MS15, MS20) previously used to differentiate P. vivax populations in Thailand [13] 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 [40]. 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 [41]. Genotyping success was defined as the presence of at least one allele at a given locus in a given sample.
Data analysis
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 [42]. 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 [14]. 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 [43] and GenAlEx version 6.5 [44]. 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 [45] 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 [11]. 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 [46]. In addition, the mutation rate of P. vivax calibrated from an eradicated European strain (5.57×10-7), which was ~3,500× lower than that estimated for P. falciparum, was also used [47].
Bottleneck analysis
A heterozygosity excess test at the population level was used to detect the recent population bottleneck using BOTTLENECK 1.2.02 [48] with the SMM [49] and the two-phase model (TPM) [50]. 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 [51] 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) [52]. 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.
Population structure
Deviations from Hardy-Weinberg equilibrium could indicate the presence of population structure or inbreeding [53]. 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 [54]. 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 [57] and DISTRUCT 1.1 [58] 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 [59] implemented in MEGA7 [60]. Minimum spanning tree of parasite genotypes constructed by the goeBURST algorithm using the Phyloviz software v1.1 (http://www.phyloviz.net/)