Long-read assembly and comparative evidence-based reanalysis of Cryptosporidium genome sequences reveal expanded transporter repertoire and duplication of entire chromosome ends including subtelomeric regions

An improved long-read genome assembly for C. parvum (IOWA-ATCC)

The CpIRef genome assembly, generated in 2004, has only 10 physical gaps of unknown size, but it has 18,558 ambiguous bases and is missing six telomeres. Alignment of 54,882,187 Illumina 100-bp paired-end reads (Supplemental Fig. S1) to this reference sequence revealed many regions that had become collapsed during assembly (Supplemental Table S2). To resolve these issues, we generated a new PacBio + Illumina + Nanopore hybrid genome assembly (Table 1; Supplemental Fig. S1) for the C. parvum strain IOWA (ATCCPRA-67DQ), CpIA. To minimize strain variation differences, we performed our analysis on the IOWA strain. However, because there is a 14-yr window of propagation between these two isolates, and cryopreservation has only recently been developed (Jaskiewicz et al. 2018), we modified the strain name to IOWA-ATCC (CpIA).

Table 1.

Comparative Cryptosporidium genome assembly statistics

Open in new tabDownload PowerPointLink to table

The new CpIA genome assembly is compared to the current CpIRef sequence and two closely related species with different host preferences and pathogenicity, C. hominis (Ch30976) and C. tyzzeri (CtUGA55) (Table 1; Šlapeta 2013; Nader et al. 2019; Sateriale et al. 2019). These particular assemblies were selected because they are the best available. The new CpIA long-read assembly increases the genome size by 19,939 bases (∼152 kb when including new proposed subtelomeric regions) and putatively identifies all 16 telomeres. There are no gaps and no ambiguous bases. As expected, the CpIA genome sequence has diverged slightly but shares 99.93% average pairwise identity with the 2004 assembly in regions that are comparable (Supplemental Table S3). The main Cryptosporidium subtyping marker, the 60 kDa surface protein (gp60 locus subtype IIa), shows four amino acid differences (two in the serine repeat region) between CpIA and CpIRef (Supplemental Fig. S2; Supplemental Methods).

Structural differences between the C. parvum IOWA assemblies are confirmed

The 2004 CpIRef genome assembly used Sanger reads combined with available HAPPY-map data to scaffold the contigs. We compared the CpIRef and CpIA assemblies to identify potential rearrangements. Inversions and relocations were detected in Chromosomes 2, 4, and 5 (Fig. 1A). These inversions may be previous assembly artifacts or represent genuine differences between the isolates. We investigated the synteny between CpIRef, CpIA, Ch30976, and CtUGA55 and observed that C. hominis and C. tyzzeri also share the Chr 4 and Chr 5 inversions. Examination of the inverted region boundaries in CpIRef revealed regions of ambiguous nucleotide bases or physical gaps (Fig. 1A). To further investigate, PCR primers were designed to test each possible inversion arrangement in genomic DNA from C. parvum KSU-1 strain 2006 and Bunch Grass farms IOWA (CpBGF) (Fig. 1B; Supplemental Fig. S3; Supplemental Table S4). The results support the revised assembly orientation. Long-read ONT data also support the CpIA assembly (Supplemental Fig. S4). Better assemblies for the other species will be needed to determine the true level of synteny across these species.

Figure 1.

Syntenic relationships between select Cryptosporidium chromosome assemblies. (A) Synteny between Chromosomes 2, 4, and 5. Vertical black lines within a chromosome represent known physical gaps. Synteny between chromosomes is shown in pink and inversions in blue. (B) PCR validation using C. parvum KSU-1 DNA (Supplemental Table S4). Lanes 1, 2, and 3 in all gels are 1-kb ladder, positive control Dnmt2 gene, and no template control, respectively. The remaining lanes test each orientation of the left (L) and right (R) inversion boundaries. Red stars indicate the location of primers designed based on the CpIA assembly, and gray stars indicate the same on the CpIRef assembly.

Open in new tabDownload PowerPointLink to figure

Consistent structural gene annotation resolves inconsistencies and improves functional annotation

