Colorectal tumor molecular phenotype and miRNA: expression profiles and prognosis.

MiRNAs regulate gene expression by post-transcriptionally suppressing mRNA translation or by causing mRNA degradation. It has been proposed that unique miRNAs influence specific tumor molecular phenotype. In this paper, we test the hypotheses that miRNA expression differs by tumor molecular phenotype and that those differences may influence prognosis. Data come from population-based studies of colorectal cancer conducted in Utah and the Northern California Kaiser Permanente Medical Care Program. A total of 1893 carcinoma samples were run on the Agilent Human miRNA Microarray V19.0 containing 2006 miRNAs. We assessed differences in miRNA expression between TP53-mutated and non-mutated, KRAS-mutated and non-mutated, BRAF-mutated and non-mutated, CpG island methylator phenotype (CIMP) high and CIMP low, and microsatellite instability (MSI) and microsatellite stable (MSS) colon and rectal tumors. Using a Cox proportional hazard model we evaluated if those miRNAs differentially expressed by tumor phenotype influenced survival after adjusting for age, sex, and AJCC stage. There were 22 differentially expressed miRNAs for TP53-mutated colon tumors and 5 for TP53-mutated rectal tumors with a fold change of >1.49 (or <0.67). Additionally, 13 miRNAS were differentially expressed for KRAS-mutated rectal tumors, 8 differentially expressed miRNAs for colon CIMP high tumors, and 2 differentially expressed miRNAs for BRAF-mutated colon tumors. The majority of differentially expressed miRNAS were observed between MSI and MSS tumors (94 differentially expressed miRNAs for colon; 41 differentially expressed miRNAs for rectal tumors). Of these miRNAs differentially expressed between MSI and MSS tumors, the majority were downregulated. Ten of the differentially expressed miRNAs were associated with survival; after adjustment for MSI status, five miRNAS, miR-196b-5p, miR-31-5p, miR-99b-5p, miR-636, and miR-192-3p, were significantly associated with survival. In summary, it appears that the majority of miRNAs that are differentially expressed by tumor molecular phenotype are MSI tumors. However, these miRNAs appear to have minimal effect on prognosis.

Introduction

MiRNAs are small, non-protein-coding RNA molecules that regulate gene expression either by post-transcriptionally suppressing mRNA translation or by causing mRNA degradation[1-6]. We know that miRNAs play a critical role in regulation of proliferation, differentiation, apoptosis, and stress response and are involved in the majority of physiological processes[7, 8]. In our previous work we have shown that colorectal carcinoma cells show widespread dysregulation from normal mucosa[9] and that differences by tumor molecular phenotype also exist[10]. Others also have explored miRNAs with tumor molecular phenotype to gain insight into unique pathways involved in colorectal carcinoma. Most work in this area has focused on microsatellite instability and CpG island methylator phenotype tumors[11], although data also suggest TP53-mutated tumors, KRAS-mutated tumors, and BRAF-mutated tumors may be uniquely associated with specific miRNAs[12-14]. Most studies to date have focused on specific miRNAs to see if they are associated with specific tumor molecular phenotype. Few have looked broadly across the range of miRNAs to identify associations with specific tumor molecular phenotype.

Survival differences by tumor phenotype also have been shown especially for microsatellite instability, CpG island methylator phenotype high, and BRAF-mutated tumors[15-18]. We have observed differences in survival for microsatellite unstable colon and rectal cancer when compared to microsatellite stable tumors, in that individuals diagnosed with colon cancer with microsatellite unstable tumors had improved survival while those diagnosed with rectal cancer and microsatellite unstable tumors had worse survival[16, 17]. It is possible that miRNAs that differ by microsatellite instability could explain differences in survival observed for colon and rectal cancer.

In this study we examine differences in expression of miRNAs in colon and rectal carcinoma based on tumor molecular phenotype. We assess unique miRNA expression associated with TP53-mutated, KRAS-mutated, BRAF-mutated, CpG island methylator phenotype, and microsatellite unstable tumors. An Agilent Array, containing over 2000 miRNAs was used to fully evaluate the range of miRNAs expressed in colorectal cancer tissue. We consider the association with survival for those miRNAs that are differentially expressed by tumor molecular phenotype to better understand the role of these miRNAs on prognosis after diagnosis with colorectal cancer.

Methods

Study Participants

The study was approved by the Institutional Review Board of the University of Utah; study participants signed informed consent. Study participants came from two population-based case-control studies that included all incident colon and rectal cancers between 30 to 79 years of age who resided along the Wasatch Front in Utah or were members of the Kaiser Permanente Medical Care Program in Northern California. Participants were white, Hispanic, or black for the colon cancer study; the rectal cancer study also included Asians and American Indians not living on reservations[19, 20]. Cases had to have tumor registry verification of a first primary adenocarcinoma of the colon or rectum and were diagnosed between October 1991 and September 1994 for the colon cancer and between June 1997 and May 2001 for the rectal cancer. Tumor tissue was obtained for 97% of all Utah cases diagnosed and for 85% of all Kaiser Permanente Medical Care Program in Northern California study participants[21]. Detail study methods have been described[22]. Individuals with familial adenomatous polyposis coli, inflammatory bowel disease, Crohn’s disease, and known Lynch Syndrome were excluded.

miRNA processing

RNA was extracted from formalin-fixed paraffin embedded tissue using Ambion’s RecoverAll Total Nucleic Acid isolation kit. Normal mucosa adjacent to the carcinoma tissue was used. The Agilent Human miRNA Microarray V19.0 was used given the number of miRNAs, its high level of reliability (repeatability coefficient was 0.98 in our data) and the amount of RNA needed to run the platform. The microarray contains probes for 2006 unique human miRNAs. A total of 100 ng total RNA was labeled with Cy3 and hybridized to the Agilent Microarray and were scanned on an Agilent SureScan microarray scanner model G2600D. Data were extracted from the scanned image using Agilent Feature Extract software v.11.5.1.1. Data were required to pass stringent quality control parameters established by Agilent that included tests for excessive background fluorescence, excessive variation among probe sequence replicates on the array, and measures of the total gene signal on the array to assess low signal. If samples failed to meet quality standards for any of these parameters, the sample was re-labeled, hybridized to arrays, and scanned. If a sample failed quality control assessment a second time the sample was deemed to be of poor quality and the individual was excluded from down-stream analysis. The Agilent platform was shown to have high repeatability (r=0.98) and good correlation with other platforms such as Nanostring[22]. Other analysis of a subset of important miRNAs showed 100% agreement in terms of directionality of differential expression and fold change (Pellatt 2016 In press).

