Eimeria zuernii (Eimeriidae: Coccidia): mitochondrial genome and genetic diversity in the Chinese yak

Background Coccidiosis caused by Eimeria zuernii (Eimeriidae: Coccidia) represents a significant economic threat to the bovine industry. Understanding the evolutionary and genetic biology of E. zuernii can assist in new interaction developments for the prevention and control of this protozoosis. Methods We defined the evolutionary and genetic characteristics of E. zuernii by sequencing the complete mitogenome and analyzing the genetic diversity and population structure of 51 isolates collected from eight yak breeding parks in China. Results The 6176-bp mitogenome of E. zuernii was linear and encoded typical mitochondrial contents of apicomplexan parasites, including three protein-coding genes [PCGs; cytochrome c oxidase subunits I and III (cox1 and cox3), and cytochrome b (cytb)], seven fragmented small subunit (SSU) and 12 fragmented large subunit (LSU) rRNAs. Genome-wide comparative and evolutionary analyses showed cytb and cox3 to be the most and least conserved Eimeria PCGs, respectively, and placed E. zuernii more closely related to Eimeria mephitidis than other Eimeria species. Furthermore, cox1-based genetic structure defined 24 haplotypes of E. zuernii with high haplotype diversities and low nucleotide diversities across eight geographic populations, supporting a low genetic structure and rapid evolutionary rate as well as a previous expansion event among E. zuernii populations. Conclusions To our knowledge, this is the first study presenting the phylogeny, genetic diversity, and population structure of the yak E. zuernii, and such information, together with its mitogenomic data, should contribute to a better understanding of the genetic and evolutionary biological studies of apicomplexan parasites in bovines. Graphical Abstract Supplementary Information The online version contains supplementary material available at 10.1186/s13071-023-05925-8.

