Assessing Onchocerca volvulus Intensity of Infection and Genetic Diversity Using Mitochondrial Genome Sequencing of Single Microfilariae Obtained before and after Ivermectin Treatment

Onchocerciasis is a neglected tropical disease targeted for elimination using ivermectin mass administration. Ivermectin kills the microfilariae and temporarily arrests microfilariae production by the macrofilariae. We genotyped 436 microfilariae from 10 people each in Ituri, Democratic Republic of the Congo (DRC), and Maridi County, South Sudan, collected before and 4–5 months after ivermectin treatment. Population genetic analyses identified 52 and 103 mitochondrial DNA haplotypes among the microfilariae from DRC and South Sudan, respectively, with few haplotypes shared between people. The percentage of genotype-based correct assignment to person within DRC was ~88% and within South Sudan ~64%. Rarefaction and extrapolation analysis showed that the genetic diversity in DRC, and even more so in South Sudan, was captured incompletely. The results indicate that the per-person adult worm burden is likely higher in South Sudan than DRC. Analyses of haplotype data from a subsample (n = 4) did not discriminate genetically between pre- and post-treatment microfilariae, confirming that post-treatment microfilariae are not the result of new infections. With appropriate sampling, mitochondrial haplotype analysis could help monitor changes in the number of macrofilariae in a population as a result of treatment, identify cases of potential treatment failure, and detect new infections as an indicator of continuing transmission.


Introduction
Over 120 million people are at risk of contracting onchocerciasis, a neglected tropical disease (NTD) caused by the filarial nematode Onchocerca volvulus and transmitted by blackflies in the genus Simulium, primarily in sub-Saharan Africa [1]. Onchocerciasis pathology includes skin depigmentation and lesions, ocular damage that can lead to blindness, and neurological-endocrine disorders such as epilepsy and Nakalanga syndrome [2][3][4]. The World Health Organization (WHO), in conjunction with government and non-government agencies, are targeting the elimination of O. volvulus transmission in endemic countries [5]. The strategy for interrupting transmission is the mass drug administration of ivermectin (MDAi). Ivermectin kills the larval stages (microfilariae) found in the skin of infected persons (microfilaricidal effect), and temporarily stops female adult worm reproduction (embryostatic effect). It is distributed to affected communities annually, biannually, or even quarterly. With each subsequent exposure to ivermectin, the fertility of female worms may be reduced, although the estimated degree of this cumulative effect ranges across studies [6][7][8][9]. The effects of long-term MDAi on skin microfilariae density and macrofilariae fertility have resulted in the WHO-verified elimination of transmission in most endemic areas in the Americas and may have also resulted in elimination in several endemic areas in sub-Saharan Africa [10][11][12][13][14]. These successes motivated shifting the goal from control as a public health problem to the elimination of onchocerciasis transmission [5,15].
One challenge facing the elimination of the transmission of O. volvulus with MDAi in sub-Saharan Africa is the variable response of adult worms to the embryostatic effect of ivermectin, with increases in skin microfilariae density post treatment beginning as early as a few weeks or as late as >1 year [16][17][18][19][20][21][22][23][24][25]. Repopulation of the skin with microfilariae 80 days post treatment, or the observation of stretched microfilariae in the uteri of adult female worms 80-90 days post treatment, has been referred to as a sub-optimal response (SOR) to ivermectin [16][17][18][19][20][21][22][23][24]26]. In addition, a sub-optimal microfilariae response (SOMR) has been described, in which the rate or extent of the initial reduction in skin microfilariae following ivermectin treatment is reduced [24,25,27]. The microfilaricidal and embryostatic effect of ivermectin may differ between geographic areas [20,24]. Thus, tracking community-specific changes in microfilaridermia and worm fertility is useful for the quantitative assessment of the effect of MDAi and/or for parameterizing models used to predict the optimal timing and duration of MDA for the successful elimination of transmission [28,29].
Current parasitological methods for assessing the effect of MDAi on an individual or community level [16,17,21,22] do not allow the estimation of the number of fertile adult worms in a host or in a community. Palpable subcutaneous nodules containing adult worms can be excised for analysis, but an unknown (and significant) number of worms are inaccessible deep in the body [30][31][32][33][34], and the number of individuals who can be sampled is limited by the fact that minor surgery is required. Estimates of the number of live female adult worms per person living in the West African savannah based on the statistical analyses of macrofilariae in excised nodules indicate high variability among hosts within the same community, ranging from approximately 4 to 177 [35,36] worms per person.
Our goal is to develop a tool that can quantitively measure progress towards reducing the number of fertile worms during MDAi (and complementary interventions). Because maternal sibling microfilariae sampled from the same person will have identical mitochondrial DNA sequences (or haplotypes) due to the strictly maternal nature of mitochondrial inheritance, we can estimate how many adult females (macrofilariae) are reproductively active by counting the number of unique mitochondrial haplotypes identified in the microfilariae sampled from the skin, assuming sister nematodes rarely infect the same person. Because a particular skin snip may not contain microfilariae offspring from all female worms in the person, the number of haplotypes detected will represent the minimum number of reproductively active female worms. However, assuming that the subset represents a random sample of the microfilariae in the host, changes over time in the number of haplotypes estimated should indicate whether MDAi is effectively reducing the number of reproductively active worms. More explicitly, the number of haplotypes present amongst the microfilariae should decline, as MDAi reduces transmission (and thus new infections) and adult females reach the end of their reproductive life span, thus providing a measure of progress towards the endpoint of zero fertile adult females.
We opportunistically utilized microfilariae collected in two studies whose focus was on onchocerciasis-associated epilepsy, one in the Logo Health zone in the Ituri province of the Democratic Republic of the Congo (DRC) and the other in Maridi county of South Sudan. People enrolled in the studies had epilepsy (potentially onchocerciasis-associated epilepsy) and were ivermectin-naïve in the Logo Health zone [37]. To estimate how many microfilariae per person would need to be sequenced to determine the minimum number of reproductively active adult female worms, we performed rarefaction and extrapolation analysis of mitochondrial haplotype diversity. We compared these estimates between sampling locations to explore differences in baseline infection intensity. We compared haplotype diversity between the samples from the two countries and between microfilariae obtained pre-and post-treatment for four people from whom both pre-and post-ivermectin microfilariae were available. Finally, we placed mitochondrial haplotype diversity identified in the microfilariae from DRC and South Sudan in the context of sequence data from parasites from other endemic areas in sub-Saharan Africa.

