Genome-wide mRNA and miRNA expression profiling reveal multiple regulatory networks in colorectal cancer

Despite recent advances in cancer management, colorectal cancer (CRC) remains the third most common cancer and a major health-care problem worldwide. MicroRNAs have recently emerged as key regulators of cancer development and progression by targeting multiple cancer-related genes; however, such regulatory networks are not well characterized in CRC. Thus, the aim of this study was to perform global messenger RNA (mRNA) and microRNA expression profiling in the same CRC samples and adjacent normal tissues and to identify potential miRNA-mRNA regulatory networks. Our data revealed 1273 significantly upregulated and 1902 downregulated genes in CRC. Pathway analysis revealed significant enrichment in cell cycle, integrated cancer, Wnt (wingless-type MMTV integration site family member), matrix metalloproteinase, and TGF-β pathways in CRC. Pharmacological inhibition of Wnt (using XAV939 or IWP-2) or TGF-β (using SB-431542) pathways led to dose- and time-dependent inhibition of CRC cell growth. Similarly, our data revealed up- (42) and downregulated (61) microRNAs in the same matched samples. Using target prediction and bioinformatics, ~77% of the upregulated genes were predicted to be targeted by microRNAs found to be downregulated in CRC. We subsequently focused on EZH2 (enhancer of zeste homolog 2 ), which was found to be regulated by hsa-miR-26a-5p and several members of the let-7 (lethal-7) family in CRC. Significant inverse correlation between EZH2 and hsa-miR-26a-5p (R2=0.56, P=0.0001) and hsa-let-7b-5p (R2=0.19, P=0.02) expression was observed in the same samples, corroborating the belief of EZH2 being a bona fide target for these two miRNAs in CRC. Pharmacological inhibition of EZH2 led to significant reduction in trimethylated histone H3 on lysine 27 (H3K27) methylation, marked reduction in cell proliferation, and migration in vitro. Concordantly, small interfering RNA-mediated knockdown of EZH2 led to similar effects on CRC cell growth in vitro. Therefore, our data have revealed several hundred potential miRNA-mRNA regulatory networks in CRC and suggest targeting relevant networks as potential therapeutic strategy for CRC.

and aging. 8,11,12,13,14,15 A variety of methods have been used ranging from miRNA microarrays to global miRNA expression profiling with deep sequencing to determine the expression pattern of miRNAs in cancer tissues. 16,17 Abnormal expression of miRNAs has been associated with a large number of human cancers. In CRC, a number of studies have reported altered miRNA expression pattern, suggesting a plausible role for the aberrant miRNA expression in CRC biology. 17,18 The natural mechanisms for the dysregulation of miRNAs are still largely unknown, although gene amplification, genomic loss, and promoter hypermethylation have been reported as the potential mechanisms in various cancers. 14,19,20 Strategies based on restoration of downregulated miRNAs or inhibition of upregulated miRNAs have opened a new area of investigating the potential therapeutic value in cancer therapy.
Previous studies have examined global messenger RNA (mRNA) and miRNA expression in CRC; 21,22,23 however, only a few have examined the global miRNA and mRNA expression profiling in the same CRC tissue and compared with matched normal tissues. 24 Thus, we performed global mRNA and miRNA expression profiling in 13 CRC specimens and their matched adjacent normal tissues obtained from Saudi patients. We identified more than 700 potential miRNA-mRNA regulatory networks in CRC controlling various key pathways relevant to CRC development, progression, and therapy failure.

