In this molecular epidemiological study, we examined P. falciparum isolates collected from 2013-2015 in Ecuador from individuals of all ages presenting with uncomplicated malaria cases confirmed by microscopy and/or PET-PCR (34). Details on the study population and data collection have been previously published (19, 35, 36). Briefly, these samples were collected during passive surveillance by the Ecuadorian Ministry of Health from consenting participants who lived in the areas where the samples were taken, were over 2 years old, agreed to participate and provided a blood sample (either venous blood or dried blood spot) and answered a basic demographic questionnaire (including places recently travelled and their address). Genomic DNA was extracted from venous blood or dried blood spots using a QIAamp DNA MINI KIT (QIAGEN, USA) as recommended by the manufacturer. The study was approved by the ethics committees at the Pontificia Universidad Católica del Ecuador (Quito, Ecuador) and The University of Melbourne (Melbourne, Australia).

Var genotyping PCR, sequencing details, and related data processing steps have been previously published (33, 37). Briefly, for each P. falciparum isolate, a single-step PCR to amplify the DBLα domains of var genes was performed using universal degenerate primers targeting homology blocks D (forward primer) and H (reverse primer) originally described in (38), with the addition of 10bp GS FLX Titanium multiplex identifier primer sequence for barcoding (39). The PCR reaction was prepared in a total volume of 40μl, containing 3μL of genomic DNA, MgCl₂ at a final concentration of 2mM, dNTPs at a final concentration of 0.07mM, each primer at a final concentration of 0.375mM and 3 units of Flexi DNA Taq polymerase (Promega). The cycling conditions were as follows: initial denaturation step of 2 min at 95°C was followed by 30 cycles of: annealing for 40 seconds at 95°C, extension for 90 seconds at 49°C, denaturation for 90 seconds at 65°C, and then a final extension step of 10 min at 65°C. The resulting individually barcoded DBLα amplicons (approximately 450-700bp in length) were pooled equimolarly and barcoded libraries were prepared using the KAPA Low-Throughput Library Preparation Kit (Kapa Biosystems). The libraries were sequenced on a MiSeq Illumina platform using the 2x300bp paired-end protocol and MiSeq Reagent kit v3 chemistry (Australian Genome Research Facility, Melbourne, Australia).

The raw illumina sequence data was then cleaned and processed using the DBLaCleaner pipeline [(37, 40), http://github.com/Unimelb-Day-Lab/DBLaCleaner]. Our customized bioinformatic pipeline has been described in detail in (37, 40). Briefly, we de-multiplexed and merged the paired reads as well as removed low-quality reads and chimeras using several filtering parameters (see Supplementary Figure S1 for details). This pipeline resulted in 2,141 cleaned DBLα sequences ( Supplementary Figure 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 (30) 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 (32) using the clusterDBLalpha pipeline (http://github.com/Unimelb-Day-Lab/clusterDBLalpha). 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). All the cleaned DBLα sequences in this study ( Supplementary Figure S1 ) have been submitted to the DDBJ/ENA/GenBank (BioProject Number: PRJNA642683). A tutorial on the DBLα processing pipelines can be found at http://github.com/Unimelb-Day-Lab/tutorialDBLalpha.

In addition, any P. falciparum isolate with low sequencing quality (< 10 DBLα types) was removed from the analysis. Therefore, from a total of 70 genotyped P. falciparum Ecuadorian isolates, 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 (30)] 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. 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 (41) to generate species accumulation curves where the number of unique DBLα types are plotted as a function of the number of sampled sequences. Plateauing of this curve indicates saturation and robust sampling of the diversity pool (i.e. sampling more sequences will not identify any new types).

