Comparative Transcriptome Analysis Reveals Relationship among mRNAs, lncRNAs, and circRNAs of Slow Transit Constipation

Background Slow transit constipation (STC) is characterized by persistent, infrequent, or incomplete defecation. Systematic analyses of mRNA, lncRNA, and circRNA expression profiling in STC provide insights to understand the molecular mechanisms of STC pathogenesis. The present study is aimed at observing the interaction of mRNAs, lncRNAs, and circRNAs by RNA sequencing in vivo of STC. Methods A rat model of STC was induced by loperamide. The expression profiles of both mRNAs and miRNAs were performed by RNA sequencing. Enrichment analyses of anomalous expressed mRNAs, lncRNAs, and circRNAs were performed in order to identify the related biological functions and pathologic pathways through the Gene Ontology (GO) database and Kyoto Encyclopedia of Genes and Genomes (KEGG) database. Results In total, 26435 mRNAs, 5703 lncRNAs, and 7708 circRNAs differentially expressed were identified between the two groups. The analyses of GO and KEGG show that (1) upregulated genes were enriched in a positive regulation of GTPase activity, cell migration, and protein binding and lipid binding and (2) GO annotations revealed that most trans-target mRNAs are involved in the regulation process of immune signal together with the proliferation and differentiation of immune cells. Additionally, the protein-protein interaction (PPI) network of differentially expressed (DE) mRNAs was constructed. Interestingly, all of the core lncRNAs and their coexpression mRNAs in this network are downregulated. Moreover, downregulated circRNAs have a set of target mRNAs related to immunoreaction, which was consistent with the overall tendency. Conclusion Our investigation enriches the STC transcriptome database and provides a preliminary exploration of novel candidate genes and avenues expression profiles in vivo. The dysregulation of mRNAs, lncRNAs, and circRNAs might contribute to the pathological processes during STC.


Introduction
Slow transit constipation (STC) is one of the refractory digestive tract diseases. Commonly, this syndrome concludes slow colonic peristalsis and delayed excretion of intestinal contents, except normal rectal discharge and normal pelvic floor function, which is the most common subtype of functional constipation [1]. Reference investigation shows that the incidence of chronic constipation is increasing, even though it has become one of the key factors that affect people's quality of life globally. The epidemiological survey of chronic consti-pation shows that the incidence of chronic constipation is 9.9% in southern China, besides the prevalence rate significantly increases with age [2]. According to the statistical data of Europe and the United States, the incidence of constipation may be higher than expected. At least 65% of constipation patients are treated with laxatives on their own [3]. Patients with chronic constipation often have progressive difficulty in defecation, abdominal distension, sometimes they will gain introverted personality even depression. The above syndrome seriously declines life quality and could lead to complications (i.e., myocardial infarction, stroke, and colorectal cancer) [4]. Meanwhile, patients with intractable constipation show obvious psychological problems, which are often accompanied by anxiety, paranoia, obsessivecompulsive disorder and social maladjustment. It affects life quality and physical and mental health; with these, it causes serious burden to the social economy [5].
Current increasing studies have shown that STC has numerous pathogenesis and complex mechanism. Therefore, it is difficult to fully clarify the pathogenesis from one facet. Additionally, the clinical treatment is not sufficient [6,7]. Long noncoding RNAs (lncRNAs), a little understood type of transcribed RNA molecules, have over 200 nucleotides in length and no significant protein-coding capacity, which have been identified as the key regulators of various biological functions [8,9]. circRNAs are a genetic element, which are evolutionarily conserved and covalently closed. Some of them are rich in eukaryotes possessing cell-specific and tissue-specific expression profiles [10].
To disclose complex and heterogenous mechanisms of STC, we used a rat model of STC induced by loperamide and simultaneously performed mRNA, lncRNA, and cir-cRNA microarray analyses to identify mRNA, lncRNA, and circRNA interactions, plus explored how these interactions influence on the pathogenesis of STC.

