Low genetic diversity and functional constraint in loci encoding Plasmodium vivax P12 and P38 proteins in the Colombian population

Plasmodium vivax is one of the five species causing malaria in human beings, affecting around 391 million people annually. The development of an anti-malarial vaccine has been proposed as an alternative for controlling this disease. However, its development has been hampered by allele-specific responses produced by the high genetic diversity shown by some parasite antigens. Evaluating these antigens’ genetic diversity is thus essential when designing a completely effective vaccine. The gene sequences of Plasmodium vivax p12 (pv12) and p38 (pv38), obtained from field isolates in Colombia, were used for evaluating haplotype polymorphism and distribution by population genetics analysis. The evolutionary forces generating the variation pattern so observed were also determined. Both pv12 and pv38 were shown to have low genetic diversity. The neutral model for pv12 could not be discarded, whilst polymorphism in pv38 was maintained by balanced selection restricted to the gene’s 5′ region. Both encoded proteins seemed to have functional/structural constraints due to the presence of s48/45 domains, which were seen to be highly conserved. Due to the role that malaria parasite P12 and P38 proteins seem to play during invasion in Plasmodium species, added to the Pv12 and Pv38 antigenic characteristics and the low genetic diversity observed, these proteins might be good candidates to be evaluated in the design of a multistage/multi-antigen vaccine.


Background
Malaria is a disease caused by protozoan parasites from the Plasmodium genus, five of which cause the disease in human beings (Plasmodium falciparum, Plasmodium vivax, Plasmodium ovale, Plasmodium malariae and Plasmodium knowlesi) [1,2]. This parasite is transmitted by the bite of an infected Anopheles female mosquito. Around 3.3 billon people are at risk of malaria annually, mainly in tropical and subtropical areas of the world, children aged less than five years and pregnant women being the most vulnerable [3]. Plasmodium falciparum is responsible for the disease's most lethal form, being predominantly found on the African continent whilst P. vivax is widely distributed around the world. Even though it has been thought that infection caused by the latter species was benign, recent studies have shown that P. vivax can cause clinical complications [4]. It has been found that 2,488 million people are at risk of becoming infected by P. vivax on the continents of Asia and America, 132 to 391 million cases occurring annually [5].
In spite of control strategies having been introduced in different countries, malaria continues to be a public health problem due to the parasite's resistance to anti-malarial treatments [6] and the vector's resistance to insecticides [7], among other causes. More effective measures have thus to be implemented for controlling such disease, including the development of an anti-malarial vaccine.
Several antigens have been characterized as promising candidates for inclusion in a vaccine [8,9], however, the genetic diversity of some of them [10][11][12][13][14][15][16][17][18] has hampered the development of such vaccine [19,20] as these genetic variations provoke allele-specific responses [21,22] making them become a mechanism for evading the immune system [23]. It has been necessary to focus vaccine development on conserved domains or antigens to avoid such responses [24], since these regions could have functional constraint and have had slower evolution [25].
Developing a multi-antigen vaccine against the parasite's blood stage has been focused on blocking all hostpathogen interactions to stop merozoite entry to red blood cells (RBC) [26]. A group of proteins anchored to the membrane via glycosylphosphatidylinositol (GPI) has been identified in P. falciparum, predominantly located in detergent-resistant membrane (DRM) domains [27,28]; they have been implicated in the parasite's initial interaction with RBC [29][30][31][32][33] and some have been considered as being candidates for being included in a vaccine [34,35]. One group of proteins belonging to the 6-cystein (6-Cys) family is particularly noteworthy among these DRMs (i.e., Pf12, Pf38, Pf41 and Pf92) as they have been characterized by having s48/45 domains (ID in PFAM: PF07422). Members of this family are expressed during different parasite stages [28,36] and some of them (e.g., Pf48/45, Pf230) have been considered as vaccine candidates for the sexual stage [36,37].
Pf12 and Pf38 are expressed during late stages of the intra-erythrocyte cycle, each having two high binding peptides, suggesting an active role during invasion of RBC [30]. Orthologous genes encoding these proteins have been characterized recently in P. vivax [38,39]. Both proteins have a signal peptide, a GPI anchor sequence and have been associated with DRMs [38,39]. Pv12 has two s48/45 domains [39] whilst Pv38 has a single domain located towards the C-terminal end [38]. These proteins have been shown to be antigenic [38][39][40], suggesting that they are exposed to the immune system, probably during P. vivax invasion of RBC.
The present study involved a population genetics analysis for evaluating the genetic diversity of pv12 and pv38 loci and the evolutionary processes generating this variation pattern; the results revealed these antigens' low genetic diversity in the Colombian population, possibly due to functional/structural constraints in s48/45 domains. Since the proteins encoded by these genes share structural characteristics with other vaccine candidates, added to the fact that Pv12 and Pv38 are targets for the immune response [38][39][40] and have conserved domains, they should be considered when designing a multistage/ multi-antigen anti-malarial vaccine.

