Profiling the lncRNA-miRNA-mRNA ceRNA network to reveal potential crosstalk between inflammatory bowel disease and colorectal cancer

Introduction

Inflammatory bowel disease (IBD) includes ulcerative colitis (UC) and Crohn’s disease (CD), and is a chronic and incurable inflammatory disorder most frequently affecting the gastrointestinal tract. The main pathological changes of IBD involve intestinal fibrosis, abscesses, fistulas, and ulcers (Baumgart & Sandborn, 2012; Danese & Fiocchi, 2011). According to a report by Kappelman in 2007 in the United States, the prevalence of CD and UC in children <20 years of age was 43 and 28 per 100,000, respectively, while among adults the prevalence was 201 and 238 per 100,000, respectively (Kappelman et al., 2007). In Europe, the prevalence of IBD ranged from 3–7 cases for CD and 4–11 cases for UC per 100,000 (Vegh, Kurti & Lakatos, 2017). The prevalence has also been reported to vary widely in Asian countries (6.5–121 per 100,000) (Ng et al., 2013; Thia et al., 2008). The etiology and pathogenesis of IBD still remain unclear, and abnormal immune regulation, gut commensal microbes, gastrointestinal environmental factors, and genetic susceptibility are commonly cited as key factors (Khor, Gardet & Xavier, 2011). A recent study has also reported that IBD is still a global disease with an accelerating incidence in newly industrialized countries (Ng et al., 2018). Patients with IBD need a lifelong course of drugs, which results in immense financial and psychological burdens to patients and their families.

IBD is a significant risk factor for colorectal cancer (CRC). In 2001, by summarizing 116 studies of the relationships between IBD and CRC, Eaden, Abrams & Mayberry (2001) concluded that the overall prevalence of CRC was 3.7% in UC patients. Another meta-analysis reported the prevalence of UC-associated CRC as 0.85%, which was elevated to nearly 14% at 30 years in those with long-standing extensive colitis (Bopanna et al., 2017). Comparable findings have been demonstrated in CD and the reported incidence was 8% at 22 years. Similarly, the cumulative risk of CRC in CD was reported as 0.3%, 1.6%, and 2.4% at 5, 15, and 25 years after diagnosis (Shawki et al., 2018). Furthermore, IBD-associated CRC patients exhibit differences both in chemotherapeutic and clinical outcomes (Axelrad et al., 2017).

The competing endogenous RNA (ceRNA) hypothesis provides a better understanding of tumorigenesis (Salmena et al., 2011; Wee et al., 2012). LncRNAs contain miRNA-response elements, which function as ceRNAs, interacting with miRNAs to indirectly regulate mRNAs. The molecular pathogenesis of IBD-associated CRC has been widely reported in recent studies. Kameyama et al. (2018) reported similar genomic alterations between samples obtained from CRC and colitis-associated colorectal cancer, except for some distinct variations in the times and frequencies. Although an expanding number of miRNAs and lncRNAs have been shown to play a role in IBD-induced CRC, few of them reported comprehensive lncRNA-mRNA-miRNA ceRNA networks during disease development (Josse & Bours, 2016; Kanaan et al., 2012; Lewis et al., 2017; Svrcek et al., 2013; Zhu et al., 2016).

In the present study, we identified a potential ceRNA network and characterized the potential mechanism during IBD-related CRC using bioinformatic analysis. We examined differentially expressed genes (DEGs) among CRC (or IBD related neoplasia), IBD and normal tissue based on three previously published gene datasets, GSE10714, GSE4183 and GSE68306, available in the Gene Expression Omnibus (GEO) database. By taking full advantage of The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov), an authoritative dataset COAD was used to identify abnormally expressed mRNAs, lncRNAs, and miRNAs in CRC compared to normal tissue. Together with four datasets, the crosstalk network among aberrant expressions of lncRNA-miRNA-mRNA was constructed to predict candidate diagnostic biomarkers and therapeutic targets for IBD-induced CRC. Further functional analyses, including Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment, protein-protein interaction (PPI) network, transcription factor (TF) enrichment, and survival analysis were performed to explore potentially important genes involved in the transformation from IBD to CRC.

