Identification of Transcriptional Markers and microRNA–mRNA Regulatory Networks in Colon Cancer by Integrative Analysis of mRNA and microRNA Expression Profiles in Colon Tumor Stroma

The aberrant expression of microRNAs (miRNAs) and genes in tumor microenvironment (TME) has been associated with the pathogenesis of colon cancer. An integrative exploration of transcriptional markers (gene signatures) and miRNA–mRNA regulatory networks in colon tumor stroma (CTS) remains lacking. Using two datasets of mRNA and miRNA expression profiling in CTS, we identified differentially expressed miRNAs (DEmiRs) and differentially expressed genes (DEGs) between CTS and normal stroma. Furthermore, we identified the transcriptional markers which were both gene targets of DEmiRs and hub genes in the protein–protein interaction (PPI) network of DEGs. Moreover, we investigated the associations between the transcriptional markers and tumor immunity in colon cancer. We identified 17 upregulated and seven downregulated DEmiRs in CTS relative to normal stroma based on a miRNA expression profiling dataset. Pathway analysis revealed that the downregulated DEmiRs were significantly involved in 25 KEGG pathways (such as TGF-β, Wnt, cell adhesion molecules, and cytokine–cytokine receptor interaction), and the upregulated DEmiRs were involved in 10 pathways (such as extracellular matrix (ECM)-receptor interaction and proteoglycans in cancer). Moreover, we identified 460 DEGs in CTS versus normal stroma by a meta-analysis of two gene expression profiling datasets. Among them, eight upregulated DEGs were both hub genes in the PPI network of DEGs and target genes of the downregulated DEmiRs. We found that three of the eight DEGs were negative prognostic factors consistently in two colon cancer cohorts, including COL5A2, EDNRA, and OLR1. The identification of transcriptional markers and miRNA–mRNA regulatory networks in CTS may provide insights into the mechanism of tumor immune microenvironment regulation in colon cancer.


Background
The tumor microenvironment (TME) is complex "ecosystems" which comprises many different cell types including stromal cells [1]. The functions of stromal cells are often altered in the TME, and the tumor stroma is associated with tumor cellular migration, neoangiogenesis, immunosurveillance evasion, and drug resistance [2]. Colorectal cancer (CRC) is one of the leading causes of cancer mortality worldwide [3]. Transcriptional signatures of CRC stromal cells have been associated with poor prognosis in CRC [4][5][6]. MicroRNAs (miRNAs) play a critical role in modulating the TME in CRC [7]. Bullock et al. showed that stromal and epithelial miRNAs play key roles during CRC progression [8]. Hiyoshi et al. revealed that the elevated expression of miRNAs-34b and -34c in stromal tissues is associated with poor prognosis in colon cancer [9]. These prior studies demonstrated the significant contribution of aberrantly expressed miRNAs in colorectal tumor stroma to CRC growth, invasion, and metastasis.
The miRNA-mRNA interaction networks in cancers have been widely explored. For example, Taguchi et al. identified the potential miRNA-mRNA interactions within multiple cancer types, including colorectal/colon cancer, hepatocellular carcinoma, non-small cell lung cancer, esophageal squamous cell cancer, prostate cancer, and breast cancer [10]. Yang et al. constructed miRNA and mRNA integrative networks as key regulators in colon adenocarcinoma (COAD) [11]. Li et al. built microRNA-mRNA interactions in human colon cancer [12]. Paul et al. identified potential miRNA-mRNA regulatory modules in CRC based on miRNA and mRNA expression data [13]. Sells et al. showed that the tissue kallikrein-related peptidase 6 (KLK6) mediated specific microRNA-mRNA regulatory networks associated with colon cancer invasion [14]. Izadi et al. investigated the conserved mRNA-miRNA interactions in colon and lung cancers [15]. Besides CRC, other cancer-associated miRNA-mRNA interactions have been explored. For example, Pham et al. identified miRNA-mRNA regulatory relationships which were consistent across different breast cancer subtypes [16]. Lou et al. constructed an miRNA-mRNA regulatory network potentially associated with glioblastoma multiforme [17]. These prior studies demonstrated that the miRNA-mRNA regulatory networks play crucial roles in the pathogenesis of cancer. However, the investigation of miRNA-mRNA regulatory networks in the TME, such as tumor stroma, remains lacking.