Ethics statement
The parasitized DNA used in this study was extracted from total blood collected from different Colombian areas (Antioquia, Atlántico, Bogotá, Caquetá, Cordoba, Chocó, Guainía, Guaviare, Magdalena, Meta, Nariño, and Tolima) from 2007 to 2010. All P. vivax-infected patients who provided blood samples were notified about the object of the study and signed an informed consent form if they agreed to participate. All procedures involved in taking blood samples were approved by Fundación Instituto de Inmunología de Colombia (FIDIC) ethics committee.
DnaSP software (v.5) [47] was used for evaluating intrapopulation genetic polymorphism by calculating: the number of polymorphic segregating sites (Ss), the number of singleton sites (s), the number of parsimony-informative sites (Ps), the number of haplotypes (H), haplotype diversity (Hd, which was multiplied by (n-1)/n according to Depaulis and Veuille [47,48]), the Watterson estimator (θw) and nucleotide diversity per site (π). DNA sequence variation was calculated using the sequences obtained from the aforementioned databases, plus the Colombian ones (worldwide isolates, global diversity) and just those obtained for the Colombian population (local diversity). The frequency for each Colombian haplotype was also estimated by count and year.
Two test families were used for evaluating the neutral molecular evolution model for the Colombian population: (1) frequency spectrum test, and (2) haplotype test. The former involved calculating Tajima's D statistics [49], Fu and Li's D* and F* [50] and Fay and Wu's H statistic [51]. Tajima's D statistic compares the difference between segregating sites and the average of nucleotide differences between two randomly taken sequences. Fu and Li's D* statistic takes the difference between the number of singleton sites and the total of mutations, whilst F* takes the difference between the number of singleton sites and the average of nucleotide differences between two randomly taken sequences. Fay and Wu's H statistic is based on the difference of the average number of nucleotide differences between pairs of sequences and the frequency of the derived variants. Fu's Fs statistic [52], K-test and Htest [48] are tests for calculating haplotype distribution. The Fs statistic compares the number of haplotypes observed to the expected number of haplotypes in a random sample. K-test and H-test [48] are based on haplotype number and haplotype diversity, respectively; these statistics are conditioned by sample size (n) and the number of segregating sites (Ss). Test significance was determined by coalescence simulations using DnaSP (v.5) [47] and ALLE-LIX software (kindly supplied by Dr Sylvain Mousset). Sites having gaps were not taken into account in any of the tests performed.
The effect of natural selection was evaluated regarding intra and interspecies; the average number of nonsynonymous substitutions per non-synonymous site (d N ) and the average number of synonymous substitutions per synonymous site (d S ) were calculated for the former by using the modified Nei-Gojobori method [53]. The significant differences between the above were determined by using Fisher's exact test (suitable for d N and d S < 10) and codon-based Z-test incorporated in MEGA software (v.5) [54]. Differences between d N and d S per site were calculated by using SLAC, FEL, REL [55], IFEL [56], MEME [57], and FUBAR [58] methods. The average number of non-synonymous divergence substitutions per nonsynonymous site (K N ) and the average number of synonymous divergence substitutions per synonymous site (K S ) were calculated using the modified Nei-Gojobori method [53], with Jukes-Cantor correction [59], to infer natural selection signals which may have prevailed during malarial parasite evolutionary history (interspecies; using Plasmodium cynomolgi (accession number BAEJ 01001076.1) and P. knowlesi (accession number NC_011 912.1) orthologous sequences). The significant differences between K N and K S were determined by using a codon-based Z-test incorporated in MEGA software (v.5) [54]. The McDonald-Kreitman test [60] was also calculated; this is based on a comparison of intraspecific polymorphism to interspecific divergence (using Plasmodium cynomolgi (accession number BAEJ01001076.1) and P. knowlesi (accession number NC_011912.1) orthologous sequences). This test involved using a web server [61], which takes Jukes-Cantor divergence correction into account [59]. All the above tests were calculated using the sequences obtained from the databases plus the Colombian ones and just those obtained for the Colombian population.
Z nS [62] and ZZ [63] statistics were calculated for evaluating the influence of linkage disequilibrium (LD) and intragenic recombination, respectively. The minimum number of recombination (Rm) events was also calculated; this included calculating effective population size and the probability of recombination between adjacent nucleotides per generation [64]. Additionally, the GARD method [65] available at the Datamonkey web server [66] was performed. These tests were performed using the sequences obtained from the Colombian population.

