Impact of microRNAs on regulatory networks and pathways in human colorectal carcinogenesis and development of metastasis

Matched miRNA and gene expression in normal colon mucosa, primary colon carcinoma and liver metastasis

Of 847 miRNAs represented in the array, 309 were considered for later analyses after filtering (see Methods). Similarly, 30% of genes with low expression profile variability were filtered out, and 15,761 genes with moderately to highly variable expression profiles were considered further.

Unsupervised hierarchical cluster analysis was performed on selected miRNAs and genes. Normal samples clustered together and were relatively well separated from T and M samples, both according to unsupervised hierarchical cluster analysis based on miRNA and on gene expression data (Figure 1A). In fact, only miRNA expression profiles separated T from M quite well. Considerable per-patient pairing of T and M samples was observed in both dendrograms, in which triplets and pairs of samples from the same patient are shown in the same color (see also Additional file 1: Figure S1). In roughly 25% of patients, the M samples were more similar to the T from which they had derived, rather than to the M samples of other patients. The similarity between T and M samples of the same patient was more evident when samples were classified according to gene expression (20% and 28% of per-patient sample pairing, in miRNA- and gene-based heatmaps, respectively).

miRNA expression variability in the two main transitions in the N-T-M progression was evaluated, by considering the number of miRNAs and genes resulting up- or down modulated (absolute fold-change > 1) or unchanged in the T vs N and in the M vs T comparisons. Figure 1B shows that most expression variation took place in the T vs N comparison, revealing two alternative patterns of miRNA (and gene) expression variation in the progression: i) miRNAs (genes) up- or down-modulated in the comparison T vs N and basically stable in M vs T (Figure 1B, left), and ii) miRNAs (genes) unchanged in T vs N and modulated in M vs T (Figure 1B, right).

Differentially expressed miRNAs

Many significantly differentially expressed miRNAs (DEMs) can be found during tumor progression. The whole set of samples with miRNA expression data was considered in unpaired tests (e.g., all T vs all N samples). Other sample subsets matched by patient were considered for paired comparisons (Additional file 1: Methods and Results and Additional file 1: Figure S2).

By considering the larger unpaired dataset, we identified 62, 63 and 11 DEMs in T vs N, M vs N and M vs T comparisons, respectively (Additional file 1: Table S1). Several miRNAs were significantly modulated in more than one contrast. Only 5 miRNAs varied merely when M and T samples were compared (miR-146a, miR-15a, miR-15b, miR-196a, miR-708). Of 53 DEMs shared by at least two comparisons, 25 were always under- and 26 over-expressed, whereas two miRNAs did not follow the same trend in the various comparisons. miR-100 and mir-99a, both putative tumor suppressors, were under- expressed in T vs N and over-expressed in M vs T.

Identification of DEM confirmed that more miRNAs are modulated in T vs N than in M vs T comparison; however, DEMs in metastasis compared with primary tumors may be of great importance, since they include miR-10b, miR-210 and miR-708, which are key regulators of several processes related to disease progression, such as DNA repair, angiogenesis, hypoxia, EMT induction, and cancer recognition by the immune system [16–18]. Additional file 1: Figure S3 shows the expression profiles of 22 miRNAs involved in EMT which were differentially expressed in T vs N and/or M vs T comparisons.

miRNA biogenesis is not impaired during cancer progression

According to our results, during tumor progression we could not detect either general impairment of miRNA biogenesis or an overall increase in miRNA expression. In fact, in both T vs N and/or M vs T contrasts, we observed a comparable number of up- and down-modulated DEMs: 29 and 33 were respectively up- and down-modulated in T vs N; as against 5 and 6 in M vs T. Regarding global miRNA expression, when we considered the distribution of all DEM expression levels measured in N, T and M, again we could not find any significant differences between the groups (mean values 6.10, 5.95 and 6.01, respectively; p-value of pairwise mean equality t-test >0.7).

qRT-PCR validation of miRNA expression

To assess the reproducibility of the identified miRNAs, we used qRT-PCR to measure the expression of 5 miRNAs, both in the same 78 samples used for miRNA expression profiling and in an independent set of 21 samples obtained from matched samples of 7 patients. We quantified three down-regulated (miR-150, miR-10b, miR-146a) and two up-regulated miRNA (miR-210 and miR-122). miR-150 was the most down-regulated miRNA in the T vs N comparison; miR-10b was significantly down-regulated in all three contrasts; miR-122 and miR-146a were differentially expressed in the M vs T comparison; and miR-210 was up-regulated in the M vs T and M vs N comparisons. For each of the above miRNAs, the Spearman rank correlation test between qRT-PCR expression estimation (${2{-}^{∆\mathrm{C}}}_{\mathrm{t}}$) and array-based expression level, in the same 78 samples, was significant (p < 0.01) (Additional file 1: Figure S4A). miR-150 and miR-146a were also tested on an independent set of 21 samples (N, T and M samples from 7 new patients) and the results were fully consistent with microarray-based observations, confirming miR-150 down-regulation in T vs N and miR-146a down-regulation in M vs T (Additional file 1: Figure S4B).