Origin of Microfilariae
Skin snips were collected from people with epilepsy in Maridi County in South Sudan and in the Logo Health Zone, Ituri Province, DRC in 2018 during epidemiological studies to investigate the association between onchocerciasis and epilepsy. Two snips were obtained from each participant at two time points: immediately before ivermectin treatment ("pretreatment") and 4 months (DRC) or 5 months (South Sudan) later ("post-treatment"). After the 24-h incubation of the snips in isotonic saline, the microfilaridermia was quantified microscopically [38][39][40] and samples were transferred to 80% ethanol and sent to La Trobe University, Victoria, Australia. In Maridi County, annual MDAi had been instituted since the early 2000s, but was interrupted for several years because of insecurity, and only reintroduced in 2017 with very low coverage (40.8%) [41]. At the time the skin snip was taken, each participant had either never taken ivermectin or had taken it only once. In the Logo Health zone, DRC, MDAi had never been instituted. A recent study concluded that today, Simulium dentulosum appears to be the main vector of human onchocerciasis in the area, and that Simulium vorax may be a secondary vector [42]. Ov16 IgG4 antibody positivity, determined with the Ov16 IgG4 Bioline rapid diagnostic test [43][44][45][46], between 2016 and 2018 was 0% (0/55) among 6-year-old and 7.1% (13/182) among 7-10-year-old children in the Logo Health zone [47]. In contrast, the Ov16 seroprevalence, determined with the same rapid test in an area in Maridi close to the blackfly breeding site, was 40% among children 3-6 years old and 66.7% among children 7-9 years old [48]. This suggests much higher ongoing O. volvulus transmission in Maridi compared to the Logo Health zone.

Origin of Adult Worms and Adult Worm Sequences
In 2016-2017, 27 adult female worms from the Centre Region and 12 from the Littoral Region in Cameroon were excised from nodules, as described in [20], and the heads (i.e., without uterine tissue) were placed in RNAlater (Invitrogen, Thermo Fisher Scientific, Waltham, MA, USA).