Results and discussion
The presence of genomic DNA (gDNA) and identification of single Plasmodium vivax strain infection An 18S subunit rRNA gene fragment was amplified from 100 samples of P. vivax collected from different areas of Colombia and stored from 2007 to 2010. Seventy-seven samples revealed an amplicon at the expected size, indicating the presence of P. vivax gDNA. A region of the pvmsp-1 gene was then amplified and digested with restriction enzymes, showing that seven of the 77 samples proving positive for P. vivax had multiple infections. Only 70 samples were thus considered for later analysis. Due to the low number of samples collected from some areas, they were grouped according to geographical localisation and epidemiological conditions (South-west: Chocó, Nariño; South-east: Caquetá, Guainía, Guaviare, Meta; Midwest: Bogota, Tolima; North-west: Atlántico, Antioquia Cordoba, Magdalena).

Genetic diversity in pv12
Seventy samples amplified a 1,200 base pair (bp) fragment corresponding to the pv12 gene (South-west n = 6; South-east: n = 20; Midwest: n = 8; North-west: n =36). These amplicons were purified and sequenced; the sequences were then analysed, compared to different reference sequences obtained from various sequencing projects [43,44] and those having a different haplotype were deposited in the GenBank database (accession numbers KF667328 and KF667329).
Four single nucleotide polymorphisms (SNP) were observed throughout the pv12 gene sequence ( Figure 1A Six haplotypes were found in pv12 ( Figure 1A and Table 1) around the world, four of which are present in Colombia at 8.7, 5.8, 10.1, and 75.4% frequency for haplotypes 2, 3, 5 and 6, respectively. Haplotypes 2, 5 and 6 were present in the different Colombian locations (Additional file 1), haplotype 6 being the most predominant per year (2007 n = 9; 2008 n = 17; 2009 n = 15; 2010 n = 29) and per location, having higher than 70% frequency ( Figure 2A and Additional file 1). The remaining haplotypes were absent or had low frequency (Figure 2A and Additional file 1). Interestingly, haplotype 3 was present in Colombia during 2009 but absent in the other years studied (Figure 2A). The percentage of samples from the South-east area (some of them presenting haplotype 3) was greater than for other years, suggesting that haplotype 3 was restricted to a particular geographical area (Additional file 1) and/or that this had very low frequency in different Colombian subpopulations. Haplotype 2 was absent from 2007 to 2008 but present between 2009 and 2010 ( Figure 2A); differently to haplotype 3, this haplotype was present everywhere, except in the Southwest location (Additional file 1). This appeared to be consistent with previous studies which have reported numerous private haplotypes in American Plasmodium vivax populations [67]. These results suggested that the Colombian population had one predominant pv12 haplotype and several low frequency alleles, which are geographically isolated or were not detected during some periods of time. Since P. vivax populations within countries seem to be strongly structured [67], new pv12 haplotypes could appear in other parasite populations.
This gene had 0.0004 ± 0.0001 global nucleotide diversity (π) and 0.0003 ± 0.0001 for the Colombian population (Table 1). This value was about 2.5 times less than that reported for its orthologue in P. falciparum (π = 0.001) [68]; however, both values were low when compared to other membrane proteins [10][11][12][13][14]17], suggesting that this gene is highly conserved in different Plasmodium species. This value places pv12 among the most conserved antigen-encoding genes characterized to date in P. vivax.

