Changes in long non-coding RNA expression profiles related to the antagonistic effects of Escherichia coli F17 on lamb spleens

Sheep colibacillosis is one of the most common bacterial diseases found at large-scale sheep farms. The aim of this study was to employ RNA-seq to screen differentially expressed (DE) long non-coding RNAs (lncRNAs) that impart antagonistic or sensitive effects on Escherichia coli F17. In this study, individuals who had antagonistic or sensitive responses to E. coli F17 were identified by feeding E. coli F17 strains to Hu lambs. The sensitive group had higher levels of intestinal bacteria than that in the antagonistic group (P < 0.05), the jejunum showed various levels of mucosal tissue damage and had a dark colour, and disintegration of part of the small intestinal villi was observed. Totals of 34 DE lncRNAs and 703 DE mRNAs in two groups were identified. qRT-PCR results for 12 randomly selected DE lncRNAs and DE mRNAs were consistent with the RNA-seq data. Gene Ontology (GO), KEGG Pathway enrichment and lncRNA-mRNA interaction analyses identified 6 co-expressed genes, namely, MYO1G, TIMM29, CARM1, ADGRB1, SEPT4, and DESI2. This is the first study that has performed expression profiling of lncRNAs in the spleen of antagonistic and sensitive lambs. The identification of DE lncRNAs can facilitate investigations into the molecular mechanism underlying resistance to diarrhoea in sheep.

Sheep colibacillosis is one of the most common bacterial diseases found at large-scale sheep farms. The traditional method of controlling the bacterial disease is by antibiotic therapy, although this approach also has several disadvantages. Detecting the expression of antagonistic genes in sheep colibacillosis provides information that may facilitate the elucidation of the molecular mechanism underlying disease resistance to Escherichia coli. In 1966, Orskov et al. 1 first reported the porcine E. coli adhesion antigen K88, which is an episome-determined antigen 2 . In addition, the morphology of the pilus situated on the surface of the bacteria has been examined by electron microscopy. To date, numerous animal-derived enterotoxigenic Escherichia coli (ETEC) pili have been identified, including K88, K99, 987P, F17 and F41, which are all vital virulence factors. Long non-coding RNAs (lncRNAs) are a class of non-coding RNAs longer than 200 nucleotides. Several studies have revealed that lncR-NAs are closely related to the development of human tumours, cardiovascular diseases, and metabolic diseases. Furthermore, there is growing evidence that lncRNAs play a significant regulatory role in anti-viral and other natural immune responses [3][4][5][6][7] . However, research on the function of lncRNAs has mainly focused on the regulation of muscle growth, testicular development, hair follicle development, and other traits in sheep, which are less well studied compared to humans [8][9][10] . A few studies have recently focused on disease prevention and control in sheep 11,12 , including disease resistance, yet our understanding of its underlying molecular mechanism is limited. In the present study, the expression levels of lncRNAs in two sheep spleen phenotypes, antagonistic or sensitive to pili of E. coli F17, were assessed by RNA-seq, and key lncRNAs were investigated by target prediction and functional annotation analysis based on the cis mechanism, a transcription activation and expression control method for adjacent mRNAs by non-coding RNAs. The results were further verified by q-PCR. Our results provide a