Animals.
Male Sprague-Dawley rats aged 6-9 weeks were purchased from the Jiangsu Laboratory Animal Centre in Suzhou, China (License approval No.: SCXK [Su] 2017-0007). They were fed at a room temperature of 26°C with 45%-55% of the humidity and 12 h light/dark cycles. The environment was clean and quiet (low noise level ≤ 60 dB) in a well-ventilated place. All experiments involved in this study were performed under the requirements of the Provision and General Recommendation of the Chinese Laboratory Association. The study was approved by the Medical Research Committee on Animal Care and Use of Suzhou TCM Hospital Affiliated to Nanjing University of Chinese Medicine.

Loperamide-Induced Constipation
Model. In this study, animals were randomized into two groups: the control group (n = 5) and Lop group (n = 5). Rats were given normal saline in the control group, while the Lop group were induced with 4 mg/kg loperamide (an antidiarrheal drug) suspension to make the STC models (oral administration, twice per day: 09:00 and 17:00) for 14 days [11].

Parameters Evaluated.
Fecal pellets of rats were collected after treatment of 24 h to collect their total number, stool weight, and water content. All measurements were performed five times.

Measuring
Intestinal Charcoal Transit Ratio. The intestinal motility by charcoal meal was assessed following Kim et al. [12]. Briefly, on day 14, all rats were fasted for 12 h but no limited water, after which, they were fed charcoal within 10% acacia gum. 0.5 hours later, the rats' intestines from pylorus to ileocecal junction were removed; then, the total length of the truncated intestine and the charcoal transport distance was measured. Finally, the rate of intestinal motility was calculated using formula (1 The blood of 5-10 mL was collected from the abdominal aorta and injected into the tube of EDTA anticoagulant. Then, it was separated 2 mL of plasma by centrifugation (3000 r/min, 30 min) to detect indexes of gastrin (GAS), motilin (MTL), substance P (SP), and 5-hydroxytryptamine (5-HT) by commercial ELISA kits (Biolegend, San Diego, USA). The procedure was strictly performed according to the kit instructions.

Histopathological and Immunohistochemical Analyses.
Colon tissue samples were collected from the sacrificed Sprague-Dawley rats. They were fixed at room temperature with 10% buffered formalin for 48 hours. The fixed colonic tissues were embedded in paraffin before being sliced into 5 μm thick sections. Then, the slices were deparaffinized and stained using hematoxylin and eosin (HE; Sigma-Aldrich Co.). After that, we analyze their histological morphology and the thickness of mucosa and muscle with the Leica Application Suite (Leica Microsystems, Switzerland). C-kit proto-oncogene protein (C-kit) and aquaporin 3 (AQP3) expression levels were analyzed with immunohistochemistry (IHC). Then, the paraffin colon sections of rats were routinely fix with 4% paraformaldehyde for 10 minutes; after rinsing with phosphate-buffered saline for 3 times (3 minutes per time), 50 μL of peroxidase blocking solution dropwise was added and incubated at room temperature for 10 minutes. Then, the slides were incubated consecutively overnight at 4°C with anti-CD117 (1 : 1000) and AQP3 (1 : 1000) primary antibody and at room temperature with HRP-conjugated anti-rabbit IgG incubation for 30 minutes. They were then stained using diaminobenzidine (DAB). By applying the high-power microscope (×200) and Image-Pro Plus 6 graphic processing software, the expression of C-kit and AQP3 per unit colon tissue area was observed and analyzed. Each section was observed with 5 fields randomly. The optical density values of C-kit and AQP3 expression of rats in both groups were statistically and, respectively, assessed. When identifying circRNA, CIRCexplorer (v2.2.3) was used to find circularizing junction and spliced sequence with the fusion junctions obtained from TopHat2. Candidates with junction reads ≥ 2 were considered bonafide circRNAs, and the expression levels of circRNAs were estimated by TPM (transcript per million).

