Nineteen new sequences of dbp regions III–VI were obtained

A single 1388 bp band was amplified. Nineteen different sequences are being reported here amongst the 23 Colombian isolates analysed, corresponding to Sal-1 strain sequence amino acid positions 582–1043 (Regions III to VI). Only one type of sequence was found for each isolate. Four sequences were identical to Sal-1 dpb, while the remaining 19 were geographically grouped and consecutively numbered within each region as follows: COLns01-COLns08 (Norte de Santander) (GenBank accessions AY341900-07), COLlla01-COLlla11 (Villavicencio, Llanos Orientales) (GenBank accessions AY341889-99). Nucleotide (ALIGN_000583) and amino acid (ALIGN_000584) alignments were reported in the EMBLALIGN database (http://www.ebi.ac.uk) (Lombard et al. 2002).

Region I comprised residues 1 to 104; region II, residues 105 to 429; region III, 430 to 595; region IV, 596 to 707; region V, 708 to 808; and region VI, residues 809 to 901.

Fig. 2 shows the amino acid alignment for P. vivax DBP regions III–VI including new Colombian isolates. One sequence (COLns08, from Norte de Santander, Colombia) had a 27 pb (9 amino acid – FAESTKSAE) insertion in the middle of the central hydrophilic region (at the beginning of region IV); this insertion has been previously reported in some Papua New Guinean (Tsuboi et al. 1994) and Korean (Suh et al. 2001) dbp sequences. The inserted sequence is a tandem repeat created by partial duplication of the sequence immediately following the insertion site and is not repeated in any other part of the gene (Tsuboi et al. 1994).

Fig. 2. Colombian isolate DBP region III–VI amino acid alignment. Grey shading indicates specific high binding affinity peptides (1640 and 1642); bold represents co-variant residues with Korean and Papua New Guinean sequences; italics show insertion in COLns08 sample. Dots indicate conserved residues and dashes represent gaps introduced for alignment.

Cluster and recombination analysis of P. vivax DBP regions

Sequences from Papua New Guinea, Salvador, Korea and Colombia were employed for cluster analysis. These sequences were divided into regions I–VI, according to the region delimitation previously reported by Fang (Fang et al. 1991) (Fig. 1B).

Fig. 3 shows the amino acid cluster trees constructed for each dbp region by the Neighbour-Joining method with 1000 bootstrap replicates. Maximum Parsimony and Minimum-Evolution methods also resulted in similar tree topology (not shown); nucleotides and amino acids show equivalent results. Non-significant bootstrap values were observed in all trees. This is explained by the scarcity of variant sites and the scattered pattern of the variation shown among dbp sequences (the majority of variant sites were singletons, see Table 1). These low bootstrap values are common to other studies with P. vivax antigens: AMA and MSP1 genes (Figtree et al. 2000), TRAP (Putaporntip et al. 2001), MSP3α (Rayner et al. 2002) and MSP1 (Putaporntip et al. 2002).

Fig. 3. Neighbour-Joining tree of amino acid sequences from Plasmodium vivax DBP regions (I–VI) from endemic populations based on pair-wise differences using P distance. The significance of the nodes was proved by Bootstrap calculated with 1000 replicates. The scale indicates 0·01 substitutions per site. (A) Tree for region I. (B) Region II. (C) Region III. (D) Region IV. (E) Region V. (F) Region VI. Colombia, COL; Norte de Santander, ns; Llanos Orientales, lla; Korea, KOR; Papua Nueva Guinea, PNG; Wosera, W; Madang, m.

Table 1. Diversity analysis of dbp gene (n, number of sequences sampled; Sites, number of sites included; V.S., number of variant positions; S, number of singleton positions, P.I., number of parsimony informative sites; π, average number of nucleotide substitutions in DNAsp 3·51; C, conservative change in amino acid; NC, non-conservative change in amino acid; d, average number of sequence diversity for amino acid with p-distance.)

No clear geographical clustering was found amongst the Colombian samples; this might indicate constant intermingling of P. vivax in eastern Colombia (Fig. 2). Colombian and Salvador sequences frequently clustered together; moreover, 4 Colombian sequences found were identical to Sal-1 sequence. These observations indicate a close relationship between Central American and Colombian strains. The sequence from Belem strain (Brazil), only available for regions I and II, clustered with Villavicencio strains (Fig. 3A, B). These observations taken together suggest continuous dispersion of strains throughout the whole of the central and Southern American regions via Colombia, constituting a bridge between these two regions.

A geographical clustering pattern was only clearly observed in region VI (Fig. 3F). Region II topology (Fig. 3B) was similar to that reported by Cole-Tobian & King (2003), showing a pattern having a low degree of mix between sequences from Papuan and American strains. Regions I, III, IV and V show more mixing among strains from different geographical locations (Fig. 3A, C, D and E respectively).

This pattern (i.e. protein regions showing geographic clustering and others that do not) highlights the modular separation between protein regions, indicating the different extent of host specialization in these modules. This tendency was shared with other P. vivax antigens such as MSP-1 (Putaporntip et al. 2002) and MSP-3α (Rayner et al. 2002). The phenomena originating these modules could also be based on genetic drift, recombination processes and different selective immune pressures.

Sequences with a characteristic insertion at the end of region III do not show a unique cluster, indicating that there are characteristics in the rest of region III sequences driving the clustering pattern. This sector may be considered as evidence of a recombination process, or may represent a convergence event between Papuan and Colombian sequences. Evidence of recombination was found by using DnaSP 3.51 software (Rozas & Rozas, 1999). There were 9 recombination events (Rm) in region II between nucleotides (667, 676) (676, 707) (716, 718) (718, 722) (734, 737) (835, 847) (903, 904) (905, 1073) (1073, 1096) and 7 Rm in regions III–VI (1426, 1525) (1550, 1714) (1714, 1738) (1898, 2068) (2092, 2321) (2321, 2497) (2500, 2678). Only region I did not present recombination sites. Regions II and III displayed the highest recombination events in the dbp gene. This agrees with the recombination rate observed in P. vivax (Putaporntip et al. 2002; Rayner et al. 2002).

Geographical clustering has been reported in other P. vivax genes. For example, some evidence has been shown in the circumsporozoite protein (CSP) (Escalante et al. 2002; Qari et al. 1992) concerning the existence of strains having a specific geographical location. Polymorphic variants have been found in tropical regions of Brazil, Southeast Asia and Papua New Guinea (Escalante et al. 2002; Qari et al. 1992), but strikingly different strains have been described in China and Korea (Kim et al. 2002; Mann et al. 1994). This wide distribution of strains can be explained by the world's population mobility favouring the spread of different parasite strains.

Diversity analysis of P. vivax DBP sequences

Table 1 shows the calculations for the different diversity parameters used, sorted by region studied and geographical location. Different diversity patterns were observed between the geographical regions analysed.

Colombian sequences showed higher nucleotide diversity in region I than Papuan and Korean samples; on the other hand, Colombian isolates showed lower nucleotide diversity in all the other regions with respect to Papuan samples, displaying the lowest nucleotide diversity values for regions II and III with respect to Papuan and Korean samples. Korean samples showed higher nucleotide diversity in region II, but the lowest nucleotide diversity in region IV with respect to Papuan and Colombian samples. Altogether, Papuan samples were more diverse at the nucleotide level than Colombian and Korean samples. The diversity pattern was more similar between Papuan and Korean samples than Colombian samples. Sequence diversity at the protein level showed a similar pattern; the only difference was that Colombian samples showed higher sequence diversity in region VI. Calculations from all sets of sequences showed that dbp regions II and VI were the most variable; regions IV and V were the least variable. Region II was divided into portions described as being crucial for erythrocyte binding (II-1, II-2, II-3) (Ranjan & Chitnis, 1999; Singh et al. 2003). The highest diversity exhibited inside region II was observed within segment II-2 (comprising Cys5 to Cys8) (data not shown); this agrees with previous studies where the specific erythrocyte binding site displays higher diversity (Ampudia et al. 1996; Tsuboi et al. 1994; Xainli, Adams & King, 2000).

All changes observed in region I corresponded to parsimonious positions, when all sets of sequences are considered, or when sorted by geographical region. Several variant positions showed co-variation, exhibiting the same amino acid substitutions (ALIGN_000584). This might indicate that the changes in this region could be the result of similar selective pressure. Moreover, most of these changes were conservative. This pattern of higher conservative than non-conservative substitutions was shared by all regions in dbp.

Papuan and Korean samples showed the highest proportion of variant sites in Region II; however, the nature of the variation differed amongst them. Papuan samples showed a higher number of singletons than parsimony positions, whereas Korean samples showed the opposite tendency. Colombian samples showed the lowest proportion of variant sites in this region, also showing more parsimonious than singleton sites, like the Korean samples. A high proportion of singletons indicated a disperse substitution pattern, possibly searching for greater sequence diversity (for immune evasion or intra-specific competence). Thus, diversity generation may play a more major role in Papuan samples than in Colombian and Korean samples. On the other hand a more localized variation pattern (estimated by the proportion of parsimonious sites) might be a response to a different type of pressure, acting in a more focused way on specific sectors of the molecule, as observed in Colombian and Korean samples.

Region III showed the lowest difference in variant positions amongst geographical samples; this homogeneity extended to singletons and parsimonious positions; more singleton positions were detected than parsimonious. This region's variability was noticeable for Papuan and Korean samples, both showing the same behaviour as Region II in Papuan samples. Papuan and Colombian Regions IV and V showed a high proportion of singletons and low parsimonious positions. This behaviour was also similar to that observed for region II in Papuan samples, but with less sequence variability in Colombian samples.

Region VI showed different behaviour; the number of singletons and parsimonious informative positions was very similar; Colombian and Papuan samples showed very similar values for variant sites and sequence diversity at nucleotide level, but at amino acid level, Colombian sample values were higher. This indicates that Colombian samples have a more intense process of protein diversification than Papuan samples, although nucleotide diversity is equivalent in both geographical regions.

The overall calculations of the variant sites for each protein region show that their values were always higher than the values from each geographical location in every case. This indicates that the majority of changes were present in positions that did not co-vary and each isolate had a unique mode of variation.

In total, 27 co-variant residues among Colombian, PNG and Korean samples were found when compared to Sal-1 (ALIGN_000584: amino acid position 46, in region I; positions 149, 226, 239–41, 245, 272, 279, 292 and 358 in region II; 476, 509, 572 and 586–94 in region III; 774 in region V; and 834, 888 and 893 in region VI).

The diversity pattern was consistent with observations regarding transmission rates seen in the different geographical regions. In other words, Papua regions suffer stronger rates of transmission and prevalence of P. vivax infection than Colombia, and this behaviour was related to high sequence diversity (Cui et al. 2003). This means that Colombia has low transmission rates and, therefore, Colombian samples show lower sequence diversity than Papuan isolates.

Tendency towards a higher number of transversion changes rather than transitions, previously described in genes encoding P. falciparum membrane proteins due to the A+T composition bias, was observed in all P. vivax dbp gene regions; transition value was 4·97 and transversion value 5·39. This tendency was maintained throughout geographical locations for all regions. Although transitional substitutions are more common as transversions (Collins & Jukes, 1994), P. vivax genome exhibiting a strong codon bias may affect the transition/transversion ratio and estimates of synonymous and non-synonymous substitutions (Escalante, Lal & Ayala, 1998).

The dbp gene has lower nucleotide diversity values (π=0·01224) than other P. vivax surface proteins, such as AMA1 (π=0·0174), MSP1 (π=0·0451) (Figtree et al. 2000), but higher than TRAP (π=0·00594) (Putaporntip et al. 2001). This may reflect a major functional restriction in DBP protein, due to its crucial role in the invasion process.

These values are considerably lower than those for P. falciparum surface proteins like MSP-1 (π=0·08792), MSP-2 (π=0·04409) and AMA-1 (π=0·01635) (Escalante et al. 1998). P. vivax dbp is also more diverse than EBA-175 (π=0·00366), its homologue in P. falciparum (Escalante et al. 1998), a value that might explain P. vivax's dependence on genetic variability for immune response evasion, as opposed to P. falciparum's reliance on alternative invasion pathways. The differences in the diversity pattern between P. vivax dbp and other non-homologous proteins in P. falciparum might reflect differences in invasion roles, virulence and number of clonal generations.

Evidence of selection in P. vivax DBP regions

Several tests were performed to determine the type of selective pressure on dbp: Nei–Gojobori's (Nei & Gojobori, 1986), and Tajima's (Tajima, 1989). Nei–Gojobori's test was computed by considering either all positions or just variable positions. All positions were used for Tajima's test. All dbp sequences used were those reported in the GenBank database, together with our dbp sequences for regions III–VI. Table 2 shows the analysis for all the studied regions (I–VI). Fu-Li's (Fu & Li, 1993) test was also performed, producing similar results to Tajima's test (data not shown).

Table 2. Nei–Gojobori's and Tajimas's tests for dpb gene (n, number of sequences sampled; Sites, number of sites included; dS, synonymous substitution; dN, non-synonymous substitution. For dN–dS in Nei–Gojobori's test and Tajima's test, positive values correspond to positive selection pressure, and negative values correspond to negative selection pressure. Bold represents significant values P<0·05. Italic and underlined letters represent significant values P<0·01.)

Tajima and Nei–Gojobori's tests showed distinct tendencies. Nei–Gojobori's test for all positions showed mainly positive selective pressure both when each geographical region was considered individually and also when all sequences were analysed altogether; this positive selection was statistically significant only in region II, except for Colombian samples (data not shown). This tendency to present more non-synonymous than synonymous substitutions has been previously observed in other reports (Ampudia et al. 1996; Tsuboi et al. 1994; Xainli, Adams & King, 2000). If only variant positions are used, more regions than just region II show positive selection for all sets of sequences (Table 2). Tajima's test showed mainly negative selection pressure in the dbp regions considered, this being significant for all sets of sequences in regions II, IV and V. Cole-Tobian et al. (2003) showed the same tendencies for region II, where opposite results were found according to the test used (Table 2). The contrasting detection of both negative and positive selection pressures in several of the molecule's regions (especially in regions II and VI) may imply a balance between a high functional restriction for maintaining structural constraint due to DBP's key role in the erythrocyte invasion process and the search for diversity as a response against host immune system pressure (as occurs in other Plasmodium surface antigens (Escalante et al. 1998; Putaporntip et al. 2001; Verra & Hughes, 2000)).

The differences between the selective pressures in samples from geographical regions showed that Colombian isolates possess a diversifying selection in more molecule regions than Papuan and Korean isolates. This might reflect a true variation between isolates, in which case, Colombian parasites suffer a stronger selective stress than Papuan and Korean parasites. Although these differences between geographical samples were also sustained by the patterns of diversity observed, these results might reflect a sample size effect and the method used for isolating dbp genes (PCR or DNA library). More sampling will be necessary for clarifying this issue.

Analysis of the 3′ cysteine region

A special approach was used for region VI to define conserved residues and physicochemical properties, due to this region's particular characteristics suggesting its possible relevance in vaccine design: the interspecies conserved cysteine arrangement; the presence of one high-activity reticulocyte binding peptide (1660) (Fig. 2) (Ocampo et al. 2002) and its diversifying and selection pattern, suggesting its exposure to the host immune system.

A clear conservation tendency could be seen in the general alignment of Colombian samples, particularly in cysteines and aromatic amino acids in region VI (Fig. 2).

Region VI presented the highest geographical diversity (Table 1), although it possessed the lowest number of variable sites, suggesting that this region was highly conserved. This conservation pattern might indicate a key role for these residues in maintaining protein tertiary structure or other undetermined functions (Michon et al. 2002).

Colombian and Papua New Guinean sequences were chosen for specifically determining which residues were conserved and some of their physicochemical properties. EBL family proteins from different species were also included in this analysis (see Materials and Methods section).

Fig. 4 shows region VI sequence's protein alignment. High conservation levels (around 67% or higher) were observed, as well as the previously described 8 cysteine residues and hydrophobic amino acids, like tyrosine, leucine, phenylalanine and isoleucine (Adams et al. 1992).

Fig. 4. Conserved carboxyl cysteine-rich alignment of the erythrocyte binding proteins. The shading indicates levels of identity among sequences. White letters on a black background represent amino acid residues that were 100% identical; white letters on a grey background represent amino acid residues that were 80% identical; and black letters on a grey background represent amino acid residues that were 60% identical. Bold letters represent conserved physicochemical properties: Small, S; Polar, P; Hydrophobicity, H; Charged, C; Aromatic, A; P. vivax Duffy Binding Protein, PVDBP; Madang Papua New Guinea, PNGm; Colombia, COL; Norte de Santander, ns; Llanos Orientales, lla; erythrocyte binding protein, EBP; erythrocyte binding antigen, EBA; erythrocyte binding like, EBL; P. cynomolgy, Pc; P. knowlesi, Pk; P. vivax, Pv; P. falciparum, Pf; P. yoelii, Py; P. reichenowi, Pr; P. berghei, Pb; P. chabaudi, Pch; P. vinckei, Pvin. Dashes represent gaps introduced to maintain alignment.

Five physicochemical properties proved to be conserved following analysis with GeneDoc software. Size was the most conserved factor (17/70 residues), followed by polarity (14/70), hydrophobicity (5/70) and aromatic residues (6/70); overall, physicochemical properties were conserved in 68% of total residues.

Kyte–Doolittle hydrophobicity variation analysis was done for every position in region VI (Fig. 5) to see the effect of amino acid substitution on the region's hydrophobic characteristics (Kyte & Doolittle, 1982). The tendency towards maintaining hydrophobic characteristics that could be observed in this region was probably related to its function.

Fig. 5. Kyte–Doolitle Hydropathy Plot of region VI (3 cys) for variant sequences from Papua New Guinea and Colombia. Score corresponds to normalized (to 1) values used in the computation. Position corresponds to amino acid numeration according to Sal-1 sequence (Fang et al. 1991).

The results from hydrophobicity analysis of variant sequences led to the conclusion that region VI was a diverse region as presented in Table 1, having a high number of non-synonymous substitutions that did not alter its general hydrophobic characteristics.

The presence of conserved hydrophobic residues makes this group's recombinant expression very difficult due to its insolubility (Fraser et al. 1997), a fact which has paralysed functional assays aimed at describing region VI's role in the P. vivax invasion process.

The role of P. falciparum EBA-175 recombinant protein excluding the carboxy-terminal cysteine-rich region (C-cys), transmembranal and cytoplasmic domains, has been indirectly evaluated in preceding studies. These studies determined that the whole protein is not necessary during the invasion process or for parasite growth, but that this domain might be involved in intracellular mobilization processes and the use of alternative invasion pathways. Such domains are probably indispensable for the protein's structural conformation, facilitating its binding to erythrocytes. This preliminary evidence suggests that these C-cys domains are involved in the invasion process, not directly related to adhesion but being capable of transmitting or receiving signals from the interaction of the binding regions (EBL) with receptors on the target cell (Reed et al. 2000).

The presence of DBL and C-cys domains in several different Plasmodium species' molecules suggests that they are performing a diverse role in receptor-ligand interactions. The residues conserved within these domains can be useful in identifying crucial determinants for cell–cell interaction during malaria.

Concluding remarks

Three forces in conflict (functional restriction, adaptive response to immune pressure, and contingency phenomena), represented in processes like diversifying selection, purifying selection, recombination and functional restriction, drive the evolution of DBP. These forces operate with different intensity depending on the considered molecule regions. Thus, two main region classes (or modules) can be delimited. The first type showed relatively low variation amongst sequence groups, showing mainly purifying selection and did not display a specific geographical clustering pattern (regions III, IV, V). The second type showed relatively high variation amongst sequence groups, diversifying selection when just variant positions were considered and purifying selection when all positions were analysed; it also displayed a geographical clustering pattern (regions II and VI belong to this second type).

Regions II and VI were associated with receptor-binding functions, a key role in the invasion process. This has been proved in region II and multiple evidence has implied a similar function in region VI. Simultaneously, exposure to the host's immune system causes selective pressure leading to an evasive response. The specific nature of this adaptation to the host drives these sectors' evolution towards geographical clustering patterns; these patterns are absent in other molecule sectors where low exposure to the host immune system and/or a strongest functional restriction (like structural constraints) produce high conservation.

Multiple and opposite selective forces acting in a relatively small region of a protein have been previously described in HLA (Suzuki & Gojobori, 1999), indicating that selective forces act on specific residues rather than whole regions; similar behaviour could be occurring within dpb. Although this might imply the existence of different functional subdomains within these sectors, more information is necessary to clarify this issue.

Recombination forces appear to be acting in all DBP and, therefore, their influence may be less important. Recognizing these different modules and understanding the factors that drive their evolution will allow the design of strategies for elucidating host–parasite interactions, in turn leading to new vaccine targets.