Clinical malaria incidence following an outbreak in Ecuador was predominantly associated with Plasmodium falciparum with recombinant variant antigen gene repertoires

Defining varcodes circulating in Ecuador

To define distinct varcodes circulating locally in Ecuador, we estimated parasite genetic relatedness by calculating the similarity index, Pairwise Type Sharing (P_TS) and a Bayesian version of the same (Bayesian pairwise type sharing; BP_TS, see Methods). BP_TS estimates include statistical uncertainty due to differences in isolate repertoire size (i.e. the number of DBLα types identified in each isolate repertoire), enabling a rigorous examination of relatedness. The median repertoire size in Ecuadorian isolates was 37 DBLα types (range: 11-43 types, Table S1), which is in line with the expected number of var genes in the Honduran laboratory reference strain HB3 (n=42 var genes excluding var2csa, ³³). Nonetheless, we used a uniform prior distribution for repertoire sizes between 40 and 50 types, which accounted for the possibility of missing data in total isolate repertoire sizes using degenerate primers for var DBLα PCR. For a further analysis to validate our PCR and sequencing methodology, including a further sub-classification of DBLα types into upsA and non-upsA groups and their respective isolate genomic proportions, see Methods and Supplementary Text 1.

We constructed genetic relatedness networks at a threshold of P_TS ≥ 0.90 to identify putatively identical genomes within the margin of error of detection of all DBLα types in an isolate. We confirmed this definition by comparing varcodes derived from P_TS to those derived from posterior mean BP_TS estimates (Supplementary Fig. S3), and by examining the lower and upper bounds of the 95% highest density posterior intervals (HDPIs). This revealed nine genetically distinct varcodes in this study with varcode4, varcode5, varcode9 seen only once and at the other extreme 36 isolates had varcode1 (Fig. 3a and Supplementary Fig. S4a). We were unable to confirm whether any of the Ecuadorian varcodes were circulating pre-outbreak as samples were not available. As expected, the outbreak varcode1 was clonal (indicated by HDPIs that included 1; Fig. S5d), as previously demonstrated by microsatellite genotyping⁸. It is worth noting that even though we only sampled 30 P. falciparum isolates in Esmeraldas City from 140 (21%) outbreak cases in 2013, the fact that the disease outbreak was clonal and represented 91% of the 154 total cases on the northwest coast and Amazon region during 2013, means that we characterized a representative sample of the cases from the outbreak and during that year.

• Download figure
• Open in new tab

Figure 3. Reconstructing the disease transmission dynamics of P. falciparum reveals signatures of local transmission and sexual recombination after the outbreak.

(a) A network visualization of the genetic relatedness of P. falciparum isolates at the threshold of PTS ≥0.90 to define varcodes (see Methods). Every node represents a P. falciparum isolate and an edge represents the P_TS value between two particular nodes/isolates. Isolates that cluster together (i.e., connected by edges) are considered to have the same varcode. Each color represents a different varcode. Two P. falciparum isolates belonging to varcode1 appear as outliers in the network due to undersampling of their DBLα types. N refers to the number of isolates. For comparison to BPTS estimates see Supplementary Fig. S3-5 and S8. (b) The number of P. falciparum isolates within each varcode (i.e. size of circle) in each year and their persistence over time. The three varcodes (varcode1, varcode2, varcode3) identified in 2013 were identified again in 2014, and varcode1 was identified again in 2015. Two varcodes identified in 2014 (varcode6, varcode7) were also identified in 2015. (c) A network visualization of the genetic relatedness of varcodes at the threshold of PTS ≥0.50 to discriminate between recombinant varcodes and genetically distinct varcodes. Isolates/varcodes that cluster together (i.e., connected by edges) represent recombinants. (d) A clustered heatmap showing the genetic profiles of each P. falciparum isolate with rows representing each isolate and columns representing each DBLα type. Black and white denote the presence and absence of each type, respectively. Isolates that clustered together were more genetically similar (i.e., the same DBLα types were present). Similarly, varcodes that clustered together were more genetically similar (e.g. recombinants). (e) The spatial distribution of varcodes in each year, with the size of the circle representing the number of P. falciparum isolates sampled in a given location in each year and the pie chart depicting the proportion of each varcode identified. (f) Spatiotemporal genetic relatedness network of varcodes between 2013-2015. Every node represents a sampling location, the size of the circle represents the number of P. falciparum isolates sampled in each location, the pie chart depicts the proportion of each varcode identified, and the weighted edges show genetically-related varcodes (P_TS ≥0.50). (Online version: GIF animation depicts the spatiotemporal transmission dynamics chronologically during and after the outbreak).

