Construction of ceRNA Network to Reveal Potential Biomarkers in Crohn’s Disease and Validation in a TNBS Induced Mice Model

Introduction

Crohn’s disease (CD) and ulcerative colitis are chronic inflammatory disorders of the gastrointestinal tract, with symptoms evolving in a relapsing and remitting manner that comprise the term inflammatory bowel disease (IBD).¹ CD is characterized by the involvement of all parts of the intestine, the most common being the terminal ileum and colon.^2,3 The complex interplay between genetic susceptibility, environmental factors, and altered gut microbiota, leading to dysregulated innate and adaptive immune responses, might result in CD.⁴ Although the mortality caused by CD is very low, the morbidity is still a serious problem, and a large number of CD patients do not respond to drug treatments and must undergo multiple operations to relieve symptoms.⁵ In addition, CD is incurable and increases the risk of lymphoma, biliary tract cancer, and colorectal cancer.⁶ In the past 10 years, CD has developed into a global disease, the prevalence on every continent is rising, adding a heavier economic burden to society. The incidence of CD in China, despie being lower than that in western countries, has increased rapidly with urbanization and industrialization. The current estimated incidence rate is 0.51–1.09 per 100,000.^7–10 However, the pathogenesis and pathophysiology of CD are not yet fully understood.¹¹ Therefore, it is crucial to explore the pathogenic mechanisms of CD.

In recent years, a new hypothesis of RNA is “competing endogenous RNA (ceRNA)”, which proposes the regulatory dialogue between different RNAs, including long noncoding RNAs (lncRNAs), microRNAs (miRNAs), transcribed pseudogenes, and circular RNAs.¹² To date, lncRNAs have been proved to be involved in many cellular processes and pathological conditions.¹³ lncRNAs can serve as miRNA sponges and compete for reducing the ability of miRNA binding to protein-coding transcripts.^14,15 At present, many studies have shown that ceRNA were implicated in different kinds of diseases, such as cancer, acute kidney injury, Alzheimer’s disease, and acute myocardial infarction.^16–18 Emerging data have suggested that ceRNA networks play a pivotal role in progression in IBD.^19,20

In this study, we obtained the data related to CD from the NCBI GEO, and compared an expression profile between CD and healthy controls. Following, mRNAs, miRNAs, LncRNAs of CD and normal samples were applied to differential expression analysis. Simultaneously, weighted correlation network analysis (WGCNA) was used to enrich modules to select mRNAs and lncRNAs most related to CD. Furthermore, miRWalk and LncBase databases were used to predict interacting mRNAs and lncRNAs. Thus, four lncRNAs, four miRNAs, and eight mRNAs were identified to construct a lncRNA–miRNA–mRNA ceRNA network in CD. Immune infiltration analysis was performed to investigate the effect of candidate mRNA biomarkers in regulation of immune cells. Furthermore, the expressions of candidates were validated in a 2,4,6-trinitrobenzene sulfonic acid (TNBS) induced colitis mice model. Finally, the predicted ceRNA network was confirmed via RNA interference (RNAi) in the HCT116 cell line. Collectively, our results provided insights into the molecular mechanism and investigated potential biomarkers of CD.

Materials and Methods

GEO Dataset Selection

A flowchart of the study design is shown in Figure 1. The expression profile datasets of lncRNAs, miRNAs, and mRNAs were downloaded from the NCBI GEO (www.ncbi.nlm.nih.gov/geo) database. lncRNAs, miRNAs, and mRNAs between active human CD tissue and matched normal tissue were included. The lncRNA expression profile data were obtained from GSE67106 (20 active CD tissues and 21 matched normal tissues from intestinal biopsies).²¹ miRNAs and mRNAs expression profiling data were obtained from GSE102134 (composed of GSE102127 and GSE 102133), GSE102127 (33 active CD tissues and six matched normal tissues from ileum biopsies) and GSE102133 (55 active CD tissues and 12 normal tissues from ileum biopsies) datasets.²² The validation set was GSE75214 (75 active CD tissues and 22 normal tissues from ileum biopsies).²³

Figure 1 Flowchart of this work.

Differential Expression Analysis and WGCNA

Differential expression analysis was performed by GEO2R, with an adjusted P-value<0.05, the log(fold-change) >1.0 or <−1.0 was defined as a DEmRNA, DEmiRNA, DElncRNA, respectively.

WGCNA was applied to assess the relative importance and module membership of lncRNAs and mRNAs.²⁴ The minimum number of module genes was set at 30. The hierarchical clustering dendrogram summarized the Gene modules with different colors.

