- Open Access
Malaria in Venezuela: changes in the complexity of infection reflects the increment in transmission intensity
Malaria Journal volume 19, Article number: 176 (2020)
Malaria incidence has reached staggering numbers in Venezuela. Commonly, Bolívar State accounted for approximately 70% of the country cases every year. Most cases cluster in the Sifontes municipality, a region characterized by an extractive economy, including gold mining. An increase in migration to Sifontes, driven by gold mining, fueled a malaria spillover to the rest of the country and the region. Here samples collected in 2018 were compared with a previous study of 2003/2004 to describe changes in the parasites population structures and the frequency of point mutations linked to anti-malarial drugs.
A total of 88 Plasmodium falciparum and 94 Plasmodium vivax isolates were collected in 2018 and compared with samples from 2003/2004 (106 P. falciparum and 104 P. vivax). For P. falciparum, mutations linked to drug resistance (Pfdhfr, Pfdhps, and Pfcrt) and the Pfk13 gene associated with artemisinin delayed parasite clearance, were analysed. To estimate the multiplicity of infection (MOI), and perform P. falciparum and P. vivax population genetic analyses, the parasites were genotyped by using eight standardized microsatellite loci.
The P. falciparum parasites are still harbouring drug-resistant mutations in Pfdhfr, Pfdhps, and Pfcrt. However, there was a decrease in the frequency of highly resistant Pfdhps alleles. Mutations associated with artemisinin delayed parasite clearance in the Pfk13 gene were not found. Consistent with the increase in transmission, polyclonal infections raised from 1.9% in 2003/2004 to 39% in 2018 in P. falciparum and from 16.3 to 68% in P. vivax. There is also a decrease in linkage disequilibrium. Bayesian clustering yields two populations linked to the time of sampling, showing that the parasite populations temporarily changed. However, the samples from 2003/2004 and 2018 have several alleles per locus in common without sharing multi-locus genotypes.
The frequency of mutations linked with drug resistance in P. falciparum shows only changes in Pfdhps. Observations presented here are consistent with an increase in transmission from the previously circulating parasites. Following populations longitudinally, using molecular surveillance, provides valuable information in cases such as Venezuela with a fluid malaria situation that is affecting the regional goals toward elimination.
Although the global malaria burden has been declining for almost two decades , Latin America is an exception to this trend. Malaria incidence dropped in the region from 2000 to 2015, resulting in a 61.2% reduction in cases. However, it started to increase from 2015  with Brazil, Nicaragua, and Venezuela driving this surge. During 2018, due to the worsening of the economic and political crisis, Venezuelan cases dramatically increased, accounting for 51% of the region reported cases, displacing Brazil, which reported 23% .
Like in other South-American countries [2,3,4], Plasmodium vivax accounts for nearly 77% of Venezuela reported cases followed by Plasmodium falciparum (17%), mixed P. vivax/P. falciparum infections (5%), and Plasmodium malariae with less than 1% . Historically, malaria incidence remained relatively focalized in the endemic areas located in the Bolívar State (bordering Brazil) until 2014 (70–80% of reported cases in the country). An economic and political crisis fueled a flow of migrants to gold mines in Bolívar State, boosting incidence locally, followed by a spillover to other areas within Venezuela, and across international borders . Furthermore, the connectivity among the countries that form the Guiana Shield (Guyana, French Guiana, Suriname, together with parts of Venezuela, and Brazil)  has the potential of increasing the regional impact of the Venezuelan situation.
Regional connectivity may not only affect malaria in the region in terms of imported cases but also because the local epidemiology in the area has allowed for the independent emergence of P. falciparum mutations linked to anti-malarial drug resistance [6,7,8,9,10,11,12,13]. Thus, such emergent drug-resistant genotypes are at risk of spilling over the region. Hence, to provide must needed baseline information, the genetic makeup of the P. vivax and P. falciparum populations circulating in Bolívar State, Venezuela was described.
The premise is that a dramatic increase in incidence, as the one observed in Venezuela, should lead to changes in the parasite population in terms of simple metrics, particularly multiplicity of infection (MOI) and linkage disequilibrium (LD) . Also, a population expansion may affect the frequency and dispersion of mutations linked to anti-malarial drugs [14, 15]. Thus, this study compares the population genetic structure and MOI in a population of the Bolívar State between samples collected in 2003/2004, before the ongoing crisis, with samples collected in 2018. All these samples were part of the malaria surveillance programme. In particular, samples from Tumeremo, the capital of the municipality of Sifontes, and surrounding villages in that municipality, were analysed. These human settlements act as a catchment area to those migrating to Brazil and/or working in gold mines.
The hypothesis was that the average multiplicity of infection should increase, and linkage disequilibrium decrease as a result of the sustained demographic expansions of both major malarial parasites. Frequency of mutations were also compared to previously used anti-malarial drugs in P. falciparum and interrogated those samples for mutations liked to the artemisinin delayed parasite clearance in the Pfk13 propeller.
Study sites, sample collection, and DNA isolation
Bolívar is the largest Venezuelan state (242,801 km2), it is in the southeastern part of the country bordering Brazil and Guyana. Bolívar accounts for approximately 70% of the Venezuelan cases, mainly among miners and indigenous groups. The state malaria cases are clustered in the Sifontes municipality, an area where gold mining is the major economic activity [5, 12]. Between June 2018 and August 2018, as part of a surveillance study, a total of 182 symptomatic volunteers infected with Plasmodium spp. were passively recruited when visiting the health posts for malaria diagnosis. All patients with confirmed malaria infection by microscopy (age range = 13–73 years) were invited to provide a blood sample for further characterization of the parasite, and they received verbal and written explanations about the study. Then, a signed informed consent (IC) form was requested for those patients that agreed to participate. For further molecular analyses to confirm diagnostic and other clinical aspects, a total of 8 mL of blood was collected by venipuncture from every patient, of which 5 mL were stored in tubes containing EDTA and in Whatman™ 903 Protein Saver Cards (GE Healthcare Life Sciences, UK). Malaria infection was determined by microscopic examination of Giemsa-stained thick blood smears (TBS). The study protocol was approved by the Ethics Committee “Complejo Hospitalario Universitario Ruiz y Páez.”
To estimate the frequency of drug-resistant mutations in P. falciparum and other analyses in both parasites (P. falciparum and P. vivax), genomic DNA was extracted from anonymized whole blood samples using QIAamp DNA Micro Kit (Qiagen, GmbH, Hilden, Germany).
Drug resistance genes and Pfk13 gene
A subset of 50 samples were genotyped for P. falciparum mutations at: (1) dihydropteroate synthase (Pfdhps) codons 436, 437, 540, 581, and 613; (2) bifunctional dihydrofolate reductase-thymidylate synthase (Pfdhfr) codons 50, 51, 59, 108, and 164; and (3) chloroquine resistance transporter (Pfcrt) codons 72, 73, 74, 75, 76, and 350. In addition, almost the complete kelch protein Pfk13 was also amplified (2164 bp out of 2181 bp) for the eighty-eight P. falciparum samples gathered in this study. The PCR primer sequences used for each gene were: (1) forward 5′-AAC CTA AAC GTG CTG TTC AA-3′ and reverse 5′-AAT TGT GTG ATT TGT CCA CAA-3′ for Pfdhps; (2) forward 5′-ATG ATG GAA CAA GTC TGC G-3′ and reverse 5′-GAA TTC TTC TAC TTG TAG GAT C-3′ for Pfdhfr; and (3) forward 5′- TTA CAT ATA ACA AAA TGA AAT TCG C-3′ and reverse 5′-TAT TGT GTA ATA ATT GAA TCG ACG-3′ for Pfcrt. In the case of Pfk13 gene, 88 samples were successfully amplified by nested PCR using 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′-ATG ACG TAT GAD AGR GAR TCN G-3′ and reverse 5′-AAT CTG GGA ACT AAT ARD GRD GG-3′ .
For Pfdhps, Pfdhfr, Pfcrt, and Pfk13 genes, PCR amplifications were carried out in 50 µl volume reaction using 5 µl of total genomic DNA, 2.5 mM MgCl2, 1X PCR buffer, 1.25 mM of each deoxynucleoside triphosphate, 0.4 mM of each primer, and 0.03 U/µl AmpliTaq polymerase (Applied Biosystems, Roche-USA). The PCR conditions for Pfdhps were a partial denaturation at 95 °C for 7 min, and 40 cycles with 30 s at 95 °C, 30 s at 50 °C and 1 min extension at 68 °C, and a final extension of 5 min at 68 °C in the last cycle. For Pfdhfr, the conditions were a partial denaturation at 95 °C for 7 min, and 40 cycles with 1 min at 95 °C, 1 min at 54 °C and 2 min extension at 72 °C, and a final extension of 10 min at 72 °C in the last cycle. For Pfk13 gene, the conditions for the primary PCR were a partial denaturation at 94 °C for 4 min, and 36 cycles with 1 min at 94 °C, 1 min at 53 °C and 2 min extension at 72 °C, and a final extension of 10 min at 72 °C in the last cycle. The same PCR conditions were used for the nested PCR but with 56 °C for the annealing temperature. In the case of Pfcrt, PCR reactions were carried out in 50 μl using TaKaRa LA Taq Polymerase (TaKaRa Mirus Bio Inc.) following manufacturers’ directions. Amplification conditions for the PCR were a partial denaturation at 94 °C for 1 min and 30 cycles with 30 s at 94 °C and 5 min at 60 °C, followed by a final extension of 10 min at 72 °C. Standard P. falciparum and P. vivax DNA positive and negative (nuclease-free dH2O) controls were included in each batch of PCR. Then, PCR products (50 μl) were excised from the gel and purified using QIAquick® Gel extraction kit (Qiagen, GmbH, Hilden, Germany). Both strands for PCR products were directly sequenced using an Applied Biosystems 3730 capillary sequencer. After careful inspections of each electropherogram, mutations associated with drug resistance (Pfdhps, Pfdhfr, and Pfcrt genes) or artemisinin delayed parasite clearance resistance were recorded. Sequences obtained in this study for kelch protein Pfk13 were deposited in GenBank under the accession numbers MT120216 to MT120303.
An alignment with 1099 complete Pfk13 sequences was performed to explore the genetic relationships among the Venezuelan alleles with those reported worldwide. The alignment included 65 Pfk13 complete gene sequences obtained here for Venezuelan samples, as well as those available in PlasmoDB and NCBI databases (N = 1034). Then, a median-joining network was estimated by using Network v22.214.171.124 (Fluxus Technologies 2011). Transversions were set equal to transitions, and the epsilon parameter 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, were 2155 bp out of 2,181 bp using the P. falciparum 3D7 strain (PF3D7_1343700) as a reference. Estimates of haplotype diversity and Tajima’s D test were obtained using DnaSP software, version 6 . Also, evidence of natural selection was explored. In particular, the average number of synonymous substitutions per synonymous site (dS) and the average number of nonsynonymous substitutions per nonsynonymous site (dN) were calculated between all pairs of sequences from the Venezuelan samples (2018). This analysis was carried out using the Nei-Gojobori method with the Jukes and Cantor corrections, as implemented in the MEGA v7 software . The null hypothesis tested was that the observed polymorphism is not under selection (dS = dN). The difference between dS and dN with its standard error (estimated using bootstrap analysis with 1000 pseudo-replicates) was calculated, and a two-tailed codon-based Z-test tested for significance on the difference between dS and dN .
Microsatellite (STRs) genotyping
Genotyping was performed using fluorescently labeled PCR primers for a set of eight standardized microsatellites for P. falciparum (POLYa, TAA60, ARA2, Pfg377, TAA109, TAA81, TAA42, and PfPK2  and nine for P. vivax (MS1, MS2, MS5, MS6, MS8, MS10, MS15 , 14.185, and 2.21, ). PCRs were performed in 15 μL reactions with 2 μL of total genomic DNA, 0.25 mM of each primer, and 7.5 μL of PCR Master Mix (Promega, USA) (it contains 0.05 U/μL of Taq DNA polymerase, 2X reaction buffer, 0.4 mM each dNTP, and 3 mM MgCl2; Promega, USA). P. falciparum and P. vivax DNA positive samples and negative (nuclease-free dH2O) controls were included in each batch of PCR, respectively. Amplification conditions for PCRs were reported elsewhere, depending on the set of primers [20,21,22]. Fluorescently labelled PCR products were separated on an Applied Biosystems 3730 capillary sequencer and scored using GeneMarker v2.6.7 (SoftGenetics LLC). All alleles were scored at a given locus if minor peaks were more than one-third the height of the predominant peak. The finding of one or more additional alleles at any locus was interpreted as a multiple (polyclonal) infection with two or more haploid genotypes in the same isolate (transmitted by one or several mosquitoes). A sample with a single infection was considered if it had only one allele per locus at all the genotyped loci. Missing data (no amplification) were reported by locus but not considered for defining multilocus genotypes.
The multiplicity of infection (MOI) and allele frequencies
MOI was defined as the average number of distinct parasite (P. falciparum or P. vivax) genotypes concurrently infecting a patient. Thus, average MOI in these P. falciparum and P. vivax populations was estimated (per loci—ignoring missing data-and on average per loci/year) as the average number of super-infections (neglecting co-infections) using a maximum-likelihood (ML) method that allows estimating profile-likelihood confidence intervals [23, 24]. Then, MOI and allele frequencies were calculated using the R-script provided in . In the case of allele frequencies, they were estimated for each locus using ML method .
Plasmodium falciparum and P. vivax population genetic analysis
A suite of approaches was used to characterize the genetic variation in the circulating P. falciparum and P. vivax and the differentiation of the parasite populations between 2003/2004 and 2018. In particular, the microsatellite data gathered in this investigation were compared against previous data reported for the same location in 2003/2004 . The genetic diversity within each parasite population was estimated using a series of summary statistics implemented in the Haplotype Analysis software v1.05 . For both parasites (P. falciparum and P. vivax), the number of different sampled multilocus genotypes (SMG), the number of unique genotypes (G), the number of private genotypes (PG), and the Nei’s index of genetic diversity (He) were estimated  on all the multi-locus genotypes that were unambiguously identified. He was estimated as
where n is obtained by taking the sum of identifiable genotypes (phased) overall samples, and pi is the relative frequency of the i-th haplotype (i = 1, …, L) in all sampled haplotypes. He gives the average probability that a pair of alleles randomly selected from the population is different. Complex infections with differences at more than two loci were not included in this analysis because the haploid genotypes could not be inferred. He was also calculated per locus. In this case, pi is the frequency of allele i, and He is the average probability that a pair of alleles randomly selected from the population is different. Non-parametric bias-corrected and accelerated bootstrap confidence intervals were estimated based on 1000 bootstrap replications using the jackknife estimate for the acceleration factor, as described in .
Two methods were used to identify the parasite clusters (P. falciparum and P. vivax) that circulated in Bolivar State, (i) principal component analyses (PCA) that do not use explicit admixture model (but that is not assumption-free), and (ii) a Bayesian model-based clustering algorithm that considers admixture as implemented in the Structure v2.3.4 software . PCAs were estimated for both population parasites (P. falciparum and P. vivax) in R-script. To include samples with missing data and samples with multiple infections, for each PCA, alleles at each locus were coded as 0–1 variable, indicating the absence and presence of the allele. Hence, the number 0–1 variables for each locus reflects the number of alleles at that locus. This allowed exploring the effect of potential amplification errors without eliminating samples.
For the Bayesian clustering approach, which assigns genotypes to K populations or clusters characterized by a set of allele frequencies at each locus, the observed genetic diversity was evaluated under different K values (K = 1 to 12), and each K value was run independently 15 times with a burn-in period of 100,000 iterations followed by 100,000 iterations. The admixture model that allows for the presence of individuals with ancestry in two or more of the K populations was used in all the analyses . The posterior probability for each number of populations or clusters (K) was computed, and the K-value that better explains the genetic data was an estimate of the number of parasites circulating clusters. Delta K values were computed using Structure Harvester  and then, CLUMPP (Cluster Matching and Permutation Program) was used to facilitate the interpretation of population-genetic clustering results . Finally, distruct v1.1 was used to display the results graphically . For this analysis, complex infections with differences at more than two loci were not included.
To evaluate the differentiation in time between the parasite populations, 2003/2004 vs. 2018, normalized fixation indexes (Fst) were estimated; their significance was assessed using a randomization test. A limitation in this analysis is that haploid genotypes cannot be inferred on samples with multiple alleles at two or more loci, so such samples were excluded. Considering that these are neutral physically unlinked loci, measurements of linkage disequilibrium (LD) between pair of loci were also estimated to detect potential clonal expansions of both parasite populations. Conditional asymmetric linkage disequilibrium (ALD) measurements were used since they consider differences in the number of alleles in each pair of loci (many pairs of loci have a different number of alleles) . ALD estimates go from 0 to 1 with 0 implying total independence and 1 complete linkage. ALD requires frequency estimates of the two-locus haplotypes observed at each pair of loci under consideration. Thus, two-locus haplotypes were inferred between all pairs of loci, and their frequencies estimated using the ML method of  using the R-script provided in . If there were multiple alleles at two loci under consideration for a given pair of samples, their haplotypes could not be phased, and just that comparison was disregarded for that pair of samples.
Finally, population genetic analyses were complemented by inferring the haplotype genealogies for P. falciparum and P. vivax populations contrasting the data obtained from 2003/2004 against those from 2018. These haplotype genealogies were estimated on eight microsatellites for P. falciparum and nine for P. vivax, using the Global Optimal eBURST algorithm , as implemented in PHYLOViZ . These analyses only included single infections and complex infections with differences at only one locus, given that the haploid genotypes can be inferred. A Minimum Spanning Tree-like (MST) structure was drawn to cluster the 198 sequence types (STs) for P. falciparum and 159 for P. vivax into a clonal complex (CC) based on their multilocus genotypes by using an extension of the goeBURST rules up to n-locus-variants-level (nLV, where n equals to the number of loci in these datasets: eight for P. falciparum and nine for P. vivax respectively). In addition, to better explore the resulting tree and the relationships between genotypes, an NLV graph was obtained using a different operation that modifies the default characteristics of the MST. This kind of NLV graph easily identifies sets of closely related nodes by relaxing the MST construction restriction, allowing the display of all possible links up to a specific threshold.
In this study, 88 patients were infected with P. falciparum and 94 with P. vivax. Out of the 88 infected with P. falciparum, 19 were co-infected with P. vivax (mixed infections) detected by PCR.
Drug resistance genes
No multiple-strain infections were detected in the electropherograms. In 2018, mutations associated with sulfadoxine-pyrimethamine (SP) resistance in Pfdhfr and Pfdhps genes revealed four multilocus genotypes: (1) Pfdhfr (50R/51I/108 N) linked to Pfdhps (437G/540E/581G) present in 58.1% of the samples (hereinafter referred to as a sextuple mutant), (2) Pfdhfr (50R/51I/108 N) linked to Pfdhps (437G/581G) in 32.6% of samples, (3) Pfdhfr (50R/51I/108 N) linked to Pfdhps (437G) in 7% of the samples, and (4) Pfdhfr (51I) linked to Pfdhps (437G) present only in 2.3% of the samples (only one patient) (Table 1). These results contrast with a previous report where the sextuple mutant was present in 90.7% of the samples, and a quadruple mutant Pfdhfr (51I/108 N) linked to Pfdhps (437G/581G) was present in 9.3% of the samples . Although, in 2003/2004 two synonymous of Pfcrt chloroquine resistance alleles, StctVMNT (91%) and SagtVMNT (9%) were circulating in the population , in 2018 only the allele StctVMNT was found in the set of samples used (Table 1).
Pfk13 population analyses
Eighty-eight Pfk13 sequences were obtained in this study for the group of samples collected in 2018. The propel region was sequenced for all 88 samples to check for mutations associated with artemisinin delayed parasite clearance phenotype. However, only 65 Pfk13 alleles were sequenced completely and used for the genetic diversity analyses. For this dataset, a very low genetic diversity was found in Pfk13 gene (π = 0.0005 ± 0.0), with only two haplotypes circulating in the area (Table 2). In addition, the polymorphism measured using the number of polymorphic (segregating) sites had an S value of 1 (Table 2). Indeed, only two sites showed substitutions in the Venezuelan sequence alignment, one was a nonsynonymous substitution outside of the propeller domain (K189T; 55.4%), and the other was a synonymous substitution (VGTT636VGTG, 1.5%) in the propeller domain (440-680 bp). More important, none mutations associated with the delayed parasite clearance phenotype were found in the Venezuelan samples.
Then, these sequences (N = 65) were analysed together with 1034 Pfk13 gene sequences from around the globe, including samples from The Americas (Table 2, Additional file 1). Thus, the K189T mutant found in Venezuelan samples was the most predominant in Africa, Brazil, and French Guyana. In contrast, the VGTT636VGTG found in one Venezuelan sample was previously reported in two sequences from India and one from Nigeria (Additional file 1). When sequences from Central and South America were compared, there were only two nonsynonymous (G496V and S679P) and four synonymous substitutions (G453G, K503K, G638G, and F673F, Additional file 1) within the propeller domain in the Haitian sequences (N = 103) and all of them in very low frequency (0.9–8.7%). It worth noticing that none of them have been associated yet with the clinical delayed parasite clearance. Interestingly, none substitutions were found in the propeller domain from other countries of The Americas.
Worldwide, a total of 124 Pfk13 haplotypes were found with a haplotype (allele) diversity (Hd) of 0.7031 (standard deviation [SD] = 0.0117) (Table 2). These results are similar to those recently reported by Pacheco et al. . 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 alleles H1 and H2 (Fig. 1) with a broad distribution (Asia, Africa, and South America). One of the alleles (H1, n = 525, 47.8%) is mostly distributed in Asia, and the other (H2, n = 282, 25.7%) in Africa. Interestingly, both haplotypes contain sequences from Venezuela and Brazil. Specifically, for Venezuelan samples, out of the 65 complete sequences, 36 (55.4%) corresponded to the H1 and 29 (44.6%) to the H2, which contains the K189T mutant (Fig. 1). A third haplotype (H3, n = 50, 4.5%), in terms of frequency, is only found in Thailand and contains the C580Y mutation found in the SEA region (Fig. 1). Nevertheless, most of the haplotypes were restricted to single countries (Fig. 1).
The multiplicity of infection and allele frequencies
Overall, for P. falciparum and P. vivax infections, the prevalence of multiple infections showed a substantial increase from 2003/2004 to 2018. Many of those P. falciparum and P. vivax multiple infections were the result of having more than one allele in one or two loci, especially in P. vivax (Table 3), so the lineage-specific genotypes were only inferred for those samples with more than one allele at only one locus. In 2018, among the P. falciparum cases, 39% of the 77 samples successfully genotyped (using microsatellites) had infections with more than one lineage in at least one locus. In contrast, 68% out of 94 P. vivax samples were found with multiclonal infections (Table 3). This difference translated into a slightly higher MOI in P. vivax (> 1.5 vs. < 1.5, Fig. 2) with overall more complex infections (several loci with up to three alleles). In addition, the estimated MOI parameter showed that the mean MOI also substantially increases for each microsatellite marker from 2003/2004 to 2018 in both parasite populations, especially for P. vivax (Fig. 2b). This increase was statistically significant for both parasite species (Additional file 2).
Population genetic analysis
Using those samples for which their haplotypes could be inferred, the number of sampled multilocus genotypes (SMG) from the human specimens, the number of distinct genotypes (G), the number of private genotypes (PG), and the Nei’s index of genetic diversity (He) were estimated for each population using Haplotype Analysis software v1.05 and are shown in Table 4. Interestingly, mean genetic diversity was high and similar in both parasites (P. falciparum-mean He: 0.939 and P. vivax-mean He = 0.938), and none of the P. falciparum and P. vivax genotypes found in 2003/2004 were recovered in 2018. Fst values were significant and relatively high for both parasites (Table 4). When He was also calculated per locus, there was some overall tendency to increase heterozygosity but not for each locus in the case of P. falciparum (only 4 loci, Fig. 3a and Additional file 3). However, for P. vivax, there were high levels of heterozygosity in 2003/2004 and 2018. There was an overall tendency to increase heterozygosity in 2018 but not for each locus (only 5 loci, see Fig. 3b). The number of alleles per locus was similar between the sampled years for both parasite populations (Additional files 4 and 5), but with a tendency to increase the number of alleles in 2018 for P. vivax with a greater number of low-frequency-variants.
In the case of the PCA analyses, for P. falciparum dataset, 64.2% of the variance between the sampled years was explained by the first 4 components (Fig. 4a, Additional file 6) and the first PC component separates both indicating a shift between the two time points (2003/2004 and 2018). In contrast, the four PCs explained around 34.5% of the variance between the two P. vivax dataset (Fig. 4b, Additional file 6). By analysing the posterior probabilities for each clustering from K = 1 to 12, the clustering obtained with K = 2 showed the highest Delta K values (Delta K = mean (|L”(K)|)/sd(L(K))) for both Plasmodium species. Thus, the Bayesian clustering using the Structure v2.3.4 software  yield two populations (K = 2) for P. falciparum and P. vivax and linked to the time of sampling showing that both parasites populations diverged in time (Fig. 4).
Although the asymmetric linkage disequilibrium (ALD) estimated for P. falciparum was stronger for 2003/2004 than in 2018, the values per locus were low for both sampled periods (Fig. 5a). However, for P. vivax, the pattern was slightly different, with strong ALD in 2003/2004 and almost absent in 2018 (Fig. 5b). Low ALD of physically unlinked loci, as used here, is contrary to the expectations under a bottleneck or clonal expansion scenarios. Lastly, the haplotype genealogies inferred for P. falciparum (Fig. 6a.1) and P. vivax (Fig. 6b.1) populations and the NLV graph (Fig. 6a2, b2) showed no clear boundaries between the sampled periods.
The increase in malaria cases in Venezuela, and the subsequent mass migration across international borders, is likely the most important event affecting the goal toward malaria elimination in Latin-America [5, 13, 35]. It is worth noting that the long-term impacts of the current situation in Venezuela, in the context of cross-border malaria, are still poorly understood. Most reports focus on counting the number of imported cases [5, 35]. Yet, data are being collected in other aspects, such as the prevalence of mutations linked to drug resistance in P. falciparum [12, 13]. Here, baseline data on the circulating parasites found in Bolívar State in 2018, the epicentre of the current malaria epidemic affecting Venezuela, were described. Although there are other factors, the surge of unplanned gold mining activities with limited access to effective anti-malarial treatment is considered the leading cause that drove the increase in malaria transmission observed in the Bolivar state [4,5,6].
Bearing in mind the limitations of the sample size, it is worth noting that 21.5% (19 of the 88 P. falciparum positive parasites) were mixed infections with P. vivax detected when using PCR. These results differ from the 6.8% reported from the same area in a recent study , but using different protocols. Thus, it may be possible that the reported national average of mixed infections (5%) is an underestimation. Considering that all those infections were diagnosed as P. vivax, inadequate anti-malarial treatment could lead to an increase in “new infections” by P. falciparum that were undetected mixed infections treated with chloroquine.
Although SP and chloroquine are no longer used to treat P. falciparum, fluctuations in the frequency of mutations linked to resistance were observed in 2018 when compared to 2003/2004. Unlike the situation found in some African countries [37, 38], the lack of Pfcrt-wild type in the P. falciparum population did not allow the recovery of sensitive chloroquine parasites even when the drug is no longer used to treat falciparum malaria. This scenario was suggested by early studies [8, 9]. The only way that chloroquine-sensitive parasites could re-emerge in this region is via their reintroduction by migration, something that has not been documented, or the emergence of a new mutation that revert to chloroquine sensitivity . It could be speculated that the P. falciparum population is still under some chloroquine drug pressure. The drug is commonly used to treat uncomplicated P. vivax, and there could be a treatment spillover due to mixed infections or inadequate diagnostic . The fixation of the Pfcrt StctVMNT can be explained by genetic drift since it was the most common allele in the region [8, 9].
The situation is different in SP. Like has been reported in Peru [40, 41], there is a decline in the highest resistant alleles. Here, however, that decline was observed on the frequency of the Pfdhps triple mutant that confers high resistance to sulfadoxine, one of the SP components. It is also worth noting the presence of a single Pfdhps mutant; this allele was not observed in an early study . Thus, the regional changes in malaria treatment seem to have affected at least the frequency of mutations linked to sulfadoxine resistance.
Mutations linked to artemisinin delayed phenotype in Pfk13 were not found in the 88 isolates genotyped for this gene; a result consistent with a recent study on samples collected on Venezuelan migrants in Brazil . This is an important finding since there is constant movement between Guyana, Brazil, and Venezuela as a result of the illegal exploitation of gold mines . Fortunately, it seems that the Pfk13 mutations reported in the region  remind contained. Still, the circulating Pfk13 alleles show polymorphism with mutations observed in low frequency in other areas. E.g., the K189T mutant found in Tumeremo has been found in Africa and neighbouring countries such as Brazil and French Guyana. Indeed, the alleles circulating in Venezuela belong to the two major alleles already reported . Such local substructure likely accounts for the positive, but no significant, Tajima’D value reported here (Table 2). The world sample with N = 1099 sequences has a significant and negative Tajima’D consistent with a population expansion or selective sweep, the later less probable simply because this sample of worldwide alleles only has a few of the mutations conferring resistance. The bias on reporting resistant alleles has to do with how the data on Pfk13 mutations is usually generated .
Nevertheless, the worldwide presence of mutations in low frequency in the Pfk13 gene could be explained, at least in part, by the P. falciparum population expansion. This demographic process could make natural selection less efficient to eliminate otherwise semi-deleterious mutations . Such a repertoire of mutations in low frequency may provide the variation for drug-mediated selective pressure to act.
The observed increase in the number of polyclonal infections (multiplicity of infection or MOI) in both parasites between 2003/2004 and 2018 correlates with the overall changes in transmission [14, 42,43,44,45,46]. Although MOI measures the frequency of superinfections and co-infections , its relationship with transmission intensity is not linear [42,43,44]. In particular, comparing MOI between regions that differ in transmission does not seem to be a sound approach if those populations also differ in other epidemiological and historical characteristics, including the genetic diversity/evolutionary history of the circulating parasites [42,43,44, 46]. Understandably, MOI is also affected by those other factors making it a poor predictor of transmission intensity in such contexts [43,44,45]. Here, however, the increase in the average MOI appears to be explained by the surge in transmission.
The differences in MOI are interesting also because the genetic diversity in Tumeremo has a moderate increase in both parasites locally. There was observable differentiation when the two different time points are compared, in both P. falciparum and P. vivax, using PCA and structure analyses indicated temporal structure (Fig. 4). However, such differentiation seems to be explained by recombination and changes in allele frequencies from a common allele pool rather than migration, as indicated by the haplotype network (Fig. 6). One of the factors that indicate this expansion from the original population are the changes in linkage disequilibrium, as it decreases in both parasites when 2018 parasites are compared with those collected in 2003/2004. Thus, as transmission intensifies and the frequency of polyclonal infections increased, also the rate of outcrossing increased, breaking the linkage that originally was sustained by inbreeding when malaria was under control in 2003/2004 among these microsatellites that are physically unlinked. MOI, in addition to changes in transmission, has been linked to disease severity for P. vivax in some settings . However, its relationship with clinical disease severity may depend on several factors, including how MOI is measured in terms of loci, sampling, clinical endpoints, and age, among other confounding factors [42,43,44,45,46]. A clinical study should provide some insights into the role of MOI in setting such as Venezuela, where transmission has reached unprecedented high levels. Nevertheless, the observed changes in MOI and linkage disequilibrium are consistent with the use of these metrics to follow transmission intensity longitudinally [45,46,47,48].
The frequency of mutations linked with drug resistance changed in SP, where the highly resistant triple mutants decreased in frequency, particularly in the case of Pfdhps gene. These have been observed in other areas in Latin America where SP is no longer used to treat P. falciparum malaria. Importantly, no evidence of Pfk13 mutations associated with artemisinin delayed clearance were found here. Overall, population-level parameters such as the average MOI and linkage disequilibrium are internally consistent and affected by the increase in malaria incidence observed in Venezuela. These findings indicate that these parameters are informative when following populations longitudinally. Patients in Venezuela harbour more complex infections nowadays, and that criteria may allow separating a potentially imported case from a Venezuelan migrant that acquired malaria locally since MOI tend to be lower in other regions of Latin America.
Availability of data and materials
Sequences were deposited in the GenBank with the accession numbers MT120216 to MT120303.
Asymmetric linkage disequilibrium
- Pfdhfr :
Bifunctional dihydrofolate reductase-thymidylate synthase
- Pfdhps :
Nonsynonymous substitutions per nonsynonymous site
Synonymous substitutions per synonymous site
Number of unique genotypes
Nei’s index of genetic diversity
Multiplicity of infection
Minimum Spanning Tree-like
Principal component analyses
Polymerase chain reaction
- Pfcrt :
Chloroquine resistance transporter
Number of private genotypes
Number of segregating sites
Sampled multilocus genotypes
Standard microsatellite loci
Thick blood smears
WHO. World Malaria Report 2019. Geneva: World Health Organization; 2019.
Chowell G, Munayco CV, Escalante AA, McKenzie FE. The spatial and temporal patterns of falciparum and vivax malaria in Perú: 1994-2006. Malar J. 2009;8:142.
Arevalo-Herrera M, Quiñones ML, Guerra C, Céspedes N, Giron S, Ahumada M, et al. Malaria in selected non-Amazonian countries of Latin America. Acta Trop. 2012;121:303–14.
Recht J, Siqueira AM, Monteiro WM, Herrera SM, Herrera S, Lacerda MVG. Malaria in Brazil, Colombia, Peru and Venezuela: current challenges in malaria control and elimination. Malar J. 2017;16:273.
Grillet ME, Hernández-Villena JV, Llewellyn MS, Paniz-Mondolfi AE, Tami A, Vincenti-Gonzalez MF, et al. Venezuela’s humanitarian crisis, resurgence of vector-borne diseases, and implications for spillover in the region. Lancet Infect Dis. 2019;19:e149–61.
Musset L, Pelleau S, Girod R, Ardillon V, Carvalho L, Dusfour I, et al. Malaria on the Guiana Shield: a review of the situation in French Guiana. Mem Inst Oswaldo Cruz. 2014;109:525–33.
Cortese JF, Caraballo A, Contreras CE, Plowe CV. Origin and dissemination of Plasmodium falciparum drug-resistance mutations in South America. J Infect Dis. 2002;186:999–1006.
McCollum AM, Mueller K, Villegas L, Udhayakumar V, Escalante AA. Common origin and fixation of Plasmodium falciparum dhfr and dhps mutations associated with sulfadoxine-pyrimethamine resistance in a low-transmission area in South America. Antimicrob Agents Chemother. 2007;51:2085–91.
Griffing S, Syphard L, Sridaran S, McCollum AM, Mixson-Hayden T, Vinayak S, et al. pfmdr1 amplification and fixation of pfcrt chloroquine resistance alleles in Plasmodium falciparum in Venezuela. Antimicrob Agents Chemother. 2010;54:1572–9.
Pelleau S, Moss EL, Dhingra SK, Volney B, Casteras J, Gabryszewski SJ, et al. Adaptive evolution of malaria parasites in French Guiana: reversal of chloroquine resistance by acquisition of a mutation in pfcrt. Proc Natl Acad Sci USA. 2015;112:11672–7.
Chenet SM, Akinyi Okoth S, Huber CS, Chandrabose J, Lucchi NW, Talundzic E, et al. Independent emergence of the Plasmodium falciparum kelch propeller domain mutant allele C580Y in Guyana. J Infect Dis. 2016;213:1472–5.
Pacheco C, Moreno J, Herrera F. A high number of pfmdr1 gene copies in P. falciparum from Venezuela. Parasitol Res. 2019;118:3085–9.
Lucchi NW, Abdallah R, Louzada J, Udhayakumar V, Oliveira-Ferreira J. Molecular surveillance for polymorphisms associated with artemisinin-based combination therapy resistance in Plasmodium falciparum isolates collected in the State of Roraima, Brazil. Am J Trop Med Hyg. 2020;102:310–2.
Escalante AA, Ferreira MU, Vinetz JM, Volkman SK, Cui L, Gamboa D, et al. Malaria molecular epidemiology: lessons from the International Centers of Excellence for Malaria Research Network. Am J Trop Med Hyg. 2015;93:79–86.
Dalmat R, Naughton B, Kwan-Gett TS, Slyker J, Stuckey EM. Use cases for genetic epidemiology in malaria elimination. Malar J. 2019;18:163.
Pacheco MA, Kadakia ER, Chaudhary Z, Perkins DJ, Kelley J, Ravishankar S, et al. Evolution and genetic diversity of the k13 gene associated with artemisinin delayed parasite clearance in Plasmodium falciparum. Antimicrob Agents Chemother. 2019;63:e02550-18.
Librado P, Rozas J. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics. 2009;25:1451–2.
Kumar S, Stecher G, Tamura K. MEGA7: molecular evolutionary genetics analysis version 7.0 for bigger datasets. Mol Biol Evol. 2016;337:1870–4.
Nei M, Kumar S. Molecular evolution and phylogenetics. New York: Oxford University Press; 2000.
Anderson TJ, Su XZ, Bockarie M, Lagog M, Day KP. Twelve microsatellite markers for characterization of Plasmodium falciparum from finger-prick blood samples. Parasitology. 1999;119:113–25.
Karunaweera ND, Ferreira MU, Munasinghe A, Barnwell JW, Collins WE, King CL, et al. Extensive microsatellite diversity in the human malaria parasite Plasmodium vivax. Gene. 2008;410:105–12.
Imwong M, Sudimack D, Pukrittayakamee S, Osorio L, Carlton JM, Day NP, et al. Microsatellite variation, repeat array length, and population history of Plasmodium vivax. Mol Biol Evol. 2006;23:1016–8.
Schneider KA, Escalante AA. A likelihood approach to estimate the number of co-infections. PLoS ONE. 2014;9:e97899.
Schneider KA. Large and finite sample properties of a maximum-likelihood estimator for multiplicity of infection. PLoS ONE. 2018;13:e0194148.
Chenet SM, Schneider KA, Villegas L, Escalante AA. Local population structure of Plasmodium: impact on malaria control and elimination. Malar J. 2012;11:412.
Eliades N-G, Eliades DG. HAPLOTYPE ANALYSIS: Software for analysis of haplotype data. Distributed by the authors. Forest Genetics and Forest Tree Breeding, Georg-August University Goettingen, Germany. 2009. http://www.uni-goettingen.de/en/134935.html.
Nei M. Analysis of gene diversity in subdivided populations. Proc Natl Acad Sci USA. 1973;70:3321–3.
Efron Bradley, Tibshirani RJ. An introduction to the bootstrap. Chapman & Hall/CRC Monographs on Statistics & Applied Probability; 1994.
Pritchard JK, Stephens M, Donnelly P. Inference of population structure using multilocus genotype data. Genetics. 2000;155:945–59.
Earl DA, vonHoldt BM. STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conserv Genet Resour. 2012;4:359–61.
Jakobsson M, Rosenberg NA. CLUMPP: a cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure. Bioinformatics. 2007;23:1801–6.
Rosenberg NA. Distruct: a program for the graphical display of population structure. Mol Ecol Notes. 2004;4:137–8.
Thomson G, Single RM. Conditional asymmetric linkage disequilibrium (ALD): extending the biallelic r2 measure. Genetics. 2014;198:321–31.
Francisco AP, Bugalho M, Ramirez M, Carriço JA. Global optimal eBURST analysis of multilocus typing data using a graphic matroid approach. BMC Bioinform. 2009;10:152.
Jaramillo-Ochoa R, Sippy R, Farrell DF, Cueva-Aponte C, Beltrán-Ayala E, Gonzaga JL, et al. Effects of political instability in Venezuela on malaria resurgence at Ecuador-Peru Border, 2018. Emerg Infect Dis. 2019;25:834–6.
Pacheco C, Moreno J, Flor Herrera F. Molecular detection and species determination of malaria parasites, Venezuela. Emerg Infect Dis. 2019;25:355–7.
Mekonnen SK, Aseffa A, Berhe N, Teklehaymanot T, Clouse RM, Gebru T, et al. Return of chloroquine-sensitive Plasmodium falciparum parasites and emergence of chloroquine-resistant Plasmodium vivax in Ethiopia. Malar J. 2014;13:244.
Sitali L, Mwenda MC, Miller JM, Bridges DJ, Hawela MB, Chizema-Kawesha E, et al. En-route to the ‘elimination’ of genotypic chloroquine resistance in Western and Southern Zambia, 14 years after chloroquine withdrawal. Malar J. 2019;18:391.
Winter DJ, Pacheco MA, Vallejo AF, Schwartz RS, Arévalo-Herrera M, Herrera S, et al. Whole genome sequencing of field isolates reveals extensive genetic diversity in Plasmodium vivax from Colombia. PLoS Negl Trop Dis. 2015;9:e0004252.
Zhou Z, Griffing SM, de Oliveira AM, McCollum AM, Quezada WM, Arrospide N, et al. Decline in sulfadoxine-pyrimethamine-resistant alleles after change in drug policy in the Amazon region of Peru. Antimicrob Agents Chemother. 2008;52:739–41.
Bacon DJ, McCollum AM, Griffing SM, Salas C, Soberon V, Santolalla M, et al. Dynamics of malaria drug resistance patterns in the Amazon basin region following changes in Peruvian national treatment policy for uncomplicated malaria. Antimicrob Agents Chemother. 2009;53:2042–51.
Pacheco MA, Lopez-Perez M, Vallejo AF, Herrera S, Arévalo-Herrera M, Escalante AA. Multiplicity of infection and disease severity in Plasmodium vivax. PLoS Negl Trop Dis. 2016;10:e0004355.
Fola AA, Harrison GLA, Hazairin MH, Barnadas C, Hetzel MW, Iga J, et al. Higher complexity of infection and genetic diversity of Plasmodium vivax than Plasmodium falciparum across all malaria transmission zones of Papua New Guinea. Am J Trop Med Hyg. 2017;96:630–41.
Koepfli C, Waltmann A, Ome-Kaius M, Robinson LJ, Mueller I. Multiplicity of infection is a poor predictor of village-level Plasmodium vivax and P. falciparum population prevalence in the Southwest Pacific. Open Forum Infect Dis. 2018;5:ofy240.
Escalante AA, Pacheco MA. Malaria molecular epidemiology: An evolutionary genetics perspective. Microbiol Spectr. 2019;7:AME-0010-2019.
Pacheco MA, Schneider KA, Céspedes N, Herrera S, Arévalo-Herrera M, Escalante AA. Limited differentiation among Plasmodium vivax populations from the northwest and to the south Pacific Coast of Colombia: a malaria corridor? PLoS Negl Trop Dis. 2019;13:e0007310.
Bei AK, Niang M, Deme AB, Daniels RF, Sarr FD, Sokhna C, et al. Dramatic changes in malaria population genetic complexity in Dielmo and Ndiop, Senegal, revealed using genomic surveillance. J Infect Dis. 2018;217:622–7.
Chenet SM, Taylor JE, Blair S, Zuluaga L, Escalante AA. Longitudinal analysis of Plasmodium falciparum genetic variation in Turbo, Colombia: implications for malaria control and elimination. Malar J. 2015;14:363.
We thank the DNA laboratory at the School of Life Sciences (Arizona State University) for technical support and Ariana C. Pacheco for the Venezuelan map design. We also thank the patients that volunteered their samples for this study.
This study was supported in part by a grant from the U.S. National Institute of Health, 2U19 AI089681. K. Schneider’s contribution is funded by the German Academic Exchange (Project-ID 57417782), the SMWK-SAB (Project-ID 100257255), and the DFG (Project-ID 656983).
Ethics approval and consent to participate
The study protocol was approved by the Ethics Committee “Complejo Hospitalario Universitario Ruiz y Páez. ”The written informed consents were administered to adult patients who met the recruitment criteria and were procuring treatment. Patients agreed on allowing genetic studies on the parasite.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Additional file 1: Table S1. (A) Complete list of P. falciparum Pfk13 gene sequences available in databases. The codons associated with artemisinin resistance located at the propeller domain are shown. Venezuelan complete sequences (N=65) are labelled in purple. In this study, only the complete sequences were used for the analyses. (B) The mutation found in the Venezuelan samples and the other countries where these mutations were also found.
Additional file 2: Table S2. (A)Plasmodium falciparum data and (B)P. vivax: One-sided bootstrap test for differences in MOI. The test statistic (difference in estimated λ), and p-value of bootstrap test based on B=10,000 bootstrap repeats using bias correction and acceleration are shown. The lower and upper bounds of the 95% bias-corrected and accelerated (BCa) bootstrap condensed interval of the test statistic based on B=10,000 bootstrap repeats are also shown. All markers showed an increase in MOI.
Additional file 3: Table S3. (A)Plasmodium falciparum and (B)P. vivax data: Two-sided bootstrap test for differences in heterozygosity per marker. The test statistic (difference in heterozygosity), and p-value of bootstrap test based on B=10,000 bootstrap repeats using bias correction and acceleration are shown. The lower and upper bounds of the 95% bias-corrected and accelerated (BCa) bootstrap condensed interval of the test statistic based on B= 10,000 bootstrap repeats, as well as the p-value for a permutation test to compare the underlying frequency distributions using the differences in heterozygosity as test statistic, are also shown. The permutation test was approximated using 10,000 permutations. Note, that the permutation test will have relatively low power because of the employed test statistic. Namely, different frequency distributions can have similar heterozygosity values. Only some markers show a difference in MOI.
About this article
Cite this article
Pacheco, M.A., Forero-Peña, D.A., Schneider, K.A. et al. Malaria in Venezuela: changes in the complexity of infection reflects the increment in transmission intensity. Malar J 19, 176 (2020). https://doi.org/10.1186/s12936-020-03247-z
- Drug resistance genes
- Pfk13 gene
- Multiplicity of infection
- Plasmodium falciparum
- Plasmodium vivax
- Transmission intensity