Mutations in pv12 appear to be selectively neutral
Several tests for evaluating the hypothesis that mutations in pv12 are neutral were performed. No significant values were found for the Tajima, Fu and Li, Fay and Wu or Fu tests (Table 2); likewise, the Colombian population's number of haplotypes (4) and haplotype diversity (0.406 ± 0.07) ( Table 2) were as expected under neutrality according to the K-test and H-test. Since neutrality could not be ruled out, the mutations or haplotypes found in pv12 could have been randomly fixed; this might explain the possible geographical isolation of haplotype 3, since different alleles could have become fixed in different populations according to the neutral hypothesis. Alternatively, the geographical isolation of haplotype 3 could have resulted from the structured P. vivax population in America, where haplotypes may have diversified in situ.

Natural selection in pv12
The gene was split into two regions: region A, nucleotides 1-546 (amino acids 1-182 including one s48/45 domain) and region B, nucleotides 547-1,095 (amino acids 183-365 including the other s48/45 domain). Synonymous substitution per synonymous site (d S ) and non-synonymous substitution per non-synonymous site rates (d N ) were calculated using the gene's total length to evaluate whether natural selection had any effect on pv12 evolution. Full length gene and split regions had non-significant values (Table 3); likewise, when d N and d S were estimated for s48/45 domains, no significant values were observed (Additional file 2), contrary to that suggested for pf12, where purifying selection action has been reported [69]. The Datamonkey server was used for calculating d N and d S rates for each codon; no selected sites were found, indicating (once more) that the gene did not appear to deviate from neutrality.
However, assessing how natural selection acts on low genetic diversity antigens is not easy [70]; the fact that Plasmodium vivax shares its most recent common ancestor with parasites infecting primates (e.g. P. cynomolgi and P. knowlesi) led to inferring patterns which may have prevailed during their evolutionary history [70,71]. When synonymous divergence substitution per synonymous site (K S ) and non-synonymous divergence substitution per non-synonymous site (K N ) rates were calculated, a significantly higher K S than K N was found (Table 4). Moreover, a sliding window for ω (d N /d S and/or K N /K S ) revealed < 1 values throughout the gene (Figure 3), which could have been a consequence of negative selection. Moreover, significant values were observed when the McDonald-Kreitman (MK) test was used for comparing intraspecific polymorphism and interspecific divergence (using all the haplotypes found for this gene): P N /P S > D N /D S (Table 5), revealing (similar to K S rates) a large accumulation of synonymous substitutions between species, which could be interpreted as negative selection. Such accumulation of interspecies synonymous substitutions suggested that evolution tried to maintain protein structure by eliminating all deleterious mutations. However, when the MK test was done with haplotypes found in Colombia (and in spite of the accumulation of synonymous substitutions between species), no significant values were observed in this population (Table 5). Although Pv12 is exposed to the immune system [39,40], it had a high level of conservation. This pattern could have been because pv12 had diverged  by negative selection, due to a possible functional/structural constraint imposed by the presence of s48/45 domains [72] which seem to play an important role during host cell recognition [30,69,72].

