Dynamic changes in genetic diversity, drug resistance mutations, and treatment outcomes of falciparum malaria from the low-transmission to the pre-elimination phase on the isolated islands of São Tomé and Príncipe

With effective vector control and case management, substantial progress has been made in the elimination of malaria on the islands of São Tomé and Príncipe (STP). During the critical period from the low-transmission to the pre-elimination phase, this study tracked the dynamic changes in the genetic diversity in Plasmodium falciparum, the distribution of antimalarial drug-resistance genes, and the treatment outcomes in patients to provide insights for the prevention of rebounded malaria in STP. Dried blood spots (DBSs) and case follow-up data were collected from malaria patients who had visited the Central Hospital between 2010 and 2016. Genomic DNA of P. falciparum was extracted from DBSs. The polymorphic regions on the genes for merozoite surface proteins 1 and 2 (msp1 and msp2) were amplied in 118 pre-treatment samples to identify the genetic diversity of the infected parasites. Antimalarial drug resistance mutations in the multi-drug resistance (pfmdr1), chloroquine resistance transporter (pfcrt), and kelch 13 (pfK13) genes were genotyped using polymerase chain reaction-restriction fragment length polymorphism (PCR-RFLP) and DNA sequencing in 111 samples. Treatment outcomes were categorized based on the parasitological results from microscopy during the 28-day follow-up after treatment. Factors related to malaria recurrence were characterized by logistic regression models using case follow-up data (total number = 7,482).

Recently, biological challenges associated with malaria have focused on the species composition of vectors and parasites (indigenous or imported), insecticide resistance in mosquitoes, and drug resistance in parasites [13][14][15]. Our previous study reported that there is only one vector species, Anopheles coluzzii, that is responsible for malaria transmission in STP, and that it has developed insecticide resistance mutation against pyrethroid after long-term application of IRS [7]. With respect to malaria parasites, P. falciparum is the major cause of malaria infections in STP [10]. The genetic structure of P. falciparum was studied in the parasite population collected at the pre-intervention stage, before 2004, by analyzing the diversity levels of microsatellite loci [16]. The authors detected differences in parasite populations across the years of 1997, 2000, and 2004, showing that local malaria control strategies could cause dynamic changes in parasite populations [16]. Over a decade later, with the expansion and transformation of the national malaria control program, changes in local parasite populations were expected but have not yet been proven. To examine the post-intervention changes in P. falciparum populations in STP, this study assessed the polymorphic genetic markers merozoite surface protein 1 and merozoite surface protein 2 (msp1 and msp2) with the aim of characterizing the genetic diversity of the parasite from the low-transmission (2010-2013) to the pre-elimination period (2014)(2015)(2016). MSP1 and MSP2 are asexual blood-stage antigens that are considered vaccine candidates, and are common markers for the identi cation of parasite populations due to their high allelic diversity [17,18]. MSP1 has three allelic families, namely, K1, MAD20, and RO33, based on the variable sequence in the block 2 region [19]. MSP2 has two allelic families, namely, FC27 and 3D7/IC, based on the central repetitive domain in the block 3 region [20]. The patterns and extent of polymorphisms on msp1 and msp2 can indicate the diversity of local malaria infections and the transmission dynamics [21], providing a better understanding of post-intervention effects on P. falciparum in the STP.
Several drug resistance surveys were conducted before 2004 in STP, of which mutations in the resistance markers including chloroquine transporter (pfcrt), multidrug resistance I (pfmdr1), dihydrofolate reductase (pfdhfr), and dihydropteroate synthase (pfdhps) were found when chloroquine and sulfadoxinepyrimethamine (SP) were used as standard treatments [22,23]. Subsequently, there were no updated reports for over a decade after the introduction of ACT as a standard treatment. ACT is a combination of a short-acting drug, an artemisinin derivative, with another long-acting partner drug, commonly amodiaquine (AQ), lumefantrine (LF), me oquine (MQ), piperaquine (PPQ), or sulfadoxine-pyrimethamine (SP) [24]. Artemisinin resistance is associated with mutations in the propeller domain of kelch 13 (pfK13) [25], while resistance to ACT partner drugs, including AQ and LF, used in STP is associated with mutations in pfmdr1 and pfcrt [26,27]. Therefore, surveillance of these drug-resistance mutations is crucial to understand the current treatment e cacy of antimalarials in STP. This study targeted malaria patients showing recurrent parasitemia after a 3-day regimen. Epidemiological analysis and genotyping of antimalarial drug-resistance markers in recurrent patients were performed to characterize the potential risk factors that may induce treatment failure. In brief, temporal changes in P. falciparum from the lowtransmission to the pre-elimination phase were of interest for this study, which aimed to address the challenges posed by the genetic patterns and antimalarial drug resistance of P. falciparum on the endemic foci of STP.

Ethics statement
The study content and methods were reviewed and approved by the Centro Nacional de Endemias of the Ministry of Health in the STP (o cial No. 19/P0CNE/2016) and the Research Ethics Committee of National Taiwan University Hospital (NTUHREC No. 201110023RD).