The expression level of miR-122, a miRNA which is highly expressed in the liver, was significantly higher in M than in T, suggesting the residual presence of normal liver tissue in metastasis samples. miR-122 expression was evaluated in samples in which normal colon mucosa and liver tissue had been mixed in different proportions. qRT- PCR analysis showed that mixing up to 95% of N with 5% of liver tissue caused a significant increase in miR-122 expression, compared with samples with 100% of N tissue. This result indicates that the presence of only 5% of liver cells in a sample significantly increases miR-122 expression.

Identification of most probable miRNA targets and definition of regulatory networks modulated in development of tumor and metastasis

We obtained a set of 77 samples (23 N, 29 T, 25 M) with matched miRNAs and gene expression data from the same biological sample. The combined analysis of target prediction and expression of 305 miRNAs with predicted targets and 12,748 target genes allowed us to reconstruct post-transcriptional regulatory networks describing the most probable regulatory interactions and circuits active in transitions characterizing tumor development and progression.

miRNAs interact in several ways with target mRNAs, and may exert non-canonical regulatory actions, but they commonly act post-transcriptionally by altering the stability of target mRNAs. The expression profile of a given miRNA was therefore expected to be inversely correlated with that of its target genes. To identify the most probable miRNA targets, we first adopted a classical approach, to enrich a large set of predicted miRNA targets in truly regulated genes, by identifying significant negative correlations between miRNAs and predicted target expression profiles [19, 20]. With the selected threshold on correlation significance (see Methods), we identified 3,078 relations, involving 117 miRNAs and 1,423 target genes (correlation value < −0.412). It should be noted that less than 1 out of 1,500 miRNA/target relations predicted by miRSVR turned out to be supported by expression data. Indeed, among supported relationships, 2,690 (87%) were based on predicted target sites, which are conserved across species, whereas only one-quarter of the original predictions involved conserved sites (see Methods).

The number of supported target genes per miRNA ranged from 1 to 216, with an average value of 26.3. About half the genes were supported targets of only one specific miRNA, whereas other genes were putatively regulated by up to 10 different miRNAs.

Within the whole group of supported target genes, of particular importance are the various subsets of genes significantly differentially expressed (DEGs), which may be viewed as the prominent effects of post-transcriptional regulatory action exerted by miRNAs.

As previously stated, in the unpaired analysis only 5 miRNAs appeared to be significantly modulated in the transition from primary tumor to metastasis, whereas the number of DEMs observed in the comparison between normal tissue and primary tumor and metastasis was higher (Figure 1C). It is reasonable to assume that DEMs, being moderately to highly expressed in all the considered samples and those most variable in each contrast, are responsible for most target repression. Thus, for the T vs N and M vs N comparisons, we focused only on DEMs up- and down-regulated in each contrast with a fold change (FC) > 3 and with an average expression over background in at least one of the two contrasted groups. The intersection between the post- transcriptional network and the results of miRNA differential expression analysis induced various sub-networks, describing post-transcriptional regulatory circuits involving those miRNAs whose expression variation may be important for tumor development and progression. Figures 2A and 3A show post-transcriptional regulatory networks with miRNAs differentially expressed in the T vs N comparison and their supported relations with target genes. Two network components are observed, involving respectively 6 up-regulated (Figure 2A) and 17 down-regulated (Figure 3A) DEMs. The component regarding the 6 up-regulated miRNAs was smaller, and a large fraction of genes appeared to be regulated by miR-182. The largest component involved 17 down-regulated DEMs, which putatively regulated a number of targets about twice as high as that observed for the whole set of miRNAs. A large proportion of the predicted target genes were indeed significantly differentially expressed. This was more evident for the up-regulated miRNA component, in which 62 down- regulated DEGs accounted for 27% of the supported targets of up-regulated DEMs. Conversely, our analysis predicted that only a few genes would be targets of more than one up-regulated miRNA. Among the interactions predicted between miRNA and target genes identified in the T vs N regulatory network we considered the relationship between miR-145 and c-Myc. qRT-PCR analysis of the entire set of samples confirmed the opposite behavior of miR-145 and c-Myc in tumour progression (Figure 3B).

