The animal study was reviewed and approved by the Animal Ethics Committee of Chungnam National University (202012A -CNU-164).

Figure 1 summarizes the experimental design. This study used 39 1-day-old male broiler chicks (Indian River; Aviagen, AL, USA). After ensuring that the chicks were free from oocysts, we randomly allocated them to 13 cages. Following a 14-day acclimation period, 21 chickens, referred to as the positive control (PC), were orally administered 1 mL of Livacox T (Biopharm, Prague, Czech Republic), which is a live attenuated vaccine. Each chicken in the PC group received a 10-fold increased dose containing 3 to 5×10³ active oocysts of Eimeria acervulina, Eimeria tenella, and Eimeria maxima for infection. The 18 chickens in the non-challenged negative control (NC) were inoculated with distilled water to simulate inoculation stress.

At 4, 7, and 21 days post-infection (dpi), birds in the PC and NC treatments were humanely euthanized. Cecum tissue samples were collected from all chickens and used for RNA extraction. At each time point, three samples from the NC group and five samples from the PC group were randomly selected for RNA sequencing.

Total RNA was extracted from cecal tissue samples using TRIzol reagent (Invitrogen, Carlsbad, CA, USA). The TruSeq Stranded mRNA sample preparation kit (Illumina, San Diego, CA, USA) was used to construct a cDNA library. Following library construction, all samples were sequenced on the Illumina NovaSeq 6000 platform, which generated 101 bp paired-end reads.

Adaptor-trimmed sequence data, which we described previously [11], were mapped to the chicken reference genome (GRCg7b, GCA_016699485.1). The reference genome was constructed using bowtie2 v2.5.1 and SAMtools v1.17. STAR v2.7.10b was used for sequence mapping to the reference genome and counting mapped reads. The gene annotation version used was GRCg7b (National Center for Biotechnology Information Annotation Release 106). The EdgeR v3.42.4 package in R [19] was used to quantify gene expression. To avoid statistical bias, only genes with a total read count >8 were used for estimating expression levels. Then, the read count was normalized using the trimmed mean of M-value method. Expression levels were compared between the PC and NC groups at every time point. The p-values of the differential expression tests were adjusted using the false discovery rate (FDR) method [20]. Then, genes located on chicken chromosome 16 were selected for the next step. Significant DEGs were defined based on the following criteria: FDR<0.05 and absolute fold change ≥1.5. A lower level of significance for differential expression was also identified within the range 0.05≤FDR<0.1 and absolute fold change ≥1.5. The final DEGs were identified at each time point. Volcano plots depicting the DEGs were generated using the ggplot2 package in R [21]. To determine the functions of the group of DEGs, gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) analyses were performed using the Database for Annotation, Visualization, and Integrated Discovery (DAVID). Gene functions were also searched in the relevant literature.

Sequencing of the 24 samples yielded an average of 40,967,230 raw reads per sample after combining forward and reverse sequences. After removing sequencing adapters and trimming low-quality reads, at least 98.29% of the reads remained in all datasets. On average, 40,455,803 reads were retained for subsequent mapping. We aligned the reads to the chicken reference genome, and successfully assembled an average of 37,600,629 reads (approximately 92.95% of the total mapped reads). Table 1 summarizes the sequence data processing.

The chicken chromosome 16 reference genome contained 309 annotated genes. Figure 2 and Table 2 present the findings for 4 dpi, in which five significant DEGs were identified as having FDR<0.05. Of these DEGs, three genes were upregulated in Eimeria-infected chickens, whereas two genes were downregulated. Two additional upregulated genes were identified among DEGs with a lower significance threshold (0.05≤FDR<0.1). Upon comparing the NC and PC samples at 7 dpi, we found four significant DEGs having FDR<0.05. Of these genes, three genes were upregulated and only one gene was downregulated in Eimeria-infected chickens. Also, there were two additional upregulated genes and one downregulated gene showing a lower significance level (0.05≤FDR <0.1). Notably, no differential expression was observed in any gene at 21 dpi. Therefore, the gene expression comparison results for 21 dpi were excluded from subsequent analyses.

In Eimeria-infected chickens, five genes had significantly higher expression at 4 dpi (FDR<0.1). These genes included MHCY36 and LOC112529954, which are located in the MHC-Y region of the chicken genome, and DMB1, TAP1, and TAPBP, which are located in the MHC-B region. Of the downregulated genes, two (LOC124417245 and LOC124417241) were in the non-MHC region. Notably, at 4 dpi, the most significant DEG was the TAP1 gene and the largest fold change was for the MHCY36 gene.

At 7 dpi, five genes were significantly upregulated (FDR<0.1): the TRIM39.2 gene in the MHC-B region and LOC121106941, LOC121106957, and two ribosomal RNA genes in the non-MHC region. Two genes (MHCYL and LOC124417244) downregulated in infected chickens were classified in the MHC-Y class I gene cluster. Among the upregulated genes, the TRIM39.2 gene had the most significantly upregulated differential expression (FDR = 0.0046). The differential gene expression results are provided in Table 2 and Figure 2.

GO analysis conducted on the DEGs at 4 and 7 dpi yielded no significant GO terms. However, KEGG analysis revealed one significant biological pathway (p = 0.03) at 4 dpi, herpes simplex virus 1 infection (gga05168) pathway, which includes the TAP1 and TAPBP genes.