Study site
The island nation STP is located in the Gulf of Guinea, 400 km off the coast of Central Africa [10]. The total area is 1,001 km 2 , and the country has nearly 200,000 residents [7]. It consists mainly of two islands. Príncipe offshore island is an autonomous region with a low resident population; while the São Tomé main island consists of six administrative districts, namely, Agua Grande (AG), Me-Zochi (MZ), Lobata (LO), Cantagalo (CT), Lemba (LE), and Caue (CU). Due to the distribution of volcanic mountains in the southwest of São Tomé from LE to CU, the population is mainly concentrated in the districts of AG, LO, and MZ, which are plains along the east coast [7]. The climate is equatorial, with an average temperature of 25-27℃, and with frequent rain except in the dry season from June to early September [7].
Genomic DNA extraction from dried blood spots DNA extracts of parasites from infected hosts were obtained from dried blood spots (DBSs), which were collected by nger-pricking patients and spotting the blood on Whatman ® 903 lter cards. This study evaluated two types of samples for analysis: (a) a total of 118 pre-treatment samples (day 0) collected from 2010 to 2016 and (b) a total of 60 pre-treatment (day 0) and 51 post-treatment (day R, R = 3-28) samples showing positive microscopic results in follow-up assessments from 2014 to 2016. The former sample set was used to identify the genetic diversity of infected parasites through analyses of the msp1 and msp2 polymorphic regions. The latter sample set was chosen from 41 recurrent and 10 non-recurrent patients to compare the genotyping results of antimalarial drug resistance markers. All samples were collected from patients with malaria who had visited the Central Hospital Ayres de Menezes (HAM) located in the capital district, AG, from 2010 to 2016. Specimen collection was approved by the Centro Nacional de Endemias of the Ministry of Health in STP (o cial No. 19/P0CNE/2016) and the Research Ethics Committee of National Taiwan University Hospital (NTUHREC No. 201110023RD). Each DBS was air-dried and sealed in a plastic bag labeled with the sampling ID and collection date. Personal information was recorded in encrypted case follow-up data.
Genomic DNA was extracted from three 6-mm circular punches on DBS samples using the Geneaid Genomic DNA Mini Kit for Tissue (Geneaid Biotech Ltd., Taipei, Taiwan). The blood discs were immersed in 1000 μL of lysis buffer and protease K followed by shaking and incubation at 60 ℃ for 1.5 h. The liquid mixture was then transferred to a DNA-binding membrane and washed according to the manufacturer's instructions. Genomic DNA was eluted using 80 μL of the preheated elution buffer. The extracted DNA was examined using the Loopamp™ Malaria Pan/Pf Detection kit (Eiken Chemical, Tokyo, Japan) to con rm the presence of P. falciparum DNA. DNA was 10-fold diluted to a volume of 30 μL for each loop-mediated isothermal ampli cation (LAMP) reaction [28]. Positive results were determined based on the presence of uorescence under UV light (365 nm). Samples con rmed to contain DNA from P. falciparum were stored at -20 ℃ for further analysis.
PCR and sequencing of msp1 and msp2 Block 2 of msp1 and block 3 of msp2 containing variable repeats were anked by primer pairs in the conserved regions [29]. The sequences were obtained from 118 pre-treatment samples by PCR with a reaction mixture containing 10 μL of 2X HotStarTaq Master Mix (Qiagen, Hilden, Germany), 10 μM of forward and reverse primers, 3 μL of the DNA template, and RNase-free water to a total volume of 20 μL. The cycling conditions for PCR included an initial denaturation at 95°C for 15 min, followed by 40 cycles of denaturation at 95°C for 30 s, annealing at 50°C for 30 s, extension at 72°C for 30 s, and a nal extension at 72°C for 10 min using the Biometra TRIO Thermal Cycler (Analytik Jena AG, Jena, Germany). The fragment lengths of the ampli ed products varied from 200 to 400 bp in msp1 and 400 to 600 bp in msp2 based on the tandem repeats in the sequences. PCR products showing a single band on 2.5% agarose gel were sent to Mission Biotech Co., Ltd., Taipei for direct sequencing. PCR products showing multiple bands on the gel were suspected to indicate polyclonal infections. The bands of different sizes were gel-puri ed and sequenced separately.
Thirty samples were further tested by nested PCR (nPCR) using family-speci c primers for two purposes [30,31]. One was to identify polyclonal infections, and the other was to discriminate new infections (due to new infectious mosquito bites) from recrudescence (incomplete clearance of asexual parasites after antimalarial treatment) in post-treatment samples [18,32]. Nested PCR was performed according to the methods described in the published protocols [30,31]. The nal nPCR products were double-checked by both traditional gel electrophoresis and capillary electrophoresis using the QIAxcel Advanced System (Qiagen, Hilden, Germany). Discrimination between new infection and recrudescence was performed in accordance with the WHO protocol [18]. Recrudescence was determined by the presence of at least one allele at each locus that was identical between the paired pre-and post-treatment samples, while new infection was determined when all alleles identi ed from post-treatment samples were different from those in the pre-treatment samples.
Genotyping of pfmdr1, pfcrt, and pfK13 polymorphisms A total of 60 pre-treatment and 51 post-treatment samples were genotyped for resistance-associated single nucleotide polymorphisms (SNPs) on pfmdr1, pfcrt, and pfK13. Detection of SNPs on pfmdr1 N86Y, Y184F and pfK13 was performed by PCR sequencing using primers designed in previous studies [33,34]. The PCR reaction mixture consisted of 10 μL of 2X HotStarTaq Master Mix, 10 μM of forward and reverse primers, 3 μL of DNA template, and RNase-free water to a total volume of 20 μL. The PCR program was as follows: initial denaturation at 95°C for 15 min, followed by 40 cycles of denaturation at 95°C for 30 s, annealing at 50°C for pfmdr1 and 54°C for pfK13, extension at 72°C for 30 s, and a nal extension at 72°C for 10 min. The amplicon sizes for pfmdr1 86+184 and pfK13 were 353 bp and 986 bp, respectively. PCR products of correct sizes were sequenced using an Applied Biosystems 3730xl DNA Analyzer (Thermo Fisher Scienti c, Waltham, MA, USA).
The pfmdr1 D1246Y and pfcrt K76T were genotyped by PCR-RFLP using methods described elsewhere [35,36]. Brie y, the PCR product of pfmdr1 1246 was digested with EcoRV (20,000 units/mL, New England Biolabs, Massachusetts, USA), and the fragment of pfcrt 76 was digested with ApoI endonuclease (10,000 units/mL, New England Biolabs, Massachusetts, USA). After the RFLP reaction, DNA bands were visualized on 2.5% agarose gels. In addition, we selected 25 samples of pfcrt 76 fragments (145 bp) to sequence and identify the haplotype of pfcrt 72-76. Genomic DNA from the P. falciparum reference clone 3D7 (Catalog Number PRA-405D, ATCC, Virginia, USA) was used as the wildtype control for the genotyping results of pfmdr1, pfcrt, and pfK13.
Sequence alignment and data analysis All sequences were aligned using Lasergene SeqMan version 7.1 and BioEdit version 7.0. The MEGA 7 tool was used to transform DNA sequences into amino acid codes, and to construct phylogenetic trees for msp1 and msp2. The protein sequence alignment of MSP1 and MSP2 was based on the sequence structures illustrated by Kang et al. [37], including family-speci c regions and tandem repeat regions. MSP1 haplotypes were categorized into three families: K1, MAD20, and RO33. MSP2 haplotypes were categorized into two families, FC27 and 3D7/IC, based on family-speci c structures.
Maximum likelihood (ML) trees for msp1 and msp2 were constructed by bootstrapping 100 times using the best substitutional model computed by the general time reversible (GTR) model. Chi-square and Fisher's exact tests were performed to test the associations of the msp1 and msp2 allelic families with antimalarial drug-resistance types and the characteristics of the human hosts.
In addition, SNPs in antimalarial drug-resistance markers were compared in 41 patients showing recurrence. Matched-pair comparisons of these markers between pre-treatment (day 0) and posttreatment (day R) samples were carried out using McNemar's test. Mixed alleles consisting of both wild and mutant types were considered as mutants in the McNemar's test.
Characterization of infective parasite density, treatment outcomes, and factors related to recurrence The malaria epidemic curve in STP ( Figure 1) shows four periods corresponding to the control, lowtransmission (ascending), low-transmission (descending), and pre-elimination phases from 2003 to 2016. Our sampling period was from 2010 to 2016 (periods B to D), which was the critical duration of the transition from the low-transmission to the pre-elimination phase in STP [7]. This study aimed to explore the temporal changes in parasites infecting hosts across these periods. The case follow-up data documented patients' ID, name, sex, age, residency (district and village), date of con rmed diagnosis, microscopic results of the initial infection, treatment, and microscopic results during follow-up (days 3, 7, 14, 21, and 28 post-treatment). Malaria cases that had been registered in the Central Hospital Ayres de Menezes (HAM) from 2010 to 2016 with complete records (excluding cases with incomplete baseline information and those lost to follow-up) were enrolled in the analysis to explore factors related to infective parasitemia and post-treatment outcomes, particularly recurrence during follow-up. First, the infective parasite density at the entry point was recorded as "grade (F1-F10)/ parasite counts per μL blood" (Table S1) [38]. Statistical analyses used log 10 -transformed parasite counts as a parameter indicating the infective parasite density. Factors related to infective parasite density at baseline were evaluated using a linear regression model. Second, treatment outcomes were classi ed based on the parasitological results irrespective of axillary temperature, which was not recorded in the case follow-up data. The de nition of each treatment outcome was based on the WHO guidelines [18,39]. All patients with complete records were classi ed into one of the following groups: adequate clinical and parasitological response (ACPR), early treatment failure (ETF), insu cient clearance on day 3 (ISC), late parasitological failure (LPF), and gametocyte carriage (GC), as shown in Figure 2. In this classi cation scheme, ISC and GC were the two additional outcomes that were not speci ed in the WHO guidelines [39], but were separately categorized in our study.
This was because, in accordance with the strict de nition of recurrent malaria, the presence of gametocytes (sexual parasitemia) only is not considered as recurrence [32]. Therefore, patients showing only gametocytes (sexual parasitemia) after treatment were classi ed into the GC group. Early treatment failure (ETF) was de ned as the presence of danger signs or higher parasitemia on day 3, i.e., ≥25% of the parasite counts on day 0. Those who showed mild symptoms and decreased parasite counts on day 3 after treatment were classi ed into the ISC group in our study.
Finally, factors related to recurrence were evaluated using logistic regression models. Two response variables were considered: "early (day 3 positive)" and "late (day 7+ positive)" parasitological failures after treatment. The host-related predictive variables were patient age, gender (male vs. female, female as reference), and residency (urban vs. other districts, the AG, MZ, and LO district, which included 75% of the citizens, were de ned as urban). Time-related predictive variables were the time point of initial infection (year was transformed to "period" as de ned in Figure 1; month was transformed to "endemic season" as de ned in Figure S1). Diagnostic and treatment variables were log 10 parasite density at enrollment, and treatment selection (quinine vs. ACT, with ACT as reference). Mediation analysis was carried out for log 10 parasite density and treatment selection for early and late recurrence. All statistical analyses were performed using R version 4.0.2.

