Population genetic structure and natural selection of Plasmodium falciparum apical membrane antigen-1 in Myanmar isolates

Background Plasmodium falciparum apical membrane antigen-1 (PfAMA-1) is one of leading blood stage malaria vaccine candidates. However, genetic variation and antigenic diversity identified in global PfAMA-1 are major hurdles in the development of an effective vaccine based on this antigen. In this study, genetic structure and the effect of natural selection of PfAMA-1 among Myanmar P. falciparum isolates were analysed. Methods Blood samples were collected from 58 Myanmar patients with falciparum malaria. Full-length PfAMA-1 gene was amplified by polymerase chain reaction and cloned into a TA cloning vector. PfAMA-1 sequence of each isolate was sequenced. Polymorphic characteristics and effect of natural selection were analysed with using DNASTAR, MEGA4, and DnaSP programs. Polymorphic nature and natural selection in 459 global PfAMA-1 were also analysed. Results Thirty-seven different haplotypes of PfAMA-1 were identified in 58 Myanmar P. falciparum isolates. Most amino acid changes identified in Myanmar PfAMA-1 were found in domains I and III. Overall patterns of amino acid changes in Myanmar PfAMA-1 were similar to those in global PfAMA-1. However, frequencies of amino acid changes differed by country. Novel amino acid changes in Myanmar PfAMA-1 were also identified. Evidences for natural selection and recombination event were observed in global PfAMA-1. Among 51 commonly identified amino acid changes in global PfAMA-1 sequences, 43 were found in predicted RBC-binding sites, B-cell epitopes, or IUR regions. Conclusions Myanmar PfAMA-1 showed similar patterns of nucleotide diversity and amino acid polymorphisms compared to those of global PfAMA-1. Balancing natural selection and intragenic recombination across PfAMA-1 are likely to play major roles in generating genetic diversity in global PfAMA-1. Most common amino acid changes in global PfAMA-1 were located in predicted B-cell epitopes where high levels of nucleotide diversity and balancing natural selection were found. These results highlight the strong selective pressure of host immunity on the PfAMA-1 gene. These results have significant implications in understanding the nature of Myanmar PfAMA-1 along with global PfAMA-1. They also provide useful information for the development of effective malaria vaccine based on this antigen.


Background
Global incidence and mortality of malaria have decreased in recent years, but malaria is still a great public health concern. In 2015, there were approximately 212 million clinical cases of malaria with an estimated 429,000 deaths [1]. The challenge for malaria control and elimination has increased due to the spread of parasites with resistance to anti-malarial drugs and insecticide-resistant Anopheles mosquitoes. The lack of effective vaccines is one of the major obstacles in combating malaria. Therefore, development of efficacious vaccine is urgently needed for malaria control. Up to date, several candidate proteins including circumsporozoite protein (CSP), Duffy-binding protein (DBP), merozoite surface protein-1 (MSP-1), apical membrane antigen-1 (AMA-1), and thrombospondin related anonymous protein (TRAP) have been tested for their potentials as candidate antigens to develop effective vaccines [2]. The majority of these antigens are expressed either in pre-erythrocytic and erythrocytic stages [2,3]. However, genetic polymorphisms identified in these parasite antigens are great hurdles to develop effective vaccines since they generate a variant-specific immune response which is less effective against parasites with other genetic variants [4][5][6].
Apical membrane antigen-1 (AMA-1) is a 83-kDa type I integral membrane protein that is mainly expressed in the merozoite and sporozoite [7,8]. Biological function of AMA-1 is not clearly understood yet, but its stage specific expression and localization suggest that the protein might play a crucial role in invasion of erythrocytes and hepatocytes by malaria parasites [8][9][10][11]. Structure analysis of AMA-1 revealed that this protein consists of a signal sequence, a cysteine-rich ectodomain, a conserved cytoplasmic region and a transmembrane region [12]. The ectodomain of AMA-1 is further subdivided into three domains, domains I, II, and III. The ectodomain AMA-1 is very immunogenic and natural immune responses against the domain have been reported in patients exposed to Plasmodium falciparum and Plasmodium vivax [13][14][15][16]. Immunization with AMA-1 can produce antibodies to effectively inhibit the invasion of erythrocyte by malaria parasite and confer protective immune responses [15,17]. As a result, AMA-1 is a leading blood stage vaccine candidate for malaria control [16,18,19]. Several studies have shown that anti-AMA-1 antibodies confer protective immunity in adults living in malaria-endemic areas [20,21]. However, these antibodies recognize either conserved or allele-specific epitopes of AMA-1, resulting limited protection against different alleles [22]. AMA-1 is known to be less variable than other malaria vaccine candidate antigens such as MSPs or CSP, but AMA-1 also possesses genetic diversity among global malaria parasites. A high rate of polymorphisms has been identified in domain I of AMA-1 and this region appears to be a major target of anti-AMA-1 protective antibodies [23][24][25][26][27][28]. To design an efficient and protective malaria vaccine, it is essential to monitor genetic variations of vaccine candidate antigens among global malaria isolates circulating in endemic areas [29].
The incidence of malaria in Myanmar has decreased in recent years. However, Myanmar still accounts for more than half of the malaria cases and approximately three quarters of the malaria deaths in the Greater Mekong Subregion [30]. Although P. falciparum cases have decreased in recent years [1,31], this species is still the most critical concern for malaria control and prevention in Myanmar with the emergence of artemisinin resistance [32]. In this study, population genetic structure and natural selection of P. falciparum AMA-1 (PfAMA-1) in Myanmar P. falciparum isolates were analysed.