We consistently annotated and compared CpIRef, CpIA, Ch30976, and CtUGA55 which have >95% genome identity to assess differences in gene content. The new annotation for each species was generated with three de novo approaches and evidence-based manual annotation. Curation of the annotation was performed pairwise between each assembly to take full advantage of syntenic regions. Data from one species could be used to assess computational predictions in others. Using this approach, fragments of genes that were previously missing in C. hominis were identified. This approach resulted in 1926 gene alterations in CpIA when compared to CpIRef, resulting in improved functional annotation. These changes increase the overall number of predicted genes, introns (100% supported by RNA-seq data), and exons (Table 2; Supplemental Results). The average mRNA length increased. These structural fixes led to the repair of the N terminus of the methylase ortholog, Dnmt2 (Supplemental Fig. S5) as well as 523 other genes and 113 fragmented genes previously annotated as pseudogenes.

Table 2.

Reannotation summary statistics

Open in new tabDownload PowerPointLink to table

Cryptosporidium has a very compact genome sequence with 76.88% covered by protein-coding sequences (CDS). As a result, RNA-seq data, which is the best evidence for annotation, contains reads that overlap adjacent genes creating false fusions of exons belonging to different genes. Available strand-specific RNA-seq was used to characterize some of these regions, but expression data were not available for all predicted genes (87% of the annotated genes were covered); thus, genes of unknown function in close proximity on the same strand remain problematic. The expression data also revealed three putative alternative spliced genes (CPATCC_0027530; CPATCC0027960; CPATCC_0035590) and 474 potential noncoding RNAs (ncRNAs), predominantly antisense lncRNAs with differential expression as reported (Li et al. 2021).

Comparative analysis reveals few gene content differences between closely related Cryptosporidium species

There is a cluster of species, C. parvum, C. hominis, C. tyzzeri, C. meleagridis, and C. ubiquitum that are highly syntenic relative to species outside of this cluster. The syntenic species are biologically distinct and largely host-adapted with the main zoonotic exceptions being C. parvum and C. ubiquitum. A synteny analysis of the clustered species and Cryptosporidium muris as an outgroup reveals high synteny (99.4%–87%) within the cluster and only 4% synteny to C. muris (Supplemental Fig. S6; Supplemental Tables S5, S14).

The consistent annotation of the species closest to CpIA (GCA_015245375.1), Ch30976 (GCA_001483515.1), and CtUGA55 (GCA_007210665.1) permitted the analysis of differences in CDS content and CNVs. Orthology analysis revealed that 94% of the genes were conserved among all species. Of the 4008 ortholog groups identified, most gene families were maintained with a similar number of paralogs (max = 6) detected in the same ortholog group, but variation was detected among singletons (Fig. 2A; Supplemental Table S6). Some of these gene differences appear to be unique to a particular species (Supplemental Table S7). Of the 224 singletons detected, we observed only 0, 1, and 1 potential truly species-specific genes in CpIA, Ch30976, and CtUGA55, respectively, following manual inspection (Fig. 2B,D). Both species-specific genes are uncharacterized proteins. The remaining 253 singletons are detected but incomplete in the fragmented assemblies of Ch30976 and CtUGA55, appearing as split genes, frameshifts, missed calls near a gap and missing subtelomeric regions, or contig break and putative false gene predictions in small contigs (Fig. 2C). The major protein-coding gene content differences between these species are gene copy number variations and not gene presence or absence.

Figure 2.

Ortholog distribution of protein-coding genes reveals few differences between the species. (A) Venn diagram of automated protein orthology assignment between CpIA, Ch30976, and CtUGA55. (B) Venn diagram of the same orthologous genes following manual investigation and removal of putative false positives. (*) The 139 genes shared between C. hominis and C. tyzzeri in A are in complex regions with repeats and gaps and do not have enough evidence to prove their uniqueness at this stage (Supplemental Table S7). (C) Example of a false positive paralog count caused by gene fragments on different scaffolds. (D) Putative unique uncharacterized gene found on Chr 8 in CtUGA55.

Open in new tabDownload PowerPointLink to figure

To identify and assess putatively overly collapsed repetitive regions within the genome assemblies analyzed in this study, that is, repetitive regions represented by only a single repeat in the assembly, we mapped Illumina reads from CpIA to the new CpIA, CpIRef, Ch30976, and CtUGA55 genome assemblies (Supplemental Table S2; Supplemental Fig. S1). Our pipeline detected 12 compressions of at least 2× read depth and >100 bp in length in the CpIRef genome assembly compared to six in the new CpIA assembly. The six compressed regions drop to four if the three putative new subtelomeric regions proposed in this study are included (see below). The Ch30976 and CtUGA55 genome assemblies contain more than 20 compressions mostly attributed to the short reads used to generate these assemblies. The CpIA collapsed regions have two hits in regions with gene annotations in Chr 1 and Chr 2. Both genic regions are composed of rRNA genes, some uncharacterized proteins, GMP synthase, aspartate-ammonia ligase, tryptophan synthase beta, and MEDLE genes, all associated with complex subtelomeric regions discussed below. The four intergenic compressions all match small simple repeat regions (Supplemental Table S8).