Results
Sequence polymorphisms ofmsp1block 2 inP. falciparumisolates from STP Twenty-two haplotypes of msp1 block 2 were identi ed in 118 samples collected between 2010 and 2016. The sequence accession nos. for the msp1 haplotypes were MW001371-MW001392. After translating nucleotides into protein sequences, nine haplotypes belonged to the K1 family, nine belonged to the MAD20 family, and four belonged to the RO33 family. The K1 and MAD20 allelic types contained tri-peptide repeat units, showing length polymorphisms from 312 to 366 bp (104-122 amino acids) in the K1 types, and 285-411 bp (95-137 amino acids) in the MAD20 types. In the K1 haplotypes, anked by family-speci c regions (E1 and E2 in Figure S2), the tri-peptide repeat region (R1) started from SAQ and terminated with SGT, and could be divided into two subgroups on the phylogenetic tree according to the central repeat region (Fig. 3). One group included KH2, 4, 5, 6, and 9, which contained SAQ and SGA-rich repeats, and was more similar to the P. falciparum NF7 clone isolated from West Africa. Another group included KH1, 3, 7, and 8 consisting of SGT and SGP-rich repeats, and was more similar to the 3D7 reference clone. Among the total 118 MSP1 sequences, K1 haplotypes were identi ed in 41 isolates (34.7%, 41/118). The prevalent haplotype in the K1 allelic family was KH1 (36.6%, 15/41), followed by KH3 (26.8%, 11/41).
In the MAD20 allelic types, the tri-peptide repeat region (R1) started with SKG, SSG, or SGG motifs, followed by the SVT motif ( Figure S2). The central repeat region majorly contained SGG and SVA-rich repeats and terminated with SGG. The predominant type in the MAD20 family was MH3 (63.5%, 33/52), which was translated from two nucleotide sequences with a synonymous SNP, represented as MH3-1 and MH3-2 on the phylogenetic tree (Fig. 3). Except for MH3, each MSP1 haplotype was translated from a single nucleotide sequence. Two subgroups were found in the MAD20 family in the phylogenetic tree. One included MH1, 5, 6, 7, and 8, and was grouped together with the HB3 reference sequence, and sequences identi ed in Tanzania and Oceania islands. Another subgroup included MH2, 3, and 4, which was similar to sequences identi ed in Southeast Asia and Ghana.
The RO33 allelic type had the least isolates (21.2%, 25/118) and haplotype diversity in STP across the sampling period. The RO33 allele did not encode any tandem repeats. The RO33 haplotypes had identical lengths (285 bp) with four nonsynonymous SNPs ( Figure S2). The predominant haplotype in the RO33 family was RH1 (84%, 21/25), isolated from samples collected after 2012 across all the districts in São Tomé. Phylogenetic analysis showed that RH1 and RH2 were grouped together with the RO33 clone isolated from Ghana, while RH3 and RH4 were grouped together with the sequence isolated from Malawi. Sequence polymorphisms ofmsp2block 3 inP. falciparumisolates from STP Twenty-one haplotypes were identi ed from 118 samples, of which 11 belonged to the 3D7/IC family (IH1-IH11, accession nos. MW001393-MW001403), and 10 belonged to the FC27 family (FH1-FH10, accession nos. MW001404-MW001413). Two tandem repeat regions of MSP2 block 3 (R1 and R2) were anked by three family-speci c regions (E1-E3, Figure S3). The R1 region in 3D7/IC haplotypes consisted of GSA-rich repeats differing in the lengths and numbers of repeat units, resulting in diverse length polymorphisms (from 430 to 600 bp). The R2 region of 3D7/IC was composed of 5-14 poly-threonine (T) residues, and only IH2 and IH9 presented additional copies of the TPA motif before the poly-threonine segments. Several nonsynonymous amino acid substitutions and deletions were found in the familyspeci c regions (E1-E3) of the 3D7/IC types. At the bottom of the phylogenetic tree of 3D7/IC (Fig. 4A), IH5 and IH11 were grouped together with reference sequences from Asian countries. The middle part of the tree showed that IH1, 2, 8, 9, and 10 were similar to sequences identi ed in Gambia and Nigeria.
Finally, IH3, 4, 6, and 7 were grouped together with the 3D7 clone and sequences identi ed from Asian and Oceanian countries.
Among the MSP2 sequences, the majority (72.9%, 86/118) belonged to the FC27-allelic type. Ten haplotypes were identi ed in the FC27 family, and FH6 and FH1 were the most prevalent types, accounting for 30% each. The FC27 haplotypes were divided into two subgroups. One subgroup included FH1, 2, and 3, containing three copies of the family-speci c repeats (32 amino acids) in the R1 region, without showing tandem repeats in the R2 region. In contrast, another subgroup included FH4 to FH10, showing tandem repeats of 12 amino acids (aa) in the R2 region, without presenting any repeats in the R1 region ( Figure S3). At the top of the phylogenetic tree (Fig. 4B), the prevalent haplotype FH6 was most similar to the reference sequences from the neighboring country, Gabon. On the other hand, haplotypes FH1 to FH3 located at the bottom of the tree were grouped together with the sequences found in Asian countries, Myanmar and India.
Temporal changes ofmsp1andmsp2alleles before and after the recent peak incidence Both MSP1 and MSP2 sequences were acquired from 118 pre-treatment samples collected from 2010 to 2016, which was a period of low transmission leading to the pre-elimination phase. During this period, a slight increase in annual malaria incidence was found in 2012 (incidence rate = 6.9%, Fig. 1). The incidence declined after 2012 and decreased to ~ 1% after 2014, reaching the pre-elimination criteria. This study tracked the temporal changes in MSP1 and MSP2 in infective parasites before and after the recent incidence peak. The results showed signi cant changes in both MSP1 and MSP2 allelic types before and after the peak incidence at 2012 (MSP1: χ² = 7.94, df = 2, p = 0.02; MSP2: χ² = 23.46, df = 1, p < 0.001), with the predominant alleles in MSP1 changing from K1 to MAD20, and those in MSP2 changing from 3D7/IC to the FC27 type (Fig. 5).
The major MSP1 allelic types before 2012 were found to be the haplotypes KH1 and KH8, which were not detected in the samples isolated after 2012. An increase in MAD20 and RO33 allelic types was found in parasites collected after 2012, particularly in the MH3 and RH1 haplotypes (Fig. 5). Haplotype MH3 accounted for 12% of the parasites isolated before 2012, but became the predominant type after 2012, accounting for 27%-35% of the prevalence. Haplotype RH1 was not identi ed in the isolates collected before 2012, but was found to constitute 20%-27% of the parasite population isolated after 2012.
Before 2012, 67% (18/27) of the parasite isolates were identi ed as the MSP2-3D7/IC allelic types, of which IH1 and IH4 were the prevalent types, and the minor types were IH3 and IH5-IH10 (Fig. 5). However, after 2012, only IH1, IH2, IH3, and a new IH11 were found in the MSP2-3D7/IC allelic family, and the prevalence of 3D7/IC alleles dropped by 50%. The predominant MSP2 types were replaced by FC27 alleles, since their prevalence increased by 50% after 2012. Haplotypes FH1 and FH6 were present throughout the sampling period in our sample set, and increased their frequencies after 2012. Haplotype FH5 was isolated only after 2014, and became the predominant type in 2016.
In our sample set, 6.8% (8/118) of the samples showed polyclonal infections by nested PCR and gel puri cation (Table S2). Two patients were infected with more than one haplotype in both MSP1 and MSP2, while the other six patients were infected with 2-3 haplotypes in MSP1, but a single haplotype in MSP2. The majority of the samples were identi ed as monoclonal infections, being infected by one haplotype in either MSP1 or MSP2. Overall, the prevalent circulating parasite strains before 2012 were MSP1-K1 and MSP2-3D7/IC, and those after 2012 were MSP1-MAD20 and MSP2-FC27. The changes occurred in the year of peak incidence, showing that the prevalent strains had changed from the lowtransmission to the pre-elimination period in STP.
The post-treatment samples were selected from 41 patients showing inadequate treatment outcomes, of which nine (22.0%, 9/41) of them had insu cient clearance on day 3, nine (22.0%, 9/41) showed early treatment failure on day 3, 13 showed late parasitological failure after day 7 (31.7%, 13/41), and ten (24.4%, 10/41) showed gametocytes during follow-up. Thirty-one (75.6%, 31/41) of these patients were treated with quinine and 10 (24.4%, 10/41) were treated with ASAQ on day 0. Five patients were originally treated with quinine and showed positive parasitemia for more than one follow-up day, and the alternative drugs were administered afterwards: ASAQ was administered to three patients, AL was administered to one patient, and both ASAQ and AL were administered to one patient after day 7. Fisher's exact test showed no signi cant differences in the pfmdr1 and pfcrt polymorphisms after treatment with either quinine or ASAQ at baseline.
Matched-pair comparisons of pfmdr1 and pfcrt polymorphisms were performed between pre-treatment (day 0) and post-treatment samples (day R) by McNemar's test, and no signi cant changes were found. More than 80% (41/51) of the parasite isolates on day R shared identical pfmdr1 and pfcrt polymorphisms with their initial infections on day 0. This was also con rmed by the genotyping results of msp1 and msp2 to discriminate recrudescence from new infections. The results showed that 83% (34/41) of the recurrent patients shared the same msp1, msp2, pfmdr1, and pfcrt alleles between pre-and post-treatment samples, indicating that the positive parasitemia at follow-up was mainly due to the incomplete clearance of the initial infections (recrudescence). However, seven patients (17%, 7/41) showed different pfmdr1 or pfcrt patterns between day R and day 0 (Table S3). Among these seven patients, ve showed substitutions one week after treatment, while two showed early substitutions at day 3. The substitutions mostly changed from pfmdr1 86Y to N86, 184F to Y184, and pfcrt 76T to K76 after treatment with quinine (6 patients) and ASAQ (1 patient) at baseline. One patient who was treated with AL by day 7 showed substitutions from pfmdr1 86Y to N86 and pfcrt 76T to K76 on days 7, 9, and 14. Moreover, through con rmation by nested PCR in these seven patients, four were con rmed to have polyclonal infections at day 0, and new infections were also detected in four recurrent cases. Thus, substitutions of pfmdr1 or pfcrt types may be due to the polyclonal infections at the beginning or new infections during follow-ups.
Finally, temporal changes were only found in pfmdr1 Y184F among the 60 pre-treatment samples from 2014 to 2016 (χ² = 8.25, df = 2, p = 0.02). The pfmdr1 184F type showed a signi cant increase of 30%-40% in 2016 ( Figure S5). In conclusion, from 2014 to 2016, pfmdr1 86Y (75%-95%) and pfcrt 76T (95%-100%), which were associated with the reduced sensitivity of AQ, the rst-line ACT partner drug used in STP, were constantly present in the local parasite population. The pfcrt 76T mutation was nearly xed in the local parasite population, and the pfcrt 72-76 mutant type was identi ed as the CVIET haplotype from sequencing results. Among these markers, pfmdr1 Y184F showed the largest variation in our sample set from 2014 to 2016, while the other markers showed the same prevalent types across the period.