Identification of Differentially Expressed Transcripts.
The mRNAs (DE mRNAs) and lncRNAs (DE lncRNAs) expressed differentially between the STC and non-STC groups were identified by DESeq2, and the different expressions of cir-cRNAs (DE circRNAs) were identified using the limma package in R. log 2 fold change | ≥1 and P value < 0.05 were indicated.

Analysis of Target
Genes Regulated by lncRNAs. The lncRNA function is mostly achieved by acting on target genes in cis or trans. cis-acting elements were DNA sequences that are adjacent to the structural portion of a gene that regulates gene expression. Hence, we selected DE mRNAs within 100 k upstream or downstream of DE lncRNAs. trans-acting factors usually are proteins binding to cis-acting elements to control gene expressions. Therefore, a coexpression network of lncRNA-mRNA was constructed according to the interregulatory correlation between DE lncRNAs and DE mRNAs. The Pearson correlation coefficient (PCC) was calculated using log 2 ðFPKMs + 1Þ, when |PCC | ≥0:98 and P value < 0.001 were considered meaningful. The network was visualized by Cytoscape.
2.11. Functional Enrichment Analysis. The enrichments of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were performed for the host genes of DE mRNAs by using DAVID (v6.8) and KOBAS (v3.0) software, respectively. P < 0:05 was recognized significant.

Effect of Loperamide on Fecal
Pellets. The STC rat model was established by loperamide inducing. The total number and weight of feces were remarkable lower in the Lop group than the control group (Figures 1(a) and 1(b)). In addition, fecal water content was decreased more than the Lop group ( Figure 1(c)). The outcomes of fecal parameter changes indicate that the in vivo model was built (Figure 1).

Effect of Loperamide on Intestinal Transit
Rate. The intestinal propelling movement of carbon ink displayed that the intestinal transit rate in the Lop group was remarkably lower than that in the control group ( Figure 1(d)).

Effects of Loperamide on Neurotransmitter Concentration.
Commonly, constipation would lead to colonic motor dysfunctions, intestinal neurological abnormalities, and disease states. The combination of myogenic, intestinal plexus, and extrinsic neurons could affect colonic motor activity [12]. Gas, MTL, SP, and 5-HT play an important role in regulating gastrointestinal motility. Loperamide can significantly reduce more Gas, MTL, SP, and 5-HT than the control group (Figures 2(a)-2(d)).

Histopathology and Immunohistochemistry
Findings. The pathological changes in the colon tissue were observed in the Lop group, with damage to colonic mucosa, thickening of the submucosal muscle, and the decrease in intestinal water ( Figure 3). A group of cells, which exists in all layers of the colon [13] and interstitial cells of Cajal (ICC), plays an important role in regulating intestinal motility. C-kit could maintain the normal phenotype of ICCs during development and maturation [14]. AQPs have recently been observed to have an integral effect on human water transport systems [15]. In the human colon, AQP3 has been demonstrated to be important to colonic water transport, and it is mainly expressed in mucosal epithelial cells [16]. In the con-trol group, C-kit-and AQP3-positive immunoreactive structures are dyed in brown. The expressions of C-kit and AQP3 decreased more in the Lop group than the control group.

Quality Test and Alignment Analysis of RNA-seq Data.
The clean data of 9 samples were obtained with an average of 56.63 million reads and 16.99G bases per sample ( Table 1). The Q30 ratio was >98%, and no GC bias was observed, suggesting sequencing clean data are qualified. Over 94% of the clean reads were perfectly mapped to the rat reference genome (rn6), and 67.83-77.50% uniquely mapped reads were obtained from the total reads from the 9 samples, which indicated an excellent performance of the sequencing reads with high credibility of the results in downstream analysis.
3.6. Characterization of Transcripts. After assembly, a total of 26435 mRNAs, 5703 lncRNAs, and 7708 circRNAs were identified in STC rats among 9 samples. Overall, the chromosomal distribution of mRNAs, lncRNAs, and circRNAs were consistent. Most transcripts were mainly located on the first 10 chromosomes (Figure 4(a)), in which the larger chromosomes contained more lncRNAs. The transcript lengths of mRNAs were longer than those of noncoding RNAs, especially the length greater than 1000 nt. 80% of the total number of noncoding RNAs was no longer than 2000 nt, while 60% of them were approximately 200 to 2000 nt in length ( Figure 4(b)). The amount of novel lncRNAs accounted for nearly half proportion of mRNAs (Figures 4(c) and 4(d)), which suggested that there were still many unknown regions in the rat genome. We also classified noncoding RNAs that were identified in this study and the numbers and proportions of different types for lncRNAs ( Figure 4(e)) and cir-cRNAs ( Figure 4(f)).