For each pairwise comparison, we calculated the HDPI width as an estimate of uncertainty (analogous to calculating the width of a frequentist confidence interval) and found there was high confidence (HDPI width ≤0.2) for 94.6% (632/668) of pairwise BP_TS estimates (Supplementary Fig. S4b). Low confidence pairwise BP_TS estimates (HDPI width>0.2; max = 0.33) were found only between the two varcode1 isolates with the smallest repertoire sizes (11 and 19 DBLα types; 36/36 comparisons). This lower statistical confidence follows naturally from under sampling of those isolates’ DBLα types. After aggregating the within-varcode pairwise comparisons, there was high confidence in the definition of all varcodes (Supplementary Fig. S5a-b).

For comparison of relatedness by microsatellite genotyping, we calculated Pairwise Allele Sharing (P_AS) using previously published microsatellite data for the same isolates^{8, 50}. Not surprisingly, 7 microsatellites did not resolve the nine varcodes but only identified four, three of which were single parasite isolates, at the threshold of P_AS ≥ 0.90 (Supplementary Fig. S6a). Because two of the seven microsatellite markers were fixed in the population^{8, 50}, we also generated the relatedness networks at the threshold of P_AS ≥ 0.80, which allowed for an error of detection of 1 microsatellite allele. However, the majority of isolates clustered at this threshold and only three varcodes were resolved (Supplementary Fig. S6b). There were several pairwise comparisons (14%, 72/541 by P_TS and 6%, 31/541 by BP_TS) where P. falciparum isolates were found to be identical by microsatellites (P_AS = 1) but much less related by var (P_TS or BP_TS ≤ 0.50, Supplementary Fig. S7). Nevertheless, the relationship between P_TS and P_AS scores was positively correlated (Pearson’s correlation coefficient = 0.757, p < 0.001, Supplementary Fig. S7a). Similarly, the relationship between BP_TS and P_AS estimates was also positively correlated (Pearson’s correlation coefficient = 0.779, p < 0.001, Supplementary Fig. S7b).

Persistent association of clinical cases with the same varcodes over time and large distances

We next explored whether parasites with the same varcodes caused clinical episodes after the outbreak. Persistent disease transmission over time (Fig. 3b) and large distances (Fig. 3e-f) was found as the same varcodes were identified post outbreak and sometimes in different sampling locations. The median time between first and last identification of the same varcodes with any clinical case was 216 days or approximately 7 months (range = 190 – 823 days) during the study period. For example, the outbreak varcode1 identified in Esmeraldas City was also identified in San Lorenzo, Esmeraldas (∼150km from Esmeraldas city) in 2013, then Cascales, Sucumbios (>300km away) in 2014, and then in Tobar Donoso, Carchi (∼150km away) in 2015 (Fig. 3b). This is in line with an earlier study using microsatellites that showed that P. falciparum isolates from the outbreak were genetically related to an Ecuadorian parasite isolate collected around 1990⁸. This result was also consistent with reports in Peru⁵¹ and Colombia^{18, 52} where identical parasites were shown to be circulating up to five and eight years after first identification, respectively.

Spatiotemporal genetic relatedness networks reveal signatures of recombination

