Twenty healthy full-sib chicken families were randomly selected from Jiangsu Jinghai Poultry Industry Group Co., Ltd., and each family contained one male and four females. Each family was raised separately, and offspring were produced through artificial insemination (the number of F1 generations in each family was ≥ 10). The F1 generation chicks of each family were housed individually in sterile animal rooms and fed antibiotic-free feed and drinking water. At 18 days of age, the F1 individuals of each family were randomly divided into two groups (the challenge group and the control group). Each bird in the challenge group was orally challenged with 3.5×10⁴ E. tenella sporulated oocysts, and the control group was fed the same amount of PBS solution. According to the description in a previous study (14), on day 4.5 post-infection, the most resistant and most susceptible families were selected from 20 families by clinical symptoms (emaciation, drooping wings, and dying death), fecal scores and cecal lesion scores, respectively.

Coccidia oocysts were a gift from the Department of Parasites, College of Veterinary Medicine, Yangzhou University, and were expanded and rejuvenated in healthy 15-day-old chicks. The oocyst collection, purification and sporulation processes were described by Jiao et al. (15).

On day 4.5 post-infection, 3 individuals from each group were randomly selected to be euthanized by cervical dislocation. The cecal tissues were collected and washed with sterilized PBS to remove the intestinal contents and mucus. Then, the treated cecal tissues were immediately stored in liquid nitrogen for future use. All the animal experimentation protocols were approved by the Animal Welfare Committee of Yangzhou University (license number: SYXK (Su) IACUC 2012-0029).

The cecal tissues (JC, JS, and JR groups) were cut into approximately 0.5 cm pieces and fixed with 4% tissue fixative. After embedding in paraffin and sectioning (5 µm), the sections were stained with hematoxylin-eosin. Then, the slices were observed and imaged under a microscope (Nikon, Japan)

Total RNA was extracted from cecum tissues with TRIzol reagent (Invitrogen, CA, USA) following the manufacturer’s recommended procedure. RNA purity and concentration were evaluated using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA), and RNA integrity was assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). If the quality of the sample RNA met the experimental requirements (RNA integrity number>=7 and 28S/18S>=0.7), then subsequent sequencing was performed. Libraries were constructed with 1 µg of RNA from each sample using the TruSeq Stranded Total RNA with Ribo-Zero Gold Kit (Illumina, Cat. RS-122-2301) according to the manufacturer’s protocol. High-throughput sequencing, including lncRNA sequencing and mRNA sequencing, was carried out on an Illumina 2500 platform, and 150 bp paired-end reads were generated (OE Biotech, Shanghai, China).

The transcriptome from each dataset was assembled independently and then merged to generate a final transcriptome using the Cufflinks 2.0 program (16). All the transcripts that overlapped with known mRNAs, other noncoding RNAs and non-lncRNAs were removed. The obtained transcripts were screened according to a length greater than 200 bp and a number of exons >2. Next, CPC (v 0.9-r2) (17), PLEK (v 1.2) (18), CNCI (v 1.0) (19), and Pfam(v 30) (20) were used to predict and analyze the transcripts with coding potential. The novel predicted lncRNAs were obtained through these processes. The characteristics (including length, type, and number of exons) of the lncRNAs were analyzed after screening.

Sequencing reads were mapped to the reference genome using HISAT2 (21). The novel predicted lncRNAs and known lncRNAs (from the NCBI and Ensembl databases) were subjected to differential expression analysis, and the expression of lncRNAs was calculated using fragments per kilobase of transcript per million mapped reads (FPKM) (22). The number of counts of each sample lncRNA was standardized using DESeq (23), and the negative binomial distribution test (NB) method was used to test the significance of the difference in the number of reads. For mRNA, HTseq-count (24) software was used to obtain the number of reads aligned to the protein-coding gene in each sample, and cufflinks was used to calculate the FPKM value of each gene. Finally, a fold change (FC) > 1 and a P value < 0.05 were set as the thresholds for significantly differential expression of lncRNAs and mRNAs.