The network involving five DEMs observed in the M vs T comparison and their supported target genes was small and consisted of five unconnected components (Additional file 1: Figure S5). No DEGs were observed among the supported targets of DEMs in the T to M transition. These results further support the concept that the main change in the transcriptome occurs early in the CRC progression.

miR-10b expression in metastases: potential association with survival

The 26 miRNAs included in the TN and MT networks (Additional file 1: Table S3) were considered for survival analysis. The expression of miR-10b measured in liver metastasis showed a statistically significant association with the survival of patients affected with stage IV CRC (hazard ratio = 1.47, 95% confidence interval = 1.23-1.75; adjusted p-value: 0.00052). The effect of miR-10b on prognosis is given in Additional file 1: Figure S6, which shows that patients with high levels of miR-10b expression in their metastatic disease have a shorter time to event (median survival: 8 months) compared with those with low levels (51 months). In our study, the expression levels of miR-10b measured in primary tumors had no significant impact on prognosis.

miR-182 controls ENTPD5

miR-182 was one of the most upregulated DEMs in the T vs N contrast. Among its predicted target genes we focused on ENTPD5, due to the involvement of the gene product in energy metabolism. Since ENTPD5 resulted significantly down-regulated in our analysis, and another study provided support for downregulation during cancer progression [21], we decided to study the relationship between miR-182 and its predicted target gene ENTPD5. To investigate the opposite behavior of the miRNA and its target gene, we performed qRT-PCR to measure miR-182 and ENTPD5 expression in a panel of five cell lines (CG-705, HT29, from a primary colorectal tumor, and MICOL-S, MICOL-14, and LoVo, from colorectal carcinoma metastases [22]). As shown in Figure 2B, all tumor cell lines showed an inverse correlation (−0.85, p-value < 0.05) between high expression of miR-182 and low expression of ENTPD5. This result not only confirms the microarray profiling data but also suggests a role of the anti-correlated relationship in conferring some advantageous properties to the tumor cells. We next wanted to provide support to the direct targeting of ENTPD5 by miR-182. To this end, we performed a Luciferase reporter assay in HEK293T cells transfected with a construct containing the firefly luciferase gene fused to the 3′-UTR of ENTPD5 (pMIR-ENTPD5). When cells were co-transfected with miR-182, a 50% reduction in luciferase expression was observed, compared with cells transfected with negative control RNA or pMIR-ENTPD5 only (Figure 2C).

miRNA modulated KEGG pathways

For the T vs N contrast, we identified pathways enriched in gene targets of DEMs, the expression of which changes in the comparisons, and in DEG targets of DEMs, with only partially overlapping results. Some pathways, such as “Cell cycle”, “Purine metabolism” and “Pathways in cancer”, were found in both cases, whereas others, such as “P53 signaling” was identified by the first, less conservative strategy. Conversely, “Wnt signaling” and “Colorectal cancer” pathways were found to be specifically enriched in DEG targets of DEMs (Additional file 1: Table S4). It should be noted that “Pathways in cancer” is a collection of many pathways, including various cancer hallmarks such as EMT, and is an indication that many important processes may be under the control of identified miRNAs.

Mixed regulatory circuits involving interplay of miRNAs and transcription factors contain many “cancer genes”

miRNAs regulate target genes mainly at the post-transcriptional level, operating in highly interconnected regulatory networks and pathways, with complex cross-talk of miRNAs and transcription factors (TFs), which are frequently master regulators of biological processes. miRNA expression can be activated or repressed by transcription factors (TFs), whereas mRNAs encoding TFs can be silenced by miRNAs. In many cases, key mixed regulatory circuits involve miRNAs, TFs and common target genes. Thus, miRNAs and TFs can form feedback or feedforward loops, cooperating to switch or tune gene expression.

To account for the interplay of miRNAs and TFs, we have recently developed a new method for the integrated analysis of target prediction, MAGIA2 [23], which allows one to dissect regulatory complexity by exploiting target predictions, in combination with information on experimentally validated TF-miRNA and TF-gene interactions, and which can identify both negative and positive expression profile correlations. By applying this method to the matched expression profiles of all genes and to the 70 DEMs in TN and/or MT comparisons, we were able to reconstruct mixed regulatory networks involving miRNAs and TFs. Our results show that 18 miRNAs are involved in the strongest 127 interactions identified as significant by MAGIA2 (Additional file 1: Figure S7A). Then, we focused on the identification of two types of triangular mixed circuits representing putative feed-forward or feed-back loops with interplay between transcriptional and post-transcriptional regulation: (i) circuits in which one TF regulates both a given miRNA and its target gene, and (ii) circuits in which one miRNA regulates both a given TF and its regulated gene. Additional file 1: Figure S7B reports 20 significant mixed circuits involving five distinct miRNAs which interact with 6 TFs and 15 non-TF protein-coding genes.