Materials & Methods

Expression profile dataset

Gene expression profiles for IBD-related CRC (GSE10714, GSE4183 and GSE68306) were obtained from the GEO database. GSE10714 and GSE4183 were normalized genome-wide gene expression profiles, evaluated by HGU133 Plus 2.0 microarrays. In GSE4183, the RNA expression profiles of 15 patients with CRC, 15 with IBD, and eight healthy normal controls were included. In GSE10714, seven patients with CRC, seven with IBD (four with CD and three with UC), and three healthy normal controls were chosen. Meanwhile, GSE68306 is a miRNA array, in which samples were grouped into normal tissue (n = 16), UC without neoplasia (UC, n = 9), and UC-associated neoplastic tissue (n = 11). The quality of normalized gene expression data was analyzed and visualized using the ggplot2 package of R software for each group and sample. Moreover, COAD RNA-seq dataset was analyzed online by GEPIA based on TCGA database which contained 463 tumor tissue and 85 normal tissue (two recurrent and one metastasis tissue have been excluded). COAD miRNA expression profile containing 451 tumor sample and eight normal samples (one recurrent and one metastasis samples have been excluded) was obtained from XENA (a portal website for TCGA) and analyzed in R. All of the following processes have been shown in flow chart (Fig. 1).

Results

Data distribution analyses and DEG screening

In GSE10714 and GSE4183, a total of 54,675 genes were detected in 17 samples. The expression values and the distributions were similar between the three groups in each dataset (Figs. S1A, and S1D for GSE10714; Figs. S1B, and S1E for GSE4183). These data distribution analyses showed that the two GEO datasets were qualified for further bioinformatics analyses. In GSE68306 miRNA array dataset, a total of 654 miRNAs were detected in 36 samples. Among them, samples Normal 12-16, IBD 8-9, and Neoplasia 9-11 were ruled out because data distribution of these samples were quite different from others in the same group (Fig. S1F). The expression density of three groups in GSE68306 is similar after filtration (Fig. S1C). Following data processing using LIMMA, we identified 2,179 DEGs, 1,291 DEGs, and 74 DE-miRNAs in the CRC/Neoplasia group vs. the IBD group for GSE10714 (Fig. 2A), GSE4183 (Fig. 2B) and GSE68306 (Fig. 2C), respectively; meanwhile, 1765 DEGs, 2265 DEGs and 72 DE-miRNAs were identified in the IBD group vs. the Normal group for GSE10714 (Fig. 2E), GSE4183 (Fig. 2F) and GSE68306 (Fig. 2G), respectively. According to GEPIA, 5,372 DEGs in the CRC group vs. the Normal group were obtained with the same cut-off criteria. Also, 341 DE-miRNAs in the CRC group vs. the Normal group were identified by R based on COAD miRNA expression profile (Fig. 2D). Repetitive genes with different probes and these predicted genes without official gene symbol or transcript annotation in NCBI or Ensembl were removed in the following analysis.

Construction of the ceRNA network

The lncRNA-miRNA-mRNA network was then constructed following the steps listed in the Materials and Methods, which contained three lncRNA, 15 miRNAs, and 138 mRNA (Fig. 3A). CeRNA is said to act as a sponge for miRNAs and free mRNA from miRNA binding. Therefore, a core ceRNA network was extracted including synergic expressed mRNA and lncRNA but divergent expressed miRNAs, and shown in Fig. 3B. These three lncRNAs in Fig. 3A were termed key ceRNAs. The interactions and nodes information in ceRNA network were shown in Tables S10 and S11, respectively.

Discussion

