Comprehensive analysis of negatively correlated miRNA-mRNA regulatory pairs associated with microsatellite instability in colorectal cancer

Abstract

BACKGROUND:

Several studies have demonstrated that microRNAs (miRNAs) and target mRNAs are associated with different frequencies of microsatellite instability.

OBJECTIVE:

The study aimed to elucidate the profiles of miRNAs and target mRNAs expression and their associations with the phenotypic hallmarks of microsatellite instability in colorectal cancers (CRC) by integrating transcriptomic, immunophenotype, methylation, mutation, and survival data.

METHODS:

Differentially expressed miRNAs (DEmiRNAs) and mRNAs (DEmRNAs) were screened out and then the miRNA-mRNA regulatory pairs were identified through two databases. We verified that the expression levels were detected in 40 microsatellite instable (MSI) and 40 microsatellite stable (MSS) CRC samples and used the logistic regression and the Cox regression method to evaluate the diagnostic and prognostic value of negative regulatory pairs respectively.

RESULTS:

The best diagnostic model that combines miR-31-5p, PLAGL2, miR-361-5p, and RAB27B, which were associated with immune microenvironment, tumor mutation burden (TMB), and overall DNA methylation, could significantly predict microsatellite instability in colon tissues. MiR-31-5p and RAB27B could also predict the overall survival of MSS CRCs.

CONCLUSION:

This study generated a predictive model of the combination of miRNAs and mRNAs to distinguish MSI versus MSS CRCs and elaborated their potential molecular mechanisms and biological functions.

Keywords

miRNA miRNA-mRNA networks microsatellite instability colorectal cancer

Figure 1.

Flow chart for identifying the miRNA-mRNA regulatory pairs and the comprehensive analysis of regulatory pairs role of microsatellite instability in colorectal cancer.

1. Introduction

Colorectal cancer (CRC) where the incidence ranked second has a third largest estimated mortality for all types of cancers worldwide in 2020 [1]. The molecular subtyping has been suggested to be important for treatment plan selection and prognosis prediction [2]. Microsatellite instability (MSI) which is one of the classicized bases of the current widely accepted consensus molecular subtype (CMS) has been observed in CRCs and endometrial cancers [3]. Microsatellites (MS) which are small tandem repeats of DNA bases and extensive throughout the human genomes can lead to genomic instability and mutation due to the DNA slippage [4]. MSI can be divided into three groups: high-frequency MSI (MSI-H), low-frequency MSI (MSI-L) and microsatellite stable (MSS) depending on the degree [4].

MSI detection is crucial for diagnosing Lynch syndrome and selecting CRC patients eligible for immunotherapy or basic chemotherapy. A panel of two mononucleotides (BAT-25 and BAT-26) and three dinucleotides (D2S123, D5S346, and D17S250) repeat microsatellites called the NCI panel for the detection of MSI in CRC is considered as a golden standard in many countries [5]. In addition, a whole range of different clinical features and molecular biomarkers to differ microsatellite instable samples from CRCs are available such as tumor infiltrating lymphocytes [6].