Characteristics of infective parasite density, treatment outcomes, and factors related with recurrence
To examine the host and temporal factors that may affect parasite density in initial infections, a linear regression model was used. The outcome variable was log 10 -transformed parasite density (parasites/µL) at the date of initial infection. The predictive host-related variables were age, gender, residential area, and time-related variables, including period and endemic season. Results showed that "age" is a signi cant factor (coe cient = -0.012, p < 0.001), which may decrease the infective parasite density with an increase in age (Table S4). The levels of parasitemia in infected hosts were signi cantly higher (coe cient = 0.40, p < 0.001) in the low-transmission period (2010-2013) than in the pre-elimination period (2014-2016).
A total of 7,482 malaria cases with complete records were enrolled in the analysis, and 86.7% (6,490/7,482) of the cases showed adequate clinical and parasitological responses (ACPR) after treatment. The ACPR rate in the ACT treatment group ranged from 88.2%-97.6%, which was slightly higher than that in the quinine treatment group (73.6%-94.0%). The proportion of cases with inadequate treatment outcomes (positive parasitemia at days 3-28 after treatment) after quinine treatment (17.3%, 847/4,906) was higher than that after ACT (5.6%, 145/2,576). In particular, patients treated with quinine were mostly hospitalized patients with higher parasitemia levels at enrollment, and their follow-up assessments showed higher probabilities of insu cient clearance on day 3 (0.7%-14.2% in patients receiving quinine and 0.3%-4.2% in those receiving ACT) and late parasitological failure after day 7 (3.6%-7.8% in patients receiving quinine and 0.8%-4.7% in those receiving ACT) (Fig. 6). In the lowtransmission period, more than half of the patients (61.5%-81.5%) were treated with quinine, while a lower proportion of patients were treated with ACT (18.5%-38.5%). After the decrease in parasitemia levels in the pre-elimination period, cases treated with ACT were greater or nearly equal to the quinine-treated numbers, and the overall ACPR rate was higher in the pre-elimination phase (2014-2016) than in the lowtransmission period (2010-2013).
Logistic regression models ( Table 2) and mediation analysis ( Figure S6) were performed to characterize the factors that were correlated with early and late parasitological failure. Results from logistic regression models showed that age was a signi cant protective factor against recurrence (OR = 0.97 for early recurrence, and OR = 0.98 for late recurrence, both p < 0.001). Thus, younger children had higher odds of recurrence than the elderly. In the low-transmission period, the odds of early and late recurrence were signi cantly higher than those in the pre-elimination period (OR = 1.63-5.82, p < 0.01). Both parasite density (log 10 transformed, OR = 1.44, p < 0.001) and treatment selection (quinine vs. ACT, OR = 1.96, p < 0.001) were signi cant risk factors for early recurrence; however, only treatment selection (quinine vs. ACT, OR = 1.91, p < 0.001) was a signi cant risk factor for late recurrence. Since treatment selection would be partially determined by the infective parasite density at enrollment, the mediation pathway of "initial parasite density→ treatment selection→ recurrence" was further tested in the model. The direct effect of parasite density on early and late recurrence was 0.46 (p < 0.001), and 0.16 (p = 0.009), respectively. After partial mediation by treatment selection, the coe cient of parasite density declined, and its indirect effect on both early and late recurrence was 0.09. Moreover, log 10 parasite density was no longer a signi cant factor for late parasitological failure after the mediation effect of treatment selection. In conclusion, quinine treatment was a signi cant risk factor correlated with early and late recurrence, and a partial mediator between infective parasite density and recurrence during follow-up.

