Study design and sample collection
In 2021, malaria patients aged 2–10 years old were recruited as part of the National Health Facility Survey (NMCP, Ministry of Health) that evaluated quality of malaria case management at public outpatient clinics in six provinces of Mozambique (Maputo, Inhambane, Manica, Zambézia, Niassa and Nampula). 40 health facilities were randomized out of all eligible centers to maximize representativity (Fig. 1). In 2022, targeted sampling was conducted at selected health facilities in Maputo province (all ages > 6 months) and Inhambane, Manica, Zambézia, Sofala and Manica (children 2–10 years old). Samples from Nampula were obtained from children aged 3 months to 5 years old at health facilities in areas where seasonal malaria chemoprevention was implemented. Samples from Tete and Cabo Delgado were obtained from the 2022 therapeutic efficacy survey (children aged 6 months to 5 years old), with an additional inclusion criterion of parasite density > 2500 parasites per microliter by microscopy20. Additional surveys were conducted in 2022 during dry season (June to October) for Maputo and Manica provinces to study the effect of malaria seasonality on molecular markers of antimalarial resistance and parasite diversity.
In all surveys, patients presenting with a confirmed diagnosis of uncomplicated P. falciparum malaria by RDT were invited to participate20. After they provided informed consent, two to four 50 µL DBS were prepared onto one or two filter papers through finger prick. DBS were identified with anonymous barcodes, air-dried during 72 h and stored at 4 °C in sealed bags with silica gel until laboratory processing.
Genomic DNA extraction and quantification
Genomic DNA was extracted from DBS samples using a Tween-Chelex based protocol as described by Brokhattingen et al.30, with some modifications. Five mm (~ 12,5 µL blood) discs were cut from each DBS into 96-well deep well plates with a manual puncher. One mL of freshly made 0.5% Tween 20® detergent diluted in PBS was added to each well plate containing a DBS punch and incubated overnight in a thermomixer at 15ºC and 300 rpm. The next morning, the supernatant was removed, 1 mL of fresh PBS was added per well and the plate was briefly vortexed and then incubated at 4ºC for 30 min. After incubation, the liquid was aspirated, and 150 µL of a solution of 10% Chelex (C7901, Merck) in molecular grade water was added. The samples were incubated at 95ºC in a water bath for 15 min with gentle vortexing every 5 min. The plate was then centrifuged for 5 min at 1500 rpm to pellet the Chelex® beads. Supernatant (approximately 130 µL) containing the eluted DNA was transferred to a new PCR 96 well plate, centrifuged again and finally 100 uL transferred to barcoded tubes placed on 96 well plate (Wilmut).
P. falciparum infection was confirmed in all DNA samples by qPCR targeting the 18 S rRNA gene on an ABI PRISM 7500 HT Real-Time System (Applied Biosystems), as previously described41. Parasite density was quantified by extrapolation to an external standard curve composed of six 1:10 dilutions of 3D7 (MRA-151, MR4, Bei Resources) cultured parasites in whole blood spotted onto filter paper (range 100.000 to 1 parasites/µl). DNA was stored at -20ºC until sequencing.
Amplicon-based sequencing
Sequencing was performed using the MAD4HatTeR multiplex amplicon sequencing panel using CleanPlex reagents and CleanMag Magnetic Beads (Paragon Genomic Inc, California, USA) as previously described22. Briefly, we used two multiplexed PCR reactions with primer pools D1.1, R1.2 (reaction 1) and R2.1 (reaction2). Together, these primer pools target 241 P. falciparum loci of 225–300 bp22. Multiplex PCR was performed with 10 cycles if parasite density ≥ 500 parasites/µL, or 20 cycles for those samples with < 500 parasites/µL. Following multiplex PCR, reactions proceeded in a single tube for each sample, where products were bead-cleaned, digested, and indexed via PCR to generate Illumina-compatible libraries. A randomly selected subset of 10 libraries from each full plate was assessed using automated capillary electrophoresis in a TapeStation 4150 (Agilent technologies, California, USA) to confirm library quality (size and concentration). Finally, libraries from each sample were pooled adjusting volumes based on parasitemia, and the pool was bead-cleaned using a 1X bead ratio to remove primer dimers. Pooled libraries were run on an agarose gel, from which the amplicon-sized band was excised. DNA was extracted using Monarch® DNA Gel Extraction Kit (New England Biolabs Inc., Massachusetts, USA), and products were quantified using a TapeStation and a Qubit fluorometer. For 2021 samples, pools contained 96 samples and were 150 paired-end sequenced in a MiSeq System with v2-300 cycles reagents (Illumina, USA) at CISM laboratory; for 2022 samples, pools of 288 samples were also sequenced with 150 paired-end reads in a NextSeq 2000 System using P1 reagents at ISGlobal laboratory.
Positive (n = 2, matching the parasitemia category of the samples in the plate) and negative (n = 2) controls were included in every library preparation plate to control for run quality and contaminations. Positive controls were prepared from P. falciparum laboratory strains 3D7 (MRA-151), HB3 (MRA-155), Dd2 (MRA-156 and MRA-1255). Cultures were synchronized in the ring stage, mixed with uninfected human whole blood to obtain a range of parasite densities (1, 10, 100, 1,000, 10,000 and 100.000 parasites/µL) and spotted onto filter paper. Negative controls were prepared from P. falciparum negative DBS.
Sequence data analysis
FASTQ files were subjected to filtering, demultiplexing and allele inference using a Nextflow-based pipeline version 0.1.8 (https://github.com/EPPIcenter/mad4hatter)22. The 3D7 genome sequence was used as reference for alternative allele calling (https://github.com/EPPIcenter/mad4hatter/blob/main/resources/v4/ALL_refseq.fa). The resulting allele tables were subsequently filtered based on read counts and coverage across loci within a sample and across samples. Alleles with fewer reads than the maximum observed reads in any locus for negative controls were removed, along with alleles with < 1% within-sample frequency.
Prevalence of antimalarial drug resistance markers was defined as the number of P. falciparum infections carrying mutant alleles (including mixed infections) out of total infections with valid allele calls per each locus. Reconstruction of pfdhps double, pfdhfr triple and pfdhfr/pfdhps quintuple haplotypes was done for samples with no mixed genotypes at selected loci and for those with allele frequency > 95%. Within-host and population level metrics of diversity were calculated from diversity locus (n = 165) as recently described30. Samples with < 50 diversity loci covered at a read depth of at least 100, and diversity loci with < 100 samples covering them with at least 100 reads were filtered out. Intra-host complexity of infection (COI) and effective COI (eCOI) were estimated from polyallelic genomics data using a Markov Chain Monte Carlo (MCMC)-based approach implemented in MOIRE v3.4.0 (https://github.com/EPPIcenter/moire). eCOI considers within-host relatedness, and can be interpreted as the expected COI if population diversity was infinite (heterozygosity = 1). Polyclonal infections were defined as having eCOI > 1.130. Wright’s inbreeding co-efficient (Fws) was calculated as 1-Fws across all diversity loci as the allele heterozygosity of the individual (HW) relative to the population42. Genetic diversity of the parasite population was measured using expected heterozygosity (HE), i.e., the probability that two randomly selected parasites carry distinct alleles at each diversity locus30. Principal Coordinates Analyses (PCoA) of in-sample allele frequencies and presence/absence of alleles were done to visualize the differences between provinces and regions. Bray-Curtis dissimilarity matrices were computed using the R package vegan (https://vegandevs.github.io/vegan/)43. A permutational Mantel test (10000 permutations) utilizing Spearman’s correlation method was conducted to evaluate the correlation between genetic (Bray-Curtis dissimilarity) and geographic distances using vegan. Geographic coordinates for each province were used to calculate a Haversine matrix.
To infer the evolutionary history of the mutant alleles in pfdhps, we focused on six specific amplicons surrounding the pfdhps gene covering positions 549,583 to 596,266 on chromosome 8 (549583–549807, 549960–550215, 550057–550318, 585331–585590, 585703–585949, and 585993–586266). Each amplicon was aligned separately with the R package msa44 and concatenated into a 1143-nucleotide sequence. Pairwise distances between infections were calculated, and an initial unrooted tree was constructed using the minimum evolution algorithm45. The best evolutionary model was identified using the Bayesian Information Criterion (BIC), and 1000 bootstrap replicates were performed for statistical support. Phylogenetic reconstruction was conducted with the Phangorn package in R46. Chi-squared tests assessed associations between pfdhps-436 codon (mutant or wild type) and pfdhps-437/540 codons, between haplotypes (wt/wt/wt, wt/mut/mut, and mut/wt/wt) and populations (region or province) and between frequency of each haplotype across regions and provinces. Cramér’s V was used to measure effect sizes. Haplotypes mix/wt/wt and mix/mix/wt were considered mut/wt/wt, and the rest of mixed haplotypes were excluded from the analysis. Single mutants for dhps437/540 were grouped into one category for the codon association analysis. The population association analysis only included haplotypes wt/wt/wt, wt/mut/mut, and mut/wt/wt.
Copy number of pfpm2
MAD2HatTeR data was screened for CNV at pfpm2 locus using read counts. A generalized additive model was constructed based on amplicon read counts for each sequencing run under the assumption that most amplicons correspond to single-copy genes. Fold change values were calculated from the difference between observed and expected read counts, and further normalized using the resulting fold change of single-copy laboratory 3D7 controls sequenced in each run to account for technical variation. Additionally, the 5% of amplicons with the highest read count variability across sequencing runs were excluded from the fold change calculation to reduce noise. Samples with normalized fold changes greater than 1.5 were flagged as potential CNV. The analysis was restricted to samples with a parasite density > 1000 parasites/µL, as lower parasitaemia samples are more prone to produce inconsistent read counts across runs22.
Samples with read fold-change > 1.5 were tested for pfpm2 CNV adapting a previously published qPCR protocol15 by using ubiquitin conjugated enzyme (pfuce) as reference gene47. Reactions were set-up in an ABIPrism 7500 Thermocycler with Power SYBR Green master mix (ThermoFisher). Amplification efficiencies calculated from standard curve build with 3D7 genomic DNA were 0.90% for pfpm2 and 0.81% for pfuce. Copy numbers were calculated using ddCt method and 3D7 (one pfpm2 copy) as calibrator sample. Samples with CNV > 1.5 by qPCR were considered pfpm2 duplications. DBS from a field isolate with confirmed four pfpm2 copies were included as positive control (kindly donated by Prof. Didier Ménard, Institute Pasteur).
Definitions and statistical analysis
Regions for geographical analysis were South (Maputo, Inhambane), Centre (Manica, Sofala, Zambézia, Tete) and North (Nampula, Niassa and Cabo Delgado). Maputo Province and Maputo City, which are administratively independent provinces, were considered as a single unit for sampling purposes and in secondary analysis at the province level. Maps were created using Python 3.9 from the public OpenStreetMap data.
The prevalence of infections carrying parasites with markers of antimalarial resistance were estimated per region, year and transmission season. Validated and associated markers of antimalarial drug resistance were considered as those included in WHO’s 2020 review report on antimalarial drug-resistance13. Rainy season was defined as the period between January 1st and May 30th, and samples collected after this date and up to October 31st were considered dry season. Malaria case data per each province and season was extracted from the Health Information System for Monitoring and Evaluation (SIS-MA, Ministry of Health, Mozambique). Health facilities reporting a minimum of 96 malaria rapid diagnostic tests per study month were included to avoid biases for health facility with sporadic testing patters (95% confidence interval and 80% power). Positivity rate was defined as the percentage of positive malaria RDTs out of the total tests conducted in each of the health facilities.
Statistical analyses were performed in Stata version 15.0 or R version 4.3.1. Chi-square test, Fisher’s exact were used to compare frequencies of resistance markers and differences in percentage of polyclonal infections by location and/or year. Kruskal–Wallis rank sum test was used for the comparison of the distribution of parasite densities, read counts, eMOI or 1-Fws between populations. Multivariate regression models for genetic metrics were build adjusting for the age (older/younger than 5 years), gender and parasite density (below or above 500 parasites/µL. Genome-wide He estimate and 95% confidence intervals were calculated by a linear mixed model fitting locus as a random effect. To compare the malaria positivity rates between the rainy and dry seasons within each province, the Wilcoxon rank sum test was employed. P-values lower than 0.05 was considered to indicate statistical significance.
Ethical considerations
Clinical-demographic data and blood samples were collected only after written informed consent was obtained from participants and/or their accompanying adults. All RDT positive individuals were treated according to national treatment guidelines. Study protocols were approved by the Mozambican National Committee for Bioethics in Health (CNBS, references number 354/CNBS/2021 and 604/CNBS/21) and Hospital Clinic de Barcelona Ethics Review Committee (ref. HCB/2022/0097).
