One hundred nineteen C. hominis Next-Generation Sequencing (NGS) genome projects with shotgun reads were selected and downloaded (Supplementary Table S1) from the Sequence Read Archive-SRA database. Genome projects based only in 454 Technology were excluded. The raw read data were directly download from the SRA database using the fastq-dump tool with the split option activated.

As outgroups, another 15 genomes were used: C. cuniculus (UKCU5: PRJNA492839, UKCU2: PRJNA315496), C. parvum (UKP2:PRJNA253836, UKP3:PRJNA253840, UKP4:PRJNA253843, UKP5:PRJNA253845, UKP6:PRJNA253846, UKP7:PRJNA253847, UKP8:PRJNA253848, UKP14:PRJNA315506, UKP15:PRJNA315507), and C. meleagridis (UKMEL1:PRJNA222838, UKMEL3:PRJNA315502, UKMEL4:PRJNA315503).

C. hominis oocysts were purified from a fecal sample collected from an HIV Colombian female patient by flotation in a saturated sodium chloride solution (Kar et al., 2011). Parasite diagnosis was previously confirmed by Kinyoun stain. Species and subtype identification was done by a nested PCR and sequence analysis of the small-subunit (SSU) rDNA gene and 60 kDa glycoprotein gene (IbA10G2), respectively. Purified oocysts were resuspended in PBS and quantified using a Neubauer chamber slide and light microscopy. The sample was stored at 4°C until DNA was extracted. DNA extraction was performed using the kit NORGEN Stool DNA Isolation (CAT 27600), following the manufacturer’s instructions, with previous freeze-thaw cycles (each for 2 min) in liquid nitrogen and a water bath at 37°C. A solution with approximately 1.85 × 10⁸ oocyst was used for DNA purification. DNA obtained was quantified by light absorption at 260 nm (NanoDropTM, Thermo Scientific) and with PicoGreen^® reagent. The genome was sequenced using Illumina NOVASEQ 6000 instrument at Macrogen (Seoul, Korea). Paired-end reads of 150 bases were generated and deposited at the SRA database under the accession number SRR14522748.

Genomic reads from each experiment were independently assembled using SPADES v3.14.1. First, reads ends were cleaned with rapifilt (homebrew program) filtering at Q30 and only keeping reads with a minimum of 50 bases. Then, assembly was performed using SPADES parameters: -careful -t 40 -m 160 -k 33,55,77,99.

Each assembly was filtered to excluded contaminating sequences using BLASTN (v2.10.1+). Only those longer than 1,000 bases and with a bit score greater than 300 were kept in each C. hominis assembled genome dataset. The bit score was obtained mapping the scaffolds to the C. parvum IOWAII genome reference (CryptoDB v 52), using BLASTN with parameters -evalue 1e-30 -num_alignments 5.

The genomes stats (Total length of sequence, Largest contig, N50 stats, Total number of sequences) were obtained using an in-house python script (see Supplementary Material). The “AltSNPs” (alternate single nucleotide polymorphisms) and “Coverage” of the assembly, were obtained by mapping the cleaned reads of each experiment against the assembled genome, using bowtie2 (version 2.4.1) with default parameters and using samtools (v1.10) view -F 3584 to keep only mapped reads. The coverage was obtained using samtools coverage and obtaining the mean and median of the mean depth of each scaffold. While the SNVs were obtained by counting the number of variants from the Variant Call Format (VCF) created as bcftools mpileup--redo-BAQ--in-BQ 30--per-sample-mF--skip-indels, then bcftools call--multiallelic-caller--variants-only-Ov and then bcftools view -i “%QUAL ≥ 30.”

The SNVs and Identity vs. C. hominis UDEA01 were obtained comparing each genome with the C. hominis UDEA01 using DNAdiff (version 1.3) and getting the TotalSNPs and AvgIdentity from the.report file. The C. hominis genomes that had less than 8,167,000 assembled bases and with more than 3,825 AltSNPs were dropped as they did not meet basic quality parameters, to obtain the selected genomes that are further analyzed. We run BUSCO v5.2.2 (Simão et al., 2015) with parameters “busco-m geno-l coccidia_odb10 -i” for each selected genomes, including outgroups.