Discussion
This study tracked the dynamic changes in genetic diversity, antimalarial drug resistance in local parasites, and treatment follow-ups in patients in STP. The major novel ndings were as follows: (a) signi cant temporal changes in circulating parasite strains in malaria cases following a slightly rebounded incidence (b) the most frequently observed pfmdr1 86 + 184 + 1246-pfcrt 76 polymorphisms in parasites isolated from STP were YFD-T during the low endemic periods (c) young children, patients with high parasitemia level at enrollment, and receiving quinine treatment were observed to be more likely to experience recurrent malaria during follow-ups. A thorough understanding of the molecular epidemiology of falciparum malaria in this study can provide useful information for strengthening national malaria control program in STP.
A previous study analyzed 180 P. falciparum isolates derived from STP in 2000, and showed that the frequency of the MSP1-K1, MAD20, and RO33 types was 50%, 44%, and 6%, respectively [17]. This was similar to the frequency detected before 2012 (K1 = 54%, MAD20 = 42%, RO33 = 4%). After the peak incidence in 2012, the frequencies of the MAD20 and RO33 allelic types increased, and even surpassed that of the original predominant type, the K1 family. Similar results were also observed for the MSP2 marker. The proportion of MSP2-3D7/IC and FC27 alleles detected from 2010 to 2011 was similar to that in 2000 [17], of which 3D7/IC was the prevalent allelic type (60%-65%) and FC27 was the minor type (35%-40%). However, after 2012, FC27 replaced 3D7/IC as the dominant type of STP. Since MSP1 and MSP2 are allele-speci c peptides that can trigger family-speci c immune responses in hosts, reduced immuno-protection to the prevalent types during the epidemics may lead to the detection of more symptomatic cases during or after the epidemic peak. In terms of epidemic expansion, the populations of certain clones expanded, which could be attributed to selective antibody or antimalarial drug pressures or a founder effect at the beginning of the epidemic [40][41][42]. Several studies also found temporal changes in the genetic diversity of P. falciparum due to transmission intensity. A reduction in the genetic diversity of P. falciparum from the Grande Comore Island was found 10 years after the introduction of ACT [43]. On the other hand, the distribution of msp1 and msp2 alleles in Djibouti, Northern East Africa, also showed changes in the prevalent type before, during, and after the 1999 outbreak, indicating the expansion of a few strains that had already existed locally during the epidemic [41]. Although this study was conducted with limited sample numbers, the haplotypes that showed expansion after 2012 in STP were MH3 and RH1 in MSP1, and FH5 and FH6 in MSP2. MSP1-MH3 and MSP2-FH6 were found before 2012 and increased in frequencies afterwards, whereas MSP1-RH1 and MSP2-FH5 were not detected in the samples collected before 2012. On the basis of these ndings, changes in the prevalent parasite strains may be due to the expansion of a few haplotypes that originally existed locally, or some new strains that may have evolved or imported into the hotspots. This could result in the resurgence of malaria-endemic areas, which is a major threat to the elimination of malaria. Therefore, it is necessary to monitor the genetic structure of the local P. falciparum population, especially during the time when sustaining the control e cacy is the top priority.
Antimalarial drug resistance is one of the key biological challenges in malaria control and elimination [14]. ACT resistance has been found to be associated with pfK13, pfmdr1, and pfcrt genetic markers [27,44]. According to the policy decision recommended by the WHO [24], ACT should be changed if the proportion of treated patients remaining parasitemic on day 3 exceeded 10%. In STP, the proportion of patients showing positive parasitemia on day 3 after treatment with ACT was only 2.5% (63/2,576), showing that the current treatment policy was acceptable. However, by monitoring target-site mutations, our study found an increase in the pfmdr1 86Y mutation rate in comparison to that in a previous study conducted in 2004, the time when ACT was introduced in STP [22]. The prevalence of the pfmdr1 86Y mutation in local parasites increased from 21-83%. This mutation is known to be associated with reduced sensitivity to AQ and chloroquine (CQ) in P. falciparum [45]. With an increase in the pfmdr1 86Y mutation and nearly xed mutation on pfcrt 76T, the local parasites may show increased tolerance to AQ, raising concerns regarding the use of AQ as an ACT partner drug in STP. Other SNPs showed the same pattern of predominance (pfmdr1 184F and pfmdr1 D1246) as shown in the parasites isolated in 2004 [22]. Overall, the prevalent haplotype of pfmdr1 86 + 184 + 1246-pfcrt 76 was YFD-T (51.4%) in STP, which was similar to the ndings on the neighboring island, Bioko island, Equatorial Guinea, where the prevalence of pfmdr1 YFD type was 59% [46]. In addition, no mutations were found in the artemisinin resistance marker pfK13 in the STP isolates. In fact, pfk13 mutations were relatively rare in P. falciparum isolated from Africa, suggesting a lower selection pressure on this genetic marker even when ACTs have been intensively used [47,48]. The contribution of pfK13 polymorphisms to artemisinin resistance remains uncertain in the African isolates.
This study found that young children had a higher risk of developing high levels of parasitemia and treatment failure, probably due to the lack of acquired immuno-protection, and the lower treatment compliance compared to adults [49]. Notably, patients treated with quinine were more likely to show insu cient clearance of parasites during follow-ups in comparison with those treated using ASAQ.
Several possible reasons can be provided for this. First, patients treated with quinine mostly had higher parasitemia levels, which may prolong the clearance time. The mean parasite clearance time after quinine treatment was approximately 4-5 days [50,51], which was slower than that with ASAQ treatment (2-3 days) [52] in general, and the difference was even greater in cases with higher parasitemia. Second, poor compliance and tolerability with the quinine regimen may occur, especially in children. The recommended quinine regimen in sub-Saharan Africa is 10 mg/kg administered three times daily for 7 days [24]. This prolonged treatment course may result in reduced patient compliance if the patients have to keep taking medicines after they are discharged from the hospital. Studies in Africa have found unacceptably high treatment failure rates for patients who did not complete the 7-day quinine regimen [53,54]. One study in Uganda showed that among children administered a 7-day course of oral quinine, 69% (18/26) of parasitological failures were cases of recrudescence, and the major cause was the poor adherence due to the care givers forgetting to administer the drugs, the drugs being vomited up, or the children feeling better [55]. This was similar to our ndings in STP, and our study results also proved that most of the recurrent cases in the quinine group involved recrudescence, with only a few new infections. The nal reason was the possible development of resistance against quinine in the local parasite population. The mechanism of quinine resistance has not been well elucidated [54]. Con icting results from the lack of resistance to varying degrees of resistance against quinine have been reported in Africa [56,57]. Notably, one study from Western Kenya showed that polymorphisms in pfmdr1 86Y, 184F, and pfcrt 76T were signi cantly associated with reduced quinine susceptibility [58]. This suggests that the parasite strains in STP may have developed reduced sensitivity in quinine owing to the high prevalence of pfmdr1 86Y, 184F, and pfcrt 76T found in our study. However, the contribution of these genetic markers differed among strains, and thus, the association with in vitro quinine susceptibility should be further assessed in the local parasite strains [58]. A study conducted in STP in early 1992-1993 showed no resistance to quinine in in vivo or in vitro tests [50]. Since then, quinine has been used as a standard treatment for over 30 years, and no updated study has been conducted to date. Therefore, the therapeutic e cacy of quinine should be carefully inspected in STP. Overall, considering the higher recurrence rate after quinine treatment and possible resistance mutations noted in this study, an alternative treatment protocol using intravenous artesunate is recommended for severe malaria in accordance with WHO guidelines [24]. An ACT-based post-discharge oral regimen for 3 days may increase patient adherence. If the local health care facilities encountered stock-outs of ACT, a 7-day course of quinine in combination with antibiotics should be strictly followed, and the follow-up assessments should be strengthened, at least for children.
Lastly, in combination with our previous ndings for vector mosquitoes, which showed increased mosquito density and knockdown resistance mutations in An. coluzzii [7], and signi cant changes were observed in the antigenic genes of P. falciparum from 2012 to 2013, which occurred during and after a recent peak incidence. This evidence showed that under a low transmission scenario, dynamic changes in vectors and parasites could result in uctuations in the malaria incidence rate. Undoubtedly, biological changes in vectors and parasites are related to local interventions, spatiotemporal evolution and adaptation, and human activities in the environment [59]. Therefore, routine surveillance of vectors, parasites, and clinical cases plays an important role in mapping malaria risk in low-transmission settings [60]. Through the integration of the case surveillance data, one can capture a more comprehensive risk characterization and disease mapping. Our previous and current studies pointed out that outdoor mosquito density and insecticide resistance in An. coluzzii, changes in antigenic and antimalarial drug resistance markers in P. falciparum, and recurrence in patients were signi cant risk factors for local malaria transmission. Moreover, surveillance data showed that malaria hotspots in STP were distributed along the primary and secondary roads around the Central Hospital (HAM), military camps, schools, markets, and airports located in the capital district of AG. These are places that aggregate many people and may intensify malaria transmission. Although the higher malaria case numbers observed in these hotspots may be due to the accessibility to hospitals, insu cient environmental management in these densely populated locations pose a great threat to the rebound of malaria. With respect to malaria elimination, these speci ed risk factors and hotspots should be addressed in the policy decision-making for national malaria control programs in STP.