Every parasite isolate was compared to every other parasite isolate in the population to determine the proportion of shared DBLα types. Theoretically, pairwise comparisons resulting in P_TS of 0 (0% shared types), 0.5 (50% shared types), and 1 (100% shared types) will reflect genetically distinct isolates with different var repertoires, isolates with recombinant or highly related var repertoires, and clones or genetically identical isolates with the same var repertoire, respectively. In practice, however, in low-transmission areas often “fixed” relatedness thresholds may obscure true relatedness estimates because of inbreeding events and thus require confirmation by measuring the uncertainty around each estimate with methods like our Bayesian inference of P_TS. We applied this approach to define “varcodes” as groups of isolates sharing ≥90% of their DBLα types (P_TS ≥ 0.90), identifying putatively identical genomes within the margin of error of detection of 1-5 DBLα types in an isolate. We confirmed our interpretations of varcodes, recombinants/highly related, and genetically distinct isolates at the thresholds of 0, 0.5 and 0.90-1, respectively, by comparing to unbiased BP_TS estimates (posterior means) and examining the HDPI.

To visualize the varcode relatedness between isolates as determined by P_TS or BP_TS, we constructed networks using the R packages ggraph (45) and tidygraph (46) where isolates are depicted as nodes and edges as the P_TS or BP_TS values at a given threshold. The R package ggspatial (47) was used to plot spatial networks using latitude/longitude coordinates for sampling location. To visualize varcode relatedness of parasites over time we used the R package gganimate (48) to construct spatiotemporal relatedness networks. We generated a clustered heatmap based on the presence/absence matrix of DBLα types to visualize the genetic profiles of each isolate in Ecuador and each country in South America using the R package pheatmap (49) and the “complete” clustering method. Unrooted neighbor-joining phylogenetic trees based on pairwise genetic distance (calculated as 1-P_TS) were constructed using the R packages ape (50) and ggtree (51, 52).

We used R version 3.5.2 (53), base R, and the R packages dplyr (54), epiR (55) for all analyses. We used 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.

Diversity of var genes was assessed by varcoding for 58 P. falciparum isolates that were collected between 2013 and 2015 from individuals of all ages experiencing clinical malaria ( Figure 2 ). These isolates represent 21% of the total cases reported in 2013, 61% in 2014, and 3% in 2015 (60% of the cases reported in January, 7% in May and 8% in November 2015). Overall, we identified 195 unique DBLα variants or types from 2,141 DBLα sequences ( Supplementary Figure S1 ) from the 58 isolates, representative of the diversity circulating in Ecuadorian P. falciparum populations, as indicated by sampling accumulation curves approaching a plateau ( Supplementary Figure S2 ).

To define distinct varcodes circulating locally in Ecuador, we estimated relatedness between two isolates’ varcodes by calculating the similarity index, Pairwise Type Sharing (P_TS) and constructed relatedness networks to identify genomes with the same varcodes (P_TS ≥0.90). We confirmed this by comparing varcodes derived from P_TS to those derived from posterior mean BP_TS estimates ( Supplementary Figure S3 ), which include statistical uncertainty, 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 36 isolates having varcode1 (representing the clonal outbreak in Esmeraldas City, salmon pink varcode Figures 3 and Figure 4A ) and at the other extreme varcode4, varcode5, varcode9 were seen only once ( Figure 4A , S4A ). Our definition of varcodes proved to predict genomic lineages (identity-by-descent ≥0.99) based on published whole genome sequence data in the case of 30 of the isolates (six varcodes) that were analyzed by both WGS (56) and varcoding. As expected, the outbreak varcode1 was clonal (indicated by HDPIs that included 1; Figure S5B ), as previously demonstrated by microsatellite genotyping (19) and more recently by WGS (56). 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 ( Figure 4B ). Overall, persistent disease transmission of the same varcodes was observed both over time ( Figure 4B ) and large distances ( Figures 3A–C ). 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. This is in line with reports in Peru (57) and Colombia (18, 58) where clonal parasites were shown to be circulating up to five and eight years after first identification, respectively.

