Analysis of aberrant miRNA-mRNA interaction networks in prostate cancer to conjecture its molecular mechanisms

Abstract

BACKGROUND:

MicroRNAs (miRNAs) capable of post-transcriptionally regulating mRNA expression are essential to tumor occurrence and progression.

OBJECTIVE:

This study aims to find negatively regulatory miRNA-mRNA pairs in prostate adenocarcinoma (PRAD).

METHODS:

Combining The Cancer Genome Atlas (TCGA) RNA-Seq data with Gene Expression Omnibus (GEO) mRNA/miRNA expression profiles, differently expressed miRNA/mRNA (DE-miRNAs/DE-mRNAs) were identified. MiRNA-mRNA pairs were screened by miRTarBase and TarBase, databases collecting experimentally confirmed miRNA-mRNA pairs, and verified in 30 paired prostate specimens by real-time reverse transcription polymerase chain reaction (RT-qPCR). The diagnostic values of miRNA-mRNA pairs were measured by receiver operation characteristic (ROC) curve and Decision Curve Analysis (DCA). DAVID-mirPath database and Connectivity Map were employed in GO/KEGG analysis and compounds research. Interactions between miRNA-mRNA pairs and phenotypic features were analyzed with correlation heatmap in hiplot.

RESULTS:

Based on TCGA RNA-Seq data, 22 miRNA and 14 mRNA GEO datasets, 67 (20 down and 47 up) miRNAs and 351 (139 up and 212 down) mRNAs were selected. After screening from 2 databases, 8 miRNA (up)-mRNA (down) and 7 miRNA (down)-mRNA (up) pairs were identified with Pearson’s correlation in TCGA. By external validation, miR-221-3p (down)/GALNT3 (up) and miR-20a-5p (up)/FRMD6 (down) were chosen. The model combing 4 signatures possessed better diagnostic value. These two miRNA-mRNA pairs were significantly connected with immune cells fraction and tumor immune microenvironment.

CONCLUSIONS:

The diagnostic model containing 2 negatively regulatory miRNA-mRNA pairs was established to distinguish PRADs from normal controls.

Keywords

miRNA miRNA-mRNA networks prostate adenocarcinoma TCGA PCR

1. Introduction

According to the global cancer statistics estimated by the World Health Organization (WHO) for male in 2020, the incidence and mortality of prostate adenocarcinoma (PRAD) ranked second and fifth respectively [1]. In the past few decades, extensive screening based on prostate-specific antigen (PSA) has largely improved early detection of prostate cancer, thus decreasing prostate cancer mortality [2]. However, PSA-based screening was not recommended in men over 70 years by US Preventive Services Task Force because of higher risk of overdiagnosis and overtreatment, which may bring about urinary incontinence and erectile dysfunction [3]. Since then, finding biomarkers with better diagnostic efficacy and prognostic evaluation is of great significance.

MicroRNAs (MiRNAs) contain approximately 22-nucleotide-long small noncoding RNAs, which play an crucial role in GENE regulatory networks by forming imperfect binding with the message RNA(mRNA) 3’-untranslated region (3’-UTR) [4]. MiRNA-mRNA interactions are widely involved in the occurrence and development of various tumors. For instance, miR-125a-3p targeting at RHOA promotes epithelial-mesenchymal transition in pancreatic cancer [5]. Mi RNA-mRNA interaction pairs, such as miR-135a-3p targeting at HOXD10, miR-371a-3p targeting at CACNA1D together with miR-106a-5p targeting at miR-106a-5p, participated in the complex molecules of prostate cancer [6]. Dysregulation of miR-193a-5p/TP73 pair and miR-188-5p/UBE2I pair also indicate tumor progression and poor prognosis in PRAD patients [7]. Thus, defining more reliable and effective biomarkers such as miRNA-mRNA pairs is promising in diagnosis and prognosis of prostate cancer.

Obtaining expression patterns from high throughput sequencing and microarrays is now recognized as a widely-applied and mature technology in acquiring more global views on miRNA or genes expression in tumor [8]. A total of 67 (20 down & 47 up) miRNAs and 351 (139 up & 212 down) mRNAs were selected by analyzing The Cancer Genome Atlas (TCGA) RNA-Seq data, 22 miRNA and 14 mRNA Gene Expression Omnibus (GEO) datasets. After screening from miRTarBase and TarBase, 8 miRNA (up)-mRNA (down) and 7 miRNA (down)-mRNA (up) pairs were identified with Pearson’s correlation analysis, the expression of which were detected in 30 paired prostate specimens by reverse transcription quantitative real-time PCR (RT-qPCR). Via performing Wilcoxon test, miR-221-3p (down)/GALNT3 (up) and miR-20a-5p (up)/FRMD6 (down) were chosen. Clinical analysis and survival analysis were performed to evaluate whether these miRNA-mRNA pairs could be indicators of diagnosis and prognosis in prostate cancer. Finally, correlation analysis was carried out between TCGA clinical phenotype data composed of immune microenvironment, tumor infiltrating immune cells, global methylation and tumor mutational burden and the two pairs of miRNA-mRNA expression data.

In the past decades, a variety of biomarkers have been employed in the diagnosis and prognosis prediction of prostate cancer. Although some researchers have performed differential expression analysis of miRNAs and mRNAs in prostate cancer, this study is the first to include 36 GEO datasets and TCGA databases to define the differential expression profiles. Our researches based on massive gene expression profiles and bioinformatic analysis and RT-qPCR may shed new light on the detailed mechanism in PRAD pathogenesis. Dysregulated miRNA-mRNA pairs may become underlying prognostic biomarkers and novel therapy target.

2. Materials and method

2.1 Data collection and processing

Figure 1.

The flow chart showed the processes for identifying the negatively miRNA-mRNA regulatory pairs and conducting extensive analysis of pairs in prostate cancer.