Patients and collection of tissue samples

For this study, 46 patients with sporadic colorectal adenocarcinomas, who underwent surgery at the University of Padova (Surgery Unit, Department of Surgery, Oncology and Gastroenterology) between March 1994 and September 2008, were selected from the institutional CRC database. Patients with a known history of a hereditary colorectal cancer syndrome were excluded. The Ethics Committee of the University Hospital of Padova approved the study. All patients provided written informed consent. Enrolled patients did not receive any neo-adjuvant treatment. Normal mucosa samples were taken at a minimum distance of 10 centimeters from the tumor site. All samples were immediately snap-frozen in liquid nitrogen and stored at −80° until use.

RNA extraction

7 μm sections from each tissue sample were prepared using a Leica CM 1950 cryostat (Leica Microsystems, Wetzlar, Germany). Hematoxylin and eosin stained sections of each specimen were prepared and re-evaluated by one experienced pathologist (G.E.); only samples with more than 80% of vital tumor tissue were considered for RNA extraction in toto. Laser microdissection was performed on a few frozen samples of primary tumours and metastases with a proportion of neoplastic cells lower than 80% using LMD-6000 Laser Microdissection System (Leica Microsystems, Wetzlar, Germany). Total RNA from samples was isolated using Trizol (Life Technology Corp, Carlsbad, CA, USA) according to the manufacturer’s instructions. RNA concentration was quantified on a NanoDrop 1000 Spectrophotometer (NanoDrop Technologies, Waltham, MA, USA). RNA quality was evaluated by RNA 6000 Nano LabChip (Agilent Technologies, Santa Clara, CA, USA) on an Agilent 2100 Bioanalyzer. Samples with RNA integrity number (RIN) < 6 were excluded.

Expression profiling

miRNA microarray hybridization was performed from total RNA with the Affymetrix GeneChip miRNA Array 2.0 (Affymetrix Santa Clara, CA, USA). Briefly, 100 ng of total RNA from each sample were labeled with the FlashTag Biotin RNA Labeling Kit (Genisphere, Hatfield, PA, USA) and samples were then hybridized according to the manufacturer’s instructions.

GeneChip Human Exon 1.0 ST (Affymetrix) microarray hybridization was performed from the same total RNA extracted for miRNA profiling. Total RNA (100 ng from each sample) was labeled with the Ambion WT expression kit (Ambion Inc, Austin, TX, USA), as provided by the manufacturer. End-labeling, hybridization, washing and scanning were performed according to the GeneChip Whole Transcript (WT) Sense Target Labeling Assay user manual (Affymetrix), and scanned with an Affymetrix GCS 3000 7G scanner.