Genetic diversity in pv38
Only 46 out of 70 samples could be amplified for the pv38 gene, giving a 1,121 bp fragment (South-west n = 6; South-east: n = 13; Midwest: n = 4; North-west: n = 23). The 46 sequences obtained from Colombian isolates were compared to and analysed regarding reference sequences obtained from different regions of the world [43,44]. Colombian sequences that have a different haplotype to that of previously reported ones can be found in GenBank (accession numbers KF667330-KF667340). Nine SNPs were observed in the pv38 gene ( Figure 1B), most of which were no-synonymous (nucleotides: 88 (R30S), 206/207 (A69V), 209 (R70L), 524/525 (T175N), 880 (M294L), and 998 (S333N)), similar to that found in Pf38 [73]. Positions 525 and 969 produced synonymous substitutions (a change in protein sequence was generated when the substitution in position 525 was accompanied with another one in position 524). The parasite population in Colombia has eight of these nine SNPs, all being informative-parsimonious sites. Similar to that reported for its orthologue in P. falciparum [73], most substitutions were found in the gene's 5′ region.
Seventeen haplotypes were identified from alignment (including sequences from different regions of the world) ( Figure 1B) Figure 2B). The absence of some haplotypes in determined years, or in some locations, could not just have been due to the low frequency which they might have had but also to the difference in the number of samples for each year (n = 6 in 2007, n = 6 in 2008, n = 8 in 2009 and n = 26 in 2010) or because American P. vivax populations appear to be structured and therefore several privative haplotypes might be found [67].
π in this gene was 0.0026 ± 0.0002 worldwide and 0.0024 ± 0.0002 in the Colombian population (Table 1), this being 1.3 times lower than that for its orthologue in P. falciparum (π = 0.0034) [68,73] showing that the pv38   Table 2). This suggested balanced ancestral polymorphism [48], being similar to that reported for the P. falciparum p38 gene which showed evidence of balanced selection in 5' region [73].

Natural selection in pv38
A modified Nei Gojobori method was used for calculating d N and d S rates for showing some type of selection in the pv38 gene. Similar to that used regarding pv12, the pv38 gene was divided into two regions: region A, covering position 1-459 (amino acids 1-153) and region B, nucleotides 460-1,065 (amino acids 154-355 including   the s48/45 domain). There were more d N substitutions in region A than d S substitutions, whilst there were more d S substitutions in region B than d N ones, even though no significant values were observed (Table 4 and Additional file 2). Selection tests by codon revealed positive selection in codon 70 and negative selection in codons 175 and 323, suggesting that the gene was influenced by selection. When the long-term effect of natural selection was explored by comparing divergence rates (K S and K N ), pv38 had a higher statistically significant K S rate than K N (Table 4) (Table 4), when intraspecific polymorphism and interspecific divergence was compared, showing P N / P S > D N /D S (p < 0.02). This result could have been the result of either a negative selection or a balanced selection [61,74]. K-test and H-test results (Table 2) and the presence of different haplotypes at intermediate frequencies ( Figure 2B) suggested that it is most probable that pv38  was influenced by balanced selection, similar to that reported for P. falciparum [73]. Such selection seemed to be domain specific. Significant values were observed for region A (p = 0.014) when intraspecific polymorphism and interspecific divergence was calculated in each region (Additional file 3), this being where most of the substitutions found became accumulated, whilst neutrality could not be ruled out for region B (p = 0.1). Functional/structural constraint due to the presence of an s48/45 domain was also probable for pv38, given this region's low diversity, two negatively selected sites and a statistically significant K S > K N .