To explore the functions of these differentially expressed lncRNAs during E. tenella infection, Gene Ontology (GO) enrichment analysis of the three categories of biological process (BP), cellular component (CP), and molecular function (MF) was performed on neighboring genes of the identified lncRNAs. Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were used to reveal enriched pathways of the differentially expressed neighboring genes. A P value ≤0.05 indicated significantly enriched GO terms and KEGG pathways.

The Pearson correlation coefficient (PCC) was used to analyze the correlation between differentially expressed lncRNAs and mRNAs and to construct a coexpression network. Coexpression relationship pairs of lncRNA-mRNA were selected according to a PCC value ≥0.85 and a P value <0.05, and a network map was then constructed using the R network package.

Studies have shown that lncRNAs regulate the expression of protein-coding genes through cis and trans methods (25). For cis-regulation prediction, the coding genes within 100 kb upstream and downstream of the differentially expressed lncRNAs and had significant co-expression with the lncRNAs were identified (P ≤ 0.05). For trans-regulating prediction, the following procedures were used: (1) based on the coexpression results, we screened out lncRNAs and genes that were not on the same chromosome as the candidate targets, and (2) the number of bases for the two nucleic acid molecules of the lncRNA and gene to directly interact was ≥ 10, and the base binding free energy was ≤ -50. Then, the lncRNA cis- and trans-regulation networks were constructed using Cytoscape software.

DF-1 cells (a chicken embryo fibroblast cell line) were cultured at 37°C and 5% CO2 in basal medium (Gibco, USA) containing 10% fetal bovine serum (Gibco, USA) and 1% penicillin and streptomycin.

For sporozoite exposure, E. tenella sporozoites were separated from sporulated oocysts through the G3 funnel method. To ensure the successful construction of the in vitro cell coccidia model, the ratio of the number of sporozoites to chicken DF-1 cells was 3:1.

Three lncRNA BTN3A2-specific small interfering RNA (siRNA) and negative controls sequences ( Table 1 ) were designed and synthesized by GenePharma Biotechnology (Shanghai, China). For overexpression, the lncRNA mature sequence ( Supplementary Text 1 ) was inserted into the pcDNA3.1(+) cloning vector to construct the overexpression plasmid: pcDNA3.1-lncRNA BTN3A2.

Chicken DF-1 cells were inoculated in a 12-well plate and grown to 70%. Then, the cells were subsequently transfected with 20 µmol/L siRNA or negative control, 0.8 µg/µl pcDNA3.1-lncRNA BTN3A2 and pcDNA3.1 using jetPRIME Transfection Reagent (Polyplus, Illkirch, France) for 24 h and then stimulated with sporozoites for another 6 h. The cells were harvested for RNA extraction, and the suspension was collected for ELISA.

Chicken cecal tissues from the JC, JR, and JS groups on day 4.5 post-infection were collected. Total RNA of tissue samples and DF-1 cells was isolated using TRIzol reagent for cDNA synthesis. Six differentially expressed lncRNAs in this study were selected for random verification of the RNA-Seq results, and β-actin was used as an internal reference gene. Primer sequences including chicken IL-6, TNF-α, IL-1β and IL-8 ( Table 2 ) were synthesized by Sangon Biotech Co., Ltd. (Shanghai, China). Then, qRT–PCR was performed using the ChamQ SYBR qPCR Master Mix Kit (Vazyme Biotech, Nanjing) according to the manufacturer’s specifications. The PCR program was run at 95°C for 30 s followed by 40 cycles of 95°C for 10 s and 60°C for 30 s. Each sample was analyzed in triplicate, and data analyses were performed using the 2−ΔΔCt method (26).

Chicken IL-6, TNF-α, IL-1β and IL-8 in the supernatant of DF-1 cells were measured by enzyme-linked immunosorbent assay (ELISA) according to the protocols of the kits (R&D Systems, USA).

All the data are presented as the mean ± SEM, and one-way ANOVA via SPSS 25.0 software was used to compare the differences in data in the different groups. A P value < 0.05 was considered statistically significant.

