The animal study was reviewed and approved by the Animal Ethics Committee of Shandong Agricultural University, China, and performed in accordance with the Committee’s guidelines and regulations (Approval No.: 2004006).

The experimental design and data for this study were obtained from our previous research [18]. The experimental design was as follow: 12 male Landrace piglets weaned at 21 days were collected and divided into the OS group and the NC group. In the OS group, piglets were intraperitoneally injected with DQ (Sigma-Aldrich, Saint Louis, USA) at 10 mg/kg body weight, while the NC group received an equivalent volume of isotonic saline. The trial lasted for seven days. Three piglets were selected from both the OS group and the NC group. After slaughter, the jejunum tissues were collected. Raw data with SRA number PRJNA661634 were downloaded with sratoolkit (v3.0.0). Six samples (three from the OS group and three from the NC group) from this dataset were used for this study. Of these, L01, L02, and L03 were the NC samples, and L04, L05, and L06 were the OS samples.

The raw reads were quality-controlled before read alignment. Fastp (v0.23.2) [19] was used to filter the raw reads further. The reads with adapters, reads with unknown sequences, and low-quality reads were removed for quality control.

Sequencing reads were aligned using HISAT2 (v2.2.1) [20] against the reference genome (Sus_scrofa.Sscrofa11.1.106.chr.gtf). To predict lncRNAs, the reads were then assembled and merged using StringTie (v2.2.1) [21]. Transcript expression levels were calculated for each sample using FeatureCounts (v2.0.3) [22] to obtain the gene counts. To reflect the expression levels of transcripts more realistically, the fragments per kilobase of transcript per million fragments mapped (FPKM) values were used as measures of transcript expression as calculated by Rstudio (v4.1.3).

The lncRNAs were screened according to their characteristics. Briefly, the merged transcripts were compared by GffCompare (v0.12.6) [23] with the Sus scrofa reference genome. i) Transcripts with class_code “i” (fully contained within a reference intron), “u” (intergenic), “x” (exonic overlap on the opposite strand), “o” (other same strand overlap with reference exons), “j” (multi-exon with at least one junction match) were selected. ii) Those transcripts that were more than 200 nucleotides in length and had an exon number greater than or equal to 2 were further considered for the identification of novel lncRNAs. iii) CPC2 (v3.0) [24] and CNCI (v2.0) [25] were used to identify whether the transcript is encoded. When the score <0, the transcript was considered incapable of encoding. The transcripts that were able to be compared to the Pfam [26] database were transcripts with a certain protein domain; they were considered to have coding ability, while the transcripts that were incomparable were considered potential lncRNAs. The criterion for Pfam domain screening was an E-value <1×10⁻⁵. Transcripts with an E-value >1×10⁻⁵ were screened and retained. iv) Those transcripts with FPKM >0.1 were retained for the identification of novel lncRNAs. The identification of novel lncRNAs was considered complete after the above screening. The merged transcripts of the six samples were annotated to the Sus scrofa genome to obtain known lncRNAs.

DESeq2 [27] (v1.34.0) was employed for performing the differential expression analysis based on the counts of genes in samples with biological replicates. Upregulated genes showed higher expression levels in the OS group compared to the NC group in the current study, whereas downregulated genes exhibited lower expression levels in the OS group compared to the NC group. Differentially expressed (DE) lncRNAs were defined as genes with a false discovery rate <0.05 and log2(fold change) ≥1.

Bedtools [28] (v2.30.0) was used to identify cis-target genes for DE lncRNAs. Neighbouring genes within 100,000 bp of DE lncRNAs were considered cis-target genes of the lncRNAs [29]. The Pearson correlation coefficient method [30] was used to predict trans-target genes. When the sample size was greater than or equal to six, the Pearson correlation coefficient method was used to analyse the correlation between lncRNAs and protein-coding genes among samples. Then, lncRNA–mRNA gene pairs with absolute correlation values (|cor|) >0.95 and p<0.05 were retained.

To further understand the function of DE lncRNAs, the gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) databases were used for enrichment analysis of target genes of DE lncRNAs. The GO database contains three main categories of functional information: the biological processes, molecular function, and cellular component categories. Each of these categories contains a hierarchical network of terms that describe different aspects of gene function. KEGG provides pathway maps, gene annotations, and other related information to study the functions of genes within specific pathways. The target genes were then uploaded into the DAVID [31] database 2021. GO terms and KEGG pathways with p<0.05 were considered significant by the DAVID database. To obtain the circos plot, the DAVID results were submitted to SangerBox (3.0).

The pre-mRNAs transcribed from genes can undergo various splicing events, where different exons are selected to generate diverse mature mRNAs. These mature mRNAs are then translated into different proteins, contributing to the diversity of biological traits. ASprofile [32] was used to classify 12 types of alternative splicing (AS) events for each sample.

Six DE lncRNAs were randomly selected to verify whether their expression was consistent with the trends of the RNA-seq results. Total RNA was extracted from piglet jejunal tissues using TRIzol (Invitrogen, Carlsbad, CA, USA). An Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) was employed to assess RNA quality. The cDNA was then synthesized. For the relative quantification of lncRNA, glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was employed as an internal reference gene. Primers of selected lncRNAs (Table 1) were designed and then synthesized by Accurate Biotechnology Co., Let, Changsha, China. The quantitative real-time polymerase chain reaction (qRT-PCR) system for lncRNA had a total volume of 20 μL. The system consisted of 10 μL of 2X SYBR Green Pro Taq HS Premix, 0.4 μL of forward primer, 0.4 μL of reverse primer, 2 μL of cDNA template, and 7.2 μL of RNase-free water. Run the formulated 20 μL system on a Light Cycler 96 real-time PCR system (Roche, Basel, Switzerland) with the following program: 95°C for 30 s followed by 40 cycles of 95°C for 5 s and 60°C for 30 s. The 2^−ΔΔCT relative quantification approach [33] was utilized to perform the quantitative analysis of the data.