We next constructed spatiotemporal varcode relatedness networks to explore signatures of local disease transmission after the outbreak. Figure 4C shows that 15 isolates with four varcodes (varcode3, varcode4, varcode6, varcode7) clustered with the outbreak varcode1 at the threshold of P_TS≥0.50 ( Figure 4C ). The remaining 7 isolates with different varcodes did not cluster in the network and were genetically distinct. These patterns were confirmed using BP_TS estimates based on posterior means ( Supplementary Figure S6A ) and their corresponding 95% HDPIs ( Supplementary Figures S5C–F , S6B, C ).

To better understand the varcode relatedness signatures in the “outbreak cluster” we examined both P_TS and BP_TS estimates. In some instances, highly related varcodes shared ~50% of the outbreak DBLα types, such as the outbreak varcode1 with varcodes3,4,6,7 (median P_TS range=37-58%, median BP_TS=45-70%, Figures S5C, D ). This level of varcode relatedness would be consistent with the generation of a new varcode through outcrossing by conventional meiosis resulting in recombinant var repertoires of the outbreak clone. When examining the pairwise combinations of varcodes 3,4,6,7, sharing of types also indicates high relatedness (median P_TS range=27-50%, median BP_TS=32-64%, Figures S5C, D ), pointing to the possibility of a prior cross followed by a backcross given the overall limited var diversity in the parasite population. However, the timing of such crosses is unknown. Interestingly, all the possibly recombinant and highly related varcodes (varcodes3,4,6,7) were identified in San Lorenzo with varcode7 being the main one, indicating this area is a transmission hot spot and a reservoir of parasites with diverse var repertoires ( Figure 3D ).

To further confirm the observed varcode relatedness signatures, we constructed a clustered heatmap to visualize the genetic profiles of each isolate based on the presence/absence of the 195 DBLα types, such that genetically similar isolates cluster together ( Figure 4D ). This analysis confirmed that in the case of isolates identified as highly related to the outbreak varcode1, a proportion of outbreak DBLα types as well as different DBLα types were present. By contrast, the genetic profiles of isolates with varcodes 2, 5, 8, and 9 had more different DBLα types to all other varcodes, i.e., parasites with genetically distinct varcodes, with sharing of only <30% of types in most instances (median P_TS range= 2-29% except for 2 pairwise comparisons 34-37%; median BP_TS range=4-30% except for 5 pairwise comparisons 31-41%, Figures S5E, F ).

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 post-outbreak, 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 highly-related/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 the cases sampled were caused by either parasites with the outbreak varcode1 (21%) or with any of the four highly-related/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 highly-related/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 sampled after the outbreak were caused by either parasites with the outbreak varcode1 (24%) or with any of the four highly-related/recombinants of varcode1 (56%), especially with varcode7 (representing 57% of the cases caused by parasites with highly-related/recombinant varcodes). Thus, disease transmission after the outbreak in 2014 and into 2015 was predominantly associated with parasites with highly related/recombinant var repertoires of the outbreak clone, with most of these cases occurring in the San Lorenzo hotspot.

We next examined the possible origins of the varcodes circulating locally in Ecuador by comparing to the only published var DBLα dataset from the region (30), comprising 128 P. falciparum isolates collected from South American P. falciparum populations in 2002-2008 from countries with higher malaria transmission compared to Ecuador ( Figure 5A ). In all countries (except some isolates in French Guiana), the number of DBLα types per isolate repertoire, a proxy for complexity of infection (COI, i.e. the number of P. falciparum genomes in an infection), was indicative of only one P. falciparum genome infecting an individual (i.e., repertoire size ≤ 50 or COI=1, Supplementary Figure S9 ; Supplementary Table S1 ). Despite the difference in sampling years, the 543 unique DBLα types identified in these 186 isolates represent the diversity circulating in South American P. falciparum populations, as indicated by DBLα sampling accumulation curves approaching saturation ( Supplementary Figure S7 ). An unrooted phylogenetic neighbor-joining tree showing pairwise genetic distance (1-P_TS) revealed distinct, country-specific clusters of isolates with related varcodes ( Figure 5B ), consistent with previous analyses demonstrating geographic population structure (30). Thus, despite the relatively low var diversity, there is sufficient resolution to differentiate between the parasite populations from different South American countries and to explore the possible geographic origins at the level of country-specific DBLα types but also the unique combinations of these types as varcodes.