Gene Set Enrichment Analysis and Protein–Protein Interaction Network Building

The metascape database was used for functional annotation of genes. Gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) pathway enrichment analyses were performed.²⁵

A protein–protein interaction (PPI) network was constructed by the STRING database with high confidence (0.7), and line color indicates the type of interaction evidence.²⁶ The hub genes were identified by Cytoscape cytohubba.

CeRNA Network Construction

The miRWalk database was used to predict the interaction between miRNA and mRNA.²⁷ The LncBase database was applied to explore the interaction between miRNA and lncRNA.²⁸ The intersection of DEmRNA, WGCNA module, and predict mRNA was selected to construct the ceRNA network, the same as the intersection of DElncRNA, WGCNA module, and predict lncRNA and DEmiRNA. The CeRNA network was visualized by Cytoscape.

Immune Infiltration Analysis

ImmuCellAI database was applied to estimate the differential of immune cell infiltration in two groups.²⁹

TNBS Induced Colitis Mice Model

TNBS colitis belongs to the group of chemically induced colitis animal models, since it was introduced in 1989. It is the most common experimental colitis model to closely mimic CD, which shares features of human CD. Twenty specific pathogen-free (SPF) female C57BL/6J mice (8–10 weeks old, weighing 20∼22 g) were purchased from Beijing Vital River Laboratory Animal Technology Co. Ltd. All animal experiments were conducted in strict accordance with the Institutional Animal Care and Use Committees and Institutional Review Board (IRB) of Tongji Hospital. All experimental procedures of the mice were specifically approved for this study by the Medical Ethics Committee of Tongji Hospital and were carried out in accordance with institutional guidelines (TJ-IRB20182677) and guide for the care and use of laboratory animals (National Academies Press, 2011). All experimental procedures of the mice were specifically approved for this study by the Medical Ethics Committee of Tongji Hospital. All procedures on mice were performed under sodium pentobarbital anesthesia, and all efforts were made to minimize suffering. The TNBS induced mice model most commonly utilized animal models of CD and has shown a significant consistency that reflects their extensive use during the last decades.³⁰ Before the induction of experimental colitis, a short fasting is needed of about 12 hours. Subsequently, the anesthesia is performed via sodium pentobarbital, then administered with 100 mg/kg TNBS in 50% ethanol, while the control group was only administered 50% ethanol. After that, the mice were held in a head down position in order to avoid expulsion of the fluid for a few minutes. All mice were euthanized on the 4th day and the colons were collected for further tests.

RNAi and Transfection

Human colon cancer cell line HCT116 cells were seeded in 6-wells plates at a density of 1.5×10⁶ cells/well and incubated at 50–60% confluence before transfection. Then, the cells were transfected with siRNA-MIR3142HG or mir-429 inhibitor or negative control (NC) (GenePharma, China) by using transfection reagent GP-transfect-Mate (GenePharma, China) according to the manufacturer’s protocol. Oligonucleotide sequences are listed in Supplementary Table S2.

qRT-PCR

TRIpure Reagent was used to extract total RNA from HCT116 cells and colon tissues of TNBS mice, and the PCR conditions were 95°C for 5 minutes, followed by 40 cycles of 95°C for 15 seconds, 60°C for 15 seconds, and 72°C for 30 seconds. The relative gene expression level was calculated using the 2-ΔΔCt method. To normalize the data, GAPDH was used as an internal reference. The sequences of the primers are shown in Table 1.

Table 1 The Sequences of Primers

Results

Differential Gene Expression and Enrichment Analysis

The expression of mRNAs from the GEO dataset was explored. We found 248 significantly upregulated DEmRNAs and 162 significantly downregulated DEmRNAs from 55 CD tissues and 12 normal tissues. Figure 2A shows the distribution of all the DEmRNAs on the two dimensions of -log10 (false discovery rate, FDR) and log2 (fold change, FC) through a volcano plot. To explore the biological function of DEmRNAs, GO analysis were used, and the top 10 biological processes (BP), cellular components (CC), and molecular functions (MF) are shown in Figure 2B. These terms indicate that DEmRNAs of CD are associated with “response to bacterium”, “regulated exocytosis”, and ”leukocyte migration” in BP. The top three terms of CC were ”extracellular matrix”, ”apical part of cell”, and ”apical plasma membrane”. In categories of MF, “receptor regulator activity”, “anion transmembrane transporter activity”, and “organic anion transmembrane transporter activity” were major terms. Several inflammatory-associated BPs were observed, including “myeloid leukocyte activation”, “regulation of inflammatory response”, “response to molecule of bacterial origin”, and ”leukocyte chemotaxis”. KEGG pathway analysis suggested that DEmRNAs were mainly enriched in “cytokine–cytokine receptor interaction”, “protein digestion and absorption”, and “chemokine signaling pathway”, as shown in Figure 2C and D.