Conclusion
Although malaria treatment e cacy remained acceptable in STP, the evidence from analyses of MSP1 and MSP2 in this study showed that the circulating strains of local parasites underwent temporal changes and an increase in drug-resistance mutations. Case surveillance data showed that patients with younger age, higher parasitemia level at enrollment, and quinine treatment were more likely to experience recurrence during follow-ups. The therapeutic e cacy of antimalarials and case surveillance should be strengthened, especially in the hotspots, to avoid resurgence of malaria transmission and the spread of resistant parasites on this island that had pre-eliminated malaria. Significant p-values (p < 0.05) are shown in Italic. The definition of each variable is described in method section. The mediation effect, in which parasitemia leads to recurrence through treatment are examined in Figure S6. Figure 1 Classi cation of periods based on the annual malaria incidence rate in STP. Malaria incidence rate substantially dropped through effective indoor residual spraying (IRS), and the introduction of artemisinin combination therapy (ACT) during the control phase from 2000 to 2007. Between 2007 and 2012, the incidence rate was controlled under 5% with the implementation of larval control by Bacillus thuringiensis israelensis (Bti), and rapid diagnostic tests (RDT) for mass screening. Subsequently, there was a slight increase in the incidence rate, and a rebounded peak incidence (6.9%) was shown in 2012. From 2012 to 2014, malaria incidence was again controlled and declined. The country then moved to the preelimination phase with incidence around 1% after 2014, and the electronic medical record (EMR) system had been executed in the health care facilities.   Period B and C are low-transmission phases, of which the latter period shows a recent peak incidence in 2012, and Period D is the pre-elimination phase (see Figure  1). MSP1 haplotypes identi ed in P. falciparum isolates were KH1-KH9 (yellow colors, K1 family), MH1-MH8 (blue colors, MAD20 family), and RH1-RH4 (green colors, RO33 family). In MSP2, 11 haplotypes were found in the 3D7/IC family (IH1-IH11, orange colors), and 10 haplotypes were found in the FC27 family (FH1-FH10, gray colors).