The indels vs C. parvum IOWAII were obtained using nucdiff (v2.0.3), with default parameters for each selected genome and the additional C. parvum, C. meleagridis, and C. cuniculus genomes. Using an in-house python script the indels were extracted from “_query_snps.gff” files. The principal component analysis (PCA) analysis was created with python’s sklearn package and plotted with matplotlib.

The selected C. hominis genomes, together with the C. parvum, C. meleagridis and C. cuniculus genomes were used to infer a phylogenetic tree. The sixty one neutrally evolving genes of C. parvum described by Nader et al. (2019) were extracted from each genome. Each gene from all genomes were aligned using MAFFT (version v7.475) with parameters--inputorder--adjustdirection--anysymbol--auto. All the aligned genes were concatenated using catsequences. These aligned sequences were used to infer a maximum likelihood tree using IQTREE2 (2.1.2 COVID-edition for Linux 64-bit) with parameters -B 5000-T AUTO -m MFP+MERGE -rcluster 10 (Lanfear et al., 2012; Lanfear et al., 2014).

Best-fit models selected according to BIC for each coding gene were: K3Pu + F + G4:cgd1_1450 + cgd4_4440 + cgd6_2560 + cgd7_1810 + cgd7_2600 + cgd8_3560,TPM3u + F +G4:cgd1_1730 + cgd1_3650 + cgd1_640 + cgd6_2720 + cgd6_5300 + cgd7_1270 + cgd7_340 + cgd7_890,TIM2 + F + G4:cgd1_2000 + cgd1_3790 + cgd2_2470 + cgd3_4230 + cgd4_2820 + cgd5_2250 + cgd5_4240 + cgd7_2340 + cgd8_2850,TN + F + G4:cgd1_3780 + cgd3_1720 + cgd4_3800 + cgd4_4360 + cgd6_4090 + cgd6_4280 + cgd6_5370 + cgd8_2080 + cgd8_830,TN + F + G4:cgd2_180 + cgd2_3630 + cgd8_1960 + cgd8_3030 + cgd8_5310,HKY + F:cgd2_2060 + cgd5_1340,TN + F + G4:cgd2_3110 + cgd3_1010,K3Pu + F + G4:cgd2_3810 + cgd2_940 + cgd3_380 + cgd4_2620 + cgd5_2890 + cgd7_1330 + cgd7_3550,TPM3 + F + G4:cgd3_2600 + cgd5_3600 + cgd6_2100,HKY + F + I:cgd3_3070 + cgd3_3650 + cgd4_2210 + cgd5_2730,TPM3 + F + G4:cgd3_3310 + cgd5_2700,HKY + F + G4:cgd4_2180 + cgd4_370 + cgd5_1860,TPM2 + F + G4:cgd8_140

We built a haplotype network using the minimum spanning network algorithm (Bandelt et al., 1999) and implemented in PopART (Leigh and Bryant, 2015). We included 95 genomes and excluded four genomes (Chom_EU_SWEH8, Chom_EU_SWEH5, Chom_AFRICA_SWEH3, Chom_AFRICA_SWEH4) because their long branches inferred in the phylogeny, resulting most likely from sequencing errors instead of true genetic variation. The dataset included 118 segregating sites across the 61 genes. The complete step by step analysis pipeline can be found in Supplementary Material.

Statistical analysis and graphics were done in R [R version 4.0.4 64, x86_64-apple-darwin17.0 (64-bit)] and R studio (Version 1.2.1335, Macintosh; Intel Mac OS X 10_16_0). Boxplots were generated in R with the function “boxplot.” Central tendency measures and inter-quartile ranges (IQR) were calculated with R functions “mean,” “median,” “quantile,” and “IQR.”

One hundred and nineteen Cryptosporidium hominis genomes were included in the present study. A new Colombian C. hominis genome with the code UDEAa567 was also generated following C. hominis UdeA01 strategy (Isaza et al., 2015): the oocysts were enriched using a flotation protocol from a stool sample from an HIV + infected individual and then the NGS sequencing was performed using an Illumina Novaseq 6,000 instrument. C. hominis UDEAa567 assembly genome size was 9,052,438 bp. with an N50 value of 93,658 bp. The mean and median read depths were 10.2X and 10.5X, respectively. The UDEAa567 has 99.93% genome identity and 2,655 Single nucleotide variants (SNVs) respect C. hominis UDEA01.

