Changes in circRNA expression profiles related to the antagonistic effects of Escherichia coli F17 in lamb spleens

Sheep colibacillosis is one of the most common bacterial diseases in large-scale sheep farms. In this study, we orally administered Escherichia coli F17 (E. coli F17) to lambs to obtain antagonistic and sensitive individuals. We used RNA-seq to screen for differential circRNAs in the spleens of both antagonist and sensitive individuals to explore the effect of circRNA on anti-diarrhoea in sheep. The results showed that 60 differentially expressed (DE) circRNAs were screened by RNA-seq in the spleen of antagonistic and sensitive lambs, among which 31 were up-regulated and 29 were down-regulated; q-PCR was used to validate the relative expression levels of six randomly selected circRNAs in antagonist and susceptible lambs and found to be consistent with the results of RNA-seq. Using Miranda analysis of circRNA-miRNA-mRNA interactions, we found a certain target relationship between 6 circRNAs, 5 miRNAs and 9 mRNAs. The relative expression levels of mRNA in antagonistic and sensitive lambs were verified by q-PCR and were consistent with the results of RNA-seq. This study explored the expression profile of circRNA in the spleen of an antagonistic and susceptible lamb with diarrhoea and found that differentially expressed circRNAs were helpful for determining how the lambs resist the pathogenesis of diarrhoea and provided a scientific basis for lambs to resist diarrhoea.

product of pre-mRNA processing; circRNA processing machinery competes with mRNA. At the same time, the classical clipping signal and clipping mechanism are also necessary for back-shear 11 . At present, studies on sheep disease resistance mainly focus on the prevention and treatment of disease 12,13 , but the important molecular mechanisms of disease resistance are seldom reported. In this study, we first screened for circRNAs that were differentially expressed in individuals that were antagonistic and sensitive to E. coli F17 fimbriae by RNA-seq. Miranda was used to analyse the circRNA-miRNA-mRNA interaction to find miRNA target genes and then verified by q-PCR. At the circRNA level, this study deepens our understanding of antagonizing E. coli F17 fimbriae in sheep and at the same time is expected to identify genes that can antagonize E. coli F17 fimbriae and solve key problems for Chinese local sheep breeding against E. coli disease. This study can be used as a foundation and can provide a theoretical basis for formulating a breeding strategy against E. coli in the future.

Results
Identification of transcripts in sheep spleens. After mapping the reference sequence, we identified 7,730 known circRNAs. The length of circRNAs primarily ranges from 200 bp to 900 bp, with an average length of 1943 bp and a GC content of approximately 43.5%. The number of exons of circRNA is mainly 2-4 ( Fig. 1). The statistics of the variable shear signal (GT-AG) of the reverse cleavage site in the circRNA sequence were calculated and graphs were drawn (Fig. 2a). The statistics for the circRNA types are shown in Fig. 2b: overlapping circRNA accounts for 92.11%, exonic circRNA accounts for 3.27%, intergenic circRNA accounts for 3.18%, intronic cir-cRNA accounts for 0.88% and antisense circRNA accounts for 0.56%. circRNAs were compared to the genomic elements to explore the distribution of circRNAs in the genome, to count the number of circRNAs predicted on each chromosome or scaffold, and to plot the results (Fig. 2c). It was found that circRNAs were primarily distributed on three chromosomes: NC_019458.2 (853), NC_019459.2 (772), NC_019460.2 (787).

Analysis and validation of DE transcripts.
We used the RPM value to estimate the expression level of cir-cRNA transcripts, and the expression level of circRNA transcripts was low (Fig. 3). We screened 31 up-regulated and 29 down-regulated DE circRNAs (Fig. 4). Differentially expressed circRNAs can be found on Table 1. To further validate the reliability of RNA-seq, 6 DE circRNAs were randomly selected, and their relative expression levels in antagonistic and sensitive lambs were confirmed by q-PCR and found to coincide with our RNA-seq results (Fig. 5), thus indicating that the RNA-seq data are reliable. Our analyses also show that high-throughput sequencing has the advantage of detecting genes with low expression levels. and GO databases showed that a total of 60 circRNAs were annotated and classified into 297 functional subclasses. The results showed the oxidation-reduction process (GO: 0055114), transport (GO: 0006810), extracellular region (GO: 0005576), focal adhesion (GO: 0005925), extracellular exosome 0005615), zinc ion binding (GO: 0008270) and seven more subclasses of circRNA functions, while the remaining functional subclass circRNA was less distributed (Fig. 6a). A comparison of the DE circRNA and KEGG PATHWAY databases showed that a total of 60 circRNAs were annotated and grouped into 73 KEGG PATHWAYS. The results showed that there were more circRNAs in three KEGG pathways, including the estrogen signalling pathway (path: ko04915), protein processing in the endoplasmic reticulum (path: ko04141) and regulation of the actin cytoskeleton (path: ko04810). The remaining KEGG pathways have fewer circRNAs (Fig. 6b).   Table 2 below. To further validate the relative expression, 9 mRNAs were selected, and their relative expression levels in antagonistic and sensitive lambs were confirmed by q-PCR. Significant differences were found between the two groups ( Fig. 7).