For both miRNAs and exon arrays, a first quality control check was performed with Affymetrix® Expression Console™ software (v.1.0) to determine the success of hybridizations. miRNA and gene expression measures were reconstructed from .cel files by using the Robust Multichip Average (RMA) method. Signals of 41 probes per gene, on average, were summarized to estimate gene expression with EntrezGene-based custom CDF (http://brainarray.mbni.med.umich.edu/Brainarray), obtaining expression profiles of 22,517 genes.

Detailed quality control of samples was carried out with R software: Normalized Unscaled Standard Error (NUSE) and Relative Log Expression (RLE) for global quality of signals in each array assessment, and MA plots before and after RMA for identification of biases associated with specific intensity classes.

miRNAs detected in fewer than 20 samples were discarded, without filtering out miRNAs undetected only in one sample class.

miRNA and gene differential expression

A Significance Analysis of Microarrays (SAM) with an unpaired two-class design was performed to identify differences in miRNA expression between groups of normal mucosa (N), primary tumor (T) and liver metastases (M) samples.

Differentially expressed genes (DEG) between groups of N, T and M samples were identified with SAM with an unpaired two-class unpaired design. The cut-off for significance (determined by tuning parameter delta) corresponded to a false discovery rate (FDR) < 0.01.

Integrated analysis of miRNAs and gene expression profiles

We integrated target predictions with correlation-based miRNA and gene expression profile, to identify those regulatory relationships significantly supported by expression data. This analysis is based on the assumption that, at least for miRNAs acting on mRNA stability, the expression profile of a given miRNA is expected to be inversely correlated with those of its true targets. First, miRNA target predictions were obtained with miRanda-miRSVR, a method able to identify miRNA-mRNA interactions involving the seed sequence as well as non-conserved and non-canonical sites [50].

Using a genome-wide approach, we considered over four million miRNA-gene relationships predicted with a “good” score, 25% of which involved evolutionarily conserved target sites. Pair-wise Pearson correlations between miRNA and predicted target gene expression profiles were then calculated and assigned a statistical significance (p-value). The FDR (q-value) was used to correct the correlation significance for multiple testing. Only miRNA-target relationships with significant correlations (FDR < 0.01) were considered to be supported by expression profiles.

Considering that expression values ranged from 0 to 15 in log2 scale and that 70% of miRNAs were weakly expressed, we then focused only on miRNAs with signal over 5, arbitrarily set to the average value of miRNA expression.

Enrichment of KEGG pathways

To identify significantly perturbed KEGG pathways for each considered comparison, we applied the statistical procedure GAGE [51], which takes into account gene expression variations in both directions, to genes turning out to be supported targets of DEMs (p-value < 0.05). In particular, two strategies were implemented to identify pathways enriched in genes targets of DEMs, the expression of which changes in the considered comparisons, and pathways enriched in DEGs targets of DEMs.

Survival analysis

We investigated the association between the expression levels of 26 DEMs (present in reconstructed post-transcriptional regulatory networks) in biopsies obtained from distinct primary [n = 26] or metastatic [n = 20] colorectal cancers and patients’ disease- specific survival (interval between diagnosis of primary or metastatic disease and death by disease or last follow-up).

Given the relatively low sample size, only univariate survival analysis was performed, and the Cox proportional hazard regression model was used, assuming a linear functional form of the covariates being assumed. The risk associated with a unit increase in miRNA levels was expressed as hazard ratio (HR) and its 95% confidence interval (CI). With Bonferroni’s p-value adjustment for multiple comparisons, the alpha level of significance was set at 0.002. In order to illustrate prognosis associated with different levels of the relevant miRNAs, Kaplan-Meier survival curves were generated after dichotomizing (high vs low categories) originally continuous covariates based on the median values of miRNA expression levels. All analyses were performed with Stata/SE software (version 11.0, StataCorp LP, College Station, TX, USA).

Quantitative RT-PCR

To confirm array data for miR-150, miR-146a, miR-10b, miR-122 and miR-210 and to validate miR-145/c-Myc and miR-182/ENTPD5 relationships, we conducted qRT- PCR experiments, as previously described [52]. Briefly, experiments were performed three times in triplicate with a LightCycler 480 II (Roche Diagnostics, Basel, Switzerland) and were analysed by the ∆∆Ct method. Quantification of the five selected miRNAs, miR-145 and c-Myc was normalized against internal housekeeping controls, selected for their minimum variability (measured as expression profile Shannon entropy). RNU44 and GAPDH were used as internal controls, respectively for miR-182 and ENTPD5 gene expression quantification.

Luciferase reporter assay

HEK293T cells were plated in 24-well plates at 8×104 cells/well and cotransfected using Lipofectamine 2000 (Invitrogen), with 500 ng pMir target vector customized for ENTPD5 3′UTR (OriGene Technologies, Rockville, MD, USA) and 250 ng pRL-TK (Promega, Milan, Italy) following the manufacturers’ instructions. For miR-182 analysis, cells were cotransfected with non-target RNA (Tema ricerca, Bologna, Italy) as negative control or miCENTURY OX miNatural for hsa miR-182 (Tema ricerca), in triplicate. For the above analyses, cell lysates were analysed 30 hours after transfection by the Dual-Glo Luciferase Assay System (Promega) and experiments were independently repeated three times.

Reconstruction of mixed miRNA-TF networks and circuits

Matched miRNA and gene expression data were analysed with MAGIA2 to construct mixed miRNA-TF networks and circuits. MAGIA2 was run on the expression data for the top 75% of genes with greatest differences in expression and of 70 DEMs in TN and/or MT comparisons. The TargetScan (default settings) method was selected for miRNA target prediction. For TFs, MAGIA2 utilizes experimentally validated TF– miRNA interactions from mirGen2.0 and TransmiR and on TF–gene interactions from the ‘TFBS conserved’ track of the UCSC human genome annotation. Pearson’s correlation was used as a profile association measure.