Bin Zhao1†
Xiusheng Qu2†Xin Lv3Qingdong Wang4Deqiang Bian5Fan Yang1Xingwang Zhao1Zhiwu Ji1Jian Ni1Yan Fu1Guorong Xin1Haitao Yu1*- 1Department of Proctology, First Affiliated Hospital of Jiamusi University, Jiamusi, China
- 2Department of Chemoradiotherapy, First Affiliated Hospital of Jiamusi University, Jiamusi, China
- 3Department of General Surgery, Samii Medical Center, Shenzhen, China
- 4Department of Anesthesiology, First Affiliated Hospital of Jiamusi University, Jiamusi, China
- 5Scientific Research Departments, First Affiliated Hospital of Jiamusi University, Jiamusi, China
Non-coding RNAs such as long non-coding RNAs (lncRNAs) and microRNAs (miRNAs) have been found to be indispensable factors in carcinogenesis and cancer development. Numerous studies have explored the regulatory functions of these molecules and identified the synergistic interactions among lncRNAs or miRNAs, while those between lncRNAs and miRNAs remain to be investigated. In this study, we constructed and characterized an lncRNA–miRNA synergistic network following a four-step approach by integrating the regulatory pairs and expression profiles. The synergistic interactions with more shared regulatory mRNAs were found to have higher interactional intensity. Through the analysis of nodes in the network, we found that lncRNAs played roles that are more central and had similar synergistic interactions with their neighbors when compared with miRNAs. In addition, known colon adenocarcinoma (COAD)-related RNAs were found to be enriched in this synergistic network, with higher degrees, betweenness, and closeness. Finally, we proposed a risk score model to predict the clinical outcome for COAD patients based on two prognostic hub lncRNAs, MEG3 and ZEB1-AS1. Moreover, the hierarchical networks of these two lncRNAs could contribute to the understanding of the biological mechanism of tumorigenesis. For each lncRNA–miRNA interaction in the hub-related subnetwork and two hierarchical networks, we performed RNAup method to evaluate their binding energy. Our results identified two important lncRNAs with prognostic roles in colon cancer and dissected their regulatory mechanism involving synergistic interaction with miRNAs.
Introduction
Colon adenocarcinoma (COAD) has emerged as one of the leading causes of cancer-related deaths worldwide, with an increasing prevalence (Dienstmann et al., 2015). COAD is a complex disease, whose initiation and progression is closely related with mechanisms of regulation of gene expression (Liu et al., 2017). In recent years, with the development of next-generation sequencing technologies, studies suggested that less than 2% of the human genome encoded protein-coding genes, whereas non-coding RNAs represented most of the human transcriptome (Tay et al., 2014; Levy and Myers, 2016). Non-coding transcripts are divided into various classes, among which long non-coding RNAs (lncRNAs) and microRNAs (miRNAs) have attracted increasing attention. Notably, lncRNAs have been implicated in a diverse range of biological processes, including proliferation, migration, or genomic stability (Mercer et al., 2009; Fatica and Bozzoni, 2014). For instance, some lncRNAs have been reported to regulate gene expression through binding to PRC2 and acting as important controllers of cellular functions (Lee, 2012). miRNAs are small non-coding RNAs that also play a key role in gene regulation (Xu et al., 2019). miRNAs regulate gene expression mainly by binding to the 3′-untranslated regions of mRNAs and leading to their degradation or inhibition, and various studies have demonstrated that aberrant gene expression is closely linked to tumorigenesis, metastasis, and specific tumor stages (Liu et al., 2009).
Previous studies have demonstrated that lncRNAs interact with miRNAs to act on biological traits (Zhang et al., 2019): for example, Kallen et al. (2013) found that the lncRNA H19 modulated miRNA let-7 by performing in vivo crosslinking combined with affinity purification experiments. In summary, lncRNAs and miRNAs could interact by regulating mRNAs, thus playing critical roles in the pathological processes involved in disease development (Liao et al., 2019). However, the biological roles and functions of lncRNA–miRNA synergistic interactions have not yet been described in COAD and should be investigated to improve the efficiency of early diagnosis and treatment in the tumorigenesis and progression of this disease.
In this study, we constructed and characterized the lncRNA–miRNA synergistic network involved in COAD by integrating the lncRNA/miRNA–mRNA regulation pairs and the expression profiles of these RNA molecules. In total, we identified 305 positive and 294 negative synergistic lncRNA–miRNA interactions with significantly shared mRNAs. We observed that some of the synergistic lncRNAs and miRNAs were significantly enriched with cancer RNAs, and COAD-related lncRNAs were more important than COAD-related miRNAs. Finally, we proposed a risk score model to predict the clinical outcome of COAD patients based on two prognostic hub lncRNAs, MEG3 and ZEB1-AS1, which were identified by univariate Cox regression analysis. The biological mechanism involving these two lncRNAs was further analyzed. For synergistic lncRNA–miRNA interactions in the hub-related subnetwork and two prognostic lncRNAs related interactions, we provided the total free energy of binding evaluated by RNAup method. Altogether, our analysis provides new insight for exploring the molecular mechanisms of lncRNA–miRNA synergistic interactions and uncovering candidate lncRNA biomarkers for the prognosis of COAD.
Materials and Methods
The lncRNA/miRNA–mRNA Regulation Pairs
The miRNA–mRNA target data were downloaded and filtered from StarBase (Li et al., 2014). We chose the miRNA–mRNA interactions which were predicted by at least three of seven target-predicting programs. These seven target-predicting programs included PITA, RNA22, miRmap, DIANA-microT, miRanda, PicTar, and TargetScan. Recent advances in high-throughput sequencing of immunoprecipitated RNAs after cross-linking (CLIP-Seq, HITS-CLIP, PAR-CLIP, CLASH, iCLIP) provide powerful ways to identify biologically relevant RNA-target interactions. To obtain the high-quality miRNA–mRNA datasets, we further selected the miRNA–mRNA interactions which were validated by at least three CLIP-seq data from above interactions as the final miRNA–mRNA interactions. Similarly, we obtained the lncRNA–mRNA interactions that have at least two interactions and supported by at least three CLIP-seq data (Supplementary Table S1). Moreover, we downloaded the experimentally validated lncRNA–mRNA interactions from LncReg and LncRNA2Target databases (Supplementary Table S1; Zhou et al., 2015; Cheng et al., 2019). Integrating the lncRNA–mRNA interactions downloaded from these three databases, we obtained the final lncRNA–mRNA interactions. This way, we obtained 1,336 lncRNA–mRNA and 202,712 miRNA–mRNA high-quality non-redundant interactions.
Expression Profiles and Clinical Data of COAD Samples
The lncRNA, miRNA, and mRNA expression data for COAD patients were downloaded from the Cancer Genome Atlas (TCGA) database (Tomczak et al., 2015). For each expression profile, RNAs with missing values in more than 30% samples were removed and each of the remaining missing value was imputed by the KNN Imputation. Then all expression values were log2 transformed to obtain the final expression profiles. We chose sample-matched miRNA and lncRNA/mRNA expression profiles for further analysis. The clinical data of COAD patients was also obtained from TCGA.
CpG sites with missing values in more than 30% samples were removed and each of the remaining missing value was imputed by the KNN Imputation.
Collection of COAD-Related lncRNAs and miRNAs
The COAD-related lncRNAs were downloaded from LncRNADisease (Bao et al., 2019) and lnc2Cancer (Gao et al., 2019). Similarly, we collected the COAD-related miRNAs from several databases, including miR2Disease (Jiang et al., 2009), HMDD (Huang et al., 2019), SM2miR (Liu et al., 2013), and OncomiRDB (Wang et al., 2014).
Construction of the lncRNA–miRNA Synergistic Interaction Network
To identify the synergistic lncRNA–miRNA interactions, we developed a four-step computational method by integrating lncRNA–mRNA interactions, miRNA–mRNA interactions, and expression profiles of lncRNA, miRNA, and mRNA in COAD samples (Supplementary Figure S1).
First, high-quality lncRNA–mRNA and miRNA–mRNA interactions were downloaded from several databases and processed to obtain the non-redundant data as described above. Second, the regulatory networks of lncRNA–mRNA and miRNA–mRNA were constructed by filtering the lncRNA/miRNA–mRNA pairs obtained from StarBase with their expression profiles. The regulatory correlation between lncRNA/miRNA and mRNA was evaluated by Pearson correlation coefficient based on the matched lncRNA/miRNA and mRNA expression profiles. The pairs with significant correlation were saved as the lncRNA–mRNA (p-adjusted < 0.05) and miRNA–mRNA (R < -0.4, p-adjusted < 0.05) interactions in COAD patients. Third, we identified the co-regulated lncRNA–miRNA pairs by evaluating the significance of their shared regulated mRNAs (hypergeometric-test, p-adjusted < 0.01). Fourth, we used Pearson correlation to evaluate the synergistic direction and synergistic power of each co-regulated lncRNA–mRNA pair, and the pairs with p < 0.05 were considered the final synergistic lncRNA–mRNA pairs. Finally, after assembling the synergistic lncRNA–mRNA pairs, we obtained the lncRNA–miRNA synergistic interaction network. Two types of nodes were involved in the network (lncRNAs and miRNAs), with positive or negative synergistic regulations.
Survival Analysis According to the Risk Score Model
We evaluated the clinical outcomes for COAD patients by our risk score model based on the expression levels of MEG3 and ZEB1-AS1. We first divided the COAD samples into training (70% of the samples) and test (30% of the samples) sets. The risk score model was constructed by considering the individual power of MEG3 and ZEB1-AS1 evaluated by the univariable Cox regression analysis and their expression levels in training samples as follows:
where expi is the expression level of MEG3 or ZEB1-AS1 and coefi is the regression coefficient of MEG3 or ZEB1-AS1 estimated by univariate Cox regression analysis. As a result, the risk score for each COAD patient was computed by the formula:
The median risk score value in the training samples was chosen as the cut-off value to classify patients into high-risk and low-risk groups from the training and test sets, respectively. Survival analyses were performed to assess the difference in clinical outcome between the high-risk and low-risk groups, and statistical significance was evaluated by a log-rank test using the R package “survival.”
Network Visualization
The networks were visualized by Cytoscape 3.3.0 (Shannon et al., 2003), including the synergistic lncRNA–miRNA network, the hub-related subnetwork and the MEG3/ZEB1-AS1-related hierarchical networks.
Results
Construction and Characterization of the Synergistic lncRNA–miRNA Network
Based on the lncRNA/miRNA–mRNA regulation pairs downloaded from databases and expression profiles of lncRNA, miRNA, and mRNA, we constructed a synergistic lncRNA–miRNA network in the context of COAD (Figure 1A). As described in the ‘Materials and Methods’ section, we constructed this network in four steps. In the first step, we obtained 455 lncRNA–mRNA and 28,639 miRNA–mRNA COAD-specific regulation pairs (Table 1). Then, we identified 1,368 co-regulated lncRNA–miRNA pairs with significantly shared mRNAs by a hypergeometric test (Table 1). Finally, Pearson correlation analysis was used to filter the co-regulated lncRNA–miRNA pairs, and we obtained 305 positive and 294 negative synergistic lncRNA–miRNA interactions, covering 88 lncRNA and 161 miRNAs (Table 1, Figure 1B and Supplementary Table S2).
Figure 1. The synergistic long non-coding RNAs (lncRNA)-microRNAs (miRNA) network in the context of colon adenocarcinoma (COAD). (A) Global view of the synergistic network. Orange and purple nodes represent lncRNAs and miRNAs, respectively. Red and green edges represent positive and negative synergistic relationships. Larger nodes represent higher degrees. (B) Statistics of nodes and edges in the synergistic network. (C) Degree distribution of the synergistic network. (D) The mean path length of the synergistic network is higher than randomization test. The arrow represents the mean path length in the real network. (E) The synergistic lncRNA–miRNA pairs with more co-regulated mRNAs show more significantly synergistic regulatory relationship. The white circles represents the outlier points in boxplots.
Table 1. The statistic of COAD specific regulation pairs.
To explore the architecture and features of the synergistic network, its topological properties were analyzed. Through the analysis of node degree distribution, we found that the majority of nodes had few synergistic interactions, while a small portion had many interactions. This fits with the power-law distribution, suggesting that the synergistic lncRNA–miRNA network was scale-free and different from randomly generated networks (R2 = 0.82, Figure 1C). Moreover, we compared the clustering coefficient between our synergistic network and random networks. The result showed that the lncRNAs and miRNAs in our network had tight synergistic interactions (p < 0.001, Figure 1D). The synergistic lncRNA–miRNA pairs share varying numbers of mRNAs. In order to explore the relationship between their synergistic intensity and the number of shared mRNAs, we compared the co-expression significance with different numbers of shared mRNAs. The result showed that those lncRNA–miRNA pairs with more shared regulated mRNAs tended to be more significantly co-expressed to achieve coordinated regulation (R = 0.97, p = 0.026, Figure 1E).
lncRNAs Are More Likely to Have Similar Synergistic Interactions With Its Neighbors
The lncRNAs and miRNAs in our network had positive and negative synergistic interactions. Therefore, we counted and computed the ratio of positive and negative miRNA neighbors for each lncRNA. The results indicated that each lncRNA had 1∼55 lncRNA neighbors, including 0∼50 positive and 0∼26 negative neighbors (Figure 2A). After calculating the neighbor ratio of each lncRNA, we found that 82.95% of the lncRNAs tended to have the same synergistic direction as most (≥ 80%) of its miRNA neighbors. Among these lncRNAs, 48 and 52% were likely to have positive and negative synergistic interactions with miRNAs, respectively (Figure 2B). Examples of such lncRNAs are MALAT1, DANCR, and AGAP2-AS1 (Figure 2C).
Figure 2. The synergistic interactions between long non-coding RNAs (lncRNAs) and microRNAs (miRNAs). (A) The number of positive and negative miRNA neighbors for each lncRNA. Dark red and dark green represent positive and negative interactions, respectively. (B) The ratio of positive and negative miRNA neighbors for each lncRNA. (C) The synergistic interactions for three lncRNAs: MALAT1, DANCR, and AGAP2-AS1. (D) The number of positive and negative lncRNA neighbors for each miRNA. (E) The ratio of positive and negative lncRNA neighbors for each miRNA. (F) The synergistic interactions for two miRNAs: hsa-miR-106b-5p and hsa-miR-20a-5p.
Similarly, we counted and computed the neighbors of each miRNA, and discovered that each miRNA had 1∼14 lncRNA neighbors with 0∼10 positive and 0∼7 negative lncRNA neighbors (Figure 2D). As opposed to lncRNAs, only 43.48% of the miRNAs had the same synergistic direction as most (≥ 80%) of its lncRNA neighbors, while other miRNAs such as hsa-miR-106b-5p and hsa-miR-20a-5p had mixed synergistic relationships with their lncRNA neighbors (Figures 2E,F).
The Cancer-Related lncRNAs in the Synergistic Network Have Centralized Roles
Many cancer-related genes have been discovered, and to explore the relationship between our synergistic network and COAD, we obtained COAD-related lncRNAs and miRNAs from databases as described in the ‘Materials and Methods’ section. In total, we obtained 116 COAD-related lncRNAs and 134 miRNAs. After mapping these known COAD-related RNAs to our synergistic network (Figure 3A), we found that the synergistic lncRNAs and miRNAs were significantly enriched with the cancer-related RNAs (p < 2.2 × 10–16 for lncRNAs and p = 6.36 × 10–10 for miRNAs, Figure 3B). Then, we compared the topological properties of cancer-related RNAs with other RNAs in the synergistic network. The results showed that COAD-related genes had significantly higher degrees, betweenness centrality, and closeness centrality than other nodes in the synergistic network (p = 7.27 × 10–10, 4.46 × 10–06, and 3.77 × 10–05 for COAD-related lncRNAs; p = 1.67 × 10–03, 6.28 × 10–03, and 3.65 × 10–04 for COAD-related miRNAs; Figures 3C–E). These results suggest that cancer-related lncRNAs and miRNAs have more centralized roles when compared to other RNAs. Moreover, COAD-related lncRNAs appear to be more important than COAD-related miRNAs.
Figure 3. Known colon adenocarcinoma (COAD) related RNAs in synergistic network. (A) Known COAD related RNAs were mapping to the synergistic network. Node color and edge color are same as Figure 1A. Nodes with dark red border represent the known COAD related RNAs. (B) Enrichment of known COAD related long non-coding RNAs (lncRNAs) and microRNAs (miRNAs) in network. P-values are computed by hypergeometric test. (C) COAD related lncRNAs and miRNAs have higher degrees than other nodes. (D) COAD related lncRNAs and miRNAs have higher betweenness than other nodes. (E) COAD related lncRNAs and miRNAs have higher closeness than other nodes. P-values are calculated based on Wilcoxon rank sum test.
Identification of Two Central and Prognostic lncRNAs
Numerous studies have reported hub genes as playing key roles in cancer (Chen et al., 2017; Das et al., 2017; Yin et al., 2019). To identify the hub lncRNAs and miRNAs in our synergistic network, we computed a degree for each node and sorted them in a descending order. Then, we chose 10% of the nodes with the highest degrees as the hub RNAs, including 18 lncRNAs and seven miRNAs (Figure 4A). When comparing the ratio of hub lncRNAs among all lncRNAs with the ratio of hub miRNAs among all miRNAs, we found that a greater ratio of lncRNAs were identified as the hub RNAs (20% vs. 4%, Figure 4B), suggesting important roles of lncRNAs in the synergistic network and in accordance with our previous results. Moreover, we found that the majority of the hub RNAs were known COAD-related RNAs, and this ratio was higher than that in non-hub RNAs (80% vs. 56%, Figure 4C). The COAD-related hub RNAs included 18 lncRNAs and seven miRNAs. Next, we extracted the edges that connected two hub RNAs and found that all hub RNAs were connected, except for one lncRNA: UCA1. The hub subnetwork is depicted in Figure 4D. In accordance with the synergistic network described above, we found that the lncRNAs and miRNAs in the hub subnetwork were likely to interact with their neighbors in similar and different directions, respectively (Figure 4D).
Figure 4. Identification and characterization of hub RNAs. (A) Degree of hub RNAs. Orange bars and purple bars represent degree of long non-coding RNAs (lncRNAs) and microRNAs (miRNAs), respectively. (B) The ratio of hub and non-hub RNAs in lncRNAs and miRNAs, respectively. (C) The ratio of known colon adenocarcinoma (COAD) related lncRNAs and mRNAs in hub and non-hub RNAs. (D) Hub subnetwork extracted from the synergistic network. Color of node and edge are same as Figure 3A. (E) The result of univariable Cox regression for each hub lncRNA. The prognostic lncRNAs are marked in red. (F) The result of univariable Cox regression for each hub miRNA.
Considering that hub RNAs have key roles in the development of cancer, in our next step we evaluated the prognostic ability of each hub RNA. Through univariate Cox regression analysis, we identified two risk lncRNAs in COAD patients, MEG3 and ZEB1-AS1 (HR = 1.52 and p = 0.03 for MEG3, HR = 3.19 and p = 2.13 × 10–05 for ZEB1-AS1, Figures 4E,F). Furthermore, we combined these two lncRNAs using the risk score model to predict the clinical outcome of COAD patients. The risk score model was constructed according to Equation 1, with 0.4933 as the regression coefficient for MEG3 estimated by the univariate Cox regression analysis, and 1.1077 as the regression coefficient for ZEB1-AS1. The median risk score of training samples was used as the cut-off value (1.22) to separate the high-risk and low-risk groups. Survival analysis revealed a significant difference in overall survival between these two groups (log-rank test p = 0.00255, Figure 5A). Furthermore, we computed the risk score of samples in the test set based on the risk score model and divided these samples into high-risk and low-risk groups. Comparing the clinical outcome of samples between the two groups, we found that low-risk samples in the test set also showed significantly better prognosis (log-rank test p = 0.021, Figure 5B). In addition, we randomly selected different sample sets (70∼90% of all COAD samples and the 70∼90% of the test samples), computed their risk scores, and compared the survival difference between high-risk and low-risk groups. The results showed that the risk score model could predict the clinical outcome of COAD patients (Supplementary Figure S2). This result illustrated the robust prognostic ability of the MEG3 and ZEB1-AS1 lncRNAs.
Figure 5. (A) Kaplan–Meier curves of overall survival in high-risk and low-risk groups in training set. (B) Kaplan–Meier curves of overall survival in high-risk and low-risk groups in test set. P-values were calculated by the log-rank test. (C) The colored boxes and circles represents expression values of MEG3 and ZEB1-AS1 in cancer and normal samples. ‘can’ represent ‘cancer’ and ‘nor’ represent ‘normal’. (D) The colored boxes and circles represents expression values of MEG3 and ZEB1-AS1 in high-risk and low-risk samples. ‘hr’ represent ‘high-risk’ and ‘lr’ represent ‘low-risk’. P-values are calculated by Wilcoxon rank sum test.
To explore the roles of MEG3 and ZEB1-AS1 during cancer progress, we compared the expression values of MEG3 and ZEB1-AS1 in cancer samples and normal samples. As a result, we found that MEG3 presented similar expression levels in both contexts (p = 0.66), while ZEB1-AS1 showed a significantly higher expression in cancer samples (p = 3.983 × 10–09, Figure 5C). In addition, we compared the expression levels of these two lncRNAs in high-risk and low-risk samples: results showed that MEG3 and ZEB1-AS1 had significantly higher expression values in high-risk samples (all p < 2.2 × 10–16, Figure 5D). We also downloaded two GEO lncRNA expression datasets associated with colon diseases, including GSE77013 and GSE67106 (Mirza et al., 2015; Padua et al., 2016). The expression of MEG3 and ZEB1-AS1 were found to be highly expressed in disease samples compared with control samples (t-test, P = 0.09, 0.016, and 1.23e-08, Supplementary Figure S3). Due to the lack of ZEB1-AS1 probes in GSE77013, we didn’t compare its expression values between disease and control samples. These results demonstrate that the high expression of these two lncRNAs is associated with colon disease. Summarizing these results, we believe that MEG3 played a role in cancer development while ZEB1-AS1 could act in both carcinogenesis and cancer development. These results were consistent with previous studies (Fu and Cui, 2017; Gong et al., 2017; Liu et al., 2019; Hu et al., 2020; Ni et al., 2020).
Hierarchical Networks for Elucidating the Biological Mechanism of MEG3 and ZEB1-AS1
To contribute to the understanding of the synergistic interactions of MEG3 and ZEB1-AS1, we proposed hierarchical models to systematically illustrate the regulatory process. As shown in Figure 6, MEG3 or ZEB1-AS1 regulate mRNAs by synergistic interactions with miRNAs, and further participate in cancer biological processes. Among the miRNAs that have synergistic interactions with MEG3, 14 out of 20 miRNAs were known COAD-related miRNAs. These miRNAs, along with MEG3, further regulate 11 mRNAs, including CASP3, CASP8, and vascular endothelial growth factor (VEGFA). Previous studies have indicated that polymorphisms in CASP3 and CASP8 are related to colon cancer (Goodman et al., 2006; Choi et al., 2012), and the VEGFA was also significantly associated with rectal cancer (Slattery et al., 2014). In this study, we uncovered their upstream regulatory factors, which were MEG3 and its synergistic miRNAs. This result might be another indication of the roles of these mRNAs in carcinogenesis. Through the integration of annotation information on mRNAs, we found that MEG3 and its synergistic miRNAs were mainly associated with cancer-related processes such as immune system development, cell development, tissue development, cell differentiation, protein metabolism, and other processes.
Figure 6. Two hierarchical models of prognostic long non-coding RNAs (lncRNAs). (A) Hierarchical model of MEG3 related synergistic regulatory interactions. (B) Hierarchical model of ZEB1-AS1 related synergistic regulatory interactions. The models are exhibited hierarchically by three levels, including synergistic lncRNA–microRNAs (miRNA) interactions, regulatory mRNAs and annotated gene ontology (GO) biological processes of mRNAs. GO terms annotated to at least seven mRNAs and five mRNAs are showed in (A,B), respectively. Orange, purple and green nodes represent lncRNAs, miRNAs and mRNAs, respectively. Gray nodes represent the GO terms.
Regarding the synergistic miRNAs related to ZEB1-AS1, 15 out of 20 were found to be known COAD-related miRNAs. ZEB1-AS1 and these miRNAs synergistically regulated six genes, such as CyclinD1 (CCND1) and zinc finger E-box binding homeobox 1 (ZEB1). CCND1 is a key cell cycle regulatory protein and its polymorphism has been found to be significantly associated with overall COAD risk (Xie et al., 2017). Li et al. (2012) reported that IL-1β may promote colon tumor growth and invasion through the activation of colon cancer stem cell self-renewal and epithelial-mesenchymal transition (EMT), and ZEB1 plays a critical role in these two processes. Moreover, we observed that the mRNAs regulated by ZEB1-AS1 and its synergistic miRNAs were annotated to cancer-related gene ontology (GO) terms such as cell differentiation, cell proliferation, and developmental processes.
In light of these results, we believe that MEG3 and ZEB1-AS1 play important roles in the initiation and progression of colon cancer, through their synergistic interactions with cancer-related miRNAs and finally regulating cancer-related mRNAs that were associated with cancer biological processes. Our results could contribute to the understanding of important roles of synergistic lncRNA–miRNA interactions in tumorigenesis, expand the complexity of the ncRNA–mRNA regulatory network, and provide potential therapeutic targets for colon cancer treatment.
Discussion
Colon adenocarcinoma is the third most common cancer worldwide and has become a widespread health issue for its highly mortality and morbidity (Bertotti et al., 2015; Marmol et al., 2017). Recent studies suggested that interactions between lncRNAs and miRNAs in the regulation of mRNA expression played important regulatory roles in the initiation and progression of COAD (Yu et al., 2017). However, the regulatory mechanisms through which lncRNA–miRNA interactions are involved in the progression of this disease are stil unclear.
Long non-coding RNAs–miRNA synergistic interactions are critical for many biological functions and exploring these interactions contributes to a further understanding of the process of tumorigenesis and development of COAD (Guil and Esteller, 2015; Wang et al., 2019a). More importantly, increasing evidence shows that the lncRNA–miRNA interaction network is implicated in several pathophysiological processes, including human cancers (Lin et al., 2019). In this work, we constructed and characterized the lncRNA–miRNA synergistic network by integrating lncRNA–mRNA interactions, miRNA–mRNA interactions, and the expression profiles of lncRNA, miRNA, and mRNA in COAD samples. The analysis of this synergistic network allowed the detection of complicated features and functions of RNA regulatory interactions and how lncRNAs and miRNAs could play regulatory roles in the tumorigenesis and progression of COAD (Wang et al., 2020). Our results indicated that the synergistic lncRNAs and miRNAs were significantly enriched with cancer-related RNAs. In addition, COAD-related lncRNAs and miRNAs had significantly higher degrees, betweenness centrality, and closeness centrality than other nodes in the synergistic network. Further analysis showed that cancer-related miRNAs, especially lncRNAs, had more centralized roles when compared with other RNAs. Altogether, our study of lncRNA–miRNA interactions could contribute with crucial information in the understanding of the regulatory mechanisms through which ncRNAs act, as well as with the identification of molecular targets for the diagnosis and treatment of COAD.
Reliable prediction of RNA–RNA binding energies is crucial. RNAup is an effective method, which involved two energy contributions, including (1) the energy necessary to ‘open’ the binding site and (2) the energy gained from hybridization. To improve the medical effectiveness of our results, we performed RNAup to compute the potential binding possibility between lncRNAs and miRNAs. The sequence of lncRNA transcripts and miRNA were downloaded from Ensemble and miRBase database, respectively (Supplementary Table S3; Kozomara et al., 2019; Yates et al., 2020). For interactions in the hub-related subnetwork and two hierarchical networks, we provided the total free energy of binding for each lncRNA–miRNA interaction (Supplementary Tables S4, S5). Based on the total free energy of binding we provided, users can acquire both direct and indirect interactions by their own cutoffs. This way, we expect to provide results that have higher potential medical usefulness.
Accumulating evidence revealed that lncRNAs acted as prognostic biomarkers and regulated cell functions in colorectal cancer (Yin et al., 2015; Wang et al., 2019b). Through analyses of node degree and univariate Cox regression analysis, we identified two important lncRNAs: MEG3 and ZEB1-AS1. We further proposed two hierarchical models to systematically illustrate the regulatory process of these two lncRNAs. In the hierarchical models, most miRNAs which have synergistic interactions with MEG3 or ZEB1-AS1 were found to be known COAD-related miRNAs. Moreover, we found that some mRNAs regulated by the lncRNAs and miRNAs were reported to be associated with COAD. Our results proposed another indication of the roles of these mRNAs in carcinogenesis. We believe that other ncRNAs and mRNAs in the hierarchical models were also COAD-related RNAs. Our results provide potential therapeutic targets for colon cancer treatment. Finally, we proposed a risk score model to predict the clinical outcome of COAD patients and demonstrated the utility of lncRNAs as promising biomarkers.
Data Availability Statement
All datasets generated in this study are included in the article/Supplementary Material.
Author Contributions
HY conceived and designed the experiments. BZ, XQ, XL, QW, DB, FY, XZ, ZJ, JN, GX, and YF collected and analyzed data. HY and BZ wrote this manuscript. All authors read and approved the final manuscript.
Funding
This work was supported by scientific research project of Heilongjiang Health and Family Planning Commission [2014-242].
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2020.572983/full#supplementary-material
References
Bao, Z., Yang, Z., Huang, Z., Zhou, Y., Cui, Q., and Dong, D. (2019). LncRNADisease 2.0: an updated database of long non-coding RNA-associated diseases. Nucleic Acids Res. 47, D1034–D1037.
Bertotti, A., Papp, E., Jones, S., Adleff, V., Anagnostou, V., and Lupo, B. (2015). The genomic landscape of response to EGFR blockade in colorectal cancer. Nature 526, 263–267.
Chen, P., Wang, F., Feng, J., Zhou, R., Chang, Y., Liu, J., et al. (2017). Co-expression network analysis identified six hub genes in association with metastasis risk and prognosis in hepatocellular carcinoma. Oncotarget 8, 48948–48958. doi: 10.18632/oncotarget.16896
Cheng, L., Wang, P., Tian, R., Wang, S., Guo, Q., Luo, M., et al. (2019). LncRNA2Target v2.0: a comprehensive database for target genes of lncRNAs in human and mouse. Nucleic Acids Res. 47, D140–D144.
Choi, J. Y., Kim, J. G., Lee, Y. J., Chae, Y. S., Sohn, S. K., Moon, J. H., et al. (2012). Prognostic impact of polymorphisms in the CASPASE genes on survival of patients with colorectal cancer. Cancer Res. Treat. 44, 32–36. doi: 10.4143/crt.2012.44.1.32
Das, S., Meher, P. K., Rai, A., Bhar, L. M., and Mandal, B. N. (2017). Statistical approaches for gene selection, hub gene identification and module interaction in gene co-expression network analysis: an application to aluminum stress in soybean (Glycine max L.). PLoS One 12:e0169605. doi: 10.1371/journal.pone.0169605
Dienstmann, R., Salazar, R., and Tabernero, J. (2015). Personalizing colon cancer adjuvant therapy: selecting optimal treatments for individual patients. J. Clin. Oncol. 33, 1787–1796. doi: 10.1200/jco.2014.60.0213
Fatica, A., and Bozzoni, I. (2014). Long non-coding RNAs: new players in cell differentiation and development. Nat. Rev. Genet. 15, 7–21. doi: 10.1038/nrg3606
Fu, J., and Cui, Y. (2017). Long noncoding RNA ZEB1-AS1 expression predicts progression and poor prognosis of colorectal cancer. Int. J. Biol. Markers 32, e428–e433.
Gao, Y., Wang, P., Wang, Y., Ma, X., Zhi, H., Zhou, D., et al. (2019). Lnc2Cancer v2.0: updated database of experimentally supported long non-coding RNAs in human cancers. Nucleic Acids Res. 47, D1028–D1033.
Gong, H., Wen, H., Zhu, X., Lian, Y., Yang, X., Qian, Z., et al. (2017). High expression of long non-coding RNA ZEB1-AS1 promotes colorectal cancer cell proliferation partially by suppressing p15 expression. Tumour Biol. 39:1010428317705336.
Goodman, J. E., Mechanic, L. E., Luke, B. T., Ambs, S., Chanock, S., and Harris, C. C. (2006). Exploring SNP-SNP interactions and colon cancer risk using polymorphism interaction analysis. Int. J. Cancer 118, 1790–1797. doi: 10.1002/ijc.21523
Guil, S., and Esteller, M. (2015). RNA-RNA interactions in gene regulation: the coding and noncoding players. Trends Biochem. Sci. 40, 248–256. doi: 10.1016/j.tibs.2015.03.001
Hu, D., Zhang, B., Yu, M., Shi, W., and Zhang, L. (2020). Identification of prognostic biomarkers and drug target prediction for colon cancer according to a competitive endogenous RNA network. Mol. Med. Rep. 22, 620–632. doi: 10.3892/mmr.2020.11171
Huang, Z., Shi, J., Gao, Y., Cui, C., Zhang, S., Li, J., et al. (2019). HMDD v3.0: a database for experimentally supported human microRNA-disease associations. Nucleic Acids Res. 47, D1013–D1017.
Jiang, Q., Wang, Y., Hao, Y., Juan, L., Teng, M., Zhang, X., et al. (2009). miR2Disease: a manually curated database for microRNA deregulation in human disease. Nucleic Acids Res. 37, D98–D104.
Kallen, A. N., Zhou, X. B., Xu, J., Qiao, C., Ma, J., Yan, L., et al. (2013). The imprinted H19 lncRNA antagonizes let-7 microRNAs. Mol. Cell. 52, 101–112. doi: 10.1016/j.molcel.2013.08.027
Kozomara, A., Birgaoanu, M., and Griffiths-Jones, S. (2019). miRBase: from microRNA sequences to function. Nucleic Acids Res. 47, D155–D162.
Lee, J. T. (2012). Epigenetic regulation by long noncoding RNAs. Science 338, 1435–1439. doi: 10.1126/science.1231776
Levy, S. E., and Myers, R. M. (2016). Advancements in next-generation sequencing. Annu. Rev. Genomics Hum. Genet. 17, 95–115.
Li, J. H., Liu, S., Zhou, H., Qu, L. H., and Yang, J. H. (2014). starBase v2.0: decoding miRNA-ceRNA, miRNA-ncRNA and protein-RNA interaction networks from large-scale CLIP-Seq data. Nucleic Acids Res. 42, D92–D97.
Li, Y., Wang, L., Pappan, L., Galliher-Beckley, A., and Shi, J. (2012). IL-1beta promotes stemness and invasiveness of colon cancer cells through Zeb1 activation. Mol. Cancer 11:87. doi: 10.1186/1476-4598-11-87
Liao, J., Wang, J., Liu, Y., Li, J., and Duan, L. (2019). Transcriptome sequencing of lncRNA, miRNA, mRNA and interaction network constructing in coronary heart disease. BMC Med. Genomics 12:124. doi: 10.1186/s12920-019-0570-z
Lin, C., Yuan, G., Hu, Z., Zeng, Y., Qiu, X., Yu, H., et al. (2019). Bioinformatics analysis of the interactions among lncRNA, miRNA and mRNA expression, genetic mutations and epigenetic modifications in hepatocellular carcinoma. Mol. Med. Rep. 19, 1356–1364.
Liu, B., Li, J., Tsykin, A., Liu, L., Gaur, A. B., and Goodall, G. J. (2009). Exploring complex miRNA-mRNA interactions with Bayesian networks by splitting-averaging strategy. BMC Bioinformatics 10:408. doi: 10.1186/1471-2105-10-408
Liu, H., Ye, D., Chen, A., Tan, D., Zhang, W., Jiang, W., et al. (2019). A pilot study of new promising non-coding RNA diagnostic biomarkers for early-stage colorectal cancers. Clin. Chem. Lab. Med. 57, 1073–1083. doi: 10.1515/cclm-2019-0052
Liu, R., Zhang, W., Liu, Z. Q., and Zhou, H. H. (2017). Associating transcriptional modules with colon cancer survival through weighted gene co-expression network analysis. BMC Genomics 18:361. doi: 10.1186/s12864-017-3761-z
Liu, X., Wang, S., Meng, F., Wang, J., Zhang, Y., Dai, E., et al. (2013). SM2miR: a database of the experimentally validated small molecules’ effects on microRNA expression. Bioinformatics 29, 409–411. doi: 10.1093/bioinformatics/bts698
Marmol, I., Sanchez-De-Diego, C., Pradilla Dieste, A., Cerrada, E., and Rodriguez Yoldi, M. J. (2017). Colorectal carcinoma: a general overview and future perspectives in colorectal cancer. Int. J. Mol. Sci. 18:197. doi: 10.3390/ijms18010197
Mercer, T. R., Dinger, M. E., and Mattick, J. S. (2009). Long non-coding RNAs: insights into functions. Nat. Rev. Genet. 10, 155–159.
Mirza, A. H., Berthelsen, C. H., Seemann, S. E., Pan, X., Frederiksen, K. S., Vilien, M., et al. (2015). Transcriptomic landscape of lncRNAs in inflammatory bowel disease. Genome Med. 7:39.
Ni, X., Ding, Y., Yuan, H., Shao, J., Yan, Y., Guo, R., et al. (2020). Long non-coding RNA ZEB1-AS1 promotes colon adenocarcinoma malignant progression via miR-455-3p/PAK2 axis. Cell Prolif. 53:e12723.
Padua, D., Mahurkar-Joshi, S., Law, I. K., Polytarchou, C., Vu, J. P., Pisegna, J. R., et al. (2016). A long noncoding RNA signature for ulcerative colitis identifies IFNG-AS1 as an enhancer of inflammation. Am. J. Physiol. Gastrointest. Liver Physiol. 311, G446–G457.
Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13, 2498–2504. doi: 10.1101/gr.1239303
Slattery, M. L., Lundgreen, A., and Wolff, R. K. (2014). VEGFA, FLT1, KDR and colorectal cancer: assessment of disease risk, tumor molecular phenotype, and survival. Mol. Carcinog. 53(Suppl. 1), E140–E150.
Tay, Y., Rinn, J., and Pandolfi, P. P. (2014). The multilayered complexity of ceRNA crosstalk and competition. Nature 505, 344–352. doi: 10.1038/nature12986
Tomczak, K., Czerwinska, P., and Wiznerowicz, M. (2015). The Cancer Genome Atlas (TCGA): an immeasurable source of knowledge. Contemp. Oncol. 19, A68–A77.
Wang, D., Gu, J., Wang, T., and Ding, Z. (2014). OncomiRDB: a database for the experimentally verified oncogenic and tumor-suppressive microRNAs. Bioinformatics 30, 2237–2238. doi: 10.1093/bioinformatics/btu155
Wang, W., Lou, W., Ding, B., Yang, B., Lu, H., Kong, Q., et al. (2019a). A novel mRNA-miRNA-lncRNA competing endogenous RNA triple sub-network associated with prognosis of pancreatic cancer. Aging 11, 2610–2627. doi: 10.18632/aging.101933
Wang, W., Xie, Y., Chen, F., Liu, X., Zhong, L. L., Wang, H. Q., et al. (2019b). LncRNA MEG3 acts a biomarker and regulates cell functions by targeting ADAR1 in colorectal cancer. World J. Gastroenterol. 25, 3972–3984. doi: 10.3748/wjg.v25.i29.3972
Wang, Y., He, R., and Ma, L. (2020). Characterization of lncRNA-associated ceRNA network to reveal potential prognostic biomarkers in lung adenocarcinoma. Front. Bioeng. Biotechnol. 8:266. doi: 10.3389/fbioe.2020.00266
Xie, M., Zhao, F., Zou, X., Jin, S., and Xiong, S. (2017). The association between CCND1 G870A polymorphism and colorectal cancer risk: a meta-analysis. Medicine 96:e8269. doi: 10.1097/md.0000000000008269
Xu, J., Meng, Q., Li, X., Yang, H., Xu, J., Gao, N., et al. (2019). Long Noncoding RNA MIR17HG Promotes Colorectal Cancer Progression via miR-17-5p. Cancer Res. 79, 4882–4895. doi: 10.1158/0008-5472.can-18-3880
Yates, A. D., Achuthan, P., Akanni, W., Allen, J., Allen, J., Alvarez-Jarreta, J., et al. (2020). Ensembl 2020. Nucleic Acids Res. 48, D682–D688.
Yin, D. D., Liu, Z. J., Zhang, E., Kong, R., Zhang, Z. H., and Guo, R. H. (2015). Decreased expression of long noncoding RNA MEG3 affects cell proliferation and predicts a poor prognosis in patients with colorectal cancer. Tumour Biol. 36, 4851–4859. doi: 10.1007/s13277-015-3139-2
Yin, K., Zhang, Y., Zhang, S., Bao, Y., Guo, J., Zhang, G., et al. (2019). Using weighted gene co-expression network analysis to identify key modules and hub genes in tongue squamous cell carcinoma. Medicine 98:e17100. doi: 10.1097/MD.0000000000017100
Yu, Y., Nangia-Makker, P., Farhana, L., and Majumdar, A. P. N. (2017). A novel mechanism of lncRNA and miRNA interaction: CCAT2 regulates miR-145 expression by suppressing its maturation process in colon cancer cells. Mol. Cancer 16:155. doi: 10.1186/s12943-017-0725-5
Zhang, W., Tang, G., Zhou, S., and Niu, Y. (2019). LncRNA-miRNA interaction prediction through sequence-derived linear neighborhood propagation method with information combination. BMC Genomics 20:946. doi: 10.1186/s12864-019-6284-y
Keywords: non-coding RNA, synergistic interaction, colon cancer, biological mechanism, prognostic biomarker
Citation: Zhao B, Qu X, Lv X, Wang Q, Bian D, Yang F, Zhao X, Ji Z, Ni J, Fu Y, Xin G and Yu H (2020) Construction and Characterization of a Synergistic lncRNA–miRNA Network Reveals a Crucial and Prognostic Role of lncRNAs in Colon Cancer. Front. Genet. 11:572983. doi: 10.3389/fgene.2020.572983
Received: 15 June 2020; Accepted: 24 August 2020;
Published: 15 September 2020.
Edited by:
Xiaofeng Dai, Jiangnan University, ChinaReviewed by:
Padhmanand Sudhakar, Earlham Institute (EI), United KingdomNathan Weinstein, Universidad Nacional Autónoma de México, Mexico
Copyright © 2020 Zhao, Qu, Lv, Wang, Bian, Yang, Zhao, Ji, Ni, Fu, Xin and Yu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Haitao Yu, MTU2OTQ1NDk4OThAMTYzLmNvbQ==
†These authors have contributed equally to this work