We next reconstructed spatiotemporal genetic relatedness networks to explore signatures of local disease transmission and parasite microevolution after the outbreak. We applied a P_TS threshold of ≥0.50 to detect recombinant and/or genetically-related varcodes. By using the date and location where each isolate was collected, we animated the network chronologically to examine spatiotemporal disease transmission dynamics during the outbreak in 2013 and post-outbreak in 2014 and 2015 (Fig. 3c, 3f, GIF). This analysis resolved parasite microevolution relevant to disease transmission after the outbreak and revealed that four varcodes (varcode3, varcode4, varcode6, varcode7) were recombinants of the outbreak varcode1 (Fig. 3c). The remaining varcodes did not cluster in the network and were genetically distinct. These patterns were confirmed using BP_TS estimates based on posterior means (Supplementary Fig. S8a) and their corresponding 95% HDPIs (Supplementary Fig. S5c-f and S8b-c). For the varcodes that were recombinants and/or genetically-related (i.e., cluster in the networks in Fig. 3c and Supplementary Fig. S8a), 70.9% (433/611) of the HDPIs were confidently above the threshold of BP_TS ≥0.50 (Supplementary Fig. S5d). For varcodes that were genetically-distinct (i.e., do not cluster in the networks), 98.4% (368/374) of the HDPIs were confidently below the threshold of BP_TS ≤0.50 (Supplementary Fig. S5f). As expected, the relatedness networks based on the 7 microsatellites did not resolve these patterns since all isolates were connected at the threshold of P_AS ≥ 0.50 and genetic relatedness was overestimated (Supplementary Fig. S9a). Indeed, we were not always able to discriminate recombinants from genetically distinct parasites since several pairwise comparisons (25%, 307/1232 by P_TS and 16%, 196/1232 by BP_TS) that resulted in P_AS ≥ 0.50 were found to be less related by var (P_TS or BP_TS ≤ 0.50) (Supplementary Fig. S7). Furthermore, the majority of recombinants were classified as essentially clonal since they clustered together at the threshold of P_AS ≥ 0.80 (allowing for differences in one microsatellite allele, Supplementary Fig. S6b).

We were able to detect that varcode3 was in fact a recombinant of the outbreak varcode1, even though it was sampled in San Lorenzo during the same time period as the ongoing outbreak in 2013 (Fig. 3c, 3e). With the available data we were unable to confirm whether varcode3 was ‘imported’ to San Lorenzo from Esmeraldas City. The most plausible explanation is that the recombination event occurred in San Lorenzo during the outbreak. Parasites with recombinant varcode3 caused clinical episodes in multiple locations in 2014 (Fig. 3e). The other recombinant varcodes (varcode4, varcode6, varcode7) were identified after the outbreak in 2014 or 2015. Whether the recombinant varcodes resulted from sexual recombination events between outbreak varcode1 and genetically distinct parasites that were already circulating at low levels and/or in asymptomatic reservoirs in Ecuador (e.g. ⁵) or those that were previously imported could not be ascertained from the current study population. Interestingly, all the recombinant varcodes were identified in San Lorenzo indicating this area is a transmission hot spot and a reservoir of parasites with diverse var repertoires (Fig. 3e).

To further visualize the genetic profiles of varcodes and confirm the observed recombination signatures, we constructed a clustered heatmap based on the presence and absence of the 195 DBLα types across all isolates, such that isolates with similar genetic profiles cluster together (Fig. 3d). This analysis confirmed that in the case of isolates identified as recombinants of the outbreak varcode1 (i.e., varcode3, varcode4, varcode6, varcode7), a proportion of outbreak DBLα types as well as different DBLα types were present. This is consistent with the generation of a new varcode through outcrossing between two genomes with different var repertoires by conventional meiosis (i.e., approximately ≥50% sharing of DBLα types). By contrast, the genetic profiles of varcodes 2, 5, 8, and 9 had different DBLα types, leading to the definition as parasites with genetically distinct varcodes. When using 7 microsatellites, these recombination signatures in the genetic profiles were not identified (Supplementary Fig. S9b).

Parasites with recombinant varcodes most frequently caused P. falciparum clinical episodes following the 2012-13 outbreak

