Aberrantly Expressed Genes and miRNAs in Slow Transit Constipation Based on RNA-Seq Analysis

Background This study aims to identify the key genes and miRNAs in slow transit constipation (STC). Methods MRNA and miRNA expression profiling were obtained. Differentially expressed genes (DEGs) and miRNAs were identified followed by the regulatory network construction. Functional annotation analysis and protein-protein interaction (PPI) network were conducted. The electronic validation was performed. Results Hsa-miR-2116-3p, hsa-miR-3622a-5p, hsa-miR-424-5p, and hsa-miR-1273-3p covered most DEGs. HLA-DRB1, HLA-DRB5, C3, and ICAM were significantly involved in staphylococcus aureus infection. The PPI network generated several hub proteins including ZBTB16, FBN1, CCNF, and CDK1. Electronic validation of HLA-DRB1, PTGDR, MKI67, BIRC5, CCNF, and CDK1 was consistent with the RNA-sequencing analysis. Conclusion Our study might be helpful in understanding the pathology of STC at the molecular level.


Introduction
Constipation is a gastrointestinal disorder characterized by infrequent bowel movements, defecation difficulty, and incomplete bowel evacuation sensation [1][2][3]. Constipation is more widespread in women and is often increased with age [4,5]. It is reported that constipation is associated with inflammatory bowel disease, irritable bowel syndrome, and relevant comorbidities [6]. Additionally, medications, mobility loss, gut aging alterations, and rectal sensory-motor dysfunction can cause constipation. It is worth mentioning that slow transit constipation (STC) is the prototype of functional constipation and frequently appears in general medicine and clinical practice (such as gastroenterology and geriatrics). Bharucha AE et al. found that STC was related to endocrine and metabolic disorders including hypercalcemia, hypothyroidism, or diabetes mellitus [7].
Recently, the underlying pathophysiology of STC is poor to predict and physiological assessment of gastrointestinal tract transit is the common diagnostic method. Taking laxative remains the primary treatment option in severe STC [3,8] and patients may take laxative for many years in many instances. Despite advances in therapy of STC, treatment failure is a common phenomenon [9]. Therefore, exploring the underlying pathology mechanism is urgent. In physiological conditions, colon motor activity increases after meals and wake up and decreases during sleep [10]. And STC is characterized by prolonged transit time of the stool through the colon on account of the reduced colonic transit rate [11]. Maybe, altered colonic function plays an important role in patients with STC.
To our knowledge, expression profiling of mRNAs and miRNAs in STC colon tissue through RNA-sequencing has not been reported. In this study, we used high-throughput RNA-sequencing to study the mRNA and miRNA expression profiling and investigated the functional significance of aberrantly expressed genes and miRNAs in STC.

Patients and Samples.
In this study, we selected patients with slow transit constipation (STC) strictly according to international or domestic standard. Those

Quality Control of Raw Sequence Data.
The raw image data of mRNAs and miRNAs derived from high-throughput RNA-sequencing was translated into raw FASTQ sequence data by using Base-Calling. Quality control of FASTQ data (Read QC) was performed by using FastQC v0.11.4. To obtain the high quality clean data, the low quality sequences (adaptor sequences, sequences with a quality score <20, and sequences with N base rate of raw reads >10%) were removed from sequencing data of mRNAs and miRNAs by using cutadapt v1.9.1.