Linkage disequilibrium (LD) and recombination
Several statistics were calculated for determining possible associations between polymorphisms and/or the presence of recombination in pv38. Z nS did not reveal statistically significant values, indicating that pv38 polymorphisms were not associated. Lineal regression between linkage disequilibrium (LD) and nucleotide distance revealed a reduction in LD as nucleotide distance increased, indicating that intragenic recombination might have led to new variations being produced. The ZZ statistic was calculated to confirm whether recombination affected pv38 evolution, showing no significant values (Table 2); however, 2 RM (minimum recombination events) were found. The GARD method (in Datamonkey web server) gave a recombination breakpoint in position 524. Prior studies have suggested that new haplotypes could be produced through recombination in spite of functional constraints [73]. Intragenic recombination could thus be one of the factors promoting diversity in the pv38 gene. Crosslinking during recombination could produce new combinations between the gene's 5′ (region A) and 3′ region (region B) as the breakpoint found in this gene was located upstream of the region encoding the s48/45 domain (region B). As only one polymorphic site was found in pv12, the aforementioned tests were not carried for this gene.
pv12 and pv38 should be considered for an antimalarial vaccine The lack of a totally effective vaccine against human malarial parasites is at least partly due to high genetic diversity found in proteins involved in red blood cell invasion. These molecules' constant exposure to the host's immune system allows the fixation of mutations generating an adaptive advantage preventing their recognition. Antigens such as pvmsp-1, pvdbp, pvmsp-3α, pvmsp-5, pvmsp-7C, pvmsp-7H, pvmsp-7I and pvama-1 have shown high genetic diversity which appears to be maintained by positive-balancing selection [10][11][12][13][14][15][75][76][77][78]; however, other antigens are highly conserved despite being exposed to the host's immune system. Surface antigens such as pvmsp-4, pvmsp-7A, pvmsp-7 K, pvmsp-8, pvmsp-10, pv230 or others in the rhoptries (pvrap-1 and pvrap-2) appear to evolve more slowly due to a possible functional constraint in their encoded proteins [70,71,[79][80][81][82]. Thus, most mutations have become eliminated from the population, maintaining a conserved protein structure, even throughout these parasites' evolutionary history [70,71]. The latter behaviour seems to have been directing pv12 and pv38 evolution, highlighting high conservation at both intra-and inter-species level due to the influence of negative selection exerted on s48/ 45 domains which are important for red blood cell recognition [30]. Although antigens having low genetic diversity are usually not immunogenic [83] nor do they induce protection-inducing responses [84], some limited polymorphism antigens have been shown to be able to induce immunogenicity and protection [85]. Therefore, pv12 and pv38 (or their s48/45 domains) should be evaluated regarding vaccine development because immune responses against 6-Cys family antigens appear to be directed against structural epitopes in s48/45 domains [86][87][88], blocking such domains should prevent invasion [30,88] and being highly conserved and having a functional constraint, allele-specific immune responses are thus avoided.

Conclusions
The p12 and p38 genes in P. vivax were seen to have low genetic diversity; the regions encoding the s48/45 domains seemed to be functionally or structurally constrained. Several members of the 6-Cys family are found on the surface of malaria parasites in every stage [28,[36][37][38][39]69] and some of them (e g, P48/45, P230) are considered to be promising (transmission-blocking) vaccine candidates [36,37,87]. Epitopes identified by monoclonal antibodies against this type of protein are structural and have been localized within s48/45 domains [86,87] which seem to be involved in host-pathogen interaction [30,72]. Since pv12 and pv38 share structural characteristics with members of the 6-Cys family, added to their antigenic characteristics [38][39][40] and the low genetic diversity found in this study, the proteins encoded by these genes or their functionally/structurally constrained (conserved) regions could be born in mind when designing a multistage, multi-antigen subunit-based anti-malarial vaccine.

Additional files
Additional file 1: pv12 and pv38 haplotypes distribution in the Colombian population. Haplotype distribution found in pv12 (A) and pv38 (B) from 2007 to 2010.
Additional file 2: Synonymous substitution per synonymous site rate (d S ) and non-synonymous substitution per non-synonymous site rate (d N ) in s48/45 domains from pv12 and pv38 genes. No statistically significant differences were found by codon-based Z-test or Fisher's exact tests. se: Standard error. pv12 s48/45 domain in region A: nucleotides 82-471; pv12 s48/45 domain in region B: nucleotides 589-906; pv38 s48/45 domain in region B: nucleotides 481-852 -: There is no s48/ 45 domain in the pv38 region. Numbering is based on the Sal-I reference sequence.
Additional file 3: McDonald-Kreitman test for evaluating the action of natural selection in pv12 and pv38 gene regions A and B. The McDonald-Kreitman test was done using sequences obtained from databases (worldwide isolates) together with Colombian ones, and just with those obtained in the Colombian population. The interspecies divergence data was obtained from comparing Plasmodium vivax sequences with two related species: Plasmodium cynomolgi and Plasmodium knowlesi. Significant values are underlined. pv12: region A, nucleotides 1-546 and region B, nucleotides 547-1,095. pv38: region A, nucleotides 1-459 and region B, nucleotides 460-1,065.