Varcoding resolved a total of 97 varcodes in South America (P_TS ≥ 0.90), ranging from nine in Ecuador and Venezuela, to 56 varcodes in French Guiana and the varcodes were country-specific since none of the isolates from different countries clustered with isolates from another country at this threshold ( Supplementary Figure S8 ). However, although identical varcodes were not identified in multiple countries, we were interested in examining varcode relatedness. To assess the overall sharing of types in the varcodes identified in Ecuador to those sampled in other South American parasite populations, we aggregated all possible DBLα types identified in the isolates comprising each Ecuadorian varcode (range: 34-47 types per varcode) and all DBLα types seen in any isolate from a given country (range: 112-249 types per country). The overall sharing of types was high, with the highest median sharing with Peru (53%, range: 26-72%) followed by Colombia (39%, range: 19-89%) and Venezuela (27%, range: 12-76%), and only 10% with French Guiana (range: 3-86%) except for varcode5 (86%) ( Figure 5D ).

Next we estimated P_TS between all possible isolate pairs to examine whether any of the Ecuadorian parasites had varcodes that were genetically related to these historical South American parasites (P_TS ≥ 0.50) and constructed regional varcode relatedness networks ( Figure 5C ). We identified a genetically related Peruvian P. falciparum isolate that clustered with the “outbreak cluster”, but especially with outbreak varcode1 (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 (13, 19). More recent whole genome sequence data also points to the outbreak varcode1 circulating in Colombia in the early 2000s (56). No other South American isolate clustered with the “outbreak cluster” based on our network analysis, suggesting that there may be other locally circulating parasites in unidentified reservoirs in Ecuador, e.g. the parental varcodes of the possible recombinants. In the case of varcode3, the WGS study showed that isolates with varcode1 and varcode3 belong to distinct genomic lineages but are identical by descent across 80% of the genome thus belong to a “super-cluster” (56), which is consistent with our BP_TS estimates of 70-84% type sharing ( Figure S5D ). Additionally, the WGS study identified varcode6 and varcode7 in Colombian isolates from 2014-2016 and in 2016, respectively, but varcode4 was not identified in any of their Colombian isolates (56).

We were also interested to better understand the origins of those varcodes that did not cluster with the “outbreak cluster” since historically they may have been imported and may represent other locally circulating parasites that were not a direct result of onward transmission after the outbreak. We found evidence of two putative importations of parasites from neighboring countries, varcode9 from Colombia (Case EC53, Table 1 ), based on both genetic and epidemiologic data, and varcode5 (Case EC40, Table 1 ) the origin of which was less clear. The epidemiological data for the putative infection location was recorded as Peru (Case EC40, Table 1 ), however based on our data this varcode is more related to historical P. falciparum populations from French Guiana and Venezuela than Peru. The remaining Ecuadorian isolates with varcodes 2 and 8 did not cluster with any other South American isolate. The epidemiological data for the putative infection location for varcode8 was recorded as Colombia (EC49 and EC50, Table 1 ), but these isolates did not cluster with Colombian isolates in our network analysis. Sharing of DBLα types in these varcodes with Colombian types was high (44% for varcode2 and 67% for varcode8, Figure 5D ), providing evidence that parasites with varcode2 and varcode8 may represent residual parasites historically imported from Colombia. This is consistent with WGS data for isolates with varcode2, demonstrating they are identical by descent (IBD>0.99) with Colombian isolates from the early 2000s (56), suggesting past rather than recent importation. It is worth noting that varcode2 was identified in the hotspot of San Lorenzo, Esmeraldas confirming that this location may act as a reservoir of parasites with diverse var repertoires, as described above.