Abstract Background Growing evidence indicates that N^7-methylguanosine (m^7G) modification plays critical roles in epigenetic regulation. However, no data regarding m^7G modification are currently available in Eimeria tenella, a highly virulent species causing coccidiosis in chickens. Methods In the present study, we explore the distribution of internal messenger RNA (mRNA) m^7G modification in sporulated and unsporulated oocysts of E. tenella as well as its potential biological functions during oocyst development using methylated RNA immunoprecipitation sequencing (MeRIP-seq) and mRNA sequencing (mRNA-seq), and the mRNA-seq and MeRIP-seq data were verified by the quantitative reverse transcription polymerase chain reaction (RT–qPCR) and MeRIP–qPCR, respectively. Results Our data showed that m^7G peaks were detected throughout the whole mRNA body, and the coding DNA sequence (CDS) region displayed the most methylation modification. Compared with unsporulated oocysts, 7799 hypermethylated peaks and 1945 hypomethylated peaks were identified in sporulated oocysts. Further combined analysis of differentially methylated genes (DMGs) and differentially expressed genes (DEGs) showed that there was a generally positive correlation between m^7G modification levels and gene transcript abundance. Unsurprisingly, the mRNA-seq and MeRIP-seq data showed good consistency with the results of the RT–qPCR and MeRIP–qPCR, respectively. Gene Ontology (GO) and pathway enrichment analysis of DEGs with altered m^7G-methylated peaks were involved in diverse biological functions and pathways, including DNA replication, RNA transport, spliceosome, autophagy-yeast, and cAMP signaling pathway. Conclusions Altogether, our findings revealed the potential significance of internal m^7G modification in E. tenella oocysts, providing some directions and clues for later in-depth research. Graphical abstract [41]graphic file with name 13071_2024_6580_Figa_HTML.jpg Supplementary Information The online version contains supplementary material available at 10.1186/s13071-024-06580-3. Keywords: Eimeria tenella, Oocysts, N^7-methylguanosine, Transcriptomics, Methylated RNA immunoprecipitation sequencing Background Chicken coccidiosis, an enteric disease caused by apicomplexan parasites belonging to the genus Eimeria, is a global problem in the poultry industry [[42]1, [43]2]. The disease has been estimated to cost about £10.4 billion at 2016 prices in poultry industry [[44]3]. To date, there are ten Eimeria species infecting chickens, including seven long-recognized Eimeria species and three additional cryptic operational taxonomic units [[45]1, [46]4]. Eimeria tenella, characterized by a unique tropism for the cecum, is among the most economically significant species [[47]5, [48]6]. E. tenella undergoes a strict fecal–oral life cycle, which features an exogenous phase (sporogony) in the environment and an endogenous phase (schizogony and gametogony) within the host [[49]7]. Exogenous and extracellular endogenous life cycle forms include unsporulated oocysts (UO), sporulated oocysts (SO), sporozoites, merozoites, and microgametes [[50]1, [51]8]. These stages exhibit distinct developmental and morphological characteristics [[52]2]. Several studies aiming to decipher the mechanisms governing the biology of different developmental stages through analysis of the genes expressed as well as the levels of gene expression were performed. For example, a previous study identified 3342 out of 7329 genes exhibiting differential expression during development through comparative transcriptome profiling of E. tenella in various developmental stages [[53]9]. The data generated by RNA sequencing (RNA-seq) of E. tenella gametocytes, merozoites, and sporozoites revealed upregulated gametocyte transcription of 863 genes [[54]10]. Notably, various RNA modifications are involved in gene expression regulation during eukaryotic development [[55]11, [56]12]. Over 170 types of RNA modifications have been discovered in eukaryotes thus far, including N^6-methyladenosine (m^6A), methylcytidine (m^5C), N^1-methyladenosine (m^1A) and N^7-methylguanosine (m^7G) [[57]13]. Of which, m^7G methylation was initially identified as a signature modification at the 5′ cap of messenger RNAs (mRNAs) [[58]14, [59]15]. Subsequently, m^7G methylation was also found at internal positions within mRNAs, ribosomal RNAs (rRNAs), and transfer RNAs (tRNAs) [[60]14]. This modification affects diverse aspects of post-transcriptional gene regulation, including RNA stability and splicing [[61]16]. At the cellular level, the involvement of m^7G methylation in the processes of differentiation of pluripotent stem cells has been reported [[62]17]. To date, however, the involvement of m^7G modification in epigenetic regulation during the developmental cycle of E. tenella remains unexplored. Hence, in the present study, we perform a comprehensive analysis of the transcriptome-wide m^7G methylation in sporulated and unsporulated oocysts of E. tenella. Methods Chickens and parasites Sporulated oocysts of E. tenella (SD-01 strain) were kindly provided by Professor Xiaomin Zhao, Shandong Agricultural University, China. The strain was maintained by subsequent passage every 6 months as described previously [[63]18]. One-day-old Hy-line Brown cocks, purchased from a commercial hatchery (Shanxi Kangmu Farm, Jinzhong, China), were reared in a coccidia-free environment. Preparation of E. tenella oocysts Hy-line Brown cocks were infected orally at 14 days old with a dose of 1 × 10^4 sporulated oocysts of E. tenella SD-01 strain. Seven days postinfection, the chickens were euthanized by cervical dislocation. Fresh oocysts (unsporulated oocysts) were recovered from the cecal contents of infected chickens using procedures described before with minor modifications [[64]19]. The purified unsporulated oocysts were divided into two parts, and each part included three individual aliquots (prepared from different pools of unsporulated oocysts). One part was immediately frozen in liquid nitrogen. The other part was sporulated for 3 days by incubation in 2.5% potassium dichromate solution and then frozen in liquid nitrogen. High‑throughput m^7G sequencing, mRNA sequencing Total RNA was extracted separately from sporulated and unsporulated oocysts with TRIzol reagent (Invitrogen, CA, USA) according to the manufacturer’s specifications. Then, rRNA was depleted using the GenSeq^® rRNA Removal Kit (GenSeq, Inc., Shanghai, China). Regarding fragmentation and immunoprecipitation (IP), the GenSeq^® m7G-IP Kit (GS-ET-004, GenSeq Inc., Shanghai, China) was used. For fragmentation, the rRNA depleted samples were treated with fragmentation buffer for 6 min at 70 °C and then halted with stop buffer. Following overnight incubation at −80 °C with PC buffer, PC enhancer, and 75% ethanol, the RNA precipitate was obtained by centrifugation at 15,000g for 25 min at 4 °C. After washing by 75% ethanol, the RNA pellet was allowed to air dry at room temperature and dissolved in nuclease-free water. A portion of the resultant RNA samples was kept as input. The remaining was decapped by treatment with tobacco acid pyrophosphatase (TAP) (M0608S, NEB, Ipswitch, MA) at 37 °C for 1 h and then subjected to immunoprecipitation. For immunoprecipitation, the porcine gastric mucin-conjugated (PGM) magnetic beads were washed with 1× IP buffer and then incubated with the anti-m^7G antibody (RN017M, MBL, Tokyo, Japan) at room temperature (RT) for 1 h. Following washing with 1× IP buffer, 50 µL of 5× IP buffer, and 250 µL of fragmented RNA were added to the PGM magnetic beads and then incubated at 4 °C for 1 h. The PGM beads were washed with 1× IP buffer, added with 55 µL of dig solution and rotated at 4 °C for 45 min. Afterward, the supernatant was transferred to a new tube with preadded MS beads and RLT buffer. Following incubation with ethanol and washing with 75% ethanol, RNA was eluted with nuclease-free water (immunoprecipitated samples) and used for subsequent library construction with the GenSeq^® Low Input RNA Library Prep Kit (GenSeq Inc., Shanghai, China). Meanwhile, RNA input samples without immunoprecipitation were used for RNA library generation with the same kit. Library sequencing was carried out on an Illumina NovaSeq platform at CloudSeq Inc. (Shanghai, China). Data analysis After 3′ adapter trimming and removal of low quality reads by using cutadapt software (v1.18) [[65]20], the clean reads of m^7G-IP and input libraries were compared with the reference genome of E. tenella (ToxoDB-60_EtenellaHoughton2021, [66]https://toxodb.org/toxo/app/downloads/Current_Release/EtenellaHough ton2021/) with the use of HISAT2 software [[67]21]. After mapping, the sequence alignment map (SAM) files were converted to a binary alignment map (BAM) file using SAMtools [[68]22] and then converted to BED files using the bamToBed command of BEDTools (v2.30.0) [[69]23]. Peak calling was performed using model-based analysis of ChIP-seq (MACS) software (v1.4.2) [[70]24]. For the analysis of distribution of m^7G sites across the mRNA landscape, MetaPlotR was used to calculate and scale the distance of peaks relative to transcriptomic features ([71]https://github.com/olarerin/metaPlotR) [[72]25]. Differentially methylated sites between the two groups were identified with diffReps (v1.55.6) with the criteria of a P-value < 0.05 and a fold change (FC) > 2 [[73]26]. Motif analysis was carried out by using discriminative regular expression motif elicitation (DREME) ([74]https://meme-suite.org/meme/tools/dreme). For the input RNA-seq data analysis, raw gene counts were obtained using HTSeq software (v0.9.1) [[75]27], followed by being normalized by the edgeR package (v3.16.5) [[76]28]. Differentially expressed genes (DEGs) were identified according to the following criteria: a fold change > 2 and a P-value < 0.05. DEGs, Differentially methylated genes (DMGs), and common genes between DMGs and DEGs were subjected to Gene Ontology (GO) enrichment analysis using the topGO package (v2.10) [[77]29], and P-values were calculated using the hypergeometric distribution method (significant enrichment was defined by a P-value < 0.05). For Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway annotation, DEGs, DMGs, and common genes between DMGs and DEGs were searched against the KEGG Database at [78]https://www.genome.jp/kegg/kaas, and enrichment P-values were calculated using Fisher’s exact test (significant enrichment was defined by a P-value < 0.05). Quantitative reverse transcription polymerase chain reaction (RT–qPCR) and methylated RNA immunoprecipitation qPCR (MeRIP–qPCR) The expression levels of four genes (ETH2_1435900, ETH2_1516900, ETH2_1361700, and ETH2_1248400) were determined by RT–qPCR using the QuantStudio 5 Dx Real-Time PCR System (Thermo Fisher Scientific, Waltham, MA). Actin (ETH2_1310200) was used as the endogenous control [[79]30]. Cycling conditions were as follows: 10 min at 94 °C, followed by 40 cycles of 10 s at 95 °C and 60 s at 60 °C. The comparative cycle threshold (CT) (2^−ΔΔCt) method was used to calculate the relative RNA expression levels [[80]31]. For validation of the m^7G sequencing results, both input and m^7G immunoprecipitation samples were subjected to RT–qPCR analysis. The relative m^7G enrichment of mRNA was calculated by normalizing to the input: % input = 2^[Ct(input) −Ct(IP)] × 1/DF × 100, where DF is the dilution factor between IP and input samples. The sequences of primers used for RT–qPCR and methylated RNA immunoprecipitation–qPCR (MeRIP–qPCR) are shown in Additional file [81]1: Table S1. Data analysis and statistical tests were carried out using Microsoft Excel and GraphPad Prism 8.0. Numerical data were shown as mean with standard deviation (SD). Differences between the two groups were analyzed by using a Student’s t-test. Statistical significance was defined as a P-value less than 0.05. Results General features of m^7G methylation in unsporulated and sporulated oocysts Through sequencing of methylated RNA immunoprecipitation sequencing (MeRIP-seq) and RNA-seq, 92,388,774–119,399,988 raw reads were generated from each MeRIP-seq dataset and 77,532,536–102,549,048 raw reads were generated from the RNA-seq dataset. After removing low-quality data and successfully mapping to the reference genome, the mean alignment rates were 84.82, 79.56, 77.20, and 75.03% of clean reads for RNA-seq and MeRIP-seq from the unsporulated and sporulated oocysts samples, respectively (Additional file [82]2: Table S2). We analyzed the distribution patterns of m^7G peak in the whole transcriptome, and the results showed that m^7G peaks in each group were highly enriched around start codon region and within CDS region (Fig. [83]1A). We found that 44.82% of m^7G peaks were enriched in CDS region, followed by 31.45% in start codon region, 13.07% in stop codon region, 9.01% in 5′UTR region and 1.65% in 3′UTR region in the SO group (Fig. [84]1B). In the UO group, 41.11% of m^7G peaks were enriched in CDS region, followed by 28.08% in start codon region, 15.72% in stop codon region, 13.59% in 5′UTR region, and 1.5% in 3′UTR region (Fig. [85]1C). Subsequently, the number of m^7G peaks per m^7G-modified gene was analyzed, and the results showed that more than 80% of genes contained one or two m^7G peaks in each group (Fig. [86]1D). In the SO group, ETH2_1256100, a gene encoding a protein containing the kringle domain, harbored the highest number of m^7G peaks (seven peaks) (Additional file [87]3: Table S3). In the UO group, ETH2_0101200, a gene encoding a protein containing the N-terminal chorein domain, harbored the highest number of m^7G peaks (12 peaks) (Additional file [88]3: Table S3). In addition, motif analysis for the two groups was performed using DREME, and the top two motifs are shown in Fig. [89]1E, [90]F. Fig. 1. [91]Fig. 1 [92]Open in a new tab Overview of m^7G modification in the two groups. A Metagene plots displaying the m^7G peak density distribution across the transcripts. Pie charts exhibiting m^7G peak distribution in different transcript segments in sporulated (B) and unsporulated (C) oocysts, respectively. D Proportion of genes harboring different numbers of m^7G peaks in the two groups. E The top two sequence motifs enriched from all the identified m^7G peaks in the SO group. F The top two sequence motifs enriched from all the identified m^7G peaks in the UO group DMGs and GO enrichment analysis Compared with unsporulated oocysts, 7799 hypermethylated peaks within 3087 genes and 1945 hypomethylated peaks within 847 genes were detected in sporulated oocysts (Fig. [93]2A, Additional file [94]4: Table S4). The DMGs were subjected to GO analysis, and the results were grouped into the following three categories: molecular function (MF), cellular component (CC), and biological process (BP). For BP, the genes with hypermethylated peaks were related to microtubule-based process, movement of cell or subcellular component, and microtubule-based movement (Fig. [95]2B). For CC, the genes with hypermethylated peaks were related to cilium, cell projection, and plasma-membrane-bounded cell projection (Fig. [96]2B). For MF, the genes with hypermethylated peaks were involved in hydrolase activity, acting on glycosyl bonds, microtubule motor activity, and binding (Fig. [97]2B). Fig. 2. [98]Fig. 2 [99]Open in a new tab Analysis of differentially m^7G-modified genes between the two groups. A Volcano plots depicting the differential peaks between the studied groups. B GO enrichment analysis of the genes with hypermethylated peaks in the SO group. C GO enrichment analysis of the genes with hypomethylated peaks in the SO group. D KEGG analysis of the genes presenting hypermethylated peaks in the SO group. E KEGG analysis of the genes presenting hypomethylated peaks in the SO group For BP, the genes with hypomethylated peaks were related to chromosome organization, DNA metabolic process, and DNA conformation change (Fig. [100]2C). For CC, the genes with hypomethylated peaks were related to chromosome, supramolecular complex, and AP-type membrane coat adapter complex (Fig. [101]2C). For MF, the genes with hypomethylated peaks were involved in tubulin binding, cytoskeletal protein binding, and microtubule binding (Fig. [102]2C). Pathway enrichment analysis of DMGs KEGG analysis was performed to further investigate the potential biological function of DMGs. The significantly enriched pathways for the hypermethylated genes included glycolysis/gluconeogenesis, biosynthesis of cofactors, phenylalanine, tyrosine and tryptophan biosynthesis, alanine, aspartate and glutamate metabolism, porphyrin and chlorophyll metabolism, tyrosine metabolism, thiamine metabolism, sulfur relay system, fat digestion and absorption, and aminoacyl-tRNA biosynthesis (Fig. [103]2D). The significantly enriched pathways for the hypomethylated genes included DNA replication, GABAergic synapse, cell cycle yeast, base excision repair, pentose and glucuronate interconversions, sphingolipid signaling pathway, nonhomologous end joining (NHEJ), cell cycle, homologous recombination, and meiosis yeast (Fig. [104]2E). Integrated analysis of MeRIP-seq and RNA-seq data Compared with the UO group, 4103 genes were differentially expressed in the SO group, including 2306 upregulated genes and 1797 downregulated genes (Fig. [105]3A, Additional file [106]5: Table S5). Clustering analysis revealed the distinct gene expression patterns (Fig. [107]3B). Based on the combined analysis of the RNA-seq and MeRIP-seq data, all DMGs with differential expression were divided into four groups, including 4387 hypermethylated and upregulated (hyper-up) genes, 571 hypermethylated and downregulated (hyper-down) genes, 16 hypomethylated and upregulated (hypo-up) genes, and 1568 hypomethylated and downregulated (hypo-down) genes (Fig. [108]3C, Additional file [109]6: Table S6). As shown in Fig. [110]3D, there was a generally positive correlation between the methylation levels and gene transcript abundance. Fig. 3. [111]Fig. 3 [112]Open in a new tab A conjoint analysis of MeRIP-seq and RNA-seq data. A Volcano plots of the DEGs between the studied groups. B Hierarchically clustered heat maps showing the expression changes between the studied groups. C Four quadrant graphs showing the DEGs with differentially methylated m^7G peaks. D Correlation analysis between mRNA expression levels and methylation levels Four genes (ETH2_1435900, ETH2_1516900, ETH2_1361700, and ETH2_1248400) were subjected to RT–qPCR analysis, and the results were consistent with the gene expression profiles in transcriptome data (Fig. [113]4A). In addition, a MeRIP–qPCR assay was performed to validate the MeRIP-seq data. The MeRIP–qPCR results were consistent with the sequencing results (Fig. [114]4B). Fig. 4. [115]Fig. 4 [116]Open in a new tab Validation of RNA-seq and MeRIP-seq data. A The mRNA levels of four genes were measured by RT–qPCR. B The m^7G levels of four genes were measured by MeRIP–qPCR. **P < 0.01, ***P < 0.001 GO analysis showed that the hyper-down genes were mainly associated with DNA conformation change, microtubule cytoskeleton, and ATP binding (Fig. [117]5A). The hyper-up genes were primarily associated with RNA metabolic process, nucleus, and DNA-binding transcription factor activity (Fig. [118]5B). The hypo-down genes were involved in DNA replication, chromosome, and tubulin binding (Fig. [119]5C). The hypo-up genes were mainly associated with tRNA aminoacylation for protein translation, integral component of membrane, and aminoacyl-tRNA ligase activity (Fig. [120]5D). Fig. 5. [121]Fig. 5 [122]Open in a new tab GO analysis of DEGs with altered m^7G-methylated peaks. A Bar plots showing the significantly enriched GO terms for the hyper-down genes. B Bar plots showing the significantly enriched GO terms for the hyper-up genes. C Bar plots showing the significantly enriched GO terms for the hypo-down genes. D Bar plots showing the significantly enriched GO terms for the hypo-up genes In addition, KEGG analysis revealed that the hyper-down genes were mainly associated with adrenergic signaling in cardiomyocytes, chloroalkane and chloroalkene degradation, naphthalene degradation, degradation of aromatic compounds, and biosynthesis of cofactors (Fig. [123]6A). The hyper-up genes were primarily associated with alanine, aspartate and glutamate metabolism, plant-pathogen interaction, RNA transport, spliceosome, carbon fixation in photosynthetic organisms, phenylalanine, tyrosine and tryptophan biosynthesis, fluid shear stress and atherosclerosis, lipid and atherosclerosis, phenylalanine metabolism, isoquinoline alkaloid biosynthesis, tropane, piperidine and pyridine alkaloid biosynthesis, circadian rhythm—plant, morphine addiction, PI3K-Akt signaling pathway, and longevity regulating pathway—multiple species (Fig. [124]6B). The hypo-down genes were involved in DNA replication, GABAergic synapse, pentose and glucuronate interconversions, sphingolipid signaling pathway, cell cycle-yeast, NHEJ, cAMP signaling pathway, pathogenic Escherichia coli infection, adrenergic signaling in cardiomyocytes, cell cycle, parathyroid hormone synthesis, secretion and action, hepatitis B, purine metabolism, autophagy yeast, and meiosis yeast (Fig. [125]6C). Fig. 6. [126]Fig. 6 [127]Open in a new tab Pathway analysis of DEGs with altered m^7G-methylated peaks. A Bubble charts showing the significantly enriched pathways for the hyper-down genes. B Bubble charts showing the significantly enriched pathways for the hyper-up genes. C Bubble charts showing the significantly enriched pathways for the hypo-down genes Discussion In the present study, we analyzed, for the first time, the internal mRNA m^7G modification pattern in unsporulated and sporulated oocysts of E. tenella. RT–qPCR and RNA-seq results were in good agreement, demonstrating that the RNA-seq data are reliable. Meanwhile, the reliability of m^7G-RIP-seq data was validated by MeRIP-PCR. Our findings showed that the overall distribution of m^7G modification sites is similar in the transcriptome of unsporulated and sporulated oocysts, providing insights into the conservation of m^7G modification at the two developmental stages. The m^7G peaks were primarily enriched in the CDS region in each group, indicating that m^7G modification within the CDS region may play an important role in post-transcriptional gene expression regulation during oocyst development. Intriguingly, a previous study analyzed the distribution patterns of m^7G peaks across mRNA transcripts in HL60 and HL60/MX2 cells, and the results also showed that the CDS region exhibited the most methylation modification [[128]32]. It was found that the m^6A peaks were also highly enriched in the CDS region [[129]33]. However, differential distribution patterns were observed. For example, m^6A enrichment around the stop codon region ranked second [[130]33], whereas m^7G enrichment around the start codon region ranked second. A start codon is related to translation initiation, while a stop codon ends it [[131]34]. This indicated that m^7G and m^6A modifications may play distinct and cooperative roles to ensure accurate mRNA translation initiation and termination during oocyst development, and further studies are needed to better understand the underlying molecular mechanisms. Motif analysis showed that the top two significantly enriched motifs in each group were similar. Strikingly, there are similarities between the motifs identified in the present study and those reported previously in mammalian cells [[132]35]. A large number of DMGs were identified in the SO group compared with the UO group, with hypermethylated genes accounting for the majority. Among them, many genes contained more than one m^7G peaks. For example, ETH2_0701600 contained two hypermethylated peaks and three hypomethylated peaks. This indicated the involvement of m^7G modification in post-transcriptional regulation of ETH2_0701600. Multiple processes in Plasmodium chabaudi (such as DNA replication, cell cycle, and microtubule-based movements) during the intraerythrocytic developmental cycle were reported to be affected by disruption of serpentine receptor 10 (PcSR10) [[133]36]. Hence, ETH2_0701600, the SR10 homolog in E. tenella, may play a regulatory role in the process of sporulation through affecting DNA replication, cell cycle, and microtubule-based movements; if this is the case, besides m^6A modification [[134]33], m^7G modification is a component in the complex but well-organized multilayered regulatory network that controls the transcriptome of E. tenella. Intriguingly, the DMGs between the two groups were enriched in GO terms of microtubule-based movements; also, significantly enriched KEGG pathways for DMGs included DNA replication and cell cycle. E. tenella genomic DNA double-strand breaks (DSBs) are mainly repaired by NHEJ pathway [[135]37]. In addition, E. tenella also relies on homologous recombination for the repair of genomic DSBs [[136]37]. In the present study, both NHEJ and homologous recombination were enriched for the genes with hypomethylated peaks. This indicated the involvement of internal m^7G methylation in the two pathways for the repair of DSBs in E. tenella. Transcriptome-wide m^6A profiling of sporulated and unsporulated oocysts of E. tenella revealed a positive correlation between m^6A modification levels of most genes and their mRNA expression levels [[137]33]. Herein, the combined analysis of DMGs and DEGs showed that m^7G modification levels of the majority of genes was also positively correlated with their mRNA expression levels. Nevertheless, the presence of the hyper-down and hypo-up genes suggested that mRNA m^7G modification in E. tenella may also negatively regulate gene expression. In the present study, KEGG analysis of the hypo-down genes showed that the significantly enriched pathways included cAMP signaling pathway, autophagy—yeast, and purine metabolism. cAMP signaling pathway and autophagy—yeast were also enriched for the DEGs with altered m^6A-methylated peaks [[138]33]. This indicated that more than one RNA modifications were involved in the regulation of the same signaling pathway in the stage conversion of E. tenella. In addition, purine metabolism is known to provide a cell with the necessary energy [[139]38]. This indicated that m^7G modification was involved in providing cellular energy. The hyper-up genes were associated with RNA transport; spliceosome; and alanine, aspartate, and glutamate metabolism. Alanine, aspartate, and glutamate are three amino acids that can be de novo synthesized by Toxoplasma gondii [[140]39]. It can be inferred that E. tenella may possess the ability to synthesize alanine, aspartate, and glutamate, and m^7G modification was involved in the regulation of amino acid metabolism. In eukaryotes, pre-mRNA splicing catalyzed by the spliceosome is essential for gene expression [[141]40]. Another important step in the control of gene expression is RNA transport [[142]41]. It can be inferred that m^7G modification occurred in more than one steps of gene expression regulation during oocyst development. Conclusions We present here the first transcriptome-wide map of internal mRNA m^7G modification in unsporulated and sporulated oocysts of E. tenella. We found different methylation features between the two stages and analyzed the potential functions of DMGs and DEGs with altered m^7G-methylated peaks. The findings provide a solid foundation for further investigation of the molecular mechanisms governing the development of E. tenella. Supplementary Information [143]13071_2024_6580_MOESM1_ESM.doc^ (32.5KB, doc) Additional file 1: Table S1. Sequences of primers used for MeRIP–qPCR and RT–qPCR analysis. [144]13071_2024_6580_MOESM2_ESM.xlsx^ (10.5KB, xlsx) Additional file 2: Table S2. A statistical summary of the raw and clean reads. [145]13071_2024_6580_MOESM3_ESM.xlsx^ (67.3KB, xlsx) Additional file 3: Table S3. The number of m^7G peaks harbored in each methylated gene. [146]13071_2024_6580_MOESM4_ESM.xlsx^ (661.1KB, xlsx) Additional file 4: Table S4. The differentially methylated transcripts between the two groups. [147]13071_2024_6580_MOESM5_ESM.xlsx^ (257.5KB, xlsx) Additional file 5: Table S5. The genes found to be differentially expressed between the two groups. [148]13071_2024_6580_MOESM6_ESM.xlsx^ (306.2KB, xlsx) Additional file 6: Table S6. The genes with significant changes in both m^7G and mRNA levels. Acknowledgements