Figure 2 DEmRNA screening and enrichment analysis. (A) Volcano plot of DEmRNA. (B) Top 10 GO enrichment. (C) Top 10 KEGG enrichment. (D) KEGG chord plot of top 10 terms.

WGCNA Was Applied to Analyze Gene Modules

Gene modules were analyzed by using the WGCNA among the top 2,000 mRNAs by variance comparison. As shown in Figure 3A, R² cut at 0.8, while 8 was chosen as a soft threshold. Eight gene modules were identified, which were used to analyze module-trait (CD and normal) co-expression similarity and adjacency (Figure 3B). Yellow and turquoise modules are shown as strongly positive (module trait correlation=0.626) and negative (−0.561) correlation to CD, respectively (Figure 3C). Moreover, the miRWalk database was used to predict mRNA targeted by DEmiRNAs (Supplementary Figure S1 and Table 2). The overlap of DEmRNAs, WGCNA modules, and DEmiRNA-predict-mRNA were further analyzed (Figure 4A), and the heatmap of these 82 mRNAs are shown in Figure 4B.

Table 2 The DEmiRNAs

Figure 3 mRNA modules analyzed by WGCNA. (A) Determination of soft-thresholding power. (B) Clustering dendrogram. (C) The module-trait relationships.

Figure 4 Identification of candidate overlapping mRNA. (A) Intersection of DEmRNAs, WGCNA modules, and miRwalk-predict. (B) Heatmap of overlapping mRNA.

PPI Analysis

To explore the functions of these mRNAs, we constructed the PPI network of proteins encoded by these 82 mRNAs via STRING with high confidence (Figure 5A), line color indicated the type of interaction evidence. Total 41 mRNAs were filtered into the PPI network, which contained 41 nodes and 57 edges. Furthermore, hub genes with highest degrees of connectivity in the PPI network were identified by Cytohubba plugin of Cytoscape. As shown in Figure 5B, the top 10 hub genes were CXCL10, STAT1, IL1B, CXCL8, PTGS2, CXCL9, ICAM1, SOCS3, CXCL5, and IDO1. Then, GO and KEGG enrichment analyses were performed for these 82 mRNAs. GO enrichment is shown in Figure 5C and Supplementary Table S1, the relationships between mRNAs and GO terms is exhibited in Figure 5D. The top 10 KEGG terms are shown in Figure 5E, which indicates these genes are associated with protein digestion and absorption, the TNF signaling pathway, the NOD-like receptor signaling pathway, complement and coagulation cascades, the IL-17 signaling pathway, the NF-kappa B signaling pathway, and the Toll-like receptor signaling pathway. The relationships between mRNAs and enriched KEGG pathways are shown in Figure 5F.

Figure 5 The overlapping mRNA function annotation. (A) PPI network. (B) hub gene identified via cytohubba. (C) KEGG enrichment. (D) GO enrichment. (E) GO chord plot of top 10 term. (F) KEGG chord plot of top 10 term. Lollipop chart.

LncRNA Module Analyzed by WGCNA

Differential expression analysis of lncRNAs from the GEO dataset was performed. As shown in Figure 6A, 424 upregulated and 137 downregulated DElncRNAs were found. We analyzed the top 2,000 lncRNAs via WGCNA. R² was set as 0.8, while the soft threshold was 9 (Figure 6B). The resulting gene dendrograms and module colors are shown in Figure 6C, in total 13 gene modules were identified. As shown in Figure 6D, the turquoise module exhibited a strongly positive correlation (module trait correlation=0.848), while the green module shows a strongly negative correlation (module trait correlation=−0.781) with CD. The overlap of DElncRNAs, WGCNA modules, and DEmiRNA-predict-lncRNA were further used to construct the ceRNA network (Figure 6E and Table 3).

Table 3 The DElncRNAs

Figure 6 DElncRNA screening. (A) Volcano plot of DElncRNA. (B) Determination of soft-thresholding power. (C) Clustering dendrogram. (D) The module-trait relationships. (E) Intersection of DElncRNA, WGCNA, and LncBase-predict.

CeRNA Network Constructed