Functional annotation identifies new protein features

Several approaches to assess function were applied including InterProScan and I-TASSER among others (Methods). As a result, 138 new C. parvum protein annotations were generated or modified. The percentage of CpIA genes annotated as uncharacterized proteins was reduced from 40% to 33% in all reannotated sequences (Supplemental Table S9). Many new features including domain and repeat content were added to 738 previously uncharacterized proteins. In addition, 729 predicted CpIA CDSs have signal peptides, and 1990 have GO assignments (Supplemental Results). Using I-TASSER protein structure searches, 1414 CDSs were further assessed for confidence, and 1008 predicted structures were assigned as high-confidence by random forest categorization. In CpIRef, 143 previously uncharacterized proteins were assigned high-confidence GO terms.

New transporter genes are identified

We further characterized transporter genes using three different prediction methods. A total of 152 proteins in CpIA and Ch30976 were identified as transporters, including 128 confident candidates and 24 putative candidates (Supplemental Table S10). This represents an increase of 53 transporters relative to the CpIRef GO annotation (CryptoDB v36) (Heiges et al. 2006). Most identifiable transporters are related to purine metabolism, peptidoglycan biosynthesis, oxidative phosphorylation, and N-Glycan biosynthesis pathways (Fig. 3). Six translocases were also identified.

Figure 3.

CpIA assembly and annotation reveal new transporters. The numbers of transporters correspond to the counts of genes encoding each type of transporter protein: (ABC) ATP-binding cassette transporter; (MFS) major facilitator superfamily; (DMT) divalent metal transporter; (AAAP) amino acid/auxin permease; (MC) mitochondrial carrier; (ZIP) zinc transporter protein; (CPA) cation/proton antiporter; (SulP) sulfate transporter; and (PUP) purine permeases.

Open in new tabDownload PowerPointLink to figure

Entire subtelomeric regions are duplicated

As shown in the read depth coverage analysis and in Table 1, Supplemental Table S2, and Supplemental Figure S1, the new CpIA assembly was able to recover ∼2.3 kb cumulative length in collapsed regions relative to CpIRef. One subtelomeric region on Chr 1 in CpIA, previously reported on Chr 5 in CpIRef (but not linked to Chr 5 in the HAPPY map), still shows signs of sequence compression suggesting that most of the genes present in this region have more than one copy (Fig. 4A; Supplemental Fig. S7). This region reveals at least 13 genes that vary in copy number between different Cryptosporidium species (Supplemental Fig. S8). The genes contained in this region are 18S rRNA, 5S rRNA, and 28S rRNA; uncharacterized proteins; a GMP synthase; an aspartate-ammonia ligase; tryptophan synthase beta; and a cluster of several MEDLE genes. Some of these genes, such as the tryptophan synthase beta and the MEDLE's are the focus of considerable research because they may be related to parasite survival and are potentially involved in parasite invasion, respectively (Sateriale and Striepen 2016; Li et al. 2017; Fei et al. 2018). The predicted number of copies of rRNAs and MEDLE's are underrepresented because they also have paralogs on Chr 2 and Chr 5, respectively. The Illumina pileup of ∼1350 reads on Chr 2, positions 681,607 to 686,953 (Supplemental Table S2; Supplemental Fig. S1) is the region where the 18S/28S rRNA gene(s) are located on this chromosome. The five 18S genes are identical, and 28S rRNAs have three gaps (Supplemental Fig. S9). Thus, alignment competition explains why the read coverage varies relative to the equivalent 18S/28S rRNA Chr 1 pileups. Regions with pileups on inner portions of Chr 5, 7, and 8 are low complexity regions composed by tandem repeats (Supplemental Table S8). In Chr 5, we have one uncharacterized protein (CPATCC_0023030), full of tandem repeats and good RNA-seq support for its expression. On Chr 7 and 8, these regions are smaller than 100 bp and do not contain any annotated genes.

Figure 4.