Results
Gene expression profiling in CRC. Global gene expression profiling was conducted on 13 colon cancer specimens and 13 adjacent normal tissues. The clinical characteristics of patients involved in current study are provided in Table 1. As shown in Figure 1a, hierarchical clustering based on differentially expressed mRNAs revealed clear separation of the two groups where the cancer tissues clustered separately from the normal tissues, except for one branch of the cancer group representing one cancer sample, which was misclassified to the normal group. Using significance analysis, 1273 up-and 1902 downregulated genes were identified (2.0 FC, Po0.02; Supplementary Table 1). Pathway analysis on the upregulated genes using GeneSpring GX revealed significant enrichment in several pathways related to cell cycle, DNA damage, matrix metalloproteases, Wnt, and TGF-β signaling (Figure 1b and Supplementary Table 2). The Wnt and TGF-β signaling pathways are illustrated in Supplementary Figures 1  and 2. To confirm their relevance to CRC, we used smallmolecule inhibitors to inhibit the Wnt (XAV939 and IWP-2) or TGF-β (SB-431542) signaling in HT115 cells, which led to significant reduction in cell viability (Figure 1c). Selected number of the upregulated (wingless-type MMTV integration site family member 2 (WNT2), matrix metallopeptidase 9 (MMP9), enhancer of zeste homolog 2 (EZH2)) and downregulated (bone morphogenetic protein 3 (BMP3)) genes from the microarray data were subsequently validated using quantitative reverse transcription-PCR (qRT-PCR) (Figure 1d). Small-interfering RNA (siRNA)-mediated knockdown of FOXM1 (forkhead box protein M1) and FOXQ1 (forkhead box protein Q1) reduced HT115 cell growth in vitro (Supplementary Figure 3), corroborating a biological relevance of the identified genes from the microarray in CRC biology. miRNA expression profiling in CRC. To identify potential miRNA-mRNA regulatory networks in CRC, we performed global miRNA expression profiling on the same 13 cancer and adjacent normal samples that were used for mRNA profiling shown in Figure 1. Using significant analysis, we identified 61 significantly downregulated and 42 significantly upregulated miRNAs (1.5-fold change, Po0.02; Table 2). Hierarchical clustering of the differentially expressed miRNAs in 13 colon cancer specimes and 13 normal tissues is shown in Figure 2a. The data revealed clear separation of the two groups. We subsequently focused on the downregulated miRNAs and their correlation with the upregulated target genes. Using TargetScan prediction feature in GeneSpring GX software (Agilent Technologies, Santa Carla, CA, USA), 16 157 genes were predicted to be targeted by the identified downregulated miRNAs (Supplementary Table 3). Pathway analysis using the predicted gene targets revealed significant     Figure 2b. We subsequently focused on the upregulated genes in CRC specimens, which could potentially be regulated by miRNAs found to be downregulated in the same specimens. Crossing the list of predicted gene targets for downregulated miRNAs with the list of the upregulated genes in CRC revealed 794 upregulated genes, which were predicted to be targeted by downregulated miRNAs in CRC (Figure 2c and Supplementary Table 4). The expression levels of selected number of the identified miRNAs from the microarray data (hsa-miR-145-5p, hsa-miR-26a-5p, and hsa-miR-30a-5p) were subsequently validated using Taqman qRT-PCR (Figure 2d).
Depletion of EZH2/PRC2 complex reduces colon cancer cell proliferation and cell migration. Among the identified upregulated genes in our study is EZH2. Elevated expression of EZH2 has been observed in different human cancers, 13,25,26,27 which we found to be upregulated in CRC in the current study as well (Figure 1d and Supplementary  Table 1). To assess the biological ramifications of EZH2 depletion on CRC cancer cells, we treated HT115, HT-29, and SW620 colon cancer cells with 3-deazaneplanocin A (DZNep), a small-molecule inhibitor known to target EZH2 protein, and assessed cell viability on days 4 and 8 posttreatment. As shown in Figure 3a, a significant dose-and time-dependent decrease in colon cancer cell viability was observed, which was associated with a reduction in EZH2 protein expression (Figure 3b, left) and a marked reduction in tri-methylated lysine 27 (H3K27-3me) (Figure 3b, right). Concordant with these data, HT115 cells treated with DZNep also exhibited marked reduction in cell migration as measured using transwell migration assay (Figure 3c). To identify the molecular pathways regulated by EZH2 in CRC, HT115 cells were treated with DZNep to induce reduction in EZH2 and subsequently examined the effects of reduced EZH2 on global gene expression using microarray analysis (Figure 3d). Pathway analysis on the differentially expressed genes revealed multiple enriched pathways including senescence and autophagy, apoptosis, and FAK (Figure 3e and Supplementary Table 5). The senescence and autophagy pathway is illustrated in Supplementary Figure 4, with all matched entities indicated.
EZH2 is regulated by several miRNAs in CRC. Our results suggest that EZH2 is involved in multiple aspects of CRC cell biology; therefore, we hypothesized that the elevated expression of EZH2 in CRC could be attributed to the downregulation of miRNAs that targets EZH2. Figure 4a illustrates the map for EZH2 3′-untranslated region (UTR) and the list of miRNAs predicted to target EZH2 based on TargetScan prediction. Among the predicted miRNAs, EZH2 was found to be regulated by six microRNAs (hsa-miR-26a-5p, hsa-Let-7b-5p, hsa-Let-7c-5p,  hsa-Let-7e-5p, hsa-Let-7g-5p, and hsa-miR-363-3p), which were downregulated in CRC (Figure 4b and Table 2). We subsequently focused on hsa-miR-26a-5p and hsa-let-7b-5p, as we previously reported those two miRNAs to target EZH2 in nasopharyngeal carcinoma. 13 The alignment between EZH2 3′-UTR and these two miRNAs is shown in Figure 4c (upper panel). Overexpression of hsa-miR-26a-5p or hsa-let-7b-5p in HT115 cells led to significant reduction in EZH2 protein levels (Figure 4c, lower panel). Interestingly, significant inverse relationship between EZH2 and hsa-miR-26a-5p (R 2 = 0.56, P = 0.0001) and hsa-let-7b-5p (R 2 = 0.19, P = 0.02) expression by microarray was observed in the 13 CRC and their matched adjacent normal tissue specimens (Figure 4d), corroborating EZH2 being relevant biological target for these two miRNAs in CRC. Exogenous expression of hsa-miR-26a-5p and hsa-let-7b-5p (Figure 4e, upper panel) led to significant reduction in cell viability, similar to those seen with EZH2 knockdown (Figure 4e, lower panel).

