{'fraction_non_alphanumeric': 0.0314698809488466, 'fraction_numerical': 0.02015853694742155, 'mean_word_length': 5.240088921822897, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
R H: a genetic metric for measuring intrahost Plasmodium falciparum relatedness and distinguishing cotransmission from superinfection

Abstract Multiple-strain (polygenomic) infections are a ubiquitous feature of Plasmodium falciparum parasite population genetics. Under simple assumptions of superinfection, polygenomic infections are hypothesized to be the result of multiple infectious bites. As a result, polygenomic infections have been used as evidence of repeat exposure and used to derive genetic metrics associated with high transmission intensity. However, not all polygenomic infections are the result of multiple infectious bites. Some result from the transmission of multiple, genetically related strains during a single infectious bite (cotransmission). Superinfection and cotransmission represent two distinct transmission processes, and distinguishing between the two could improve inferences regarding parasite transmission intensity. Here, we describe a new metric, RH, that utilizes the correlation in allelic state (heterozygosity) within polygenomic infections to estimate the likelihood that the observed complexity resulted from either superinfection or cotransmission. RH is flexible and can be applied to any type of genetic data. As a proof of concept, we used RH to quantify polygenomic relatedness and estimate cotransmission and superinfection rates from a set of 1,758 malaria infections genotyped with a 24 single nucleotide polymorphism (SNP) molecular barcode. Contrary to expectation, we found that cotransmission was responsible for a significant fraction of 43% to 53% of the polygenomic infections collected in three distinct epidemiological regions in Senegal. The prediction that polygenomic infections frequently result from cotransmission stresses the need to incorporate estimates of relatedness within polygenomic infections to ensure the accuracy of genomic epidemiology surveillance data for informing public health activities.


Introduction
The past two decades have seen an increase in the collection and application of pathogen genetics in public health surveillance programs. Genomic epidemiology surveillance is now used for a wide variety of viral, bacterial, and eukaryotic pathogens (1)(2)(3)(4)(5). For Plasmodium falciparum, a major focus of genomic epidemiology surveillance has been to identify genetic markers associated with parasite transmission. Collectively, these surveys have revealed sev-eral genetic markers associated with parasite transmission (6)(7)(8)(9). Chief among them is the frequency of clonal (genetically identical) parasites in the population, the frequency of multistrain (polygenomic) infections, and the number of strains per infection [commonly referred to as complexity of infection (COI)] (10)(11)(12).
Each of these metrics represents different aspects of parasite transmission. The increased frequency of clonal parasites is hypothesized to reflect decreases in transmission due to a reduc- Under standard assumptions of superinfection, individuals in high-transmission settings are exposed to multiple infectious bites and have multiple opportunities for outcrossing. This results in parasite populations characterized by high proportion of polygenomic infections, high COIs, and few clones. Individuals in low-transmission settings are exposed to fewer infectious bites, resulting in fewer polygenomic infections and low COIs. Low transmission also limits outcrossing and promotes the transmission of clonal parasite strains. (B) Cotransmission weakens the relationships proposed in panel (A) because polygenomic infections can also result from a single bite. Cotransmission is initiated when a mosquito ingests multiple parasite strains from an initial polygenomic infection. The number of strains present in the initial mosquito blood feed depends on the COI and strain proportions in the initial host. These parasites undergo sexual outcrossing to produce genetically related sporozoites that are transmitted into a new host. Cotransmissions do not reflect multiple exposures and instead represent a single exposure event.
tion in parasite genetic diversity and a reduction in outbreeding opportunities. The frequency of polygenomic infections and their COIs are hypothesized to reflect the frequency of superinfection, or the repeated infection of individuals from multiple, infectious mosquito bites (Fig. 1A). If true, genetic metrics could augment traditional, but labor-intensive estimates of transmission, such as the entomological inoculation rate (EIR, the number of infectious mosquito bites per individual) (13,14). To date, COI is one of the most frequently reported genetic metrics of parasite transmission (5,12,(15)(16)(17)(18)(19), and a variety of statistical tools have been developed to estimate COI from various sources of genotypic data (20)(21)(22). COI is positively associated with transmission intensity and typically ranges between 1 and 5, but extreme values >10 have occasionally been reported (23).
However, polygenomic infections are not always the result of superinfection. Whole genome sequencing analyses have detected genetically related parasites in both single-strain (monogenomic) and polygenomic infections (10,(23)(24)(25)(26). These genetically related parasites are created when a mosquito vector bites a polygenomic infection and ingests multiple parasite strains that then mate and undergo sexual reproduction to complete the sexual stage of their life cycle (24,27) (Fig. 1B). Cotransmission occurs when the mosquito vector bites a new host and transmits multiple parasite strains that are oftentimes genetically related. The genetic composition of the resulting polygenomic infection reflects the assortment and segregation of genomes during meiosis, the sampling of gametocytes by the mosquito vector, the injection of sporozoites into the human host, and the invasion of sporozoites into the liver (27). Thus, failing to accurately distinguish superinfections from cotransmissions could lead to inaccurate estimates of malaria transmission based on COI alone.
One way of distinguishing superinfection from cotransmission is to examine the relatedness of coinfecting strains within polygenomic infections. Superinfection assumes these strains are ran-domly drawn from the population, while cotransmission assumes they are genetically related. While several metrics have been developed to assess parasite relatedness, they have largely required whole genome sequencing and/or the establishment of genomic phase (22,28). Here, we developed a new metric, R H , that utilizes intrahost heterozygosity to assess the relatedness of malaria infections and distinguish cotransmission from superinfection. Both superinfection and cotransmission contribute to the complexity of parasite population genetics, and distinguishing the two could lead to greater insight into how transmission and parasite genetics are linked. As a proof of concept, we used R H to quantify the genetic relatedness of parasites genotyped with a 24 single nucleotide polymorphism (SNP) barcode (29). Using R H , we show that cotransmission is prevalent in three epidemiologically distinct regions of Senegal with moderate-to-low transmission intensity.

Genomic epidemiology of Kédougou, Thiès, and Richard Toll
A total of 1,758 malaria infections were collected and genotyped using a TaqMan-based molecular barcode ( Fig. 2; see the "Materials and methods" section). These samples were collected from three regions of Senegal: Kédougou, Thiès, and Richard Toll. Kédougou has the highest intensity of transmission (Fig. 2B), with incidences of 300 to 500 + cases per thousand reported by the National Malaria Control Program (PNLP, Programme National de Lutte contre le Paludisme) (30)(31)(32)(33)(34)(35) between 2014 and 2020. Thiès and Richard Toll have lower transmission intensity and incidences of <10 and <1 case per thousand. Based on incidence data, we expected Kédougou to have the least amount of clonal population structure, the highest proportion of polygenomic infections, and the highest COI. The allele frequencies observed for each of the assays used in the molecular barcode are presented in Fig. S1. Polygenomic fraction and COI were elevated in Kédougou ( Fig. 3 and Fig. S2A). The inverse-variance-weighted averages for polygenomic fraction (the proportion of infections polygenomic) were 0.53 (0.50, 0.57) for Kédougou, 0.27 (0.24, 0.31) for Thiès, and 0.33 (0.28, 0.38) for Richard Toll. There was no statistically significant change associated with the sampling year in Richard Toll between 2012 and 2015 (ordinary least-squares regression, slope = −0.01, Pvalue = 0.601 for Richard Toll). For Thiès and Kédougou, we found a small, but statistically significant increase in the polygenomic fraction from 2015 to 2020 at a rate of 0.03 (0.010, 0.04) per year (P-value = 0.01) for Thiès and 0.020 (0.001, 0.035) per year for Kédougou (P-value = 0.04).
As expected, we observed the least amount of clonal population structure in Kédougou ( Fig. 3B and C). We detected a total of 7 clonal haplotypes in Kédougou, 23 in Richard Toll, and 38 in Thiès. We found no evidence of persistent clonal haplotypes in Kédougou (Fig. 3C). Two of the clonal haplotypes in Richard Toll and six of the clonal haplotypes in Thiès were detected in multiple years ( Fig. 3F and I). Most clonal haplotypes were detected twice, but some haplotypes were detected five or more times per year in Thiès and Richard Toll.
Our point estimates for the unique monogenomic fraction were 0.95 (0.92, 0.98) for Kédougou, 0.70 (0.66, 0.75) for Thiès, and 0.86 (0.83, 0.90) for Richard Toll (Fig. 3B, E, and H and Fig.  S2B). These were calculated as the inverse-variance-weighted average across all years. We found little evidence that the unique monogenomic fraction changed significantly over time in Kédougou between 2015 and 2019 or Richard Toll between 2012 and 2015 (ordinary least-squares regression, P-value >0.23). Interestingly, we detected a small increase of 0.07 (0.038, 0.102) per year (P-value = 0.004) in the unique monogenomic fraction in Thiès, but the increase was negligible when compared with

Verifying the accuracy of R H for identifying cotransmission and superinfection based on samples genotyped with the 24 SNP molecular barcode
The main goal of this study was to determine if genetic metrics that measure intrahost heterozygosity and inbreeding could distinguish superinfection from cotransmission. Theoretically, an R H > 0 reflects cotransmission, an R H = 0 reflects superinfection with COI = 2, and an R H < 0 reflects superinfection with COI ≥ 3. However, it was unclear whether the 24 SNP molecular barcodes would bias R H estimates and result in errors in its interpretation. To address this, we simulated polygenomic infections genotyped with the 24 SNP molecular barcode under different hypotheses of superinfection and cotransmission ( "Materials and methods" section, Supplementary Material).
These simulations showed that an R H quantified from the 24 SNP molecular barcode could be used to distinguish superinfection from cotransmission (Fig. 4A). Superinfection with a COI = 2 resulted in a theoretical expectation (R H ) of zero. Superinfection with COI > 2 resulted in negative R H values of −0.5 for a COI = 3 infection and −1.0 for a COI = 4 infection. Cotransmission resulted in positive R H values even when the COI of the initial infection was high. For cotransmission chains originating from a COI = 2 superinfection, R H increased from 0.40 to 060 and 0.70 following the first, second, and third cotransmission event.
However, these simulations also showed that individual superinfection and cotransmission events could result in positive and negative R H values (Fig. 4A). To address this, we used the Bayes factor to define a classification scheme that would (1) distinguish cotransmission from superinfection (Fig. S3C), (2) identify serial cotransmission ( Fig. S3D) or (3) estimate the COI of superinfections ( Fig. S3E and F). Based on these results, we defined cotransmission as R H > 0.30, superinfection with COI = 2 as -0.1 < R H < 0.3, and superinfection with COI > 2 as R H < -0.1 (Table 1).
Of key interest was to determine whether using a threshold of R H > 0.30 accurately discriminated cotransmission from superinfection. With this threshold, the true cotransmission identification rates were 0.66 (0.63, 0.69), 0.85 (0.83, 0.87), and 0.93 (0.92, 0.95) for the first, second, and third cotransmission events following a COI = 2 superinfection (Fig. 4B). These rates decreased to 0.47 Table 1. Classification of polygenomic infections as cotransmission or superinfection using R H .)

R H Classification
The boundaries and classification notations are based on the standards for Bayes factor classification described in ref. (36)  We next evaluated the performance of this threshold using three diagnostic tests designed to evaluate how accurately a given threshold distinguishes cotransmission and serial cotransmissions from COI = 2 superinfections (Fig. S4). These tests showed that the optimum threshold differs for cotransmission and serial cotransmissions and depends on the COI of the originating polygenomic infection (Fig. S4B to D). In general, the optimal threshold is lower for cotransmission events originating from high COI infections. In the context of Senegal, where the average COI was 2, our tests show that a cotransmission threshold of R H > 0.30 is near the optimum threshold needed to accurately detect single and serial cotransmission events (Fig. S4B). Note that the true optimum is somewhat uncertain, as identifying it requires priors for the relative proportion of cotransmissions, and serial cotransmissions in the population.

R H versus F WS
We compared R H with another commonly used estimate of intrahost heterozygosity, F WS (37, 38) ( "Materials and methods" section, Supplementary Material). However, because the original definition of F WS reported in (37) required the use of reads generated from next generation sequencing data, we could only compare R H and F WS for the 2020 samples where barcode and Illumina short-read whole genome sequencing data were available for both monogenomic and polygenomic infections. Barcode-derived R H estimates were consistent with their corresponding whole genome sequencing-derived F WS estimates (R-squared = 0.66, ordinary-least squares) (Fig. S5).
To examine the theoretical differences between R H and F WS , we simulated barcode-derived R H and F WS in populations with different levels of clonal parasite sharing. For these simulations, we relied on the F WS described in (39)(40)(41) that normalizes intrahost heterozygosity with the expected, population-level heterozygosity but relaxes the requirement that heterozygosity be estimated from sequence counts. In populations with no clonal parasites, R H and F WS can both be used to identify cotransmissions and evaluate COI (Fig. S6). However, as the frequency of clones increases and the total heterozygosity of the population declines, the expected F WS values for simulated cotransmissions and superinfections shift and deviate from the expectations observed in nonclonal populations (Fig. S6). Simulated R H values were unaffected. The first cotransmission always had an expected R H of 0.5, while COI = 2 and COI = 3 superinfections always had an expected R H of 0.0 and −0.5, respectively (Fig. S6).

R H detects widespread cotransmission in Kédougou, Thiès, and Richard Toll
We next calculated R H for all the polygenomic infections collected from Kédougou, Thiès, and Richard Toll. Individual polygenomic samples had a wide range in individual R H values, with some samples with R H values close to one and some samples with negative R H values (Fig. 5A to C). Despite this variation, the average R H ( R H ) in each population was positive but inversely associated with regional incidence estimates. The R H was 0.14 (0. 12 Based on the R H thresholds for classifying superinfection and cotransmission defined in Table 1, we found that most polygenomic infections were the result of cotransmission (Fig. 5D to F). Despite its relatively high incidence, we found that 0.43 (0.33, 0.54) of the polygenomic infections from Kédougou were the most likely result of cotransmission. This was smaller than either Thiès or Richard Toll, where the proportion was 0.53 (0.38, 0.69) and 0.52 (0.45, 0.60), respectively. These cotransmitted polygenomic infections were classified as either COI = 1 or COI = 2 by the THE REAL McCOIL (Fig. S7). Approximately half of the THE REAL McCOIL COI = 2 samples had R H values consistent with cotransmission.
These elevated R H values could also result from a technical inability to accurately detect heterozygous sites (29). From the laboratory-generated 3D7, Dd2, and TM90C6B mixtures, we knew that heterozygous site detection was reduced in polygenomic mixtures with extreme strain ratios (Supplementary Material). We adjusted the heterozygosity of each polygenomic barcode based on the strain proportions inferred from the cycle thresholds (CTs) reported by the TaqMan-based barcode assays (Supplementary Material Figs. S8 and S9). This adjustment resulted in reduced R H values but had no effect on our overall results (Fig. S10).

Discussion
To date, most genetic epidemiology analyses have ignored cotransmission and interpreted polygenomic fraction and COI under simple assumptions of superinfection. However, it is uncertain how well this assumption holds as transmission falls in moderateor near-elimination settings where superinfection is expected to be comparatively rare. Resolving this conundrum requires collecting additional data regarding polygenomic relatedness and cotransmission rates across diverse epidemiological settings.
To enable broad-scale genetic surveillance of polygenomic relatedness and cotransmission in parasite populations, we developed a new metric, R H , that normalizes the observed polygenomic heterozygosity with the expected heterozygosity of a COI = 2 superinfection as a possible means for distinguishing superinfection from cotransmission in parasite populations. Polygenomic infections with less heterozygosity than a typical COI = 2 superinfection indicate cotransmission, while polygenomic infections with more indicate a superinfection whose COI is greater than 2.
The normalization of observed polygenomic heterozygosity to the expected heterozygosity of a COI = 2 superinfection is critical and differentiates R H from previous metrics of intrahost heterozygosity such as F WS (37-41). F WS normalizes intrahost heterozygos- ity to the total, population-level heterozygosity of the population (H S , Eq. S10). Unlike F WS , the sampling of genetic clones is excluded from R H because they are indistinguishable from COI = 1 monogenomic infections and do not directly contribute to outcrossing or cotransmission. In highly diverse populations with few genetic clones, F WS behaves similarly to R H . However, F WS estimates are sensitive to the presence of parasite clones in the population. This complicates F WS based inferences of either cotransmission or COI because these inferences would need to be recalibrated to consider the frequency of clones in each population.
In contrast, R H is unaffected by the frequency of parasite clones, and its predictions remain consistent across populations with different clonalities. This property makes R H particularly relevant for examining parasite populations in historically low-transmission settings such as Southeast Asia and South America, where population fragmentation and clonal parasite populations are the norm (37,42,43). In fact, the frequency of parasite clones in these regions can be as extreme as the simulated populations used in Fig. S6; 60% of the P. falciparum cases collected in Quibdó, Columbia, in 2001 showed evidence of being infected by the same parasite strain (44). While the advantages of R H are clearest in low-transmission settings with highly clonal parasite populations, it can also be advantageous in moderate-and high-transmission settings to account for parasite clonality arising from falling transmission intensity or focal transmission (45).
Another advantage of R H compared to other metrics used to assess relatedness or inbreeding (22,28,37) is that R H does not require whole genome sequencing or that the genomes within polygenomic infections be phased. In fact, barcode-derived R H were consistent with whole genome sequence-derived estimates of F ws despite relying on a small number of SNPs present throughout the genome. This greatly expands our ability to assess superinfection and cotransmission in epidemiological settings using accessible genotyping technologies like the 24 SNP TAQman-based molecular barcodes where whole genome sequencing is unfeasible due to logistical or financial constraints.
Estimating R H from the 24 SNP molecular barcode also avoids some of the issues involving microsatellites (39). Use of an SNP barcode rather than microsatellites based on di-and trinucleotide repeats avoids the potential problems of replication slippage and unequal crossing over, which can alter the number of repeating units (46). SNP markers are not prone to such events, although they may change due to random mutation or gene conversion; these events are infrequent.
Using R H , we provide new evidence that cotransmission is ubiquitous among symptomatic infections reporting to clinics in three sites spread across Senegal. These results are consistent with our previous analyses of whole-genome sequences of polygenomic infections in Thiès, which also suggest that cotransmission occurs at high frequency (24). Here, we also found evidence of high rates of cotransmission in both Richard Toll and Kédougou. These results were unexpected for Kédougou, where the incidence was 300 to 500 + cases per thousand per year, as the prevailing consensus was that polygenomic infections in moderate-to hightransmission settings were the result of superinfection. However, our results are consistent with previous reports of cotransmission in polygenomic infections collected from Malawi , where transmission is high (25,26). Altogether, our results showed that at nearly half of the sampled polygenomic infections are the result of a single infectious bite. The true cotransmission rates are likely higher, based on the theoretical misclassification rates reported in Fig. 4. Regardless, these results imply that a significant fraction of polygenomic infections result from single infectious bites and that superinfection is a rarer than expected in populations with incidences less than 500 cases per thousand per year. Whether cotransmission is equally prevalent in asymptomatic infections is a topic for future research.
Our results present a problem for genetic epidemiology analyses, as they show that a large fraction of polygenomic infections do not reflect multiple bites. Accurate EIR inference will likely require model-based approaches that incorporate all available information regarding superinfection, cotransmission, host immunity, and age (47), and any potential observation or sampling biases, such as focal transmission heterogeneity, that could affect COI estimation. An important application of R H would be to enable calibration of superinfection and cotransmission rates in future model-based inferences of EIR. Superinfection and cotransmission represent two distinct types of transmission whose impact on parasite genetics may be useful for helping genomic epidemiology models disentangle different aspects of parasite transmission. However, these models should also assess the effects of host immunity on COI as immunity could also cause intrahost heterozygosity to be depressed and be conflated with cotransmission.
Interestingly, initial genomic surveillance with the 24 SNP molecular barcode suggests that metrics of parasite clonality, COI, and polygenomic fraction did not distinguish Thiès from Richard Toll, both of which have incidences less than 10 cases per thousand per year. While this could reflect a limitation of the 24 SNP molecular, we note that whole genome sequences of polygenomic infections in Thiès also reported high cotransmission frequencies (24). Alternatively, importation from migration and migrant workers in Richard Toll could be obscuring the genetic signals associated with reduced transmission intensity (48)(49)(50). High levels of importation could explain the high unique monogenomic fraction relative in Richard Toll and the anomalously low R H values observed in Thiès during 2018. Additional analyses throughout different epidemiological clines will be valuable for determining which genetic signals are most strongly associated with transmission intensity and whether anomalies like the low R H in Thiès 2018 reflect other epidemiological factors such as migration.
To guide our understanding of epidemiology and malaria population genetics, genomic epidemiology analyses should identify and separate cotransmission from superinfection polygenomic infections. R H represents a significant advancement in malaria genomic epidemiology by explicitly measures the genetic impact of sexual reproduction and outcrossing when evaluating the relationship between malaria population genetics and transmission. Understanding the differences between these two types of polygenomic infections and their relationship with parasite transmission is relevant for future model-based estimates of parasites transmission and operational decision-making by national malaria control programs.

Sample collection
Samples were collected through passive-case detection from febrile patients reporting to health posts or clinics during the malaria transmission season in Senegal (September to December), or actively detected in households in response to a case de-tected at the Richard Toll clinics. Patients over 6 months of age with fever within the past 24 hours of visiting the clinic with no history of antimalarial use were diagnosed with malaria using microscopy or rapid diagnostic tests (RDTs). Filter papers spotted with blood were collected from malaria-positive patients in Thiès and Kédougou. RDTs were collected from malaria-positive patients in Richard Toll as previously described (49). Ethical approval for the study was obtained (IRB Protocols: 16,330 and 17 to 1288 from Harvard T.H. Chan School of Public Health).

Molecular barcode genotyping and SNP calling
The molecular barcode consists of a series of 24 neutral SNPs spread across the malaria genome that are genotyped using a panel of TaqMan-based quantitative PCR genotyping assays. Nucleic acid material was extracted from either filter paper or RDT material and preamplified using the methods described in (49,51). A description of the reagents and methods used for the TaqMan assays is in (29).
Criteria for calling homozygous and heterozygous sites were determined using a set of laboratory-generated mixes containing 3D7, Dd2, and TM90C6B DNA mixed with proportions ranging from 1 : 1 to 1 : 5 (Supplementary Material Table S1). These criteria were based on the normalized CTs calculated by the Applied Biosystems ViiA 7 Real-Time PCR System (v1.2). Baselines for allele 1 and allele 2 were determined by the software, and the RN that determines CTs was set to 20.
Homozygous sites were denoted by their allelic identity (A, T, C, and G) and identified as sites with (1) no heterogeneous amplification, (2) heterogeneous amplification where the difference in CT between allele 1 and allele 2 ( CT) was greater than 8 , or (3) heterogeneous amplification, where one of the alleles had a CT > 38. For homozygous sites with heterogeneous amplification, the allelic identity was determined by the allele with the smaller CT value.
Heterozygous sites were denoted as "N" and identified as sites with heterogeneous amplification and (1) a CT < 8 and (2) CT < 38 for both alleles. Barcodes with two or more heterozygous sites were considered polygenomic. Sites without amplification were considered missing and denoted as "X." Barcodes with more than two missing sites were excluded from analysis. These thresholds were identified by benchmarking their accuracy in detecting heterozygous and homozygous sites using a set of lab-generated DNA mixtures generated from DNA isolated from 3D7, Dd2, and TM90C6B P. falciparum strains (Supplementary Material Figs. S11 to S13).

Whole genome sequencing
Monogenomic and polygenomic samples collected from 2020 were submitted for whole genome sequencing. Selective whole genome amplification was performed on extracted DNA and used to construct libraries with a NEBNext Ultra II library Kit for Illumina short-read sequencing. Variant-calling and read alignment was performed following the best practice standards set by the Pf3k consortium. Briefly, short-reads were aligned to the P. falciparum 3D7 reference genome (PlasmoDB v. 28) using BWA-mem and Picard Tools, Variants were called using HaplotypeCaller in GATK v3.5.
Analyses were carried out using a set of ∼150,000 SNPs identified from a set of 1,328 monogenomic samples obtained from the Pf3k database (52). These sites were chosen to (1) exclude the core chromosomal regions, (2) reside in nonoverlapping 2 kb windows whose average intrahost heterozygosity < 0.03, (3) be farther than seven base pairs from indels, and (4) have an average intrahost heterozygosity < 0.04 across all samples.

Inverse-variance weighting
To correct for the varying sample sizes observed in each year, we used the inverse-variance-weighted average, which is defined aŝ where n years is the number of sample years, y i is the average of the ith year, and σ 2 i is the variance of the ith year.
) where H mono,y sample is the expected heterozygosity obtained from sampling 200 monogenomic pairs from the monogenomic barcodes present in each sample year. A sampled monogenomic pair was treated as an independent one and sampling was proportional to the frequency of each monogenomic barcode haplotype in each sample year. Monogenomic pairs consisting of identical monogenomic barcode haplotypes were excluded because these would not be recognized as polygenomic infections. H poly,sample is the observed heterozygosity of the polygenomic sample. This pro-cess was repeated 200 times and the average was used as the point estimate of R H,sample .
The average R H ( R Hregion ) for each region was calculated through bootstrapping. Bootstrapped estimates of R Hyear, region were obtained by calculating the average from 200 randomly samplied (with replacement) R H,sample estimates observed in the region during the specified sample year. This calculation was repeated 200 times to obtain the final bootstrapped distribution of R Hyear, region . The average and SDs of these bootstrapped distribution were then used to evaluate the inverse-variance-weighted average.

Simulating superinfection and cotransmission
Superinfection was simulated as the random sampling of unique monogenomic barcodes. Sites that were concordant across all sampled barcodes were considered homozygous. Sites that were discordant were considered heterozygous. Sites that were missing or heterozygous in the sampled monogenomic barcodes were excluded when quantifying heterozygosity and R H . Our superinfection simulations assume mass action and do not incorporate spatial or temporal transmission heterogeneity.
Cotransmission was simulated using a model that simulates the mating and sexual recombination of parasites as they are sampled and deposited by the mosquito vector (27). The model generates identity-by-descent maps that determine which portions of the genomes in cotransmitted strains were inherited from either of the parental strains present in the original polygenomic infection. All cotransmission chains were initiated from a superinfection with two to five randomly sampled monogenomic barcodes that were assumed to be unrelated.
The parameters used in the model were described in (27). Briefly, oocyst counts were drawn from a modified two-parameter Weibull distribution defined by f oocystcount (scale, shape) = 1 + round down (Weibull (scale, shape)) , ( (4)) where the scale = 2.5 and shape = 1. The round down in Eq. 4 indicates that value drawn from the Weibull distribution is round down to the nearest integer. This results in a distribution whose median is 2, mean is 3, and interquartile range is between 1 and 4 . This distribution reflects the low oocyst intensities observed in naturally infected mosquitoes (54)(55)(56)(57). Model predictions were robust to different assumptions regarding oocyst counts (Supplementary Material, Figs. S14 and S15).
Infected hepatocyte counts were drawn from a lognormal distribution (mean = 1.8, SD = 0.8) whose sampled values were rounded down to the nearest integer unless the sample was less than 1 (58). In this case, the sampled value was set to 1. Uneven strain proportions were not simulated because their primary effect is on the probability of maintaining cotransmission chains (27).
Based on the identity-by-descent maps generated by the model, the barcodes of cotransmitted strains were generated by (i) determining which sites of the barcode were inherited from the first or second strain of the original infection, and (ii) copying the allelic identities of those parental barcode sites. After each cotransmission event, we sampled two unique, cotransmitted strains to generate polygenomic barcodes representing cotransmitted polygenomic infections with COI = 2.