Results
Histological observation and comparison of the number of intestinal bacteria in the antagonistic and sensitive groups. Based on faecal morphology 13 , the experimental subjects were divided into the antagonistic group (12p, 13p, 14p) and the sensitive group (15p, 16p, 17p). In the sensitive groups, the bacterial count was arranged from 4.7*10 8 to 1.9*10 9 , the mean of the bacterial count was 1.22 *10 9 , and yet the count dropped 5.1*10 6 to 9.0*10 7 in the antagonistic groups, the mean of the bacterial count was 3.37*10 7 (Table 1). Compared to the antagonistic group, the number of intestinal bacteria in the sensitive group was significantly higher (P < 0.05, Fig. 1). The jejunum mucosal tissues of lambs in the sensitive group were damaged, dull in colour, lysed, and exhibited large lacuna that could be observed at the submucosal layer. The intestinal villi had disintegrated and were highly vascularized (i.e., capillaries), and the intestinal mucosa showed severe damage, thus making it difficult to prepare histological sections for assessment (Fig. 2).
Summary of RNA sequencing of spleens in sheep. cDNA libraries of the lamb spleens from the antagonistic and sensitive groups were constructed. Sequencing was performed using the Illumina HiSeq 2500 platform. The antagonistic and sensitive groups generated a total of 354,943,820 and 370,616,990 raw reads, with GC contents of 48.33% and 49.67%, respectively. The valid reads in the clean reads were mapped to the O. aries v4.0 reference genome, and more than 73.5% of the reads were mapped to the genome. The reads mapped to multiple locations of the reference sequence were less than 4.5% and more than 70% of the reads were uniquely mapped to the reference sequence. Approximately 35% of the reads mapped to the positive and negative chains of the genome. In addition, the number of reads that were mapped to the exonic regions (~60%) was higher than those that were mapped to intergenic and intron regions by annotation analysis. These results indicate that the matching efficiency of our de novo assembly is high, and most reads mapped to the exonic region ( Table 2).

Identification of transcripts in sheep spleens.
After mapping the reference sequence, we identified 1,988 lncRNAs and 38,843 mRNAs from 42,460 compiled transcripts. The length of the lncRNAs was mainly distributed within the range of 200 bp-5,000 bp (Fig. 3a), and the average length was 2,124 bp. Additionally, the  Table 1. Comparison of the number of intestinal bacteria between the antagonistic and sensitive lambs. Note: "ng" means "no growth".

Analysis and validation of DE transcripts.
The expression levels of lncRNA and mRNA transcripts were estimated using the FPKM values. We found that the expression level of the lncRNA transcripts was relatively low (Fig. 4). A total of 14 upregulated and 20 downregulated DE lncRNAs and 370 upregulated and 333 downregulated DE mRNAs were screened under conditions of P < 0.05 and |log 2 (fold change)| > 1 (Fig. 5). To further verify the reliability of our RNA-seq data, a total of 12 DE lncRNAs and DE mRNAs were randomly selected. Their relative expression levels in the antagonistic and sensitive lambs were confirmed by q-PCR ( Fig. 6) and were found to coincide with our RNA-seq results (Fig. 7), thus indicating that our RNA-seq data were reliable. Our analyses also showed that high-throughput sequencing has the advantage of detecting genes that are expressed at relatively very low levels (0 < FPKM < 1). Comparison between the DE lncRNAs and the entries in the KEGG Pathway database indicated that a total of 34 lncRNAs could be annotated and classified into 149 KEGG pathways. The number of DE lncRNAs in the top 30 KEGG Pathways is shown in Fig. 8. The number of DE lncRNAs in the thyroid hormone signalling pathway (path:ko04919), spliceosome (path:ko03040), leukocyte transendothelial migration (path:ko04670), neurotrophin signalling pathway (path:ko04722), lysosome (path:ko04142), MAPK signalling pathway-yeast (path:ko04011), sphingolipid signalling pathway (path:ko04071), phagosome (path:ko04145), and oxidative phosphorylation (path:ko00190) categories were relatively higher.

LncRNAs and their adjacent coding genes.
We searched for all the coding genes within the 100-kb flanking regions of the lncRNAs, and genes that were significantly co-expressed with the lncRNAs as indicated by Pearson correlation calculations were identified. These co-expressed genes that were adjacent to the lncRNAs were presumed to be regulated by lncRNAs. Therefore, we identified 6 genes may be regulated by their associated lncRNAs (Table 3).