Discussion
Although a number of previous studies have examined mRNA and/or miRNA expression in CRC, 21,22,23,28,29 only few studies have examined global mRNA and miRNA expression in the same clinical samples, 24,30 none so far has been conducted in this geographical region. Therefore, the strength of our approach is that it enables the identification of deregulated mRNA-miRNA networks in the same biological specimens. Our data revealed more than 700 potential miRNA-mRNA regulatory networks in CRC and thus provide circumstantial evidence for the involvement of miRNAs in the pathogenesis of colorectal cancer. In addition, our data provide a comprehensive molecular profiling of CRC in Saudi Arabia and the Middle East.
Several of the deregulated mRNAs and miRNAs identified in the current study have been reported previously, suggesting a common underlying molecular mechanisms leading to CRC pathogenesis regardless of ethnicity. For instance, hsa-miR--135b, hsa-miR-223, hsa-miR-18a, hsa-miR-17, hsa-miR-31, and hsa-miR-21 were upregulated in our study, and were also reported by a previous study. 18 Similarly, we found hsa-miR-375, hsa-miR-195, hsa-miR-378, hsa-miR-143, hsa-miR--145, hsa-miR-29c, hsa-miR-1, hsa-miR-30c, hsa-miR-30e, hsa-miR-26a, hsa-miR-100, and hsa-miR-338-3p to be downregulated in CRC in our data, which were also reported to be downregulated in colon cancer patients from Northern Europe. 18 Our data revealed several additional novel miRNAs, which have not been reported previously (Table 2), possibly because of a more comprehensive coverage of miRNAs in the miRNA microarray chips that cover 1205 human miRNAs and used in our study. Our gene expression data revealed multiple deregulated pathways in colon cancer such as cell cycle, DNA damage response, Wnt signaling, and matrix metalloproteases signaling, which is concordant with previous studies implicating these pathways in CRC. 31,32 We found that pharmacological inhibition of Wnt or TGF-β signaling impaired colon cancer cell proliferation in vitro (Figure 1c), which suggests a biological relevance for these pathways in CRC. Several of the identified pathways in CRC were found to be among the predicted targets for miRNAs identified in our study, which suggest a plausible role for the identified miRNAs in the pathogenesis of CRC. We have chosen EZH2, which is a member of the polycomb gene (PcG) family, as it has been implicated in the pathogenesis of a number of other cancer types. 13,[25][26][27] EZH2 is the catalytic subunit of the polycomb repressive complex 2 (PRC2), which is responsible for methylation of lysine 27 on histone H3. This epigenetic modification of H3 is necessary for gene repression through the PRC2 complex. Our current study suggests that EZH2 has a role in the pathogenesis of CRC. Pharmacological inhibition of EZH2 led to significant decrease in H3K27-3me, significant decrease in cell viability, and migration in CRC cells. In addition, siRNA-mediated knockdown of EZH2 exhibited profound effects on colorectal cancer cell growth in vitro.
In silico prediction has identified several potential miRNAs targeting EZH2 in colon cancer cells, and forced expression of hsa-miR-26a-5p and hsa-let7b-5p phenocopied the effects of EZH2 depletion in CRC cells, supporting a role of the two miRNAs in regulating EZH2 expression in colorectal cancer. Our data are concordant with our previous publication implication hsa-miR-26a-5p and hsa-let-7 family in regulating EZH2 in nasopharyngeal carcinoma. 13 Interestingly, we observed significant inverse relationship between EZH2 and hsa-miR-26a-5p and hsa-let-7b-5p expression in CRC (Figure 4d), corroborating the biological relevance of this regulatory network in this disease.
In our current study, we have validated one regulatory network for its relevance in CRC cell biology. However, we provided information regarding several other potential regulatory networks ( Supplementary Tables 3 and 4) in CRC that remain to be investigated.