To normalize the analytical conditions for all C. hominis genomes, the processing started with the original raw read data, and then they were all assembled using the same strategy. The C. parvum IOWA II genome (version 52) was downloaded from CryptoDB. All the read sequences used (C. hominis, C. parvum, C. cuniculus, C. meleagridis) are publicly available at the NCBI/SRA website (see Materials and Methods).

The C. hominis genome projects were executed in the last 10 years with different Illumina instruments, different library preparation kits, and were conducted in different laboratories around the world. A small set of these genomes were generated on the ION TORRENT platform. To assess the quality of the assemblies, their general metrics were compared in Figure 1: total assembled bases, median sequencing depth, assembly N50 values, genome nucleotide identity with C. hominis UDEA01, alternate allele count, BUSCO genome completeness, and BUSCO single-copy genes (Supplementary Table S1). The Cryptosporidium assembled genome sizes ranged between 1,148,020 and 10,587,153 Mbp. Despite the broad range of the assembled bases metric, a very narrow dispersion range (an IQR of 21,974 bases) was observed around the median, 9,073,830 bp (Figure 1A). By contrast, the observed median sequencing depth and its IQR were both very variable. The depth ranged from 1.1× to 714×, with a median value of 183× and an IQR of 197 (Figure 1B).

The assembly N50 value also showed a very broad range, from 1,302 to 1,106,561 bp. The median N50 value was 495,148 bases with an IQR of 597,468 (Figure 1C). Furthermore, the global nucleotide identity of all genomes against the C. hominis UDEA01 reference was calculated. They have a median identity of 99.9% and an IQR of 0.02% (Figure 1E). Only 10 genomes had a global identity value below 99.7% (Chom_EU_4127, Chom_AFRICA_SWEH7, Chom_EU_SWEH9, Chom_EU_6940, Chom_EU_6946, Chom_EU_6925, Chom_EU_6934, Chom_EU_4120, Chom_EU_4128, Chom_EU_4118). These genomes are poor-quality references since they are outliers in either assembled genome size or alternative allele count metrics (Supplementary Table S1).

To assess the possible Cryptosporidium coinfections, the reads of each isolate were mapped to its respective assembly and alternate alleles were counted. The observed median number of genomic positions with an alternate allele was 83, with an IQR of 196 (Figure 1D). Eighty-five percent of the C. hominis isolates showed alternate allele counts below 470. Then, to estimate the genetic variations among the C. hominis genomes, we counted the number of SNVs against the C. hominis UDEA01 reference. This analysis showed a median value of 3,324 SNVs with an IQR of 831 (Figure 1F). Ninety percent of the C. hominis genomes had less than 4,000 SNVs regardless its continental origin.

We analyzed the C. hominis genomes with BUSCO using the coccidia ODB10 database. The BUSCO analysis shows the same median and IQR values for both completeness and single-copy metrics, 98.6 and 0.2%, respectively (Figures 1G,H). The duplicated BUSCO median value was 0 with an IQR of 0. All C. hominis genomes, except 1 (Chom_EU_6925), showed duplicated values below 1%. Finally, the BUSCO fragmented metrics showed a median value of 0.4% with an IQR of 0.2% (Supplementary Table S1).

The genomic descriptive measurements showed outliers in several statistic indices, indicating the presence of poor-quality genomes. To filter out them, we set two thresholds: assembled genome size and alternate allele count. The first parameter is directly related to the coverage achieved for each genome and failing to reach a minimum threshold will affect the assembly size and also the accuracy of the assembly. The second parameter, alternate allele count, might indicate the presence of more than one Cryptosporidium genome in the human stool sample. The proposed thresholds rationale are: 1) Assembly size > 8,166,447 bp., which corresponds to 90% of the median value of C. hominis assembled genome size. We think that a representation of 90% of the expected assembly genome size is enough for subsequent comparative or phylogenetic analysis and is equivalent to a BUSCO genome completeness score of >90% (Supplementary Table S1); and 2) AltSNVs threshold of 3,992 SNVs, which corresponds to the percentile 90th of the observed HQ SNVs of the analyzed C. hominis genomes (Figure 2). Applying these two thresholds, 20 genomes were excluded for further analysis, 18 with European origin and 2 with African origin (Chom_EU_SWEH9, Chom_EU_UKH6, Chom_EU_SWEH13, Chom_EU_SWEH11, Chom_EU_H58, Chom_EU_6946, Chom_EU_6945, Chom_EU_6942, Chom_EU_6940, Chom_EU_6939, Chom_EU_6935, Chom_EU_6934, Chom_EU_6925, Chom_EU_6919, Chom_EU_4128, Chom_EU_4127, Chom_EU_4120, Chom_EU_4118, Chom_AFRICA_SWEH7, and Chom_AFRICA_SWEH2). After filtering, the number of genomes per continent was Africa, 33; America, 6; Asia, 29; European Union, 49; Oceania, 2.