DNA Extraction, Amplification and Sequencing
To minimize the transfer of host skin cell debris during DNA extraction, the microfilariaecontaining ethanol solution was washed with Milli-Q ® water (Merck, KGaA, Darmstadt, Germany) by carefully transferring the ethanol solution containing the microfilariae, but without skin snip, into a 15 mL falcon tube, filling to 15 mL with Milli-Q water, and centrifuging at 500× g for 5 min. Approximately 14 mL was aspirated carefully from the tube, leaving the microfilarial pellet in approximately 1ml of now diluted ethanol solution. This was transferred into a clean reusable glass petri dish (glass because microfilariae are less adhesive to glass compared to plastic). Each microfilaria was picked by aspiration using a pipette (0.1-2 µL) under a dissecting microscope into 20 µL lysis buffer (10mM Tris-HCl, pH 8.0; 1mM EDTA, pH 8.0; add 1% Tween ® 20 (Sigma-Aldrich, Burlington, MA, USA) and Proteinase K 300 µg/mL (recombinant, PCR Grade, Sigma-Aldrich) added just prior to use). The transfer of any visible cell debris was avoided while picking microfilariae. Tubes each containing a microfilaria in lysis buffer were incubated at 55 • C for 2 or 3 h, followed by heat inactivating the solution at 80 • C or 85 • C for 20 min. DNA from adult worms from Cameroon was extracted using the Isolate II Genomic DNA kit (Bioline, London, UK) as per the manufacturer's instructions.
Quantitative PCR (qPCR) was performed on the microfilarial lysates to confirm the presence of mitochondrial DNA prior to whole genome amplification. qPCR reactions were performed using 2 µL of 1:5 diluted microfilarial lysates, 5 µL of SsoAdvanced Universal SYBR Green master mix (Bio-Rad Laboratories, Hercules, CA, USA), 2 µL of nuclease-free water, and 0.5 µL each of 10 µM forward (SP-Ov-mt-10062: 5 -attggtgaccaataaccttca-3 ) and reverse (ASP-Ov-mt-10062: 5 -ttgattcaatatcagggacgta-3 ) primers. A synthesized 68-bp oligonucleotide (ttg att caa tat cag gga cgt ata ttt cgt caa tct gag ttg act ttg aag gtt att ggt cac caa t) was used as a positive control and standard and HPLC water as a negative control for all reactions. Oligonucleotides were synthesized by Integrated DNA Technologies (Redwood City, CA, USA). qPCR assays were run on a CFX Real-Time System (Bio-Rad Laboratories), with an initial denaturing step of 3 min at 98 • C followed by 40 cycles (of 98 • C 10 s, 54 • C 15 s, 72 • C 15 s, read plate) including melt curve analysis at 65 • C to 95 • C for 5 s with an increment of 0.2 • C for 5 s. Each of the diluted lysates was assayed in duplicate, and the standards in triplicate. The overall statistics for the qPCR runs were assessed using the CFX Maestro Software (Bio-Rad Laboratories). The cycle threshold (Cq) values for each sample were determined as positive if the Cq < 30, and negative if Cq > 30 [51,52].
Microfilariae are relatively small in size (250-360 × 5-9 µm [29]), with low yield of DNA from a single microfilaria. A minimum amount of starting DNA is prerequisite to the successful generation of Illumina sequencing libraries. Whole genome amplification (WGA) was used as an intermediate step on microfilarial DNA lysates to generate high yields of amplified DNA as required for library construction. To generate a sufficient quantity of microfilarial DNA for Illumina library construction, high-fidelity, multiple displacement whole genome amplification (WGA) was performed on each microfilaria using 2 µL lysate as the starting material and processed according to the manufacturer's instructions (REPLIg, QIAGEN GmbH, Hilden, Germany). The WGA reactions were performed at 30 • C for 16 h and heat-inactivated at 65 • C for 10 min. Concentrations of the WGA DNA were determined using a Qubit 2.0 fluorometer (Thermo Fisher Scientific, Waltham, MA, USA).
Sequencing libraries were constructed on samples whose DNA concentrations were above 1.0 ng/µL, using the Nextera DNA Flex/Illumina DNA prep library construction kits and barcoded using unique dual indexes, according to the manufacturer's instructions (Illumina, San Diego, CA, USA). An Agilent TapeStation analysis was performed on libraries using D1000 ScreenTapes (Agilent Technologies, Santa Clara, CA, USA) to ensure that the libraries were within the selected size range of 400-650 bp. The Qubit fluorometer was used to quantify the libraries for pooling. The final 4nM pooled library was spiked with 1% PhiX control (Illumina) and run on a NovaSeq SP, 300 cycles (i.e., 150 bp paired reads) at the Australian Genomic Research Facility (Melbourne, Australia) or on a NextSeq 500, 300 cycles, at the La Trobe Genomics Platform (Bundoora, Australia).

Variant Calling
Raw reads were trimmed using trimmomatic v.0.32 [53]. Trimmed reads were competitively mapped to the O. volvulus nuclear and mitochondrial genomes v.4 [54], the Wolbachia bacterial endosymbiont (GenBank accession number NZ_HG810405.1), and Homo sapiens (genome GRCh38.p13) using bwa v.0.7.17 [55,56]. Reads that mapped to H. sapiens were counted for quality control assessment (Supplementary Table S1) and discarded from the data using Unix command line programs grep and awk. Samtools v.1.9 [57] was used to remove secondary and supplementary reads, and only reads that mapped uniquely to O. volvulus were retained. If a microfilaria was sequenced in more than one experiment, the mapped reads were combined. Depth was assessed using bedtools v.2.26.0 [58] and coverage for each O. volvulus chromosome at a depth of at least 20× was estimated using a custom Perl script (Supplementary Table S1). Samples that did not meet a minimum coverage of 85% at a depth of 20× were not included in the downstream analyses. RepeatScout v1.0.5 [59] was used to identify repetitive regions across O. volvulus chromosomes.
Variants were also called using freebayes v.1.0.2 [60], and filtered for quality and depth using vcftools v.0.1.13 [61]. Variant calls were normalized using bcftools v.1.2 [62] and haplotypes were simplified using the vcfallelicprimitives function in vcflib [63]. Bcftools was then used to find the intersection between the calls made by GATK and those made by freebayes. Finally, vcftools was used to further filter the dataset to remove individuals with less than 99% of variants called, to remove sites with missing data, and to identify and remove indels, singletons, and repeat regions.

Data Analysis
The R package adegenet v.2.1.5 was used to produce principal component analysis (PCA) plots of the genetic diversity within people at each of the different collection time points and to cluster the worms [64,65]. Discriminant analysis of principal components (DAPC [66]), which maximizes the differences between predefined groups, was used to estimate the probability that a microfilaria would be correctly identified to person or to date of collection based on genotype. The number of PCs to use in DAPC was estimated using cross-validation.
SNPEff v.4.1l [72] was used to assess whether variants were in coding or noncoding regions and to calculate their effects on known genes (e.g., amino acid changes).
PGDSpider v.2.1.1.3 for Windows [73] was used to convert vcf to nexus format, which was then edited to add a "trait" block indicating person or country of origin. PopArt v.1.7 [74] was used to generate haplotype networks using the TCS approach [75], calculate diversity statistics, and to perform an analysis of molecular variance (AMOVA; [76]).