Prior to the outbreak there had been a steady decline in clinical cases, however, increased incidence of disease occurred after the outbreak. Therefore, we analyzed trends in the epidemiology of P. falciparum cases that occurred post-outbreak to understand if there were any key risk factors. Of the 25 individuals in our study that had a clinical P. falciparum episode after the outbreak in 2014 or 2015, we had age data for 18 individuals (72%, Table 1). There was no significant difference (Kruskal-Wallis test, p = 0.65) in the median age of individuals experiencing clinical episodes caused by parasites with the outbreak varcode1 (19 years, range = 19 – 57 years, N = 3 patients), a recombinant of varcode1 (i.e., varcodes 3, 4, 6, or 7) (25 years, range = 17 – 58 years, N = 13 patients), or by a different varcode (34.5 years, range = 32 – 37 years, N = 2 patients). A greater diversity of varcodes (1-9 inclusive) was associated with incidence of clinical malaria post-outbreak than during the 2013 outbreak (varcodes 1,2,3). Indeed in 2014, 79% of thecases sampled were caused by either parasites with the outbreak varcode1 (21%) or with any of the four recombinants of varcode1 (58%). The trend was similar in 2015 with 83% of the cases sampled caused by either parasites with the outbreak varcode1 (33%) or with any of the four recombinants of varcode1 (50%), although our sampling of reported cases in 2015 after January was limited. Importantly, overall we found that 80% of the cases we sampled after the outbreak were caused by either parasites with the outbreak varcode1 (24%) or with any of the four recombinants of varcode1 (56%), especially with varcode7. Thus, disease transmission after the outbreak in 2014 and into 2015 was predominantly associated with parasites with recombinant var repertoires of the outbreak clone.

Comparative analyses to South America: elucidating possible origins of Ecuadorian varcodes

We next examined the possible origins of the varcodes circulating locally in Ecuador by comparing to published data²⁸ from 128 P. falciparum isolates collected from neighboring South American P. falciparum populations (Fig. 4a). When we combined the sequences from Ecuadorian isolates with those previously sequenced, we identified 543 unique var DBLα types in South America (compared to 458 in ²⁸). This is largely representative of the var diversity circulating in South American P. falciparum populations, as indicated by sampling accumulation curves approaching saturation (Supplementary Fig. S10). Varcoding resolved a total of 97 varcodes in South America (P_TS ≥ 0.90), with those identified in each country ranging from 9 in Ecuador and Venezuela, to 56 varcodes in French Guiana (Supplementary Fig. S11). The number of var DBLα types per isolate repertoire was low in all countries and indicative of only one P. falciparum genome infecting an individual (i.e. repertoire size ≤ 60, Fig. 4b, Supplementary Table S1). We found that the median number of var DBLα types in each isolate ranged from 36.5 in Venezuela to 48.0 in French Guiana and the maximum number of var DBLα types per isolate was 42 or 43 in all countries except French Guiana that appeared to have more multi-genome infections (max = 92). The median repertoire size was in the range of the number of var genes seen in the genome of the Honduran laboratory reference strain HB3 (N=42 excluding var2csa, ^{32, 33}). An unrooted phylogenetic neighbor-joining tree revealed distinct clusters of genetically-related isolates in South America generally clustering by country (Fig. 4c), consistent with previous analyses in the region demonstrating geographic population structure²⁸.

• Download figure
• Open in new tab

Figure 4. Relatedness networks in South America provide insights into the origins of the Ecuadorian varcodes and diversity patterns in South America.

(a) A map showing the study sites across South America. The dates of sample collection in these countries occurred from 2002 to 2008, around five to thirteen years prior to the sample collection of the Ecuadorian isolates. N refers to the number of isolates. (b) Repertoire sizes are shown for each P. falciparum isolate, from which multiplicity of infection can be inferred (i.e., ∼40-60 var DBLα types per P. falciparum genome). (c) An unrooted neigh14bor-joining tree shows the relatedness patterns among South American P. falciparum isolates. (d) A network visualization of the genetic relatedness of South American varcodes at the threshold of PTS ≥0.50. Due to the differences in time points and sample collections, for this analysis a genetically-related parasite did not necessarily reflect a recent transmission event. Every node represents a P. falciparum isolate and an edge represents the PTS value between two particular nodes/isolates. Isolates/varcodes that cluster together (i.e., connected by edges) are genetically-related. (e) Relatedness of Ecuadorian varcodes to P. falciparum populations from South America. The relatedness of each varcode was measured by first concatenating all possible var DBLα types that were identified in the P. falciparum isolates comprising each varcode as well as concatenating all var DBLα types that were identified in each country. Then PTS was calculated between each varcode and each country. The color gradient denotes the PTS value for a particular comparison, with darker shades showing higher relatedness and the different colors corresponding to the different country comparisons. N refers to the number of var DBLα types. (f) A clustered heatmap showing the patterns of diversity across South America, with a row representing all the types identified in a particular country, and columns representing each of the 543 DBLα types. Black and white denote the presence and absence of each type, respectively. The total number of unique var DBLα types identified in each country ranged from 112 in Colombia to 249 in French Guiana. Countries that cluster together share more var DBLα types.