Discussion
Due to the limited efficiency of traditional molecular biology methods for the detection of circRNAs, circRNA has long been regarded as a product of the abnormal splicing of RNA 7 . In recent years, with the rapid development of bioinformatics and high-throughput sequencing technology, a large number of circRNAs have been identified in eukaryotes 6,14,15 , and they may play an important role in the regulation of gene expression 9,16,17 . It has been found that most circRNAs contain miRNA binding sites that act as efficient competitive endogenous RNAs and efficiently adsorb miRNAs to regulate the target genes of miRNA 18,19 , which have the following four biological functions: miRNA sponge effect 17,18 , protein translation template 20,21 , regulation of gene transcription 22,23 and regulation of competition for linear RNA production 8,10 . However, investigations of lamb diarrhoea relative to circRNAs are limited. Hu sheep are a unique breed with high fecundity and a strong adaptability to warm-wet climates; they can be kept indoors all year round. This study provides the first overview of circRNAs in relation to diarrhoea in sheep, as well an investigation into their possible roles in disease resistance.
In the present study, we found that the expression level of circRNA was very low: 5 miRNAs binding to 6 cir-cRNAs were identified by the Miranda software, and the adjacent genes of 4 circRNAs were identified as Btnl 1, GSTM1, NRAMP2 and B2M.
Btnl 1 is a key suppressor of T-cell activation and immune diseases 24 . The mechanism of action of Btnl 1 is different from those of Btnl 2 and BtnlA1, which directly inhibit T cell activation through anti-receptor binding [24][25][26][27] . The study found that the Btnl gene may be a new important local regulator of intestinal inflammation 28 . GSTM1 encodes the glutathione-S-transferase (GST) M1 enzyme, which is involved in the detoxification of various carcinogens of lung cancer 29 and plays a key role in protecting cells from oxidative stress 30 . NRAMP2 is a metal transporter protein, and in the absence of manganese, NRAMP2 is involved in the regeneration of Mn in the Golgi and promotes plant root growth 31 . B2M encodes the beta chain of the major histocompatibility complex (MHC) class I molecule and is up-regulated in inflammatory and tumour cells 32 .
We also used Miranda software to predict the 3 miRNA target genes and verified that they were significantly differentially expressed in antagonistic and sensitive groups, namely, NEB, UBE3B, ADGRF2, LAMA1, LTF, MGAT5, TLN2 and SLC25A29.
NEB encodes nebulin, a large protein component of the cytoskeletal matrix that coexists with myofilaments in skeletal muscle. Mutations in the NEB gene are the most common causes of myotubes, accounting for approximately 50% 33 . UBE3B is a ubiquitin ligase (UBE3), and its unique combination of E2-binding enzymes provides . Differentially expressed circRNAs in antagonistic and sensitive lambs Note: Gray represents circRNAs that have no significant differences; Red represents significantly upregulated circRNAs; Green represents significantly downregulated circRNAs; Blue indicates that the difference multiple is more than 2 times, but the circRNA is not significant in the difference significance test. The horizontal axis is the display of log2 FoldChange, and the vertical axis is the display of log10 Pvalue. Continued high substrate specificity, which is required to target specific protein degradation 34 . ADGRF2 is a member of the adhesion G protein-coupled receptor family and plays an important role in adhesion in the cell-cell and cell-matrix 35 . LAMA1 mutations may be related to Poretti-Boltshauser syndrome, and studies have shown that LAMA1 deficiency can lead to cytoskeletal changes 36 . LTF is a member of the transferrin family, whose protein products initiate host defence against a wide range of microbial infections and antigenic activity 37 . The protein encoded by MGAT5 belongs to the glycosyltransferase family and is one of the most important enzymes involved in the regulation of the biosynthesis of glycoprotein oligosaccharides. Changes in oligosaccharides on cell surface glycoproteins cause significant changes in cell adhesion or migration behaviour; increased enzyme activity is associated with the development of invasive malignancies 38 . Talin is a large adapter protein that links the integrin family of adhesion molecules to F-actin; Talin 1 is required for integrin-mediated cell adhesion, and TLN2, like Talin 1, is considered to be unique. Transmembrane receptors bind to form new connections between the extracellular matrix and the actin cytoskeleton 39 . SLC25A29 encodes a nuclear-encoded mitochondrial protein that is a member of the large family of solute transporter family 25 (SLC25) mitochondria. The primary physiological role of SLC25A29 is to introduce basic amino acids into the mitochondria for mitochondrial protein synthesis and amino-acid degradation 40 . GO is a bioinformatics tool that is widely used to study the functional relationship of genes. GO and KEGG pathway analyses of 60 DE circRNAs showed that the relevant circRNAs may potentially participate in the process of pili adhesion to intestinal mucosa. However, the role of these pathways in disease resistance remains largely unknown.
We found that a total of 60 circRNAs were significantly differentially expressed between antagonistic and sensitive groups, with 31 up-regulated and 29 down-regulated. In addition, we identified a total of 1,942 new circRNAs in both groups. To further verify the results of RNA-Seq, the expression levels of the six circRNAs were verified by q-PCR, and the results were consistent.
We studied the expression profiles of circRNAs in the spleens of antagonistic and sensitive sheep that developed diarrhoea to further understand their regulatory role in disease resistance in sheep. We found that circRNAs were differentially expressed in spleen tissues of antagonistic and sensitive lambs. Our research may help to determine how lambs resist the mechanism of diarrhoea. In addition, further studies of these circRNAs can provide a scientific basis for lambs to resist diarrhoea.  Table 1. Differentially expressed circRNA.