Sequencing Results
Out of 804 sequenced in total, 225 microfilariae from 10 people from the DRC and 211 microfilariae from 10 people from South Sudan successfully passed the filtering criteria and were included in the analyzed data (Table 1). Of the successfully sequenced microfilariae from the DRC, 207 were collected before treatment (D0); 15 and 3 microfilariae collected four months post treatment, were, respectively, from two persons from whom we do not have pre-treatment mtDNA sequences. The dataset from South Sudan includes 163 microfilariae collected before treatment and 48 collected around five months post treatment, with four people having parasites successfully sequenced from both pre-and post-treatment microfilariae samples.

Genetic Variation of Microfilariae within People
After merging and filtering, 439 variants (places in the genome where the base differs from the reference sequence) in the 13,744 bp mitochondrial genome were called, of which 119 were singletons (i.e., found in only one worm). Because singletons could be due to experimental error (e.g., the introduction of variation during whole genome amplification) rather than biological variation, we excluded these when determining the number of mitochondrial haplotypes represented in each sample, leaving 320 variants.
The number of mitochondrial haplotypes identified for the microfilariae from a person represents the minimum number of reproductively active female worms in that person (Table 1), since there may be unidentified haplotypes from microfilariae either not present in the skin snip, lost during transfer of the skin snip to ethanol, or not successfully sequenced. Table 1. Number of Onchocerca volvulus mitochondrial haplotypes estimated per person in the study areas in the Democratic Republic of Congo (DRC) and South Sudan. The mean counts of microfilariae from two skin snips, the number of microfilariae successfully sequenced and in the final dataset, and the number of haplotypes before and after treatment are given. The total number of parasite haplotypes per person are also indicated, as this number will be different from the sum of the pretreatment and post-treatment haplotypes when haplotypes are found in both samples. Totals for each country represent the total number of unique pre-and post-treatment haplotypes used in the analysis; the total number in parentheses indicates the number of unique haplotypes across all worms, as some haplotypes were found in more than one person. A rarefaction curve indicates how the number of haplotypes identified changes given the number of successfully sequenced microfilariae, and the shape of the curve indicates the extent to which the number of microfilariae sequenced effectively captures the genetic diversity within the population investigated, i.e., the microfilariae within an individual person (host). The rarefaction curves for samples from 12 individuals (those with microfilariae counts ≥10; Figure 1) suggest that, overall, the number of microfilariae successfully sequenced captured the genetic diversity in the parasite population in the DRC study area (green) to a larger extent than that in the South Sudan study area (purple). Extrapolated estimates for the total number of haplotypes predicted per person further suggest that the predicted number of reproductively active female worms (observed+unobserved) is greater in hosts from the South Sudan study area (average across people based on Chao estimate: 29.4 haplotypes per person; based on abundance-based accumulation: 36.4) when compared to the DRC (9.18 and 13.4, respectively), with highly variable associated standard errors ( A rarefaction curve indicates how the number of haplotypes identified changes given the number of successfully sequenced microfilariae, and the shape of the curve indicates the extent to which the number of microfilariae sequenced effectively captures the genetic diversity within the population investigated, i.e., the microfilariae within an individual person (host). The rarefaction curves for samples from 12 individuals (those with microfilariae counts ≥10; Figure 1) suggest that, overall, the number of microfilariae successfully sequenced captured the genetic diversity in the parasite population in the DRC study area (green) to a larger extent than that in the South Sudan study area (purple). Extrapolated estimates for the total number of haplotypes predicted per person further suggest that the predicted number of reproductively active female worms (observed+unobserved) is greater in hosts from the South Sudan study area (average across people based on Chao estimate: 29.4 haplotypes per person; based on abundancebased accumulation: 36.4) when compared to the DRC (9.18 and 13.4, respectively), with highly variable associated standard errors ( Table 2).  We used several different methods for assessing and comparing mtDNA genetic diversity in the microfilariae from the DRC and South Sudan locations. Because the sampling of the parasite population was non-random (i.e., the repeated sampling of sibling microfilariae in each host is non-representative of the population of parasites across all hosts), we estimated genetic diversity statistics using the set of unique maternal haplotypes rather than individual microfilariae sequences. Nucleotide diversity (π: the average number of pairwise nucleotide differences per site) of haplotypes from the DRC was 0.0198, and that of haplotypes from South Sudan was similar at 0.0212. Tajima's D (the difference in South Sudan. This indicates that in the DRC, most of this variation was partitioned within individual people (90.59%), with minor variation among people (9.41%). In South Sudan, this pattern of variation was even stronger (within individuals: 96.73%; between individuals 3.27%). In other words, the genetic differentiation between the microfilariae present in different people was marginally higher in the DRC than in South Sudan, but in both populations, nearly all of the genetic diversity was within the individual hosts (as in other parasites, such as the filarial nematode that causes lymphatic filariasis, Wuchereria bancrofti [77,78]). We visualized the relationships among the microfilariae using haplotype networks, which indicate how similar specific haplotypes are to each other. There are fewer haplotypes observed across the 225 microfilariae from 10 people from the DRC (Figure 2a) than observed across the 211 microfilariae from 10 people in South Sudan (Figure 2b), and the network is correspondingly simpler. There are two haplotypes that form central nodes in the South Sudan haplotype network that are shared by microfilariae from multiple people (Figure 2b).
We also explored whether genetic diversity was clustered by infrapopulation (the microfilariae within an individual host). The PCA of the samples from DRC suggests that there is considerable overlap in genetic diversity across microfilariae from different people (Figure 3a), but with only a few shared haplotypes (Figure 2a). DAPC is an analytical approach that maximizes differentiation among groups. A DAPC of the mitochondrial variants from DRC (using the first 80 PCs as determined by cross-validation) indicated that worms sampled from different people could be genetically differentiated (Figure 3b). While there are a few shared haplotypes, the proportion of individual microfilariae that could be correctly assigned to their host, based on their mitochondrial genotype, is 0.8756. The assignment proportion does not necessarily serve as a metric for similarity among worms within a person; rather, it indicates how different the worms sampled are between people. The higher the assignment proportion, the more distinct the parasite genotypes present in different people. We also explored whether genetic diversity was clustered by infrapopulation (the microfilariae within an individual host). The PCA of the samples from DRC suggests that there is considerable overlap in genetic diversity across microfilariae from different people (Figure 3a), but with only a few shared haplotypes (Figure 2a). DAPC is an analytical approach that maximizes differentiation among groups. A DAPC of the mitochondrial variants from DRC (using the first 80 PCs as determined by cross-validation) indicated that worms sampled from different people could be genetically differentiated (Figure 3b). similarly suggests overlap in genetic diversity across the microfilariae in different hosts (Figure 3c), with a couple of haplotypes found in several people (Figure 2b). DAPC (using the first 60 PCs) was less able to discriminate among the parasites from different people compared to the DAPC performed on the DRC samples (Figure 3d), and the assignment proportion for individual microfilariae to the infected person from which they came was only 0.641, indicating that discriminating infrapopulations based on maternal mitochondrial haplotype was not feasible. We focused further analysis on the 98 microfilarial haplotypes in the four people from South Sudan from whom we had both pre-and post-ivermectin treatment samples. It was not possible to discriminate haplotypes based on whether they were collected prior to or following treatment: haplotypes found before treatment (gray) were also found posttreatment (black) and some of the haplotypes were found in more than one person. Figure  4a visually demonstrates this genetic overlap between microfilariae haplotypes, and Figure 4b indicates the inability of the discriminant function to differentiate between preand post-treatment haplotypes: each haplotype identified (columns) has a similar probability of being assigned to either the pre-treatment (gray) or post-treatment (black) The PCA of mitochondrial genetic variation in the South Sudan microfilariae similarly suggests overlap in genetic diversity across the microfilariae in different hosts (Figure 3c), with a couple of haplotypes found in several people (Figure 2b). DAPC (using the first 60 PCs) was less able to discriminate among the parasites from different people compared to the DAPC performed on the DRC samples (Figure 3d), and the assignment proportion for individual microfilariae to the infected person from which they came was only 0.641, indicating that discriminating infrapopulations based on maternal mitochondrial haplotype was not feasible.
We focused further analysis on the 98 microfilarial haplotypes in the four people from South Sudan from whom we had both pre-and post-ivermectin treatment samples. It was not possible to discriminate haplotypes based on whether they were collected prior to or following treatment: haplotypes found before treatment (gray) were also found posttreatment (black) and some of the haplotypes were found in more than one person. Figure 4a visually demonstrates this genetic overlap between microfilariae haplotypes, and Figure 4b indicates the inability of the discriminant function to differentiate between pre-and posttreatment haplotypes: each haplotype identified (columns) has a similar probability of being assigned to either the pre-treatment (gray) or post-treatment (black) category based on their genotype. When discriminating based on human host and the day sampled, the post-treatment haplotypes largely overlapped with the pre-treatment haplotypes sampled from the same individual, although haplotypes from different individuals were moderately differentiated (Figure 4c,d).
Pathogens 2022, 11, x FOR PEER REVIEW 11 of 20 category based on their genotype. When discriminating based on human host and the day sampled, the post-treatment haplotypes largely overlapped with the pre-treatment haplotypes sampled from the same individual, although haplotypes from different individuals were moderately differentiated (Figure 4c,d).

Genetic Variation in the African Context
Maternal haplotype sequences from the DRC and South Sudan mf data (i.e., retaining only one copy of a haplotype if it was found in multiple people and/or at multiple time points) were combined with sequences from Guinea (n = 4), Sierra Leone (n = 3), Liberia (n = 1), Mali (n = 10), Côte d'Ivoire (n = 13), Ghana (n = 189), Benin (n = 1), Cameroon (n = 39), and Uganda (n = 2) [49,50] for a total of 438 O. volvulus. After filtering, there were 851 single-nucleotide polymorphic variants called in total; of these, 427 were singletons. For downstream analyses, the singletons were removed.
Of the 424 polymorphic variants in the 13,744 bp mitochondrial genome, 139 variants were nonsynonymous (i.e., resulting in a change in the amino acid in the gene product), 208 were synonymous (i.e., a nonsynonymous:synonymous ratio across all genes of 0.668), 4 were nonsense mutations that would truncate the predicted gene product, and the remaining were in non-coding regions. The overall nucleotide diversity was 0.0111. The reduction in nucleotide diversity compared to that estimated from microfilariae from the DRC and South Sudan was largely driven by low nucleotide diversity among the 221 West African worms (π = 0.0064). Most of the variance in this diversity was within countries rather than between them: a simple non-nested AMOVA resulted in a fixation index of ФST = 0.0255 (p = 0.009; significance based on 1000 permutations), where only 2.55% variance was among countries and 97.45% was within country.

Genetic Variation in the African Context
Maternal haplotype sequences from the DRC and South Sudan mf data (i.e., retaining only one copy of a haplotype if it was found in multiple people and/or at multiple time points) were combined with sequences from Guinea (n = 4), Sierra Leone (n = 3), Liberia (n = 1), Mali (n = 10), Côte d'Ivoire (n = 13), Ghana (n = 189), Benin (n = 1), Cameroon (n = 39), and Uganda (n = 2) [49,50] for a total of 438 O. volvulus. After filtering, there were 851 single-nucleotide polymorphic variants called in total; of these, 427 were singletons. For downstream analyses, the singletons were removed.
Of the 424 polymorphic variants in the 13,744 bp mitochondrial genome, 139 variants were nonsynonymous (i.e., resulting in a change in the amino acid in the gene product), 208 were synonymous (i.e., a nonsynonymous:synonymous ratio across all genes of 0.668), 4 were nonsense mutations that would truncate the predicted gene product, and the remaining were in non-coding regions. The overall nucleotide diversity was 0.0111. The reduction in nucleotide diversity compared to that estimated from microfilariae from the DRC and South Sudan was largely driven by low nucleotide diversity among the 221 West African worms (π = 0.0064). Most of the variance in this diversity was within countries rather than between them: a simple non-nested AMOVA resulted in a fixation index of Φ ST = 0.0255 (p = 0.009; significance based on 1000 permutations), where only 2.55% variance was among countries and 97.45% was within country.
Based on PCA, the haplotype diversity observed in South Sudan and DRC contained some sequence variation that overlapped with mitochondrial sequences obtained from adult female worms from West and Central Africa, and some sequence variation that was quite distinct (Figure 5a; consistent with the haplotype network presented in Figure S1). The results of a DAPC including only parasites from those countries where at least 10 parasites were sampled (Cameroon, Côte d'Ivoire, DRC, Ghana, Mali, and South Sudan), using the first 70 PCAs as determined through cross-validation, was visually consistent with the results of the PCA (Figure 5b). Assignment proportions (the proportion of parasites correctly assigned to their country of origin) were higher for parasites from Cameroon (0.90), DRC (0.94), and South Sudan (0.81) than parasites from Ghana (0.73), Côte d'Ivoire (0.08), and Mali (0.40) (Figure 5c). The very low proportions for Côte d'Ivoire and Mali could be due to fewer haplotypes being available compared to the number of sequences from countries with higher assignment proportions (i.e., an analytical bias). Consideration has to be given to the possibility that parasites in countries in West Africa are more closely related overall (due to the migration of hosts or vectors and/or the more recent divergence of those populations in evolutionary history) and thus less able to be discriminated using mitochondrial data alone.
was among countries and 97.45% was within country.
Based on PCA, the haplotype diversity observed in South Sudan and DRC contained some sequence variation that overlapped with mitochondrial sequences obtained from adult female worms from West and Central Africa, and some sequence variation that was quite distinct (Figure 5a; consistent with the haplotype network presented in Figure S1). The results of a DAPC including only parasites from those countries where at least 10 parasites were sampled (Cameroon, Côte d'Ivoire, DRC, Ghana, Mali, and South Sudan), using the first 70 PCAs as determined through cross-validation, was visually consistent with the results of the PCA (Figure 5b). Assignment proportions (the proportion of parasites correctly assigned to their country of origin) were higher for parasites from Cameroon (0.90), DRC (0.94), and South Sudan (0.81) than parasites from Ghana (0.73), Côte d'Ivoire (0.08), and Mali (0.40) (Figure 5c). The very low proportions for Côte d'Ivoire and Mali could be due to fewer haplotypes being available compared to the number of sequences from countries with higher assignment proportions (i.e., an analytical bias). Consideration has to be given to the possibility that parasites in countries in West Africa are more closely related overall (due to the migration of hosts or vectors and/or the more recent divergence of those populations in evolutionary history) and thus less able to be discriminated using mitochondrial data alone.

Discussion
Our results showed that the O. volvulus mitochondrial DNA haplotype diversity varied between people and that the number of haplotypes in a single person was in the range of 2-26. The number of haplotypes indicates the minimum number of adult female worms that produced the microfilariae in our analyses, and is within the range of previous estimates of the number of female adult worms per person in West African populations [35,36]. Furthermore, our analyses showed that: (1) the number of haplotypes identified increased with the number of microfilariae successfully sequenced (Figure 1), indicating that data from a higher number of microfilariae than was available to us is required to estimate the number of haplotypes within one host as well as within the human population to which the host belongs; (2) the total number of microfilariae that need to be successfully sequenced to capture the genomic diversity of their parents will vary depending on the number of adult parasites within a person and geographic location (Figure 1, Table  2); and, thus, that (3) the number of haplotypes we identified underestimates the number

Discussion
Our results showed that the O. volvulus mitochondrial DNA haplotype diversity varied between people and that the number of haplotypes in a single person was in the range of 2-26. The number of haplotypes indicates the minimum number of adult female worms that produced the microfilariae in our analyses, and is within the range of previous estimates of the number of female adult worms per person in West African populations [35,36]. Furthermore, our analyses showed that: (1) the number of haplotypes identified increased with the number of microfilariae successfully sequenced (Figure 1), indicating that data from a higher number of microfilariae than was available to us is required to estimate the number of haplotypes within one host as well as within the human population to which the host belongs; (2) the total number of microfilariae that need to be successfully sequenced to capture the genomic diversity of their parents will vary depending on the number of adult parasites within a person and geographic location ( Figure 1, Table 2); and, thus, that (3) the number of haplotypes we identified underestimates the number of reproductively active female worms in most people, more so in those from South Sudan than those from DRC. If a rarefaction curve begins to asymptote, then the microfilariae genotyped is more likely to be representative of the genetic diversity within that person. In South Sudan, the number of microfilariae per person that would need to be successfully sequenced to capture the genetic diversity in that parasite population appears to be higher than for people in the DRC (Figure 1). These differences were observed despite similar numbers of microfilariae sequenced from a similar number of hosts, which further suggests that the intensity of infection, or worm burden per person, is higher in the area sampled in South Sudan than in DRC. The differences in infection intensity may be related to differences in prevalence; in the cohorts studied here, 36.5% of people recruited for the study based in Ituri compared to 84.9% of participants from South Sudan were microfilaridermic [38,40]. However, because the participants did not represent a random sampling of onchocerciasis-infected people in either location, sampling parasites from additional people would be required to confirm what would be a significant epidemiological difference between the two areas. Differences in onchocerciasis prevalence are driven by differences in annual biting rates, which are, in turn, driven by ecological conditions favorable for blackfly breeding, the amount of time people spend in areas with different biting rates, and the presence and effectiveness of interventions [79][80][81].
The statistically significant fixation index (Φ ST from AMOVA) of haplotypes from the DRC and South Sudan demonstrates that each individual infected person carries a genetically variable population of worms. In population genetic terms, each infected person carries their own infrapopulation that is drawn from a much larger pool of parasites present in the population as a whole (the metapopulation). Most of the genetic variation is partitioned within people, which means that the differentiation between infrapopulations is moderate, and stronger in the DRC than South Sudan. The difference between the DRC and South Sudan infrapopulation differentiation is driven by fewer haplotypes (i.e., fewer reproductively active females) detected per host from the DRC, and thus that each sampled host contained a smaller proportion of the total parasite population's genetic diversity than hosts in South Sudan. Regardless, the higher within-than between-infrapopulation diversity supports the hypothesis that the parasite populations in each transmission zone-the geographic area over which parasites are transmitted and thus able to interbreed-are large, and the probability of parasites that have identical-by-descent mitochondrial haplotypes being transmitted to one person is low.
The O. volvulus populations sampled in DRC and South Sudan harbor populationassociated genetic variants that can be used to discriminate many of the parasites collected from these locations from parasites found elsewhere in sub-Saharan Africa. As shown in similar analyses of mitochondrial DNA data from O. volvulus [49,78], there is a widespread haplotype found in nearly all countries, which causes the overlap in genetic diversity observed ( Figure 5; central node in Figure S1). This is consistent with large population sizes maintaining ancestral variation in mitochondrial genomes, which are under selective constraint to retain function.
For the four hosts from South Sudan from whom we had mitochondrial DNA sequences from microfilariae collected before as well as five months post-ivermectin treatment, we found that the haplotypes of the microfilariae collected post treatment were not distinct from those taken before treatment (Figure 4b). This shared genetic diversity indicates that the microfilariae present before and five months after ivermectin treatment are likely the offspring of those mothers that recovered fertility relatively rapidly. However, some haplotypes identified post treatment were not found in the pre-treatment sample. New post-treatment infections can be discounted as a source of microfilariae in the skin only five months post treatment, given that it takes 12-18 months for L3 larvae transmitted to mature into a reproductive adult worm and that it might take up to 3 years for the microfilariae of that adult to be sufficiently numerous to be detectable in skin snips [32,[82][83][84]. An alternative source for post-treatment mitochondrial haplotypes not detected in a pretreatment sample could be the random assortment of unique haplotypes into offspring from a heteroplasmic mother (i.e., a female worm with more than one mitochondrial haplotype). A likely interpretation is that sampling was not sufficient to detect all of the genetic diversity actually present within the host: some genotypes are missing from the pre-treatment sample because of under-sampling rather than true absence. We tested this hypothesis by producing a rarefaction curve, which indicates how the detection of unique mitochondrial sequences (or haplotypes) found in each person in both the DRC and South Sudan microfilariae increases with additional sequenced microfilariae (Figure 1). Under the assumption that post-treatment microfilariae would be a genetic subset of the pre-treatment worms, the results are consistent with the hypothesis that there was insufficient sampling of microfilariae from people in the South Sudan cohort to capture the genetic diversity present at the time of ivermectin treatment. Extrapolation may be useful for guiding researchers in determining when additional sequencing is required; while the extrapolated number of haplotypes varies depending on the assumptions made about detection of rare haplotypes [70,71], we found that in as many as three people (M206, M224, M238), only about half of the mitochondrial diversity was predicted to have been sampled ( Table 2).
The analysis of parasite mitochondrial data to track changes in the number of reproductively active female worms could be useful for elimination programs and researchers working where there has been persistent transmission despite MDAi with high compliance and community coverage. High vector-biting rates can sustain transmission, even when skin microfilariae load is low. However, ongoing transmission might also be driven by parasites transmitted via vector or human movement from other onchocerciasis endemic areas or by SOR to ivermectin (reviewed in [85]). Changes in the estimated number of reproductively active worms over time could serve as an indicator of whether years of programmatic MDAi are successfully reducing the worm burden, or whether new imported infections (represented as new haplotypes) are identified. In this latter case, the genotyping of worms from areas with ongoing transmission that are likely to be connected by movement of people (as in [86,87]) or vectors (including long-distance, wind-assisted migration [88][89][90][91]) could help determine whether transmission among endemic areas might be contributing to persistent prevalence (see, e.g., [92,93] for modeling that explores the impact of movement of people or vectors on prevalence). We have shown here that mitochondrial haplotypes can discriminate among parasites from different countries, and, thus, may be informative where cross-border transmission might occur. Nuclear genotypes derived from the whole-genome sequencing of adult worms were able to discriminate O. volvulus collected in forest vs. savannah bioclimes within West Africa [49]; thus, the further development of cost-effective approaches for the nuclear genotyping of microfilariae would be useful for identifying transmission occurring between endemic areas.
In areas where SOR is suspected or has been demonstrated, monitoring changes in the number of reproductive SOR macrofilariae can help inform decisions about whether to deploy certain alternative and/or complementary interventions. SOR macrofilariae resume reproduction as soon as a few weeks after ivermectin treatment. Assuming there is no difference in the growth, survival, and probability of transmission of microfilariae that are the offspring of SOR and non-SOR macrofilariae, and assuming that the early resumption of reproduction is a heritable trait [20], the prevalence of reproducing SOR parasites will increase over time. This would jeopardize onchocerciasis control and elimination efforts and would require alternative intervention strategies. To date, the methods for identifying SOR parasites involve counting the microfilariae in skin snips and/or evaluating developmental stages in the uteri of macrofilariae excised from palpable subcutaneous nodules (where a fraction of the macrofilariae reside) soon (typically around 3 months) after ivermectin treatment [16][17][18][19][20][21][22][23][24]26]. Neither method can identify the number or percentage of SOR macrofilariae. A genetic approach applied to longitudinal, post-treatment samples of skin microfilariae or parasites in the vectors could indicate whether or not there is an increase in the number of SOR adult female worms producing those microfilariae repopulating the skin early, and, thus, whether alternative intervention strategies should be considered.
A major challenge for developing a genetics-based tool is the expense and sampling effort required to sufficiently capture the genetic diversity of reproducing parasites. A costeffective approach would minimize the amount of effort needed to obtain microfilariae from skin snips, be able to sequence DNA from parasite pools, and ideally be suitable for use with blackfly pools. Long-read sequencing technologies (such as Oxford Nanopore) applied to amplified targets in the mitochondrial or nuclear genome would have the advantage of sequencing haplotypes even from DNA pools. Since national NTD programs are currently not set up for sequencing and population genetic analysis, collaboration with a research institution in the country, or building the required infrastructure and personnel capacity within the programs, would be needed.
Finally, the methods we present here may be useful, and easiest to implement, in studies evaluating the efficacy of new drugs, when it is important to know whether microfilariae appearing in the skin are due to macrofilariae reproductively active before treatment and that were not sterilized or killed by the new drug or whether they are due to new infections acquired after treatment [94]. In drug efficacy trials, four skin snips are usually taken before treatment and at varying time points after treatment. The mitochondrial haplotypes of the microfilariae sampled before and after treatment (timing dependent on the prior knowledge about the effect of the drug on the parasites) can be compared to estimate the effect of a drug on the burden of reproductively active female worms and on individual female worm fertility. Depending on the extent to which sampling captures haplotype diversity, samples taken >1 year after treatment could provide insight into the probability that the haplotypes only identified in post-treatment samples are due to new infections rather than due to parasites not affected by or having recovered from the effect of the drug. This is particularly important when the studies are conducted in areas with high transmission that could result in post-treatment infection.  NIH and ERC had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results. WHO/TDR, via co-author A.C.K., had a role in the conceptualization of the research, the writing of the manuscript, and in the decision to publish the results.

Data Availability Statement:
The data presented in this study are openly available from the National Center for Biotechnology Information Short Read Archive (https://www.ncbi.nlm.nih.gov/sra/) under study accession number PRJNA981587.