To examine whether any of the parasites with Ecuadorian varcodes were genetically related to other South American parasites, we constructed regional spatiotemporal relatedness networks and applied our threshold of P_TS ≥ 0.50 (Fig. 4d). We identified a genetically-related Peruvian P. falciparum isolate that clustered with the outbreak and recombinant varcodes (P_TS = 0.66-0.75 with varcode1 isolates). This is consistent with previous analyses using microsatellites showing the outbreak source was possibly a residual parasite lineage circulating in Peru in 1999-2000 and in Ecuador in the early 1990s^{8, 13}. The fact that these isolates did not cluster with anything else suggests the outbreak may have been caused locally due to circulating parasites in unidentified reservoirs. Varcode5, which was sampled in the Amazon had a very different genetic profile to the other varcodes (≥83% unique types, Fig. 3d). It is worth noting that the parasite isolate with varcode5 also had a different genetic background with respect to its microsatellite haplotype as compared to the rest of the Ecuadorian isolates (3 of 7 unique microsatellite alleles). The epidemiological data for the putative infection location was recorded as possibly imported from Peru (EC40, Table 1). However, our analysis would suggest that this varcode is more related to P. falciparum populations from French Guiana and Venezuela. Nonetheless, Peruvian isolates also clustered with P. falciparum isolates from French Guiana even though they did not cluster with varcode5 and we found that 53% of the types identified in varcode5 were also identified in Peruvian parasites (Fig. 4e). We were able to determine a putative importation of a Colombian parasite (varcode9, EC53, Table 1, Fig. 4d). The varcodes 2 and 8 did not cluster with any other South American isolate suggesting parasites with these varcodes may be local to Ecuador. The epidemiological data for the putative infection location for varcode8 was recorded as Colombia (EC49 and EC50, Table 1). Although it did not cluster with Colombian P. falciparum isolates in our network analysis, we found that 67% of the types identified in varcode8 (collapsing all 3 isolates) were also identified in Colombian parasites (Fig. 4e). For varcode2, we found that 44% of the types identified in varcode2 were also identified in Colombia (Fig. 4e), providing evidence that parasites with varcode2 and varcode8 may represent residual parasites previously or recently imported from Colombia.

Conservation of individual DBLα types in space and time

Of significance, we found the same var DBLα types (56% of types identified in Ecuador, N=110/195, Fig. 4f) and varcodes that were genetically related to other South American parasites and sampled up to 13 years before (Fig. 4d-e) were contributing to local transmission in Ecuador. Indeed, conservation of the DBLα types in all varcodes was observed in space and time to varying degrees (Fig. 4e). We also found conservation of DBLα types on a continental scale with the number of the same types identified in any two or three countries being 124 (23%) and 59 (11%), respectively (Fig. 4f). Moreover, even DBLα types that were identified in relatively low frequencies in Ecuador (i.e., identified in a few P. falciparum isolates) were also identified in the South America dataset, often at relatively high frequencies in other countries (Supplementary Fig. S12). These patterns are consistent with some parasite varcodes being putative importations and/or highly-related to parasites in other countries, as described previously.

No evidence of recombination in any DBLα type over the course of the outbreak