Resolution of repetitive subtelomeric regions found on Chr 1 identifies missing telomeres on Chromosomes 7 and 8. (A) Illumina reads from CpIA are mapped to the CpIA Chr 1 long-read assembly subtelomeric region to identify read pileups and estimates of sequence copy number by normalizing against the average genomics Illumina read depth. Vertical gray areas indicate regions with annotated genes. Annotated genes are represented below the shaded regions, the 5.8S rRNA is present but not indicated. (B) Subtelomeric variaton observed on different CpIA chromosomes is supported by CpBGF ONT long reads. Individual ONT long reads provide evidence of at least four different yet related subtelomeric regions that extend into the chromosomes that were missing telomeres in CpIA (Chr 7 and Chr 8) in addition to Chr 1. The white and black reference bar above each collection of annotated ONT reads identify the resolved subtelomeric regions (white) and linkage to existing assembly (black). The penultimate read on the Chr 7 3′ end panel indicates a unique region of insertion (nucleotide positions 1,191,705–1,217,462). This region contains mostly uncharacterized proteins and two transferases. Each ONT read is annotated as indicated in the key shown in A.

Open in new tabDownload PowerPointLink to figure

Because there is an apparent compression in a subtelomeric region assembly with no gaps and good PacBio long-read coverage, we hypothesized that these extra copies might derive from additional copies of this region. The CpIA assembly was only missing three telomeric regions, both ends of Chr 7 and one telomere of Chr 8. Using existing PacBio long reads we were able to identify a few reads that extended into rRNA regions on the chromosomes missing telomeres. We attempted reassembly with only PacBio reads and we could not resolve the missing regions. Thus, we generated very deep (2260×) ONT single-molecule reads from CpBGF, (ATCC DNA was not available, only 143 SNVs are detected between the strains, of which 108 are indels) (Supplemental Table S11). The ONT reads revealed related, yet unique subtelomeric regions linked to the chromosomes missing their telomeres, in addition to Chr 1 (Fig. 4B). We found good ONT long-read support for these regions (Supplemental Fig. S7). Each distinct subtelomeric region begins with chromosome-specific sequences followed by a conserved ribosomal RNA cluster that is followed by the duplicated subtelomeric region and telomere. There are many ONT and PacBio reads that link the unique chromosomal regions and the beginning of the subtelomeric gene families but only a few span the entire chromosome end. We also note that there is slight variation observed among the reads for each subtelomeric region distal to the rRNA cluster.

New positively selected genes are identified in C. parvum