Totally 22 gene expression profiles and 14 miRNA expression profiles were downloaded from the Gene Expression Omnibus database (GEO, http://www.ncbi. nlm.nih.gov/geo/) with the key word “(prostate OR prostatic) AND (cancer OR tumor OR carcinoma OR adenocarcinoma) NOT (cell line)”. We chose “Homo sapiens” in “organism” and “expression profiling by array” & “Non-coding RNA profiling by array” in “study type”. TCGA RNA-sequencing profiles together with relevant clinical information of PRAD were obtained from the GDC data portal in the National Cancer Institute (https://portal.gdc.cancer.gov/). As is displayed in Fig. 1, dysregulated miRNAs of PRAD mentioned in Human MicroRNA Disease Database (HMDD) [9], Database of Differentially Expressed MiRNAs in human Cancers (dbDEMC) [10] and miRCancer [11] were also taken into consideration. The microarray expression profiles possessed from GEO datasets were analyzed via GEO2R, an online appliance for comparison by the GEOquery and limma R packages (http://www.ncbi.nlm.nih.gov/geo/geo2r/). The RNA-sequencing data derived from TCGA were processed with R software with psych [12], a data analysis package based on R language.

2.2 Gene oncology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) pathway enrichment analysis and connectivity map analysis

In this study, Gene oncology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses of DE-miRNAs and DE-mRNAs were conducted by DAVID-mirPath [13] and the cluster Profiler tool on Hiplot (https://hiplot.com.cn), respectively for miRNAs and mRNAs. The DE-mRNAs lists made up of dysregulated tags were uploaded to the cMap [14] (https://clue.io/query), which collected gene expression patterns of cells treated with about 5,000 small molecule compounds. Normalized connectivity source (norm_cs) was employed to measure the relevance between dysregulated mRNAs and agents and the cut-off criterion was the negative normalized connectivity score (norm_cs) $<$ $-$ 1.35 and the significant negative log10Q-value (fdr_q_nlog10) $>$ 2.

2.3 Forecasting of miRNA-mRNA interaction

We combined miRTarBase (http://miRTarBase.cuhk .edu.cn/) with TarBase (http://www.microrna.gr/tarbase), both of which are databases collecting experimentally verified miRNA-mRNA networks, to predict the target genes of DE-miRNAs. Next, the negative correlation between miRNA and mRNA expression level was assessed by Pearson’s correlation using psych R-based package. The miRNA-mRNA interactions were plotted with Cytoscape v3.8.2, an opensource software for complex network.

2.4 Sample collection

30 paired Formalin-Fixed and Paraffin-Embedded (FFPE) samples of PRAD tissues and adjacent normal prostate tissues were taken from patients undergoing surgery at the First Affiliated Hospital with Nanjing Medical University, the clinical features of which were presented in Table 1. All patients did not receive any chemotherapy or radiotherapy before their surgery. The study was conducted in accordance with the guidelines of the Hospital Ethics Committee and approved by the Institutional Review Boards of the First Affiliated Hospital of Nanjing Medical University (ID: 2016-SRFA-148). This study was conducted in accordance with the Declaration of Helsinki. All patients were informed about the study and invited to participate. All patients signed an informed consent statement before starting the study.

Table 1
Clinicopathological and molecular features of prostate cancer patients

	Prostate cancer ( $n=$ 30)	Rate (%)
Age(year)
Mean (SD)	66.7 (7.07)
Median (min, max)	66 (56, 82)
Cancer location
Left prostate	6	20
Right prostate	3	10
Both sides	21	70
Gleason score
6	5	16.67
7	14	47
8	2	6.7
9	4	13.3
NA	5	16.7
WHO/ISUP grade
1	5	16.67
2	6	20.0
3	8	26.67
4	2	6.7
5	4	13.3
NA	5	16.7
Nearby lymph node
Yes	6	20
No	24	80

2.5 RNA isolation and quantitative reverse transcription qPCR

FFPE tissues were stored at 20 ${}^{\circ}$ C. RNAprep Pure FFPE Kit (TIANGEN Biotech, Beijing, China, DP439) was applied in isolating total RNA from FFPE prostate tissue samples according to the instructions. The concentrations of entire RNA were measured via NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific). To verify the DE-miRNAs and DE-mRNAs of bioinformatics analysis, we conducted reverse transcription with PrimeScript RT reagent Kit (Takara, Tokyo, Japan, RR037A) and performed RT-qPCR with SYBR Premix Ex Taq II (Takara, Tokyo, Japan, RR820A). All kits above were used following the manufactures’ interactions. We performed the PCR detection on the qTOWER ${}^{3}$ 84 (Analytik Jena) plate in a volume of 10 $\mu$ L with the temperature setting 95 ${}^{\circ}$ C for 20 s, followed by 40 cycles of 10 s at 95 ${}^{\circ}$ C and 20 s at 60 ${}^{\circ}$ C in the final volumes of 10 $\mu$ L. GAPDH and 18S were set as reference genes of mRNA; RUN6B was recognized as internal reference of miRNA [15, 16]. All the sequences of primer used in this study were acquired from Primer Bank [17] (https://pga.mgh.harvard.edu/primerbank/) and listed in Table S1.

2.6 Estimation of miRNA-mRNA interactions with tumor-related phenotypes

CIBERSORT was adopted to estimate the proportion of 22 immune cell types in tumor samples by conducting deconvolution analysis with gene expression data (https://cibersort.stanford.edu/index.php/) [18]. ESTIMATE, a tool for forecasting tumor purity, was employed to generate the Stromal Score, immune Score and estimate Score of TCGA-PRAD samples [19]. The discrepancy between these immune cells in TCGA-PRAD and normal tissues were compared by the Wilcoxon test via psych package. Correlation heatmap in Hiplot (https://hiplot.com.cn) was applied in evaluating the relevancy between the expression of miRNA-mRNA pairs and tumor phenotypes.

2.7 Diagnostic effectiveness evaluation and overall survival analysis

The gene expression pattern and survival information of PRAD patients downloaded from TCGA database was analyzed to estimate the prediction probability. The Kaplan-Meier analysis and the receiver operation characteristic (ROC) curve were performed by Hiplot (https://hiplot.com.cn) and GraphPad Prism v8.0.1 respectively. The patients with PRAD were divided into two groups (high vs. low) according to the median expression level. We employed rmda (https://cran.r-project.org/web/packages/rmda/index.html) for the Decision Curve Analysis (DCA).

2.8 Statistically analysis

SPSS Statistics 26.0 (IBM, New York, USA), GraphPad prism 8.0.1 and R language v4.0.1 (https://www.r-project.org/) were adopted in data analyses The criterion for DE-mRNAs and DE-miRNAs were $p$ -value $<$ 0.05 and $|\text{log2(FC)}|>$ 0.58. In Pearson’s correlation between miRNA and mRNA $p$ -value $<$ 0.05 was considered significant. The diagnostic value of miRNA-mRNA interactions was assessed by the area under the ROC curves (AUC), where $p$ -value $<$ 0.05 was of statistical significance. DCA based on Binary Logistic Regression was employed to evaluate clinical diagnosis models. The expressions of miRNA and gene were calculated using the 2 ${}^{-\Delta\Delta\text{CT}}$ method ( $\Delta\text{CT}=\text{Ct}_{\text{miRNA or mRNA}}-\text{Ct}_{\text{normalizer}}$ ; Ct: the threshold cycle) [20]. The inter-group significance was examined by Wilcoxon test All derived data were shown as the median with 95% CI. The $p$ -value $<$ 0.05 denotes statistical significance. Hiplot, Cytoscape, GraphPad prism 8.0.1 together with R language v4.0.1 were applied in plotting.

3. Results

3.1 Identification of DE-mRNAs and DE-miRNAs

Figure 2.

(A) The circular bar-plot exhibited the detailed information of 36 GEO datasets. Among the 22 datasets using mRNA microarray, 20 were drawn from tissues, 1 from cultured tumors cells and normal tissues (GSE34312) and 1 from circulating tumor cells (GSE153514) and the remaining 14 miRNA expression profiles included 2 datasets from urine, 7 datasets from tissue, 4 datasets from serum and 1 dataset from plasm. (B) The negative regulatory relationship between 67 microRNAs (miRNAs) and 351 mRNAs were plotted. 118 miRNA (up)-mRNA (down) pairs and 13 miRNA (down)-mRNA (up) pairs were selected via TarBase and MiRTarBase. Orange nodes represent the upregulated miRNAs/mRNAs in PRADs versus NCs, while blue nodes represent the downregulated miRNAs/mRNAs in PRADs versus NCs.

Table 2

Information pertaining to the selected GEO datasets for prostate cancer

	Specimen type	GEO accession	Number of tumor	Number of non-tumor	Total number	Country	Platform
mRNA array	Tissue	GSE133906	4	4	8	USA	GPL10332
		GSE103512	50	7	57	USA	GPL13158
		GSE62872	264	160	424	USA	GPL19370
		GSE73397	3	3	6	China	GPL20952
		GSE140927	4	4	8	China	GPL22120
		GSE89194	49	49	98	USA	GPL22571
		GSE28615	2	1	3	Germany	GPL4133
		GSE46602	36	14	50	Denmark	GPL570
		GSE55945	13	8	21	USA	GPL570
		GSE45016	10	1	11	Japan	GPL570
		GSE32448	40	40	80	USA	GPL570
		GSE26910	12	12	24	Italy	GPL570
		GSE69223	15	15	30	Germany	GPL570
		GSE6956	69	20	89	USA	GPL571
		GSE54809	19	11	30	USA	GPL6244
		GSE86260	7	7	14	Australia	GPL6244
		GSE28204	4	3	7	China	GPL6480
		GSE32571	60	38	98	Germany	GPL6947
		GSE6919	90	81	171	USA	GPL8300
		GSE68882	23	3	26	USA	GPL91
	Cells cultured from tumors	GSE34312	10	10	20	USA	GPL6884
	and normal tissues
	Circulating tumor cells	GSE153514	9	6	15	Spain	GPL10332
miRNA array	Urine	GSE138740	146	89	235	USA	GPL21572
		GSE159177	277	188	465	USA	GPL21572
	Tissue	GSE45604	50	10	60	Spain	GPL14613
		GSE14857	11	10	21	Germany	GPL6955
		GSE76260	32	32	64	Italy	GPL8179
		GSE64318	27	27	54	USA	GPL8227
		GSE23022	20	20	40	Germany	GPL8786
		GSE36802	21	21	42	USA	GPL8786
		GSE46738	53	4	57	Brazil	GPL8786
	Serum	GSE112264	809	241	1050	Japan	GPL21263
		GSE113486	40	100	140	Japan	GPL212633
		GSE113740	25	10	35	Japan	GPL212633
		GSE16512	6	7	13	USA	GPL8686
	Plasma	GSE113234	60	27	87	Italy	GPL19730

As is exhibited in Fig. 2A, totally 14 miRNA expression microarray GEO datasets were selected out, of which 7 were drawn from tissues, 4 from serum, 2 from urine and 1 from plasm. Among these 14 databases, 108 decreased miRNAs were mentioned in 3 or more datasets and 166 increased miRNAs were mentioned in 4 or more datasets. As for mRNA expression profiling by microarray, 22 GEO databases matched the inclusion criteria, specimen of which are tumor tissues, cells cultured from tumors and circulating tumor cells. After comparing to normal samples, 184 upregulated genes appeared in 7 or more databases and 215 down-regulated mRNAs appeared in 9 or more databases. The detailed information of 36 GEO datasets is exhibited in Table 2. After integrating all the databases including GEO datasets, TCGA, dbDEMC, HMDD and miRCancer, 351 DE-mRNAs (139 up and 212 down) together with 67 DE-miRNAs (47 up and 20 down) were finally determined.

3.2 GO and KEGG pathway enrichment analysis and Identification of candidate small molecular compounds

Outcomes of GO/KEGG analysis based on 47 upregulated miRNAs and 20 downregulated miRNAs were separately revealed in Table S2 and Table S3, ranked with counts as the primary keyword and $p$ -value as the secondary keyword. The top 10 of KEGG pathways, top10 of GO biological process, top 2 of GO molecular function and top 2 of GO cellular component were presented in Fig. S1A and S1B. Besides, we also generated GO/KEGG analysis of 351 dysregulated genes. However, no result was put out for upregulated genes, so only the results of downregulated genes were displayed in Fig. S1C. According to the outcomes of cMap database concentrated on prostate cancer cell lines (LNCAP and PC3), 60 agents were selected out as the potential drugs for PRAD. The 20 compounds with smallest norm_cs were presented in Fig. S2A. cMap target analysis indicated a total of 138 drug-target genes of 60 compounds demonstrated in Fig. S2B.

Figure 3.

A total of 15 negatively miRNA-mRNA regulatory pairs validated by RT-qPCR in 30 paired FFPE samples were presented.

3.3 Identification of miRNA/mRNA reverse regulation pairs in PRAD

Table 3
Pearson’s correlation analysis of miRNA-mRNA pairs in prostate cancers in TCGA

miRNA (down)	mRNA (up)	$p$ -value	$r$ -value
hsa-miR-133b	ERG	0.0282	$-$ 0.0896
hsa-miR-205-5p	ERBB3	0.0009	$-$ 0.1705
hsa-miR-221-3p	GALNT3	0.0051	$-$ 0.1458
	TSPAN13	0.0000	$-$ 0.2696
hsa-miR-222-3p	GALNT3	0.0001	$-$ 0.1927
	TSPAN13	0.0000	$-$ 0.2870
hsa-miR-378a-3p	P4HB	0.0212	$-$ 0.1225
miRNA (up)	mRNA (down)	$p$ -value	$r$ -value
hsa-miR-17-5p	FRMD6	0.0001	$-$ 0.2000
	CAV1	0.0183	$-$ 0.1251
hsa-miR-182-5p	PDK4	0.0245	$-$ 0.1199
hsa-miR-183-5p	FBXO32	0.0197	$-$ 0.1237
	GAS1	0.0254	$-$ 0.1192
hsa-miR-20a-5p	CAV1	0.0114	$-$ 0.1331
	FRMD6	0.0348	$-$ 0.1134
hsa-miR-93-5p	FRMD6	0.0000	$-$ 0.2369

To explore the negative regulatory relationships between miRNA and mRNA, 67 DE-miRNAs (47 up and 20 down) were queried in TarBase and MiRTarBase. As is presented in Fig. 2B, 118 miRNA (up)-mRNA (down) pairs and 13 miRNA (down)-mRNA (up) pairs were screened after intersections of 351 DE-mRNAs and mRNAs targeted by DE-miRNAs. Pearson’s correlation conducted in TCGA database selected out 8 miRNA (up)-mRNA (down) pairs and 7 miRNA (down)-mRNA (up) pairs (shown in Fig. 3) with significant negative correlation listed in Table 3.

3.4 Evaluation of miRNAs and mRNAs expression in FFPE PRAD samples via RT-qPCR

Table 4
Pearson’s correlation analysis of miRNA-mRNA pairs in FFPE prostate cancer samples

miRNA (up)	$p$ -value	Fold change	mRNA (down)	$p$ -value	Fold change	Pearson’s correlation
		( $2^{-\Delta\Delta\text{CT}}$ )			( $2^{-\Delta\Delta\text{CT}}$ )	$p$ -value	$r$ -value
17-5p	0.251	1.332	CAV1	0.006	0.472	0.698	$-$ 0.100
			FRMD6	0.009	0.723	0.439	$-$ 0.157
182-5p	0.091	1.500	PDK4	0.794	1.112	0.876	0.039
183-5p	0.182	2.100	FBXO32	0.0002	0.3423	0.234	$-$ 0.284
			GAS1	0.449	1.979	0.410	$-$ 0.167
20a-5p	0.0002	2.029	CAV1	0.006	0.472	0.399	$-$ 0.189
			FRMD6	0.032	0.723	0.023	$-$ 0.446
93-5p	0.1445	2.307	FRMD6	0.032	0.723	0.410	$-$ 0.166
miRNA (down)			mRNA (up)			$p$ -value	$r$ -value
133b	0.0767	1.470	ERG	0.114	1.668	0.844	$-$ 0.051
205-5p	0.0889	2.180	ERBB3	0.665	1.649	0.942	0.020
221-3p	0.0119	0.523	GALNT3	0.016	2.197	0.019	$-$ 0.507
			TSPAN13	0.379	1.721	0.399	$-$ 0.197
222-3p	0.211	1.752	GALNT3	0.016	2.197	0.399	$-$ 0.175
			TSPAN13	0.379	1.721	0.600	0.120
378a-3p	0.0194	0.417	P4HB	0.0002	3.436	0.909	$-$ 0.027

Figure 4.

The expression of 10 DE-miRNAs and 10 DE-mRNAs were confirmed by ploy(A) tailing RT-PCR. (A) The mRNA expression levels of FRMD6, CAV1 and FBXO32 were downregulated, while the expression levels of GALNT3 and P4HB were upregulated. The miRNA expression levels of miR-205-5p and miR-20a-5p climbed while miR-378a-3p and miR-221-3p declined in PRADs. (B) Pearson’s correlation analysis of miRNA-mRNA regulatory pairs in 30 paired samples. MiR-221-3p/GALNT3 and miR-20a-5p/FRMD6 were the two negatively regulated pairs displayed in figure. Data are presented as mean $\pm$ SEM. ${}^{*}p<$ 0.05, ${}^{**}p<$ 0.01, ${}^{***}p<$ 0.001 and ${}^{****}p<$ 0.0001 (Wilcoxon test). $p$ -values are listed in Table 4.

Figure 5.

Receiver operating characteristic (ROC) curves and decision curve analysis (DCA) were adopted to measure the ability of the complex predictive model including 4 signatures (miR-221-3p, GALNT3, miR-20a-5p and FRMD6) to distinguish PRAD samples from normal samples. (A) The ROC curves of 15 models, including 1 model containing all four signatures, 4 models containing three signatures, 6 models containing two signatures and 4 models containing one signature, were drawn with data of TCGA database. (B) The decision curve analysis (DCA) of 15 models were conducted in TCGA database and validation cohort containing 30 paired PRAD tissues by qRT-PCR.

In this study, we measured the expression of 10 DE-mRNAs and 10 DE-miRNAs in 30 paired PRAD and normal tissue by poly(A) RT-qPCR. As is displayed in Fig. 4A, GALNT3 ( $p=$ 0.16, FC $=$ 2.197) and P4HB ( $p=$ 0.0002, FC $=$ 3.44) were remarkably overexpressed in PRAD tissues and FRMD6 ( $p=$ 0.032, FC $=$ 0.72), CAV1 ( $p=$ 0.006, FC $=$ 0.47) together with FBXO32 ( $p=$ 0.0002, FC $=$ 0.34) were downregulated in PRAD samples. Prominent decline of miR-221-3p ( $p=$ 0.0119, FC $=$ 0.52) and miR-378a-3p ( $p=$ 0.0194, FC $=$ 0.42) can be observed in PRAD samples comparing to normal controls, while expression of miR-20a-5p ( $p=$ 0.0002, FC $=$ 2.03) climbed in PRAD tissues, exhibited in Fig. 4B. As for the Pearson’s correlation analysis, miR-221-3p/GALNT3 ( $p=$ 0.019, $r=$ $-$ 0.507) and miR-20a-5p/FRMD6 ( $p=$ 0.023, $r=$ $-$ 0.446) turned out to be statistically significant, shown in Fig. 4C and 4D respectively. The full outcomes are listed in Table 4. As is presented in Fig. S3, IHC images downloaded from the Human Protein Atlas (HPA) database confirmed the higher level of GALNT3 expression and the lower level of FRMD6 expression for PRAD cells than for normal prostate cells.

3.5 Estimation of miRNA-mRNA pairs’ predictive value in PRAD

In present study, we estimated the prognostic value of a panel containing 2 miRNA-mRNA pairs: miR-221-3p/GALNT3 and miR-20a-5p/FRMD6. Receiver operating characteristic (ROC) curves and decision curve analysis (DCA) of 15 models based on logistic regression were displayed in Fig. 5. The areas under the curve (AUC) of the complex model (miR-221-3p, GALNT3, miR-20a-5p, FRMD6) calculated in TCGA-PRAD were 0.98 (95% CI: 0.9656 to 0.9936, $p<$ 0.001), the highest AUC of the 15 models. All curves were plotted in Fig. 5A. We carried out DCA in TCGA database and validation cohort containing 30 paired PRAD tissues to verify the diagnostic applicability of 15 models mentioned above. The complex model may obtain higher diagnostic efficiency than other 14 models if the threshold probability was between 0–0.70 in TCGA testing cohort, while in external validation cohort, the threshold probability was between 0.33 and 0.5 (Fig. 5B).

3.6 Comparing miRNA-mRNA expression level between clinical subgroups in PRAD

Figure 6.

Expression level of miRNA-mRNA pairs were compared in groups divided according to clinical features of PRAD patients in TCGA database. (A) Expression level of has-miR-221-3p was higher in group “age $\leqslant$ 60” than in group “age $>$ 60”. (B) MiR-221-3p, GALNT3 and FRMD6 were upregulated in high Gleason score subgroup, while miR-20a-5p was downregulated in high Gleason score subgroup.

Figure 7.

The picture showed profiles of association between immune-related phenotype and miRNA-mRNA pairs. (A) The bar plot of 22 immune cells differentially infiltrated fraction in PRADs versus NCs. (B) The association between cell fraction of 9 differently infiltrated immune cells and expression level of miRNA-mRNA pairs by Pearson’s correlation. (C) The relevance between expression level of miRNAs and targeting mRNAs and global methylation and four tumor microenvironment factors.

The number of clinical cases we collected was so few that we had to employ information from TCGA-PRAD to compare expression level of miRNAs and mRNA between several subgroups. In Fig. 6A, miR-221-3p expression was higher in patients with diagnosis age younger than 60 ( $p=$ 0.0412, FC $=$ 1.214). No differences were observed in GALNT3, FRMD6, and miR-20a-5p expression between patients older than 60 and those younger than 60. The expression of GALNT3, FRMD6 and miR-221-3p declined in high Gleason score group ( $p<$ 0.001, FC $=$ 0.745; $p<$ 0.001, FC $=$ 0.596; $p<$ 0.001, FC $=$ 0.726), while miR-20a-5p increased in high Gleason score group ( $p<$ 0.001, FC $=$ 1.243).

3.7 Correlation analysis between tumor-related phenotypes and signatures in PRAD

It is presented in Fig. 7A that nine subsets of immune cells showed differentiation between PRAD samples and control tissues, with entire outcome listed in Table S4. We also investigated the correlation between expression of miRNA-mRNA pairs and the 9 categories of immune cells described above. MiR-20a-5p and target gene FRMD6 were associated with T cell regulatory (Treg), neutrophils, monocytes and resting mast cells (Fig. 7B). Besides, miR-221-3p/GALNT3 pairs could interact with M1 macrophages. As is revealed in Fig. 7C, level of the negatively regulatory pair miR-221-3p/GALNT3 and FRMD6 are relevant with tumor purity, ESTIMATE Score, Immune Score, and Stromal Score in tumor tissue, while there was no correlation between miR-20a-5p and these tumor microenvironment evaluation index.

3.8 Overall survival analysis

Only TCGA data was employed in survival analysis due to lack of sufficient clinical information from our FFPE tissue samples and the GEO datasets. Nevertheless, too few deaths were recorded during follow-up of TCGA database (only 10 in 380 cases), thus no statistical significance estimated by univariate and multivariate cox regression analysis was observed between age, tumor stage, miR-221-3p/GALNT3 or FRMD6 expression level and overall survival (OS). As is disclosed in Fig. S4, miR-20a-5p level was significantly correlated with the OS (HR: 1.000437, 95% CI: 1.000018 to 1.000856, $p=$ 0.041).

4. Discussion

MicroRNAs capable of regulating mRNA expression by incompletely binding to mRNA 3’-UTR are of significance in tumor [21]. Numerous researches have partly defined the microRNA signatures in PRAD diagnosis and prognosis [22]. Despite plenty of studies, the role of miRNA-mRNA interactions in PRAD has not been fully identified. In this study, we identified 351 DE-mRNAs (139 up and 212 down) and 67 DE-miRNAs (47 up and 20 down) by combining 36 GEO datasets, TCGA database and several cancer-related databases Then, 118 miRNA (up)-mRNA (down) pairs and 13 miRNA (down)-mRNA (up) pairs were confirmed after intersections between 351 DE-mRNAs, TarBase and miRTarBase. Next, 15 negatively correlated miRNA-mRNA pairs ( $p$ -value $<$ 0.05) were filtered out with Pearson’s correlation of TCGA data. We further explored the expression level of 15 miRNA-mRNA regulatory pairs in 30 paired prostate tissues by poly(A) qRT-PCR. Via Wilcoxon test and Pearson’s correlation, miR-221-3p/GALNT3 and miR-20a-5p/FRMD6 were chosen.

We carried out GO/KEGG analyses and cMap analysis with candidate DE-miRNAs and DE-mRNAs. Pathway enrichment analysis demonstrated that DE-miRNAs and DE-mRNAs of PRAD were associated with lipid metabolism. As is reported, lipid metabolism is involved in cancer originated from the exocrine glands of the reproductive system like breast cancer and prostate cancer [23, 24]. Notably, low carbohydrate or low ketogenic diets are introduced for dietary strategies of PRAD patients owing to the effect of insulin resistance and adipokine [25]. Besides, 7 compounds (KPT-330, phenytoin, chlorpyrifos, montelukast, ketorolac, mitoxantrone and sorafenib) have been reported to exert negative effects on PRAD development in cell experiments, animal experiments or clinical trials [26, 27, 28, 29, 30, 31, 32].

Finally, miR-20a-5p/FRMD6 and miR-221-3p/GALNT3, two significantly negatively correlated miRNA-mRNA pairs, were incorporated into the logistic regression model. Among the 15 models randomly combining 4 signatures, the complex model (miR-221-3p & GALNT3 & miR-20a-5p & FRMD6) showed the best predictive value. For lack of GEO datasets containing both miRNA and mRNA expression profile, RNA-sequencing data from TCGA was employed in ROC and DCA. Insufficient miRNA and mRNA expression pattern put off the verification of their diagnosis potential. Experiments with large sample size are required. In the future, we will conduct further experiment based on more freshly frozen tissues and FFPE samples to verify the prognostic efficiency of 2 miRNA-mRNA pairs.

MiR-20a-5p is a member of the miR-17-92, a highly conserved cluster reported to be overexpressed in colorectal cancer (CRC) [33], breast cancer (BC) [34] and PRAD [35] in plenty of studies. Likewise, miR-221-3p belonging to the miR-221/222 cluster is relevant with cervical squamous cell carcinoma (CSCC) [36], epithelial ovarian cancer (EOC) [37] and prostate cancer cells [38]. It is possible that miR-221-3p may serve as novel biomarkers in diagnosis and prognosis of pancreatic cancer [39], gastric cancer [40] and diffuse large B cell lymphoma (DLBCL) [41]. As is reported in previous studies, the significantly downregulated gene FRMD6, short for 4.1 Ezrin Radixin Moesin (FERM)-Containing Protein, is capable of serving as a novel tumor suppressor gene and prognostic biomarker candidate in prostate cancer [42]. GALNT3, Psolypeptide N-Acetylgalactosaminyltransferase 3, was found to be efficient predictors for discriminating PRAD patients from normal controls [43].

By analyzing the overall survival (OS) data of TCGA, we observed that high miR-20a-5p expression indicated poor OS in prostate cancer ( $p=$ 0.046 in K-M survival curve), which corresponds with the research of Stoen et al. [44]. Besides, it has been reported that upregulation of FRMD6 could function as tumor suppressor in prostate cancer [42]. Similarly, miR-221 has been suggested as prognostic predictor of cancer-related death [45]. Even if neither FRMD6 nor miR-221-3p showed significant relevancy with OS in Kaplan-Meier (K-M) survival curves, it is unreasonable to completely deny the prognostic value of FRMD6 and miR-221-3p in PRAD. Insufficient clinical information of follow-up study blocked detailed analysis, and the prognostic value of these two miRNA-mRNA pairs need to be further acknowledged.

In addition, miR-20a-5p/FRMD6 and miR-221-3p/GALNT3 may affect tumor immune microenvironment by influencing infiltration of immune cells including regulatory T cells (Treg), neutrophils and mast cells, thus interact with cancer-related inflammation. High level of miR-221-3p was reported to promote the transformation of M2-macrophages to a pro-inflammatory status, leading to immune hyperfunction [46].

MiR-20a-5p also has significant impact on pro-inflammatory macrophages and processes fibrosis in tumor tissue [47]. The specific effects of miRNA-mRNA interaction networks on the tumor immune microenvironment and their molecular mechanisms need to be further studied.

In this study, we only investigated dysregulation of miRNAs and mRNAs in tissues, which is not convenient enough in clinical application. Besides, lack of further research for the detailed mechanisms of the miRNA targeting mRNA level influence cancer cells is also weakness of our study. Thus, studies on body fluid and specific mechanisms deserve further exploration.

5. Conclusions

In summary, the present study tried to figure out the essential pairs in PRAD miRNA-mRNA network via bioinformatics analysis and RT-qPCR. Our research may provide novel biomarkers for diagnosis and prognosis of PRAD in near future.

Author contributions

Conception: Tongshan Wang, Feng Gao and Wei Zhu.

Interpretation or analysis of data: Shuang Peng and Cheng Liu.

Preparation of the manuscript: Shuang Peng, Cheng Liu and Xingchen Fan.

Revision for important intellectual content: Shuang Peng, Jingfeng Zhu and Shiyu Zhang.

Supervision: Xin Zhou.

Supplementary data

The supplementary files are available to download from https://dx-doi-org.web.bisu.edu.cn/10.3233/CBM-220051.

sj-pdf-1-cbm-10.3233_CBM-220051.pdf - Supplemental material

Supplemental material, sj-pdf-1-cbm-10.3233_CBM-220051.pdf

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. This study was supported by the Graduate Research and Practice innovation Plan of Graduate Education Innovation Project in Jiangsu Province [Grant number: JX10213729].

References

Sung

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.

Carlsson

S.V.

and Vickers

A.J.

, Screening for Prostate Cancer, Med Clin North Am 104(6) (2020), 1051–1062.

US Preventive Services Task Force Grossman

D.C.

et al., Screening for Prostate Cancer: US Preventive Services Task Force Recommendation Statement, JAMA 319(18) (2018), 1901–1913.

Fabian

M.R.

et al., Regulation of mRNA translation and stability by microRNAs, Annu Rev Biochem 79 (2010), 351–79.

et al., Long noncoding RNA NORAD, a novel competing endogenous RNA, enhances the hypoxia-induced epithelial-mesenchymal transition to promote metastasis in pancreatic cancer, Mol Cancer 16(1) (2017), 169.

Zhang

et al., Exploration of the molecular mechanism of prostate cancer based on mRNA and miRNA expression profiles, Onco Targets Ther 10 (2017), 3225–3232.

Zhang

et al., Identification of a novel microRNA-mRNA regulatory biomodule in human prostate cancer, Cell Death Dis 9(3) (2018), 301.

Summerer

et al., Targeted high throughput sequencing of a cancer-related exome subset by specific sequence capture with a fully automated microarray platform, Genomics 95(4) (2010), 241–6.

Huang

et al., HMDD v3.0: a database for experimentally supported human microRNA-disease associations, Nucleic Acids Res 47(D1) (2019), D1013–D1017.

10.

Yang

et al., dbDEMC 2.0: updated database of differentially expressed miRNAs in human cancers, Nucleic Acids Res 45(D1) (2017), D812–D818.

11.

Xie

et al., miRCancer: a microRNA-cancer association database constructed by text mining on literature, Bioinformatics 29(5) (2013), 638–44.

12.

Revelle

W.R.

, psych: Procedures for personality and psychological research, Psychology (2017).

13.

Vlachos

I.S.

et al., DIANA-miRPath v3.0: deciphering microRNA function with experimental support, Nucleic Acids Res 43(W1) (2015), W460–6.

14.

Subramanian

et al., A Next Generation Connectivity Map: L1000 Platform and the First 1,000,000 Profiles, Cell 171(6) (2017), 1437–1452.e17.

15.

Murthi

et al., GAPDH, 18S rRNA and YWHAZ are suitable endogenous reference genes for relative gene expression studies in placental tissues from human idiopathic fetal growth restriction, Placenta 29(9) (2008), 798–801.

16.

Tang

et al., Selection and validation of reference genes for quantitative expression analysis of miRNAs and mRNAs in Poplar, Plant Methods 15 (2019), 35.

17.

Spandidos

et al., PrimerBank: a resource of human and mouse PCR primer pairs for gene expression detection and quantification, Nucleic Acids Res 38(Database issue) (2010), D792–9.

18.

Newman

A.M.

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

19.

Yoshihara

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

20.

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–8.

21.

Afonso-Grunz

and Muller

, Principles of miRNA-mRNA interactions: beyond sequence complementarity, Cell Mol Life Sci 72(16) (2015), 3127–41.

22.

Filella

and Foj

, Prostate Cancer Detection and Prognosis: From Prostate Specific Antigen (PSA) to Exosomal Biomarkers, Int J Mol Sci 17(11) (2016).

23.

Stoykova

G.E.

and Schlaepfer

I.R.

, Lipid Metabolism and Endocrine Resistance in Prostate Cancer, and New Opportunities for Therapy, Int J Mol Sci 20(11) (2019).

24.

Wang

et al., JAK/STAT3-Regulated Fatty Acid beta-Oxidation Is Critical for Breast Cancer Stem Cell Self-Renewal and Chemoresistance, Cell Metab 27(1) (2018), 136–150.e5.

25.

Kaiser

et al., The evolving role of diet in prostate cancer risk and progression, Curr Opin Oncol 31(3) (2019), 222–229.

26.

Gravina

G.L.

et al., KPT-330, a potent and selective exportin-1 (XPO-1) inhibitor, shows antitumor effects modulating the expression of cyclin D1 and survivin [corrected] in prostate cancer models, BMC Cancer 15 (2015), 941.

27.

Abdul

and Hoosein

, Inhibition by anticonvulsants of prostate-specific antigen and interleukin-6 secretion by human prostate cancer cells, Anticancer Res 21(3B) (2001), 2045–2048.

28.

Svensson

R.U.

Bannick

N.L.

Marin

M.J.

Robertson

L.W.

Lynch

C.F.

and Henry

M.D.

, Chronic Chlorpyrifos Exposure Does Not Promote Prostate Cancer in Prostate Specific PTEN Mutant Mice, J Environ Pathol Toxicol Oncol 32(1) (2013), 29–39.

29.

Tang

et al., Montelukast inhibits hypoxia inducible factor-1alpha translation in prostate cancer cells, Cancer Biol Ther 19(8) (2018), 715–721.

30.

Cormio

et al., Combined perianal-intrarectal (PI) lidocaine-prilocaine (LP) cream and lidocaine-ketorolac gel provide better pain relief than combined PI LP cream and periprostatic nerve block during transrectal prostate biopsy, BJU Int 109(12) (2012), 1776–80.

31.

Tannock

I.F.

et al., Docetaxel plus prednisone or mitoxantrone plus prednisone for advanced prostate cancer, N Engl J Med 351 (2004), 1502–12.

32.

Fan

et al., Clofoctol and sorafenib inhibit prostate cancer growth via synergistic induction of endoplasmic reticulum stress and UPR pathways, Cancer Manag Res 10 (2018), 4817–4829.

33.

Meng

W.J.

et al., MicroRNA Expression Profile Reveals miR-17-92 and miR-143-145 Cluster in Synchronous Colorectal Cancer, Medicine (Baltimore) 94(32) (2015), e1297.

34.

Moi

et al., Differential expression of the miR-17-92 cluster and miR-17 family in breast cancer according to tumor type; results from the Norwegian Women and Cancer (NOWAC) study, J Transl Med 17(1) (2019), 334.

35.

Liu

et al., The SOX4/miR-17-92/RB1 Axis Promotes Prostate Cancer Progression, Neoplasia 21(8) (2019), 765–776.

36.

X.G.

et al., Cancer-derived exosomal miR-221-3p promotes angiogenesis by targeting THBS2 in cervical squamous cell carcinoma, Angiogenesis 22(3) (2019), 397–410.

37.

et al., MiR-221-3p targets ARF4 and inhibits the proliferation and migration of epithelial ovarian cancer cells, Biochem Biophys Res Commun 497(4) (2018), 1162–1170.

38.

Krebs

et al., miR-221-3p Regulates VEGFR2 Expression in High-Risk Prostate Cancer and Represents an Escape Mechanism from Sunitinib In Vitro, J Clin Med 9(3) (2020).

39.

Zhou

et al., Plasma miRNAs in diagnosis and prognosis of pancreatic cancer: A miRNA expression analysis, Gene 673 (2018), 181–193.

40.

Zhang

et al., Combined Detection of Serum MiR-221-3p and MiR-122-5p Expression in Diagnosis and Prognosis of Gastric Cancer, J Gastric Cancer 19(3) (2019), 315–328.

41.

Larrabeiti-Etxebarria

et al., Systematic Review of the Potential of MicroRNAs in Diffuse Large B Cell Lymphoma, Cancers (Basel) 11(2) (2019).

42.

Haldrup

et al., FRMD6 has tumor suppressor functions in prostate cancer, Oncogene 40(4) (2021), 763–776.

43.

Mengual

et al., Using gene expression from urine sediment to diagnose prostate cancer: development of a new multiplex mRNA urine test and validation of current biomarkers, BMC Cancer 16 (2016), 76.

44.

Stoen

M.J.

et al., Overexpression of miR-20a-5p in Tumor Epithelium Is an Independent Negative Prognostic Indicator in Prostate Cancer-A Multi-Institutional Study, Cancers (Basel) 13(16) (2021).

45.

Kneitz

et al., Survival in patients with high-risk prostate cancer is predicted by miR-221, which regulates proliferation, apoptosis, and invasion of prostate cancer cells by inhibiting IRF2 and SOCS3, Cancer Res 74(9) (2014), 2591–603.

46.

Quero

et al., miR-221-3p Drives the Shift of M2-Macrophages to a Pro-Inflammatory Function by Suppressing JAK3/STAT3 Activation, Front Immunol 10 (2019), 3087.

47.

et al., miR-20a-5p/TGFBR2 Axis Affects Pro-inflammatory Macrophages and Aggravates Liver Fibrosis, Front Oncol 10 (2020), 107.

Analysis of aberrant miRNA-mRNA interaction networks in prostate cancer to conjecture its molecular mechanisms

Abstract

BACKGROUND:

OBJECTIVE:

METHODS:

RESULTS:

CONCLUSIONS:

Keywords

1. Introduction

2. Materials and method

2.1 Data collection and processing

2.3 Forecasting of miRNA-mRNA interaction

2.4 Sample collection

Table 1 Clinicopathological and molecular features of prostate cancer patients

2.6 Estimation of miRNA-mRNA interactions with tumor-related phenotypes

2.7 Diagnostic effectiveness evaluation and overall survival analysis

2.8 Statistically analysis

3. Results

3.1 Identification of DE-mRNAs and DE-miRNAs

Table 3 Pearson’s correlation analysis of miRNA-mRNA pairs in prostate cancers in TCGA

Table 4 Pearson’s correlation analysis of miRNA-mRNA pairs in FFPE prostate cancer samples

3.6 Comparing miRNA-mRNA expression level between clinical subgroups in PRAD

3.8 Overall survival analysis

4. Discussion

5. Conclusions

Author contributions

Supplementary data

sj-pdf-1-cbm-10.3233_CBM-220051.pdf - Supplemental material

Footnotes

Acknowledgments

References

Table 1
Clinicopathological and molecular features of prostate cancer patients

Table 3
Pearson’s correlation analysis of miRNA-mRNA pairs in prostate cancers in TCGA

Table 4
Pearson’s correlation analysis of miRNA-mRNA pairs in FFPE prostate cancer samples