Discussion
Due to the rapid development of transcriptome analysis, lncRNAs have recently been considered as a novel modulator of cell development 14 . The lncRNAs that we have identified in this study are mainly associated with cancer, such as prostate cancer 15 , gastric cancer 16 , lung cancer 17 and breast cancer 18 , as well as reproduction [19][20][21][22] . However, investigations of lamb diarrhoea in relation to lncRNAs are limited. The Hu sheep is a unique breed with high fecundity and strong adaptability to warm-wet climates and can be kept indoors all year. This study has provided the first overview of lncRNAs in relation to diarrhoea in sheep, as well an investigation into their possible roles in disease resistance. Major economic losses in sheep farms are often due to diarrhoea. In this study, We mainly study the immune status of lambs, and the spleen is the largest immune organ; thus, we selected the spleen as the research object. we found that the expression level of lncRNAs was lower than mRNAs in lamb spleens (Fig. 4), which was in agreement with the results from sheep testicular tissues 9 , and the average lengths of lncRNAs and mRNAs in sheep were longer than those in pigs (1,713 bp and 1,983 bp, respectively) 20 . We searched for all the coding genes within the 100-kb flanking regions of the lncRNAs and identified intersecting genes that were significantly co-expressed with the lncRNAs as indicated by Pearson correlation calculations 23 . We also identified the following 6 genes as being co-expressed with the lncRNAs: myosin IG (MYO1G), translocase of inner mitochondrial membrane 29 (TIMM29), co-activator associated arginine methyltransferase 1 (CARM1), adhesion G protein-coupled receptor B1 (ADGRB1), septin 4 (SEPT4), and desumoylating isopeptidase 2 (DESI2).
MYO1G plays an important role in maintaining cell stiffness in B-cell lymphocytes. The deletion of the myo1g gene results in a reduction in cell stiffness, which in turn affects cell adhesion, proliferation, phagocytosis, and endocytosis in B-cell lymphocytes 24 . Investigations of TIMM29 are limited. TIMM29 has been identified as the first specific component of the mammalian TIMM22 protein complex and plays an important role in the assembly   of the TIMM23 protein 25,26 . CARM1, a member of the protein arginine methyltransferase (PRMT) family, is an enzyme with a highly conserved domain with methyltransferase activity. CARM1 knockout mice die at birth 27 , indicating that CARM1 is essential to postnatal survival. It was later discovered that CARM1 inhibition promotes HIV-1 activation 28 . ADGRB1 is a member of the transmembrane protein-adhesion G protein coupled receptor (aGPCR) family, which is characterized by a conserved GAIN domain that has autologous proteolytic activity that can cleave the receptor near the first transmembrane domain. Studies have shown that the new N-terminal  stalk, which is revealed by GAIN domain cleavage, can directly activate aGPCRs as a tethered agonist 29 . Septins are a highly conserved cytoskeletal family with GTPase activity. The tumour suppressor SEPT4 is a member of the septin family that can induce cancer cell apoptosis 30 . Mutations in the SEPT4 gene in mice can lead to disorders involving the annuli of spermatozoa 31 and adjacent cortical structures, thereby causing low sperm motility, ultimately leading to infertility 32 . DESI2 is a pro-apoptotic gene; in vitro experiments have shown that its overexpression induces apoptosis in pancreatic cancer and other tumour cells, which can effectively inhibit the proliferation of some cancer cells. Gene therapy using DESI2 and IP10 significantly inhibits tumour growth and effectively prolongs the survival of tumour-bearing mice 33,34 .   . Relative expression levels of DE lncRNAs and mRNAs between antagonistic and sensitive lambs were confirmed by q-PCR. Note: "**" means highly significant correlation; "*" means significant correlation; "ns" or "no SuperiorScript" means no significant correlation. The same as below. A total of 703 mRNAs and 34 known lncRNAs were differentially expressed between the antagonistic and sensitive groups, including 14 upregulated lncRNAs and 20 downregulated lncRNAs. In addition, the present study identified 1,942 novel lncRNAs. We searched for all the coding genes within the 100-kb flanking regions of the lncRNAs and identified genes shared between the two experimental groups that were significantly co-expressed with lncRNAs as indicated by Pearson correlation calculations. We have determined that 6 genes may be regulated by their associated lncRNAs. To validate our RNA-seq results, q-PCR was performed to verify the expression levels of the 12 known lncRNAs and mRNAs. The final results coincided with our RNA-seq data.
GO is a bioinformatics tool that has been extensively utilized in studying the relationship of various gene functions. GO analysis indicated that 16 out the 34 DE lncRNAs were enriched with the protein binding (GO: 0005515) category. Moreover, KEGG Pathway analysis showed that the sphingolipid signalling (path: ko04071), axon guidance (path: ko04360), and glycosylphosphatidylinositol (GPI)-anchor biosynthesis (path: ko00563) pathways may be important KEGG pathways of genes co-expressed with DE lncRNAs, and the related lncRNAs may be potentially involved in fimbriae adhesion to the intestinal mucosa. However, the role of these pathways in disease resistance remains largely unknown.
In this study, the expression profiles of lncRNAs in the spleens of lambs that were antagonistic or sensitive to diarrhoea were investigated to further understand the regulation of lncRNAs in sheep disease resistance. Several differentially expressed lncRNAs in lamb spleens between the antagonistic and sensitive groups were identified, and we found that 6 genes (MYO1G, TIMM29, CARM1, ADGRB1, SEPT4, and DESI2) may be regulated by their associated lncRNA. Our study may help elucidate the mechanism underlying resistance to diarrhoea in lambs. Further investigations of these sheep lncRNAs in relation to diarrhoea-resistance are warranted.     colony-forming units (CFUs)·mL −1 ], as well as ad libitum access to drinking water. Stool features of the experimental lambs were recorded daily ( Table 5). Lambs that exhibited diarrhoea for two days were classified as antagonistic and sensitive and were euthanized. The intestinal tissues were collected in 4% paraformaldehyde. The liver, spleen, duodenum, jejunum, and ileum of each lamb were collected and immediately frozen in liquid nitrogen until RNA extraction.

