Transmission dynamics of co-endemic Plasmodium vivax and P. falciparum in Ethiopia and prevalence of antimalarial resistant genotypes

Ethiopia is one of the few African countries where Plasmodium vivax is co-endemic with P. falciparum. Malaria transmission is seasonal and transmission intensity varies mainly by landscape and climate. Although the recent emergence of drug resistant parasites presents a major issue to malaria control in Ethiopia, little is known about the transmission pathways of parasite species and prevalence of resistant markers. This study used microsatellites to determine population diversity and gene flow patterns of P. falciparum (N = 226) and P. vivax (N = 205), as well as prevalence of drug resistant markers to infer the impact of gene flow and existing malaria treatment regimes. Plasmodium falciparum indicated a higher rate of polyclonal infections than P. vivax. Both species revealed moderate genetic diversity and similar population structure. Populations in the northern highlands were closely related to the eastern Rift Valley, but slightly distinct from the southern basin area. Gene flow via human migrations between the northern and eastern populations were frequent and mostly bidirectional. Landscape genetic analyses indicated that environmental heterogeneity and geographical distance did not constrain parasite gene flow. This may partly explain similar patterns of resistant marker prevalence. In P. falciparum, a high prevalence of mutant alleles was detected in codons related to chloroquine (pfcrt and pfmdr1) and sulfadoxine-pyrimethamine (pfdhps and pfdhfr) resistance. Over 60% of the samples showed pfmdr1 duplications. Nevertheless, no mutation was detected in pfK13 that relates to artemisinin resistance. In P. vivax, while sequences of pvcrt-o were highly conserved and less than 5% of the samples showed pvmdr duplications, over 50% of the samples had pvmdr1 976F mutation. It remains to be tested if this mutation relates to chloroquine resistance. Monitoring the extent of malaria spread and markers of drug resistance is imperative to inform policy for evidence-based antimalarial choice and interventions. To effectively reduce malaria burden in Ethiopia, control efforts should focus on seasonal migrant populations.


Introduction
Despite considerable progress towards malaria control, two-thirds of the population in Ethiopia, i.e., approximately 66 million people, reside in areas of low or high malaria transmission [1]. Apart from human factors such as population mobility, urbanization, and agricultural development, emergence of drug resistant parasites and insecticide resistance present a major hurdle to malaria control programs in Ethiopia and worldwide [2]. Reports of emerging Plasmodium vivax resistance to chloroquine (CQ) in Ethiopia threaten the efficacy of P. vivax treatment [3][4][5][6]. Also, the well-documented emergence of P. falciparum resistance to artemisinin in Southeast Asia may endanger current malaria treatment programs in Ethiopia, given that both CQ and sulfadoxine-pyrimethamine (SP) resistance originated in Southeast Asia and spread quickly to East Africa [7]. Thus, knowing how malaria parasites spread as well as monitoring prevalence of drug resistant markers in high-risk areas are important to informing antimalarial interventions.
While P. vivax is the most widespread human malaria parasite, it is rare in Africa where P. falciparum predominates. Due to its low prevalence in the continent, little is known about the transmission patterns of P. vivax in Africa. Ethiopia is unique in that P. vivax is co-endemic with P. falciparum at approximately equal case incidence rates. Other African countries with significant P. vivax infections are Eritrea, Sudan, and Madagascar [1]. Although Ethiopia carries a substantial malaria burden, information on the transmission dynamics and spread of drug resistance across the country is scarce.
P. vivax and P. falciparum exhibit different biological and epidemiological features. Compared to P. falciparum, P. vivax has a broader temperature tolerance, an early onset of gametocyte development, and a dormant life cycle stage, the hypnozoite, in the host liver that can cause relapse. Relapse infections may present opportunities for P. vivax to exchange and genetic structures, as well as antimalarial drug efficacy in Ethiopia. Monitoring for markers of antimalarial drug resistance is imperative to informing public health interventions.