Tumor molecular phenotype

We have previously assessed TP53 and KRAS mutations, the CpG island methylator phenotype using the classic panel[18], and microsatellite instability on the mononucleotide repeats at BAT26 and TGFβRII and a panel of 10 tetranucleotide repeats that were correlated highly with the Bethesda Panel[23]; our study was done prior to the development of the Bethesda Panel. The classic CpG island methylator phenotype panel consisted of MLH1, p16, and MINT1, MINT2, and MINT31. Tumors were scored as CpG island methylator phenotype high if two or more of the CpG islands were methylated; otherwise they were classified as CpG island methylator phenotype low. This panel was run prior to the advent of more recent panels[24, 25]. Additionally, the BRAF V600E mutation was assessed in colon carcinomas[15].

Tumor Registry Data

Both Utah and California are members of the National Cancer Institute funded Surveillance, Epidemiology, and End Disease Program. In addition to disease stage data, they provided follow-up data on all study participants which included total number of months survived, date of death or date of last follow-up, and cause of death. Follow-up was complete through 2006.

Statistical Methods

Of the 2006 unique human miRNAs assessed, 1278 were expressed in colorectal carcinoma tissue. To normalize differences in miRNA expression that could be attributed to the array, amount of RNA, location on array, or other factors that could erroneously influence expression, total gene signal was normalized by multiplying each sample by a scaling factor[26] (http://genespring-support.com/files/gs_12_6/GeneSpring-manual.pdf), which was the median of the 75th percentiles of all the samples divided by the individual 75th percentile of each sample_ENREF_11. We limited our analysis to miRNAs that were expressed in at least five carcinoma tissue samples and also had the mutated tumor phenotype of interest. Analysis ranged from 983 miRNAs assessed with TP53 in colon cancer to 694 miRNAs assessed with MSI in rectal cancer.

We assessed differences between tumor phenotype mutated and non-mutated, microsatellite unstable vs microsatellite stable and CIMP high vs. CIMP low samples for carcinoma tissue. We used log base 2 transformed miRNA expression levels using the significance analysis of microarrays technique implemented in the R package siggenes[27]; p-values were based upon 1000 permutations. To adjust for multiple comparisons, we applied a False Discovery Rate level of significance of 0.05 based on Benjamini and Hochberg[28] as implemented in siggenes. For those miRNAs that were significantly differentially expressed, we report the mean level of expression and the fold change (on non-log-transformed data) between the tumor phenotypes. We present statistically significant data where the fold change was 1.50 or higher or 0.67 or lower.

We further assessed the impact of differentially expressed miRNAs between mutated and non-mutated tumors on survival for rectal cancer and colon cancer separately. We calculated Cox proportional hazard ratios and corresponding 95% confidence intervals adjusting for age, sex, and AJCC stage. We also adjusted for microsatellite instability status since microsatellite instability has been shown to be associated with survival. Our statistical endpoint for our survival analysis was CRC cancer-specific survival based on months between diagnosis date and date of death or lost to follow-up. Individuals dying of other causes or who were lost to follow-up were censored at the time of their death or date of last contact. Because several miRNAs were infrequently expressed we calculated p values based on 1000 permutations for the HR using the R survival package[29, 30]. For miRNAs that were infrequently expressed, which we defined as expressed in less than 50% of individuals, we calculated the HR based on any expression. SAS 9.4 (SAS Institute, Cary, NC) was used to calculate the Cox proportional hazard ratios with the unit of change being the interquartile range. Adjustment for multiple comparisons was done by controlling the familywise error rate using the Sidak method[31]; we accounted for both the tests performed on the more commonly expressed miRNAs as well as on miRNAs infrequently expressed when calculating the familywise error rate.

Bioinformatics

Experimentally verified miRNA-target genes were identified for miRNAs associated with colon microsatellite instability tumor phenotype. Target genes were identified using the Homo sapiens download from miRTarBase v6.0 (http://mirtarbase.mbc.nctu.edu.tw/)[32]. R was used to filter the miRTarBase repository for genes associated with the miRNAs of interest, as well as compute the number of miRNAs each gene was associated with. The list of unique genes associated with these miRNAs comprised 9,404 target genes, ranging from genes having one miRNA association to 29 associations. Genes that were associated with seven miRNAs or more were kept for functional analysis; this number was chosen as it created a sub-list of 1,048 genes, which is suitable for gene set enrichment analysis. This subset of target genes was used as input to Ingenuity Pathway Analysis (http://www.ingenuity.com/)[33, 34] for functional analysis. Four genes were not mapped by Ingenuity Pathway Analysis: C1ORF21, C11ORF57, C5ORF51, and NDUFA4P1. We performed a core analysis, using Ingenuity Pathway Analysis knowledgebase of genes only and direct relationships only, included all sources, only experimentally verified results, mammalian species, and all tissues and mutations in the analysis with stringent filters. We exported the canonical pathways, where the enrichment score was corrected for multiple comparisons using the Benjamini Hochberg correction.

Results

A total of 1893 individuals had miRNA assessed; 60.8% were colon carcinomas and 39.2% were rectal carcinomas (Table 1). The mean age of study participants was 64.2 years and 45.8% of participants were female. For the 1862 individuals with AJCC stage available, 30.0% were Stage 1, 26.3% were Stage 2, 29.4% were Stage 3 and 14.3% were Stage 4. At the end of five years, 550 of 1893 participants had died of CRC.

Table 1

Description of Study Population

		Overall	TP53-mutated		KRAS-mutated		CpG Island Methylator Phenotype High		Microsatellite Instability		BRAF-mutated
		(N = 1893)	(N = 864)	0.01	(N = 569)	%	(N = 354)	(%)	(N = 170)	(%)	(N = 90)	(%)
Study Center	Kaiser Permanente	1144	507	0.44	336	0.29	222	0.19	107	0.09	58	0.05
	Utah	749	357	0.48	233	0.31	132	0.18	47	0.06	32	0.04
Tumor Site	Colon	1150	500	0.43	347	0.30	275	0.24	154	0.13	90	0.08
	Rectal	743	364	0.49	222	0.30	79	0.11	16	0.02	NE2
Age (mean)		64.2	63.8		64.0		66.2		66.8		67.6
Sex (% F)		45.8	44.0		45.5		54.5		63.0		54.4
AJCC Stage	1	559	274	0.49	139	0.25	85	0.15	34	0.06	15	0.03
	2	489	195	0.40	160	0.33	93	0.19	59	0.12	23	0.05
	3	548	258	0.47	168	0.31	123	0.22	50	0.09	35	0.06
	4	266	123	0.46	89	0.33	50	0.19	10	0.04	16	0.06
5 Year Survival	Alive	1081	487	0.45	307	0.28	194	0.18	107	0.10	49	0.05
	Dead Colorectal Cancer	550	253	0.46	194	0.35	117	0.21	23	0.04	33	0.06
	Dead Other	198	93	0.47	57	0.29	34	0.17	20	0.10	7	0.04
	Lost to Follow-up	62	31	0.50	10	0.16	9	0.15	4	0.06	1	0.02

% refers to percentage of row

NE indicates not evaluated

Several miRNAs were differentially expressed between mutated and non-mutated carcinomas using a false discovery rate of 0.05; those with a fold change of 1.5 or higher or 0.67 or lower in miRNA expression between mutated and non-mutated carcinomas are shown, although many more miRNAs were statistically significantly different in level of expression (Table 2). There were more differentially expressed miRNAs in TP53-mutated tumors than for either KRAS-mutated, CpG island methylator phenotype high, or BRAF-mutated tumors. The majority of miRNAs that were dysregulated were different between colon and rectal tumor sub-site. MiR-4700-3p was down-regulated in TP53-mutated tumors for both colon and rectal cancer. The majority of miRNAs dysregulated for TP53-mutated tumors were up-regulated, were infrequently expressed, and had low-levels of expression when expressed. Unlike for TP53, KRAS differentially expressed miRNAs were only observed for rectal cancer. The majority of miRNAs associated with KRAS-specific mutations were down-regulated in KRAS-mutated tumors, and were more frequently expressed than those differentially expressed for TP53. However, the three miRNAs that were up-regulated in KRAS rectal tumors were seldom expressed and had low levels of expression. Differentially expressed miRNAs for CpG island methylator phenotype high and BRAF-mutated tumors were only observed for colon cancer. MiR-196a-5p and miR-196b-5p were down-regulated in both CpG island methylator phenotype high and BRAF-mutated tumors.

Table 2

Differences in miRNA expression by tumor molecular phenotype with a false discovery rate set at 0.05

miRNA	% Expressing	Mean Expression	% Expressing	Mean Expression	Siggenes p value	Fold Change (mutated/not mutated)
Colon Cancer	TP53 Not Mutated		TP53 Mutated
hsa-let-7b-3p	10.6	2.1	15.4	4.22	0.00148	2.01
hsa-let-7d-3p	8.7	2.45	14	4.45	0.00282	1.81
hsa-miR-1226-3p	0.2	0.36	2.4	2.13	0.00165	5.89
hsa-miR-1227-3p	4.2	1.3	7.8	2.82	0.00174	2.17
hsa-miR-124-3p	20.9	1.33	15	0.85	0.01242	0.64
hsa-miR-1296	3.7	0.92	7.4	2.41	0.00154	2.61
hsa-miR-133b	12.4	1.47	19.4	2.34	0.00074	1.59
hsa-miR-1976	0.2	0.53	2.6	9.05	0.00059	17.07
hsa-miR-3190-5p	0.3	1.32	2.6	4.61	0.00238	3.49
hsa-miR-328	12.7	2.37	17.4	4.89	0.00213	2.06
hsa-miR-34b-5p	2.5	3.61	1.2	1.53	0.01979	0.42
hsa-miR-3613-3p	0.3	0.42	2.4	4.7	0.00306	11.27
hsa-miR-4469	27.1	1.7	20.8	1.12	0.00334	0.66
hsa-miR-449b-3p	4.5	2.71	8.8	5.04	0.0031	1.86
hsa-miR-4695-3p	2.7	1.51	5.6	3.5	0.00273	2.32
hsa-miR-4700-3p	13.9	2.06	9.8	1.35	0.0479	0.65
hsa-miR-4763-5p	13.6	1.54	18.2	3.03	0.00211	1.96
hsa-miR-513a-3p	13.1	0.84	9.2	0.49	0.02807	0.59
hsa-miR-625-5p	11.1	1.77	6.2	1.09	0.00439	0.62
hsa-miR-6509-3p	0.2	0.71	2.2	4.91	0.00289	6.88
hsa-miR-6511a-3p	3.9	2.24	7.4	4.81	0.00198	2.15
hsa-miR-885-5p	7.9	4.16	12.4	7.44	0.00246	1.79
Rectal Cancer
hsa-miR-127-5p	6.7	0.23	11.3	0.37	0.02098	1.65
hsa-miR-192-3p	23.9	1.42	36.5	2.21	0.00006	1.56
hsa-miR-3180	4.8	0.86	11.8	1.54	0.00022	1.79
hsa-miR-4700-3p	9.6	3.03	5.5	0.48	0.00375	0.16
hsa-miR-590-5p	11.8	1.13	20.6	1.81	0.00159	1.6
Rectal Cancer	KRAS Not Mutated		KRAS Mutated
hsa-miR-192-3p	32.9	2.1	21.6	1.1	0.00009	0.53
hsa-miR-196b-5p	77.9	16.36	73.4	10.54	0.00062	0.64
hsa-miR-204-3p	4.7	0.4	10.8	0.67	0.01043	1.66
hsa-miR-324-5p	48.6	4.37	32	2.74	0.00001	0.63
hsa-miR-374b-5p	38.2	6.09	30.2	3.57	0.00127	0.59
hsa-miR-4255	3.3	0.34	7.2	0.56	0.04833	1.66
hsa-miR-4290	28.7	3.92	18.5	2.48	0.00311	0.63
hsa-miR-518e-5p	0.2	1.23	2.7	2.05	0.03663	1.66
hsa-miR-532-5p	47.7	5.03	42.3	3.27	0.00491	0.65
hsa-miR-552	13.6	4.19	6.8	2.31	0.00009	0.55
hsa-miR-652-3p	25.6	1.88	17.1	1.21	0.00323	0.64
hsa-miR-660-5p	14.7	3.11	7.7	2.05	0.0002	0.66
hsa-miR-7-5p	54.3	6.59	45	4.44	0.0002	0.67
Colon Cancer	CpG Island Methylator Phenotype Low		CpG Island Methylator Phenotype High
hsa-miR-196a-5p	66.3	6.85	50.9	4.6	<0.00001	0.67
hsa-miR-196b-5p	82.7	18.03	68	9.4	<0.00001	0.52
hsa-miR-31-5p	9.9	7.31	32.7	12.12	<0.00001	1.66
hsa-miR-4267	3.7	0.59	7.6	1.37	0.026	2.33
hsa-miR-513a-3p	10.1	0.58	15.6	0.99	0.03441	1.71
hsa-miR-548d-5p	11.4	0.97	7.3	0.38	0.00389	0.39
hsa-miR-590-5p	23.6	2.28	16.4	1.14	0.0008	0.5
hsa-miR-625-5p	6.6	0.96	13.8	2.33	0.00063	2.43
Colon Cancer	BRAF Not Mutated		BRAF Mutated
hsa-miR-196a-5p	63.5	6.48	42.2	3.18	0.00042	0.49
hsa-miR-196b-5p	79.7	16.39	67.8	7.28	0.00027	0.44

Of the 94 differentially expressed miRNAs between microsatellite unstable and microsatellite stable colon tumors, all but eight were down-regulated (Table 3). The majority of the down-regulated miRNAs were frequently expressed. Additionally, the levels of expression appeared to be greater than those observed for most of the other tumor molecular phenotype differentially expressed miRNAs. Likewise, we observed that 41 miRNAs were significantly differentially expressed between microsatellite unstable and microsatellite stable rectal carcinomas with a fold change of 1.5 or greater or 0.67 or lower; almost all were downregulated (Table 4). Twenty-eight of the 41 miRNAs that were dysregulated between microsatellite unstable and microsatellite stable rectal carcinomas also were dysregulated between microsatellite unstable and microsatellite stable colon carcinomas. For the most part, the fold change observed between microsatellite unstable and microsatellite stable miRNA expression levels for rectal carcinomas was larger than observed for colon carcinomas.

Table 3

Associations between miRNAs differentially expressed between microsatellite stable and microsatellite unstable colon tumors with a false discovery rate set at 0.05

	Microsatellite Stable		Microsatellite Unstable
miRNA	% Expressing	Mean Expression	% Expressing	Mean Expression	Siggenes p value	Fold Change1
hsa-let-7a-5p	99.0	276.83	98.7	183.31	0.000	0.66
hsa-let-7e-5p	94.1	24.07	90.3	15.57	0.000	0.65
hsa-let-7f-5p	98.5	155.97	97.4	96.19	0.000	0.62
hsa-let-7g-5p	97.1	53.66	95.5	34.33	0.000	0.64
hsa-miR-100-5p	70.4	8.80	50.0	4.79	0.000	0.54
hsa-miR-125b-5p	97.3	53.77	96.8	30.40	0.000	0.57
hsa-miR-127-3p	15.0	2.70	5.2	1.15	0.000	0.43
hsa-miR-1291	76.8	5.92	59.1	3.86	0.000	0.65
hsa-miR-130a-3p	49.3	5.88	25.3	2.79	0.000	0.47
hsa-miR-143-3p	63.8	7.07	39.0	2.68	0.000	0.38
hsa-miR-145-5p	99.8	134.19	99.4	77.28	0.000	0.58
hsa-miR-148a-3p	79.4	17.83	60.4	8.17	0.000	0.46
hsa-miR-151a-3p	51.1	4.83	23.4	2.33	0.000	0.48
hsa-miR-151a-5p	91.2	14.36	80.5	6.85	0.000	0.48
hsa-miR-151b	83.5	6.95	66.9	3.51	0.000	0.51
hsa-miR-15a-5p	61.8	6.82	48.1	3.68	0.000	0.54
hsa-miR-17-5p	98.1	53.02	98.1	31.47	0.000	0.59
hsa-miR-183-5p	24.3	4.88	9.7	2.51	0.000	0.51
hsa-miR-185-5p	12.3	3.49	5.8	2.23	0.001	0.64
hsa-miR-1913	14.1	1.93	6.5	0.93	0.001	0.48
hsa-miR-192-3p	21.1	1.58	8.4	0.53	0.000	0.33
hsa-miR-192-5p	97.2	87.65	91.6	50.24	0.000	0.57
hsa-miR-193b-3p	84.6	8.50	64.9	4.59	0.000	0.54
hsa-miR-194-5p	97.3	84.55	92.9	54.19	0.000	0.64
hsa-miR-195-5p	44.1	3.23	20.8	1.12	0.000	0.35
hsa-miR-196a-5p	64.3	6.72	42.2	3.36	0.000	0.50
hsa-miR-196b-5p	81.0	17.51	51.9	4.72	0.000	0.27
hsa-miR-1973	99.8	793.25	100.0	512.18	0.000	0.65
hsa-miR-199a-3p	94.1	38.36	91.6	21.28	0.000	0.55
hsa-miR-199a-5p	89.5	17.08	81.8	8.86	0.000	0.52
hsa-miR-199b-5p	37.9	3.72	12.3	1.62	0.000	0.43
hsa-miR-19b-3p	91.9	25.37	87.7	14.21	0.000	0.56
hsa-miR-200a-3p	93.5	27.08	90.9	17.29	0.000	0.64
hsa-miR-200b-3p	98.4	151.75	98.1	97.72	0.000	0.64
hsa-miR-203a	66.4	11.82	56.5	7.58	0.000	0.64
hsa-miR-204-3p	15.0	1.12	9.7	0.71	0.044	0.64
hsa-miR-20a-5p	97.6	61.83	96.8	33.17	0.000	0.54
hsa-miR-20b-5p	86.7	15.29	79.9	7.72	0.000	0.50
hsa-miR-214-3p	88.7	11.97	81.2	7.47	0.000	0.62
hsa-miR-215	95.6	41.56	89.6	26.31	0.000	0.63
hsa-miR-221-3p	84.6	12.18	70.8	7.46	0.000	0.61
hsa-miR-223-3p	89.5	22.79	94.8	36.61	0.000	1.61
hsa-miR-224-5p	67.6	13.82	17.5	4.55	0.000	0.33
hsa-miR-23a-3p	99.1	157.85	98.7	93.21	0.000	0.59
hsa-miR-23b-3p	98.2	65.43	96.1	35.36	0.000	0.54
hsa-miR-26a-5p	98.7	96.40	98.1	55.09	0.000	0.57
hsa-miR-26b-5p	89.3	18.41	84.4	11.01	0.000	0.60
hsa-miR-27a-3p	97.6	49.86	96.8	30.52	0.000	0.61
hsa-miR-27b-3p	95.5	32.03	92.2	17.97	0.000	0.56
hsa-miR-28-5p	20.7	1.67	5.2	0.42	0.000	0.25
hsa-miR-29a-3p	98.1	97.72	98.1	54.89	0.000	0.56
hsa-miR-29b-3p	93.7	22.49	89.6	13.68	0.000	0.61
hsa-miR-29c-3p	85.5	18.35	77.9	10.14	0.000	0.55
hsa-miR-30a-5p	32.7	2.24	11.7	0.53	0.000	0.24
hsa-miR-30b-5p	94.4	20.87	90.3	10.81	0.000	0.52
hsa-miR-30c-5p	79.0	6.79	66.9	3.53	0.000	0.52
hsa-miR-30d-5p	99.8	32.98	99.4	21.46	0.000	0.65
hsa-miR-30e-5p	46.6	3.03	30.5	1.58	0.000	0.52
hsa-miR-31-5p	14.5	10.05	37.0	15.68	0.000	1.56
hsa-miR-330-3p	46.5	2.99	61.0	4.69	0.000	1.57
hsa-miR-3609	45.9	2.57	18.8	0.81	0.000	0.31
hsa-miR-361-3p	11.5	1.36	4.5	0.56	0.000	0.41
hsa-miR-361-5p	88.3	11.15	70.8	4.32	0.000	0.39
hsa-miR-362-5p	33.4	3.63	13.0	2.34	0.000	0.64
hsa-miR-3651	99.0	55.13	98.7	35.97	0.000	0.65
hsa-miR-3653	94.1	13.36	87.0	7.12	0.000	0.53
hsa-miR-365a-3p	63.7	7.91	37.0	4.18	0.000	0.53
hsa-miR-374b-5p	29.4	6.75	6.5	1.99	0.000	0.29
hsa-miR-424-5p	33.3	4.08	13.6	2.50	0.000	0.61
hsa-miR-425-5p	81.2	9.86	74.7	5.34	0.000	0.54
hsa-miR-4284	100.0	1352.71	100.0	873.15	0.000	0.65
hsa-miR-429	66.3	12.28	49.4	6.86	0.000	0.56
hsa-miR-4492	12.2	1.20	20.1	1.97	0.013	1.65
hsa-miR-455-3p	14.8	5.00	7.8	2.24	0.001	0.45
hsa-miR-4681	3.5	0.70	8.4	1.45	0.024	2.06
hsa-miR-483-3p	41.5	9.94	33.1	3.80	0.002	0.38
hsa-miR-497-5p	23.3	1.51	7.1	0.44	0.000	0.29
hsa-miR-505-5p	8.5	0.50	3.9	0.32	0.006	0.63
hsa-miR-520d-3p	12.4	1.62	22.1	3.05	0.014	1.88
hsa-miR-532-3p	27.9	3.10	6.5	1.57	0.000	0.51
hsa-miR-532-5p	41.1	5.74	7.8	1.62	0.000	0.28
hsa-miR-590-5p	23.0	2.14	5.8	0.34	0.000	0.16
hsa-miR-625-5p	7.6	1.27	16.9	2.47	0.003	1.95
hsa-miR-6515-5p	27.4	1.14	35.7	1.79	0.010	1.56
hsa-miR-652-3p	19.5	2.50	5.8	1.12	0.000	0.45
hsa-miR-664a-3p	39.5	4.27	18.2	2.27	0.000	0.53
hsa-miR-664a-5p	93.0	11.72	81.8	6.73	0.000	0.57
hsa-miR-664b-3p	99.7	31.27	100.0	16.29	0.000	0.52
hsa-miR-664b-5p	99.7	59.42	98.7	30.99	0.000	0.52
hsa-miR-92a-3p	99.8	107.15	99.4	52.95	0.000	0.49
hsa-miR-92b-3p	11.2	1.03	8.4	0.48	0.042	0.47
hsa-miR-98-5p	34.5	3.59	9.1	1.29	0.000	0.36
hsa-miR-99a-5p	49.8	5.52	22.1	2.07	0.000	0.38
hsa-miR-99b-5p	44.0	5.10	31.8	3.33	0.001	0.65

Fold Change is Microsatellite unstable/microsatellite stable

Table 4

Associations between miRNAs differentially expressed between microsatellite unstable and microsatellite stable rectal tumors with a false discovery rate set at 0.05

	Microsatellite Stable		Microsatellite Unstable
miRNA	% Expressing	Mean Expression	% Expressing	Mean Expression	Siggenes p value	Fold Change1
hsa-let-7a-5p	99.9	272.16	100	160.81	0.010	0.59
hsa-let-7f-5p	99.4	149.04	100	81.94	0.010	0.55
hsa-let-7g-5p	99	52.05	100	27.61	0.010	0.53
hsa-miR-100-5p	80.8	9.71	50	2.56	0.000	0.26
hsa-miR-125b-5p	98.9	62.09	100	28.37	0.010	0.46
hsa-miR-1281	100	42.21	100	63.57	0.010	1.51
hsa-miR-145-5p	100	133.78	100	47.27	0.000	0.35
hsa-miR-151a-5p	95.9	14.13	62.5	4.25	0.000	0.3
hsa-miR-151b	88.7	6.81	50	2.12	0.000	0.31
hsa-miR-196b-5p	77.6	14.95	37.5	1.33	0.000	0.09
hsa-miR-199a-5p	93.8	18.26	75	6.39	0.010	0.35
hsa-miR-221-3p	90.6	11.27	62.5	4.57	0.010	0.41
hsa-miR-23a-3p	100	164.28	100	83.89	0.000	0.51
hsa-miR-23b-3p	99.3	68.68	100	29.5	0.000	0.43
hsa-miR-24-3p	100	97	100	57.89	0.000	0.6
hsa-miR-26a-5p	100	97.29	100	40.86	0.000	0.42
hsa-miR-26b-5p	94.2	16.93	62.5	6.01	0.010	0.35
hsa-miR-27b-3p	98.5	31.49	87.5	14.04	0.010	0.45
hsa-miR-29c-3p	90.9	15.69	56.3	5.74	0.010	0.37
hsa-miR-30b-5p	97.6	21.37	75	8.42	0.010	0.39
hsa-miR-30c-5p	88.5	7.56	62.5	2.7	0.010	0.36
hsa-miR-30d-5p	99.9	31.73	100	18.7	0.000	0.59
hsa-miR-3150b-5p	97.1	12	100	18.48	0.000	1.54
hsa-miR-331-3p	98.8	13.8	75	6.55	0.010	0.47
hsa-miR-361-5p	93.5	11.38	56.3	2.74	0.000	0.24
hsa-miR-3653	97.8	13.88	75	6.64	0.010	0.48
hsa-miR-365a-3p	74.8	8.09	37.5	3.27	0.010	0.4
hsa-miR-424-3p	100	37.2	100	23.46	0.010	0.63
hsa-miR-425-5p	89.8	10.87	62.5	5.37	0.010	0.49
hsa-miR-4532	100	296.27	100	446.24	0.030	1.51
hsa-miR-4655-3p	98.5	14.89	100	22.35	0.010	1.5
hsa-miR-4674	91	5.41	100	10.55	0.020	1.95
hsa-miR-4769-3p	92.3	7.77	100	12.21	0.000	1.57
hsa-miR-520e	94.7	9.87	100	15	0.000	1.52
hsa-miR-5703	100	197.43	100	309.88	0.000	1.57
hsa-miR-636	99.7	8.97	100	13.59	0.010	1.51
hsa-miR-664b-3p	99.9	33.48	100	19.84	0.000	0.59
hsa-miR-664b-5p	99.9	58.55	100	35.64	0.000	0.61
hsa-miR-877-3p	100	24.65	100	37.46	0.000	1.52
hsa-miR-92a-3p	100	116.97	100	58.86	0.000	0.5
hsa-miR-940	100	591.89	100	895.9	0.010	1.51

Fold change is based on microsatellite unstable/microsatellite stable

We assessed if miRNAs significantly differentially expressed with a fold change of 1.5 or greater (or conversely 0.67 or less) were associated with survival (Table 5). Associations with survival could give insight into functionality as well as potentially explain differences in survival patterns for those with microsatellite unstable and microsatellite stable tumor phenotypes[16, 17]. There were few miRNAs that were associated with survival after adjustment for multiple comparisons along with age, sex, disease stage at time of diagnosis and even fewer significant associations after adjustment for microsatellite instability tumor phenotype. Prior to adjustment for microsatellite instability status, there were four miRNAs associated with microsatellite instability status and one miRNA associated with CpG island methylator phenotype status in colon tumors that was associated with colorectal cancer survival. For rectal cancer, there was one miRNA associated with KRAS, three associated with microsatellite instability, and one associated with TP53 mutations in rectal tumors that also were associated with colorectal cancer survival after diagnosis with rectal cancer. After adjustment for microsatellite instability, there was one differentially expressed miRNA between colon BRAF-mutated and non-mutated tumors and two miRNAs differentially expressed between colon microsatellite unstable and microsatellite stable associated with survival after diagnosis with colon cancer. Additionally, one miRNA associated with microsatellite unstable and microsatellite stable in rectal tumors and one miRNA associated with TP53 were associated with colorectal cancer survival after diagnosis with colon cancer. Interestingly, miR-196b-5p associated with BRAF mutations in colon cancer was only associated with colorectal cancer survival after adjustment for microsatellite instability status. For those that remained significant after adjustment for microsatellite status, the hazard ratio changed from 1.49 to 1.70 for miR-31-5p after adjustment for microsatellite instability status; most of the other hazard ratio only changed marginally if at all.

Table 5

Associations between miRNAs dysregulated for colon and for rectal cancer and colorectal cancer survival

Survival	Original Tumor							P-value
Association	Phenotype	miRNA	% Expressing	25th %ile1	75th %ile1	Hazard Ratio2	(95% Confidence Intervals)	Raw2	Adjusted3
Colon	CIMP and MSI4	hsa-miR-31-5p	14.8	0.00	0.00	1.49	(1.16,  1.91)	0.002	0.01
Colon	MSI	hsa-miR-193b-3p	86.8	1.92	3.40	1.28	(1.14,  1.44)	0.0001	0.009
Colon	MSI	hsa-miR-214-3p	90.9	2.36	3.99	1.26	(1.11,  1.43)	0.0003	0.03
Colon	MSI	hsa-miR-28-5p	21.0	0.00	0.00	1.62	(1.26,  2.08)	0.0002	0.02
Colon	MSI	hsa-miR-99b-5p	47.3	0.00	2.54	1.55	(1.25,  1.92)	<.0001	0.005
Rectal	KRAS	hsa-miR-552	9.7	0.00	0.00	0.47	(0.28,  0.77)	0.003	0.03
Rectal	MSI	hsa-miR-3653	94.8	3.28	4.19	0.80	(0.71,  0.90)	0.0006	0.02
Rectal	MSI	hsa-miR-520e	91.4	3.00	3.68	1.21	(1.07,  1.37)	0.001	0.04
Rectal	MSI	hsa-miR-636	99.4	3.00	3.44	1.25	(1.11,  1.41)	0.0005	0.02
Rectal	TP53	hsa-miR-192-3p	23.1	0.00	1.64	0.65	(0.48,  0.88)	0.005	0.03
Additional Adjustment for MSI
Colon	BRAF	hsa-miR-196b-5p	76.8	1.58	4.24	0.83	(0.71,  0.97)	0.02	0.04
Colon	CIMP and MSI	hsa-miR-31-5p	14.8	0.00	0.00	1.70	(1.32,  2.19)	<.0001	0.004
Colon	MSI	hsa-miR-99b-5p	47.3	0.00	2.54	1.59	(1.28,  1.97)	<.0001	0.003
Rectal	MSI	hsa-miR-636	99.4	3.00	3.44	1.23	(1.09,  1.39)	0.001	0.04
Rectal	TP53	hsa-miR-192-3p	23.1	0.00	1.64	0.65	(0.48,  0.89)	0.007	0.03

Expression values for the 25th and 7th are expressed as log2 transformed values that were used in the Cox proportional hazard models

p values adjusted for age, sex, and AJCC Stage

Adjusted for multiple comparisons

CIMP=CpG Island Methylator Phenotype and MSI = Microsatellite Instability status

Assessment of canonical pathways associated with differentially expressed miRNAs between microsatellite unstable and microsatellite stable tumors showed hundreds of pathways significantly enriched by genes regulated by differentially expressed miRNAs (Supplemental Table 1). Targeted pathways include those related to molecular mechanisms of cancer, estrogen-mediated signaling, cell cycle regulation, PI3K/AKT-signaling, PTEN-signaling, TP53-signaling, IGF-1-signaling, and TGFβ-signaling. The top 50 pathways are displayed in Figure 1.

Figure 1

Discussion

There were few differentially expressed miRNAs by TP53, KRAS, CpG island methylator phenotype, and BRAF molecular phenotype for either colon or rectal carcinoma. Conversely, 94 miRNAs were differentially expressed between microsatellite unstable and microsatellite stable tumors for colon carcinomas and 41 miRNAs were differentially expressed between microsatellite unstable and microsatellite stable tumors for rectal carcinomas. For the most part, miRNA expression was down-regulated in microsatellite unstable tumors. Evaluation of those miRNAs that were associated with specific tumor molecular phenotype with survival showed few associations. Those that were associated with survival were originally identified as being differentially expressed between microsatellite stable and microsatellite unstable tumors. MiRNAs that influenced colon cancer survival increased the likelihood of dying, while those associated with rectal cancer improved survival.

In interpreting the findings from this study, several things should be considered. We only reported associations for miRNAs that were both statistically significant and had a fold change of 1.5 or greater or 0.67 or lower. This excluded numerous associations that were statistically significant, and it also excluded associations that had lower differences in expression between mutated and not-mutated tumors. While this is a common practice, it is not without problems. First, we do not know what level of fold change is biologically meaningful. Second, the level of miRNA expression can impact the fold change. This is most evident in our study when looking at fold changes for TP53-mutated, KRAS-mutated and CpG island methylator phenotype high tumors. The larger fold changes are often associated with low levels of miRNA expression and miRNAs expressed in few individuals; this can yield a large fold change that may not be biologically meaningful. In our data, larger fold changes associated with MSI and MSS tumors may be more meaningful given miRNA expression levels for that specific tumor phenotype were more frequently expressed and had higher levels of expression. Additionally we are looking at miRNAs that uniquely associated with tumor phenotype, thus some miRNAs, such as miR-21 being significantly differentially expressed in our overall colorectal cancer data[22], was not associated with any one tumor molecular phenotype.

Comparison of our data to the literature on miRNA and tumor molecular phenotype is difficult given most studies that suggest associations between specific tumor phenotype and miRNAs have been conducted on cell lines or cells from cultured tissue and were limited to targeted miRNAs. Those conducted in populations have generally studied few people and often evaluated expression only compared to normal mucosa, often do not report specificity to a given phenotype, and may have samples limited to specific disease characteristics such as Crohn’s Disease. Despite these differences in methodology, there are suggestions that certain miRNAs can be associated with tumor molecular phenotype. This study gives us an opportunity to test some of the previously hypothesized associations. A study of six patients with either Crohn’s Disease or ulcerative colitis and expression of 88 miRNAs identified six differentially expressed miRNAs (miR-122, mIR-214, miR-372, miR-15b, let-7e, and miR-17) among TP53-mutated tumors[35]. Others have found miR-34 associated with the TP53 network[36, 37]; miRNA-16[38], miR-1915[39], miR-221[40], and miR-148b[41] also have been linked to TP53 mainly in cell lines or tissue cultures. We only found miRNA-34b-5p associated with TP53-mutated colon tumors in our data regardless of level of fold change. It is believed that miR-34 is involved in apoptosis and cell-cycle arrest[42].

We examined TP53-mutated tumors compared to non-mutated tumors previously in a small set of cases and this study replicates some of the previously reported associations[10]. For those that did replicate, which included miR-135b for rectal carcinoma and miR-224, miR-17, miR-1226, miR-532-5p, miR-17, miR-574-5p, miR-424, and miR-16 for colon cancer, the direction of association was the same as previously reported. Inconsistent nomenclature between an earlier versions of the Agilent array we used previously in our pilot study and to the much larger platform used in this study, prohibited comparison of some miRNAs. Likewise, we adjusted for more comparisons because of the larger platform in the current study and many of the previous associations that did not replicate were marginally significant with adjusted p values between 0.04 and 0.05. We identified 22 miRNAs associated with TP53-mutated colon tumors and five miRNAs associated with TP53-mutated rectal tumors in this study. With one exception the differentially expressed miRNAs for colon and rectal TP53-mutated tumors were different.

Several miRNAs also have been examined in conjunction with KRAS-mutated and non-mutated tumors. Among these, let-7, miR-100, miR-126, miR-143, miR-145, miR-200c, miR-221, miR-222, miR-224, miR-345, miR-4689, miR-96-5p have been linked to KRAS-mutations mainly in studies of cell lines and cell cultures[43-53]. MiR-31-5p has been associated with both KRAS and BRAF mutations[14, 54]; miRNA-31-5p was seen as the most up-regulated miRNA associated with the BRAF V600E mutation in a study of 760 miRNA in 29 CRC cases[14]. This miRNA also has been associated with CIMP status and serrated polyps[54]. BRAF also has been identified as a direct target of miR-378-5p[55]. In our data, miR-31-5p was significantly upregulated for CpG island methylator phenotype high colon carcinomas with a fold change of 1.66 and in BRAF-mutated tumors, however the fold change of 1.17, although statistically significant, may have less importance biologically. On the miRNAs previously linked to CpG island methylator phenotype, we observed statistically significant small changes in expression levels resulting in small fold changes for miR-145-5p (fold change=0.90), miR222-3p (fold change=1.30) and miR-224-5p (fold change=0.91). Our previous study identified no miRNAs that were associated with KRAS in colon carcinomas and 11 that were associated with rectal carcinoma when compared to non-KRAS-mutated tumors[10]. In this study only two miRNAs, miR-572 and miR-638, which although statistically significant in our data had low fold changes (fold change 1.16 and 1.11 respectively), showed any association with KRAS-mutated tumors. None of the other miRNAs linked in the literature to KRAS were associated with KRAS-mutated tumors in our data. Likewise, miR-378-5p was not associated with colon BRAF-mutated tumors. We replicated our previously identified association between CIMP with miR-31 and miR-492 from our smaller pilot data, however the fold change for miR-492 was small (fold change=1.26). We currently did not observe any CpG island methylator phenotype specific associations for rectal tumors, whereas previously we identified eight miRNAs that were associated with CpG island methylator phenotype high tumors, although the adjusted p values were marginal at 0.047.

Most previous work with miRNAs and tumor molecular phenotype has been with microsatellite unstable tumors. In our data the majority of associations also were observed for differences in miRNA expression between microsatellite unstable and microsatellite stable carcinomas; this was true for both colon and rectal cancer. For the most part, miRNAs were down-regulated in microsatellite unstable tumors. Others have shown that expression of miR-92, miR-223, miR-155, miR-196a, miR-31, and miR-26b were significantly different in a group of 23 microsatellite stable and 16 microsatellite unstable tumors[56]. They also showed that miR-31 and miR-223 were overexpressed in colorectal cancer in patients with hereditary non-polyposis colorectal cancer syndrome and that the oncogenic miR-17-92 family were significantly up-regulated in microsatellite stable cancers[56]. Similar findings have been reported for individuals with Lynch syndrome[57]. Others have reported that miR-484 expression was significantly decreased in microsatellite unstable colorectal cancers[58]. MiR-155 has been associated with MSI in several studies[58-61]. Several studies also have shown that miR-21 was over-expressed in microsatellite unstable tumors[58, 59]. In our previous work we also evaluated microsatellite unstable vs microsatellite stable tumors and observed eight miRNAs that were significantly up or down-regulated[10]; all but one of those miRNAs (miR-552) was significantly associated with microsatellite unstable carcinomas in this larger study.

Microsatellite unstable carcinomas appear to have the most unique associations with miRNA expression of any tumor molecular phenotypes examined. Given the number of miRNAs that were differentially expressed between microsatellite unstable and microsatellite stable carcinomas for both colon and rectal cancer, we further evaluated canonical pathways that were influenced by genes regulated by miRNAs specific to microsatellite unstable carcinomas. While many of the pathways that were enriched by genes associated with miRNAs differentially expressed between microsatellite unstable and microsatellite stable carcinomas were directly cancer related, such as Molecular Mechanisms of Cancer, others were involved in estrogen signaling, TP53-signaling (noted by IPA as p53), PTEN-signaling IGF1 (noted by IPA as IGF-1) signaling, TGFβ-signaling, STAT3 pathway, (MAPK7) ERK5-Signaling and WNT/β catenin -signaling. This reflects the wide range of pathways and mechanisms that could be associated with microsatellite unstable tumors. Many of these pathways have been previously identified as being associated with miRNAs[62]. Unfortunately, many of the bioinformatics databases available to assess functionality are limited in their specificity. Thus many of these pathways could be associated with miRNAs and colorectal as well as with microsatellite unstable tumors specifically.

For those miRNAs that were differentially expressed by tumor molecular phenotype, we tested their association with survival. The motivation for this was the observation that microsatellite unstable tumors are associated with survival differently than microsatellite stable tumors and that these associations differ for colon and rectal cancer[16, 17]. For the most part, we observed few associations between differentially expressed miRNAs by tumor phenotype and impact on survival. Those associated with colon microsatellite unstable carcinomas had the greatest impact on survival. Differences in associations between colon and rectal cancer could be from importance of different pathways associated with miRNAs and microsatellite unstable carcinomas for colon and rectal cancer. Some of these miRNAs have previously been associated with either advanced stage or survival. Others have reported a direct association between miR-31 with advanced disease stage and worse survival[14, 62]. Our replication of previously reported associations with disease stage and survival using half of the existing dataset did not show an association with miR-31-5p after adjusting for disease stage while miR-99b-5p was associated with worse survival among individuals diagnosed with colon cancer[63]. Given the infrequent expression of miR-31-5p in colorectal carcinomas we had limited power to evaluate associations with the smaller sample size. Adjustment for microsatellite instability status attenuated some of the associations with survival, while others became stronger. We presented both adjusted and not-adjusted microsatellite instability associations. It is unclear however if microsatellite instability is a confounder of associations in that it is associated with both the exposure and the outcome, or if it is in the causal pathway given that these miRNAs were associated with different levels of expression between microsatellite stable and microsatellite unstable carcinomas which could lead to tumor phenotype which then influenced survival.

This study has strengths, including the use of the Agilent platform with over 2000 miRNAs. We have previously reported on the comparability of the Agilent to Nanostring and qPCR, with excellent results[22]. The reliability of the Agilent platform in our data was 0.98 and the overall platform had relatively good performance when compared with Nanostring[22]. Comparison of expression levels of a set of specific miRNAs with qPCR showed extremely high correlation in terms of directionality of expression and fold change observed between carcinoma and normal mucosa. Additionally, the study is large and contains detailed information on multiple tumor molecular phenotypes as well as survival. In this study, we applied an FDR of 0.05 to take into consideration multiple comparisons being made. However, ideally another large dataset would replicate these findings. In previous analyses, we have split our dataset to incorporate both discovery and testing; we chose not to do that in these analyses to maintain more power to identify associations for those tumor molecular phenotypes that are rarer, such as MSI. We thus encourage others to replicate these findings in a more targeted approach that focuses on specific miRNAs, especially those that have prognostic implications.

We believe that this study makes a unique contribution to the literature. First, it is large so that associations are measured with more precision. Second, we are able to incorporate discovery of new associations as well as test those previously identified as being associated with specific miRNAs. In doing this we have adjusted for multiple comparisons, where others looking at candidate miRNAs have not, and thus, many previous associations may have been false positives that do not withstand the adjustment for multiple comparisons. This is not unlike the transition from candidate SNP and candidate pathways to genome-wide association studies, where more significant associations are identified when fewer adjustments for multiple comparisons are made.

In summary, our data suggest that most unique associations between tumor molecular phenotype and miRNAs are with microsatellite unstable vs. microsatellite stable tumors. While several miRNAs were differentially expressed between these tumor phenotypes, few of the miRNAs were associated with survival. Microsatellite unstable colon carcinomas have previously been associated with better survival, while microsatellite unstable rectal tumors were associated with worse survival. Interestingly, miRNAs in this study uniformly show that higher levels of expression increase risk of dying for colon cancer, but improve survival if diagnosed with rectal cancer. We encourage others with similar datasets to confirm these associations.