Methods
HE staining. The jejunum tissue was washed with 0.9% normal saline and fixed in 4% paraformaldehyde for 48 h at room temperature and then used in histological analysis. Next, 7 μm-thick sections were stained with haematoxylin-eosin and the morphology of the jejunum epithelia was assessed under a microscope.
Library construction and sequencing. RNA

Identification of lncRNAs and mRNAs.
The raw data were filtered to eliminate low-quality reads.
Clean reads mapped to the reference genome (Ovis aries v4.0) were selected for de novo assembly. Coding and non-coding RNA candidates from the unknown transcripts were categorized using four coding potential analysis methods, CPC 35   were employed to compare different expression levels in the same transcript in both samples, two criteria were selected: 1) fold change, which is the change in the expression level of the same transcript in both samples; and 2) the p value or false discovery rate (FDR) (adjusted p-value). The FDR error control method was used in correcting the p-value multiple hypothesis test 41 .

GO and KEGG pathway analyses.
After screening for differentially expressed transcripts, functional annotation was performed using GO enrichment analysis. Enrichment analysis employed counting the number of transcripts in each GO term, followed by Fisher's exact test to assess statistical significance (p < 0.05). KEGG 42 is the main public database used in pathway analysis, which was followed by Fisher's exact test to assess statistical significance (p < 0.05). The analysis of different transcripts was used to identify enriched pathways.
Prediction of the target genes of DE lncRNAs. The target genes of DE lncRNAs were predicted by calculating the Pearson correlation coefficients and P values among multiple genes. The |correlation| ≥ 0.7 and P ≤ 0.05 were used to filter the transcripts 43 . The DE lncRNAs associated with adhesion were selected, and the target genes of all DE lncRNAs were predicted by cis-acting 44 .

Verification of the expression level of DE lncRNAs.
To verify whether the screened DE lncRNAs play a role in the process of antagonism, q-PCR was used to detect the expression levels of 12 DE lncRNAs and DE mRNAs in the lamb spleens between the antagonistic and sensitive groups. The relative expression of each RNA was normalized to that of GAPDH using the 2 −∆∆Ct method 45 , and the primers used in the amplification of the lncRNAs are shown in Table 6.
Statistical analysis. All data were analysed by SPSS (version 20.0), and the relative expression levels of different transcripts were analysed by ANOVA. Tukey's test was used for multiple comparisons. Statistical significance was determined when p < 0.05. Each group contained three samples, and each experiment was repeated thrice.

Data Availability Statement
We guarantee that our data is valid.