Blood samples and ethics
Blood samples used in this study were obtained from P. falciparum infected patients who resided in rural villages located in Naung Cho, Shan State, and Pyin Oo Lwin and Tha Beik Kyin, Mandalay Division, Upper Myanmar between July to December 2015 (Fig. 1). Malaria transmission was heterogeneous and seasonal in these areas and most malaria cases occurred with the peak during and just after the rainy season. The major ethnic groups in these regions were Bamar and Shan. Communitybased survey was performed in 825 residents in 7 villages of the areas and P. falciparum infection was confirmed by Giemsa-stained thick and thin blood smear examination. The average age of patients was 34.7 years old ranged from 13 to 60 years. Two or three ml of venous blood was collected from the P. falciparum infected patients in EDTA tubes before drug treatment for further molecular analysis. All P. falciparum positive samples were further confirmed by polymerase chain reaction (PCR) targeting 18S ribosomal RNA (rRNA) gene [31,33]. The use of blood samples in this study was approved by the Ministry of Health, Myanmar (Approval Number: 97/Ethics 2015). It was also approved by Biomedical Research Ethics Review Board of Inha University College of Medicine, Republic of Korea (Approval Number: INHA 15-013). Written consent was obtained from each individual prior to blood collection.

Amplification and sequencing analysis of PfAMA-1
Genomic DNA of parasite was extracted from 200 µl of whole blood sample of patients using QIAamp DNA Blood Kit (Qiagen, Hilden, Germany) following the manufacturer's instruction. Full-length PfAMA-1 gene was amplified by PCR. PfAMA-1 gene-specific primers were designed based on complete coding sequences of PfAMA-1 of P. falciparum isolates 3D7 (GenBank Accession No.: U65407) and FC27 (GenBank Accession No.: M27133). Forward and reverse primers were 5′-GTACTTGTTATAAATTGTACA-3′ and 5′-TTTT AGCATAAAAGAGAAGC-3′, respectively. Thermal cycling parameters for PCR were as follows: one cycle of an initial denaturation at 94 °C for 5 min, 30 cycles of 94 °C for 1 min, annealing at 55 °C for 1 min and extension at 72 °C for 2 min, followed by a final extension step at 72 °C for 10 min. Nested PCR with primers 5′-ACAAAAATGAGAAAATTATACTGC-3′ and 5′-TTAATAGTATGGTTTTTCCATCAGAAC-3′ was performed with the same amplification condition, if necessary. To minimize nucleotide mis-incorporation into sequences during PCR amplification, Ex Taq DNA polymerase (Takara, Otsu, Japan) with proof-reading activity was used in all PCR. Each PCR product was analysed by electrophoresis on 1% agarose gels. The resulting PCR product was cut from the gel followed by gel purification.  (Hd), nucleotide diversity (π), and average number of pair-wise nucleotide differences within a population (K) were estimated using DnaSP ver. 5.10.00 [34]. Value of π was calculated to estimate step-wise diversity throughout the full-length PfAMA-1 based on a sliding window of 100 bases with a step size of 25 bp. Values of non-synonymous (dN) and synonymous (dS) substitutions were estimated and were compared using Z test (P < 0.05 was considered significant) in MEGA4 program [35] based on method of Nei and Gojobori [36] with Jukes and Cantor correction. The Tajima's D value [37] and Fu and Li's D and F values [38] were analysed using DnaSP ver. 5.10.00 to evaluate neutral theory of evolution [34]. Recombination parameter (R), which included the effective population size and probability of recombination between adjacent nucleotides per generation, and minimum number of recombination events (Rm) were analysed by DnaSP ver. 5.10.00 [34]. Linkage disequilibrium (LD) between different polymorphic sites was computed based on the R 2 index using DnaSP ver. 5.10.00 [34].