Identification of DEGs and DEmiRs between CTS and Normal Stroma
We used Network Analyst [19] to identify the DEGs between CTS and normal stroma by a meta-analysis of two CTS gene expression profiling datasets (GSE31279 and GSE35602). Each dataset was normalized by quantile normalization. The batch effects from both datasets were removed using the ComBat method [35]. We used the R package "limma" to identify the DEGs between CTS and normal stroma, and the Cochran s combination test to perform the meta-analysis [36]. The false discovery rate (FDR) [37] was used to adjust for multiple tests. We identified the DEGs using a threshold of absolute combined effect size (ES) > 1. 24 and FDR < 0.05. We used miRNet [18] to identify DEmiRs between CTS and normal stroma by analyzing the preprocessed GSE35602-GPL8227. The dataset was normalized by quantile normalization, and the R package "limma" was utilized to identify the DEmiRs between CTS and normal stroma using a threshold of |log2(FC)| > 1.0 and adjusted p-value < 0.05. The "FC" denotes the fold change of mean expression levels between two compared groups.

Identification of Upstream Transcription Factors and Target Genes of DEmiRs
The upstream transcription factors (TFs) regulating DEmiRs were predicted using FunRich [28], a tool for analyzing functional enrichment and interaction network of genes and proteins. The most significant (with the smallest p-values) 10 upstream TFs for the upregulated and downregulated DEmiRs were identified, respectively. We used miRNet [18], DIANA-microT-CDS [38], and DIANA-TargetScan [22] to identify target genes of DEmiRs. The online tool "Calculate and draw custom Venn diagrams" (http://bioinformatics.psb.ugent.be/webtools/Venn/) was used to identify common genes between the target genes of DEmiRs and the DEGs between CTS and normal stroma.

KEGG Pathway Enrichment Analysis of DEmiRs
We identified KEGG [39] pathways that were significantly associated with the DEmiRs by DIANA-mirPath (version 3) [21]. In this tool, we chose to analyze the validated miRNA-gene interactions using archive in microT web server v5.0 [38]. A threshold of adjusted p-value < 0.05 (Fisher's exact test) was used to define the statistical significance.

Construction of Protein-Protein Interaction (PPI) Network of DEGs
We constructed a PPI network of the DEGs between CTS and normal stroma by STRING [26]. The hub genes (with no less than 10 edges connected to other nodes) in the PPI network were identified using the node explorer module of Network Analyst [19].

Comparisons of the Expression Levels of Hub Genes between Colon Cancer and Normal Tissue and Between Colon Tumor Stromal Fibroblast and Normal Fibroblast
We compared the expression levels of the hub genes in the PPI network which were also the targets of the DEmiRs between TCGA COAD and normal tissue, and between colon tumor stromal fibroblast and normal fibroblast. A threshold of Student s t test p-value < 0.05 and |FC| > 1 was utilized to denote the statistical significance.

Associations of the Expression Levels of the Hub Genes Targeted by DEmiRswith Survival Prognosis in Colon Cancer
We analyzed the associations of the expression levels of the hub genes targeted by DEmiRs with survival prognosis in colon cancer using PrognoScan [24] and GEPIA [23]. The log-rank test p-value < 0.05 denoted the statistical significance.

Associations of the Expression Levels of the Hub Genes Targeted by DEmiRs with Immune Signature Enrichment Levels in Colon Cancer
We analyzed the correlations of the expression levels of three hub genes (COL5A2, EDNRA, and OLR1) with the abundance of immune signatures in TCGA COAD using TIMER [25]. Moreover, we compared the ratios of pro-inflammatory cytokines to anti-inflammatory cytokines between the colon cancers with high expression levels of the hub genes (expression levels > median) and the colon cancers with low expression levels of the hub genes (expression levels < median) in the TCGA COAD dataset using Student s t test. We defined the ratio of pro-/anti-inflammatory cytokines in a tumor sample as the ratio of average expression levels (base-2 log transformed) of their marker genes. The pro-inflammatory cytokines represent the immune-stimulatory signature with marker genes IFNG, IL-1A, IL-1B, and IL-2, and the anti-inflammatory cytokines represent the immune-inhibitory signature with marker genes IL-4, IL-10, IL-11, and TGFB1 [40].

Identification of Food and Drug Administration (FDA)-Approved Drug-Hub Gene Interaction
We identified the drugs that target the hub genes using DGIdb [27]. DGIdb collects drug-gene interaction data from 30 disparate sources, including ChEMBL, DrugBank, Ensembl, NCBI Entrez, PharmGKB, PubChem, clinical trial databases, and literature in NCBI PubMed. The drug-gene interactions supported by at least one database and/or PubMed reference were identified. From the identified drug-gene interactions, we selected the drugs that have been approved by the FDA.

Identification of DEmiRs and Their Target DEGs
We identified seven downregulated and 17 upregulated DEmiRs in CTS relative to colon normal stroma (Table 1). A literature review showed that all seven miRNAs downregulated in CTS have been associated with colon cancer and other cancer types. For example, hsa-mir-135b-5pwas significantly downregulated in CTS with the highest fold change of expression levels (FDR = 0.009, FC > 7). This miRNA has been associated with gastric carcinogenesis [41]. Hsa-miR-214-3p was significantly downregulated in CTS (FDR = 0.02, FC > 3) and was related to colon cancer risk [42]. Hsa-mir-224-5p was related to the pathogenesis of COAD [11]. Hsa-mir-495-3p might play an important role in the tumorigenesis of glioma [43]. Interestingly, our results showed that hsa-miR-21-5p, hsa-mir-21-3p, and hsa-mir-409-3p were downregulated in CTS, while their expression levels were increased in CRC tissues [44][45][46]. Many of the significantly upregulated DEmiRs in CTS have been associated with CRC. For example, the expression levels of hsa-mir-375 were 5.74 times higher in CTS than in normal stroma (p = 0.03). This miRNA was also differentially expressed between CRC and normal tissue [47]. Hsa-mir-192-3p was enriched in different extracellular vesicle subtypes derived from CRC cell lines [48]. The elevated expression of hsa-mir-215-5pwas associated with CRC metastasis [49]. Hsa-mir-378a-3pwas differentially expressed between human colorectal carcinoma and normal colonic mucosa [50]. The analysis of circulating miRNA expression profiles revealed that hsa-mir-498was upregulated in colon cancer patients [46]. Therefore, a majority of the aberrantly expressed DEmiRs in CTS relative to normal stroma identified by the bioinformatics approach have been related to the pathogenesis of colon cancer.

Identification of Pathways Significantly Associated with DEmiRs
We identified 25 KEGG [39] pathways significantly associated with the downregulated DEmiRs in CTS using DIANA-miRPath [21]. These pathways are mainly involved in cancer (pathways in cancer, pancreatic cancer, and proteoglycans in cancer), cellular development (cell adhesion molecules, Wnt, TGF-β, FoxO, adherens junction, Gap junction, and ErbB pathways), immune regulation (cytokine-cytokine receptor interaction), and metabolism (glycosphingolipid biosynthesis-lacto and neolacto series, biosynthesis of unsaturated fatty acids, and ubiquitin mediated proteolysis) (Supplementary Table S4). Previous studies have shown that some of these pathways are significantly associated with colon cancer or CTS. For example, the Wnt signaling, cell adhesion molecules, adherens junction, gap junction, and cytokine-cytokine receptor interaction pathways have been associated with CTS [6]. Genetic variation in the TGF-β signaling pathway may lead to an increased risk of CRC [57]. ErbB signaling pathway is involved in cell proliferation, invasion, metastasis, and neoangiogenesis in colon cancer [58]. Moreover, we identified 10 KEGG pathways that were significantly associated with the upregulated DEmiRs in CTS (Supplementary Table S5). Many of these pathways have been associated with colon cancer. For example, the extracellular matrix (ECM)-receptor interaction pathway is highly enriched in CTS [6]. Dysregulation of proteoglycans is associated with colon cancer prognosis [59].

Identification of Differentially Expressed Hub Genes Targeted by DEmiRs
We identified 460 DEGs between CTS and normal stroma by a meta-analysis of two GEO microarray datasets (GSE31279 andGSE35602). These DEGs included 160 downregulated and 300 upregulated genes in CTS (Supplementary Tables S6 and S7). Based on the 460 DEGs, we constructed a PPI network to identify hub genes (Supplementary Table S8) and found eight hub genes (degree ≥ 10) targeted by the DEmiRs (Figure 1A). All the eight genes (COL5A2, EDNRA, OLR1, TGFBI, MET, TNFRSF11B, TWIST1, and WNT5A) were upregulated in CTS and were targeted by the downregulated DEmiRs ( Figure 1B). COL5A2 is predominantly expressed in COAD stroma [60]. A meta-analysis of gene expression profiles in 950 cancer cell lines revealed that OLR1 was upregulated in 20% of CRC cells [61]. Interestingly, we found that all eight hub genes (COL5A2, EDNRA, TGFBI, MET, OLR1, TNFRSF11B, TWIST1, and WNT5A) were also upregulated in TCGA COAD samples versus normal samples (Student′s t test, p < 0.001) (Figure 2). It indicates that stromal components may have a contribution to COAD transcriptome. However, none of the hub genes showed significant expression differences between colon tumor stromal fibroblast and normal fibroblast (Supplementary Table S9), suggesting that these hub genes are specifically expressed in other CTS cells. Interestingly, we found that all eight hub genes (COL5A2, EDNRA, TGFBI, MET, OLR1, TNFRSF11B, TWIST1, and WNT5A) were also upregulated in TCGA COAD samples versus normal samples (Student s t test, p < 0.001) (Figure 2). It indicates that stromal components may have a contribution to COAD transcriptome. However, none of the hub genes showed significant expression differences between colon tumor stromal fibroblast and normal fibroblast (Supplementary Table S9), suggesting that these hub genes are specifically expressed in other CTS cells.

The Hub Genes Are Negative Prognostic Factors in CRC
Survival analyses using Prognoscan [24] showed that three hub genes (COL5A2, EDNRA, and OLR1) had a negative expression correlation with survival prognosis (overall survival (OS) and/or disease-free survival (DFS)) in CRC ( Figure 3A). Moreover, the elevated expression of COL5A2, EDNRA, TGFBI, and OLR1 was associated with a worse DFS prognosis in TCGA COAD ( Figure 3B). Previous studies have shown that some of the hub genes are associated with prognosis in colon cancer. For example, the depressed expression of EDNRA is associated with a better prognosis in colon cancer [62]. Taken together, these data demonstrate that the upregulated hub genes regulated by the downregulated DEmiRs in CTS are negative prognostic factors in colon cancer.

The Hub Genes Are Negative Prognostic Factors in CRC
Survival analyses using Prognoscan [24] showed that three hub genes (COL5A2, EDNRA, and OLR1) had a negative expression correlation with survival prognosis (overall survival (OS) and/or disease-free survival (DFS)) in CRC ( Figure 3A). Moreover, the elevated expression of COL5A2, EDNRA, TGFBI, and OLR1 was associated with a worse DFS prognosis in TCGA COAD ( Figure 3B). Previous studies have shown that some of the hub genes are associated with prognosis in colon cancer. For example, the depressed expression of EDNRA is associated with a better prognosis in colon cancer [62]. Taken together, these data demonstrate that the upregulated hub genes regulated by the downregulated DEmiRs in CTS are negative prognostic factors in colon cancer.   . Kaplan-Meier curves shows that the elevated expression of four hub genes (COL5A2, EDNRA, OLR1, and TGFBI) is associated with a worse DFS prognosis in the TCGA colon cancer cohort by GEPIA [23] (log-rank test, p < 0.05).

The Elevated Expression of Hub Genes is Associated with the Immunosuppressive TME in Colon Cancer
The tumor-infiltrating lymphocytes (TILs) level is an independent predictor of survival in CRC [63]. We analyzed the correlations between the expression levels of three hub genes and the levels of TILs in TCGA COAD. The three hub genes included COL5A2, EDNRA, and OLR1 which were identified as negative prognostic factors in both TCGA COAD [33] and two microarray CRC cohorts (GSE17536 [30,31] and GSE14333) [32]. Interestingly, we found that the elevated expression of the three genes was consistently correlated with the upregulation of immunosuppressive signatures in colon cancer, such as tumor-associated macrophages (TAMs), M2 macrophages, regulatory T cells (Tregs), and T cell exhaustion (Table 3). Moreover, the expression levels of the three genes exhibited significant positive correlations with the expression levels of immune-inhibitory marker genes, including PD-L1, PD-L2, TGFB1, and TGFBR1 ( Figure 4A). Interestingly, we found that the ratios of pro-/anti- . Kaplan-Meier curves shows that the elevated expression of four hub genes (COL5A2, EDNRA, OLR1, and TGFBI) is associated with a worse DFS prognosis in the TCGA colon cancer cohort by GEPIA [23] (log-rank test, p < 0.05).

The Elevated Expression of Hub Genes Is Associated with the Immunosuppressive TME in Colon Cancer
The tumor-infiltrating lymphocytes (TILs) level is an independent predictor of survival in CRC [63]. We analyzed the correlations between the expression levels of three hub genes and the levels of TILs in TCGA COAD. The three hub genes included COL5A2, EDNRA, and OLR1 which were identified as negative prognostic factors in both TCGA COAD [33] and two microarray CRC cohorts (GSE17536 [30,31] and GSE14333) [32]. Interestingly, we found that the elevated expression of the three genes was consistently correlated with the upregulation of immunosuppressive signatures in colon cancer, such as tumor-associated macrophages (TAMs), M2 macrophages, regulatory T cells (Tregs), and T cell exhaustion (Table 3). Moreover, the expression levels of the three genes exhibited significant positive correlations with the expression levels of immune-inhibitory marker genes, including PD-L1, PD-L2, TGFB1, and TGFBR1 ( Figure 4A). Interestingly, we found that the ratios of pro-/anti-inflammatory cytokines were significantly lower in the colon cancers highly expressing the hub genes (expression levels > median) than in those lowly expressing the hub genes (expression levels < median) in the TCGA COAD cohort (Student's t test, p < 0.001) ( Figure 4B). Altogether, these results suggest that the elevated expression of these hub genes is associated with the immunosuppressive TME in colon cancer that could explain their oncogenic function to a certain degree. inflammatory cytokines were significantly lower in the colon cancers highly expressing the hub genes (expression levels > median) than in those lowly expressing the hub genes (expression levels < median) in the TCGA COAD cohort (Student's t test, p < 0.001) ( Figure 4B). Altogether, these results suggest that the elevated expression of these hub genes is associated with the immunosuppressive TME in colon cancer that could explain their oncogenic function to a certain degree.  The correlation analyses were performed using TIMER [25]. The Spearman's correlation test p-values (p) and correlation coefficients (cor) are shown. RSEM: RNA-Seq by Expectation Maximization [64]. (B) The ratios of pro-/antiinflammatory cytokines are significantly lower in the colon cancers highly expressing the hub genes (expression levels > median) than in those lowly expressing the hub genes (expression levels < median) in the TCGA COAD cohort (Student′s t test, p < 0.001). The ratio of pro-/anti-inflammatory cytokines in a tumor sample is defined as the ratio of average expression levels (log2-transformed) of their marker genes. The pro-inflammatory cytokines represent the immune-stimulatory signature with marker genes IFNG, IL-1A, IL-1B, and IL-2, and the anti-inflammatory cytokines represent the immune-inhibitory signature with marker genes IL-4, IL-10, IL-11, and TGFB1 [40]. ***, p < 0.001. Table 3. Expression correlations between three hub genes and marker genes of immune-inhibitory signatures evaluated by TIMER [25].  . The correlation analyses were performed using TIMER [25]. The Spearman's correlation test p-values (p) and correlation coefficients (cor) are shown. RSEM: RNA-Seq by Expectation Maximization [64]. (B) The ratios of pro-/anti-inflammatory cytokines are significantly lower in the colon cancers highly expressing the hub genes (expression levels > median) than in those lowly expressing the hub genes (expression levels < median) in the TCGA COAD cohort (Student s t test, p < 0.001). The ratio of pro-/anti-inflammatory cytokines in a tumor sample is defined as the ratio of average expression levels (log2-transformed) of their marker genes. The pro-inflammatory cytokines represent the immune-stimulatory signature with marker genes IFNG, IL-1A, IL-1B, and IL-2, and the anti-inflammatory cytokines represent the immune-inhibitory signature with marker genes IL-4, IL-10, IL-11, and TGFB1 [40]. ***, p < 0.001. Table 3. Expression correlations between three hub genes and marker genes of immune-inhibitory signatures evaluated by TIMER [25].

Identification of Candidate Drugs Targeting Hub Genes
Using DGIdb [27], we identified eight FDA-approved drugs that potentially target the protein products of two hub genes (EDNRA and COL5A2) ( Table 4). Table 4 shows that the drugs targeting EDNRA are mostly antagonists. EDNRA (Endothelin receptor) antagonists ambrisentan, bosentan, and macitentan have been evaluated as hepatobiliary transporter inhibitors and substrates in sandwich-cultured human hepatocytes [65]. The efficacy of macitentan, a dual endothelin receptor antagonist, has been evaluated in treating the brain metastases of breast and lung cancers in mice [66]. The use of macitentan, a dual endothelin receptor antagonist, in combination with temozolomide, may lead to glioblastoma regression and long-term survival in mice [67]. Nevertheless, to our knowledge, inhibition of EDNRA or COL5A2 has not been tested for colon cancer treatment. Our data suggest that these genes could be promising targets for development of anticancer drugs in treating colon cancer patients.

Discussion
For the first time, we built the miRNA-mRNA regulatory network based on the aberrantly expressed miRNAs and genes in CTS. We identified seven downregulated and 17 upregulated DEmiRs in CTS relative to normal stroma by a meta-analysis. A recent study has revealed that TFs can modulate the expression of miRNAs, and they have a synthetic effect on the expression of their common target genes in colon cancer [68]. Thus, we identified upstream TFs that potentially regulated the DEmiRs. Notably, SP1, a zinc finger TF, potentially regulated a majority of the DEmiRs. Pathway enrichment analysis revealed that these DEmiRs were significantly associated with the pathways involved in cancer, cellular development, immune regulation, and metabolism. Furthermore, we identified DEGs between CTS and normal stroma and constructed a PPI network of these DEGs. From the PPI network, we identified eight hub genes which were significantly upregulated in CTS and were targeted by the downregulated DEmiRs. Among the eight hub genes, three (COL5A2, EDNRA, and OLR1) showed significant inverse expression correlations with OS and/or DFS in three CRC cohorts ( Figure 3). Moreover, the elevated expression of the three prognostic hub genes (PHGs) was significantly correlated with immunosuppressive signatures in colon cancer (Figure 4). Network analysis showed that many immune-inhibitory signature genes interacted with the PHGs (Figure 5A). For example, OLR1 interected with CCL2 and CD68. These results suggest that the negative prognostic impact of PHGs in colon cancer may result from the suppressive antitumor immune microenvironment. Collectively, the expression pattern, prognostic effect, immune-inhibitory activity demonstrated the oncogenic effect of COL5A2, EDNRA, and OLR1 in colon cancer.
A major limitation of this study is that the CTS-associated miRNA-mRNA regulatory network identified by in silico analysis has not been proved by experimental validation. Thus, although our findings could provide potentially useful biomarkers for colon cancer diagnosis and prognosis, as well as therapeutic targets, further experimental and clinical validation is needed to translate these findings into clinical application.

Conclusions
The identification of transcriptional markers and miRNA-miRNA regulatory networks in CTS may provide new biomarkers for colon cancer diagnosis, prognosis, and treatment, as well as insights into the mechanism of tumor immune microenvironment regulation in colon cancer.  Table S6. Upregulated genes in CTS versus colon normal stroma. Table S7. Downregulated genes in CTS versus colon normal stroma. Table S8. The hub genes (degree ≥ 10) in the protein-protein interaction network of DEGs between CTS and normal stroma. Table S9. Comparisons of the expression levels of eight hub gene between colon tumor stromal fibroblast and normal fibroblast.
Author Contributions: M.N.U. conceived the research, performed data analyses, and wrote the manuscript. M.L. performed data analyses and helped prepare the manuscript. X.W. conceived the research, designed analysis strategies, and wrote the manuscript. All the authors read and approved the final manuscript.
Funding: This work was supported by China Pharmaceutical University (grant numbers 3150120001 and 2632018YX01 to X.W.).

Conflicts of Interest:
The authors declare no conflict of interest.