This study was conducted according to the principles expressed in the Declaration of Helsinki. Before blood collection, the study protocol, potential risks, and benefits were explained to the participants, and written informed consent was obtained from all adult participants and from the parents or legal guardians of children. Blood was collected following institutional ethical guidelines reviewed and approved by the Ethics Committee at the National Institute of Parasitic Diseases, Chinese Center for Disease Control and Prevention (no. 20120826).

A total of 180 blood samples of malaria patients infected with P. vivax were collected from the CMB area in YP (China) during 2009–2011 (n = 77) and 2014–2016 (n = 103) ( Figure 1 ). All samples were stained by Giemsa and verified by microscopic examination then confirmed for single infection with P. vivax by nested PCR (Zhou et al., 2014). Genomic DNA was extracted from whole blood using the DNeasy Blood & Tissue Kit (Qiagen, Germany) as previously reported (Chen et al., 2017; Kassegne et al., 2020).

To identify Plasmodium species, nested PCR was performed as previously described (Zhou et al., 2014). In the first PCR, the DNA fragment was amplified by 2×Taq PCR MasterMix (Tiangen Biotech, Beijing, China) using the primers rPLU1 5′-TCAAAGATTAAGCCATGCAAGTGA-3′ and rPLU5 5′-CCTGTTGTTGCCTTAAACTTC-3′. One microliter of template DNA was added to a 20-μl PCR mixture consisting of 0.4 μM of each primer, 10 μl 2×Taq PCR MasterMix (Tiangen Biotech, Beijing, China) containing 0.1 U Taq polymerase/μl, 500 μM deoxynucleotide triphosphates (dNTP), 3 mM MgCl₂, 100 mM KCl, and 20 mM Tris–HCl, pH 8.3. The cycling parameters to amplify the fragments were as follows: initial denaturation at 94°C for 5 min; 30 cycles of 94°C for 30 s, 55°C for 30 s, and 72°C for 1 min; followed by a final extension at 72°C for 5 min. One microliter of the first PCR product was used in the second amplification. Conditions and concentrations used for the second amplification were identical to those used for the first, except that rVIV1/rVIV2 (rVIV1 5′-CGCTTCTAGCTTAATCCACATAACTGATAC-3′, rVIV2 5′-ACTTCCAAGCCGAAGCAAAGAAAGTCCTTA-3′) were used as primers and amplification was performed over 35 cycles. The size of the DNA target amplified by these outer primers is about 1,600–1,700 bp and that by the inner primers is 121 bp. PCR products were qualitatively analyzed on 2% agarose gel and were sent to Beijing Genomics Institution (BGI, Shenzhen, China) for sequencing.

The PvAMA-1 fragment was amplified by PrimerSTAR Max DNA Polymerase (Takara, Japan) using the primers SeqF1 5′-CCCTACCAGCGGCTACTTC-3′ and SeqR1 5′-CGTTTGCTTGGCCAACTC-3′. All PCR amplifications were performed in a 50-μl PCR reaction volume containing 0.4 μM of each primer pair, 1.5 mM MgCl₂, 1× PCR buffer (50 mM KCl, 10 mM Tris–Cl, pH 8.3), 0.2 mM dNTPs, and 0.5 unit of DNA polymerase. The cycling parameters to amplify the fragments were as follows: initial denaturation at 98°C for 5 min, 35 cycles of denaturation at 98°C for 10 s, annealing at 69°C for 15 s, extension at 68°C for 1 min 30 s, and a final extension at 68°C for 5 min. The PCR fragments were about 1,900 bp in size and contained all the gene fragments of PvAMA-1 (1,689 bp). The PCR products were qualitatively analyzed on 1% agarose gel and were sent to Beijing Genomics Institution (BGI, Shenzhen, China) for sequencing. All unique mutations were carefully checked, and ambiguous bases were confirmed by resequencing. Meanwhile, a set of 73 sequences of Pvama-1 (1,290 bp) isolates from Myanmar was downloaded from NCBI (KX495505–KX495577) (Zhu et al., 2016).

