ABSTRACT
Mutations in the Plasmodium falciparum k13 (Pfk13) gene are linked to delayed parasite clearance in response to artemisinin-based combination therapies (ACTs) in Southeast Asia. To explore the evolutionary rate and constraints acting on this gene, k13 orthologs from species sharing a recent common ancestor with P. falciparum and Plasmodium vivax were analyzed. These comparative studies were followed by genetic polymorphism analyses within P. falciparum using 982 complete Pfk13 sequences from public databases and new data obtained by next-generation sequencing from African and Haitian isolates. Although k13 orthologs evolve at heterogeneous rates, the gene was conserved across the genus, with only synonymous substitutions being found at residues where mutations linked to the delayed parasite clearance phenotype have been reported. This suggests that those residues were under constraint from undergoing nonsynonymous changes during evolution of the genus. No fixed nonsynonymous differences were found between Pfk13 and its orthologs in closely related species found in African apes. This indicates that all nonsynonymous substitutions currently found in Pfk13 are younger than the time of divergence between P. falciparum and its closely related species. At the population level, no mutations linked to delayed parasite clearance were found in our samples from Africa and Haiti. However, there is a high number of single Pfk13 mutations segregating in P. falciparum populations, and two predominant alleles are distributed worldwide. This pattern is discussed in terms of how changes in the efficacy of natural selection, affected by population expansion, may have allowed for the emergence of mutations tolerant to ACTs.
INTRODUCTION
Antimalarial resistance is one of the greatest threats to malaria control and its eventual elimination (1, 2). The emergence of a delayed parasite clearance phenotype following artemisinin-based therapies in Plasmodium falciparum was first reported in 2007 and 2008 in western Cambodia and was later reported in Myanmar, Thailand, Vietnam, Lao People’s Democratic Republic (Lao PDR), and South China (2–7). The occurrence of this phenotype has raised concerns that it may spread to Africa or lead to resistance, following the previous path of chloroquine (CQ) and sulfadoxine-pyrimethamine (SP) resistance (8–9). However, there is currently no evidence that this phenotype has spread outside of Southeast Asia (SEA), nor is there evidence of the emergence of resistant parasites (2).
Artemisinin delayed parasite clearance was first defined and confirmed using in vitro ring-stage survival assays (RSAs) (1, 10, 11). A series of genome-wide association studies have shown that this phenotype is associated with nonsynonymous mutations in the propeller domain of a P. falciparum kelch gene (codons 441 to 726) found on chromosome 13 (Pfk13) (12, 13). To date, 14 independent Pfk13 mutations have been associated with clinical delayed parasite clearance (2), and of these, only 5 have been validated by in vivo and in vitro experiments (2): N458Y, Y493H, R539T, I543T, and C580Y. Furthermore, there is evidence of the independent emergence of the C580Y lineage across SEA (14, 15). The most frequent and widely detected Pfk13 mutations found in Cambodia, Vietnam, and Lao PDR are Y493H, R539T, I543T, P553L, and C580Y (5, 7). The mutations F446L, P553L, N458Y, R561H, P574L, and C580Y are the most frequently observed in western Thailand, Myanmar, and China (7, 16–18). Only two mutations, P553L and C580Y, are found across SEA, with the latter mutation being more widespread in eastern Thailand (16–18).
In contrast, a high number of single nonsynonymous mutations, including some associated with delayed parasite clearance, have been found at a low frequency in Africa (19). However, many of the mutations reported in Africa have not expanded to the local parasite populations (2, 19–23). In the case of South America, only Pfk13 C580Y mutant parasites have been found in Guyana, with the flanking microsatellite profiles being different from those observed in SEA, suggesting an independent origin (24). Studies in other countries in South America have searched for polymorphisms in the Pfk13 gene, but no evidence of these mutations has been found to date outside Guyana (25, 26).
Although our understanding of which specific mutations within the PfK13 domain are associated with the delayed parasite clearance phenotype is still incomplete (especially in the context of the genetic background of parasites outside of SEA), the identification of polymorphisms in the Pfk13 gene has been used as a valuable genetic marker for the rapid detection, tracking, and monitoring of this phenotype. Here, to better understand the evolution of the k13 gene, its divergence across the genus Plasmodium and its genetic variation in P. falciparum populations were studied. In particular, we explored the mode and rate of evolution of the k13 gene by analyzing its orthologs across the Plasmodium genus and 982 Pfk13 sequences, including those in public databases, along with new data obtained by next-generation sequencing of African and Haitian isolates.
RESULTS
Phylogenetic-based analyses of the k13 gene.This investigation reports 10 new k13 orthologous gene sequences, obtained by nested PCR, from macaque and gibbon parasite strains and chimpanzees infected with Plasmodium reichenowi, Plasmodium billcollinsi, and Plasmodium gaboni. The sequences obtained here were aligned with those available in public databases (see Materials and Methods and Data Set S1 in the supplemental material) (27). The alignment constructed with the Plasmodium species orthologous k13 genes showed that the encoded protein is remarkably conserved across the genus. Estimates of the evolutionary divergence among orthologous k13 gene sequences were low at both the nucleotide and amino acid levels compared to that of other genes (28). The k13 gene was almost identical among species belonging to the subgenus Laverania. For example, the genetic distance between P. falciparum and Plasmodium praefalciparum from gorillas was 0.00, and that between P. falciparum and P. reichenowi from chimpanzees was only 0.01 ± 0.00 (Data Sets S2 and S3). Furthermore, all k13 orthologous genes displayed a high degree of conservation in the kelch propeller domain among Plasmodium spp. (Fig. 1). In particular, the number of base substitutions per site between the two major malarial parasite k13 orthologs (Pfk13 and Plasmodium vivax k13 [Pvk13]) was low (0.09 ± 0.01) (Data Set S2) and showed synonymous substitutions between P. falciparum and P. vivax, with the exception of V568I versus V568G in Pfk13 and Pvk13, respectively (Fig. 1; Data Sets S1 and S4).
Alignment of the Plasmodium K13 kelch propeller domain. The Pfk13 mutations that have been associated with clinical delayed parasite clearance (2) are shown. The five mutations validated in vivo and in vitro are indicated with an asterisk. The Haemoproteus tartakovskyi (an avian haemosporidian parasite) k13 gene was included as an outgroup.
To compare all Plasmodium species k13 orthologous genes considered in this study (including different strains of some parasites), a Bayesian phylogeny was estimated using the nucleotide sequence alignment, given the high conservation at the protein level (Fig. 2). As expected, the phylogeny was consistent with earlier reports using putatively neutral genes, genes encoding merozoite antigens, and/or mitochondrial genomes (27–29). Overall, the Laverania subgenus, which includes P. falciparum, formed a monophyletic group, as did P. vivax and its closely related nonhuman primate malaria parasites from SEA and all rodent malaria parasites. All mammal Plasmodium species shared a common ancestor with avian Plasmodium species (Fig. 2).
Bayesian phylogenetic hypothesis of Plasmodium parasites based on the k13 gene (42 sequences; 1,985 bp, excluding gaps). The values at the nodes are posterior probabilities (see Materials and Methods). Haemoproteus tartakovskyi (an avian haemosporidian parasite) k13 gene was included as an outgroup.
Then, using the Bayesian tree topology, the relative evolutionary rates of the k13 gene were estimated in a maximum likelihood framework with the RelTime method (30). This revealed a pattern in which there was a faster k13 evolutionary rate in malaria parasites infecting Southeast Asian macaques and their relatives, including the human malaria parasite P. vivax (Fig. 3). The k13 relative evolutionary rates were approximately 12 times higher in malaria parasites infecting Southeast Asian macaques and their relatives than in parasites from other mammalian hosts and African ape parasites. Differences in rates have also been observed in other genes in this clade, which includes the human parasite P. vivax (31, 32). Interestingly, Pfk13 also showed a k13 relative rate of evolution faster than and similar to that in SEA parasites (Fig. 3). The acceleration in the P. falciparum k13 evolutionary rate, however, did not translate into the gene, being under diversifying natural selection when P. falciparum diverged from its closest related species (see below).
Relative evolutionary rates of the Plasmodium k13 gene. Branches are colored from high (green) to low (purple) according to their relative rates with respect to the rate at the root, which was set to 1, as estimated by RelTime without calibration constraints.
The relative evolutionary rate differences in the k13 gene between clades were further examined by estimating synonymous and nonsynonymous substitution rates per codon per million years (My) using a codon model implemented in the CodonRates (v1.0) program (33). In this analysis, the k13 nonsynonymous substitution rate (0.022 ± 0.012 per codon per My) was slightly higher (on average) than the synonymous substitution rate (0.012 ± 0.006 per codon per My) across the genus. The concordance statistic (S) showed that synonymous and nonsynonymous substitution rate changes over time were correlated (P < 0.5), suggesting that there was no pattern indicative of positive selection acting on this gene. The branch leading to P. falciparum (27, 28) from the common ancestor with its closest sister taxon (Plasmodium praefalciparum) showed a synonymous substitution rate of 0.005 per codon per My, while the nonsynonymous substitution rate was 0.001 per codon per My.
To further assess the role that selection may have had during the divergence of the k13 gene, we utilized methods implemented in the HyPhy package (34) (see Materials and Methods). A test that examines each individual site independently for evidence of episodic diversification (MEME [35]) found no evidence of episodic positive/diversifying selection acting on any site of the gene. In addition, a branch site-unrestricted statistical test for episodic diversification (BUSTED) (36), found no evidence (likelihood ratio test, P = 0.997, which is ≥0.05) of gene-wide episodic diversifying selection in any of the branches of the k13 phylogeny. Finally, RELAX (37), a test for identifying whether the strength of natural selection (positive or negative) has been relaxed or intensified along a specified set of test branches (e.g., P. falciparum), was not significant (K = 0.74, P = 0.665, likelihood ratio = 0.19). These results showed no evidence of positive diversifying selection acting on the k13 gene either during the evolution of the sampled Plasmodium species or in the lineages of P. falciparum or P. vivax.
Plasmodium falciparum population analyses.Using an alignment with 982 complete k13 gene sequences from around the globe (519 from public databases and 463 reported in this study; see Materials and Methods and Data Set S5), low genetic diversity was found in the Pfk13 gene, as measured by π (0.0005 ± 0.0; Table 1). In contrast, the polymorphism measured using the number of polymorphic (segregating) sites had an S value of 146. The number of singleton variable sites, considering two variants, was 108. It is worth noting that no mutations associated with the delayed parasite clearance phenotype were found in our samples from Africa and Haiti. In the case of P. vivax, all 56 Pvk13 sequences available from the databases were identical.
Polymorphism found in k13 gene of Plasmodium falciparum
Evidence of natural selection acting on the Pfk13 polymorphism was explored by estimating the average pairwise number of synonymous substitutions per synonymous site (dS) and the average pairwise number of nonsynonymous substitutions per nonsynonymous site (dN). This revealed a significant (P < 0.05) excess of nonsynonymous over synonymous substitutions only in Asia, the region where mutations associated with the artemisinin delayed parasite clearance phenotype were first reported (Table 1).
A total of 147 Pfk13 haplotypes were found with a haplotype (allele) diversity (Hd) of 0.742 (standard deviation [SD] = 0.012). The uneven sampling effort did not allow us to estimate the haplotype relative frequencies per country. Overall, the median-joining network for the Pfk13 gene, estimated using Network software (Fluxus Technologies, 2011), showed two star-like shapes with at least two most predominant and common haplotypes (haplotype 1 [H1] and H2; Fig. 4A) shared by the three continents (Asia, Africa, and South America) with a broad distribution. Interestingly, one of the haplotypes (H1, n = 442, 45%) is mostly distributed in Asia, and the other (H2, n = 223, 22.7%) is mostly distributed in Africa. A third haplotype (H3, n = 50, 5%), in terms of frequency, is only found in Thailand (Fig. 4A). A noteworthy finding was that haplotype H3 contains the C580Y mutation found in the SEA region. In addition, a less frequent haplotype was found only on the China-Myanmar border (H4, n = 28, 2.9%) and had the F446I mutation, one of the most common mutations found in China and on the China-Myanmar border (Fig. 4B). However, most of the haplotypes were restricted to single countries (Fig. 4A).
Median-joining network of the P. falciparum k13 haplotypes. (A) Median-joining network of the P. falciparum k13 haplotypes using 982 cosmopolitan complete sequences. All the k13 gene sequences obtained here for P. falciparum population samples (n = 463), as well as those available in the PlasmoDB and NCBI databases (n = 519), were included. Branch lengths are proportional to divergence; node sizes are proportional to the total haplotype frequencies. The network shows 147 haplotypes found in 982 sequences with a haplotype diversity (Hd) of 0.742 (SD = 0.012). Every color corresponds to a different geographic origin. Lines separating haplotypes represent mutational steps. The four most frequent haplotypes are indicated. (B) Frequencies of candidate/validated mutations in the propeller domain of the Pfk13 gene that have been associated with clinical delayed parasite clearance by country (Table 2). Frequencies were estimated using an alignment constructed with 3,916 partial (n = 2,934) and complete (n = 982) k13 gene sequences. C-M, China-Myanmar.
Finally, using 3,912 partial and complete Pfk13 sequences (Data Set S5), we estimated the frequencies (by country) of candidate/validated mutations associated with a clinical delay in parasite clearance (2). Those mutations, found in the propeller domain of a Pfk13 gene, are shown in Fig. 4B (see also Table 2). As previously reported, mutations F446I and C580Y were the most frequent mutations in this data set and were mostly found in Asia (Fig. 4B; Table 2; Data Set S6).
Number of Pfk13 sequences with mutations (candidates/validated) associated with clinical delayed parasite clearance (per codon and per country)a
DISCUSSION
Unlike other genes with mutations associated with antimalarial drug resistance, such as P. falciparum dhfr (Pfdhfr) (38, 39), the Pfk13 gene shows several single mutations at a low frequency, and many of these have no clear associations with the delayed parasite clearance phenotype (2). However, similar to genes such as Pfdhfr (28, 38), the Pfk13 mutations associated with this phenotype appeared solely in the P. falciparum lineage, indicating that they are all of recent origin (40). Although they evolve at heterogeneous rates, as indicated by the relative rate analysis, it is evident that orthologous k13 genes are highly conserved across Plasmodium species (in several species of human, nonhuman primate, other mammalian, and avian parasites). Evidence specific to the P. falciparum lineage can be found in the low nonsynonymous substitution rate estimated in this study. This finding suggests that the k13 gene was under a strong constraint (negative selection) from undergoing natural variation during its evolution in the genus Plasmodium.
There was no evidence of fixed nonsynonymous differences between P. falciparum and its closest related Laverania species (P. praefalciparum and P. reichenowi). Together, this evidence indicates that all Pfk13 mutations, including those associated with delayed parasite clearance, are of recent occurrence. The presence of several Pfk13 mutations found at a low frequency translates into the extraordinary haplotype diversity found in this study. A naive interpretation of this pattern could find it to be contradictory to the strong conservation observed among Plasmodium species at the protein level (indicative of negative selection). However, this finding is consistent with the proposed historical population expansion of P. falciparum following expansion of the human population (41). Population expansions can increase the number of segregating mutations (even semideleterious mutations) by making negative selection less effective (42, 43). Nonetheless, as parasite populations expand, negative selection becomes more efficient in eliminating those semideleterious mutations and maintaining them at a low frequency. The spatial differential distribution of those Pfk13 mutations at a low frequency is also expected because P. falciparum is not a single population, so there are local population expansions where different mutations occur. This pattern of population expansions is summarized in Tajima’s D results, where the low value of π and a high number of segregating sites (mostly singletons) lead to a significant negative estimate for the statistic D (−2.683, P < 0.001). Thus, considering the parasite demography, the Pfk13 gene is expected to have several recent geographically distinct mutations at a low frequency, with each individual parasite lineage carrying one or very few of those segregating sites (42, 43). Although in P. vivax there is also evidence of a recent population expansion (44, 45), no segregating mutations were found in Pvk13. Given the limitation of our sample (n = 53), we cannot discern whether this is a real pattern or a sampling artifact.
The data are consistent with both selection and demography driving the dynamic of Pfk13 mutations conferring artemisinin-based combination therapy (ACT) tolerance (46). The parasite population expansion (demography) relaxes purifying (negative) selection, allowing for semideleterious mutations to segregate (41, 43, 46, 47), and the introduction of ACTs selects for mutations conferring delayed parasite clearance/tolerance out of those several variants present at a low frequency, even if such variants have a fitness cost in the absence of drug pressure (8, 46–48). Under such dynamics, as the data in Africa seem to indicate, an asymptomatic untreated reservoir will favor ACT-sensitive parasites, keeping a negative selective pressure on those semideleterious mutations conferring tolerance to the drugs (47–50). Consequently, the Pfk13 mutations present today are unlikely to increase in frequency in regions such as Africa, where their relative advantage is still overcome by their fitness cost in the absence of treatment (e.g., subclinical or asymptomatic infections at a high frequency [49, 50]), and parasites with those mutatiosn are still susceptible to aggressive/monitored drug treatment (8, 46). The interplay between the putative fitness cost of these mutations and the absence of drug pressure in asymptomatic individuals could be facilitated by the higher population sizes of P. falciparum in Africa, making differences in selection regimens prevail over genetic drift in the dynamics of these mutations (46, 48). Under the current scenario, interventions directed to reduce the asymptomatic reservoir, such as targeted mass drug administration (TMDA) with ACTs, will remain effective because there are not yet Pfk13 mutations conferring complete resistance (46, 51).
Although the observed patterns are robust, it is worth noting that many studies in P. falciparum directly report the frequency of mutations related to ACT-delayed parasite clearance. This practice not only may lead to underreporting of mutations occurring at a low frequency but also limits our capacity to track geographic changes in haplotype frequencies or infer other patterns. As an example, the frequency of the C580Y mutation is high in western Cambodia (52, 53), but the deposited sequence data are scarce. Whether there is an expansion of a particular haplotype in the parasite population cannot be addressed across studies by using mutation frequency data.
The proposed dynamics explain why the actions deployed thus far have limited the expansion of these mutations from Asia, and there has been only a single haplotype of C580Y identified in Guyana. Whether the spread will remain limited (17, 24) is a matter that requires continued monitoring. Importantly, reporting novel mutations is critical since they may change the parasite fitness under drug pressure. Thus, considering that all observed mutations in Pfk13 are of recent origin, understanding their evolutionary dynamics at the early stages of drug resistance could further inform containment strategies wherever putative resistance mutations start to increase in frequency.
MATERIALS AND METHODS
Samples/sequences from the genus Plasmodium.Orthologous k13 genes were sequenced from macaque and gibbon parasite laboratory strains provided by the Centers for Disease Control and Prevention (CDC): Plasmodium inui (OS strain), Plasmodium hylobati, Plasmodium fieldi (N-3 strain), Plasmodium simiovale, and Plasmodium gonderi. Information on these species and strains can be found elsewhere (54). Archived remainders of blood samples from chimpanzees infected with Plasmodium reichenowi (Shegue and 20 isolates [28]), Plasmodium billcollinsi (the Ithaito isolate [28]), and Plasmodium gaboni (Shegue and 39 isolates [28]) were also used in this study (see Data Set S1 in the supplementary material). The chimpanzees were housed at the Jane Goodall Institute (JGI) Tchimpounga Chimpanzee Rehabilitation Center in the Democratic Republic of Congo (DRC), and their blood samples were collected throughout 2009 and 2010 by the veterinary staff as part of the chimpanzees’ routine health examinations, following standards approved by the Pan Africa Sanctuary Alliance. More information on these species and ethical guidelines can be found in reference 28.
In addition, all orthologous k13 gene sequences from different nonhuman Plasmodium species and Plasmodium vivax with publicly available genomes were retrieved from the PlasmoDB (v38) (55) and MalAvi (56) databases (Data Sets S1 and S4).
Plasmodium falciparum samples/sequences.Plasmodium falciparum samples from Africa (n = 244) and Haiti (n = 103) and samples from U.S. malaria cases without travel histories (n = 116) reported to the CDC malaria surveillance network from 2014 to 2017 (n = 463) were used (Data Set S5). Genomic DNA isolated from P. falciparum strains 7G8 and HB3 was used as a positive control, and DNA isolated from whole blood from individuals without malaria was used as a negative control. The use of U.S. surveillance P. falciparum samples was approved by the Office of the Associate Director of Science, Center for Global Health, Centers for Disease Control and Prevention, as research not involving human subjects (protocol 2017-192). In addition, the methodology employed here for the amplification of the k13 gene was also approved by the same office as a surveillance activity (protocol 2017-309).
The data were expanded by retrieving complete and partial k13 gene sequences for P. falciparum from publicly available genomes/sequences found in the PlasmoDB (55) and NCBI (57) databases. All sequences were analyzed together with the new data obtained here from African and Haitian isolates (Data Set S5).
DNA extraction and Sanger sequencing for nonhuman primate parasites.Genomic DNA was extracted from whole blood of the nonhuman primate Plasmodium species using a DNeasy blood and tissue kit (Qiagen, GmbH, Hilden, Germany). Orthologous k13 genes from macaque and gibbon parasites were amplified by nested PCR using external primers forward 5′-C ATT TCC AAC TTC TCC GTC-3′ (58) and reverse 5′-TAT ATT TGC TAT TAR NAC NGA RTG-3′ and internal primers forward 5′-A TGA CGT ATG ADA GRG ART CNG-3′ and reverse 5′-AA TCT GGG AAC TAA TAR DGR DGG-3′. The primary PCR amplifications were carried out in a 50-μl reaction volume using 5 μl of total genomic DNA, 2.5 mM MgCl2, 1× PCR buffer, 1.25 mM each deoxynucleoside triphosphate, 0.4 mM each primer, and 0.03 U/μl AmpliTaq polymerase (Applied Biosystems, Roche, USA). The primary PCR conditions were a partial denaturation at 94°C for 4 min and 36 cycles with 1 min at 94°C, 1 min at 54°C, and a 2-min extension at 72°C, and a final extension of 10 min at 72°C was added in the last cycle. Then, the nested PCR mixtures were also made in 50-μl reaction volume using 1 μl of the primary PCR mixture, 2.5 mM MgCl2, 1× PCR buffer, 1.25 mM each deoxynucleoside triphosphate, 0.4 mM each primer, and 0.03 U/μl AmpliTaq polymerase. The nested PCR conditions were a partial denaturation at 94°C for 4 min and 25 cycles with 1 min at 94°C, 1 min at 56°C, and a 2-min extension at 72°C, and a final extension of 10 min at 72°C was added in the last cycle. Then, orthologous k13 genes from chimpanzee parasites, which share a recent common ancestor with P. falciparum, were amplified by using the same protocol described above, but with external primers forward 5′-GAA TTT TTC TAT NAC RTA YGA DAG-3′ and reverse 5′-ATT TGC TAT TAR NAC NGA RTG NCC-3′ and internal primers forward 5′-A TGA CGT ATG ADA GRG ART CNG-3′ and reverse 5′-AA TCT GGG AAC TAA TAR DGR DGG-3′. Then, two independent PCR products (50 μl) were excised from the gel (bands of approximately 2 kbp), purified using a QIAquick gel extraction kit (Qiagen, GmbH, Hilden, Germany), and cloned into the pGEM-T Easy vector system (Promega, USA). To detect mixed infections, at least two independent PCR products were cloned, and four clones were sequenced from each sample. Both strands for all the k13 orthologous genes were sequenced using an Applied Biosystems 3730 capillary sequencer. All the sequences obtained were identified as the k13 gene using analysis with the BLAST program (59). The sequences were submitted to GenBank and may be found under accession numbers MK103884 to MK103893.
DNA extraction from human samples and next-generation sequencing.DNA from human samples was extracted using the whole-blood protocol for the QIAamp DNA blood kit (Qiagen, GmbH, Hilden, Germany). Then, an end-to-end Illumina targeted amplicon deep sequencing (TADS) protocol and bioinformatics pipeline for molecular surveillance of drug resistance in P. falciparum, called malaria resistance surveillance (MaRS), was used to amplify the Pfk13 gene (60). First, the gene was amplified using a New England BioLabs (NEB) High Fidelity PCR kit (catalog no. 51104; New England BioLabs, USA) according to the manufacturer’s instructions with a 50.0-μl master mix preparation using 5× GC buffer. Next, 10.0 μl of each PCR product for each sample was used to generate the next-generation sequencing library. Pooled library fragment sizes were checked using an Agilent D5000 ScreenTape station (catalog numbers G2940C, 5067-4626, and G2940CA; Agilent Technologies), and the pooled library concentration was measured using a Qubit (v3.0) fluorometer (catalog numbers Q33216 and Q32853; Life Technologies Corp.). Lastly, a MiSeq desktop sequencer (Illumina Inc.) was used to obtain the sequences, and a consensus sequence for each sample was assembled using a pipeline. Additional details about the MaRS protocol are provided in Text S1 in the supplemental material of reference 60, including primer sequences, thermocycling conditions, and library amplification, purification, quantification, and normalization. The Pfk13 sequences obtained here were submitted to the GenBank database and may be found under accession numbers MK103421 to MK103883.
Phylogenetic-based method analyses of the k13 gene.First, a nucleotide sequence alignment was performed by using 42 sequences (1,985 bp, excluding gaps), including the entire orthologous k13 gene sequences obtained here, the sequences of the different haplotypes for some well-identified parasite species (e.g., P. reichenowi and P. gaboni), as well as sequences available in the PlasmoDB and MalAvi databases (Fig. 1 and 2; Data Set S1). Then, a second alignment based on 28 orthologous k13 gene sequences selected from different species of the first alignment was performed (duplicate haplotypes were not included; Fig. 3). All the alignments were generated using the ClustalX (v2.0.12) and Muscle programs implemented in the SeaView (v4.3.5) program (61) with manual editing. In both alignments, the Haemoproteus tartakovskyi (an avian haemosporidian parasite) k13 gene was included as an outgroup (56).
Using both alignments, phylogenetic relationships were inferred using the Bayesian method implemented in MrBayes (v3.2.6) software with the default priors (62). A general time-reversible model with gamma-distributed substitution rates and a proportion of invariant sites (GTR + Γ + I), which was the model with the lowest Bayesian information criterion (BIC) scores, as estimated by MEGA (v7.0.26) software, was used (63). Bayesian support for the nodes was inferred in MrBayes by sampling every 1,000 generations from two independent chains lasting 4 × 106 Markov chain Monte Carlo (MCMC) steps. The chains were assumed to have converged once the average standard deviation of the posterior probability was below 0.01 and the value of the potential scale reduction factor (PSRF) was between 1.00 and 1.02 (62). As a burn-in, 25% of the sample was then discarded once convergence was reached. The phylogenetic trees were visualized and edited using the FigTree (v1.4.3) program (64).
To detect changes in the rate of evolution, relative evolutionary rates of the k13 gene were estimated on the second alignment using a maximum likelihood framework with the RelTime method (30), as implemented in the command line version implemented in MEGA (vX) software (65). The substitution model was the same as that used for Bayesian analyses (GTR + Γ + I). Then, the tree was visualized in the FigTree (v1.4.3) program and colored using the relative evolutionary rate estimates of the k13 gene (Fig. 3). To quantitate the rate of evolution of this gene, absolute rates of synonymous and nonsynonymous substitutions per codon per My for the k13 gene were obtained using the CodonRates (v1.0) program in a Bayesian framework (33). In addition, by disentangling synonymous from nonsynonymous substitution rates, this analysis provides a first assessment of how natural selection may have affected the evolution of the gene. Because this method requires a tree topology and time calibration priors, including one for the in-group root node (the Plasmodium clade in this case), the Bayesian methods program BEAST (v2.4.1) (66) was used to estimate the origin of the Plasmodium species included in this investigation.
In BEAST, the relaxed clock (67) was applied with the GTR + Γ + I evolutionary model with a Yule model as the tree prior. Two independent chains lasting 108 MCMC steps were run until convergence. One scenario, using a combination of two calibration constraints that consider information on vertebrate hosts’ fossils (see the Paleobiology Database at http://www.paleodb.org/; last accessed 13 November 2017) was explored, and the results were compared with those of earlier studies based on the mitochondrial genomes (29). In all cases, uniform priors were used to calibrate divergences. Therefore, no particular time point within a given time interval was favored, and no punctual calibrations were used. The first calibration constraint considers that African parasites found in Mandrillus spp. (Plasmodium gonderi) diverged from those Plasmodium species found in Southeast Asia macaques (a monophyletic group) when Macaca branched from Papio (29, 68, 69). This event assumed that the divergence at that node took place at any time between 6 and 14.2 My ago (Ma) using a uniform prior. This time interval includes molecular estimates for the same Papio-Macaca divergence event (29, 69), and it is also consistent with the age of older fossils reported for Macaca spp. (see the Paleobiology Database at http://www.paleodb.org/; last accessed 13 November 2017). This calibration constraint considers the fact that P. gonderi has also been seen in Chlorocebus (a member of the Cercopithecini). The second fossil-based calibration constraint is the minimum of 23.5 Ma for the human-Macaca split (70), which is assumed to be the minimum time when the monophyletic group that includes P. malariae (a parasite found in humans) and all the Asian parasites originated. This event is used to inform a calibration prior with a maximum of 65 Ma that allows this parasite clade to be as old as the origin of primates (71). More details about this method can be found elsewhere (29, 69).
Given that the time estimated here for the avian Plasmodium-primate/rodent Plasmodium split (45.36 Ma; 95% credibility intervals, 33.23 to 61.10 Ma) was comparable to previous estimates also obtained using BEAST (45.31 Ma; credibility intervals, 38.95 to 51.96 [29]), the absolute rates of synonymous and nonsynonymous substitutions per codon per My were then estimated by using CodonRates (33) with a calibration prior of 45.4 ± 10 Ma for the in-group root (genus Plasmodium), as estimated from BEAST. This analysis was performed on a tree topology estimated with MrBayes that used Haemoproteus tartakovskyi as an outgroup.
Finally, to elucidate selective pressures during the evolution of the k13 gene, several methods were performed with the HyPhy package (34) implemented in the DataMonkey web server (http://www.datamonkey.org/; last accessed August 2018). Evidence of both episodic and pervasive positive selection at the level of an individual site was identified using a mixed-effects model of evolution, MEME (35), which detects sites evolving under positive selection under a proportion of branches. To investigate whether the k13 gene has experienced positive selection at least at one site on at least one branch, a branch site-unrestricted statistical test for episodic diversification (BUSTED) was used. This provides a gene-wide (not site-specific) test for positive selection (36). Sites positively selected were identified at a significance level (P value) of <0.05. In addition, to test whether the strength of natural selection has been relaxed or intensified along a specified set of test branches, the RELAX method (37) was performed. RELAX is most useful for identifying trends and/or shifts in the stringency of natural selection on a given gene.
Plasmodium falciparum population analyses.To study the genetic relationships among worldwide haplotypes of the k13 gene, a third alignment was performed using 982 cosmopolitan complete sequences, including all the k13 gene sequences obtained here for P. falciparum population samples (n = 463), as well as those available in the PlasmoDB and NCBI databases (n = 519). Evidence of natural selection was explored by estimating the average number of synonymous substitutions per synonymous site (dS) and the average number of nonsynonymous substitutions per nonsynonymous site (dN) between a pair of sequences under the Nei-Gojobori method with the Jukes and Cantor corrections, as implemented in the MEGA (v7) program (63). The difference between dS and dN with its standard error (estimated using bootstrap analysis with 1,000 pseudoreplications) was calculated, and significance was tested by a two-tailed codon-based Z-test on the difference between dS and dN, as described by Nei and Kumar (72). The null hypothesis was that the observed polymorphism is not under selection (dS = dN).
The alignments were used to estimate median-joining networks using Network (v4.6.1.0) software (Fluxus Technologies, 2011). Transversions were set equal to transitions, and the epsilon parameter was set equal to 0 with only one round of star contraction, which collapses star-like structures in the network into single nodes. The total number of sites included in this analysis, excluding gaps or missing data, was 2,137 bp out of 2,181 bp using the P. falciparum 3D7 strain as a reference. Estimates of nucleotide diversity, haplotype diversity, and Tajima’s D test parameter were obtained using DnaSP (v5.50) software (73). Lastly, the frequency of candidate/validated mutations located in the propeller domain of a P. falciparum kelch (Pfk13) gene that have been associated with clinical delayed parasite clearance was estimated by country, as well as by continent, using a fourth alignment constructed with 3,916 partial (n = 2,934) and complete (n = 982) k13 gene sequences (Table 2; Data Sets S5 and S6).
ACKNOWLEDGMENTS
We thank Michael Li for the silhouettes design and the DNA laboratory at the School of Life Sciences (Arizona State University) for technical support.
This study was supported in part by grants from the U.S. National Institutes of Health (2U19 AI089681-08 and R01AI130473). Computer resources were supported by a grant from the Pennsylvania Department of Health (TU-420721).
We declare that no conflicts of interest exist.
FOOTNOTES
- Received 5 December 2018.
- Returned for modification 24 February 2019.
- Accepted 28 April 2019.
- Accepted manuscript posted online 13 May 2019.
Supplemental material for this article may be found at https://doi.org/10.1128/AAC.02550-18.
- Copyright © 2019 American Society for Microbiology.