According to the above analyses, miRNA-mRNA and lncRNA-miRNA relationship pairs were revealed. To explore how lncRNA regulates mRNA by binding to miRNA in CD, based on relationship pairs of miRNA-mRNA and lncRNA-miRNA obtained, the lncRNA-miRNA-mRNA ceRNA network was constructed via Cytoscape. The upregulation and downregulation of these nodes were also determined, and the ceRNA network is shown in Figure 7A. According to competition rule in the ceRNA network, the final regulation network, which included four lncRNAs, four miRNAs, and eight mRNAs, was obtained (Figure 7B).

Figure 7 The construction of (A) ceRNA network, red color represents upregulate and blue represents downregulate, circle shape represents mRNAs, triangle represents miRNAs, and rhombus represents lncRNAs. (B) Sankey diagram of final ceRNA network.

Receiver Operating Characteristic (ROC) Analysis and Validation

ROC curves of the gene in the final regulation network are shown in Figure 8A and B. In GSE75214 (tissues from ileum biopsies), except CTSS, AUCs of those gene were 0.75–1. AUCs of LPL, STAT1, TRIB2, and PTGS2, even beyond 0.9, which suggest those genes might have outstanding value for predicting the responsiveness in CD patients.

Figure 8 Validation and immune infiltration analysis. (A) and (B) ROC of mRNAs in the ceRNA network. (C) The differential of immune infiltration between CD and normal controls. (D) Correlation matrix of immune cell proportions and PTGS2. (E) Lollipop chart of correlation between immune cell and PTGS2. *p<0.05, **p<0.01, ***p<0.001. Abbreviation: ns, not significant.

Immune Infiltration Analysis

The ImmuCellAI database was applied to assess the composition of immune cells in CD and normal tissues. As shown in Figure 8C, CD tissues showed a higher fraction of Exhausted T-cells, Type 1 regulatory T-cells, Monocytes, Natural Killer (NK) cells, and a lower fraction of effector memory T-cells, Gamma delta T-cells, and T-cells CD8. As shown in Figure 8D and E, in CD tissues, the expression of PTGS2 exhibited high correlation with critical immune cells in the pathogenesis of CD, like Natural regulatory T-cells, Type 1 regulatory T-cells, Th1 helper cells, and Th2 helper cells (p<0.01).

TNBS Induced Mice Model Validation

The expression of PTGS2 and mir-429 were measured by qRT-PCR in an experimental colitis mice model. Compared with the control group, PTGS2 was high expressed, while mmu-mir-429 was low expressed in the experimental colitis group (Figure 9A and B). The raw data are available in Table S3.

Figure 9 The relative expression validated via qRT-PCR. (A) The relative expression of PTGS2 in TNBS-induced mice and control mice. (B) The relative expression of mmu-mir-429 in TNBS-induced mice and control mice. (C) The relative expression of PTGS2 in HCT116 cells transfected with NC, si-MIR3142HG, and mir-429 inhibitor. ***p<0.001.

RNAi and Transfection

To validate the MIR3142HG/mir-429/PTGS2 ceRNA network, a series of experiments were designed and conducted by using siRNA and miRNA inhibitors. Results of qRT-PCR demonstrated that knockout of MIR3142HG decreased the mRNA level of PTGS2, whereas inhibition of mir-429 increased the mRNA level of PTGS2 in HCT116. Furthermore, the co-transfection of si-MIR3142HG and mir-429 inhibitor could eliminate effects (Figure 9C).

Discussion

CD is a chronic inflammatory disease of the gastrointestinal tract, due to the increasing incidence worldwide, understanding the molecular mechanism, identifying diagnostic biomarkers as well as treatment of CD are of urgent concern.⁴ Recent reports have described lncRNAs which act as ceRNAs or miRNA sponges by competing for binding to shared miRNAs that are important post-transcriptional regulators of gene expression.¹² CeRNAs regulatory network has been reported to be involved in development and progression of many kinds of diseases. Exploring key ceRNA networks will benefit in CD therapy. For instance, Lnc-ITSN1-2 serve as ceRNA for IL-23R via sponging mir-125a in IBD.¹⁹ In addition, circRNA_103765 functions as a ceRNA to promote DLL4-mediated cell apoptosis by decoying the miR-30 family in CD.³¹ Therefore, the construction of the ceRNA to understand regulation mechanisms and explore potential biomarkers is essential.