The new gapless CpIA genome assembly and annotation presented an opportunity to revisit the prediction of genes evolving under positive selection in this species. We performed a single-nucleotide variant (SNV) analysis using 136 different C. parvum WGS data sets obtained from the NCBI GenBank database (https://www.ncbi.nlm.nih.gov/genbank/) (Supplemental Table S12) using the new CpIA assembly and annotation. A total of 24,407 positions were found to contain at least one high-confidence biallelic variant. Multiallelic calls were removed to guard against mixed infections. The biallelic variants reflect 3892 genes, 342 of which show a πN/πS ratio of nonsynonymous/synonymous rates of >1.5 (Supplemental Table S13). Of the 342, 17 genes were previously identified and 145 are classified as uncharacterized proteins, 105 of which are annotated as having a signal peptide or being secreted. All previously identified genes evolving under positive selection were detected, including Insulinase-like protein (CPATCC_0017080), an uncharacterized secreted protein (CPATCC_0010380), gp60 (CPATCC_0012540), and others (Strong et al. 2000; Sanderson et al. 2008; Nader et al. 2019; Zhang et al. 2019). Of the top 10 genes by πN/πS ratio, nine appear to be new to this study. Gene family members such as MEDLEs, FLGN, and SKSR were also detected, but new members of each of these families are identified as also evolving under positive selection. Because the putative new subtelomeric repeats (Fig. 4) were not included in these analyses (they were identified in a different strain), evolution of the MEDLE genes may be an overestimation. A family of WYLE proteins (Sanderson et al. 2008) is also identified as being positively selected.

Sample DNA sources

C. parvum IOWA-ATCC (CpIA) DNA from oocysts/sporozoites was purchased from the ATCC. The source was the University of Arizona, Sterling Parasitology Laboratory. It is a GP60 subtype (IIa) like the current C. parvum IOWA II reference (CpIRef) genome sequence. C. parvum IOWA DNA was also prepared from oocysts obtained in 2018 from Bunch Grass Farms, Deary, Idaho, referred to as CpBGF in this study. C. parvum KSU-1 genomic DNA was also prepared in 2006 from oocysts obtained from Steve Upton. Public sequence data were accessed from the NCBI BioProject database (https://www.ncbi.nlm.nih.gov/bioproject/) under accession numbers PRJNA252787, PRJEB3213, PRJNA388495, and PRJEB10000. Accession numbers for the 136 C. parvum sequences used for evolutionary analysis are detailed in Supplemental Table S12.

C. parvum IOWA-ATCC sequencing and genome assembly

PacBio RSII and Illumina HiSeq 2000 sequencing were performed at the Wellcome Sanger Institute, United Kingdom. The CpIA reads were first assembled using the PacBio open source SMRTlink v6.0 from nine PacBio SMRT cells, with ∼75× mean genome coverage. The resulting assembly was then submitted to the accuracy improver tool Sprai 0.9.9.23 (https://sprai-doc.readthedocs.io/en/latest/index.html), and then gaps were filled using PBJelly 15.24.8 (English et al. 2014) with PacBio reads and IMAGE 2.4.1 (Swain et al. 2012) with Illumina reads. A manual inspection and improvement using GAP5 (Bonfield and Whitwham 2010) was applied, and the final scaffolded genome assembly was polished with Illumina reads using iCORN2 0.95 (Otto et al. 2010) and Pilon 1.22 (Walker et al. 2014).

ONT single-molecule long-read sequencing was performed on DNA from CpBGF (ATCCPRA-67DQ was out of stock) following the recommended R9.4.1 flow cell protocol. MinION ONT sequencing was performed at the Georgia Genomics Bioinformatics Core at the University of Georgia, using an R.9.4 flow cell and the Ligation Sequencing kit (SQK-LSK109). The ONT long reads generated >2000× coverage of the C. parvum genome. This high coverage complemented the PacBio data to confirm and resolve several complex regions (Supplemental Methods). The final assembly was submitted along with CpIRef, Ch30976, and CtUGA55 to QUAST v.5.02 (Gurevich et al. 2013) to compare and evaluate the quality of CpIA. All sequencing statistics are in Supplemental Table S12.

Cryptosporidium genome reannotation

Genome annotation was generated with ab initio prediction using GeneMark-ES 4.57 (Lomsadze et al. 2005); evidence-trained predictions were made using SNAP/MAKER (Cantarel et al. 2008; Johnson et al. 2008) and AUGUSTUS (Stanke and Morgenstern 2005). For training, we used publicly available data from each respective species: RNA-seq, ESTs, previously predicted proteins, and MassSpec proteomics data when available. In parallel we also generated transcriptome assemblies using HISAT2 v.2.1.0 (Kim et al. 2015) and StringTie v.1.3.4 (Pertea et al. 2015), and noncoding RNA predictions were generated for C. parvum as described (Li et al. 2021). Manual curation of all genes in the context of existing molecular evidence was performed using WebApollo2 (Lee et al. 2013).

We performed comparative genome annotation using the Artemis Comparison tool 17.0.1 (Carver et al. 2005). OrthoFinder v.2.3.7 (Emms and Kelly 2015) was used to detect paralogs, orthologs, and singletons. All singletons were then manually compared using MCScanX 0.8 (Wang et al. 2012) and JBrowse (Buels et al. 2016) to verify their uniqueness and assess the contribution of sequence gaps or misassembly to the findings. We considered the following error types: split genes caused by frameshifts or early stop codons, lack of stranded RNA-seq to confirm the gene model, and the presence of a gapped region in the genome assembly. All genes that did not fall into one of these categories were identified as unique.

Functional annotation

Following structural annotation, the predicted protein sequences were used to search the Swiss-Prot, Trembl, and the NCBI nonredundant protein database with BLASTP using an e-value threshold at the superfamily level of 1 × 10⁻⁶. Protein structure similarity was explored using I-TASSER (Roy et al. 2010) as in Ansell et al. (2019) and Supplemental Methods. Blast2GO (Conesa et al. 2005) version 4.1.9 was used to assign Enzyme Code (EC) and Gene Ontology (GO) terms. We compared the existing protein product names to the new functional results. Some structural information, such as protein domain and repeat pattern content were added to some uncharacterized proteins, and nomenclature errors were corrected according to the NCBI guidelines.

Transporter prediction

Predicted proteins were submitted to four different transporter prediction methods: (1) BLASTP against TCDB (Saier et al. 2009) with a threshold e-value of 1 × 10⁻⁵ cutoff; (2) TMHMM (Server v. 2.0) (Krogh et al. 2001) and SignalP (Server 4.1) (Bendtsen et al. 2004) to reduce false positives from the TCDB BLASTP results. Transporter candidates with no transmembrane domains or candidates with only one transmembrane prediction while having signal peptides predicted were removed. (3) TransAAP (Ren et al. 2007), a Transporter Classification tool at TransportDB v2.0 (http://www.membranetransport.org/transportDB2/index.html), was used to provide information about potential transporter identity and substrate; and (4) a structural proof for candidate transporters using Phyre2.0 (Kelley et al. 2015). Final candidate transporters were checked according to the preceding results as well as annotations obtained from InterProScan 5.44 (Jones et al. 2014).

Comparative analyses

Structural variation sites were calculated using NucDiff v2.0.3 (Khelik et al. 2017), and the major inversions observed between CpIRef and CpIA were verified by PCR. Primers to test both ends of each orientation were designed using Primer3 v0.4.0 (Supplemental Table S4; Untergasser et al. 2012). Primers designed to be specific and conserved among the species were tested using an in silico PCR amplification tool (San Millán et al. 2013). These regions contain repeats, so the amplicons range from 2 to 9 kb to avoid them. PrimeSTAR GXL DNA polymerase (TAKARA) was used with Long-PCR conditions: initial 3 min hot start at 98°C, 35 cycles of 10 sec denaturation at 98°C; 15 sec primer annealing at 55.4°C; and 10 min elongation at 68°C; followed by 10 min elongation at 72°C. PCR products were separated in a 0.8% agarose gels and stained with ethidium bromide.

The consistency of annotation and potential gene family CNVs were determined with OrthoFinder v.2.2.7. CNVs were also determined by aligning Illumina sequence reads from each species studied to the new CpIA genome sequence to check for variations in read depth coverage. Alignment was performed using BWA-MEM 0.7.17 (Li and Durbin 2009) with default options, and the alignment depth per base was calculated using BEDTools genomecov 2.29.2 (Quinlan and Hall 2010) and SAMtools depth 1.6 (Li et al. 2009). Read depth coverage plots were generated using the reshape R package (Wickham 2007; R Core Team 2011). To avoid the interference of multiply-mapped regions, only mapped reads were kept for this analysis. For plotting purposes, the Ch30976 genome was scaffolded using the CpIA chromosomal genomic structure using RagTag v.2.0.1 (Alonge et al. 2019).

Resolving the structure of repetitive subtelomeric regions

Following the CNV analysis, the sequence content of the subtelomeric compressed regions and their CpIA assembly noncompressed chromosomal sequence boundaries containing at least 10 genes of Chromosomes 1, 7, and 8 were used to build a BLAST database. We then used this database and BLASTN 2.10.0 (Camacho et al. 2009) to detect CpBGF ONT reads capable of aligning to both subtelomeric and chromosomal boundary regions. The few reads meeting these criteria were evaluated and visualized by aligning all subtelomeric ONT reads to the unique pre-subtelomeric regions of Chromosomes 1, 7, and 8 using the Geneious mapper 2019.1.3 (https://www.geneious.com) with medium sensitivity and minimap2 v.2.22 (Li 2018). Finally, the longest ONT reads were polished with Illumina reads and annotated as previously described for gene content analysis.

Variant analysis, selection prediction, and populational analysis

Illumina sequence reads from 136 different isolates of C. parvum from different geographical locations (Supplemental Table S8) as well as CpBGF were aligned against the CpIA reference genome sequence using BWA-MEM. The BAM files were parsed to select uniquely mapped reads using Picard (http://broadinstitute.github.io/picard/) and then submitted to the GATK 3.8 Haplotypecaller (McKenna et al. 2010). The results were filtered by mapping quality greater than 40 and depth coverage greater than 10. Because mixed infections exist, we restricted analysis to biallelic sites. The individual VCF files were combined into one GVCF file using the GATK tool GenotypeGVCF. After selecting only SNVs from this data, the combined GVCF file was annotated using SnpEff v.4.3 (Cingolani et al. 2012). The SNV variants from the combined annotated GVCF file had their πN/πS ratio (Nei and Gojobori 1986) estimated using SNPGenie 1.0 (Nelson et al. 2015). To avoid noise in the data and identify top candidates, genes with ratios >1.5 were detected and denoted as evolving under positive selection in the C. parvum population analyzed. The higher threshold of 1.5 was chosen based on known genes evolving under positive selection, such as gp60 (Strong et al. 2000) and Insulinase-like (Zhang et al. 2019).