In the current study, we used bioinformatic analysis to identify IBD-related CRC mRNAs and lncRNAs based on the published database obtained from the GEO and TCGA. By constructing a lncRNA-miRNA-mRNA network, a total of three lncRNAs, 15 miRNAs, and 138 mRNAs were identified as potential mediators in the pathophysiological processes of IBD-related CRC. Moreover, we also performed functional analyses to investigate the potential biological functions and pathways involved in the disease progression from IBD to CRC.

According to our GO annotation results, DEGs included in ceRNA network, could regulate the IBD-related CRC through involvement in the regulation of the immune response, response to lipopolysaccharide (LPS), and inflammation. The immune response was shown to play an important role in the pathogenesis of IBD. LPS was well known as a TLR4 activator (Chowdhury, Sacks & Sheerin, 2006). As an important mediator involved in the host immune response, TLR4 is overexpressed in inflammation-associated colorectal neoplasia tissues, and contributes to the production of inflammatory factors in chronic colitis. Moreover, TLR4-deficient mice exhibit lower rates of colonic polyps (Fukata et al., 2007). Together, these results showed a critical role for the TLR4-mediated immune response in the tumorigenesis of colitis-associated CRC. Conversely, Lowe et al. reported that mice with colitis-associated colorectal cancer deficient in TLR2 easily developed significantly more and larger colorectal tumors compared with wild-type mice, together with increased expressions of IL-6, IL-17A, and phospho-STAT3 (Lowe et al., 2010). Inflammation was also widely reported in IBD-related CRC. Necrosis and regeneration of intestinal mucosa were the main pathological changes of IBD. Longer-duration chronic inflammation in IBD could lead to the abnormal regulation of the immune system (Kuraishy, Karin & Grivennikov, 2011). A number of proinflammatory pathways have been recently identified in inflammation-associated tumor development (Long, Lundsmith & Hamilton, 2017). Specific molecular signaling pathways such as NF-κB, has been implicated in the carcinogenesis of IBD-related CRC. NF-κB has been identified as a key upstream TF in regulating DEGs from ceRNA and NF-κB related pathway, a downstream part of LPS/TLR4, is highly enriched and popped up in KEGG pathway enrichment analysis. Deletion of IKKβ in intestinal epithelial cells decreased the incidence of colitis-associated cancer in the colitis-associated mouse cancer model (Greten et al., 2004). IL-6/STAT3 was also demonstrated to be involved in metabolic reprogramming, to promote the progression from chronic colitis to CRC (Qu et al., 2017), and the COX-2/PGE2 pathway was also shown to contribute to inflammation and carcinogenesis in IBD. Overexpression of COX-2 was observed in UC-associated neoplasia at both the gene and protein levels (Agoff et al., 2000). Peroxisome proliferator-activated receptor δ (PPARδ) was shown to exert an important role in chronic inflammation and pathophysiological processes of CRC via the COX-2/PGE2 pathway (Wang & DuBois, 2010; Wang et al., 2014). In vivo studies showed that PPAR δ was necessary for the production of proinflammatory mediators, such as chemokines, cytokines, COX-2, and PGE2, which led to colitis-associated tumor growth (Wang et al., 2014). In general, these immune response and inflammatory response-related hub genes identified in our study, including CXCL10, CSF1, CXCL1, CXCL9, IL1R1, IL1RN, CXCL11, FCGR1A, CXCL2 and GBP2 (centrality degree >10), deserve further research and discussion. Interestingly, all these genes are downregulated in CRC vs. IBD group, while they are all higher expressed in CRC and IBD compared to normal tissue. It seems that they are dynamically regulated according to the necessity of tumorigenesis process. Further research need be done to indicate the potential function of these genes in IBD related CRC.