Following the previous work published by Nader et al. (2019), 61 genes identified as neutrally evolving in C. parvum and C. hominis were used for the phylogenomic analysis. After filtering out the low-quality C. hominis genomes, ninety-nine genomes were kept for the phylogenetic reconstruction. As outgroups, 15 genomes were used: two C. cuniculus (UKCU5 and UKCU2), 8 C. parvum parvum (UKP2, UKP3, UKP4, UKP5, UKP6, UKP7, and UKP8), three C. parvum anthroponosum (UKP14 and UKP15), and three C. meleagridis (UKMEL1, UKMEL3, and UKMEL4).

A maximum-likelihood tree was constructed encompassing 114 tips. The total matrix of the 61 aligned genes encompasses 140,425 sites. BIC partitions were reduced to 13 (see Materials and Methods). The consensus tree has a Log-likelihood of −266,302.24. The basal branches recreate the previously described topology for the analyzed species, with C. meleagridis defined as the outgroup. The C. parvum subspecies clades C. parvum parvum and C. parvum anthroponosum are reciprocally monophyletic with 100% bootstrap support. The clade C. hominis comprises two well-supported lineages with strong continental bias: One lineage representing isolates from Asia and Africa (100% nodal support), and the other lineage representing genomes isolated from European or American infected humans (97% nodal support). The two isolates from Oceania are part of the European-American clade (Figure 3A). From now on, these two lineages will be referred to as the C. hominis Euro-American lineage (EA) and the Afro-Asian lineage (AA). The consensus tree shows very well-supported basal nodes but low support at the terminal branches due to low divergence among closely related genomes (Figure 3B). It is also noteworthy that the C. hominis EA lineage is populated by Ib genotypes, mostly IbA10G2/IbA12G2, except for C. hominis UdeA01, which has been classified with the genotype Ie.

Four terminal branches showed exceptional longer branch lengths, one in the AA lineage, and three in the EA lineage. It is also particular to see that these four samples come from the same study and were sequenced with the same instrument, the ION TORRENT.

To further confirm the lineage segregation in C. hominis, haplotype analysis was performed based on the 61 neutrally evolving genes. Haplotype network is congruent with the two main lineages (Euro-American and Afro-Asian lineages, Figure 3C). However, we did not find structure within these two lineages, and haplotype diversity appears to be higher in the Afro-Asian lineage.

Indel events are common in Cryptosporidium parasites (Arias-Agudelo et al., 2020). To assess the genomic deletion and insertion profiles within C. hominis species, genomic alignment of the assemblies against the C. parvum IOWA II reference was performed and analyzed. Insertion and deletion events were treated independently to give more clarity.

Deletion events were more common than insertions and spanned from 1 to 78 bases. The most frequent loss harbored a single base and was followed in frequency by the loss of three bases (Figure 4A). All the C. hominis genomes displayed a set of 3,958 shared deletion events. Furthermore, there are 609 specific deletion events in C. hominis EA lineage, and 2,338 specific deletion events in AA lineages.

Like the deletions, the insertion events spanned from 1 to 68 bases, being, in order, more frequent those that involve one and three bases (Figure 4B). Most of the insertion events were shared between both lineages (n = 3,272). Nonetheless, each lineage displayed specific insertion events. The AA lineage has 1,996 common insertions, while the EA lineage has 380.

Indel events may help to understand the relationship among the different parasite clades. With the aim to test the informative signal of the indel event profile of each genome, a Principal component analysis (PCA) analysis of the studied genomes was performed. As can be seen in Figures 4C,D, deletions or insertion event profiles segregate the parasites according to species and the lineage of C. hominis.