All 111 specimens were successfully subtyped and 1–11 sequence types were found at the 12 gene loci. Among them, GP60 and DZ-HRGP had relatively higher sequence polymorphism, with 11 and 7 subtypes being identified, respectively. In contrast, all 111 specimens had the same sequence at the normally polymorphic RPGR locus. Altogether, these C. parvum IId specimens formed 25 MLST subtypes (Table S1).

Sequence data of all 12 loci were concatenated to make a multilocus gene of 7209 bp in length. The genetic diversity of sequences was analyzed using DnaSP 5.10.01. They had 135 polymorphic (segregating) sites, and 12 haplotypes with a haplotype diversity Hd = 0.698 ± 0.039, nucleotide diversity Pi = 0.00692, and average number of nucleotide differences k = 47.296. The genetic diversity of Swedish isolates was higher than that of Chinese isolates (Table S2). Gene flow and genetic differentiation analysis showed significant differences among China, Egypt, and Sweden populations (P < 0.001), with the pairwise fixation index (F_ST) of 0.811 to 0.919. Thus, the value of Nm representing gene flow was only 0.05 by both N_ST and F_ST estimates based on sequence data.

To estimate genetic diversity level of C. parvum IId subtypes, the GP60 locus representing the most polymorphic sequence identified thus far was used in this study. 43 GP60 sequences belonging to 15 subtypes were obtained from the GenBank database, including IIdA14G1, IIdA15G1, IIdA16G1, IIdA17G1, IIdA18G1, IIdA19G1, IIdA20G1, IIdA21G1, IIdA22G1, IIdA23G1, IIdA24G1, IIdA25G1, IIdA26G1, IIdA28G1, and IIdA29G1. DNA polymorphism was analyzed using DnaSP 5.10.01. The nucleotide diversity (Pi) of GP60 sequences was 0.00267 in Western Asia (Iran, Kuwait, and Jordan), which was higher than that in other areas, including European countries, other Asian countries (China, India, and Malaysia), Egypt, and Australia (Table S3).

Neighbor-joining analysis of the concatenated sequences grouped the 25 MLST subtypes into three distinct lineages (Fig. 1). The lineages occupy strikingly distinct geographic distributions; lineage 3 occurred almost exclusively in China, with only Egyptian isolates (MLST-6) being included. In contrast, the Sweden isolates consisted of lineages 1 and 2. Likewise, the phylogenetic network analysis also produced the same result (Fig. S1).

To examine LD within and between loci, we analyzed pairwise associations of polymorphism across the concatenated sequences. LD was estimated using the ZnS statistic, which calculates the average correlation among alleles in all pairs of polymorphic sites. ZnS was elevated not only for all isolates from the three countries but also for isolates within China or Sweden (Table 1). These high ZnS values indicated the C. parvum IId subtypes had a clonal population structure, suggesting the absence of recombination among lineages, which was supported by the D′ statistic (Fig. S2). In addition, the standardized index of association (I^S_A) analysis also confirmed this conclusion (Table S4).

A Bayesian statistical approach was used to infer population substructure in allelic variation using the software STRUCTURE. The 111 C. parvum IId isolates derived from the admixture of three ancestral lineages (K = 3) (Fig. S3). Cryptosporidium parvum IId isolates from cattle, hamsters, and horses in China and cattle in Egypt consisted of ancestral lineage III; whereas the 43 C. parvum IId isolates from cattle in Sweden formed ancestral lineages I and II (Fig. 2). Generally, Chinese isolates and Swedish isolates had different ancestral lineages. In contrast, the isolates from Egypt mainly showed the admixture of lineages I and III (Fig. 2).

A total of 111 C. parvum IId specimens were used in this study including Chinese isolates from dairy cattle, horses, and hamsters and Swedish and Egyptian isolates from dairy cattle (Table S1). Most of the specimens were identified in previous studies41628, whereas the 10 bovine isolates from Egypt and six horse and five bovine isolates from China were diagnosed as C. parvum IId subtypes by DNA sequence analysis of the GP60 gene37.

Twelve microsatellite, minisatellite, and SNP loci were used in the present study. Primers and amplification conditions used in nested-PCR analysis of MSC6-7 (serine repeat antigen), RPGR (retinitis pigmentosa GTPase regulator), MSC6-5 (hypothetical trans-membrane protein), DZ-HRGP (hydroxyproline-rich glycoprotein), ZPT (zinc protease telomerase), CP56 (56 kDa trans-membrane protein), Chom3T (T-rich gene fragment), CP47 (47 kDa protein), HSP70 (70 kDa heat shock protein), Mucin1 (mucin-like protein), ML2 (possible sporozoite cystein-protein), and GP60 genes were previously described1. KOD-FX-Neo amplification enzyme (Toyobo Co. Ltd, Osaka, Japan) was used in PCR. 400 ng/μl of non-acetylated bovine serum albumin (Solarbio Co. Ltd, Beijing, China) was used in primary PCR to neutralize PCR inhibitors. The secondary PCR products were examined by agarose gel electrophoresis and visualized after GelRed™ (Biotium Inc., Hayward, CA) staining. The secondary PCR products were sequenced on an ABI 3730 DNA Analyzer (Applied Biosystems, Foster City, USA), using the secondary primers and the Big Dye Terminator v3.1 Cycle Sequencing kit (Applied Biosystems). The sequence accuracy was confirmed by bi-directional sequencing and by sequencing a new PCR product if necessary.

Sequence alignment was made using the program ClustalX 1.83 (http://www.clustal.org/). Neighbor joining trees were constructed using the program Phylip 3.69, based on the evolutionary distances calculated by Kimura-2-parameter model. The reliability of these trees was assessed using the bootstrap analysis with 1,000 replicates. The phylogenetic network analysis was done using the median-joining algorithm as implemented in Network 4.611 based on the polymorphism of multilocus sequence data from 111 C. parvum IId isolates38. DNA polymorphism, pairwise F_ST values for the different populations, and gene flow were calculated using DnaSP 5.10.0139.

Sequences from 12 loci were combined in a single contig and analyzed for linkage disequilibrium (LD) across the entire composite sequence. LD analysis was conducted using software DnaSP 5.10.0139. D′ was calculated for all pairs of sites. Both the two-tailed Fisher's exact test and the χ² test were used to determine significance of the associations between polymorphic sites. The average LD was estimated by using the ZnS statistic, which averages LD over all pairwise comparisons for S polymorphisms in N sequences. Alternatively, as a confirmation for the ZnS statistic, the standardized index of association (I^S_A) by Habould and Hudson (2000) was also used to assess the LD across all loci. The index and its probability under a null model of complete panmixia were calculated using LIAN version 3.540 with hypothesis testing by a parametric method.

The subpopulations of C. parvum IId subtypes were determined based on the allelic data using STRUCTURE 2.3.3 by K-means partitional clustering and the admixture model41. Twenty simulation runs were conducted for each of K = 2 to 10 using a length of burn-in of 10⁴ and 10⁵ replicates of Markov chain Monte Carlo simulation. The average membership coefficients for the 20 simulation runs of a given K value were generated by CLUMPP v1.1.2 (windows version)42 and a graphical presentation of the average membership coefficients for each isolate was generated in Microsoft Excel. The most appropriate number of K was calculated using an ad hoc statistic-based approach implemented in software Structure Harvester v0.6.93, as described previously43.

Representative nucleotide sequences were deposited in the GenBank under accession numbers KC823118 to KC823157.