Amino acid polymorphisms in Myanmar PfAMA-1 compared to global PfAMA-1
When patterns of amino acid polymorphisms identified in Myanmar PfAMA-1 were compared to those of PfAMA-1 from other countries, similar but not identical polymorphic patterns were observed. Most amino acid changes identified in global PfAMA-1 were found in domains I and III. Amino acid changes were also found in domain II, 5ʹ-terminal (5ʹ-T), and 3ʹ-terminal (3ʹ-T) regions (Fig. 3). Overall patterns of amino acid changes observed in global PfAMA-1 were similar to each other. However, frequencies of each amino acid change differed by country. The most notable amino acid changes identified in Myanmar PfAMA-1 were that three amino acid changes (R39H, Y175D and P330S) were fixed in all Myanmar PfAMA-1 sequences, even though their frequencies were also high in PfAMA-1 of other global isolates. Frequency (%)  were higher than those of Asian and Pacific PfAMA-1. Nucleotide diversity across the full-length PfAMA-1 from different geographical areas also slightly differed by geographical areas ( Table 2). The level of nucleotide diversity across PfAMA-1 in Myanmar P. falciparum isolates (π = 0.01064) was lower than that for PfAMA-1 from Ghana (π = 0.01405), PNG (π = 0.01226), Tanzania (π = 0.01361), or Solomon Islands (π = 0.01179), but similar to that for PfAMA-1 from Thailand (π = 0.01103), Philippines (π = 0.01136), or Vanuatu (π = 0.01012). A sliding window plot of π revealed that PfAMA-1 from different geographical areas shared highly similar patterns of nucleotide diversity through their sequences with two peaks at domains I and III. The maximum diversity was found at domain I (Fig. 4a). All PfAMA-1 sequences from different countries showed positive Tajima's D values, suggesting that balancing selection of global PfAMA-1 ( Table 2). A sliding window plot analysis also showed that global PfAMA-1 had similar pattern for Tajima's D across the gene, albeit some differences were identified among PfAMA-1 with different geographical origins (Fig. 4b).

Recombination and linkage disequilibrium
The minimum number of recombination events between adjacent polymorphic sites (Rm) for Myanmar PfAMA-1 was estimated to be 24. The R values between adjacent sites (Ra) and per gene (Rb) were 0.0141 and 26.3, respectively. Possible recombination events were also identified in global PfAMA-1. The highest R values were predicted for African PfAMA-1 (Ghana and Tanzania), while the lowest R values were predicted for Vanuatu PfAMA-1 (Table 3). LD index (R 2 ) for global PfAMA-1 also decreased with increasing distance across the gene (Fig. 5).

Haplotype network analysis
Haplotype network analysis of PfAMA-1 haplotypes from global P. falciparum populations showed a dense network with extremely complex relationships (Fig. 6)

Nucleotide differentiation among global PfAMA-1
To investigate the degree of gene flow and genetic differentiation among global PfAMA-1 sequences, Fst values were evaluated for full-length PfAMA-1 sequences from different geographical P. falciparum populations deposited at GenBank (Table 4). Fst values between different geographical PfAMA-1 populations varied from 0.00019 (between Ghana and Tanzania) to 0.14141 (between Vanuatu and Tanzania).

Association between natural selection and host immune pressure
The selective pressure of host immunity on PfAMA-1 was evaluated by analysing genetic polymorphisms in predicted RBC-binding sites, B-cell epitopes and IUR regions. Results showed that most amino acid changes caused by SNPs were found in predicted RBC-binding sites, B-cell epitopes, or IUR regions of PfAMA-1 (Fig. 7a). Of 83 amino acid changes, 66 were found in global PfAMA-1 compared to PfAMA-1 of 3D7 (Gen-Bank Accession No.: U65407) at predicted RBC-binding sites, B-cell epitopes, or IUR regions. Among 51 commonly identified amino acid changes in global PfAMA-1, 43 were found in predicted RBC-binding sites, B-cell epitopes, or IUR regions. The 23 less commonly changed amino acids in global PfAMA-1 sequences were also found in predicted RBC-binding sites, B-cell epitopes, or IUR regions. Eight of 11 predicted B-cell epitopes were polymorphic. In particular, high levels of π were predicted for B-cell epitopes 3, 4, 5, 8, and 10. These regions contained major polymorphic amino acid residues in global PfAMA-1 (Fig. 7b). Tajima's D values for the predicted B-cell epitopes 3, 4, and 8 were all positive, indicating a decrease in population size and/or balancing selection. Meanwhile, Tajima's D values for the predicted B-cell epitopes 5 and 10 were negative (Fig. 7b). Common amino acid changes found in global PfAMA-1 were located at C1-L covering region, which is located near the hydrophobic pocket of PfAMA-1. Meanwhile, only a few less frequent amino acid changes were identified in loop II of PfAMA-1.