Ethics statement
Scientific and ethical clearance was obtained from the institutional scientific and ethical review boards of Jimma and Addis Ababa Universities in Ethiopia and University of California, Irvine, USA. Written informed consent/assent for study participation was obtained from all consenting heads of households, parents/guardians (for minors under age of 18), and each individual who was willing to participate in the study.  Table). This area encompasses an elevation gradient from ca. 50m in the basin to over 2,500m in the highlands west of the Great Rift Valley. Finger-prick blood samples were collected from malaria symptomatic (who has fever with axillary body temperature > 37.5˚C and with confirmed asexual stages of malaria parasite based on microscopy) or febrile patients visiting the health centers or hospitals at each of the study sites. Thick and thin blood smears were prepared for microscopic examination and three to four spots of blood, equivalent to~50 μl, from each individual were blotted on Whatman 3MM filter paper. Parasite DNA was extracted from dried blood spots by the Saponin/Chelex method [34]. Nested and quantitative PCR were performed to identify and confirm parasite species of the infected samples [35]. A total of 226 and 205P. falciparum and P. vivax samples (ranged from 18-58 samples per site) were included in microsatellite analyses.

Data analyses
Linkage disequilibrium and genetic diversity. All analyses were performed separately on the P. falciparum and P. vivax datasets. To examine whether the microsatellite loci represent an independent set of markers of the parasite genome, linkage disequilibrium (LD) was tested by Fisher's exact test for each pair of loci with GenePop version 4.2, using the Markov chain method with 100 batches and 10,000 iterations per batch [39]. Significance values were adjusted by sequential Bonferroni correction for multiple comparisons. In addition, mulltilocus LD was assessed among the parasite samples for each of the populations using the webbased LIAN 3.5 [40]. The standardized index of association (I A S ), which measures the strength of linkage disequilibrium and views as a function of the rate of recombination among samples, was calculated with 10,000 random permutations of the data.
The percentage of polyclonal infections (i.e., samples with more than one allele at any given locus) as well as multiplicity of infections (MOI: the number of genetically distinct clones The first three axes that contain nearly 90% of the total variation are shown. Locations of the studied sites from different parts of Ethiopia are presented in the map below as well as S1 Table. https://doi.org/10.1371/journal.pntd.0005806.g001 present within a host) were estimated for each of the study sites, respectively, for P. falciparum and P. vivax. For each sample, MOI was scored as the maximum number of alleles observed when all loci were taken into account and the average MOI was calculated for each population.
Genotypic variation was calculated in GenoDive version 2.0b27 [41]. We first calculated genetic distances using the method of Smouse and Peakall, a squared Euclidean distance based on the number of times a certain allele was found in the two individuals [42]. The minimal distance class was set as threshold to identify the following: (i) the number of multilocus genotypes (G); (ii) Simpson's diversity index (D), also known as Nei's genetic diversity corrected for sample size that ranges from zero (where two randomly chosen individuals in a population share a single genotype) to one (where individuals have different genotypes); and (iii) genotype evenness (E) that ranges from zero (where one or a few genotypes dominate in a population) to one (where all genotypes are of equal frequency in a population). In addition, the number of effective alleles and expected heterozygosity were estimated for each study site.
Population structure and isolation-by-distance A model-based Bayesian method implemented in STRUCTURE v2.3.4 was performed to examine partitioning of individuals to genetic clusters [43]. The number of clusters (K) was determined by simulating a range of K values from 1 (no genetic differentiation among all sites) to 6 (all sites were genetically differentiated from one another). The posterior probability of each value was then used to detect the modal value of ΔK, a quantity related to the second order rate of change with respect to K of the likelihood function [44]. Posterior probability values were estimated using a Markov Chain Monte Carlo (MCMC) method. A burn-in period of 500,000 iterations followed by 10 6 iterations of each chain was performed to ensure convergence of the MCMC. Each MCMC chain for each value of K was run ten times with the 'independent allele frequency' option that allows individuals with ancestries in more than one group to be assigned into one cluster. Individuals were assigned into K clusters according to membership coefficient values (Q) ranged from 0 (lowest affinity to a cluster) to 1 (highest affinity to a cluster). The partitioning of clusters was visualized with DISTRUCT [45]. Neighboring-joining trees were constructed using T-REX [46,47] to show the genetic relatedness among P. falciparum and P. vivax samples. The squared Euclidean distance, which is based on the number of times a certain allele found in two individuals [48], was used for tree constructions. The resulted trees were visualized in FigTree v1.4.2.
An F ST analysis was conducted using θ, an F ST -estimator in SPAGeDi v1.2e [49]. F ST values were tested for significance using 10,000 permutations. Genetic differentiation among sites was displayed by multidimensional scaling plot based on the estimated D S values (an analog of F ST ) in R v3.3.0. Furthermore, an analysis of molecular variance (AMOVA) was used to determine the hierarchical distribution of genetic variance within and among populations, as well as among regions (north, east, and south of Ethiopia) using GENALEX [50]. The relationships between genetic distances (D S values) and Euclidean geographical distance (estimated from spatial coordinates using R for multivariate and spatial analysis; [51]) were examined by Mantel tests (10,000 randomizations) and reduced major axis (RMA) regression in the Isolation By Distance v3.23 [52].

Bottlenecks and migration rates
Signature of genetic bottleneck was detected with BOTTLENECK v1.2.02 [53]. Only sites with a sample size of 20 or above were included for statistical significance. Two tests were performed using three different mutation models: the infinite alleles model (IAM), the stepwise mutation model (SMM), and a combination of those two extreme hypotheses, the two-phase model (TPM). First was the overall distribution of allele frequency classes. Second was the Wilcoxon-signed rank test to compare the number of loci that present a heterozygosity excess to the number of such loci expected by chance only.
Frequency of gene flow among populations was estimated for each parasite species by a maximum-likelihood analysis implemented in Migrate-N v2.4.4 [54]. Parameters including Θ (defined as 4N e μ, where N e is the effective population size and μ is the mutation rate per generation and site) and M (m/μ, where m is the immigration rate scaled by mutation rate) were estimated. Four independent runs were conducted with the Brownian motion model using 10 short chains with 5,000 sampled genealogies and three long chains with 50,000 sampled genealogies to obtain the mean and range of Θ and M values. In addition, we inferred migration rate using a Bayesian approach implemented in BayesAss v3 [55], which is not dependant on the assumption of equilibrium and can be used with populations that are not in migrationdrift or Hardy-Weinberg equilibrium. A MCMC algorithm was used to estimate the posterior probability distribution of the proportion of migrants from one population to another. We performed the analyses with 9×10 6 iterations, with a burn-in of 10 6 iterations, and a sampling frequency of 2,000 to ensure the parameters of the model were converged. The correlation between migration rate and geographical distance was tested for all pairs of populations.

Landscape genetics
To test for the effects of landscape factors on gene flow between populations, we performed a landscape genetic analysis as follows. First, we created landscape resistance surfaces based on our predictions of the factors influencing gene flow of Plasmodium species, specifically factors influencing vector ecology (land cover and precipitation), human movement (distance to roads), and Plasmodium biology (elevation, which is tightly correlated with temperature). A resistance surface is a spatial layer in which each cell in a grid is assigned a value that represents the degree to which that cell constrains gene flow or movement [56]. These values were often based on numerous assumptions about relationships between a landscape or environmental feature and the ability of a given organism to move through that feature. Landscape resistance surfaces were derived from publicly available data: NASA MODIS MCD12Q2 for land cover (forest, shrubland, woody savanna, savanna, grassland, cropland, and sparsely vegetated) [57,58]; NASA SRTM v4.1 for elevation [59]; WorldClim v1.4 for precipitation [60]; and Roads Africa shapefile in ArcGIS for distance to roads. All raster files were resampled to a resolution of 1km in ArcGIS 10. Next, we used ResistanceGA to optimize landscape resistance surfaces based on our genetic data [61]. ResistanceGA uses a genetic algorithm to unbiasedly assign landscape resistance values to continuous or categorical data. Circuitscape v.4.05 was used to measure resistance distance between populations [62]. Circuitscape relies on electrical circuit theory to predict landscape connectivity and incorporates all possible pathways between populations into the resistance distance measure. To test the fit of resistance surfaces in relation to the genetic data, linear mixed effects models with the maximum likelihood population effects (MLPE) were fit in lme4 [63]. Finally, Akaike information criterion with a penalty for extra parameters (AICc) was calculated from the linear mixed effect model and used as the means for model selection.

Pfmdr1 and Pvmdr1 gene copy estimation
The pfmdr1 gene copy number of P. falciparum was assessed by real-time PCR. Genomic DNA of P. falciparum clones 3D7 (which has a single copy of pfmdr1) was used as a calibrator and pfβ-tubulin, a house-keeping gene, was used as an internal control. The primers for pfmdr1 and β-tubulin were described previously [65]. For P. vivax, the Salvador I strain was used as a calibrator and the pvaldolase gene, which is known to be a single copy gene in P. vivax, was used as an internal control using the published primers [66].
Amplification was performed in triplicate in a total volume of 20 μl containing 10μl of SYBR Green PCR Master Mix, 0.75 μl of each of the sense and anti-sense primers (10 μM), 20 ng of genomic DNA and 3.5 μl of water. Reaction was performed in CFX96 Touch™ Real-Time PCR Detection System (Bio-Rad), with an initial denaturation at 95˚C for 3 min, followed by 45 cycles at 94˚C for 30 sec, 55˚C for 30 sec, and 68˚C for 1 min with a final 95˚C for 10 sec. This was then followed by a melting curve step of temperature ranged from 65˚C to 95˚C with 0.5˚C increment to determine the melting temperature of each amplified product. A negative control with no template was used in each run. Each sample was run in triplicates and the C t values and melting temperature were recorded at the end of the reactions. The average and standard deviation of the three C t values were calculated, and the average value was accepted if the standard deviation was lower than 0.32. The 2 ΔΔCt±SD method for relative quantification was used to estimate the gene copy number [66] and the results were expressed as the N-fold copy number of the targeted gene in relative to the reference. Fisher's exact test (given small sample size) was used to test for significant differences in mutation prevalence and gene copy number among the study populations. All statistical analyses were performed in R (R Core Team 2013).  Table 1). Compared to P. falciparum (8.8%; 20/226; Table 1), P. vivax indicated a lower rate of polyclonal infections (4.3%; 9/205). Polyclonal samples were observed in all sites for P. falciparum, with the highest rate of polyclonal infections in the southern lowlands (HA: 16.7% and JM: 11.8%). For P. vivax, polyclonal infections ranged from 5.3% in Bure (BU) to 0% in Shewa Robit (SR), despite a slightly smaller sample size. Likewise, MOI for P. falciparum from all sites (mean MOI = 1.10; Table 1) was significantly higher than that of P. vivax, (mean MOI = 1.04, P<0.01), indicative of a higher complexity within P. falciparum infections.

Linkage disequilibrium and complexity of infections
Among all the polyclonal infections, 18 were bi-clonal of which two equally dominant alleles were detected in a single locus. We separated the genotypes of the two strains and included them in the analyses. For the 11 samples that showed >1 alleles in two or more loci, we were unable to confidently differentiate the genotypes of the different strains and thus these samples were discarded in the analyses (S3 Table).

Genetic diversity comparison
Both P. falciparum and P. vivax revealed similar levels of allelic and genotypic diversity (Table 2). Nevertheless, genotypic evenness in P. vivax from the highlands (E = 0.75; BU and MA) was significantly lower than the other samples (P<0.01; two-tailed t-test; Table 2). This suggested a less even distribution of genotypes in these populations and that some genotypes were more common than the others. AMOVA indicated that most of the genetic variation was within populations in both P. falciparum and P. vivax (>84%; S4 Table). In P. falciparum, a greater proportion of variation was found among regions (13% among the north, east and south Ethiopia) than among populations in a region (3%); whereas in P. vivax, the proportion of variation among regions (8%) and populations (7%) were comparable and significant.

Genetic clustering of samples
Populations in north and east Ethiopia were slightly differentiated from those in the south (Fig  1). This pattern was shown in both P. falciparum and P. vivax, though populations were more scattered in P. falciparum. The two northern populations (BU and MA) were genetically close to samples in the eastern Rift Valley (SR and ME; Fig 1).
In P. falciparum, three most probable genetic clusters were detected by STRUCTURE analyses (Fig 2), but these clusters did not clearly represent geographical regions. The red cluster was most apparent in the northern (MA and BU) and eastern (SR and ME) populations, but less significant in the southern populations (JM and HA). All three genetic clusters were found in sites BU and ME, but the blue cluster was almost absent in MA and SR (Fig 2). In P. vivax, samples from north and east Ethiopia constituted predominantly the purple cluster, contrast with those from the south that constituted an admixture of the purple and yellow clusters. Unlike P. falciparum, the genetic composition between the northern and eastern populations was largely similar. Neighbor-joining trees did not indicate clear distinction among the P. falciparum and P. vivax population samples (Fig 3). Both trees had relatively short internodes but long terminal branches, which suggested that the parasite lineages were rapidly diverged from one another and that frequent gene exchange occurred among the populations. In P. falciparum, there were a number of subclades where parasites from the same population were genetically closely related (e.g., subclades I, II, and III in Fig 3A). However, in P. vivax, samples from different sites were clustered together in the same clade without clear differentiation (Fig 3B).

Distance and landscape factors
Mantel tests indicated no significant association between geographical and genetic distances among populations of P. falciparum (R 2 = 0.14, P>0.05) and P. vivax (R 2 = 0.09, P>0.05), respectively. These results suggest that parasite gene flow was not limited by geographical distance. Further, for both P. falciparum and P. vivax, we found that none of the tested landscape factors explained pairwise genetic distance (F ST ) among populations more than the Euclidean distance alone based on AICc (Table 3). These results indicated that the differences in land cover, elevation, precipitation, and distance to roads (a proxy for accessibility) did not significantly influence parasite gene flow and that our study populations were clearly connected (S1 Fig).

Demographic change and migrations
All populations of P. falciparum showed a normal L-shape distribution in allele frequency ( Table 4), suggesting that these populations did not experience a recent severe bottleneck. In P. vivax, allele frequency was shown with a shifted mode in site SR (east Ethiopia), indicative of a significant genetic bottleneck. In addition, a significant excess of heterozygosity was  Table. https://doi.org/10.1371/journal.pntd.0005806.g002 Transmission dynamics of P. vivax and P. falciparum in Ethiopia observed in this population as well as other northern and eastern populations (ME and BU) under the IAM and SMM mutation models, suggestive of a deviation in the mutation-drift equilibrium. Based on Migrate-N analyses, both P. falciparum and P. vivax showed a relatively small effective population size (Θ = 0.1-0.6 and 0.17-0.88, respectively; S5 Table), which suggested that the effect of drift was unequivocally as significant as migration. Given that most values of M (m/μ) were >1, the effect of migration (m) was larger than the effect of mutation (μ). For P. falciparum, the effective number of migrants per generation N e m ranged from 0.12-8.58. The greatest migration was observed between the north and east Ethiopia (e.g., from BU to SR: M 30.12 and N e m 8.58) and these north-east migrations were frequent and mostly bidirectional. While migrations between the north-east and south Ethiopia appeared to be less significant, Transmission dynamics of P. vivax and P. falciparum in Ethiopia these migrations in most cases were bidirectional (Fig 4). For P. vivax, the effective number of migrants per generation N e m ranged from 0.12-3.58. The greatest migration was between the northern populations (e.g., from MA to BU: M 22.58 and N e m 3.58) followed by the migration from the north and east to south Ethiopia (e.g., from ME to JM: M 18.77 and N e m 2.02; from MA to JM: M 16.60 and N e m 1.96). Interestingly, the migrations to the south were primarily unidirectional (Fig 4). BayesAss analyses supported the estimates of migration rates from Migrate-N. No significant correlations were found between migration rate and geographical distance in P. falciparum (R 2 = 0.13, P = 0.09) and P. vivax (R 2 = 0.02, P = 0.38; S2 Fig).

Resistance gene marker polymorphisms
Among the 226 P. falciparum and 204 P. vivax samples, we successfully amplified and obtained complete resistance gene data in 199 (88%) and 185 (90%) of the samples (S7 Table). Samples with incomplete data were excluded in the analyses.
By contrast, the pattern of mutation in genes pfdhps and pfdhfr that associated with SP resistance appeared to vary among geographical regions in Ethiopia (Fig 5). For instances, 62% (40/65) of P. falciparum in the northern populations had the mutant 396K of pfdhps, which was significantly higher than that in the eastern (8/72 = 11%) and southern populations  Transmission dynamics of P. vivax and P. falciparum in Ethiopia Amplification and sequencing of the entire pfK13 indicated that this gene was highly conserved among all the P. falciparum samples. No polymorphisms were detected at the codon positions putatively related to artemisinin resistance (Fig 5; S7 Table). For P. vivax samples from the north, east and south Ethiopia, all revealed a similar pattern of mutations in pvcrt-o and pvmdr1, the genes that associated with chloroquine resistance. Sequences of these two genes were highly conserved among the samples. Almost all had the wild type genotype except pvmdr1 codon 976 where over 50% of the samples had the mutant 976F genotype (north: 35/52 = 67.3%; east: 20/38 = 52.6%; south: 52/93 = 55.9%; Fig 5). qPCR data indicated that less than 4% of the samples had two or more copies of the pvmdr1 gene.

Discussion
Ethiopia is a unique malaria endemic country in sub-Saharan Africa where P. falciparum and P. vivax coexist. The present study showed that both species revealed similar population structure. The northern and eastern populations of P. falciparum and P. vivax were genetically closely related and slightly distinct from the southern populations. While parasite gene flow was most frequent between the northern highlands and the highland-fringe areas along the Rift Valley, the southern basin populations were not excluded. We did not find a significant association between any landscape factors and population genetic structure. This result may be caused by the chosen metric for human movement (distance to roads), which does not fully capture the seasonal migration patterns that coincide with harvest season (i.e., September to November). The seasonal migration from highland to lowland areas for agricultural harvest appears to most closely reflect that of the observed gene flow patterns. Since the end of civil war in 1991, seasonal migration of Ethiopians from highland and highland-fringe areas to lowlands has increased due to the growth of large-scale agricultural development projects [67]. Sugarcane production and coffee plantations in the lowlands are important sources of employment and income to people who live in the highlands where agricultural activities are scarce. It is possible that agricultural employments in the eastern Rift Valley and southern basin areas elicit seasonal human migrations and consequently enable spread of P. falciparum and P. vivax  Table. https://doi.org/10.1371/journal.pntd.0005806.g004 Transmission dynamics of P. vivax and P. falciparum in Ethiopia across broad areas without landscape or distance barriers [68,69]. For example, the lowland districts of which populations can increase by 20-30% as a result of the arrival of tens of thousands of farmers during the harvest season [70,71]. Hence, in addition to maintaining control efforts in the community, seasonal migrant populations should not be ignored in order to The proportion of wild type (white) and mutants (gray) for each codon position was indicated in each column of the respective resistant gene associated with the antimalarial drug (bottom bar). GCN denotes gene copy number of pfmdr1 and pvmdr1 (white: single copy; gray: ! two copies). Significant difference was detected in the mutation frequency of codons in pfdhps and pfdhfr among the study sites, as indicated by a small letter above the column. , which is about 50km away from our study site Jimma in southern Ethiopia. The moderate and comparable genetic diversity observed in both species in Ethiopia contrasts with those reported in South Pacific and Southeast Asia where P. vivax showed a higher microsatellite diversity and less fragmented gene pool than the sympatric P. falciparum [9][10][11][12]. In countries or areas where malaria transmission is intense and stable, relapse could be common. This may, in turn, facilitate recombination and local spread of P. vivax leading to higher within-population diversity and lower among-population differentiation compared to P. falciparum. By contrast, our study sites in Ethiopia represent areas of low to moderate transmission settings where the relapse rate of P. vivax is largely unknown [35,75]. Malaria there displays a strong seasonal pattern with a lag time varying from a few weeks at the beginning of the rainy season to more than a month at the end of the rainy season [76]. The rainy season is relatively short in the highland and highland-fringe areas, and thus transmission season is usually short-lived [77]. These apparent seasonal and/or landscape differences not only influence the behavior and distribution of vector mosquitoes and the length of the parasite life cycle [78], but also the patterns of human movement and settlement that can in turn determine transmission dynamics of malaria. For instances, although P. vivax infections can occur periodically throughout the year, human migration is seasonal in Ethiopia and this could limit the spread of relapse P. vivax. Also, the long arid or semi-arid season particularly in the highlands may constrain gametocyte development of P. vivax or local transmission even when relapse occurs. A recent study in southern Ethiopia indicated an approximately 9.4% (ranged from 6.4-13.6% by sites) of P. vivax patients showed recurrent infection by day 28 [5,6]. However, because relapses can occur as early as 21 days following initial treatment, it is unclear how many of the recurrent cases were due to relapse. Given that our P. falciparum and P. vivax samples were collected at the same time during the peak transmission season, the lack of contrast between P. falciparum and P. vivax population structure and diversity may suggest similar demographic factors and a less significant impact of relapse. Future study should investigate the incidence of relapse among transmission settings and how such influences parasite population diversity.
The overall polyclonal infection rates observed in P. falciparum (8.8%) and P. vivax (4.3%) were low in our study area. The proportion of polyclonal P. vivax infections was similar to that reported in areas of low endemic setting such as Central China (2-19%) [79], but considerably lower than hypo-endemic areas such as Vietnam (71.4%; [80]), Sri Lanka (68%; [81]), Colombia (60-80%; [82]), and the Amazon Basin in Brazil (50%; [83]). Likewise, the polyclonal rate of P. falciparum was comparable to low transmission areas such as Haiti (12.9%) [72] and southern China (10-23%) [84] that are approaching elimination phase, but lower than endemic countries in West Africa such as Gambia and Senegal (36-50%) [85,86], East Africa such as Kenya (70-90%) [87], Papua New Guinea (39-45%) [88], and Southeast Asia such as Malaysia (65%) [89] and Cambodia (47%) [9].The higher polyclonality in P. falciparum than in P. vivax may imply a potentially large P. falciparum reservoir present in asymptomatic hosts that emerges during the transmission/rainy season [90]. By contrast, because P. vivax infections can occur periodically throughout the year in Ethiopia, our samples may reflect only a fraction of the existing P. vivax gene pool. Although relapse and early production and circulation of gametocytes in P. vivax infection can lead to increased opportunities for recombination, drug-sensitive clones (both blood-stage parasites or hypnozoites) could be eliminated during the non-transmission season by antimalarial treatment and resulted in reduced withinhost diversity [8,82]. The lack of demographic data and drug use history of our patients limits our exploration of immunity and other factors on the observed MOI. It is important to note that microsatellites have the potential to underestimate polyclonality due to their lower sensitivity and specificity in detecting minority alleles compared to amplicon deep sequencing [91]. Despite the concern of relapse by P. vivax, the asymptomatic P. falciparum and P. vivax reservoirs that remain undetected during non-transmission season or in low endemic areas could pose a long-term impact on local transmission. Ethiopia adopted AL as the first-line treatment for uncomplicated P. falciparum in 2004 in response to increased resistance to CQ and SP [2]. Unlike the Greater Mekong Subregion of SE Asia where delayed ACTs response has been reported, AL is shown to be highly effective in clearing parasite and fever within three days of drug treatment in Ethiopia [92]. The high efficacy of AL concords with our finding of predominantly wild type pfK13 genotypes in P. falciparum from the northern, eastern, and southern populations. Although CQ and SP have not been used for P. falciparum treatment in the last decade, the high prevalence of mutations in pfcrt 76, pfmdr1 184, pfdhps 396, and pfdhfr 51 and 108 across Ethiopia suggested that strong selection may still exist possibly by the use of CQ for treating P. vivax malaria, as well as SP for intermittent preventive treatment (IPT) on pregnant women as part of the antimalarial schemes in sub-Saharan Africa [93]. It is concerning that the high prevalence of SP resistance mutations observed in the present study may indeed influence the outcome or effectiveness of IPT [94]. Another explanation for the high frequency of CQ and SP resistance mutations is the spread of resistance parasites from one population to another. Microsatellites indicated a substantial admixture of P. falciparum genotypes between north and east Ethiopia. It is possible that resistant parasites can spread via frequent human movements and become locally selected.
The observations of an alarming level of CQ resistance prevalence in Papua New Guinea [95], India [96], SE Asia [97,98], and South America [99] after a decade-application deepen the concern for the appearance of CQ resistance to P. vivax in Ethiopia. Recent reports of therapeutic failure of CQ (5.76-13%) [3][4][5][6] as well as high rates (9-32%) of recurrent infections subsequent to CQ usage in different parts of Ethiopia [6,100] have threatened the efficacy of P. vivax treatment. The highly conserved sequences of pvcrt-o as well as the predominantly wild type T958M, F1076L, and single copy of pvmdr1 suggest that these attributes may not be relevant to CQ resistance. By contrast, we detected a high prevalence of pvmdr1 976F mutation, which is considerably higher than that reported in India (22%) where the observed resistance genotypes were confirmed by in vitro drug sensitivity testing [101]. Thus, one possible explanation for the high 976F prevalence in our study area is that this mutation may associate with emerging CQ resistance. Such association merits further clinical observations and/or in vitro testing to confirm its functional significance. Moreover, it is yet unclear whether our P. vivax samples were relapse or recrudescent infections. In Cambodia, 89% of the P. vivax samples that were isolated from patients with recurrent/relapse infections within a 42-day follow-up had pvmdr1 976F mutation [102]. This reiterates the importance of distinguishing recrudescent from relapse infections in order to clarify the implications of the observed mutations and accurately elucidate resistance prevalence of P. vivax. Given that chloroquine monotherapy has been the recommended regimen for P. vivax malaria in Ethiopia for the past decades [2], selection as well as parasite gene flow may explain the emergence and spread of resistance genotypes across the country. Alternative P. vivax treatment regimes such as ACT or CQ in combination with primaquine are suggested to prolong the efficacy of CQ and prevent/reduce relapse in Ethiopia.
In summary, P. falciparum and P. vivax revealed moderate levels of genetic diversity and similar population structure in Ethiopia. Human migrations may promote parasite gene flow while environmental heterogeneity and geographical distance did not appear to be a major gene flow barrier. To effectively reduce malaria burden in Ethiopia, control efforts should focus on seasonal migrant populations. Unconstrained parasite gene flow may partly explain similar patterns of resistance marker prevalence across the country. Our findings are paramount to monitoring the emergence and spread of antimalarial drug resistance and offer evidence-based guidelines to existing treatment regimes.
Supporting information S1