Images of chicken cecal tissues at 4.5 days post-infection in the JC, JR and JS groups are shown in Figure 1 . Bleeding spots appeared on the mucosal surface of the cecal tissue in the JR group. The cecum wall was slightly swollen and thickened, and the contents were relatively normal ( Figure 1B ). However, the chicken cecal tissues in the JS group had severe intestinal bleeding, and the intestinal wall was significantly swollen ( Figure 1C ). The cecal tissues of the control group are shown in Figure 1A . In addition, the pathological section analysis of the intestinal tissue in this study showed that the structure of each layer of the cecal tissue in the JC group was clear, the mucosal epithelium was intact, and the epithelial cells were tightly arranged ( Figure 1D ). The cecal tissue of the JR group chickens exhibited more mucosal necrosis ( Figure 1E , black arrows) with moderate infiltration of lymphocytes and heterophils ( Figure 1E , blue arrows). The cecal tissue in the JS group showed extensive mucosal necrosis and ulceration ( Figure 1F , black arrows), a large number of lymphocytes ( Figure 1F , blue arrows) and massive multifocal bleeding ( Figure 1F , red arrows). Therefore, these data indicate that the E. tenella infection model in chickens of resistant and susceptible families was successfully constructed.

After clean reads were compared to the chicken reference genome, a total of 13,809 lncRNAs were identified in the sequencing results from 9 samples. According to the source classification, the proportion of identified lncRNAs included exonic lncRNAs, intronic lncRNAs, intergenic lncRNAs and antisense lncRNAs as shown in Figure 2A . Moreover, the length distribution of the identified lncRNAs were shown in Figure 2B , and most of the lncRNAs were over 2000 bp in length. LncRNAs were distributed on every chromosome of chickens, and a large number of lncRNAs were distributed on chromosome 1, followed by chromosome 2 ( Figure 2C ).

To screen out the key lncRNAs during E. tenella infection, we analyzed the differential expression of the lncRNAs based on an FC > 1 and P value <0.05. A total of 192 lncRNAs were differentially expressed between the JS and JR groups, of which 114 were upregulated and 78 were downregulated; 263 differentially expressed lncRNAs were identified in the JS and JC groups, namely, 121 upregulated and 142 downregulated lncRNAs. Additionally, in the JR and JC groups, 109 differentially expressed lncRNAs were found, with 62 upregulated and 47 downregulated ( Table 3 ; Supplementary Table 1 ; Figures 3A–C ). Hierarchical clustering analysis revealed that the differentially expressed lncRNAs exhibited different expression patterns among the control group and the treatment groups ( Figures 3D–F ).

To predict the potential functions of these differentially expressed lncRNAs, GO terms and KEGG pathways enrichment were analyzed according to the neighboring genes of the differentially expressed lncRNAs. The GO analysis results showed that the top 30 terms enriched in the JS vs JC groups mainly included bacterial defense response, circulating immune complexes and antigen binding. Similarly, the immune response, B cell receptor signaling pathway and immunoglobulin receptor binding were significantly enriched in the JS vs JR groups, while signaling receptor binding and transcription factor binding were significantly enriched in the JR vs JC groups ( Figure 4 ; Supplementary Table 2 ).

The top 20 enriched KEGG pathways in the three groups that may participate in the coccidial immune response process mainly included the NF-kappa B signaling pathway, B cell receptor signaling pathway and natural killer cell mediated cytotoxicity in the JS vs JC and JR vs JS groups and the MAPK signaling pathway in the JR vs JC groups ( Figure 5 ; Supplementary Table 3 ).