The GC content of each sample was greater than 50.28%, Q20% was above 97.71%, Q30% was above 93.50%, and N% was below 0.01% after quality control (Table 2). This indicated that the clean reads obtained were suitable for subsequent transcriptome analysis. In total, 95.33% to 96.67% of the clean reads of individual samples were mapped to the Sus scrofa reference genome. A total of 84.05% to 87.18% of the mapped reads were considered unique mapped reads, while only 4.96% to 7.18% of the mapped reads were multiple mapped reads (Table 3). The mapping results indicated that the quality and reliability of the sequencing data were relatively high.

Overall, 78,398 transcripts were used as input. A total of 790 candidate novel lncRNAs were finally identified (Figure 1A). A total of 13,053 transcripts were kept by filtering and retaining transcripts with class_codes “i”, “u”, “x”, “o”, “j”. A total of 83.35% of the initial transcripts were filtered; this was the step that filters out the most transcripts in identifying novel lncRNAs. After removing transcripts that were less than 200 nt in length and had an exon numbers less than 2, 12,702 transcripts remained. Only 0.48% of transcripts were further filtered out in this process. A total of 2,486 transcripts were retained through CNCI filtering; 2,162 transcripts were retained through CPC2 filtering; and 2,138 transcripts were retained through Pfam protein structural domain analysis. The results of the three software programs were intersected, and 1,179 transcripts were preserved (Figure 1B). A total of 9.07% of transcripts were further filtered out in this process. Transcripts with low expression were then removed. A total of 0.49% of transcripts were further filtered out in this process. Finally, 790 transcripts were retained. The 790 transcripts obtained were considered to be novel lncRNAs. The merged transcripts were compared with the Sus scrofa reference genome, and the transcripts with low expression were removed, resulting in 5,289 known lncRNAs. Ultimately, a total of 6,079 lncRNAs were identified.

The exon number distribution statistics of lncRNAs showed that the highest number of exons was 1, 2 (Figure 1C). Most lncRNAs were distributed on chromosome 1 in jejunal tissue, followed by chromosome 6 and chromosome 13 (Figure 1D). The length distribution statistics of lncRNAs showed that the most common length bin of the lncRNAs in the samples was 3,001 nt to 6,000 nt, followed by 200 nt to 3,000 nt and 6,001 nt to 9,000 nt (Figure 1E).

After differential expression analysis in the OS group and NC group, 79 DE lncRNAs (34 lncRNAs were upregulated and 45 lncRNAs were downregulated) were screened (Figure 2A; Supplementary Table S1). The heatmap showed that the expression levels of DE lncRNAs (Supplementary Table S2) were significantly different between the OS group and the NC group (Figure 2B).

Of the 79 DE lncRNAs, 61 lncRNAs had cis-target genes within 100,000 bp upstream and downstream, in a total of 135 targets. A total of 12,980 lncRNA–mRNA pairs were predicted to be associated by the Pearson correlation coefficient method. The top 50 lncRNA-mRNA gene pairs (Supplementary Table S3) with the most significant p-values were shown in Figure 2(C). In this figure, the DE lncRNAs point towards their target mRNAs, representing the differential expression of lncRNAs that regulated their target genes in trans.

In the present study, the target genes of DE lncRNAs were significantly enriched in 68 BP, 40 CC, and 48 MF categories (Supplementary Table S4). The top 20 BP terms mainly included B-cell activation, cargo loading into COPII-coated vesicle, cell cycle, DNA recombination, regulation of growth, tissue development, and other terms (Figure 3A). The top 20 CC terms mainly included nucleus, cytosol, nucleoplasm, cytoplasm, kinetochore, chromatin, and other terms (Figure 3B). The top 20 MF mainly included RNA polymerase II core promoter proximal region sequence-specific DNA binding, chromatin binding, methyl-CpG binding, protein kinase A binding, and other terms (Figure 3C).

A total of 4,645 target genes of DE lncRNAs were significantly enriched in 96 signaling pathways (Supplementary Table S5). The top 20 KEGG pathways mainly included metabolic pathways, cell cycle, chemical carcinogenesis-receptor activation, mismatch repair, DNA replication, chemical carcinogenesis -ROS, and other pathways (Figure 3D). Among the 96 KEGG pathways associated with OS were chemical carcinogenesis-ROS, colorectal cancer, Parkinson disease, insulin pathway, cellular senescence, Alzheimer’s disease, and other pathways.

AS is an important mechanism of gene expression regulation, playing a crucial role in the normal functioning and disease development of organisms. In this research, a total of 12 AS events were detected (Figure 4). In this study, it was found alternative 5’ first exon (TSS) and alternative 3’ last exon (TTS) accounted for 67.84% to 72.13% of the AS events. This suggests that TSS and TTS were the two most prevalent AS events observed in this research.

To validate the accuracy of RNA-seq, six DE lncRNAs (three upregulated lncRNAs and three downregulated lncRNAs) were randomly selected to perform qRT-PCR (Figure 5A). The correlation and p-value were obtained after linear fitting of the log₂ (fold change). The results showed that there was a significant correlation between the RNA-seq and qRT-PCR (r = 0.87; p = 0.03 for DE lncRNAs) (Figure 5B). This also indicated the accuracy of high-throughput sequencing results.