Discussion
Population genetic structure and natural selection of PfAMA-1 among Myanmar P. falciparum isolates were analysed in this study. Most of amino acid changes identified in Myanmar PfAMA-1 were clustered at domains I and III, consistent with previous reports on PfAMA-1 from other endemic areas [40,[42][43][44][45]. Overall distribution patterns and frequencies of amino acid changes found in Myanmar PfAMA-1 were similar to those of global PfAMA-1 analysed in this study, but several differences among and between global PfAMA-1 were also identified. The implication of these substantial geographic differentiations is currently unclear. Considering that only a limited number of PfAMA-1 sequences in each geographical area were analysed in this study, these amino acid changes and their different frequencies  may not be statistically significant. Indeed, many amino acid changes found in global PfAMA-1 were commonly detected and evenly distributed among global PfAMA-1 sequences, although their frequencies differed among and between populations. Further examination of PfAMA-1 nucleotide and amino acid variations in diverse PfAMA-1 populations with a larger number of global PfAMA-1 sequences would be necessary to better understand the polymorphic nature of PfAMA-1. It is also worthy to mention that Myanmar PfAMA-1 sequences analysed in this study were from P. falciparum isolates that collected in restricted areas of Myanmar. Therefore, nation-wide analysis of PfAMA-1 in P. falciparum isolates collected from various regions of Myanmar is also needed to clearly understand the genetic diversity and population structure of Myanmar PfAMA-1. The π value for global PfAMA-1 varied, ranging from 0.01012 (Vanuatu) to 0.01405 (Ghana). The π values in African PfAMA-1 sequences (Ghana and Tanzania) were relatively higher than those in Asian or Pacific PfAMA-1 sequences. Nucleotide diversity was not evenly distributed through Myanmar PfAMA-1. The 5ʹ-T, domain II and 3ʹ-T showed low levels of nucleotide diversity, indicating that these regions might be more conserved with lower frequency of polymorphisms. Meanwhile, much higher value of nucleotide diversity was observed at domain I and domain III. Highly similar patterns of nucleotide diversity were also detected in global PfAMA-1 analysed in this study, strongly indicating that global PfAMA-1 might share highly similar nucleotide diversity across this gene. The dN − dS value for Myanmar full-length PfAMA-1 was positive, implying balancing selection in this gene. However, the dN − dS value for 5ʹ-T was negative, suggesting negative or purifying natural selection in this domain. Tajima's D value (the statistic value predicting departure from neutrality) for Myanmar PfAMA-1 also suggested that this gene might have evolved under balancing selection. In sliding plot analysis of Tajima's D, general patterns of Tajima's D values across Myanmar PfAMA-1 were similar to those of global PfAMA-1, although differences among and between PfAMA-1 from different countries were also found. Regardless of theses slight differences, domains I and III and the junction between domains II and III shared very similar patterns of Tajima Table 3 Comparison  of  recombination  events  between global PfAMA-1 The R and Rm were estimated excluding the sites containing alignment gaps or those segregating for three nucleotides. The R was computed using R = 4Nr, where N is the population size and r is the recombination rate per sequence (per gene) n, number of isolates; Ra, recombination parameter between adjacent sites; Rb, recombination parameter for entire gene; Rm, minimum number of recombination events between adjacent sites suggest that domains I and III of Myanmar PfAMA-1 are highly polymorphic with strong selection force acting on these domains, similar to other global PfAMA-1. Meiotic recombination between the adjacent polymorphic sites is one of main forces that drive allelic diversity in Plasmodium parasites [40]. Substantial levels of recombination events in PfAMA-1 from different geographical isolates have been reported previously [40,43,[45][46][47][48][49]. Recombination events within Myanmar PfAMA-1 and decline of LD index R 2 with the increase of nucleotide distance of the sequences evidenced meiotic recombination in Myanmar PfAMA-1, consistent with previous studies      Fst is a measure of population substructure and is one of the most useful tools for analysing the overall genetic differentiation among populations. Fst values at each locus are considered as no differentiation (0), low genetic differentiation (0-0.05), moderate differentiation (0.05-0.15), or high differentiation (0.15-0.25) [50]. Fst values for Myanmar PfAMA-1 in comparison with PfAMA-1 from other different geographical populations showed moderate levels of genetic differentiation, ranging from 0.058281 to 0.09740. Fst values between PfAMA-1 populations belonging to the same geographical area were relatively low. In particular, the value between Ghana and Tanzania was very low (0.00019), suggesting very low genetic differentiation among these two populations. However, Vanuatu showed relatively high Fst values for Ghana (0.12425) and Tanzania (0.14141), providing evidence for geographical isolation followed by population subdivision in these regions. Overall Fst values between global PfAMA-1 were in moderate differentiation ranges, indicating limited genetic differentiation of PfAMA-1 among global parasite populations.
The significant ratio of dN/dS and evidence of positive natural selection identified in PfAMA-1 of Myanmar and the other seven P. falciparum populations suggest that amino acid replacements are generally favored in PfAMA-1 and this adaptive evolution is presumably due to host immune pressure. Previous studies have suggested that polymorphic amino acid residues identified in PfAMA-1 are located on one side of the molecule, thereby exposing them to the exterior environment more easily with more access to host immune responses [40,43,51]. A recent study has provided information on potential RBC-binding regions, B-cell epitopes, and IURs across the ectodomain of PfAMA-1 [40]. To assess the association between natural selection and host immune pressure for PfAMA-1, genetic polymorphisms in predicted RBC-binding sites, B-cell epitopes, and IUR regions were analysed. Most amino acid changes found in global PfAMA-1 were predicted to be localized at the predicted RBC-binding sites, B-cell epitopes, or IUR regions of PfAMA-1. Eight of 11 predicted B-cell epitopes were polymorphic. In particular, B-cell epitopes 3, 4, 5, 8, and 10 (corresponding to domain I, junction between domain II and III, and domain III) showed high levels of nucleotide diversity in global PfAMA-1 analysed. Tajima's D values for these predicted B-cells epitopes also suggested that these epitopes were under natural selection. The C1-L cluster, which located near the hydrophobic pocket of domain I in PfAMA-1, affects binding capacity of inhibitory monoclonal antibody, thereby contributing to antibody escape [52][53][54][55]. Global PfAMA-1 analysed in this study had clustering amino acid polymorphisms, including tri-and hepta-morphic changes in C1-L region. This supports the notion that strong natural selection force which generates genetic diversity of PfAMA-1, thereby contributing to host immune escape acts on this region [21,54,55]. Meanwhile, only a few less frequent dimorphic amino acid changes were identified in loop II of global PfAMA-1, the target of 4G2 inhibitory antibody [56], coinciding with results of a previous study [40]. High degree of amino acid sequence conservation of the loop II among global PfAMA-1 further supports the hypothesis that amino acid sequences in this region are subjected to strong natural selection. Therefore, this region might be useful as a component of PfAMA-1based vaccine.
To develop a globally effective malaria vaccine based on PfAMA-1, it would be important to understand whether PfAMA-1 haplotypes with distinct clusters are antigenically different. Haplotype network analysis of 517 global PfAMA-1 sequences suggested that a total of 174 distinct haplotypes were identified among global population. No haplotype that fully covers haplotypes from all geographical areas analyzed in this study was identified. PfAMA-1 based vaccines currently studied are based on sequences of most extensively studied forms of P. falciparum, 3D7 and FVO [57][58][59]. However, only 13.7% of global PfAMA-1 haplotypes analysed in this study was identical to the FVO. The frequency of global PfAMA-1 haplotypes with identical sequence to 3D7 was much lower. These raised a concern that vaccine based on these two forms of haplotypes may not offer effective protection for global P. falciparum isolates. A recent report suggested that PfAMA-1 sequence variations may not necessarily be strong predictors of antigenic differences or the level of cross-inhibitory antibody activity because not all polymorphic residues contribute equally to antibody binding and escape [60]. Antigenic diversity of PfAMA-1 is limited and a vaccine targeting a small number of PfAMA-1 alleles might be sufficient for coverage against naturally-circulating P. falciparum populations in different endemic areas [27]. The results of this study also suggested that the genetic diversity is relatively limited among global PfAMA-1, even though substantial geographic differentiations were also identified among the populations. However, global PfAMA-1 is under natural selection and high level of meiotic recombination that can produce new alleles occurs in the gene population. And therefore, a multi-allele approach for developing PfAMA-1 based malaria vaccines should be considered to maximize vaccine efficacy.