The fact that the outbreak varcode1 was clonal provided a unique opportunity to examine the conservation of its 47 individual DBLα types after transmission to multiple human hosts over the course of the outbreak and up to two years after (N=36 infections). DBLα types are very variable, with the average pairwise identity of the sequences encoding this domain being approximately 42%³³. The cut-off of 96% sequence identity in our pipeline has been routinely used to identify the same DBLα types within the limits of sequencing errors and PCR artifacts (^{27, 53}, Methods). Our data show that the same DBLα types were conserved in isolates with varcode1 over the course of the outbreak (N=30 infections) and up to two years after the outbreak in the other six cases with varcode1 (Supplementary Fig. S13a). In total, we report 1,283 independent observations of the same 41 DBLα types identified in 36 cases with varcode1 assessed over 2 to 3 years (any type was identified in a median of 34 infections, ranging from 1 to all 36 infections) (Supplementary Fig. S13a). Given our threshold of P_TS ≥0.90 to identify putatively identical var DBLα repertoires, only 10 of the 47 DBLα types found in the outbreak clone were not found in the majority of the isolates with varcode1 (median frequency = 1, max = 6, Supplementary Fig. S13a). Four of these were identified in the cases with recombinant varcodes (median frequency = 17.5, range = 6-19), two of which were also identified in the Peruvian isolate that clustered with the outbreak/recombinants in Fig. 4d, suggesting they are also conserved. The other 6/10 were singletons in the entire dataset (including the South America data) and could be a consequence of the PCR with degenerate primers and/or low-quality DNA.

Among the 15 cases with recombinant varcodes (varcodes 3,4,6, and 7), we report 339 independent observations of 39 of the 47 outbreak DBLα types. There were 26-27 outbreak types shared between varcode1 and the recombinant varcodes 3,4, and 7, and 17 types between varcode1 and varcode6 (Supplementary Fig. S13b). These findings highlight the conservation of DBLα types of the outbreak clone in space and time, a pattern further supported by the fact that 31 of 47 outbreak DBLα types were found among previously recorded South American isolates (in 2002 to 2008²⁸; see Supplementary Fig. S13c). This evidence points to conservation of these types despite the many transmission events that must have occurred between 2002 and 2013-15 across South America. We found that the highest number of outbreak DBLα types were identified in Peru (N=28, 60%) followed by Colombia (N=14, 30%) and Venezuela (N=8, 17%), and the lowest in French Guiana (N= 2, 4%, Supplementary Fig. S9c).

Ethics statement

The P. falciparum isolates collected from individuals of all ages presenting with uncomplicated malaria cases were obtained from the malaria surveillance protocol approved by the Ethical Review Committee of Pontificia Universidad Católica del Ecuador. Ethical approval for this study was obtained from Pontificia Universidad Católica del Ecuador (Quito, Ecuador) and University of Melbourne (Melbourne, Australia). Written informed consent was provided by study participants and/or their legal guardians.

Study design, study area and sample collection

In this molecular epidemiological study, we examined samples that were collected passively from 2013 to 2015 in endemic areas of Ecuador, namely the northwest coast and the Amazon region, as part of ongoing malaria surveillance by the Ecuadorian Ministry of Health and patients showing at the local clinics. During 2013, the majority of samples were collected from an ongoing outbreak that occurred in Esmeraldas city, Esmeraldas province. Details of the sample collection are described in Sáenz et al⁸. In 2014 and 2015, the “post-outbreak” samples were mostly collected from across the northwest coast with few from the Amazon region, reflective of the epidemiology of clinical P. falciparum cases in Ecuador (Northwest coast: 20 and 173 cases; Amazon: 11 and 16 cases in 2014 and 2015, respectively). Details of the sample collection are described in Vera-Arias et al ⁵⁰. It is worth noting that although cases were reported throughout 2015, samples were not available for many of them due to logistical challenges. Briefly, subjects who were diagnosed with falciparum malaria either by light microscopy or a rapid diagnostic test were asked if they wanted to participate in the study and signed an informed consent form. The participants had to: 1) live in the areas where the samples were taken, 2) be over 2 years old, 3) agree to participate, 4) sign an informed consent form, 5) give a blood sample, 6) be able to understand the informed consent. From each participant a blood sample (either venous blood or dried blood spot) was collected as well as their responses to a basic demographic questionnaire (including places recently travelled and their address). The positivity of all samples was confirmed in the laboratory by microscopy and PET-PCR ⁷⁶.

DNA extraction

Genomic DNA was extracted from venous blood or dried blood spot samples using a QIAamp DNA MINI KIT (QIAGEN, Germantown, MD, USA) using the protocol recommended by the manufacturer.

Microsatellite genotyping