In this study, we identified DEmRNAs, DEmiRNAs, and DElncRNAs between CD and normal ileum tissues from GSE102134 and GSE67106, respectively. Then, WGCNA was performed to select enrich modules related to CD. Furthermore, DEmiRNAs were used to predict interacted mRNAs and lncRNAs in the miRwalk and LncBase database. The intersection of DEmRNAs, mRNAs in WGCNA module and predicted mRNAs defined as intersection-mRNAs, intersection-lncRNAs were obtained similarly like intersection-mRNAs. Finally, intersection-mRNAs, intersection-lncRNAs, and DEmiRNA were used to construct the ceRNA regulatory network. According to competition rule in the ceRNA network, a final regulation network, which included four lncRNAs, four miRNAs, and eight mRNAs was obtained. To access the possibility of eight mRNAs in ceRNA act as CD biomarkers, validation was carried out by ROC analysis in GSE75214 from ileum tissues. The AUCs of LPL, STAT1, TRIB2, and PTGS2, even beyond 0.9, suggested those genes would be promising biomarkers in CD.

PTGS2 encodes an enzyme, which converts arachidonic acid to prostaglandin endoperoxide H2 in humans. It is one of the key factors in the cellular response to the inflammation and epithelial regeneration.³² PTGS2 is involved in the production of prostaglandin associated with the mediation of inflammation and exacerbate inflammation.^33,34 It has been reported to be significantly upregulated in inflamed mucosa of CD patients, and resulted in activation of the prostaglandin D2 (PGD2) metabolic pathway. PGD2 is a major inflammatory mediator, which participated in inflammation response in CD.³⁵ The results of immune infiltration analysis indicated that PTGS2 might be involved in immunomodulatory CD, due to the high correlation with critical immune cells. Moreover, intestinal immune cell trafficking has been identified as a central event in the pathogenesis of CD.³⁶ The expression of PTGS2 shows a positive correlation with regulator T-cells and T-cells in CD4, which were reported to be associated with a more complicated disease course with strictures or penetration, and a trend was observed towards more abdominal surgeries.³⁷ The increase in the number of macrophages and monocytes secretes more pro-inflammatory factors, maintains or promotes the pathological environment, and forms a vicious circle of inflammation in CD patents, which is positively correlated with the expression of PTGS2. Consistent with a previous study, the fraction of T-cells CD8 was decreased in CD patients, which is a negative correlation with PTGS2.³⁸ These results indicated that PTGS2 can be the immune-related biomarker candidate of CD.

The expression of PTGS2 was regulated via the MIR3142HG/mir-429/PTGS2 ceRNA network constructed in this work. LncRNA MIR3142HG functions as a posttranscriptional modulator via interacting with mir-429, which was reported to regulate inflammation response.³⁹ A recent study showed that downregulated mir-429 induces myristoylated alanine-rich protein kinase C substrate upregulated, and may act as a protection mechanism to repair active inflammation and inhibit the cellular breakdown of the epithelial barrier in a mice model.⁴⁰ PTGS2 (AUC 0.927) plays a key role in the cellular response to inflammation, which is typically induced by inflammation and potentially impairs the reestablishment of a functional epithelial barrier in IBD.³²

TNBS induced colitis offers an alternative method to study the pathobiology of CD, since it resembles human CD. Therefore, we established a TNBS induced colitis model and validated the expression of mir-429 and PTGS2 in the mice model. The expression of PTGS2 was significantly upregulated while mir-429 was downregulated in the TNBS mice model, which were consistent with the previous bioinformatics analysis. Furthermore, to validate the MIR3142HG/mir-429/PTGS2 ceRNA network predicted in this study, we conducted RNAi. As previously analyzed, knockout of MIR3142HG decreased the mRNA level of PTGS2, whereas inhibition of mir-429 increased the mRNA level of PTGS2 in HCT116, while the co-transfection of si-MIR3142HG and mir-429 inhibitor could rescue the mRNA level of PTGS2. It demonstrated that the MIR3142HG/mir-429/PTGS2 ceRNA network may play a role in development of CD.

Conclusion

In summary, we identified a ceRNA network in CD, which included four lncRNAs, four miRNAs, and eight mRNAs. It will be of benefit in understanding the regulatory mechanism in CD. Furthermore, four mRNAs in this ceRNA network showed promising diagnostic accuracy (AUC >0.9), used to investigate the effect in regulation of immune cells. PTGS2 exhibited a high correlated with immune cells, which were related to the development of CD. These results reveal PTGS2 may be a potential diagnostic biomarker and the MIR3142HG/mir-429/PTGS2 ceRNA network was involved in the pathogenesis of CD. However, this study still has some limitations, and, owing to the lack of lncRNA and miRNA datasets of CD, the lncRNAs and miRNAs identified in the ceRNA network cannot be validated. We cannot measure the expression of MIR3142HG in the TNBS induced mice model, due to the lack of sequence of lncRNA MIR3142HG of mice in the NCBI database. Despite these drawbacks, our results still benefit in further understanding the pathogenesis and explore the immune-related diagnostic biomarker of CD.