Introduction
Theileria annulata, a haemoprotozoan parasite, causes bovine tropical theileriosis (BTT), a tick-borne disease in cattle that leads to leukoproliferation (Salim et al., Reference Salim, Chatanga, Jannot, Mossaad, Nakao and Weitzman2019). The disease is a critical disease for livestock in Asia and North Africa and is transmitted through Hyalomma ticks (Ali, Reference Ali2022). T. annulata are tick-transmitted blood parasites that hinder the livestock sector and cause major economic damage in tropical and subtropical areas globally (Poklepovich et al., Reference Poklepovich, Mesplet, Gallenti, Florin-Christensen and Schnittger2023). Only in India an economic loss of $800 million due to infection caused by T. annulata has been reported (Dandasena et al., Reference Dandasena, Bhandari, Sreenivasamurthy, Murthy, Roy, Bhanot, Arora, Singh and Sharma2018). Economic assessments have estimated that T. annulata infection resulted in financial losses of around 598 133 USD within just 2 years in endemic regions of Turkey (Inci et al., Reference Inci, Ica, Yildirim, Vatansever, Cakmak, Albasan, Cam, Atasever, Sariozkan and Duzlu2007). In China, the commonest theilerioses are caused by T. annulata and T. orientalis. Research carried out over many years has shown that the 2 protozoans are distributed mainly in Northern China, but they also occurred in Southern China. T. annulata is the most troublesome, affecting large numbers of cattle, especially those of exotic breeds (Luo and Lu, Reference Luo and Lu1997). Furthermore, in China, bovine theileriosis has been frequently reported in major cattle-producing regions, including the northwest, northeast and northern areas. Despite annual vaccination campaigns and other control measures implemented by veterinary authorities, infection rates remain high, indicating persistent endemicity and suboptimal efficacy of current interventions. The disease leads to substantial economic losses by causing fever, anaemia, reduced productivity and in severe cases, death in affected animals. Infected cattle often exhibit impaired growth and milk yield, along with increased treatment costs and herd management burdens. The Food and Agriculture Organization (FAO) of the United Nations has estimated that tick-borne theileriosis results in approximately 7 billion USD in economic losses globally each year. Given China’s large cattle population and the widespread presence of tick vectors, T. annulata infection represents a major barrier to sustainable livestock development (Zhao et al., Reference Zhao, Xie, Du, Guo, Li and Guo2017a).
T. annulata is a protozoan parasite with a complex life cycle (Elati, Reference Elati2024). The parasite’s life cycle contains a sexual phase in tick vectors and asexual developmental stages in mammalian hosts. After invading host leukocytes, the parasites reprogram many of the host’s signalling pathways and induce uncontrolled proliferation via poorly understood molecular mechanisms. Massive invasion of red bloods cell by T. annulata merozoites causes anaemia that aggravates pathology of the disease (Liu et al., Reference Liu, Guan and Yin2022). Ticks become infected by ingesting these parasitized red blood cells (RBCs) during a subsequent blood meal (Mehlhorn, Reference Mehlhorn1984; Elati, Reference Elati2024; Liu). Within bovine leukocytes, T. annulata establishes itself in the cytoplasm as a multinucleated schizont, the pathogenic form of the parasite (Woods et al., Reference Woods, Perry, Brühlmann and Olias2021). By hijacking host signalling networks, the schizont promotes continuous host cell division and suppresses programmed cell death (Heussler et al., Reference Heussler, Sturm and Langsley2006; Woods et al., Reference Woods, Perry, Brühlmann and Olias2021). Through a process known as merogony, the schizont differentiates into merozoites inside the leukocyte, and merozoites are subsequently liberated upon RBC lysis (Schmuckli-Maurer et al., Reference Schmuckli-Maurer, Shiels and Dobbelaere2008).
To gain deeper insights into how T. annulata modulates host cellular functions across the schizont and piroplasm stages, we integrated transcriptomic and metabolomic datasets to identify stage-specific host responses and metabolic rewiring associated with infection, generating testable hypotheses about processes that may support parasite survival and pathogenicity. Over the past decade, significant efforts have been made to elucidate host–pathogen interactions during T. annulata infection (Durrani et al., Reference Durrani, Weir, Pillai, Kinnaird and Shiels2012; Kühni-Boghenbor et al., Reference Kühni-Boghenbor, Ma, Lemgruber, Cyrklaff, Frischknecht, Gaschen, Stoffel and Baumgartner2012; Haidar et al., Reference Haidar, Echebli, Ding, Kamau and Langsley2015; Zhao et al., Reference Zhao, Guan, Liu, Liu, Li, Yin and Luo2017b). Studies have shown that schizont-stage parasites reprogram host cell signalling pathways to evade apoptosis and promote uncontrolled proliferation, effectively transforming bovine leukocytes into immortalized cells (Dobbelaere and McKeever, Reference Dobbelaere and McKeever2002; Dobbelaere and Küenzi, Reference Dobbelaere and Küenzi2004; Durrani et al., Reference Durrani, Weir, Pillai, Kinnaird and Shiels2012; Woods et al., Reference Woods, Perry, Brühlmann and Olias2021; Tajeri and Langsley, Reference Tajeri and Langsley2025). These effects are driven by complex host transcriptional responses, including the upregulation of anti-apoptotic genes, cell cycle regulators and inflammatory mediators (Dobbelaere and Küenzi, Reference Dobbelaere and Küenzi2004; Kinnaird et al., Reference Kinnaird, Weir, Durrani, Pillai, Baird and Shiels2013; Ahlawat et al., Reference Ahlawat, Choudhary, Arora, Kumar, Kaur and Chhabra2023; Aktas et al., Reference Aktas, Eren, Kucukler, Eroglu, Ilgun, Yanar and Aydin2023). However, most of these findings are based largely on transcriptomic data, which provides only a partial view of the cellular response (Metheni et al., Reference Metheni, Echebli, Chaussepied, Ransy, Chéreau, Jensen, Glass, Batteux and Langsley2014, Reference Metheni, Lombès, Bouillaud, Batteux and Langsley2015; Haidar et al., Reference Haidar, Metheni, Batteux and Langsley2019). Metabolomics, which captures downstream functional outputs of gene and protein activity, has been investigated in T. annulata (Zhao et al., Reference Zhao, Li, Liu, Guan and Dan2022); however, stage-resolved, host-focused metabolomic profiling and its integration with transcriptomics remain limited, which our study addresses (Metheni et al., Reference Metheni, Echebli, Chaussepied, Ransy, Chéreau, Jensen, Glass, Batteux and Langsley2014). The metabolic alterations in infected bovine erythrocytes and lymphocytes, and their potential contribution to parasite survival and immune evasion, are poorly understood (Zhao et al., Reference Zhao, Li, Liu, Guan and Dan2022). This knowledge gap is particularly pronounced for the erythrocytic stage, where piroplasms reside in RBCs – cells that lack transcriptional activity – posing a challenge for traditional gene expression profiling of host cells (Sae-Lee et al., Reference Sae-Lee, McCafferty, Verbeke, Havugimana, Papoulas, McWhite, Houser, Vanuytsel, Murphy, Drew, Emili, Taylor and Marcotte2022; Tajeri et al., Reference Tajeri, De Laté, Hemmink, Vrettou, Langsley and Morrison2025b). Furthermore, previous studies have largely focused on either the schizont stage in isolation, without a comprehensive comparison of host responses across different stages of infection. This fragmented approach limits our ability to identify stage-specific signatures and molecular pathways that could be exploited for targeted interventions.
In this study, we analysed 3 contrasts: (1) purified parasite stages – schizonts vs piroplasms; (2) infected vs uninfected bovine peripheral blood mononuclear cells (PBMCs) and (3) infected vs uninfected bovine erythrocytes. These contrasts enable separation of host-derived responses from parasite-stage signatures, delineating stage-associated transcriptional and metabolic patterns across host cells and purified parasite stages. Moreover, through the integration of transcriptomic data, which captures changes in gene expression, with metabolomic data, representing downstream metabolic alterations, we identified stage-associated interactions and candidate regulatory pathways that may be engaged during infection. This integrated approach highlights candidate molecular targets and putative metabolic vulnerabilities for future experimental validation and potential diagnostic or therapeutic exploration. Overall, our findings provide a comprehensive resource of stage-associated host and parasite profiles in BTT and illustrate the utility of combining transcriptomics with metabolomics to characterize parasite–host interactions.
Materials and methods
Sample collection
Blood samples for the multi-omics experiment were collected from 6-month-old cattle experimentally infected with T. annulata (Xinjiang Kashi strain; TaXJS) at 3–4 weeks post-infection, and from 1 age-matched uninfected cattle from the same herd, which served as a control.
The experimental infection was conducted via intravenous injection of bovine blood containing piroplasms of the TaXJS strain, which had been cryopreserved in the Vectors and Vector-borne Diseases Laboratory at the Lanzhou Veterinary Research Institute (Lanzhou, Gansu) (Ma et al., Reference Ma, Liu, Li, Xiang, Wang, Liu, Li, Yin, Guan and Luo2020). The TaXJS strain was originally isolated from naturally infected cattle and is characterized as a naturally occurring buparvaquone-resistant strain (Ma et al., Reference Ma, Liu, Li, Xiang, Wang, Liu, Li, Yin, Guan and Luo2020; Zhang et al., Reference Zhang, Zhao and Cao2022). The T. annulata (TaXJS) parasites were obtained from the same laboratory.
Peripheral blood samples were collected from each animal and processed immediately to isolate erythrocytes (RBCs) and PBMCs. For the infected animal and the uninfected control animals, parasite-infected erythrocytes were isolated directly from the blood collected from infected cattle, while erythrocytes from the uninfected cattle served as controls. Isolation of erythrocytes and PBMCs was carried out, followed by repeated washing with phosphate-buffered saline (PBS) to remove serum proteins and contaminants. Purified cell pellets from the infected animal and the uninfected control animal were aliquoted to generate 3 technical replicates per compartment for each condition (i.e. n = 3 technical replicates per group, as shown in the principal component analysis [PCA] plots) and stored at −80 °C until further transcriptomic and metabolomic profiling. Matched aliquots from the same blood draw were used for RNA-seq and metabolomic profiling, ensuring that the 2 omics layers are animal-matched.
Procurement of schizonts stage of T. annulata
T. annulata (Xinjiang Kashi strain, TaXJS)-infected bovine immune cells (leukocytes, B cells and dendritic cells) (Rashid et al., Reference Rashid, Guan, Luo, Zhao, Wang, Rashid, Hassan, Mukhtar, Liu and Yin2019; Liu et al., Reference Liu, Rashid, Wang, Liu, Guan, Li, He, Yin and Luo2020). TaXJS were cultured in RPMI 1640 medium (WILBER, China) supplemented with 10% FBS (WILBER, China). Schizonts were purified, as previously described by Wiens (Wiens et al., Reference Wiens, Xia, Von Schubert, Wastling, Dobbelaere, Heussler and Woods2014). Briefly, TaXJS cells were incubated for 16 h with nocodazole (Sigma-Aldrich, USA) to depolymerize microtubules. The cells were then treated with trypsin (WILBER, China)-activated aerolysin on ice. After removing excess aerolysin, cells were exposed to a temperature of 37 ℃ to stimulate toxin-mediated permeabilization of the host cell plasma membrane. Permeabilization was monitored using Trypan blue exclusion. Schizonts were separated from host cell debris using Percoll gradient centrifugation. Three biological replicates were prepared and extracted in parallel to minimize handling-related variation.
Purification of piroplasms stage of T. annulata
Based on previous relevant literature (Rodriguez et al., Reference Rodriguez, Bueninc, Vega and Carson1986; Sugimoto et al., Reference Sugimoto, Sato, Kawazu, Kamio and Fujisaki1991; Blackman, Reference Blackman1994), purified piroplasms of T. annulata were obtained from infected bovine erythrocytes. When the parasitaemia of T. annulata reached approximately 15%, blood was collected from the jugular vein of cattle into anticoagulant tubes and continuously agitated to prevent coagulation. The blood was then aliquoted into 50 mL centrifuge tubes and centrifuged. The supernatant was discarded, and the pellet was washed with PBS. After another centrifugation, the white cell layer was removed, and the remaining cells were resuspended in 2 times volume of PBS. This suspension was then passed through a leukocyte filter. The filtered RBC suspension was again aliquoted into 50 mL centrifuge tubes and centrifuged, discarding the supernatant. The sedimented RBCs were resuspended in 3 times the volume of 7% glycerol (prepared in PBS), mixed thoroughly and left at room temperature. After centrifugation, the supernatant was discarded, and the cell pellet was transferred to an Erlenmeyer flask. Twenty times the volume of PBS was added to the flask, which was rapidly agitated to lyse the cells. The lysate was aliquoted into 50 mL centrifuge tubes and centrifuged again. The supernatant was discarded, and the pellets were combined into a 50 mL centrifuge tube. The pellet was washed several times with PBS until the supernatant was clear. Finally, the purified piroplasms were aliquoted into 1.5 mL centrifuge tubes for further use. Three biological replicates were prepared and extracted in parallel to minimize handling-related variation.
Transcriptomic sequencing and analysis
Total RNA was extracted from each samples using Trizol reagent (Invitrogen, USA) according to the manufacturer’s protocol. RNA integrity was evaluated using the Bioanalyzer 2100 system (Agilent Technologies, CA, USA). Only RNA samples with sufficient integrity and purity were used for subsequent library construction. Messenger RNA (mRNA) was purified from total RNA using poly-T oligo-attached magnetic beads. After fragmentation, the first strand of cDNA was synthesized using random hexamer primers, followed by second-strand synthesis using dUTP to ensure strand specificity (Trapnell et al., Reference Trapnell, Williams, Pertea, Mortazavi, Kwan, Van Baren, Salzberg, Wold and Pachter2010; Bray et al., Reference Bray, Pimentel, Melsted and Pachter2016). The resulting cDNA was subjected to end repair, A-tailing, adapter ligation, size selection, USER enzyme digestion, amplification and purification (Parkhomchuk et al., Reference Parkhomchuk, Borodina, Amstislavskiy, Banaru, Hallen, Krobitsch, Lehrach and Soldatov2009; Pertea et al., Reference Pertea, Pertea, Antonescu, Chang, Mendell and Salzberg2015). Library quality was assessed using Qubit and real-time PCR for concentration measurement and a Bioanalyzer for fragment size distribution (Goldstein et al., Reference Goldstein, Cao, Pau, Lawrence, Wu, Seshagiri and Gentleman2016). Qualified libraries were pooled in appropriate proportions according to the effective concentration and desired sequencing depth, and then subjected to high-throughput sequencing using the Illumina NovaSeq 6000 platform (PE150) with the NovaSeq Reagent Kit. Sequencing followed the ‘sequencing-by-synthesis’ (SBS) principle, and the resulting base calls were processed into raw reads. Raw sequencing data in fastq format were processed using fastp to remove adapter-containing reads, reads with poly-N and low-quality reads (Garber et al., Reference Garber, Grabherr, Guttman and Trapnell2011; Liao et al., Reference Liao, Smyth and Shi2014; Patro et al., Reference Patro, Mount and Kingsford2014; Bray et al., Reference Bray, Pimentel, Melsted and Pachter2016). Clean reads were aligned to the Bos taurus reference genome ARS-UCD1.2 (NCBI Assembly ID: GCF_002263795.1) and to the Theileria annulata reference genome of the Ankara strain (NCBI Assembly ID: GCF_000003225.3) using HISAT2 (version 2.2.1) in a reference-based manner. For host-focused differential expression analyses, only RNA-seq libraries from PBMC and TaXJS leukocyte samples were used. For TaXJSM-infected erythrocyte libraries, reads mapping to the T. annulata genome were retained to define the piroplasm (parasite) transcriptome, whereas the very low number of reads mapping to the bovine genome were considered residual erythrocyte RNA and were not included in host differential expression analyses. Transcripts were assembled using StringTie (v1.3.3b) in a reference-based manner. Gene expression levels were quantified using featureCounts (v1.5.0-p3), and normalized as FPKM (Fragments Per Kilobase of transcript per Million mapped reads). Differential gene expression analysis was conducted using the DESeq2 R package (v1.20.0) (Mortazavi et al., Reference Mortazavi, Williams, McCue, Schaeffer and Wold2008; Love et al., Reference Love, Huber and Anders2014; Shen et al., Reference Shen, Park, Lu, Lin, Henry, Wu, Zhou and Xing2014), with adjusted P-value ≤ 0.05 and |log2FoldChange| ≥ 1 set as the significance thresholds. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses of differentially expressed genes (DEGs) were performed using the clusterProfiler R package (Kanehisa and Goto, Reference Kanehisa and Goto2000; Anders and Huber, Reference Anders and Huber2010; Robinson et al., Reference Robinson, McCarthy and Smyth2010; Young et al., Reference Young, Wakefield, Smyth and Oshlack2010). Terms with adjusted P-values less than 0.05 were considered significantly enriched.
Metabolomic profiling and analysis
Metabolomic profiling was performed on erythrocytes and leukocyte samples from T. annulata-infected and uninfected cattle, as well as on purified piroplasms and schizonts of T. annulata. For the metabolomic experiments, 3 pairwise contrasts were analysed: schizont-enriched vs piroplasm-enriched preparations (n = 3 vs 3), PBMC vs TaXJS leukocytes (n = 3 vs 3) and RBC vs TaXJSM erythrocytes (n = 4 vs 3). Samples were placed in Eppendorf tubes, resuspended in prechilled 80% methanol by vortexing, briefly melted on ice for 30 s, sonicated for 6 min and centrifuged at 5000 rpm at 4 °C for 1 min. The resulting supernatants were freeze-dried and reconstituted in 10% methanol. The final solutions were injected into the LC-MS/MS system for further metabolomic analysis (Want et al., Reference Want, O’Maille, Smith, Brandon, Uritboonthai, Qin, Trauger and Siuzdak2006, Reference Want, Masson, Michopoulos, Wilson, Theodoridis, Plumb, Shockcor, Loftus, Holmes and Nicholson2013; Sellick et al., Reference Sellick, Hansen, Stephens, Goodacre and Dickson2011; Yuan et al., Reference Yuan, Breitkopf, Yang and Asara2012). Metabolite separation and detection were performed using a Vanquish UHPLC system (Thermo Fisher Scientific, Germany) coupled with an Orbitrap Q Exactive™ HF-X mass spectrometer (Thermo Fisher Scientific, Germany). Samples were loaded onto a Hypersil GOLD column (100 × 2.1 mm, 1.9 µm). The mobile phases consisted of solvent A (water with 0.1% formic acid) and solvent B (acetonitrile with 0.1% formic acid). Metabolites were eluted using a 12 min linear gradient from 2% to 98% B at a flow rate of 0.2 mL min−1, followed by re-equilibration to the initial conditions. The mass spectrometer operated in both positive and negative ionization modes, with a spray voltage of 3.5 kV, capillary temperature of 320 °C, sheath gas flow rate of 35 psi, auxiliary gas flow rate of 10 L min−1, S-lens RF level of 60 and auxiliary gas heater temperature of 350 °C. The raw data generated by UHPLC-MS/MS were processed using Compound Discoverer 3.3 (Thermo Fisher Scientific) for peak alignment, peak picking and quantification (Barri and Dragsted, Reference Barri and Dragsted2013). To monitor analytical reproducibility and instrument stability, a pooled quality control (QC) sample was prepared by mixing equal aliquots from all study samples. QC injections were placed at the beginning of each batch and interspersed with study samples at regular intervals to condition the LC column and assess signal stability. Instrument performance was evaluated by calculating Pearson correlation coefficients between QC injections based on log-transformed peak intensities; pairwise QC correlations were consistently high (all r values > 0.95), indicating good stability of the LC-MS system across the run. Study samples belonging to the same biological contrast were acquired within the same analytical batch, and within each batch, replicates from the same group were injected consecutively rather than in a fully randomized order. Metabolite annotation was conducted by querying mzCloud, KEGG, Human Metabolome Database (HMDB) and LIPIDMaps. For each detected feature, an annotation confidence level (Levels 1–3) was assigned based on the information available (accurate mass, isotope pattern, MS/MS spectral match and, where available, retention time match to reference standards), following the reporting recommendations of the Metabolomics Standards Initiative. In brief, Level 1 annotations correspond to metabolites confirmed with authentic standards (matched m/z, MS/MS spectrum and retention time), Level 2 annotations correspond to putatively identified compounds based on MS/MS and/or database spectral matches and Level 3 annotations correspond to putatively characterized compound classes or tentative metabolite assignments. The annotation level for each metabolite is reported in the supplementary metabolite intensity tables (Supplementary data 2). To minimize biological misinterpretation, compounds whose tentative annotations correspond to pharmaceuticals, cosmetic ingredients or other xenobiotics (for example minoxidil, buspirone or aleuritic acid) were treated as low-confidence (Level 3) matches that may reflect database-driven spectral similarity or exogenous contamination. PCA and partial least squares discriminant analysis were performed using the metaX R package (Wen et al., Reference Wen, Mei, Zeng and Liu2017). To control technical variability, metabolites with a coefficient of variation (CV) > 30% across pooled QC injections were removed prior to downstream statistical analysis. Univariate statistical analysis (Student’s t-test) was applied to evaluate group differences at the metabolite level. At this stage, P-values were treated as nominal and were not adjusted for multiple testing. Differentially abundant metabolites (DAMs) were defined using a combined set of criteria (VIP > 1, fold change (FC) > 2 or FC < 0.5, i.e. |log2FC| > 1, and nominal P-value < 0.05). Functional classification and metabolic pathway annotation of DAMs were performed using the KEGG database. Pathway enrichment analysis was conducted based on the hypergeometric test, and pathways with P-values < 0.05 were considered significantly enriched. The full processed metabolite intensity matrices for positive and negative ionization modes, including feature identifiers, m/z values, retention times, ionization mode, annotation details and confidence levels (Levels 1–3), are provided as supplementary data.
Integrated multi-omics analysis
To explore the correlation between transcriptomic and metabolomic alterations during T. annulata infection, an integrated multi-omics analysis was performed. First, significantly DEGs and DAMs were identified based on the thresholds described above. Genes and metabolites were mapped to the KEGG database to identify shared or associated biological pathways (He et al., Reference He, Zhao, Lu, Wang, Liu, Zeng and Zhang2018). For correlation analysis, Pearson correlation coefficients were calculated between the expression levels of DEGs and the relative abundances of DAMs across matched sample groups. Correlation heatmaps were generated using the corrplot package in R, and only gene-metabolite pairs with |r| > 0.6 and P < 0.05 were considered statistically significant. These correlations were visualized to highlight potentially functionally associated gene-metabolite modules. To further explore the biological relevance of these associations, KEGG-based pathway enrichment analysis was performed using both gene and metabolite identifiers. Joint pathway enrichment was conducted using MetaboAnalyst 5.0 (www.metaboanalyst.ca), integrating transcriptomic and metabolomic datasets to identify significantly enriched metabolic and signalling pathways (P < 0.05, impact value > 0.1). The multi-omics integration enabled a systems-level understanding of host–pathogen interactions, linking transcriptomic changes in infected erythrocytes and lymphocytes to downstream metabolic rewiring.
Results
Transcriptomic and metabolomic features of 2 developmental stages of T. annulata
To help the reader follow the experimental design, we first summarize the 3 core contrasts analysed in this study: (1) purified parasite stages (schizonts vs piroplasms); (2) uninfected vs infected bovine leukocytes (PBMC vs TaXJS) and (3) uninfected vs infected bovine erythrocytes (RBC vs TaXJSM), each profiled by transcriptomics and/or metabolomics as detailed below. To define parasite stage-resolved molecular features, we profiled schizont-enriched material (from infected leukocytes) and piroplasm-enriched material (from infected erythrocytes) (Fig.1). In schizont-enriched preparations, we detected expression of 1537 parasite genes and identified 1622 metabolites. In the piroplasm stage, 1553 parasite genes and 1193 metabolites were detected from purified piroplasms. Comparative analyses between SCHZ and PIRO delineated stage-associated transcriptomic and metabolomic patterns, providing a stage-resolved resource for T. annulata across its intracellular life cycle.
Stage-specific transcriptomic and metabolomic profiling of Theileria annulata during schizont and piroplasm development. The life cycle of T. annulata in bovines includes 2 intracellular stages: schizonts in leukocytes and piroplasms in erythrocytes. A total of 1537 parasite genes and 1622 parasite-derived metabolites were detected in purified schizonts, with 11 038 host genes and 2627 host metabolites identified in infected leukocytes. In the piroplasm stage, 1553 parasite genes and 1193 parasite metabolites were detected from purified piroplasms, and 9450 host genes and 2926 host metabolites were identified in infected erythrocytes.
Stage-specific transcriptomic and metabolomic remodelling across 2 developmental stages of T. annulata
To profile parasite stage-specific molecular signatures, we analysed schizont-enriched (from infected leukocytes) and piroplasm-enriched (from infected erythrocytes) preparations and integrated transcriptomic with metabolomic readouts (Fig. 2). PCA on parasite transcript counts revealed clear separation between schizont and piroplasm samples, indicating distinct stage-associated expression patterns. At the transcript level, the T. annulata reference genome encodes 3717 (Fig.2A). Differential expression analysis identified 125 genes enriched in schizonts and 104 genes enriched in piroplasms (Fig. 2A; thresholds in Methods). At the metabolite level, 310 DAMs were observed between stages (156 higher in schizonts; 154 higher in piroplasms). These data are consistent with stage-associated metabolic differences; however, potential contributions from sample handling cannot be excluded (Fig. 2B). Differential metabolites were predominantly lipids and lipid-like molecules (31.29%) and organic acids and derivatives (27.74%), followed by organoheterocyclic compounds (11.61%) and benzenoids (8.06%), consistent with alterations in membrane-related and central metabolic pathways (Fig. 2C). Hierarchical clustering highlighted divergent expression modules across stages (Fig. 2D). GO enrichment of stage-enriched genes pointed to terms related to transcriptional regulation, catalytic/ATPase/peptidase activity, membrane components and proteasome complexes (Fig. 2E). A volcano plot (Fig. 2F) illustrates representative stage-enriched transcripts (e.g. TA21045, TA11405, TA17125 for schizonts; TA03870, TA13515, TA09450 for piroplasms). Many of these map to subtelomeric or hypothetical families whose specific functions remain poorly characterized, and we therefore do not ascribe direct roles in metabolite processing. KEGG enrichment of differential metabolites showed greater representation of glutathione metabolism, arginine biosynthesis, oxidative phosphorylation and the TCA cycle in schizont samples, along with carbon metabolism, cysteine/methionine metabolism and amino acid biosynthesis (Fig. 2H; Table 2). A circular heatmap (Fig. 2G) visualizes metabolite patterns, with schizont samples showing relatively higher abundance of features linked to nucleotide and amino acid metabolism and redox-related processes, whereas piroplasm samples display a narrower and less diverse metabolite profile. We note that stage-enrichment procedures may introduce stress; matched processing and biological replicates were used, and interpretations remain correlative. Representative differential metabolites include L-Valyl-L-phenylalanine, conjugated linoleic acids, N-dodecanoyl-N-methylglycine and dimethisterone in schizonts, and indole, 5-hydroxytryptophan, ophthalmic acid and other stress-associated metabolites in piroplasms (Fig. 2I).
Summary of changes in transcript and metabolite abundance in Theileria annulata in development. (A) Venn diagrams of differentially expressed transcripts (left) and metabolites (right) between the schizont and piroplasm stage. (B) Summary of the number of significant changes in transcripts and metabolites between the schizont and piroplasm stage. (C) Classification of 310 differential metabolites between the schizont (SCHZ) and piroplasm (PIRO) stages according to HMDB Class I chemical categories (chemical superclasses), showing the number and proportion of metabolites in each category. (D) Heatmap of transcriptomic profiles in parasites at the 2 developmental stages; SCHZ represents schizont-stage parasites, and PIRO represents piroplasm-stage parasites. (E) KEGG pathway enrichment analysis of differentially expressed transcripts between the 2 stages. (F) Volcano plot of transcriptomic differences between the 2 stages. (G) Heatmap of metabolite expression profiles at different growth stages of schizont and piroplasm stages of T. annulata and the circle graph. The primary class, metabolite number and percentage split sectors in the circular heatmap. The dendrograms in the second round denote the overall similarity of metabolite expression profiles by primary classes and the groups. P-values below 0.05. (H) KEGG pathway enrichment of differentially abundant metabolites between schizont and piroplasm stages. (I) Volcano plot showing significantly altered metabolites, with key stage-specific compounds labelled. SCHZ, schizont-enriched samples; PIRO, piroplasm-enriched samples. In heatmaps (D, G), the colour scale indicates row-scaled transcript or metabolite abundance from low to high (as shown by the accompanying colour bar). In volcano plots (F, I), points are colour-coded by regulation status (increased, decreased and non-significant features). In bubble plots (E, H), dot colour represents −log10(P-value) and dot size reflects the number of transcripts or metabolites in each pathway.
Transcriptomic alterations in host leukocytes and erythrocytes before and after T. annulata infection
To characterize host responses, we restricted host RNA-seq analyses to the leukocyte compartment. Peripheral blood mononuclear cells from the uninfected animal (PBMC) and TaXJS-infected leukocytes from the infected animal (TaXJS) were subjected to RNA-seq of bovine leukocytes detected expression of 11 038 bovine genes; differential expression between infected (TaXJS) and uninfected PBMCs was assessed (criteria in Methods; full lists in Table S1). Metabolomic profiling of PBMCs identified 2627 metabolites, of which 1111 were differentially abundant (595 up, 516 down). PCA separated infected from control leukocyte samples (Fig. 3A), and hierarchical clustering of DEGs showed clear group stratification (Fig. 3D). GO enrichment implicated immune-related processes (e.g. cytokine production, immune response), mitochondrial organization and regulation of oxidative stress (Fig. 3E). KEGG enrichment further supported immune-signalling involvement (e.g. NF-κB, T-cell receptor, Toll-like receptor, JAK-STAT) (Fig. 3E). These observations are consistent with published reports of schizont-associated transcriptional reprogramming in leukocytes and indicate concomitant metabolic remodelling; we interpret these patterns as associations, not proof of causality.
Transcriptomic comparison between healthy and Theileria annulata-infected bovine lymphocytes and erythrocytes. (A) Principal component analysis (PCA) of transcriptomic profiles from healthy peripheral blood mononuclear cells (PBMCs) and red blood cells (RBCs), and those infected with T. annulata (TaXJS and TaXJSM). (B) Venn diagram showed the overlapped identified the number of shared and unique transcripts between healthy PBMCs and TaXJS-infected lymphocytes. (C) Venn diagram showed the number of shared and unique transcripts between healthy RCs and TaXJSM-infected erythrocytes. (D) The terms enriched in significantly genes in infected T. annulata cattle lymphocyte compared with healthy cattle. (E) The top 20 enriched KEGG pathways of the genes whose expression was significantly upregulated in infected cattle lymphocyte compared with healthy cattle lymphocyte. The orange dots represent significant KEGG enrichment. (F) The terms enriched in genes significantly in infected T. annulata cattle erythrocyte compared with healthy cattle erythrocyte. (G) The top 20 enriched KEGG pathways of the genes whose expression was significantly downregulated in infected cattle erythrocyte compared with healthy cattle erythrocyte. The orange dots represent the significantly enriched KEGG pathways. The green, blue and yellow boxes represent the BP, CC and MF GO terms, respectively, in (D) and (F). (H) KEGG analyses for the DEGs in TaXJS v PBMC and TaXJSM v RBC. Different colours show the values of FDR. PBMC, peripheral blood mononuclear cell; RBC, red blood cell; TaXJS, T. annulata-infected leukocytes; TaXJSM, T. annulata-infected erythrocytes. In the GO plots (D, F), green, blue and yellow boxes indicate BP, CC and MF terms, respectively, and orange dots mark significantly enriched KEGG pathways. In the KEGG FDR heatmap (H), the colour scale denotes FDR values from low to high as indicated by the colour bar.
RNA-seq libraries were also generated from uninfected (RBC) and infected (TaXJSM) erythrocytes; however, because mature bovine erythrocytes are anucleate, we did not perform host differential-expression analysis in RBCs. Instead, (1) parasite transcripts were quantified from the TaXJSM (piroplasm) stage and (2) host-level differences in the erythrocyte compartment were evaluated by metabolomics. Metabolomic profiling of erythrocytes detected 2926 metabolites, with 726 differentially abundant (259 up, 467 down). PCA indicated separation of uninfected vs infected erythrocyte samples (Fig. 3A). In the transcriptomic domain, signals attributed to the erythrocyte compartment refer to parasite transcripts recovered from TaXJSM, not to host RBC gene regulation. Accordingly, we do not infer host pathways that are incompatible with anucleate erythrocytes (e.g. p53 signalling): any such KEGG labels associated with the TaXJSM dataset reflect parasite-stage annotations rather than host RBC biology. Accordingly, throughout the manuscript, ‘erythrocyte transcriptomic signals’ refer exclusively to the piroplasm (parasite) compartment rather than to host RBC gene regulation. Functionally, the erythrocyte compartment lacks classical immune signalling capabilities, and the metabolomic contrast indicates a more restrained and less diverse metabolic profile relative to leukocytes. Because purification can introduce handling-related stress, we cannot fully exclude potential contributions from sample processing; paired workflows and biological replicates were used to mitigate this.
Venn analyses showed that PBMC vs TaXJS shared 12 235 transcripts, with 2707 and 1281 unique to PBMCs and TaXJS, respectively (Fig. 3B). For the erythrocyte compartment, RBC vs TaXJSM co-detected 10 540 transcripts, with 2639 and 1093 unique to RBC and TaXJSM, respectively (Fig. 3C); given the anucleate nature of mature RBCs and the dominance of parasite-derived RNA in TaXJSM, these counts are interpreted cautiously and are not used for host DE calling, but rather to illustrate the detectable piroplasm transcript repertoire. To compare pathway-level trends across contrasts, we generated a KEGG FDR heatmap (Fig. 3H): immune system pathways were prominently enriched in TaXJS vs PBMC, whereas metabolic pathways-notably oxidative phosphorylation, amino acid biosynthesis and redox-related processes were selectively enriched in TaXJSM vs RBC. Together, these results emphasize stage- and cell-type associated host responses to T. annulata infection immune-focused programs in leukocytes and metabolism-focused shifts in the erythrocyte compartment while avoiding causal claims and clarifying species/context for each analysis.
Dynamic changes in host-cell metabolic profiles during T. annulata infection
To investigate the impact of T. annulata infection on host-cell metabolic reprogramming, we performed untargeted metabolomic profiling on uninfected PBMCs and red blood cells (RBCs), alongside their infected counterparts (XJS and XJSM, respectively). PCA revealed clear separations between infected and uninfected groups. In leukocytes, infected XJS samples were distinctly separated from PBMCs, with PC1, PC2 and PC3 explaining 59.61%, 15.59% and 14.49% of the total variance, respectively. Similarly, infected erythrocytes (XJSM) were well-distinguished from uninfected RBCs, with the top 3 components accounting for 50.35%, 21.00% and 17.42% of the variance (Fig. 4A, 4D), indicating that T. annulata induces profound metabolic divergence in both immune and non-immune host cells. Metabolite classification analysis revealed infection-induced alterations in the distribution of major metabolite categories. In both T. annulata-infected leukocytes and erythrocytes, the relative proportions of organic acids and derivatives, lipids and lipid-like molecules, and nucleosides/nucleotides were significantly shifted compared to their respective healthy controls (Fig 4B, 4E). These compositional changes suggest infection-induced disruptions in central metabolic processes, including energy production, nucleotide metabolism and lipid remodelling. KEGG annotation of significantly changed metabolites further indicated that most differential compounds mapped to broad metabolic pathways and amino acid metabolism. In infected leukocytes, 227 differential metabolites were annotated to general metabolic pathways and 140 to amino acid metabolism. In infected erythrocytes, similar enrichment was observed, with additional emphasis on lipid-associated processes such as sphingolipid metabolism and bile acid biosynthesis (Fig. 4C, 4F). Differential metabolite analysis identified 1111 significantly altered metabolites in XJS vs PBMC (595 upregulated and 516 downregulated) and 726 in XJSM vs RC (259 upregulated and 467 downregulated). Volcano plots illustrated extensive metabolic reprogramming in both cell types, with infection-specific elevation of thymidine, glycohyodeoxycholic acid and valmorin B in leukocytes, and upregulation of 3-hydroxydecanoic acid and taurodeoxycholate in erythrocytes (Fig. 4G, 4J). These profiles highlight distinct sets of metabolite markers associated with intracellular parasitism in different host environments. Clustering analysis of the top altered metabolites further confirmed robust separation between infected and uninfected samples across replicates. Heatmaps revealed consistent accumulation or depletion patterns in selected amino acid derivatives and xenobiotic-like features (for example cleistanthin B and glycyl-methionine) in TaXJS, and orotidine, indole-5-carboxylic acid and 3-hydroxylaurate in TaXJSM (Fig. 4H, 4K). Pathway enrichment analysis based on KEGG revealed distinct metabolic responses between cell types following T. annulata infection. In infected leukocytes (Fig. 4I), significantly enriched pathways included tryptophan metabolism, pyrimidine metabolism, glycerophospholipid metabolism and beta-alanine metabolism. These pathways are closely associated with immune activation, nucleotide biosynthesis and membrane remodelling-hallmarks of lymphocyte functional reprogramming during infection. In contrast, infected erythrocytes exhibited significant enrichment in pathways such as xenobiotic metabolism by cytochrome P450, bile secretion, glutathione metabolism, ferroptosis and neuroactive ligand-receptor interaction (Fig. 4L). These findings are further supported by comprehensive pathway enrichment statistics across developmental stages and ion modes (Table S3). These changes reflect adaptations to oxidative stress, limited biosynthetic capacity and the immunologically inert nature of RBCs, which may facilitate parasite persistence and nutrient scavenging. Together, these findings underscore the profound metabolic rewiring triggered by T. annulata in a cell type-specific manner, with lymphocytes engaging in biosynthesis and immunometabolism, while erythrocytes activate detoxification and redox-balancing pathways to accommodate parasite survival.
Metabolomic comparison between healthy and Theileria annulata-infected bovine lymphocytes and erythrocytes. (A) Principal component analysis (PCA) of the metabolomic profiles from healthy bovine lymphocytes (PBMCs) and T. annulata-infected lymphocytes (XJS), including quality control (QC) samples. Each dot represents a biological replicate. (B) Pie chart showed the classification of metabolites identified in PBMCs and XJS. (C) KEGG pathway enrichment analysis of differential metabolites between PBMCs and XJS. (D) PCA of metabolomic profiles from healthy bovine erythrocytes (RC) and T. annulata-infected erythrocytes (XJSM), including QC samples. (E) Pie chart showing the classification of metabolites identified in RBCs and XJSM. (F) KEGG pathway enrichment analysis of differential metabolites between RBCs and XJSM. (G) Volcano plot of differential metabolites between PBMCs and XJS. Red and blue dots indicate significantly upregulated and downregulated metabolites, respectively; grey dots represent non-significant metabolites. (H) Heatmap showing the abundance of selected differential metabolites in PBMC and XJS groups. (I) KEGG bubble plot showing significantly enriched pathways based on differential metabolites between PBMCs and XJS. Bubble size reflects the number of metabolites involved; colour scale indicates adjusted P-values. (J) Volcano plot of differential metabolites between RBCs and XJSM. (K) Heatmap of selected differential metabolites in RC and XJSM groups. (L) KEGG bubble plot of enriched metabolic pathways for differentially abundant metabolites between RC and XJSM. PBMC, peripheral blood mononuclear cell; XJS, TaXJS-infected leukocytes; RC, red blood cell; XJSM, TaXJSM-infected erythrocytes; QC, pooled quality control sample. In PCA plots (A, D), point colours distinguish the indicated sample groups and QC injections, with each point representing 1 biological replicate. In heatmaps (H, K), the colour scale represents relative metabolite abundance (low to high, as shown in the colour bar). In volcano plots (G, J), red and blue dots indicate significantly increased and decreased metabolites, respectively, whereas grey dots denote non-significant metabolites. In KEGG bubble plots (C, F, I, L), bubble size reflects the number of metabolites and bubble colour encodes statistical significance (−log10(P-value)).
Integrative analysis of host transcriptomic and metabolomic alterations during T. annulata infection
To gain insight into host cellular responses to T. annulata infection, we performed integrative transcriptomic and metabolomic analyses comparing infected bovine leukocytes (TaXJS) with uninfected PBMCs, and infected erythrocytes (TaXJSM) with uninfected RBCs. Differential expression analysis revealed extensive reprogramming at both the transcript and metabolite levels. In TaXJS cells, 6480 genes and 516 metabolites were significantly downregulated, while 4358 genes and 595 metabolites were upregulated (Figure 5A). Similarly, in TaXJSM cells, 5282 genes and 467 metabolites were downregulated, with 4168 genes and 299 metabolites upregulated (Fig. 5E), indicating that infection induced robust changes in host cell biosynthetic and metabolic programs. To identify shared patterns of variation across omics layers, we conducted co-inertia analysis, which revealed strong coupling between DEGs and metabolites in both comparisons (Fig. 5B, F). In TaXJS cells, 5-hydroxytryptophol, L-malate, N-acetyl-D-glucosamine and 5-oxoproline contributed strongly to the joint variance, suggesting disturbances in amino acid turnover and mitochondrial function. In TaXJSM cells, glutathione, nicotinamide, indolelactic acid and taurodeoxycholate were among the top co-contributing metabolites, implicating redox imbalance, bile acid metabolism and microbial by-products in infected erythrocyte physiology. Integrated KEGG pathway enrichment analysis further delineated infection-induced remodelling. In TaXJS, pathways enriched by both transcriptomic and metabolomic changes included tryptophan metabolism, ferroptosis, pyrimidine metabolism, serotonergic synapse, and vitamin digestion and absorption (Fig. 5C). These pathways are associated with immune activation, redox signalling and nucleotide turnover. In contrast, TaXJSM cells showed enrichment in glutathione metabolism, xenobiotic metabolism via cytochrome P450, amino sugar and nucleotide sugar metabolism, and bile secretion (Fig. 5G), indicating that detoxification, antioxidant defenses and membrane remodelling were dominant metabolic features in infected erythrocytes. To illustrate specific infection-driven alterations in canonical pathways, we mapped expression and abundance changes onto KEGG reference diagrams. In TaXJS, transcripts encoding components of mitochondrial oxidative phosphorylation, including Complex I (Ndufa1, Ndufb4), Complex III (Uqcrq) and Complex V (Atp5g1, Atp5h), were broadly downregulated (Fig. 5D), suggesting impaired energy production and mitochondrial dysfunction. Conversely, in TaXJSM, the glutathione metabolism pathway was extensively disrupted at both the transcriptional and metabolite levels (Figure 5H), with reductions in GSR, GPX and NADPH-related enzymes, along with decreased levels of reduced glutathione and its precursors, reflecting weakened antioxidant capacity under parasitic stress. These results demonstrate that T. annulata infection elicits distinct, cell type-specific molecular programs: lymphocytes mount an immune-metabolic response while undergoing mitochondrial suppression, whereas erythrocytes exhibit redox and xenobiotic stress adaptations. This integrative multi-omics view underscores the parasite’s ability to exploit divergent host cell environments through tailored metabolic manipulation.
Integrated transcriptomic and metabolomic analysis of healthy and Theileria annulata-infected bovine lymphocytes and erythrocytes. (A) Bar chart showing the number of differentially expressed genes (DEGs) and differential metabolites (DMs) between T. annulata-infected leukocytes (TaXJS) and healthy bovine leukocytes (PBMC). A total of 4558 upregulated and 6480 downregulated genes, and 595 upregulated and 516 downregulated metabolites were identified. (B) O2PLS (2-way orthogonal partial least squares) integration analysis of transcriptomic and metabolomic data between TaXJS and PBMC. The top 10 contributing genes and metabolites are annotated. Yellow triangles indicate genes; green circles indicate metabolites. The distance from each point to the origin reflects its weight in the integration analysis – the farther from the origin, the greater its contribution. Closer points indicate stronger correlations. (C) KEGG enrichment bubble plot of co-enriched pathways from transcriptomic and metabolomic data (positive and negative ion modes) between TaXJS and PBMC. The x-axis represents the ratio of DEGs or DMs in a given pathway to the total number of genes/metabolites annotated in that pathway. The y-axis lists the significantly co-enriched KEGG pathways. Dot size represents the number of DEGs or DMs; colour indicates statistical significance (−log10(P-value)). (D) Pathway visualization using KEGG Pathview for the oxidative phosphorylation pathway between TaXJS and PBMC, incorporating transcriptomic and metabolomic changes. (E) Bar chart showing the number of DEGs and DMs between T. annulata-infected erythrocytes (TaXJSM) and healthy bovine erythrocytes (RC). A total of 4498 upregulated and 7629 downregulated genes, and 259 upregulated and 467 downregulated metabolites were identified. (F) O2PLS integration analysis between TaXJSM and RC. The top 10 genes and metabolites contributing most to the joint variation are labelled. Interpretations follow panel (B). (G) KEGG enrichment bubble plot showing jointly enriched pathways of DEGs and DMs between TaXJSM and RC, with annotations and colour/size interpretations as in panel (C). (H) KEGG Pathview visualization of the glutathione metabolism pathway between TaXJSM and RC, integrating both transcriptomic and metabolomic changes. DEG, differentially expressed gene; DM, differentially abundant metabolite. In bar charts (A, E), separate bars indicate upregulated and downregulated features as labelled on the axes. In the O2PLS loading plots (B, F), yellow triangles indicate genes and green circles indicate metabolites. In KEGG bubble plots (C, G), bubble size reflects the number of DEGs or DMs per pathway and bubble colour indicates statistical significance (−log10(P-value)). In Pathview diagrams (D, H), coloured boxes represent genes or metabolites with relative increases or decreases between infected and control groups, as indicated by the colour scale bar.
Discussion
The apicomplexan parasite, T. annulata (Witschi et al., Reference Witschi, Xia, Sanderson, Baumgartner, Wastling and Dobbelaere2013), is a tick-borne intracellular pathogen of cattle in tropical and subtropical regions, undergoing proliferative development within both leukocytes and erythrocytes (Kinnaird et al., Reference Kinnaird, Logan, Kirvar, Tait and Carrington1996; Ivan Morrison, Reference Ivan Morrison2009; Jia et al., Reference Jia, Zhao, Xie, Li, Wang and Zhang2020), which can cause acute and often fatal disease. Beyond its economic impact on livestock health, T. annulata is particularly notable for its unique capacity to transform host leukocytes, a feature that distinguishes it from many other protozoan parasites (Kinnaird et al., Reference Kinnaird, Logan, Kirvar, Tait and Carrington1996; Ivan Morrison, Reference Ivan Morrison2009). This transformation, driven by the schizont stage of the parasite, induces a cancer-like phenotype in infected leukocytes. Moreover, the ease of maintaining and passaging these transformed cells in vitro makes T. annulata an attractive model for studying host–pathogen interactions and leukocyte transformation (Elati et al., Reference Elati, Tajeri, Obara, Mhadhbi, Zweygarth, Darghouth and Nijhof2023). However, according to current reports, most studies have focused on either a single time point or a single host cell type, and thus offer only a partial view of the transcriptional and metabolic changes occurring during infection. A more comprehensive understanding of the dynamic interplay between parasite development and host responses, particularly across different parasite stages and infection states, remains limited. To help fill this gap, our study investigated the host–parasite interaction across distinct developmental stages of T. annulata and corresponding host infection states. We performed an integrative multi-omics analysis combining transcriptomics and metabolomics to profile the stage-specific host responses in both bovine PBMCs and erythrocytes. By capturing molecular features of both the schizont and piroplasm stages, we provide a multilayered landscape of host reprogramming during infection. A key strength of this work lies in its parallel analysis of leukocytic and erythrocytic responses to T. annulata, which revealed striking differences in transcriptional and metabolic remodelling. Infected leukocytes exhibited immune activation, oxidative stress adaptation and mitochondrial suppression, whereas infected erythrocytes displayed redox imbalance, lipid remodelling and detoxification responses despite lacking nuclei and organelles. These contrasting host cell reactions underscore the parasite’s remarkable plasticity and its ability to tailor survival strategies to diverse intracellular niches. Our findings offer a systems-level view of T. annulata-host interactions and shed light on the parasite’s regulatory mechanisms and parasitic adaptation, paving the way for deeper understanding and potential intervention strategies.
The integration of transcriptomic and metabolomic datasets across schizont (SCHZ) and piroplasm (PIRO) stages reveals a functionally bifurcated biological program in T. annulata that corresponds with its proliferative and transmissive roles, respectively. Rather than merely indicating transcriptional or metabolic differences, these data collectively underscore a stage-specific strategic reallocation of biosynthetic investment and parasite–host interaction capacity. Consistent with these roles, the schizont stage demonstrates broader transcriptional output and more active metabolic engagement, particularly in pathways related to protein biosynthesis, glycolysis and lipid metabolism. In contrast, the piroplasm stage exhibits reduced parasite transcript diversity and a simplified metabolic profile. Given that mature bovine erythrocytes are anucleate, erythrocyte-compartment RNA is parasite-derived, consistent with limited biosynthetic activity.
Based on the subset of stage-specific genes (Table 1), the schizont stage is characterized by selective and dominant expression of nuclear-targeted effectors, such as TA20085 (Tashat1 protein), TA03155 (Tash1-like protein) and subtelomeric variable secreted proteins (SVSP family members TA09790 and TA09800). These genes are highly expressed in the schizont stage and are typically associated with parasite–host interactions that may support leukocyte transformation and parasite persistence (Schmuckli-Maurer et al., Reference Schmuckli-Maurer, Casanova, Schmied, Affentranger, Parvanova, Kang’a, Nene, Katzer, McKeever, Müller, Bishop, Pain and Dobbelaere2009; Li et al., Reference Li, Liu, Zhao, Ma, Liu, Li, Guan, Luo and Yin2021, Reference Li, Liu, Zhao, Ma, Guo, Liu, Li, Guan, Luo and Yin2022). In line with earlier work, our data show that these loci are strongly downregulated in piroplasms, consistent with a primary role in schizont-mediated host reprogramming. To date, however, no direct involvement of these proteins in metabolic pathways has been demonstrated, reinforcing the view that they mainly act as effectors of host-cell signalling, immune evasion and leukocyte transformation rather than as core metabolic enzymes.
Stage-specific transcriptionally expressed genes of Theileria annulata schizont and piroplasm stages
Note: Genes marked ‘Yes’ indicate stage-specific transcriptional expression in either the schizont or piroplasm stage, while ‘None’ indicates that the gene is not expressed in that stage. Log₂Fold Change and P-values were calculated based on normalized gene expression data using DESeq2, comparing schizont vs piroplasm stages. Gene expression levels were normalized by transcript counts using the DESeq2 median-of-ratios method. Genes with adjusted P-values < 0.05 were considered significantly differentially expressed.
In stark contrast, the piroplasm-specific gene set (e.g. TA03870, TA21045, TA17105) encodes hypothetical or stress-related proteins (Pieszko, Reference Pieszko2015), many of which are associated with conserved subtelomeric DNA sequences and show analogies to transmission-stage markers in Plasmodium (Sargeant et al., Reference Sargeant, Marti, Caler, Carlton, Simpson, Speed and Cowman2006). Several of these loci, including TA03870 and TA17105, were also reported among the most strongly upregulated genes during merozoite-to-piroplasm differentiation in Ankara-derived lines, indicating substantial overlap between our in vivo TaXJS dataset and previous in vitro transcriptome analyses. This expression profile suggests that while schizonts actively drive host manipulation, piroplasms prioritize immunoevasion, minimal metabolic footprint and preparation for vector-stage development. Nevertheless, some parasite genes deviate from this overall pattern. For example, the Ta9 family member TA15705 shows significantly higher transcript abundance in piroplasms than schizonts in our dataset, whereas Ta9 transcripts and protein have previously been detected predominantly in schizont-transformed leukocytes and linked to AP-1 activation and Hck-dependent proliferation (Unlu et al., Reference Unlu, Tajeri, Bilgic, Eren, Karagenc and Langsley2018; Tajeri et al., Reference Tajeri, Shiels, Langsley and Nijhof2025a). Several non-mutually exclusive factors could underlie this discrepancy, including differences in parasite strain background (field-derived, buparvaquone-resistant TaXJS vs laboratory Ankara isolates), in vivo vs in vitro culture conditions, stage-enrichment strategies and potential cross-contamination between fractions, or copy-specific regulation within the multicopy Ta9 family. Given these caveats, we interpret the apparent piroplasm enrichment of TA15705 as a hypothesis-generating observation that requires targeted validation (e.g. copy-specific qRT-PCR and protein localization across synchronized stages) rather than as definitive evidence that Ta9 is predominantly expressed in piroplasms.
On the metabolic axis, the biochemical partitioning supports this regulatory divergence. The classification reveals a clear bias towards lipid and organic acid derivatives in schizonts (Fig. 2C), suggesting membrane biogenesis, organelle expansion and active synthesis-processes consistent with intracellular parasitism and transformation. This is reinforced by the heatmap clustering in Figure 2G, where classes such as ‘lipids and lipid-like molecules’ and ‘organic acids and derivatives’ dominate in the schizont (SCHZ) branch. Importantly, Table 2 reveals KEGG pathway-level distinctions that mirror these roles. The TCA cycle (P = 0.028) and aromatic amino acid biosynthesis pathways are significantly upregulated in schizonts, reflecting heightened energy metabolism and protein synthesis. Additional near-significant enrichments in arginine biosynthesis, purine metabolism and pyruvate metabolism align with increased redox regulation and nucleotide turnover in rapidly proliferating stages. These patterns resonate with prior metabolic profiling of T. annulata-infected cells (Zhao et al., Reference Zhao, Li, Liu, Guan and Dan2022), which noted that proline, arginine and L-carnitine availability are essential for survival and transformation reversal.
Enrichment factor, P-value and number of differential metabolites for each metabolic pathway at schizont and piroplasm stages of Theileria annulata (combined ion modes)
Note: KEGG pathway enrichment analysis based on differential metabolites identified between schizont and piroplasm stages of Theileria annulata. Values are derived from combined positive and negative ion modes (‘all’). Enrichment factor represents the ratio of observed to expected metabolites in each pathway. P-values were calculated using a hypergeometric test. Pathways with P < 0.05 were considered significantly enriched.
In contrast, the piroplasm stage does not merely downregulate schizont-associated pathways; it selectively shifts towards maintaining minimal yet specific metabolic activity. Several metabolites, including uric acid, ophthalmic acid and 5-hydroxytryptophan, were significantly elevated in the piroplasms compared to the schizont stage (Figure 2I). These metabolites are compatible with stress-associated chemistry or catabolic turnover; while a contribution from stage-enrichment/purification cannot be excluded, the overall pattern is consistent with a shift towards detoxification and homeostatic maintenance rather than anabolic biosynthesis (Gallego-Delgado et al., Reference Gallego-Delgado, Ty, Orengo, Van De Hoef and Rodriguez2014; Vincent et al., Reference Vincent, Daly, Courtioux, Cattanach, Biéler, Ndung’u, Bisser and Barrett2016; Servillo et al., Reference Servillo, Castaldo, Giovane, Casale, D’Onofrio, Cautela and Balestrieri2018). Collectively, these profiles are consistent with schizonts showing broader transcript output and greater engagement of biosynthetic/central-carbon pathways with stage-enriched host-interaction effectors, whereas piroplasms exhibit reduced transcript diversity and a compact, maintenance-oriented metabolic profile emphasizing redox/homeostatic buffering. This division reflects a broader paradigm across apicomplexa, wherein developmental transitions are orchestrated by modular transcriptional and metabolic rewiring in response to environmental and host-derived cues, allowing stage-specific adaptation (Gubbels et al., Reference Gubbels, Coppens, Zarringhalam, Duraisingh and Engelberg2021).
To gain deeper insights into host–parasite dynamics, we further explored transcriptional and metabolic reprogramming in bovine PBMCs and erythrocytes upon infection by T. annulata. Our integrated transcriptomic and metabolomic analyses demonstrated distinct cell type-specific response patterns, highlighting how the parasite precisely adjusts its infection strategy according to the host’s metabolic environment. In the leukocytes infection model (TaXJS vs PBMC), transcriptomic enrichment analysis revealed strong activation of immune and inflammatory pathways, including antigen processing and presentation, cytokine–cytokine receptor interaction, and NF-κB signalling, alongside hematopoietic cell lineage remodelling (Fig. 3D, 3E). GO enrichment analyses revealed robust activation of pathways linked to cytokine signalling, immune response and inflammation, including NF-κB signalling, graft-vs-host disease and hematopoietic cell lineage pathways. These findings are consistent with previous studies demonstrating that T. annulata drives malignant transformation of host leukocytes through pro-inflammatory and immunomodulatory signalling cascades (Chaussepied et al., Reference Chaussepied, Janski, Baumgartner, Lizundia, Jensen, Weir, Shiels, Weitzman, Glass, Werling and Langsley2010; Durrani et al., Reference Durrani, Weir, Pillai, Kinnaird and Shiels2012; Rchiad et al., Reference Rchiad, Haidar, Ansari, Tajeri, Mfarrej, Ben Rached, Kaushik, Langsley and Pain2020; Ahlawat et al., Reference Ahlawat, Choudhary, Arora, Kumar, Kaur and Chhabra2023; Elati et al., Reference Elati, Tajeri, Obara, Mhadhbi, Zweygarth, Darghouth and Nijhof2023). Correspondingly, metabolomic data demonstrated significant enrichment of amino acid metabolism – especially tryptophan metabolism – as well as enhanced glycolysis and glycerophospholipid metabolism following infection (Fig. 4G–I). These findings reinforce the notion that T. annulata profoundly reshapes host metabolic landscapes to support its proliferation and survival, corroborating previous observations of parasite-driven metabolic reprogramming targeting energy supply, membrane remodelling and signal regulation (Zhao et al., Reference Zhao, Li, Liu, Guan and Dan2022). Distinct from PBMCs, infected erythrocytes (TaXJSM vs RBC), despite their absence of nuclei and mitochondria, exhibited evident transcriptional signals and specialized metabolic adaptations. Transcriptome analysis revealed significant enrichment of ribosome-related pathways, aminoacyl-tRNA biosynthesis and nucleotide metabolism, indicative of translational activity and protein synthesis in infected erythrocytes (Fig. 3F, 3G). Meanwhile, metabolomics highlighted significant changes in glutathione metabolism, bile acid secretion, arachidonic acid metabolism and ferroptosis-associated pathways (Fig. 4L). Crucial metabolites such as taurodeoxycholate, glutathione and indolelactic acid were markedly elevated in infected erythrocytes, implicating their roles in antioxidant stress responses, detoxification and maintaining parasite viability (Becker et al., Reference Becker, Tilley, Vennerstrom, Roberts, Rogerson and Ginsburg2004; Georgiou-Siafis and Tsiftsoglou, Reference Georgiou-Siafis and Tsiftsoglou2023; Pawłowska et al., Reference Pawłowska, Mila-Kierzenkowska, Szczegielniak and Woźniak2023; Mohammed et al., Reference Mohammed, Kuraa, Rushdi and Abou El-Ella2025). The conspicuous enrichment of glutathione metabolism pathways, in particular, likely represents a parasite-induced protective mechanism safeguarding erythrocytes against oxidative damage, thereby ensuring stable, long-term parasitic residence(Mehlotra, Reference Mehlotra1996).
T. annulata infection induces cell type-specific yet convergent metabolic adaptations centred on energy production and oxidative stress mitigation (Fig. 5). In infected lymphocytes (TaXJS vs PBMC), oxidative phosphorylation pathways showed broad downregulation at the transcript level. Multiple components of the host mitochondrial electron transport chain including Complex I (NDUFS1, NDUFS2), Complex III (UQCRC1, CYC1), Complex IV (COX6C, COX7C) and ATP synthase subunits (ATP5C1, ATP5H) were reduced in TaXJS compared with PBMC (Fig. 5D). This pattern is consistent with parasite-induced impairment of host mitochondrial respiratory capacity and a shift towards alternative energy-producing pathways during leukocyte transformation. These observations reflect parasite-induced reprogramming of host mitochondrial metabolism, potentially supporting the bioenergetic demands associated with leukocyte transformation (Woods et al., Reference Woods, Perry, Brühlmann and Olias2021). Concurrent enrichment of glutathione metabolism and NADPH-generating pathways further indicates a host compensatory response to infection-associated oxidative stress and mitochondrial dysfunction (Diez et al., Reference Diez, Traikov, Schmeisser, Adhikari and Kurzchalia2021; Kapoor, Reference Kapoor2022). In infected erythrocytes (TaXJSM vs RBC), the integrated transcriptomic-metabolomic analysis pointed to an active glutathione-based antioxidant system in the infected erythrocyte compartment (Figure 5H). Within the TaXJSM transcriptome, genes annotated to glutathione metabolism including glutathione peroxidase (GPX, EC: 1.11.1.9), glutathione reductase (GSR, EC: 1.8.1.7) and the rate-limiting enzyme glutamate cysteine ligase (GCLC, EC: 6.3.2.2) together with enzymes sustaining NADPH production, tended to show higher expression compared with the RBC compartment (Becker et al., Reference Becker, Rahlfs, Nickel and Schirmer2003; Gupta et al., Reference Gupta, Pandey, Kumar and Tripathi2015). At the metabolite level, reduced glutathione and its precursor γ-glutamylcysteine, as well as several related redox-active intermediates, accumulated in infected erythrocytes (Müller, Reference Müller2015). Taken together, these changes are consistent with strengthened glutathione turnover and redox-recycling capacity in the infected erythrocyte compartment, helping to buffer parasite-induced oxidative stress and preserve erythrocyte integrity within the metabolically constrained and immunologically ‘silent’ environment of the RBC (Pawłowska et al., Reference Pawłowska, Mila-Kierzenkowska, Szczegielniak and Woźniak2023).
Collectively, these integrated results underscore a strikingly cell-type-specific transcriptional and metabolic adaptation of host cells during T. annulata infection (Lüder et al., Reference Lüder, Stanway, Chaussepied, Langsley and Heussler2009). In lymphocytes, the parasite primarily exploits host immune signalling pathways, amino acid metabolism and bioenergetics to promote cellular transformation. Conversely, in erythrocytes, it predominantly reinforces redox homeostasis and detoxification pathways, ensuring sustained parasitism. This finely tuned adaptation to distinct host cell environments represents a fundamental mechanism underlying the successful persistence of T. annulata within its host, offering valuable insights for future therapeutic target identification. Despite the novel insights presented here, our study has several limitations. First, while the integrative omics approach provided mechanistic clues, functional validation of key pathways and molecules remains necessary. Follow-up experiments using CRISPR gene editing, siRNA knockdown or small-molecule inhibitors could help clarify the causal roles of candidate genes in parasite survival and host modulation. Second, our analyses were conducted on a single field-derived TaXJS isolate, which is a naturally buparvaquone-resistant strain (Ma et al., Reference Ma, Liu, Li, Xiang, Wang, Liu, Li, Yin, Guan and Luo2020; Zhang et al., Reference Zhang, Zhao and Cao2022). It therefore remains possible that some parasite and host transcriptional or metabolic signatures described here are influenced by the resistance-associated genetic background and the specific experimental setting, rather than being fully generalizable to buparvaquone-sensitive T. annulata strains or to all field situations. Accordingly, the representativeness of TaXJS for the broader natural T. annulata population, especially outside its endemic region, remains to be determined, and our findings should be interpreted with this limitation in mind. Comparative multi-omics analyses across matched sensitive and resistant isolates will be required to disentangle resistance-specific changes from core infection programmes. Moreover, our multi-omics profiling was restricted to peripheral blood (leukocytes and erythrocytes) and does not capture tissue-specific immune-metabolic interactions or inter-animal variability. Third, although we used pre-chilled solvents and processed infected and control samples in parallel, we cannot fully exclude that sample handling prior to metabolite quenching (e.g. density separation, repeated washing and pelleting steps) may have altered labile metabolites or induced transient stress responses. Thus, some metabolite differences observed between groups might partially reflect handling-related effects rather than purely in vivo steady-state metabolism. Fourth, as this was an untargeted LC-MS/MS experiment, individual metabolite identifications remain subject to ambiguity, particularly for features whose best database matches correspond to pharmaceuticals or cosmetic ingredients. Although we applied a tiered confidence scheme (Levels 1–3) and focus our interpretation on chemically plausible metabolites and pathway-level signatures, some singleton annotations should still be considered tentative and are reported mainly for completeness rather than mechanistic inference. Future studies using larger animal cohorts, multiple tissues and longitudinal sampling are warranted to assess tissue-specific dynamics and systemic immune responses. In addition, time-series sampling combined with single-cell resolution analyses may uncover cell-to-cell heterogeneity and dynamic regulation upon infection (Panagiotou, Reference Panagiotou2024). Finally, the metabolomic analyses were based on a limited number of replicates per contrast (3–4 per group) and on nominal P < 0.05 thresholds without formal multiple-testing correction at the individual-metabolite level. Although we mitigated technical variation by retaining only metabolites with QC CV < 30% and by focusing on concordant changes within pathways, the reported DAMs should therefore be considered exploratory and require validation in larger cohorts. In addition, because metabolomic samples from the same group were acquired in contiguous blocks within each batch rather than in a fully randomized order, subtle confounding between biological group and run order cannot be completely excluded, although QC-based monitoring did not indicate major drift.
Conclusions
This study represents the systematic integration of transcriptomic and metabolomic analyses to investigate the host response to T. annulata infection across 2 developmental stages and distinct host cell types. Our multi-omics approach reveals the complex, cell-type-specific transcriptional and metabolic reprogramming induced by the parasite, highlighting key regulatory pathways involved in immune signalling, redox balance and energy metabolism. These findings significantly expand the current understanding of host–pathogen interactions in the context of T. annulata infection. Importantly, the identification of critical metabolic and signalling pathways provides a robust scientific basis for elucidating parasite survival strategies and host manipulation mechanisms. This work offers a valuable resource for future functional studies and lays the foundation for the development of novel therapeutic strategies or vaccines targeting apicomplexan parasites.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0031182026101619. The RNA-seq data for the 2 developmental stages of a naturally occurring buparvaquone-resistant strain of Theileria annulata have been deposited in the NCBI Gene Sequence Archive (GSA) under accession numbers PRJNA1227583 and PRJNA1278822. The metabolomics datasets are available in the MetaboLights database under accession numbers MTBLS12427, MTBLS12428, MTBLS12429 and MTBLS12430.
Data availability statement
The data supporting the findings of the study are available within the article and its supplementary materials.
Acknowledgements
We gratefully acknowledge the Central Instrument Room at Lanzhou Veterinary Research Institute, Chinese Academy of Agricultural Sciences, for their indispensable support in providing instrumentation and technical assistance throughout this study.
Author contributions
YC, JC and SZ collected the samples and performed the data analysis; YC and JW conceived the study and designed the research; YC wrote the article; QR, JL and QN provided constructive suggestions; GG and HY revised the article.
Financial support
The study was financially supported by the National Natural Science Foundation of China (no. 31972701), the Science Fund for Creative Research Groups (22JR5RA024), the Special Project (22CX8NA011) of Gansu Province, the Innovation Program of the Chinese Academy of Agricultural Sciences (CAAS ASTIP-2021-LVRI), NBCITS (CARS-37), the National Parasitic Resources Center (NPRC-2019–194-30) and the grants from the National Key Research and Development Program of China (2024YFD1800100).
Competing interests
The authors declare that they have no competing interests.
Ethical standards
Animal (cattle) use was approved by the Experimental Animal Ethics Committee, Lanzhou Veterinary Research Institute, Chinese Academy of Agricultural Sciences, China (permit no. LVRIAEC-2024–073) and followed the Regulations for the Administration of Affairs Concerning Experimental Animals in China.
Consent for publication
Not applicable.
Open access