The same P. falciparum isolates in this study were previously genotyped for seven putatively neutral microsatellite markers (TA1, POLYα, PfPK2, TA109, 2490, C2M34, C3M39) as described in^{8, 50}. These previously published microsatellite genotyping data^{8, 48, 75} were used in the present analyses.

Var genotyping and related data processing

For var genotyping, the DBLα domains of antigen-encoding var genes were amplified and sequenced on a MiSeq 2×300bp paired-end Illumina platform as described in ^{29, 60}. We obtained 6,929,470 raw sequence reads, which we cleaned and processed using the DBLaCleaner pipeline (⁶⁰, http://github.com/Unimelb-Day-Lab/DBLaCleaner). Our customized bioinformatic pipeline has been described in detail in⁶⁰. Briefly, we de-multiplexed and merged the reads as well as removed low-quality reads and chimeras using several filtering parameters (see Supplementary Fig. S1 for details). This resulted in 2,141 cleaned DBLα sequences (Supplementary Fig. S1). To identify distinct or unique DBLα types (i.e., unique genetic variants), we clustered the DBLα sequences from Ecuadorian P. falciparum isolates with 5,699 previously published²⁸ DBLα sequences from other South American countries (Colombia (N=21 isolates), French Guiana (N=76 isolates), Peru (N=21 isolates), and Venezuela (N=10 isolates) at the standard 96% sequence identity²⁷ using the clusterDBLalpha pipeline (http://github.com/Unimelb-Day-Lab/clusterDBLalpha). All the cleaned DBLα sequences in this study have been submitted to the DDBJ/ENA/GenBank (BioProject Number: PRJNA642683).

We further curated our dataset by translating the DBLα types into amino acid sequences using the classifyDBLalpha pipeline (http://github.com/Unimelb-Day-Lab/classifyDBLalpha) and removing any DBLα types that were non-translatable (N=4). In addition, any P. falciparum isolate with low sequencing quality (< 10 DBLα types) was removed from the analysis. Varcoding was completed for a total of 70 P. falciparum Ecuadorian isolates, however, 12 P. falciparum isolates with low DNA quality and/or poor sequencing quality were removed and we obtained var DBLα data for 58 isolates (82.9%). A total of 543 unique DBLα types identified in the 186 South American P. falciparum isolates (N = 58 Ecuadorian P. falciparum isolates from this study, N = 128 South American P. falciparum isolates previously published in ²⁸) were used for subsequent var analyses at the regional-level, and only the 195 unique DBLα types identified in Ecuador were used for Ecuador-specific analyses.

A further validation of our sequencing and genotyping methodology was undertaken by analyzing the isolate genomic proportions of DBLα types classified as upsA and non-upsA using our classifyDBLalpha pipeline²⁹ and are described in Supplementary Text 1.

Genetic relatedness analyses

Statistical analysis

We used R version 3.5.2 for all analyses⁸⁸. The package dplyr⁸⁹ was used for further data curating. We used base R and the R package epiR⁹⁰ to perform chi-squared tests for univariate analyses of categorical variables to compare proportions and for non-parametric tests to compare distributions of continuous variables between two groups (Mann-Whitney U test) or among k groups (Kruskal-Wallis test) with a Bonferroni correction for multiple comparisons. To evaluate how well we sampled the true pool of var DBLα diversity (i.e. the true number of genetic variants circulating) in Ecuador and in South America, we used the R package vegan⁹¹ to generate species accumulation curves.

Data availability

The cleaned DBLα sequences generated in this study from Ecuadorian P. falciparum isolates have been submitted to the DDBJ/ENA/GenBank (BioProject Number: PRJNA642683). The python code for the DBLαCleaner pipeline is available on GitHub at https://github.com/Unimelb-Day-Lab/DBLaCleaner. The python code for the clusterDBLalpha pipeline is available on GitHub at https://github.com/Unimelb-Day-Lab/DBLalpha. The python code for the classifyDBLalpha pipeline is available on GitHub at https://github.com/Unimelb-Day-Lab/classifyDBLalpha. All other deidentified data and analysis code are available on the open-source GitHub repository at https://github.com/shaziaruybal/varcode-ecuador.