To further reveal the function of lncRNAs during E. tenella infection, coexpression networks of the differentially expressed lncRNAs and mRNAs in the JS, JR and JC groups were constructed ( Figure 6 ; Supplementary Table 4 ). In the JS vs JC groups, TCONS_00025088 had a significant negative coexpression relationship with the immune-related genes LRRC19 and F2RL1. TCONS_00009847 was positively coexpressed with SPP1 and RORC, and TCONS_00036436 was negatively coexpressed with F2RL2 and BTN3A3. Moreover, XR_003072630.1 and TCONS_00025084 were coexpressed with DTX1 and FABP1 in the JR vs JS groups, and ENSGALT00000093041 was coexpressed with IL22RA2 in the JR vs JC groups.

For cis-targeted gene prediction, 20 mRNAs were cis-regulated by 24 lncRNAs in the JS vs JC groups, while 13 mRNAs were regulated by 12 lncRNAs in the JS vs JR groups, 5 mRNAs were regulated by 8 lncRNAs in the JR vs JC groups ( Supplementary Figure S1 ). Interestingly, we found that BTN3A2, an important gene involved in immune and inflammatory responses, was a potential cis-target gene of ENSGALT00000093710 (termed lncRNA BTN3A2), and there was also a significant co-expression relationship (R=0.965, P <0.001). Moreover, this lncRNA was significantly increased in the cecum tissues after E. tenella infection ( Figure 7C , P <0.05). Therefore, based on these evidences, we selected lncRNABTN3A2 for further functional experiments. In addition, the trans-targeted genes of the differentially expressed lncRNAs were also identified in Supplementary Table 5 .

To verify the accuracy of the high-throughput sequencing results, we randomly selected six differentially expressed lncRNAs for qRT-PCR verification. The qRT-PCR results are consistent with the transcriptome sequencing results, which verifies the sequencing results of this study ( Figure 7 ).

After coccidia infection or sporozoite stimulation, lncRNA BTN3A2 was significantly increased in chicken cecal tissue and DF-1 cells ( Figures 7C , 8A ). This suggests that lncRNA BTN3A2 may be involved in the immune response to coccidial infection. Then, an in vitro model of E. tenella infection and the lncRNA BTN3A2 overexpression vector (pcDNA3.1-lncRNA BTN3A2) was constructed in this study ( Figure 8B ). The image of the fluorescent plasmid vector after transfection for 24 h showed that the pcDNA3.1 plasmid was successfully transferred into DF-1 cells ( Figure 8C ). As respect, the expression of lncRNA BTN3A2 in DF-1 cells was significantly increased when the cells were transfected with pcDNA3.1-lncRNA BTN3A2 compared to the NC group ( Figure 8D ). DF-1 cells were then stimulated with sporozoites for another 6 h, and the expression and protein levels of four inflammatory cytokines were detected. The results showed that sporozoite stimulation significantly increased the expression of IL-1β, IL-6, TNF-α and IL-8 compared with the no-stimulation group. Moreover, the pcDNA3.1-lncRNA BTN3A2 further inhibited mRNA expression ( Figures 8E–H ) and the release ( Figures 8I–L ) of inflammatory cytokines, including IL-1β, IL-6, TNF-α and IL-8, in sporozoite-stimulated DF-1 cells.

To further confirm the effect of lncRNA BTN3A2 on the regulation of the inflammatory response, we silenced lncRNA BTN3A2 and then detected the release of inflammatory cytokines in sporozoite-stimulated DF-1 cells. After transfection with three interfering sequences of lncRNA BTN3A2 into DF-1 cells, siRNA-lncRNA BTN3A2-1000 had the highest interfering efficiency and significantly reduced the expression of lncRNA BTN3A2 ( Figure 9A ). DF-1 cells were then transfected with siRNA-lncRNA BTN3A2 to examine the production of inflammatory cytokines in sporozoite-stimulated DF-1 cells. The data showed that knocked down of lncRNA BTN3A2 further promoted the production of IL-6, IL-1β, TNF-α and IL-8 cytokines compared with siR-NC group ( Figures 9B–I ), which produced an opposite trend to that of lncRNA BTN3A2 overexpression. Taken together, these results further demonstrated that lncRNA BTN3A2 inhibits the host immune response against E. tenella infection, thereby avoiding excessive inflammation and immune responses.