Materials and Methods
Ethics statement. The clinical study and collection of tissue samples were approved by Institutional Research Ethics Board at the King Saud University College of Medicine (Riyadh, Riyadh, Saudi Arabia).
Patient and tissue collection. Tissue specimens from 13 fresh-frozen consecutive sporadic CRCs matched with their adjacent normal mucosa were obtained from previously untreated patients who underwent surgical resection at the King Khaled University Hospital (Riyadh, Saudi Arabia). Tumor and their paired normal mucosa were selected by an experienced pathologist and specimens were snap frozen in liquid nitrogen and stored at − 80°C until use. Clinical information of the patients is provided in Table 1.
Tissue preparation and RNA isolation. Tissues were ground to powder using a mortar and pestle in the presence of liquid nitrogen. RNA was isolated from 100 to 300 mg of tissue per sample using the Total Tissue RNA Purification Kit from Norgen-Biotek Corp. (Thorold, ON, Canada). The resulting RNA was quantified using NanoDrop 2000 (Thermo Scientific, Wilmington, DE, USA) and the RNA quality and integrity was confirmed using gel electrophoresis.
Gene expression profiling. Total RNA was extracted as described above using Total RNA Purification Kit (Norgen-Biotek Corp.) according to the manufacturer's instructions. One hundred and fifty nanograms of total RNA was labeled and then hybridized to the Agilent Human SurePrint G3 Human GE 8 × 60 k v16 microarray chip (Agilent Technologies). All microarray experiments were conducted at the Microarray Core Facility (Stem Cell Unit, King Saud University College of Medicine). Normalization and data analyses were conducted using GeneSpring GX software (Agilent Technologies). Pathway analysis were conducted using the Single Experiment Pathway analysis feature in GeneSpring 12.0 (Agilent Technologies) as described before. 33,34 Twofold cutoff with Po0.02 was used. miRNA expression profiling. miRNA expression profiling was conducted on the same 13 RNA samples used for gene expression profiling. Two hundred nanograms of the extracted total RNA was used for RNA labeling and hybridization on to the Agilent Human SurePrint G3 8 × 60k v16 miRNA microarray chip according to the manufacturer's protocol. Data were subsequently normalized and analyzed using GeneSpring GX software (Agilent Technologies). A fold-change of 1.5 with Po0.02 was used as cutoff to determine the differentially expressed miRNA in cancer versus normal tissues. Target prediction was conducted using a built-in feature in GeneSpring GX based on TargetScan database.