Figure 5.
Relative expression levels of DE circRNAs between antagonistic and sensitive lambs Note: "**" means highly significant correlation; "*" means significant correlation; "ns" or "no SuperiorScript" means no significant correlation. The same as below.
SCieNtifiC Identification of circRNAs. The circBase 43 database contains only circRNA sequences of human, mouse, nematode, lagomorph, and coelacanth. Since sheep are not included, we used CIRI 44 de novo prediction of cir-cRNA. According to the position of circRNA on the genome, circRNAs can be classified into the following five categories: exonic circRNA, intronic circRNA, antisense circRNA, sense overlapping circRNA, and intergenic circRNA.

Different expression analysis.
After obtaining differentially expressed circRNAs, we performed Gene Ontology (GO) and KEGG pathway significance analyses of the source genes. DESeq 45 is suitable for experiments with biological duplication. Differential expression analysis can be performed between sample groups to obtain the circRNA list for the difference between the two biological conditions. For experiments without biological duplication, edgeR 46 differential expression analysis was used to obtain a list of circRNAs that were differentially expressed between the two samples.

GO and KEGG pathway analyses.
After screening for differentially expressed transcripts, functional annotation was performed using GO enrichment analysis. Enrichment analysis involved counting the number of transcripts in each GO term, followed by Fisher's exact test to assess statistical significance (p < 0.05). KEGG 47 is the main public database used in pathway analysis, which was followed by Fisher's exact test to assess statistical significance (p < 0.05).
circRNA-miRNA-mRNA interaction studies. As a miRNA target molecule, circRNAs are regulated by miRNAs. Because circRNAs contain multiple miRNA binding sites, analysis of circRNA-miRNA interactions can help elucidate the function and mechanism of circRNA acting as a sponge. For animals, we used Miranda 48,49 to predict circRNAs that bind to miRNAs and the target genes of miRNA and elucidated the function of this portion of circRNAs based on the functional annotations of miRNA target genes.  Table 2. Prediction of the target relationship of circRNA-miRNA-mRNA. Verification of the expression level of DE circRNAs. To verify whether the screened DE circRNAs play a role in the process of antagonism, q-PCR was used to detect the expression levels of DE lncRNAs and target genes of DE miRNAs in 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 50 , and the primers used in the amplification of the lncRNAs are shown in Table 3.
Statistical analysis. All data were analysed by SPSS (version 22.0), and the relative expression levels of various transcripts were analysed by one-way ANOVA. Statistical significance was determined when p < 0.05. Each group contained three samples, and each experiment was repeated three times.

Data Availability Statement
We guarantee that our data is valid.