Differential Transcription Expression Profile Analyses.
Correlation coefficient and principal component analysis (PCA) of all transcripts could represent the degree of similarity between groups. The analysis results showed that within group were highly correlated, while between groups were clearly distinguished ( Figure 5).
To identify the different expression of mRNAs (DE mRNAs), lncRNAs (DE lncRNAs), and circRNAs (DE cir-cRNAs), P ≤ 0:05 and |log 2 fold change | ≥1 were used as the threshold. Volcano plot showed a set of significantly differentially expressed transcripts were aptly delimited both in the STC and non-STC groups (Figures 6(a), 6(c), and 6(e)). A total of 917 DE mRNAs were determined, including 298 upregulated and 619 downregulated mRNAs. DE lncRNA analysis indicated 419 deregulated lncRNAs, among those, significantly upregulated and downregulated lncRNAs were 96 and 323, respectively. For circRNAs, out of 112 deregulated, 52 upregulated and 60 downregulated were observed. Unsupervised cluster analysis of the DE mRNAs, DE lncRNAs, and DE circRNAs also showed an obvious expression patterns between two groups (Figures 6(b), 6(d), and 6(f)). All differentially expressed transcripts were listed in Table S1.    5 BioMed Research International metabolism were significantly changed, which indicated that the activation or inhibition of related gene expressions might accelerate the pathological process.

Protein-Protein Interaction (PPI) Network of DE mRNAs.
To investigate the important role of protein interactions in STC rats, a PPI network analysis were performed using the STRING based on the top 300 DE mRNAs ( Figure 8). We found that 64 proteins formed a complex functional network, and proteins which had high connectivity with other proteins, including Dync1h1 (degree = 13), Cd19 (degree = 9), Ptprc (degree = 9), Hsp90aa1 (degree = 9), Dync1i2 (degree = 8),  4 5 6 7 8 9 10 11 12 13 14 15 16 17 18     BioMed Research International Ccr7 (degree = 8), and Fn1 (degree = 8) were identified hub proteins. Hub proteins were widely involved in important cell cycle or immune-related pathways (Table 2), which indicated that dysregulation in one of these proteins can have significant effects on STC.