We also identified three key ceRNAs, DPP10-AS1, DIO3OS and LINC01106, which were not ascertained in previous studies. Among them, LINC01106 is the only ceRNA related to overall survival of CRC patients. LIN01106 contained a transcript of 2204 bp with two exons located on chromosome 2q13. However, few studies have been conducted on the function of LINC01106. Our results showed intimate interactions between LINC01106 and miRNAs, such as MIR651, MIR193A, MIR140, MIR214, MIR34A, and MIR130A. Among them, MIR193A is a popular miRNA in scientific research and stands out for its divergent expression change in CRC vs. IBD relative to LINC01106 as well as its mRNA targets in the core ceRNA network. MIR193A, located on 17q11.2, has been widely found to be dysregulated in tumor progression. Lower expression of MIR193A was observed in human aldosterone-producing adrenocortical adenomas, and was functioned as a tumor suppressor by targeting and restraining CYP11B2 expression both at mRNA and protein level (Zhang et al., 2018). Another group demonstrated that MIR193A played an important role in the regulation of gastric cancer. Both gastric cell lines and human gastric tumors contained decreased levels of MIR193A. In addition, an in vitro study showed that upregulating MIR193A decreased tumor cell proliferation and migration, at least partly, by directly targeting cyclin D1 (CCND1) and ETS proto-oncogene 1 (ETS1) expression (Chou et al., 2018). Moreover, LINC00152, a lncRNA overexpressed in gastric cancer, was proved to sponge MIR193A to increase MCL1 expression thereby promoting cells proliferation (Huang et al., 2018). In ovarian cancer, MIR193A is silenced by methylated modification and down regulation of MIR193A expression promotes tumor aggressiveness by losing function of targeting GRB7 and MAPK/ERK pathways (Chen et al., 2018). Similar results were gained by another research that MIR193A is downregulated in metastatic prostate cancer and is involved in HOTAIR/EZH2/miR-193a feedback loop to exert a tumor suppressive function in prostate cancer (Ling et al., 2017). But, contradictorily, MIR193A increases the ability of chemoresistance in prostate cancer cell. MIR193A is noticed to be upregulated in prostate cancer tissues and cell lines, with significant suppression of cell apoptosis induced by oxidative stress by targeting Bach2 (Yang et al., 2017). The role of MIR193a in colorectal cancer and inflammatory bowel disease has also been broadly explored. In 2014, researchers demonstrated that MIR93A is downregulated in CRC to promote lymph node metastasis and result in poor survival (Zhang et al., 2014). In addition, MIR193A was evidenced to target and down regulate colonic PepT1, which reduced intestinal inflammation in response to microbiota (Dai et al., 2015). Another study confirms the tumor suppressor roles of MIR193A by perturbing KRAS function in colorectal adenocarcinoma, including reducing cell proliferation, increasing apoptosis, and inhibiting epithelial-mesenchymal transition (Mamoori et al., 2018). Moreover, MIR193A is revealed to be downregulated in ulcerative colitis neoplasia and promote tumorigenesis by releasing its control for IL17RD expression (Pekow et al., 2017). Considering the important roles of miRNAs in the development of CRC, we propose that our future research should therefore focus on the LINC01106-associated pathway and potential miRNA targets, especially MIR193A.

Our current study is based on published GEO datasets and TCGA COAD dataset. A limitation of the research is that we failed to enroll relevant mRNA, lncRNA, and microRNA by comparing IBD-induced CRC and IBD directly since it’s hard to find qualified data. We are trying to collect such kind of samples from clinical patients and will do the validation research in the future.

Conclusions

The present study identified a series of IBD-related mRNAs and lncRNAs that were potential mediators in the disease development from IBD to CRC. Using a series of bioinformatics analyses, lncRNA LINC01106, DPP10-AS1, and DIO3OS were identified, which were strongly associated with IBD-related CRC. By constructing a ceRNA network with lncRNAs, miRNAs, and mRNAs, we identified the links between these three RNA species. Further functional analyses suggested a potential mechanism through which these ceRNAs mediated tumorigenesis. Our results emphasized the important role of ceRNAs in the pathogenesis of IBD-related CRC, although more research is necessary to confirm our findings.

Supplemental Information