Identification of DEGs in STC.
According to human UCSC genome reference annotation, the alignment between cleaned mRNA sequencing reads was aligned to the human genome (GRCh38.p7 assembly) via Tophat (http://tophat .cbcb.umd.edu/). The fragment was performed to assemble and relative expression of the reads with the normalized RNA-sequencing fragment was counted by using Cufflinks (http://cufflinks.cbcb.umd.edu/). Fragments per kilobase of exon per million mapped reads (FPKM) were calculated by Cuffdiff (http://cufflinks.cbcb.umd.edu/) to determine the transcription abundance of mRNAs. Paired t-tests were performed to test the expression between STC and controls. DEGs in STC compared to controls were identified with value < 0.05 and abs (log2 (fold change)) > 1. The heat map of top 100 DEGs in STC was obtained by heatmap.2 (http://www.bioconductor.org/packages/release/bioc/html/ heatmaps.html).

Network Construction of Differentially Expressed miRNA
Targets. Identifying target genes is an important step in studying the function of miRNA in tissues. In this study, the DEGs whose expression was inverse with corresponding differentially expressed miRNAs were obtained and checked by six bioinformatic algorithms (RNA22, miRanda, miRDB, miRWalk, PICTAR2, and Targetscan). Those targets recorded by more than 4 algorithms were selected. And verified targets were also obtained by miRWalk (http://www .umm.uni-heidelberg.de/apps/zmf/mirwalk/). According to the differentially expressed miRNA-target pairs, the STCspecific differentially expressed miRNA-targets interaction network was established by Cytoscape software (http://www .cytoscape.org/).

Functional Annotation of Differentially Expressed miRNA
Targets. In order to study the biological function of differentially expressed miRNA-targets, the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were performed by using the online software GeneCodis (http://genecodis.cnb.csic.es/analysis). And the threshold of FDR<0.05 was set as the criteria of statistical significance.

Protein-Protein Interaction Analysis.
In order to understand the protein interaction between DEGs in STC, BioGRID database was utilized to select interacting protein pairs. And the protein-protein network (PPI) was visualized by Cytoscape software (http://cytoscape.org/). In the network, nodes represent proteins and edges represent interaction between two proteins.

Validation the Expression of DEGs and Differentially
Expressed miRNAs by GEO. In order to validate the identified DEGs and differentially expressed miRNAs in STC, the Gene Expression Omnibus (GEO, http://www.ncbi.nlm .nih.gov/geo) database was used for electronic validation. We compared expression levels of selected DEGs and hsa-miRNAs between STC cases and normal controls and the difference of expression levels was showed by boxplots.

DEGs and Differentially Expressed miRNAs in STC.
A total of 464 DEGs (357 upregulated and 107 upregulated genes) and 10 differentially expressed miRNAs (4 upregulated and 6 downregulated miRNAs) were obtained. Heap maps of top 100 DEGs and all differentially expressed miRNAs were displayed in Figures 1 and 2, respectively. The top 20 DEGs were shown in Table 2.

Differentially Expressed miRNAs-Targets Network.
A total of 87 differentially expressed miRNAs-target pairs including 25 upregulated miRNA-downregulated target pairs and 47 downregulated miRNA-upregulated target pairs were obtained. Among which, there were 17 differentially expressed miRNAs-target pairs that have been confirmed by miRWalk. Based on the STC-specific differentially expressed miRNAs-target interaction network (Figure 3), hsa-miR-2116-3p (degree=27), hsa-miR-3622a-5p (degree=16), hsa-miR-424-5p (degree=15), and hsa-miR-1273-3p (degree=10) were the top four differentially expressed miRNAs covered most DEGs.  Figure 1: Hierarchical clustering analysis based on the expression profile of the top 100 DEGs in STC compared to normal tissues. Orange and light blue color mean the group of STC and normal tissues, respectively. Red color represents the relative expression level of genes which was higher than mean, and green color represents the relative expression of genes which was lower than mean. Orange and light blue color mean the group of STC and normal tissues, respectively. Red color represents the relative expression level of genes which was higher than mean, and green color represents the relative expression of genes which was lower than mean. Figure 3: STC-specific differentially expressed miRNA-targets interaction network. Ellipses were used to represent the DEGs and red and green color were used to represent up-and downregulation in STC, respectively. Rectangles were used to represent differentially expressed miRNAs and rose and blue color were used to represent up-and downregulation, respectively.

Validation the Expression of DEGs and Differentially
Expressed miRNAs by GEO. Based on the RNA-sequencing data, six downregulated DEGs (HLA-DRB1, PTGDR, MKI67, BIRC5, CCNF, and CDK1) in STC were selected to perform the expression validation by GEO dataset (Figure 10). Different expression levels of these genes and miRNAs were analyzed and depicted through box-plots. HLA-DRB1, PTGDR, MKI67, BIRC5, CCNF, and CDK1 were significantly downregulated in the case group compared to the normal control group, which was consistent with our RNA-sequencing data.

Discussion
STC is an intestine disease that influences the health and life quality of patients. Therefore, understanding the pathology mechanism will make a contribution to reducing the incidence rate of STC. In our study, we identified several DEGs and differentially expressed miRNAs in STC, which may play roles in the molecular level of STC.
In top 20 DEGs, we found two cell proliferation-related genes including MKI67 and BIRC5 in colon tissue of STC. Marker of proliferation Ki-67 (MKI67) is a cellular proliferation marker that is detectable in the cell nucleus during active phases of the cell cycle except resting cells [12]. It has been demonstrated that MKI67 plays an important role in colon cancer [13]. Herein, we also found that MKI67 was one of target genes of hsa-miR-1291 in STC. It has been indicated that hsa-miR-1291 can promote apoptosis in esophageal squamous cell carcinoma by regulating mucin1 [14]. Additionally, it is reported that hsa-miR-1291 is upregulated in colon cancer [15]. Baculoviral IAP repeat containing 5 (BIRC5, also called survivin) is a proliferation marker and belongs to the members of the IAP family of proteins [16,17]. BIRC5 controls a number of cellular processes including cell apoptosis [18]. It is found that BIRC5 is highly expressed in colon cancer cell lines [19]. And increased expression of BIRC5 is related to poor clinical outcome in patients with colon cancer patients [20]. In addition, it is differentially expressed in irritable bowel syndrome rectosigmoid mucosa [21]. In our study, we found downregulated expression of cell proliferation markers of BIRC5 and MKI67 in the colon tissue of STC, which indicated that BIRC5 and MKI67 may be involved in the cell proliferation and apoptosis in the development of STC.
Microbiota infection is common in different intestine disease. Postinfectious irritable bowel syndrome (PI-IBS) presents a bout of gastroenteritis caused by viral, parasitic, or bacterial infections [22]. It is worth mentioning that staphylococcus aureus infection is a potential pathogenic factor in PI-IBS [23][24][25][26]. It is reported that inflammatory bowel disease causes inflammatory obstruction of the gastrointestinal tract, resulting in symptoms such as constipation [27]. This suggested the correlation between inflammation and constipation. In this study, we found that staphylococcus aureus infection was the most significantly enriched pathway

KEGG Pathways
−log(FDR) of DEGs in STC. And four DEGs including HLA-DRB1, HLA-DRB5, C3, and ICAM1 were involved in this signal pathway. Major histocompatibility complex, class II, DR beta 1 (HLA-DRB1) has been demonstrated to be associated with ulcerative colitis and Crohn's disease [28,29]. Furthermore, HLA-DRB1 is differentially expressed in precancerous colorectal adenomas mucosa [30]. Another major histocompatibility complex family member, major histocompatibility complex, class II, DR beta 5 (HLA-DRB5) is found significantly upregulated in human colorectal cancer epithelial cell and has prognostic value in the prediction of colon cancer metastasis [31,32]. Herein, we found downregulated expression of both HLA-DRB1 and HLA-DRB5 in STC. Interestingly, they were under the regulation of hsa-miR-939-5p. It is reported that low expression of hsa-miR-939 is significantly associated with shorter distant metastasis-free survival of colon cancer and has a prognostic value in the development of colon cancer [33]. It is demonstrated that complement C3 (C3) is locally synthesized and secreted in the human intestine [34]. It is first found that the C3 mRNA expression is significantly upregulated in the inflamed mucosa of patients with inflammatory bowel disease [35]. It is worth mentioning that the expression of C3 is downregulated in the distal colon of the rat with constipation and decreased after the laxative treatment [36]. Intercellular adhesion molecule 1 (ICAM1) is found to be involved in the colon cancer cell proliferation [37]. In this study, we found upregulated expression of both C3 and ICAM in STC. In addition, the PPI network of HLA-DRB1, HLA-DRB5, C3, and ICAM1 showed that these genes had a significant association with proteins encoded by other differentially expressed genes. Our result suggested that HLA-DRB1, HLA-DRB5, C3, and ICAM1 may play a proinflammatory role in the development of STC.
In addition, we found that C3 and ICAM were commonly regulated by hsa-miR-1273 g-3p in STC. Hsa-miR-1273 has been regarded as a functionally uncharacteristic miRNAs in colon cancer exosomes [38]. In addition, another two upregulated DEGs including F2R like thrombin or trypsin receptor 3 (F2RL3) and syndecan 2 (SDC2) were also target genes of hsa-miR-1273 g-3p in STC. F2RL3 (also called PAR4) is a protease that plays roles in antinociceptive effects in visceral sensitivity through protease-activated receptor [39]. The expression of F2RL3 is downregulated in the colon of patients with irritable bowel syndrome patients [40,41], while, it is significantly upregulated in patients with colon cancer and ulcerative colitis [42,43]. SDC2 functions as a cell surface receptor and plays roles in interacting with extracellular matrix components and cell-to-cell signaling molecules [44]. It is found that SDC2 is low methylated in peripheral blood DNA of colon cancer and serves as a bloodbased diagnostic marker of colon cancer [45]. In a word, our result suggested that the expression of C3, ICAM1, F2RL3, and SDC2 under the regulation of hsa-miR-1273 g-3p may be important molecules in the process of STC.
In the PPI network, we found several high degree proteins encoded by DEGs, such as zinc finger and BTB domain containing 16 (ZBTB16), fibrillin 1 (FBN1), cyclin F (CCNF), and cyclin dependent kinase 1 (CDK1). Promyelocytic leukemia zinc finger protein (PLZF), encoded by ZBTB16 gene, is involved in different growth regulatory and differentiation pathways. It is significantly upregulated from normal colonic mucosa to aberrant crypt foci and colon cancer [46]. FBN1 is an extracellular matrix glycoprotein and high methylation in colonic tissue of patients with colon cancer [45]. The genetic mutation of FBN1 is also associated with colon cancer [47]. It is noted that FBN1 has been identified as a promising early detection biomarker of colon cancer and colon adenomas [48]. In this study, the expressions of ZBTB16 and FBN1 were upregulated and under the regulation of hsa-miR-4473 in STC. Our result suggested that both ZBTB16 and FBN1 could be involved in the growth regulatory of STC.
CCNF contributes to the G2 to M phase transition in the cell cycle. It is a key gene in colorectal adenoma-tocancer progression [49]. It is reported that CCNF is expressed in colon cancer tissue samples and overexpressed as an epigenetic clock gene in colon adenoma tissue samples [50,51]. It is reported that the interaction between CDK1/cyclin B1 complex and ZNFX1 antisense RNA 1 (ZFAS1) leads   to cell cycle progression and cell apoptosis suppression in the progression of colon cancer [52]. It is found that CDK1 is remarkably downregulated in colon cancer cell lines [53]. Significantly, it is regarded as a promising prognostic biomarker of the metastasis risk in stage II colon cancer [54]. In this study, CCNF and CDK1 were downregulated and regulated by hsa-miR-424-5p in STC. It is found that hsa-miR-424 is significantly upregulated in colon cancer [55]. In addition, hsa-miR-424 is a promising biomarker for the early diagnosis of colon cancer [56]. Beside CCNF and CDK1, prostaglandin D2 receptor (PTGDR) is also the target gene of hsa-miR-424 in STC. It is found that the expression of PTGDR is lower in colon cancer [57]. It is proposed that hypermethylation of PTGDR may make a contribution to reducing gene expression and developing colon cancer [58].
In a word, our study suggested that CCNF, CDK1, and PTGDR may be involved in the cell cycle progression of STC.
It is reported that hsa-miR-2116 presents remarkable differential expression when there are responses to bacterial infection in human host cells [59]. Herein, we found downregulated expression of hsa-miR-2116-3p in the colon tissue of STC. Furthermore, sphingosine-1-phosphate receptor 1 (S1PR1) and serpin family E member 1 (SERPINE1) were also target genes of hsa-miR-2116-3p. In mice, genetic deletion of S1pr1 will increase vascular permeability in the colon and enhance bleeding in colitis [60]. Additionally, it is found that targeting SphK1/S1P/S1PR1 may be a potential therapeutic option to prevent the development from colitis to colon cancer [61]. SERPINE1 is an inhibitor of plasmin action and has been indicated to promote cancer invasion [62]. It is found that SERPINE1 is upregulated in colon cancer and has been considered as a biomarker associated with poor prognosis in colon cancer [63]. Hamfjord et al. first found the expression of hsa-miR-3622a-5p in colon cancer tissues [64]. In this study, we found downregulated expression of hsa-miR-3622a-5p in colon tissue of STC. And adrenoceptor alpha (ADRA1D) was one of target genes of hsa-miR-3622a-5p in STC. It is found that ADRA1D is differentially expressed in stage IV compared to stages II and III colon cancer [65]. It is reported that hsa-miR-98 is upregulated in active ulcerative colitis [66]. Furthermore, the upregulation of hsa-miR-98 could be considered as the molecular marker for the colon in active ulcerative colitis patients [66]. Our result indicated that these hsa-miRNAs and target genes may play an important role in the development of STC. Beside the abovementioned hsa-miRNAs, hsa-miR-5100 and hsa-miR-4792 were also found in the colon tissue of STC. However, there were no previous reports about the relationship between hsa-miR-5100, hsa-miR-4792, and intestine disease including STC. Therefore, further research is needed to study the function of hsa-miR-5100 and hsa-miR-4792 in STC.

Conclusion
In summary, we identified 464 DEGs and 10 differentially expressed miRNAs in STC compared to normal tissues. PPI network generated some hub proteins such as ZBTB16, FBN1, CCNF, and CDK1. Target DEGs of differentially expressed miRNAs were significantly enriched in staphylococcus aureus infection involving four DEGs (HLA-DRB1, HLA-DRB5, C3, and ICAM), which suggested that inflammation may be a factor for STC. Our study was the first article for exploration of gene expression profiling in colon tissue of STC. Our findings may provide the fundamental work for understanding the molecular pathology of STC. However, there are limitations to our study. Firstly, the sample size in the RNA-sequencing was small and large numbers of samples of STC are needed for further research. Secondly, these dysregulated genes, miRNAs, and related signaling pathways in STC were identified and the potential molecular mechanism was not studied. Therefore, in vivo and in vitro experiments including animal model and cell culture were necessary to elucidate the underlying pathological mechanism of STC in the future work. Thirdly, qRT-PCR is needed to further validate the expression of several differentially expressed mRNAs and miRNAs.

Data Availability
The data used to support the findings of this study are included within the article.