Functional Identification of DE lncRNAs.
Previous studies reported that lncRNAs cis-regulated neighboring and overlapping genes, or coexpression with their target genes by trans-factors [17,18]. To elucidate the potential function of lncRNAs, DE mRNAs within a 100 kb of upstream or downstream in the DE lncRNAs were considered cis-target, while targets in trans were predicted via calculating the expressed correlation. A total of 109 possible cis-regulatory relationships were identified (Table S2). Among them, the proportion of downregulated lncRNAs exceeded 80%, and the number of downregulated mRNAs also accounted for more than half. Target genes were widely involved in immune-related biological processes (Figure 9(a)), which indicated that suppression of immune function in STC rats was mediated by lncRNAs.
Based on the absolute Pearson correlation coefficient over 0.98 and P values less than 0.001, we identified 904 pairs of coexpressed lncRNAs and mRNAs; all of them showed strong positive correlations (Figure 9 Cmahp, degree = 37). Similar to cis action, the relationships in the trans-regulatory network were mostly in the state of cosuppression, and GO annotations showed that the vast majority of the trans-target mRNAs are involved in the regulation process of immune signal and the proliferation and differentiation of immune cells (Figure 9(a)). Interestingly, all the core lncRNAs and their coexpression mRNAs in this network are downregulated (Figure 9(b)). It suggested that these lncRNAs might be the key factors to reduce the immunity.
3.11. CeRNA Network Analysis. circRNAs generally act as the sponge of miRNAs. By competing with mRNAs in a sequence complementary manner to binding miRNAs, cir-cRNAs relieve miRNAs from inhibiting mRNAs translation, thereby exerting regulatory functions on protein coding genes [17,18]. Therefore, we predicted the circRNAs-miRNAs-circRNAs binding relationship and constructed networks between DE circRNAs and DE mRNAs that were identified in this study.
The co-activated network contained 22 circRNAs and 178 mRNAs, and the cosuppressed network contained 30 cir-cRNAs and 425 mRNAs (Figures 10(a) and 10(b)). More than half of DE mRNAs had a potential regulatory relationship with circRNAs, indicating that circRNAs were generally involved in the regulation of related pathological processes. GO enrichment analysis of the target mRNAs in the coactivated network showed that many biological processes related to the nervous system and intercellular signal transduction were activated, such as regulation of establishment or maintenance of cell polarity, negative regulation of neuron death,  (Figure 10(c)). These functions were particularly more dominant than enriched functional of all upregulated mRNAs (Figure 10(a)), which implied that they were more regulated by circRNAs. A significant set of target mRNAs were related to immunoreaction (Figure 10(d)) in downregulated circRNAs, which was consistent with the overall trend of the functions for downregulating mRNAs or lncRNAs. It was shown that immune and inflammatory responses were universally regulated in STC rats.