China has the largest population of the yak in the world, with 16 million yaks accounting for > 95% of the global population [6].Recent epidemiological evidence, however, showed that E. zuernii is commonly prevalent in yak farms and is becoming the leading cause of weight reduction in young individuals across the Qinghai-Tibetan Plateau of China [5,[7][8][9][10].Such situations highlight the significance and necessity of the diagnosis and control of E. zuernii.Current diagnosis and identification of this protozoan parasite typically rely on morphological characteristics.Unfortunately, E. zuernii has similar oocyst morphotypes and overlapping biological features with other coccidian species.In such a context, obtaining an efficient approach to identify E. zuernii infections becomes urgent for clinical diagnosis, epidemiological investigation, and control, and achieving this goal is foreseeable only through the utilization of molecular approaches.Although the internal transcribed spacer 1 (ITS-1) region and small subunit ribosomal DNA (SSU) of the nuclear ribosomal DNA (rDNA)-based PCR amplifications were applied to differentiate between bovine coccidian species, including E. zuernii [11][12][13], cumulative studies demonstrated that the mitogenomic DNA (mtDNA) seems to be more powerful in delimiting individual Eimeria species because of its rapid evolution rate, absence of recombination, and matrilineal inheritance [14].Furthermore, Awadi et al. showed that the Eimeria mitogenomes have the ability to illustrate their evolution and host adaptations [15].Therefore, it is reasonable that the mitochondrial (mt) datasets are widely used for species-specific identification in bovine coccidiosis [16][17][18][19][20]. Unfortunately, until now, no complete information on the mitogenome of E. zuernii has been characterized.Additionally, a fundamental understanding of the genetic structure and diversity of E. zuernii would also assist in the development of a long-lasting control strategy.
Several in-depth studies have proven their significance and effectiveness in controlling protozoan Plasmodium spp.and Toxoplasma spp. in humans as well as Eimeria spp. in chickens by defining their genetic diversity and population architecture [21][22][23].Of course, numerous factors, including the effective population size, distance between populations, host dispersal ability, evolutionary history, host population structures, and the complexity of the life cycle, may affect the genetic makeup of E. zuernii populations [24][25][26].However, little is known about the genetic diversity and population structure of E. zuernii so far.Herein, we decoded the complete mitogenome of E. zuernii isolated from Chinese yaks and further determined its population genetic diversity by sequencing the cox1 gene of 51 E. zuernii isolates that were collected from different geographical yak breeding parks in Sichuan, China.These results would add novel insights to phylogenetic, population genetic, and molecular epidemiological studies on this yak parasite.

Sample collection
In November 2022, as part of ongoing surveys of coccidiosis in yaks, 456 Eimeria oocyst-positive fecal samples were obtained from eight Yak Modern Industrial Parks in Sichuan, China: 43 in Xinlong, 102 in Hongyuan, 25 in Ya'an, 89 in Ruoergai, 32 in Luhuo, 48 in Litang, 31 in Daofu, and 86 in Ganzi (Fig. 1).For species identification, about 500-g feces from each sampling were mashed and suspended with running water, followed by filtration through a series of sieves (300, 150, and 80 mm; Endecotts, London, UK) and culture in a 5% (w/v) potassium dichromate solution at 28 °C for 72 h under forced aeration.Sporulated oocysts were speciated according to morphometric keys [27].Of 456 Eimeria isolates, 51 were preliminarily identified as pure infections with E. zuernii (Additional file 1: Table S1), and their oocysts were further subjected to enrichment with saturated sodium nitrate flotation and molecular identification by PCR amplification using the species-specific primers (forward: 5′-CCC ACT ACA TCC AAC CTC CTG-3′; reverse: 5′-GCG TTC GGA AAT CTG ATG GT-3′) [9] and sequence comparisons with the targeted ITS-1 regions.All isolates were determined to share > 99.2% identities with the query sequence of E. zuernii in GenBank (Additional file 2: Fig. S1; accession nos.OR351547-OR351597 vs. AB769665).

Sequencing and assembly
After species identification and oocyst quantitation using a McMaster chamber among these E. zuernii isolates, one sample from Ganzi was selected for mitogenome sequencing because of its high number of oocysts per gram (OPG = 3.8 × 10 6 ).About 2.5 million E. zuernii oocysts were harvested and purified by using a combination of saturated sodium nitrate flotation and discontinuous sucrose gradient centrifugation.Following vortex breakdown, the genomic DNA was extracted from oocysts using the Genomic DNA Kit (Tiangen, Beijing, China).Based on the quality and quantity assessment, ~ 2.5 µg genomic DNA was fragmented to construct a 250-bp paired-end (PE) library, followed by sequencing on an Illumina HiSeq X-TEN platform (Ber-ryGenomics, Beijing, China).The clean reads (~ 1.8 Gb) were assembled with MITObim v1.9.1 [28].The genome was validated by PCR amplifications using four overlapping fragments (Table 1) and then annotated using MITOS2 (http:// mitos2.bioinf.uni-leipz ig.de/ index.py) and a whole genome-guided gene alignment against Eimeria mephitidis (accession no.KT203398).Based on the E. zuernii mitogenome sequenced here, a pair of primers, cox1-F (5′-TTG GTT GGA CTC TAT ACC CTC-3′) and cox1-R (5′-AGA TAG TAC AAA ATG AAA ATG AGC -3′), were designed to amplify the cox1 gene (776 bp) from the remaining 50 E. zuernii isolates.PCR

Sequences analysis
Nucleotide composition and codon usage of the E. zuernii mitogenome were measured with MEGA X (www.megas oftwa re.net).The base skewnesses of different mitogenomic regions of E. zuernii were determined using the formulas AT skew = (A − T)/(A + T) and GC skew = (G − C)/(G + C) [29].The multi-alignments of nucleotide and amino acid sequences of protein-coding genes (PCGs) obtained from the sequenced mtDNA and those of other available Eimeria parasites (Additional file 3: Table S2) were achieved using MEGA X.Based on the alignments, the similarities and divergences of PCGs were determined using DNAstar v5.02 (DNAStar Inc., Madison, WI), genetic distances were calculated based on Kimura-2-parameter (K2P) [30] with MEGA X, and the ratio of the nonsynonymous substitution (Ka) and synonymous substitution (Ks) of each PCG was calculated with KaKs_Calculator [31].

Phylogenetic and population structure analyses
The phylogenetic position of E. zuernii and relationships between different host-originated eimerian parasites in the genus Eimeria were inferred from the complete mitogenomes of E. zuernii and other congeneric species (Additional file 3: Table S2), using Isospora sp.(accession no.KP658103) as the outgroup.After nucleotide alignments using MAFFT v7.450 [32] and filtration of the ambiguous regions using Gblocks v0.91b (http:// molev ol.cmima.csic.es/ castr esana/ Gbloc ks.html), a 6552-bp multi-sequence dataset was built to reconstruct the phylogenetic tree with three algorithms, including maximum parsimony (MP), maximum likelihood (ML), and Bayesian inference (BI).For the MP analysis, the tree was constructed by PAUP v4.0b10 [33] using equally weighted parsimony and heuristic searches and tree-bisectionreconnection (TBR) branch-swapping, and the optimal topology was obtained using the Kishino-Hasegawa method with 1000 replicates.The ML analysis was implemented using PHYML v3.1 [34] under the optimal evolutionary model "TIM + F + I + G4" that was selected with the "Auto" option on the W-IQ-TREE web server (http:// iqtree.cibiv.univie.ac.at) using an ML + rapid bootstrap algorithm with 1000 replicates.The BI analysis was achieved by MrBayes v3.2.7 [35] using the optimal evolutionary "CAT + GTR + G" model chosen by ModelFinder [36] and four independent Markov chain Monte Carlo (MCMC) chains running for 1,000,000 generations; after sampling a tree every 0.1% generation, a consensus tree was obtained and visualized in Treeview X (https:// www.linux links.com/ treev iewx/).In parallel, 50 E. zuernii isolates sampled from eight different geographical yak breeding parks were subjected to an assessment of population diversity using the cox1-based mt dataset.Diversity indices for these E. zuernii isolates, including the number of haplotypes, haplotype diversity, and nucleotide diversity, were calculated in DnaSP v6.12 (http:// www.ub.edu/ dnasp/).A phylogeny of haplotypes was estimated by MP tree using MEGA X, with the significance of each node estimated using 10,000 bootstrap replicates of the dataset.The neutrality indices, including Fu's Fs, Tajima's D, and the pairwise fixation index (Fst), were used to estimate the size variation among populations using Arlequin v3.5.2.2 (http:// cmpg.unibe.ch/ softw are/ arleq uin35/).The median-joining networks among sequences of populations from eight geographic ranges were drawn using Population Analysis with Reticulate Trees (Pop-ART; http:// popart.otago.ac.nz) to reflect whether the genetic differentiation between populations was associated with geographical isolations.

Results and discussion
Eimeria zuernii mitogenome feature The E. zuernii linear mitogenome was 6176 bp in size and encoded three PCGs, cytb, cox1, and cox3, as well as seven interspersed small subunit (SSU) and twelve interspersed large subunit (LSU) rDNA fragments (Fig. 2).
Likewise, the E. zuernii mitogenome was biased towards AT (64.52-67.37%)with T as the most favored base and G the least favored.Across the mitogenome, there were three overlapping regions located between cytb and cox1 genes (4 bp), between LSUF and LSUG rRNA genes (7 bp), and between LSUA rRNA and cox3 genes (6 bp), respectively.Moreover, 18 intergenic spacers were also observed, ranging in size from 9 to 221 bp, in the E. zuernii mitogenome (Additional file 4: Table S3).

Comparative mitogenomics among Eimeria
To understand the evolutionary divergence of Eimeria spp., the interspecific variations between and among E. zuernii and congeneric species were determined using the nucleotide and amino acid sequences of single and concatenated mt PCGs.The interspecific variations for each PCG ranged from 11.4% to 29.4% at the nucleotide level and 6.4% to 23.5% at the amino acid level.Furthermore, on the basis of the concatenated PCGs, the interspecific variations ranged from 14.0% to 23.2% for nucleotide sequences and 9.0% to 16.3% for amino acid sequences (Additional file 6: Table S5).It was obvious that the highest sequence variation was found in the cox3 gene while the lowest sequence divergence was found in the cytb gene, suggesting cytb was the most conserved PCG among Eimeria mitogenomes, which is consistent with previous findings in species of Babesia, Theileria, and Plasmodium [44,45].Moreover, it also seemed noteworthy that regardless of nucleotide or amino acid levels, E. zuernii always shared the lowest sequence variation with Eimeria mephitidis, to some extent, suggesting their closer genetic similarity than other Eimeria species.

Evolutionary and phylogenetic analyses
To measure the interspecific genetic similarity between E. zuernii and other Eimeria spp., single or concatenated mt PCGs were also used to calculate genetic distances among Eimeria based on the K2P model [30].As shown in Fig. 4a     b Phylogenetic tree inferred by MP, ML, and BI algorithms using complete mitogenome datasets of E. zuernii and other congeneric species.Eimeria species, together with their branches, are marked with different colors according to host origins as follows: bird in light blue, chicken in dark blue, horse in light green, mephitis in dark red, mouse in yellow, mustela in purple, rabbit in green, turkey in blue, and yak/cattle in red.Eimeria zuernii sequenced in this study is indicated in bold font along with a star.The color boxes at each node showed bootstrap values for MP/ML/BI E. maxima.Furthermore, among genetic distance structures, the cox3-based K2P values were all obviously larger than those of the cytb, cox1, and concatenated PCGs, in agreement with the aforementioned result in which cox3 was determined to be the most divergent gene among Eimeria PCGs.In parallel with genetic distance analysis, the selective pressure placed on mt PCGs of Eimeria during evolution was determined by qualification of the Ka, Ks, and Ka/Ks ratios (Additional file 7: Table S6).Our results showed that all Ka/Ks values were < 1, suggesting that three PCGs were subjected to negative or purifying selection, to a certain extent, in agreement on the conservation and validity of the genus Eimeria in the course of evolution.
The phylogenetic relationships among eimerian parasites were inferred based on the complete mitogenomic dataset (containing 4315 conserved and 2237 variable sites) of Eimeria spp.available in GenBank (Additional file 3: Table S2).As shown in Fig. 4b, three phylogenetic trees (MP/ML/BI) consistently supported a paraphyletic relationship among species of Eimeria, and this relationship can be marked according to their host origins, except for bird-infecting species, consistent with previously conducted phylogenetic studies [46][47][48][49].Notably, E. zuernii was placed close to E. mephitidis and together grouped with Eimeria cf.ictidea and Eimeria furonis as a clade (Clade I), with high statistical support (all statistical values ≥ 99 or 0.99).Furthermore, this clade was more closely related to Eimeria species of birds (Clade II) than to Eimeria species of rabbits (Clade III), in accordance with recent mtDNA-based phylogenetic conclusions [18,19,41,50].It appeared that these three clades occupied the main part of the phylogenetic tree, although their topological relationships with Eimeria species of mice and horses remained to be determined because of the poor bootstrap supports observed here (Fig. 4b).The previous mt cox1-based phylogenetic analysis showed that the E. zuernii-contained clade was more closely related to mouse Eimeria spp.than to bird Eimeria spp.; however, the nuclear 18S-based phylogeny indicated a closer relationship between E. zuernii-contained clade and rabbit Eimeria spp.than between E. zuernii-contained clade and bird Eimeria spp.[51].This discordance in the mt and nuclear phylogenies is surprising and might be due to a greater rate of nucleotide change in the mitogenomes of eimerian parasites than that seen in nuclear-encoded sequences, as reported in helminths [52].Therefore, future studies using mitogenomic data from more widespread species or isolates of Eimeria from herbivorous and omnivorous animals worldwide are required to determine the evolutionary relationships within the entire genus Eimeria.

Population structure
The number of haplotypes, haplotype diversity value Hd, and nucleotide diversity π of each E. zuernii population were calculated and are shown in Table 2.The Xinlong population had the highest Hd value (0.867) and π (0.00326), while the Ya'an population had the lowest Hd value (0.500) and π (0.00064).Ganzi and Litang populations had the same Hd value of 0.800, whereas the Litang population had a lower π value than the Ganzi population (0.00163 vs. 0.00232).Likewise, Daofu and Ruoergai populations had the same Hd value of 0.667, whereas the Daofu population had a lower π value than of the Ruoergai population (0.00086 vs. 0.00143).For the Luhuo population, the respective Hd value was 0.750 and π was 0.00138, and for the Hongyuan population, the Hd value was 0.644 and π was 0.00097.Combined, these results showed that these eight E. zuernii populations had high haplotype diversities and low nucleotide diversities.Similar phenomena were seen in other protozoan parasites with large standing population sizes and extremely high fecundities, including Trichomonas vaginalis, Giardia duodenalis, and Theileria annulate [53][54][55], which might reflect a high matrilineal effective population size for E. zuernii and also signify the occurrence of expansion of a low effective population size of E. zuernii after a period In addition, the negative values of Fu's Fs were observed in Hongyuan (Fs = − 0.046), Xinlong (Fs = − 0.071), and Litang (Fs = − 1.350) populations, in contrast with the positive values of Fu's Fs found in Daofu (Fs = 0.172), Ya'an (Fs = 0.172), Luhuo (Fs = 0.330), Ganzi (Fs = 0.469), and Ruoergai (Fs = 0.551) populations (Table 2).Likewise, the negative value of Tajima's D was also observed in the Litang (D = − 0.185) population, in contrast with the positive values of Tajima's D found in Luhuo (D = 1.449),Ganzi (D = 1.573), and Ruoergai (D = 1.754) populations (Table 2).Although the differences among populations in both neutrality tests were not significant, the negative values of the two neutrality indices in the Litang population indicated that the population had undergone expansion, while positive values in the Luhuo, Ganzi, and Ruoergai populations suggested these populations experienced a bottleneck [56][57][58].Furthermore, significantly high Fst values were also observed between populations (Additional file 8: Table S7).The highest Fst value was seen between Hongyuan and Daofu populations (Fst = 0.92157, P < 0.05), followed by Fst values between Hongyuan and Ya'an populations (Fst = 0.91220, P < 0.01), between Hongyuan and Litang populations (Fst = 0.89787, P < 0.05), between Luhuo and Daofu populations (Fst = 0.88665, P < 0.05), between Ruoergai and Daofu populations (Fst = 0.88406, P < 0.05), between Ya'an and Luhuo populations (Fst = 0.88360, P < 0.05), and between Ya'an and Ruoergai populations (Fst = 0.88066, P < 0.01).Such high Fst values indicated a possible geographical isolation of E. zuernii.Indeed, when taking the yak's features, including the free-ranging, high-altitude grazing, and self-reproduction model, as well as the captive history of individuals sampled here into account, it is possible that the lack of translocations and/or introductions of these host populations between parks contributes to this population isolation of E. zuernii.Interestingly, such high Fst values were also observed in other protozoan parasites, such as G. duodenalis and Cryptosporidium spp., which were sampled from self-breeding dairy farms [55,57].To some extent, this high Fst value supports the phenomenon that the low gene flow and high genetic diversity of the bovine protozoa could be accelerated in an independent self-service farming mode.Of course, the phenomenon remains further validated when more additional protozoan population genetic data become available, especially those from hosts bred in an independent self-service farming model.
Based on 51 E. zuernii isolates (including the isolate for mitogenome sequencing in this study) and three reference isolates from GenBank, the mt cox1 sequences containing 31 variable sites (3.99%) defined 24 haplotypes of E. zuernii.The ML-based phylogenetic tree revealed that Ganzi population-specific haplotypes Hap_23 and Hap_24 closely clustered with Hap_1, and then three haplotypes together were paraphyletic with haplotypes of Daofu (Hap_21 and Hap_22), Ya'an (Hap_9 and Hap_10), LuHuo (Hap_14-16), and Xinlong plus Hongyuan (Hap_2-8) (Fig. 5a).It was also clear that Litang population-specific haplotypes Hap_17-20 were more closely related to Ruoergai population-specific haplotypes Hap_11-13 than others.These relationships among these haplotypes were further confirmed using the Network analysis.As shown in Fig. 5b, the network map revealed a radial-shaped clustering (RSC) with Hap_0 as the center, which was hypothetically inferred as the ancestral haplotype.It seemed that the haplotypes in each geographic population were highly specific.Although Hap_1 was shared by Ganzi, Daofu, and Yangzhou populations, Hap_12 was shared by Ruoergai and Litang populations, and Hap_7 was shared by Hongyuan and Xinlong populations.Moreover, E. zuernii from LuHuo and Ya'an appeared to form their own populations, and the latter was inferred as a possibly unsampled or extinct population because of the presence of a median vector (Fig. 5b).In general, network analysis is regarded as a better approach for representing genealogical relationships at a population level than traditional phylogenetic analysis because this approach is able to consider several factors related to intraspecific gene evolution, including the persistence of ancestral haplotypes, the existence of multiple descendant haplotypes, and low levels of sequence variation [58,59].Combined, the phylogenetic and network results suggested a low genetic structure among yak E. zuernii populations although they were prevalent in different geographic ranges of China.

Conclusions
In this study, we sequenced the entire mitogenome of E. zuernii from the yak and characterized its genetic diversity across eight geographic ranges in China.Comparative mitogenomics determined both cytb and cox3 genes as the most and least conserved PCGs, respectively, within the genus Eimeria.Evolutionary and phylogenetic analyses indicated that E. zuernii shared the closest relationship with E. mephitidis, and all species of Eimeria can be grouped into three clades according to their host origins.Further cox1-based genetic structure inference revealed 24 haplotypes of E. zuernii with high haplotype diversities and low nucleotide diversities, suggesting a low genetic structure and rapid evolutionary rate as well as a previous expansion event among E. zuernii populations.Future additional samplings could possibly uncover the geographic structure of E. zuernii in more detail.

Fig. 1
Fig. 1 Map of sampling localities for yak Eimeria isolates included in the present study.The number of isolates infected with E. zuernii (yellow bars) and other Eimeria species (green bars) from each locality are shown, respectively

Fig. 2 Fig. 3
Fig. 2 Linear organization of E. zuernii mtDNA.The order and transcriptional direction of three coding regions for cytb, cox1, and cox3 are indicated by purple arrows and fragments of seven SSU (green) and 12 LSU (light blue) rDNAs are detected between the protein-coding regions and identical to those of other congeneric species.The nomenclature of these ribosomal fragments follows the convention of Feagin et al. (2012) [40].Positive and negative GC-skew and GC content are indicated in dark green, purple, and black, respectively, across the whole genome

Fig. 5
Fig.5 Phylogeny and network map of Eimeria zuernii haplotypes inferred on the basis of 54 cox1 data.a Topology tree of 24 haplotypes of E. zuernii.b Network map of E. zuernii haplotypes.Eimeria zuernii isolates from different geographic origins are shown in dark blue (Hongyuan), blue (Xinlong), light blue (Ya'an), gray (Guelph), pink (North America), dark red (Yangzhou), red (Ganzi), yellow (Daofu), purple (Luhuo), green (Ruoergai), and dark green (Litang).The small black dot represents the median vector, indicating unsampled or extinct haplotypes.The dotted cycle denotes the hypothetically ancestral haplotype.The proportion of haplotype frequencies is shown with sizable coils.The network branches linking the cycles show the relationships between the haplotypes

Table 1
Primer pairs for PCR amplification and their positions in the Eimeria zuernii mitogenome

Table 2
Summary of the genetic diversity of eight geographic populations of E. zuernii on the basis of the cox1 dataset