All sequences were assembled and aligned by DNAMAN. After all the sequences were aligned, data were trimmed to 1,689 bp, which contains the full coding sequence of PvAMA-1. We performed a multiple sequence alignment of the PvAMA-1 gene sequences using MEGA6 (Tamura et al., 2013). PvAMA-1 gene sequence PVP01_0934200.1 (1,689 bp) obtained from PlasmoDB (http://PlasmoDB.org) was used as reference.

After all the sequences were aligned, three domains of PvAMA-1 were also divided. These included domain I (462 bp, nucleotides 280–741), domain II (297 bp, nucleotides 793–1,089), and domain III (192 bp, nucleotides 1,162–1,353). To investigate the genetic polymorphism and selection of PvAMA-1 in different time periods in the CMB area, we used DnaSP (Librado and Rozas, 2009) to calculate the number of haplotypes (H), the mean value of nucleotide differences (k), nucleotide diversity (π), and haplotype diversity (Hd) as previously described (Murhandarwati et al., 2020), with a sliding window of 100 bp and step size of 25 bp for nucleotide diversity (π). KaKs_Calculator 2.0 software (Wang et al., 2010) was used to calculate the non-synonymous (Ka) and synonymous (Ks) substitution rates. The Ka and Ks values and Ka/Ks ratios were calculated based on a model-averaged method (Wang et al., 2010). Ka/Ks calculation was used to estimate the selection pressure of PvAMA-1 gene pairs. The algorithm was NG (Nei and Gojobori, 1986) and YN, which is an alternative model for NG (Yang and Nielsen, 2000). Additionally, Tajima’s D test was performed in DnaSP to evaluate the neutrality theory of evolution (with Fu and Li’s test as a double check). The probability of recombination between adjacent nucleotides per generation was calculated using DnaSP. The linkage disequilibrium (LD) between different polymorphic sites was computed based on the D and R ² indices. The calculations performed through DnaSP were based on the default parameters.

The recombination region was calculated by Recombination Detection Program v.4.101 (Martin et al., 2015) with the MaxChi method (Maynard Smith, 1992). After removing the recombination block, a haplotype network based on the Pvama-1 sequence was constructed using the NETWORK software Version 10200 with the median-joining method (Bandelt et al., 1999). To assess allele ancestry, STRUCTURE software was used to assess clustering of isolates under the ancestry model “Use Population Information to test for migrants” (Evanno et al., 2005). Six iterations for the numbers of clusters (K) from three to eight were run, each with a burning period of 5,000 steps and 10,000 Markov chain Monte Carlo iterations. The best K was selected as previously described (Evanno et al., 2005). The sequences of nine additional Pvama-1 populations (Myanmar, KX495505–KX495577; Thailand, FJ784891–FJ785121; Korea, KM230319–KM230384; Sri Lanka, EF218679–EF218701; Iran, JX624732–JX624760; Papua New Guinea (PNG), KC702402–KC702503; India, MH657021–MH657120; Venezuela, EU346015–EU346087; and Brazil, MH049550–MH049589) were analyzed together (Gunasekera et al., 2007; Ord et al., 2008; Putaporntip et al., 2009; Arnott et al., 2013; Zakeri et al., 2013; Kang et al., 2015; Zhu et al., 2016; Bittencourt et al., 2020; Kale et al., 2021). Here, the India population (Kale et al., 2021) should be considered as a long-distance geography distribution outgroup since Kale et al. had exhibited the genetic difference between Myanmar and South Asia. Before the Network and Structure analyses, all the sequences of PvAMA-1 were analyzed and cut to 1,290 bp according to the fragment of PvAMA-1 isolates from Myanmar (Zhu et al., 2016) through MEGA6 and using Arlequin3.5 to analyze the molecular variance (AMOVA) to evaluate fixation (F _ST) (Excoffier and Lischer, 2010). Meanwhile, the migration rate between YP and Myanmar populations was estimated by Migrate-n version 4.4.3 (Beerli, 2006).

Of the 152 Pvama-1 sequences from the CMB area, there were 76 haplotypes, giving an overall haplotype diversity (Hd) of 0.958 ( Table 1 ). The average number of pairwise nucleotide differences (k) for the entire 1,689 bp in different time periods including the full sampling period (referred to as Total), 2009–2011, and 2014–2016 were 12.191, 12.933, and 10.814, respectively ( Table 1 ). The number of haplotypes was 76, 54, and 23 for Total, 2009–2011, and 2014–2016, respectively ( Table 1 ). A total of 61 single nucleotide polymorphisms (SNPs) were detected, including 54 SNPs in 2009–2011 and 47 SNPs in 2014–2016 ( Table 1 ). Nucleotide diversity (π) of the YP samples from 2009 to 2011 and from 2014 to 2016 was 0.00766 and 0.00640, respectively, and that of the total samples was 0.00722 ( Table 1 ). Nucleotide diversity (π) of the YP samples in three domains was 0.01701, 0.00482, and 0.00355, respectively, in which the Pvama-1 domain I showed the highest genetic variation. Recently, Kale et al. (2021) reported the genetic diversity of PvAMA-1 in India and confirmed that the high genetic variation was observed in Pvama-1 domain I. π values for 2009–2011 and 2014–2016 were ranging from 0.000 to 0.03771 (350–380 bp) and from 0.000 to 0.02533 (350–380 bp), respectively ( Figure 2A ).

In order to assess the neutral evolving of PvAMA-1 and its three domains, Tajima’s D test was performed. Tajima’s D value of the full PvAMA-1 fragment was 0.205 for Total, 0.321 for 2009–2011, and 0.311 for 2014–2016 (P > 0.1; Table 1 ). However, the Tajima’s D values obtained for the three domains were obviously different ( Table 1 ). The Tajima’s D value of domain I was 0.979 for Total, 0.991 for 2009–2011, and 0.877 for 2014–2016. The Tajima’s D value of domain II was −0.257, 0.448, and −0.457 for Total, 2009–2011, and 2014–2016, respectively. The Tajima’s D value of domain III was 0.456 for Total, −0.463 for 2009–2011, and 0.618 for 2014–2016 (P > 0.1; Table 1 ). These indicate that Tajima’s D test showed opposite selection directions for different domains of PvAMA-1 across different times ( Figure 2B ). Domain II (from nucleotide positions 795 to 1,089) showed positive Tajima’s D value for 2009–2011 and negative value for 2014–2016, which are consistent with the results in Fu and Li’s tests ( Table 1 and Figures 2C, D ). For the domain III fragment, Tajima’s D value was negative for 2009–2011 and positive for 2014–2016, which were in partial deviation from the analysis results of Fu and Li’s tests with values all positive ( Table 1 and Figures 2C, D ). Meanwhile, a sliding window plot depicted significant positive values (nt 401–500, Tajima’s D: 2.233, P < 0.05) in domain I of YP samples, suggesting positive diversifying selection in this region.

The Ka/Ks in PvAMA-1 for the Total was 1.89662 (YN, P < 0.05). The Ka/Ks values were 3.97228 for 2009–2011 and 3.89083 for 2014–2016 (YN, P < 0.05), suggesting a significant positive selection for PvAMA-1 of P. vivax populations in the CMB area during these times ( Table 1 ). The LD index (R ²) also declined with distance, suggesting that intragenic recombination may also contribute to the PvAMA-1 diversity ( Figure 3 ). Furthermore, we performed the LD test on different time periods and found the 2009–2011 subpopulation under similar pattern for LD SNP pairs with the whole population ( Figure 3 ). In contrast, the 2014–2016 subpopulation showed more LD SNP pairs than the 2009–2011 subpopulation did, which may suggest less recombination ratio from the founder effect ( Figure 3 ). Meanwhile, we analyzed the change in the frequency of haplotypes over time, and the result showed that the frequency of haplotypes decreased significantly in 2015 and 2016 samples ( Figure 4A ).

A total of 40 amino acid mutation sites were found in YP samples. Among them, 37 and 29 amino acid mutation sites were found in 2009–2011 and 2014–2016 samples, respectively. In the Total samples, the mutation frequencies of 19 mutation sites were less than 10%; for 11 mutation sites, 15%–45%; and for 10 mutation sites, more than 50% ( Table 2 ). More importantly, the mutation frequencies of D107N/A, R112K/T, K120R, E145A, E277K, and R438H were more than 70% ( Table 2 ). Some mutation sites with low frequency in the 2009–2011 samples disappeared in the 2014–2016 samples ( Table 2 , Table S1 ). There were also three new mutation sites in the 2014–2016 samples, in which the mutation frequencies of N316T, M319I, and K336R were 0.66%, 7.24%, and 0.66%, respectively ( Table 2 ). Furthermore, there were 19 mutation sites in domain I of PvAMA-1, in which the most mutated amino acid sites were distributed. The number of amino acid mutation sites in domains II and III were seven and three, respectively.

Here, we also analyzed the distribution of these mutation sites of YP samples in other groups ( Table S1 ). In the 37 mutation sites of YP samples, the lowest coincidence rate was found in the Korean samples and the highest coincidence rate in the Iran and India samples ( Table S1 ). At the same time, we also calculated the mutation ratio of these mutation sites of YP samples in the Indian population ( Figure 4B ). It showed that the relatively preserved amino acid changes found in the YP PvAMA-1 were well-conserved in India ( Figure 4B ).

A total of 889 Pvama-1 sequences, containing YP, Myanmar, Thailand, South Korea, Papua New Guinea, Sri Lanka, Iran, India, Venezuela, and Brazil populations, were analyzed and cut to 1,290 bp through MEGA6. The level of genetic differentiation of Pvama-1 was estimated by F _ST values. In general, the F _ST values of <0.05, 0.05–0.15, 0.15–0.25, and >0.25 indicate little, moderate, great, and very great genetic differentiation, respectively (Balloux and Lugon-Moulin, 2002). The YP, Sri Lanka, and Brazil isolates showed great differentiation with F _ST values of 0.15269 and 0.24836, respectively. The very great differentiation was found between YP isolates and Korea isolates (the value of F _ST is 0.3761) ( Table 3 ). The differentiation between YP isolates and other isolates was moderate ( Table 3 ).

One recombination event of YP samples was detected by Recombination Detection Program v.4.101, containing a 602-bp fragment, located in 493–1,094 bp of Pvama-1. After removing the recombination block, a total of 889 sequences from the global populations, length 688 bp, were analyzed by the NETWORK software. Network analysis results for PvAMA-1 populations showed that there was an obvious cluster of populations of global samples, except those from Korea and PNG which were partially separated ( Figure 5A ). Among them, the haplotype sharing ratio of YP samples was the highest. The YP population has shared haplotypes with all other populations except Brazil. Twenty-seven of the 51 YP haplotypes (52.9%) were shared with other populations, of which 59.3% (16/27), 44.4% (12/27), and 44.4% (12/27) were identical to some of the Pvama-1 haplotypes observed in Thailand, Myanmar, and India populations, respectively. Haplotype 286 is the predominant haplotype in the YP population, with a frequency of 17.1%. Haplotype 17 is shared by most populations, consisting of YP, Myanmar, Thailand, Sri Lanka, Iran, and India populations. The haplotype network, drawn by excluding the 102 singletons from the analysis, showed that clusters from the Asian populations, Oceania, and South American overlapped ( Figure 5A ).

In addition, 77 unique haplotypes of CMB isolates were identified in 225 sequences (688 bp) from the 2009–2011, 2014–2016, and Myanmar populations. The PvAMA-1 gene information that was used for comparison in this study was from P. vivax isolates, which had been collected along the CMB area in Laiza, northeast Kachin State, Myanmar, in 2011–2012 (Li et al., 2013). Network analysis results of PvAMA-1 showed an obvious cluster for the three populations ( Figure S1A ). A total of 26.2% 2009–2011 haplotypes (11/42) were shared with the Myanmar (37) population, compared with 38.9% 2014–2016 haplotypes (7/18). Among them, the 2014–2016 population had the least haplotypes, which indicates less recombination in the CMB area.

Through the analysis by MEGA6, a total of 165 SNPs from the 10 populations (YP, Myanmar, Thailand, South Korea, Papua New Guinea, Sri Lanka, Iran, India, Venezuela, and Brazil) were detected. The results of principal component analysis (PCA) showed similar results with the network analysis, where only the Korea and PNG populations were partially separated from the global samples ( Figures 5B, C ).

Furthermore, the structure analysis results of global populations suggested 10 groups with optimized clusters value of K = 7 ( Figure 5D ). The results of structure analysis show that the YP and Thailand samples have the most K components (n = 7), followed by Myanmar and India (n = 6), Sri Lanka and Iran (n = 5), and PNG (n = 4) and the least K components (n = 3) in samples from Korea, Venezuela, and Brazil ( Figure 5D ). The distribution of K components in YP and Myanmar samples was significantly different (χ² = 70.207, P < 0.0001), although they all belong to the CMB area. Only the proportions of K2 component in YP and Myanmar samples are similar (21.71% and 19.17%, respectively) (χ² = 0.191, P = 0.66209). The K components of Myanmar samples were less than those in YP samples ( Figure 5D ), indicating that the YP population was more abundant than the Myanmar population. In addition, the migration rate between YP and Myanmar populations was estimated by Migrate-n version 4.4.3 (Beerli, 2006). The mean Θ values for YP and Myanmar populations were 0.06478 and 0.03100, respectively. The result showed asymmetric gene flow between YP and Myanmar populations, with the number of migrant individuals per generation (Nm) from Myanmar to YP (15.090) being higher than that from YP to Myanmar (11.189).

Meanwhile, the 2009–2011, 2014–2016, and Myanmar samples were analyzed by STRUCTURE software using no admixture model. A total of 62 SNPs from the three populations were detected by MEGA6. The structure results suggested three groups among the CMB samples with optimized clusters value of K = 7 ( Figure S1B ). The results of the structure analysis for the three populations show that the 2014–2016 samples generally lacked some K components that were commonly found in other samples. The distribution of K5 was particularly prominent as shown in Figure S1B . The distribution of K5 in 2009–2011 and 2014–2016 was significantly different (χ² = 15.439, P < 0.0001). Meanwhile, Myanmar samples also have fewer K components than those found in the 2009–2011 samples.