Discussion
STC is an intestinal disease that affects the health and quality of life of patients. To reduce the incidence of STC, it is essen-tial to understand the underlying pathophysiological mechanism. In this study, a group of differentially expressed transcripts and some key factors were identified that may regulate STC at the molecular level using a whole-transcriptome sequencing analysis in the colon tissue of STC rats.
There are larger amounts of deregulated mRNAs in STC rats than normal rats (Figure 9). According to enrichment analysis, the downregulated prominent functions are immune cell receptor and inflammatory pathways with its related signaling pathways (Figure 10(b)). Previous studies have shown that inflammation is closely related to gastrointestinal motility disorders; moreover, the release of inflammatory mediators leads to changes in gastrointestinal motility [19]. Our results are consistent with this, suggesting that inflammation may play an important role in STC   (Figure 10(a)). Disorders of energy metabolism strongly affect the homeostasis of intestinal environment, which can accelerate a series of diseases progression such as intestinal inflammation and neoplastic pathology [20]. The release of inflammatory mediators can induce acute inflammatory cell infiltration and promote NF-κB activity [21,22], related pathways, which can be observed from upregulated mRNAs. The mRNAs included positive regulation of GTPase activity, cell migration and I-kappaB kinase/NF-kappaB signaling, phagosome, and glycerolipid metabolism (Figure 10 A).
Constipation symptoms are also inseparable with enteric nervous system dysfunction and changes of neurotransmitters that regulate intestinal motility [23]. In our study, the nervous system development, immunological synapse, and various cell receptor signaling pathway are downregulated ( Figure 10). Interstitial cells of Cajal (ICC) are involved in intestinal neuromuscular signal transmission, which play an important role in regulating intestinal smooth muscle    Table S1 and Table S2. contraction and controlling gastrointestinal tract movements. Loss and dysfunction of ICC related to inflammation are the causes of STC [24][25][26].
In addition, we observed a large number of transcripts, which are related to encoding tumor necrosis factor ligands (Tnfsf8, Tnfsf9, and Tnfsf11) and receptors (Tnfrsf13c, Tnfrsf26, and Tnfrsf9), are downregulated in STC rats (Table S1). Tumor necrosis factor alpha (TNFα) and its related transcripts, usually acting as pathogenic proinflammatory cytokines, have been shown imbalanced in inflammatory intestinal tissues and peripheral lymphocytes [27,28]. TNF-α is also known to be a particularly toxic agent to mitochondria, render whose dysfunction, and lead to cytokine storm of inflammatory activity [29]. Therefore, it is entirely plausible that the downregulation of TNF-α-related transcripts through affecting the mitochondrial function to shift energy production and activity of the inflammatory tissue in STC.
Among PPI network, the hub proteins are encoding products from vital mRNAs, and they were speculated to be important factors regulating STC (Figure 8). Dync1h1 encodes a member of the dynein cytoplasmic heavy chain family; its dysfunction can cause gut dysmotility [30]. Cd19 is a B cell coreceptor that is important for B cell development; its low expression can negatively affect the intestinal physiology under steady-state conditions [31]. Hsp90aa1 functions as a molecular chaperone and assists in the assembly, folding, and degradation of target proteins. The loss of Hsp90aa1 function is often accompanied by inflammation and other diseases [32]. Fn1 encodes fibronectin, and overexpression can cause fibrosis in various organs or tissues [33]. It is speculated that Fn1 may reduce the efficiency of intestinal transit and increase the risk of STC.
There are thousands of noncoding RNAs in human genome, which are widely involved in biological processes, such as imprinting and chromosomal conformation [34], coordination of cell status and differentiation [35,36], enzymatic regulation [37], and disease [38]. Thus, the dysregulation of noncoding RNAs is closely related to diseases. However, there are still many gaps between the current understanding of the expression pattern and molecular process of noncoding RNAs in a specific gut environment. In our research, we use rats as animal models and focused on the regulation of two important noncoding RNAs, lncRNAs, and circRNAs.
Firstly, we annotate all expressed transcripts at the whole transcriptome level and describe their characterization ( Figure 4). Nearly 90% of the coding genes are already recorded in related database, while half of the lncRNAs are recently identified. It indicated that there are still many unknown noncoding functional regions that need to be elucidated in the rat genome. Transcriptome-wide studies have shown that the expression of noncoding RNAs is usually more specific and strictly regulated than protein coding genes [39]. Thus, we speculate that these novel transcripts may be somehow related to disease process and partially specifically expressed in STC. Overall, our results further provided reference insights for the transcript expression of rats under STC status and enriched the rat genome resources.
The complex interaction of noncoding RNAs and coding genes are widely involved in the regulation of intestinal homeostasis and physiological processes [40]. By predicting the regulatory network of lncRNAs or circRNAs for coding genes, we observed many noncoding RNAs are involved in the pathological process of STC (Figures 9 and 10). The most prominent functions of target genes are immune and inflammatory responses, followed by the nervous system and intercellular signal transduction. This is consistent with the overall function of mRNAs, suggesting that a considerable portion of expression of the coding genes is regulated by noncoding RNAs, both at the transcriptional and posttranscriptional levels. Among networks, a group of key noncoding RNAs were identified. They are the core regulator that has more potential coding gene targets.
The current study illuminated the expression profiles of mRNAs, lncRNAs, and circRNAs in the loperamideinduced constipation in rats. We identified mRNAs, lncRNAs, and circRNAs with differential expression between the Lop group and the control group and elucidated the characteristics of mRNAs, DE lncRNAs, and regulatory functions of DE mRNAs. Besides, we tapped several core regulators that may contribute to the maintenance of intestinal transit. Our findings may provide useful insights into the molecular mechanisms underlying the development of STC. Further research is required to investigate the functions of mRNA, lncRNA, and circRNA identified in the present study.

Data Availability
The authors declare that the data supporting the findings of this study are available within the article and its Supplementary Information File or from the corresponding authors upon request.

Conflicts of Interest
The authors declare that they have no conflicts of interest.