MicroRNAs (miRNAs) had been reported that were critical for the progression of cancers by influencing proliferation, invasion and metastasis [7]. The interrelationship between MSI and miRNA plays a key role in the pathogenesis of CRC. Serval studies have demonstrated some miRNAs and their target mRNAs associated with CRC using bioinformatics analysis or microarray expression profile or a multi-step approach such as mir-29a/KLF4 [8], miR-182/SEMA6D [9], miR-183-5p/PNOC and miR-6730-5p/NUSAP1 [10]. By Integrative analysis of paired miRNA-mRNA expression profiles from CRC samples, we identified some miRNAs that may have a driving role at the CRC transcriptional level. An overview of the workflow steps is provided in Fig. 1. In our study, gene and miRNA profiling data were collected from the Gene Expression Omnibus (GEO) database and the Cancer Genome Atlas (TCGA). Expression of a selected panel of miRNA-mRNA pairs was evaluated by poly(A) reverse transcriptase quantitative (real-time) polymerase chain reaction (RT-qPCR) assay of formalin-fixed paraffin-embedded (FFPE) samples obtained from MSS and MSI (including MSI-L and (MSI-H) CRC. We identified the key miRNAs and their target mRNAs with survival analysis, functional enrichment and provided new insights into miRNA and target-mRNA expression in relation to tumor phenotypes which contain immune microenvironment, tumor infiltrating immune cells, global methylation, tumor mutational burden and genetic mutation.

Our research performed extensive analysis of miRNA -mRNA regulatory pairs in microsatellite instability in CRC to yield new sights in the underlying mechanism in CRC tumorigenesis associated with microsatellite instability. Combination of bioinformatic analysis and qRT-PCR provided with convenience in identifying dysregulated miRNA-mRNA regulatory pairs to improve therapeutic strategies for MSI CRCs and MSS CRCs.

2. Materials and methods

2.1 Data acquisition and processing

TCGA colon adenocarcinoma (COAD) and rectal adenocarcinoma (READ) miRNA and mRNA sequencing expression profiles and associated clinicopathological data were downloaded from the GDC data portal. There is no apparent difference between colon and rectal samples validated by Principal component analysis (PCA) and merging samples is no need to adjust [11]. So, we combined TCGA-COAD and TCGA-READ samples into a single TCGA-CRC cohort. TCGA-CRC samples were classified into MSI (MSI-H and MSI-L) and MSS groups based on data provided by Isidro Cortes-Ciriano et al. [12]. Then we searched MSI relevant gene microarray expression datasets from GEO and differentially expressed miRNAs in published articles. High-throughput TCGA-CRC data and microarray expression GEO datasets were analyzed by “R-edgeR” and GEO2R tool respectively.

2.2 Function annotation and pathway enrichment analysis of miRNA/target-gene pairs

Firstly, we screened the target genes of potential DE-miRNAs using DIANA-TarBase [13] and miRTarBase [14], which contain experimentally verified miRNA/ target-gene pairs. Then the expression correlations between DEGs and DEMs were evaluated using Pearson’s correlation. Visualization of the miRNA-mRNA regulatory network was conducted using Cytoscape software (v3.8.0) [15]. Gene ontology (GO) functional analysis and a Kyoto Encyclopedia of Genes Genomes (KEGG) pathways analysis against the DEGs and DEMs in the network were performed by using clusterProfiler tool in Hiplot web platform [16].

2.3 Connectivity map analysis of potential compounds capable of targeting the differentially expressed genes

Connectivity Map (cMap) analysis was conducted by querying up-regulated genes tag and down-regulated genes tag in CRCs versus NCs for finding perturbagens that give rise to similar expression signatures via clue.io software platform. The negative normalized connectivity score ranging from $-$ 3 to $-$ 1.8 indicated a suppressive effect of compound on the up-regulated genes. Negative log10 transformed FDR $q$ -values (fdr_q_nlog10) $>$ 2 was set as the filter condition.

2.4 Survival analysis

Univariate and multivariate Cox regression analyzed by “R-survival” [17]. Based on the cut-off point of miRNA or mRNA expression by “R-survminer”, patients also could be classified into high- expression and low- expression. The Kaplan-Meier analysis (KM) and the receiver operation characteristic (ROC) curve were carried out using the “R-KMsurv”.

2.5 Evaluation of interactions of miRNA-mRNA pairs and tumor-relative phenotypes and gene mutation status

The fraction of 22 infiltrating immune cell types was calculated using CIBERSORT (https://cibersort.stanford.edu/index.php/) [18]. The differences of these immune cells between MSI TCGA-CRC and MSS TCGA-CRC were compared via the Wilcoxon rank-sum test. ESTIMATE software based on the mRNA-seq data was used to estimate the stromal and immune levels for TCGA-CRC samples [19]. The methylation levels of the CpG sites in TCGA-COAD and TCGA-READ samples were obtained using UCSC Xena platform detected by the Illumina Infinium HumanMethylation450 BeadChips platform. The sum of the methylation levels of all 485,577 CpG sites in each sample was calculated as overall DNA methylation level. Tumor mutation burden (TMB) was used to measure the total number of somatic variants per megabase (MB) of the genome using the “R-maftools” [20]. TCGA-CRC samples were grouped into wild-type and mutated groups based on somatic mutation datasets obtained from UCSC Xena.

2.6 Sample collection and RNA isolation

Five MSI markers of the NCI panel were detected by PCR. 80 FFPE CRC tissue samples collected at surgical resection were classified as 40 MSI samples that one or more markers showed allelic shift and 40 MSS samples that none showed altered allelic size in tumor DNA. All samples used in this study were collected with patients’ consent. The present study has been approved by the institutional ethics committee and the patients written informed consent has been obtained (ID: 2016-SRFA-148). The clinical characteristics were showed in Table 1. Total RNA was extracted from FFPE colon tissue samples using RNAprep Pure FFPE Kit (TIANGEN) according to the manufacturer’s instructions. RNA concentrations were measured with NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific).

Table 1
Clinicopathological and molecular features of colorectal cancer patients

	MSI	MSS	$P$ value
Age (year), mean (SD)	57.33 $\pm$ 13.38	61.93 $\pm$ 13.04	0.123
Gender, $n$ (%)			0.648
Female	17 (42.5)	15 (37.5)
Male	23 (57.5)	25 (62.5)
Tumor location, $n$ (%)			0.061
Left colon	7 (17.5)	11 (27.5)
Right colon	25 (62.5)	23 (57.5)
Rectum	8 (20)	6 (15)
Grade, $n$ (%)			0.141
G1	0 (0)	1 (2.5)
G2	23 (57.5)	29 (72.5)
G3	17 (42.5)	9 (22.5)
TNM stage, $n$ (%)			0.085
I	4 (10)	9 (22.5)
II	23 (57.5)	15 (37.5)
III	12 (30)	11 (27.5)
IV	1 (2.5)	5 (12.5)
KRAS, $n$ (%)			0.264
Wild-type	18 (45)	21 (52.5)
Mutation	22 (55)	19 (47.5)
BRAF V600E, $n$ (%)			0.117
Wild-type	36 (90)	36 (90)
Mutation	4 (10)	4 (10)

2.7 Quantitative reverse transcription PCR (qRT-PCR) validation

Selected DEGs and DEMs were verified by qRT-PCR using PrimeScript RT reagent Kit (Takara) and SYBR Premix Ex Taq II (Takara) after adding a poly(A) tail to RNA by Poly(A) Polymerase Kit (Takara) according to the manufacturer’s protocol. The PrimeScript RT reagent Kit and SYBR Premix Ex Taq II kit contain the commercial Uni-RT primer and Uni-Reverse primer. The PCR reactions were carried out in final volumes of 10 $\mu$ L on the qTOWER [3] 84 (Analytik Jena). The sequences of PCR primers are listed in the Table S1. RUN6B (U6), GAPDH and 18S rRNA were considered as reference genes for normalization, and the comparative cycle threshold (2 ${}^{\rm-\Delta\Delta Ct}$ ) method was used to analyze the relative expression of miRNAs and genes by Livak KJ et al. [21].

2.8 Statistical analysis

IBM SPSS Statistics v.26 software (IBM Corporation, Armonk, NY, USA), R language v3.6.3 or GraphPad Prism software were used to analyze the data. Continuous variables were reported using the mean and standard deviation (SD). Student’s t tests were performed, and $p$ -values were calculated. MiRNA and gene with a $|\log_{2}\rm FC|>1$ and $P<$ 0.05 were considered to be statistically significant. The association between the expression of miRNAs and genes was analyzed by Pearson’s correlation in MSI and MSS CRC tissues. The diagnostic value of miRNA-mRNA pairs was assessed by the area under the ROC curves (AUC), decision curve analysis (DCA) and clinical impact curves which were used to evaluate the predictive accuracy and clinical value of the predictive models. Pearson’s correlation method was used to calculate correlation between DEGs or DEMs and all tumor-related phenotypes. Survival curves plotted by the KM method were analyzed by the log-rank test and $p<$ 0.05 was regarded as statistically significant. Cox regression analysis were performed and the hazard ratio (HR) and 95% confidence interval (CI) were calculated to identify statistically significant DEGs or DEMs associated with overall survival (OS).

3. Results

3.1 Identification of differentially expressed miRNAs (DEMs) and genes (DEGs) in colorectal cancer

There were 11 gene expression microarray datasets from the same manufacturer and 2 miRNA expression microarray datasets and DEMs mentioned in 15 articles related with MSI CRC filtered out shown in Fig. 2A, and listed in Tables 2 and 3.

Figure 2.

A: The circular-barplot showing the basic information of GEO datasets (GEO accession, source name and experiment type). A total of 15 microarray datasets from tissue were included in the study. 38 DEMs and 108 DEGs were screened in microsatellite instable colorectal cancers (MSI CRCs) versus microsatellite stable colorectal cancers (MSS CRCs). B: 38 microRNAs (miRNAs) to 108 mRNAs network visualized by Cytoscape. 65 miRNA (up) - mRNA (down) pairs and 74 miRNA (down) - mRNA (up) pairs were screened out by miRtarbase and Tarbase which contain experimentally validated miRNA-mRNA regulatory pairs. Blue dot represents the upregulated miRNAs/mRNAs in MSI CRCs versus MSS CRCs, while orange dot represents the downregulated miRNAs/mRNAs in MSI CRCs versus MSS CRCs.

Table 2

Information pertaining to the selected GEO datasets for colorectal cancer

[l]Study type(array)	GEO accession	Platform	Group
			MSI-H	MSI-L	MSS
Gene expression	GSE11543	GPL80	13	/	38
	GSE13067	GPL570	11	/	63
	GSE13294	GPL570	78	/	77
	GSE18088	GPL570	19	/	34
	GSE2138	GPL96	20	/	20
	GSE26682	GPL96	17	24	99
		GPL570	18	23	119
	GSE29638	GPL5175	10	1	34
	GSE35896	GPL570	5	/	56
	GSE4554	GPL570	33	/	51
	GSE75317	GPL570	11	/	48
		GPL5175	24	/	182
	GSE92921	GPL570	5	/	53
microRNA expression	GSE30454	GPL8179	28	/	19
microRNA expression	GSE136194	GPL16384	9	/	9

Table 3

Differentially expressed miRNAs mentioned in articles about microsatellite instability in colorectal cancer

PMID	Study type	Source	Samples		Regulation (MSI vs. MSS)
			MSI-H	MSS	Up	Down
17716371	Microarray	Tissue	16	23	hsa-miR-155, hsa-miR-223	hsa-miR-17-5p, hsa-miR-20, hsa-miR-25, hsa-miR-92-1, hsa-miR-92-2, hsa-miR-93-1, hsa-miR-106a
18676867	qRT-PCR microarray	Tissue	12	37	/	hsa-miR-142-3p, hsa-miR-212-3p, hsa-miR-144-3p, hsa-miR-151a-3p
20413677	qRT-PCR	Tissue	22	33	hsa-miR-155, hsa-miR-31, hsa-miR-233, hsa-miR-26b	hsa-miR-196a, hsa-miR-92, hsa-miR-145, hsa-let-7a
21844009	qRT-PCR microarray	Tissue	33	19	hsa-miR-622, hsa-miR-938, hsa-miR-615-5p, hsa-miR-1184	hsa-miR-17-5p, hsa-miR-192-3p, hsa-miR-337-3p
19922656	qRT-PCR microarray	Tissue	12	68	hsa-miR-625, hsa-miR-31	hsa-miR-224, hsa-miR-181c, hsa-miR-196b, hsa-miR-552, hsa-miR-592, hsa-miR-328, hsa-miR-133a, hsa-miR-9, hsa-miR-1
31784758	qRT-PCR microarray	Tissue	57	44	hsa-miR-181a-2-3p	/
30786859	HiSeq	Tissue	18	30	hsa-miR-7641, hsa-miR-6087, hsa-miR-625-3p	hsa-miR-335-5p, hsa-miR-26a-5p
29735329	qRT-PCR microarray	Tissue	67	54	hsa-miR-20b, hsa-miR-23a, hsa-miR-30e-5p, hsa-miR-141, hsa-miR-597, hsa-miR-601	hsa-miR-181a, hsa-miR-135a, hsa-miR-302c, hsa-miR-484
28734015	qRT-PCR microarray	Tissue	27	notfound	hsa-miR-21-5p, hsa-miR-21-3p	/
28549913	qRT-PCR	Tissues	20	86	hsa-miR-186-5p	hsa-miR-148a-3p
31087138	qRT-PCR	Tissue	23	74	/	hsa-miR-552
28624221	qRT-PCR microarray	Tissue	48	243	hsa-miR-1290	/
21694772	Microarray	Tissue	64	158	hsa-miR-31-5p	hsa-miR-224-5p, hsa-miR-552, hsa-miR-592
30119209	qRT-PCR microarray	Tissue	28	Notfound	/	hsa-miR-92a
27198570	Microarray	Tissue	170	Notfound	hsa-miR-31-5p, hsa-miR-330-3p, hsa-miR-223-3p, hsa-miR-520d-3p, hsa-miR-625-5p	hsa-let-7a-5p, hsa-miR-148a-3p, hsa-miR-221-3p, hsa-miR-361-5p, hsa-miR-92a-3p, hsa-miR-483-3p

Table 4

Pearson’s correlation analysis between miRNAs and mRNAs in TCGA

miRNA(up)	Target gene (down)	$p$ -value	Adjust $p$ -value	$r$ value
hsa-miR-155-5p	CAB39L	0.0124	0.0195	$-$ 0.1281
hsa-miR-223-3p	CFTR	0.0236	0.0327	$-$ 0.1161
hsa-miR-31-5p	PLAGL2	0.0001	0.0002	$-$ 0.2012
	FARP1	0.0308	0.0407	$-$ 0.1050
miRNA (down)	Target gene (up)	$p$ -value	Adjust $p$ -value	$r$ value
hsa-miR-1-3p	RAB27B	0.0297	0.0389	$-$ 0.1115
hsa-miR-148a-3p	TRIB2	0.0342	0.0439	$-$ 0.1087
hsa-miR-196b-5p	TFAP2 A	$<$ 0.0001	0.0001	$-$ 0.2134
hsa-miR-361-5p	RAB27B	0.0001	0.0002	$-$ 0.2058
hsa-miR-92a-3p	DUSP4	0.0001	0.0003	$-$ 0.1972
	PLK2	0.0395	0.0498	$-$ 0.1057

Table 5

The expression level of miRNAs and mRNAs detecting in MSI CRC and MSS CRC samples by RT-qPCR

miRNA (up)	$p$ -value	Fold change (MSI vs. MSS)	Targetgene (down)	$p$ -value	Fold change (MSI vs. MSS)	Pearson’s correlation
						$p$ -value	$r$ value
hsa-miR-155-5p	0.022	1.779	CAB39L	0.003	0.326
hsa-miR-223-3p	0.056	1.073	CFTR	0.386	0.371
hsa-miR-31-5p	0.003	1.876	PLAGL2	0.002	0.074	$<$ 0.001	$-$ 0.354
hsa-miR-31-5p	0.003	1.876	FARP1	0.066	0.391
miRNA (down)	$p$ -value	Fold change (MSI vs. MSS)	Targetgene (up)	$p$ -value	Fold change (MSI vs. MSS)	Pearson’s correlation
						$p$ -value	$r$ value
hsa-miR-1-3p	0.696	0.147	RAB27B	0.0070	1.579
hsa-miR-148a-3p	0.021	1.755	TRIB2	$<$ 0.0001	2.843
hsa-miR-196b-5p	0.425	0.381	TFAP2 A	$<$ 0.0010	2.508
hsa-miR-361-5p	$<$ 0.001	0.359	RAB27B	0.0070	1.579	$<$ 0.001	$-$ 0.3270
hsa-miR-92a-3p	0.210	0.243	DUSP4	0.3610	0.827
			PLK2	$<$ 0.0010	2.718

Figure 3.

GO and KEGG pathway analysis showing the associated function of the target genes of dysregulated miRNAs in colorectal cancer. (A) The 19 enriched KEGG pathways and the combination of GO terms including the top 15 GO biological processes, 2 cellular components and 6 molecular functions for 21 upregulated miRNAs ordered by $-\log_{10}$ $P$ -value in CRC. (B) The top 20 enriched KEGG pathways and the combination of GO terms including the top 20 GO biological processes, 7 cellular components and 6 molecular functions for 17 downregulated miRNAs ordered by $-\log_{10}$ $P$ -value in colorectal cancer. (C) The top 10 enriched GO terms and KEGG pathways for downregulated mRNAs ordered by adjust $p$ -value in MSI CRCs versus MSS CRCs. There is no result of upregulated mRNAs using R-clusterProfiler by Hiplot platform.

3.2 Analysis of function enrichment analysis and pathway analysis

We performed the enrichment analysis of 38 DEMs by DAVID-mirPath. Targets of up-DEMs and down-DEMs enriched in KEGG pathways and GO terms were listed in Table S2 and S3 respectively, and the top 20 ordered by $-\log_{10}$ $P$ -value were shown in Fig. 3A and Fig. 3B respectively. We also analyzed 34 down-DEGs using clusterProfiler via Hiplot platform separately shown in Fig. 3C. There were no results of up-DEGs evaluating in Hiplot. Dysregulated miRNAs were frequently associated with tumor-related signaling pathways, innate immune response, and toll-like receptors signaling pathway.

3.3 CMap analysis of dysregulated genes in CRC

We employed cMap to find potential compounds that can disturb the dysregulated gene expression pattern. After the query for dysregulated genes tags, 35 compounds, of which the significantly normal_ncs and fdr_q_nlog10 in CRC cell lines (HCT116, HT1299, HT29 and SW480) were identified as the potential drugs for CRC shown in Fig. S1 A. CMap mode of action (MoA) for 35 compounds revealed 30 mechanisms of action shared by the above compounds shown in Fig. S1B.

3.4 Screening of negative miRNA/mRNA regulatory pairs associated with MSI CRC

Firstly, the experimentally validated target mRNAs of 38 DEMs were selected by miRtarbase and Tarbase. There were 65 miRNA(up) $-$ mRNA(down) pairs and 74 miRNA(down) $-$ mRNA(up) pairs screened out after intersections of DEGs between 108 DEGs and two databases shown in Fig. 2B. Ten negatively correlated pairs were listed in Table 4. AUC of 131071 alternative models were listed in Table S4.

Figure 4.

Validating the expression of 8 differentially expressed miRNAs and 9 differentially expressed mRNAs by RT-qPCR. (A) Two up-DEMs in MSI CRC were miR-155-5p (FC $=$ 1.779) and miR-31-5p (FC $=$ 1.876), while miR-148a-3p (FC $=$ 1.755) and miR-361-5p (FC $=$ 0.359) were downregulated. (B) Two down-DEGs were CAB39L (FC $=$ 0.326) and PLAGL2 (FC $=$ 0.074), while four up-DEGs were RAB27B (FC $=$ 1.579), TRIB2 (FC $=$ 2.843), TFAP2A (FC $=$ 2.508), and PLK2 (FC $=$ 2.718). (C) Pearson’s correlation analysis of miRNA-mRNA regulatory pairs in 80 tissue samples. Two negatively correlated miRNA-mRNA regulatory pairs were plotted. Data are presented as mean $\pm$ SEM. * $p<$ 0.05, ** $p<$ 0.01 and *** $p<$ 0.001 (Student’s $t$ -test). $P$ -values are listed in Table 5.

3.5 Validation of miRNAs and mRNAs expression in CRC tissues by RT-qPCR

MiR-155-5p and miR-31-5p were upregulated in MSI CRCs, while miR-148a-3p and miR-361-5p were downregulated. CAB39L and PLAGL2 were downregulated, while RAB27B, TRIB2, TFAP2A, and PLK2 were upregulated. By the Pearson’s correlation analysis, we identified that miR-31-5p and miR-361-5p were associated with PLAGL2 and RAB27B respectively shown in Fig. 4C and Table 5.

3.6 Evaluation of predictive value of miRNA-mRNA regulatory pairs in CRC

The logistic regression analysis was used to evaluated the predictive value of a panel including 2 miRNA-mRNA regulatory pairs: miR-31-5p/ PLAGL2 and miR-361-5p/ RAB27B in TCGA and validated in 80 samples. ROC curves were plotted in Fig. 5A and D. We also used DCA to verify the clinical applicability of 15 alternative models. They all showed that the complex model based on 4 signatures (the red line shown in Fig. 5C and F) could predict the CRC much better than other 14 models. Clinical impact curves were plotted to demonstrate the significance of the complex predictive model shown in Fig. 5B and E.

3.7 Association analysis of clinical pathological features and miRNA/mRNA expression level in CRC

Clinical pathological data of CRCs were summarized in Table 1. Because of no enough clinical information, we evaluated miRNAs and mRNAs expression levels in multiple subgroups in TCGA shown in Fig. S5A–C. ROC curve was used to find the best cutoff which was as the basis for grouping of expression levels of miRNAs and mRNAs shown in Fig. S2. MiR-31-5p and RAB27B were associated with gender, and miR-361-5p and PLAGL2 were dysregulated in different sites. As shown in Fig. S5C, we found that the expression level of miR-30-5p was associated with the stage of MSI CRCs and MSS CRCs.

Table 6
Cox regression analysis for the prediction of overall survival in MSS CRCs

	Univariate cox regression		Multiple cox regression
	Hazard ratio (95% confidence interval)	$P$ value	Hazard ratio (95% confidence interval)	$P$ value
TCGA-CRC ( $n=$ 273)
Gender	2.463 (1.290–4.442)	0.018	2.068 (1.098–4.536)	0.028
Age	1.718 (1.120–2.911)	0.025	1.806 (1.115–2.923)	0.025
Tumor stage	1.400 (0.541–3.905)	0.520
Tumor site	3.859 (2.122–7.019)	0.839
KRAS	1.130 (0.576–2.216)	0.722
BRAF	0.458 (0.188–1.112)	0.084
NRAS	0.452 (0.105–1.940)	0.285
HRAS	0.201 (0.047–0.872)	0.017	0.192 (0.043–0.863)	0.014
RAS	0.889 (0.457–1.731)	0.729
RAB27B	1.156 (0.557–2.400)	0.697
Hsa-miR-31-5p	1.632 (1.052–2.532)	0.029	1.604 (1.029–2.501)	0.037
Hsa-miR-361-5p	1.086 (0.558–2.113)	0.809

MSS: microsatellite stable; CRC: colorectal cancer; HR: hazard ration; CI: confidence interval.

Figure 5.

Receiver operating characteristic (ROC) curves, decision curve analysis (DCA) and Clinical impact curves of the complex predictive model including 4 signatures (miR-31-5p, PLAGL2, miR-361-5p and RAB27B) to distinguish MSI CRCs from MSS CRCs. ROC curves of the complex predictive model in TCGA training cohort (A) and external validation cohort (D). The X axis shows the high-risk threshold, and the Y axis represents the standardized net benefit. Clinical impact curves of the complex model including the above mentioned 4 signatures for predicting MSI status in TCGA training cohort (B) and external validation cohort (E). The number of high-risk patients (blue dotted line) and the number of high-risk patients with events (red solid line) are plotted. The decision curve analysis of 15 models, of which a complex predictive model containing all 4 signatures (miR-31-5p, PLAGL2, miR-361-5p and RAB27B) (red solid line), 4 models containing 3 signatures (green dotted lines), 6 models containing 2 signatures (magenta dotted lines) and 4 models containing 1 signature (blue dotted lines) in TCGA training cohort (C) and external validation cohort (F).

We further studied the expression level of miR-31-5p and miR-361-5p in MSI versus MSS in 32 cancer types based on TCGA data that miR-31-5p was found upregulated in MSI in 6 cancer types (fold change $>$ 2), while miR-361-5p was downregulated in 5 cancer types (fold change $<$ 0.6) shown in Fig. 3A.

And then we found BRAF gene mutation and RAS gene family (HRAS, KRAS and NRAS) mutation were associated with the microsatellite instability listed in Table S5. The level of miR-31-5p was found to be higher in BRAF-mutated CRCs, and its target gene PLAGL2 was downregulated in BRAF-mutated CRCs shown in Fig. S6A. MiR-31-5p by targeting PLAGL2 and miR-361-5p by targeting RAB27B also played an important role in RAS gene family mutation shown in Fig. S6B.

3.8 Analysis of overall survival

Because of no enough OS data in MSI TCGA-CRC data, the HR of different clinical features only in MSS CRC was estimated by univariate and multivariate cox regression analysis. Gender, age, and the high expression of miR-31-5p could be as independent risk factors for MSS CRC listed in Table 6. KM curves showed that the miR-31-5p, miR-361-5p, and RAB27B have impact on the OS of MSS CRCs shown in Fig. S4B.

Figure 6.

The profiles of immune-related subtype distribution pattern in MSI CRCs ( $n=$ 55) and MSS CRCs ( $n=$ 272) in TCGA cohort. (A) The bar plot visualizing 22 immune cell subset proportions in each sample. Data are presented as mean. (B) The bar plot of 22 immune cells differentially infiltrated fraction in MSI CRCs versus MSS CRCs. (C) The Kaplan–Meier curves of 4 immune cells which affect the overall survival of CRCs. (D) Correlation heatmap displaying 4 differentially infiltrated immune cells which were correlated with two miRNA-mRNA regulatory pairs. Association between expression of miRNAs/target mRNAs and the proportions of 22 subsets of immune cells by Pearson’s correlation. Association between expression of miRNAs/target mRNAs and four tumor microenvironment factors (stromal score, immune score, ESTIMATE score and tumor purity) (E, G), overall methylation (H), and tumor mutation burden (I). (F) Correlation heatmap displaying 3 significantly different tumor microenvironment factors (immune score, ESTIMATE score and tumor purity) which were correlated with two miRNA-mRNA regulatory pairs.

3.9 Analysis of tumor-related phenotypes associated with signatures

We applied an established computational method (CIBERSORT) to bulk gene expression profiles of CRC to infer the proportions of 22 subsets of immune cells shown in Fig. 6A and listed in Table S6. We found that 4 of 8 types which were differentially expressed in CRCs versus NCs were associated with OS of CRC. We further investigated associations between each cell type and miRNA/ target mRNA expression. The levels of miR-31-5p and target gene PLAGL2, and miR-361-5p and target gene RAB27B were significantly correlated with CD8 ${}^{+}$ T cells. In our study, two pairs could interact with DNA methylation, tumor immunity and inflammation in different groups of CRCs shown in Fig. 6E and G–I. Two pairs were associated with tumor mutation burden shown in Fig. S7.

4. Discussion

We reported that MSI CRCs and MSS CRCs were different in many molecular aspects, such as the miRNA/mRNA expression profiles, overall methylation, BRAF and RAS gene family mutations in our study. Indeed, the differentially expressed miRNAs and mRNAs were associated with the above-mentioned aspects after our analysis.

Function analysis and cMap analysis were conducted that candidate DEMs and DEGs which were screened from TCGA, GEO, and results in published articles were associated with classic cancer-related signaling pathways. Notably, fatty acid metabolism which is the enriched KEGG pathways has been paid enough attention to so far. Fatty acid synthase (FASN) overexpression is associated with MSI-H, independent of sex, tumor location, p53, or KRAS/BRAF status [22]. Then, the connectivity map analysis guided us to 35 potential drugs such as MEK inhibitor. The result also provided evidence that 108 dysregulated genes could represent the molecular background of MSI or MSS.

In this study, we adopted the strict criterion to identify vital miRNA-mRNA regulatory pairs. Briefly, the target mRNA of DEMs should be differentially expressed, significant negative correlated with DEM regulator by Pearson’s correlation, and further validated by two databases containing experimentally validated miRNA-target interactions, miRTarbase and Tarbase. Then, we detected the expression level of 10 negatively correlated miRNA-mRNA regulatory pairs in colon tissues. Ultimately, two pivotal pairs (miR-31-5p/PLAGL2 and miR-361-5p/RAB27B) were considered for inclusion in the logistic regression model. ROC curves and DCA indicated that the complex model (miR-31-5p $+$ PLAGL2+ miR-361-5p $+$ RAB27B) was the best predictive model when compared with other 14 combinations. The finding of miR-31-5p overexpressing in MSI CRC was consistent with other studies’ results [23, 24, 25, 26, 27] and miR-361-5p was also found to be downregulated in MSI CRC by Martha et al. [38]. Many studies have demonstrated the oncogenic properties of miR-31-5p in cancers [28, 29, 30, 31]. It also regulates chemosensitivity in hepatocellular carcinoma [32], malignant pleural mesothelioma [33]and triple-negative breast cancer [34]. MiR-361-5p may affect the epithelial-to-mesenchymal transition (EMT) [35] and inhibit cell metabolism, proliferation, and invasion [36, 37, 38]. Our previous research found cytokine genes such as chemokines (CXCL8, CXCL16) and interferons (IFNA1, IFNB1, IFNG, IFNL1) were overexpressed in human normal bronchial epithelial cells transfected the RNA precursor of miR-361-5p. We speculate that overexpression of miR-361-5p can be associated with the suppression of inflammation. Microsatellite instability may be possible to be correct by increasing the expression level of miR-361-5p, manipulating tumor microenvironment and attenuating inflammation.

We summarized the clinical pathological features in TCGA and miR-31-5p was found to be differentially expressed in different tumor stages of MSI CRC and MSS CRC. RAS gene family (HRAS, KRAS and NRAS) mutation status and BRAF mutation status were confirmed to interact with the expression level of miR-31-5p, miR-361-5p and their target genes. Interestingly, the high expression level of miR-31-5p might be a risk factor in MSS CRC patients and it could be a prognostic biomarker of MSS CRC. MiR-31-5p by targeting PLAGL2 and miR-361-5p by targeting RAB27B could have an impact on tumor microenvironment by interacting with tumor-related inflammation and infiltration of macrophages and CD8 ${}^{+}$ T cells. CD8 ${}^{+}$ T cells density is found to be correlated positively with the frameshift mutations in gene sequences [39, 40]. It supports the view that MSI-H CRCs may be more immunogenic than MSS CRCs and therefore MSI-H CRC patients show a better prognosis [41]. Especially, RAB27B could lead to high tumor mutation burden.

In our study, the differences in miRNA and targeted mRNA expression between MSI CRCs versus MSS CRCs were reported and bio-informatic analysis on these significantly differentially expressed molecules had revealed the potential association with cell cycle, cell proliferation, cell motility, and innate immune response especially toll-like receptors signaling pathway based on function annotation and pathway enrichment analysis. Phenotype analysis of tumor microenvironment has confirmed that two miRNA-mRNA regulatory pairs which could distinguish MSI CRCs versus MSS CRCs might play an important role in tumor-related immunity. Additional studies on the molecular and biological functions of these differentially expressed miRNAs and genes, and the possibility of therapeutic targets for regulating the tumor microenvironment to enhance immunogenicity that could improve the sensitivity of MSI CRCs to immunotherapy will be required to substantiate this hypothesis.

5. Conclusion

This study established the miRNA-mRNA regulatory network and generated a predictive model of the combination of miRNAs and mRNAs to distinguish MSI CRCs versus MSS CRCs. Real-time detecting the expression levels of miRNAs and target mRNAs during the development of the colorectal cancer may expand the population of CRC who may be earlier access to benefit from immune checkpoint inhibitor-based immunotherapy.

Author contributions

Conception and design: Wei Zhu, Wenfang Cheng and Tongshan Wang.

Administrative support: Wei Zhu.

Provision of study materials or patients: Wei Zhu and Guoxin Song.

Collection and assembly of data: Cheng Liu and Xuan Zou.

Data analysis and interpretation: Cheng Liu, Xingchen Fan, Shuang Peng, and Shiyu Zhang.

Manuscript writing: All authors.

Final approval of manuscript: All authors.

Supplementary data

The supplementary files are available to download from http://dx.doi.org/10.3233/CBM-210408.

Footnotes

Acknowledgments

We thank openbiox community and Hiplot team (https://hiplot.com.cn) for providing technical assistance and valuable tools for data analysis and visualization.

References

Sung

Ferlay

Siegel

R.L.

Laversanne

Soerjomataram

Jemal

, et al., Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries, CA Cancer J Clin 71(3) (2021), 209–249.

Gendoo

D.M.

and Haibe-Kains

, MM2S: Personalized diagnosis of medulloblastoma patients and model systems, Source Code Biol Med 11 (2016), 6.

Guinney

Dienstmann

Wang

, A. de Reynis Schlicker

Soneson

, et al., The consensus molecular subtypes of colorectal cancer, Nat Med 21(11) (2015), 1350–1356.

Luo

Huang

Luo

and Zhu

, Microsatellite instability: A review of what the oncologist should know, Cancer Cell Int 20 (2020), 16.

Baudrin

L.G.

Deleuze

J.F.

and How-Kit

, Molecular and Computational Methods for the Detection of Microsatellite Instability in Cancer, Front Oncol 8 (2018), 621.

Greenson

J.K.

Bonner

J.D.

Ben-Yzhak

Cohen

H.I.

Miselevich

Resnick

M.B.

, et al., Phenotype of microsatellite unstable colorectal carcinomas: Well-differentiated and focally mucinous tumors and the absence of dirty necrosis correlate with microsatellite instability, Am J Surg Pathol 27(5) (2003), 563–570.

Lopez-Serra

and Esteller

, DNA methylation-associated silencing of tumor-suppressor microRNAs in cancer, Oncogene 31(13) (2012), 1609–1622.

Tang

Wang

Chen

, et al., Identifying microRNA-mRNA regulatory network in colorectal cancer by a combination of expression profile and bioinformatics analysis, BMC Syst Biol 6 (2012), 68.

Zhang

Chen

Zhang

and Zhang

, Construction of an miRNA-mRNA regulatory network in colorectal cancer with bioinformatics methods, Anti-cancer drugs 30(6) (2019), 588–595.

10.

J Jiang

Xie

Liu

, et al., Establishment of a Colorectal Cancer-Related MicroRNA-mRNA Regulatory Network by Microarray and Bioinformatics, Frontiers in genetics 11 (2020), 560186.

11.

Yoshihara

Shahmoradgoli

Martnez

Vegesna

Kim

Torres-Garcia

, et al., Inferring tumour purity and stromal and immune cell admixture from expression data, Nat Commun 4 (2013), 2612.

12.

Corts-Ciriano

Lee

J.J.

Jain

Jung

Y.L.

Yang

, et al., PCAWG Structural Variation Working Group, Park PJ and PCAWG Consortium, Comprehensive analysis of chromothripsis in 2,658 human cancers using whole-genome sequencing, Nature genetics 52(3) (2020), 331–341.

13.

Karagkouni

Paraskevopoulou

M.D.

Chatzopoulos

Vlachos

I.S.

Tastsoglou

Kanellos

, et al., DIANA-TarBase v8: A decade-long collection of experimentally supported miRNA-gene interactions, Nucleic Acids Res 46(D1) (2018), D239–D245.

14.

Huang

H.Y.

Lin

Y.C.

Huang

K.Y.

Shrestha

Hong

H.C.

, et al., miRTarBase 2020: Updates to the experimentally validated microRNA-target interaction database, Nucleic Acids Res 48(D1) (2020), D148–D154.

15.

Shannon

Markiel

Ozier

Baliga

N.S.

Wang

J.T.

Ramage

, et al., Cytoscape: A software environment for integrated models of biomolecular interaction networks, Genome Res 13(11) (2003), 2498–2504.

16.

Hiplot: A free and comprehensive web platform for scientific data visualization, Openbiox Community (2021), Unpublished.

17.

Wang

Hang

Zou

and Mao

J.H.

, A novel gene expression-based prognostic scoring system to predict survival in gastric cancer, Oncotarget 7(34) (2016), 55343–55351.

18.

Newman

A.M.

Liu

C.L.

Green

M.R.

Gentles

A.J.

Feng

, et al., Robust enumeration of cell subsets from tissue expression profiles, Nat Methods 12(5) (2015), 453–457.

19.

Yoshihara

Shahmoradgoli

Martnez

Vegesna

Kim

Torres-Garcia

, et al., Inferring tumour purity and stromal and immune cell admixture from expression data, Nat Commun 4 (2013), 2612.

20.

Mayakonda

Lin

D.C.

Assenov

Plass

and Koeffler

H.P.

, Maftools: Efficient and comprehensive analysis of somatic variants in cancer, Genome Res 28(11) (2018), 1747–1756.

21.

Livak

K.J.

and Schmittgen

T.D.

, Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method, Methods 25(4) (2001), 402–408.

22.

Ogino

Kawasaki

Ogawa

Kirkner

G.J.

Loda

and Fuchs

C.S.

, Fatty acid synthase overexpression in colorectal cancer is associated with microsatellite instability, independent of CpG island methylator phenotype, Hum Pathol 38(6) (2007), 842–849.

23.

Kubota

Taniguchi

Nyuya

Umeda

Mori

Fujiwara

, et al., Upregulation of microRNA-31 is associated with poor prognosis in patients with advanced colorectal cancer, Oncol Lett 19(4) (2020), 2685–2694.

24.

Earle

J.S.

Luthra

Romans

Abraham

Ensor

Yao

, et al., Association of microRNA expression with microsatellite instability status in colorectal adenocarcinoma, J Mol Diagn 12(4) (2010), 433–440.

25.

Sarver

A.L.

French

A.J.

Borralho

P.M.

Thayanithy

Oberg

A.L.

Silverstein

K.A.

, et al., Human colon cancer profiles show differential microRNA expression depending on mismatch repair status and are characteristic of undifferentiated proliferative states, BMC Cancer 9 (2009), 401.

26.

Oberg

A.L.

French

A.J.

Sarver

A.L.

Subramanian

Morlan

B.W.

Riska

S.M.

, et al., MiRNA expression in colon polyps provides evidence for a multihit model of colon cancer, PLoS One 6(6) (2011), e20465.

27.

Slattery

M.L.

Herrick

J.S.

Mullany

L.E.

Wolff

Hoffman

M.D.

Pellatt

D.F.

, et al., Colorectal tumor molecular phenotype and miRNA: Expression profiles and prognosis, Mod Pathol 29(8) (2016), 915–927.

28.

Liu

and Sai

, High miR-31-5p expression promotes colon adenocarcinoma progression by targeting TNS1, Aging (Albany NY) 12(8) (2020), 7480–7490.

29.

Zhao

Han

Zhang

Wang

and Xu

, Increased expression of microRNA-31-5p inhibits cell proliferation, migration, and invasion via regulating Sp1 transcription factor in HepG2 hepatocellular carcinoma cell line, Biochem Biophys Res Commun 490(2) (2017), 371–377.

30.

Liu

Cai

Dong

Jiang

Yang

, et al., Circ-BPTF promotes bladder cancer progression and recurrence through the miR-31-5p/RAB27A axis published correction appears in Aging (Albany NY), Aging 13(2) (2021), 3162–3164.

31.

Y.L.

Liang

Shi

G.Q.

Cao

Qiu

Yuan

, et al., LINC00689 participates in proliferation, chemoresistance and metastasis via miR-31-5p/YAP/β-catenin axis in colorectal cancer, Exp Cell Res 395(1) (2020), 112176.

32.

Que

K.T.

Zhou

You

Zhang

Zhao

X.P.

Gong

J.P.

, et al., MicroRNA-31-5p regulates chemosensitivity by preventing the nuclear location of PARP1 in hepatocellular carcinoma, J Exp Clin Cancer Res 37(1) (2018), 268.

33.

Moody

H.L.

Lind

M.J.

and Maher

S.G.

, MicroRNA-31 Regulates Chemosensitivity in Malignant Pleural Mesothelioma, Mol Ther Nucleic Acids 8 (2017), 317–329.

34.

Shen

Lei

and Du

, miR-31-5p may enhance the efficacy of chemotherapy with Taxol and cisplatin in TNBC, Exp Ther Med 19(1) (2020), 375–383.

35.

Zhang

Wei

Liu

and Qu

, MicroRNA-361-5p inhibits epithelial-to-mesenchymal transition of glioma cells through targeting Twist1, Oncol Rep 37(3) (2017), 1849–1856.

36.

Sun

J.J.

Chen

G.Y.

and Xie

Z.T.

, MicroRNA-361-5p Inhibits Cancer Cell Growth by Targeting CXCR6 in Hepatocellular Carcinoma, Cell Physiol Biochem 38(2) (2016), 777–785.

37.

Zhang

and Guo

, MiR-361-5p inhibits glycolytic metabolism, proliferation and invasion of breast cancer by targeting FGFR1 and MMP-1, J Exp Clin Cancer Res 36(1) (2017), 158.

38.

Song

Guo

Zhang

Zheng

Lin

, et al., MiR-361-5p inhibits colorectal and gastric cancer growth and metastasis by targeting staphylococcal nuclease domain containing-1, Oncotarget 6(19) (2015), 17404–17416.

39.

Tougeron

Fauquembergue

Rouquette

, F. Le Pessot Sesbo

Laurent

, et al., Tumor-infiltrating lymphocytes in colorectal cancers with microsatellite instability are correlated with the number and spectrum of frameshift mutations, Mod Pathol 22(9) (2009), 1186–1195.

40.

Maby

Tougeron

Hamieh

Mlecnik

Kora

Bindea

, et al., Correlation between Density of CD8+T-cell Infiltrate in Microsatellite Unstable Colorectal Cancers and Frameshift Mutations: A Rationale for Personalized Immunotherapy, Cancer Res 75(17) (2015), 3446–3455.

41.

Phillips

S.M.

Banerjea

Feakins

S.R.

Bustin

S.A.

and Dorudi

, Tumour-infiltrating lymphocytes in colorectal cancer with microsatellite instability are activated and cytotoxic, Br J Surg 91(4) (2004), 469–475.

Comprehensive analysis of negatively correlated miRNA-mRNA regulatory pairs associated with microsatellite instability in colorectal cancer

Abstract

BACKGROUND:

OBJECTIVE:

METHODS:

RESULTS:

CONCLUSION:

Keywords

2. Materials and methods

2.1 Data acquisition and processing

2.2 Function annotation and pathway enrichment analysis of miRNA/target-gene pairs

2.3 Connectivity map analysis of potential compounds capable of targeting the differentially expressed genes

2.4 Survival analysis

2.5 Evaluation of interactions of miRNA-mRNA pairs and tumor-relative phenotypes and gene mutation status

2.6 Sample collection and RNA isolation

Table 1 Clinicopathological and molecular features of colorectal cancer patients

2.8 Statistical analysis

3. Results

3.1 Identification of differentially expressed miRNAs (DEMs) and genes (DEGs) in colorectal cancer

3.3 CMap analysis of dysregulated genes in CRC

3.4 Screening of negative miRNA/mRNA regulatory pairs associated with MSI CRC

3.6 Evaluation of predictive value of miRNA-mRNA regulatory pairs in CRC

3.7 Association analysis of clinical pathological features and miRNA/mRNA expression level in CRC

Table 6 Cox regression analysis for the prediction of overall survival in MSS CRCs

4. Discussion

5. Conclusion

Author contributions

Supplementary data

Footnotes

Acknowledgments

References

Table 1
Clinicopathological and molecular features of colorectal cancer patients

Table 6
Cox regression analysis for the prediction of overall survival in MSS CRCs