Sage Journals: Discover world-class research

Abstract

BACKGROUND:

MicroRNAs regulating mRNA expression by targeting at mRNAs is known constructive in tumor occurrence, immune escape, and metastasis.

OBJECTIVE:

This research aims at finding negatively regulatory miRNA-mRNA pairs in esophageal squamous cell carcinoma (ESCC).

METHODS:

GENE expression data of The Cancer Genome Atlas (TCGA) and GEO database were employed in differently expressed RNA and miRNA (DE-miRNAs/DE-mRNAs) screening. Function analysis was conducted with DAVID-mirPath. MiRNA-mRNA axes were identified by MiRTarBase and TarBase and verified in esophageal specimen by real-time reverse transcription polymerase chain reaction (RT-qPCR). Receiver operation characteristic (ROC) curve and Decision Curve Analysis (DCA) were applied in miRNA-mRNA pairs predictive value estimation. Interactions between miRNA-mRNA regulatory pairs and immune features were analyzed using CIBERSORT.

RESULTS:

Combining TCGA database, 4 miRNA and 10 mRNA GEO datasets, totally 26 DE-miRNAs (13 up and 13 down) and 114 DE-mRNAs (64 up and 50 down) were considered significant. MiRTarBase and TarBase identified 37 reverse regulation miRNA-mRNA pairs, 14 of which had been observed in esophageal tissue or cell line. Through analysis of RT-qPCR outcome, miR-106b-5p/KIAA0232 signature was chosen as characteristic pair of ESCC. ROC and DCA verified the predictive value of model containing miRNA-mRNA axis in ESCC. Via affecting mast cells, miR-106b-5p/KIAA0232 may contribute to tumor microenvironment.

CONCLUSIONS:

The diagnostic model of miRNA-mRNA pair in ESCC was established. Their complex role in ESCC pathogenesis especially tumor immunity was partly disclosed.

Keywords

Biomarkers miRNA-mRNA interaction esophageal squamous cell carcinoma TCGA GEO PCR

1. Introduction

As is estimated by World Health Organization, esophageal cancer became the seventh most common cancer and the sixth leading cause of cancer-associated deaths globally in 2020, making esophagus cancer a worldwide health problem [1]. Esophagus tumors are mainly divided into esophageal adenocarcinoma (EA) and esophageal squamous cell carcinoma (ESCC) according to histological origin [2]. Esophageal squamous cell carcinoma accounts for the majority of newly diagnosed esophageal cancer in the world [3]. Despite the swift development in diagnosis and treatment of ESCC, the overall 5-year survival rate ranged from 20% to 35% [4]. Thus, it is of significance to define novel biomarker for better diagnosis and prognosis prediction.

MicroRNAs (MiRNAs) are a serious of small non-coding RNAs involved in tumor occurrence, progression, metastasis and immune escape [5]. By incompletely binding with 3’-untranslated regions of message RNA (mRNA), miRNAs usually act as suppresser in mRNA translation and stability [6]. Dysregulation of miRNAs exert profound effects on expression of proto-oncogene and anti-oncogene, which is generally accepted as the cause of cancer [7]. Besides, miRNAs can serve as biomarkers in prognosis prediction. For instance, a risk model including miR-4521, miR-3682-3p, miR-1269a, and 7 target mRNAs were established to estimate prognosis in esophagus cancer patients [8]. Therefore, biomarkers such as miRNA/mRNA signature and miRNA targeting therapies are promising in future cancer research.

Contemporarily, databases obtaining gene expression patterns and non-coding RNA profiling from both high throughput sequencing and microarrays provide an overview of mRNA and miRNA abundance in tumor. Bioinformatic analysis based on large sample points out the direction for further experimental research. Previous studies absorbed a few Gene Expression Omnibus (GEO) datasets to define differently expressed miRNAs (DE-miRNAs) and miRNA-mRNA axes for ESCC. Regulatory networks such as miR-139-5p targeting at RAP1B and miR-1 targeting at FN1 were identified as crucial axes in ESCC [9, 10]. However, few studies contained both large amount of GEO datasets and The Cancer Genome Atlas (TCGA) database. In present study, we collected miRNA and mRNA expression profile of esophageal carcinoma and normal esophagus tissue from 14 Gene Expression Omnibus (GEO) datasets (http://www.ncbi.nlm.nih.gov/geo/) and TCGA-ESCC database (https://portal.gdc.cancer.gov/) to select out differently expressed miRNA (DE-miRNA) and differently expressed mRNA (DE-mRNA). TarBase and MiRTarBase were employed to estimate interactions between miRNA and mRNA. Then, we explored the expression level of these candidate miRNA-mRNA axes in 30 paired esophagus specimens via reverse transcription quantitative polymerase chain reaction (RT-qPCR). Clinical subgroup analysis and survival analysis were conducted with the chosen miRNA-mRNA pair to measure its predictive value. MiRNA-mRNA axis may provide a new idea for selecting biomarkers valuable in tumor diagnosis and prognosis.

2. Materials and methods

2.1 Data acquisition from databases

Table 1
Information pertaining to the selected GEO datasets for esophagus cancer

	Specimen type	GEO accession	Number of tumors	Number of non-tumors	Total number	Country	Platform
miRNA array	Tissue	GSE43732	119	119	238	China	GPL16543
		GSE114110	30	10	40	China	GPL24967
		GSE6188	153	104	257	China	GPL4508
		GSE13937	44	44	88	USA	GPL8835
mRNA array	Tissue	GSE33426	59	12	71	USA	GPL571
		GSE17351	5	5	10	USA	GPL570
		GSE100942	4	4	8	China	GPL570
		GSE77861	7	7	14	USA	GPL570
		GSE38129	30	30	60	USA	GPL571
		GSE89102	5	5	10	China	GPL16956
		GSE26886	9	19	28	Germany	GPL570
		GSE44021	40	40	40	USA	GPL96
		GSE161533	28	28	56	China	GPL570
		GSE23400	51	51	102	USA	GPL96

Figure 1.

The flow chart showed the steps for identifying the negatively miRNA-mRNA regulatory pairs and conducting extensive analysis of pairs in esophageal squamous cell carcinoma cancer.

The whole process of this study is displayed in Fig. 1. Esophageal carcinoma RNA-Sequencing data and clinical information of The Cancer Genome Atlas (TCGA) were obtained from the Genomic Data Commons portal in the National Cancer Institute (https:// portal.gdc.cancer.gov/). Via setting the search keywords as ‘(esophagus OR esophageal) AND (tumor OR cancer OR squamous cell carcinoma OR carcinoma) NOT (cell line)’, we downloaded 10 mRNA expression datasets and 4 miRNA expression datasets from the GEO database (http://www.ncbi.nlm.nih.gov/geo/). All the datasets were mRNA or miRNA expression profiles measured by microarray in esophageal squamous cell carcinoma tissue and normal tissue from homo sapiens (listed in Table 1). Differently-expressed miRNAs collected in open-access databases including Human MicroRNA Disease Database (HMDD) (https:// www.cuilab.cn/hmdd), miRCancer (http://mircancer. ecu.edu/) and Database of Differentially Expressed MiRNAs in human Cancers (dbDEMC) (https://www. biosino.org/dbDEMC/index) were also involved [11, 12, 13].

2.2 Analysis of miRNA and mRNA expression

MRNA-sequencing and miRNA-sequencing profile gained from TCGA were processed with psych, an R-based package. Wilcoxon rank-sum test were performed between ESCC and normal tissue to calculate fold change (FC) of differently expressed mRNAs (DE-mRNA) and miRNAs (DE-miRNA). GEO2R (https://www.ncbi.nlm.nih.gov/geo/geo2r/), an online data processing tool, was employed in analyzing miRNA and mRNA expression profiles downloaded from GEO database $P$ -value $<$ 0.05 and $|\text{log}_{2}\text{FC}|>$ 0.58 were considered as statically significant.

2.3 Gene oncology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses

Table 2
The sequences of primers for candidate miRNAs and targeted mRNAs

	Name	Forword primer sequences (5’ $\to$ 3’)	Reverse primer sequences (5’ $\to$ 3’)
miRNA	miR-19b-3p	TGTGCAAATCCATGCAAAA	Universal reverse primer
	miR-155-5p	TTAATGCTAATCGTGATA	Universal reverse primer
	miR-93-5p	CAAAGTGCTGTTCGTGCAGGT	Universal reverse primer
	miR-21-5p	TAGCTTATCAGACTGATGTT	Universal reverse primer
	miR-106b-5p	TAAAGTGCTGACAGTGCA	Universal reverse primer
	miR-143-3p	TGAGATGAAGCACTGTAGCT	Universal reverse primer
	miR-145-5p	GTCCAGTTTTCCCAGGAATC	Universal reverse primer
	miR-30a-5p	TGTAAACATCCTCGACTGGA	Universal reverse primer
	miR-497-5p	CAGCAGCACACTGTGGTTT	Universal reverse primer
	miR-1-3p	TGGAATGTAAAGAAGTATGT	Universal reverse primer
	miR-29c-3p	TAGCACCATTTGAAATCGG	Universal reverse primer
	miR-195-5p	TAGCAGCACAGAAATATTG	Universal reverse primer
mRNA	COL1A1	GAGGGCCAAGACGAAGACATC	CAGATCACGTCATCGCACAAC
	MMP1	GGGGCTTTGATGTACCCTAGC	TGTCACACGCTTTTGGGGTTT
	ADAM12	AACCTCGCTGCAAAGAATGTG	CTCTGAAACTCTCGGTTGTCTG
	KIAA0232	TGATGAAAGCTCTGGTATCGAGA	GGAGACCCCTCTAACTTTTTGTG
	SASH1	CGGTCCCAGATCGAAGAGTC	GTTCTTTCGGAAGTTCTGCCA
	VAT1	GGAGGCTACGACAAGGTGAAG	CATGAGGTCTGCGAAGTTGAG
	PPL	GCACCAATGAGCTGTACTGG	GCTGGGGTAGTCGAGGTTG
	SOX4	AGCGACAAGATCCCTTTCATTC	CGTTGCCGGACTTCACCTT
	BIRC5	AGGACCACCGCATCTCTACAT	AAGTCTGGCTCGTTCTCAGTG
	MXD1	CGGGCTCATCTTCGCTTGT	GATTTGGTGAACGGCTTTTCTG
Reference mRNA	GAPDH	GGAGCGAGATCCCTCCAAAAT	GGCTGTTGTCATACTTCTCATGG
	18sRNA	CGCAGCTAGGAATAATGGAATAGG	GCCTCAGTTCCGAAAACCAA
	U6	CGATAAAATTGGAACGATACAGA	ATTTGGACCATTTCTCGATTTGT

GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses of DE-mRNA was carried out with the cluster Profiler tool on Hiplot (https://hiplot.com.cn). DAVID-mirPath was employed to conduct GO/KEGG analysis of upregulated and downregulated DE-miRNA in ESCC [14] We set the screening criteria as $P$ -value $<$ 0.05 and counts $>$ 2.

2.4 Forecasting the correspondence between miRNA and mRNA

MiRTarBase v.8 (http://miRTarBase.cuhk.edu.cn/) and TarBase v.8 (http://www.microrna.gr/tarbase) collected miRNAmRNA interaction networks verified by experiments. MiRNA-mRNA axes collected in both MiRTarBase and TarBase were used for screening potential target mRNAs. Next, we searched each miRNA-mRNA axis in PubMed to confirm whether these reverse regulation phenomenon between miRNA and mRNA can be observed in esophageal tissue esophageal cell line or esophageal cancer cell line. The miRNA-mRNA interactions were drawn with Cytoscape v3.8.2, an open-source software for complex network [15].

2.5 Sample collection, RNA isolation and reverse transcription quantitative polymerase chain reaction (RT-qPCR)

Totally 30 paired ESCC and normal esophagus Formalin-Fixed and Paraffin-Embedded (FFPE) tissue were collected. All patients underwent esophageal surgery at First Affiliated Hospital of Nanjing Medical University before receiving any chemotherapy or radiation therapy. This study was approved by the Institutional Review Boards of the First Affiliated Hospital of Nanjing Medical University (ID: 2016-SRFA-148). All patients were informed about this research and signed a consent statement before starting the study.

All FFPE specimen were stored at 20 ${}^{\circ}$ C. To isolate total RNA from FFPE esophageal samples, we used RNAprep Pure FFPE Kit (TIANGEN Biotech, Beijing, China, DP439) after consulting the manual. Then, a poly(A) tail was added to all RNA with Poly(A) Polymerase Kit (Takara). RT-qPCR was performed with Prime Script RT reagent Kit (Takara) and SYBR Premix Ex Taq II (Takara) following the instructions. Entire PCR process were carried out in final volumes of 10 $\mu$ L on the qTOWER384 (Analytik Jena) plate with the temperature set as 95 ${}^{\circ}$ C for 20 s, followed by 40 cycles made up of 10 s at 95 ${}^{\circ}$ C and 20 s at 60 ${}^{\circ}$ C. We used GAPDH and 18S as reference gene and U6 as miRNA’s internal control [16, 17]. All sequences of primer were obtained from Primer Bank (https://pga.mgh.harvard.edu/primerbank/) and listed in Table 2.

2.6 Estimation of miRNA-mRNA axes’ predictive value in ESCC and survival analysis

SPSS Statistics 26.0 (IBM, New York, USA) was employed to perform binary logistic regression with expression data of selected miRNA-mRNA pairs. Receiver operating characteristic (ROC) curves and decision curve analysis (DCA) were respectively drawn by GraphPad Prism and rmda, an R-based package to estimate the prognostic value of miRNA-mRNA pairs. Owing to the lack of GEO survival data, survival information including survival data were downloaded from TCGA. High expression group and low expression group were distinguished by median. Survival analysis curves were protracted with Hiplot (https:// hiplot.com.cn).

2.7 Prediction of pertinence between miRNA-mRNA pairs and Tumor phenotypes

We used CIBERSORT (https://cibersort.stanford.edu/ index.php/) to predict the proportion of 22 types of immune cells in both ESCC tissue and normal esophageal tissue with RNA-sequencing data acquired from TCGA. The Wilcoxon-Mann-Whitney (WMW) test conducted by psych was processed to compare differences of immune cell infiltration between cancer and normal control. Correlation between miRNA-mRNA pairs and infiltrating immune cell were evaluated by correlation heatmap in Hiplot (https://hiplot.com.cn) with Spearman’s correlation. We applied ESTIMATE in forecasting tumor purity and calculating the Stromal Score, Immune Score and ESTIMATE Score of TCGA-ESCC samples. Tumor mutational burden (TMB) was measured by Masked Somatic Mutation data (varscan. Somatic. Maf) obtained with the ‘maftools’, an R-based package. We also employed correlation heatmap in Hiplot (https://hiplot.com.cn) to forecast the correlation between miRNA-mRNA axes and cancer phenotypes.

2.8 Statistically analysis

The expression data of miRNA and mRNA gained by RT-PCR were calculated by the $2^{-\Delta\Delta\text{CT}}$ method [18]. Variance analysis of expression profile between normal tissue and esophagus tumor was conducted by paired T-test. The correlations between miRNA and mRNA were assessed by Spearman’s correlation with psych. The $P$ -value $<$ 0.05 represented statistical significance. SPSS Statistics v26.0(IBM, New York, USA), GraphPad prism v8.0.1 and R language v4.0.1 were applied in data analyses. Cytoscape, Hiplot, GraphPad prism and R language were adopted in plotting figures.

3. Results

3.1 Identification of DE-mRNA and DE-miRNA via GEO and TCGA database

Figure 2.

(A) The circular bar-plot exhibited the information of 14 GEO datasets. All 10 miRNA datasets and 4 mRNA datasets were drawn from ESCC tissues. (B) The negative regulatory relationship screened by TarBase and MiRTarBase between 26 microRNAs (miRNAs) and 114 mRNAs were plotted by Cytoscape. 10 miRNA (up)-mRNA (down) pairs and 27 miRNA (down)-mRNA (up) pairs were selected. Orange nodes represent the upregulated miRNAs/mRNAs in ESCCs versus NCs, while blue nodes represent the downregulated miRNAs/mRNAs in ESCCs versus NCs. (C) Totally 14 negatively miRNA-mRNA regulatory axes were validated by RT-qPCR in 30 paired FFPE samples.

Table 3

Pearson’s correlation analysis of miRNA-mRNA pairs in paired FFPE ESCC samples

		Fold change			Fold change	Pearson’s correlation
miRNA (up)	$p$ -value	(2^ $-\Delta\Delta$ CT)	mRNA (down)	$p$ -value	(2^ $-\Delta\Delta$ CT)	$p$ -value	$r$ -value
hsa-miR-106b-5p	0.034	13.978	KIAA0232	$<$ 0.0001	0.518	0.042	$-$ 0.251
hsa-miR-21-5p	0.111	187.668	SASH1	$<$ 0.0001	0.437	0.816	$-$ 0.086
hsa-miR-93-5p	0.027	74.715	SASH1	$<$ 0.0001	0.437	0.491	0.086
			VAT1	0.001	0.560	0.967	$-$ 0.009
hsa-miR-155-5p	0.092	71.147	PPL	0.006	0.656	0.092	$-$ 0.344
hsa-miR-19b-3p	0.169	192.548	MXD1	0.133	2.151	0.639	0.127
miRNA (down)			mRNA (up)			$p$ -value	$r$ -value
hsa-miR-143-3p	0.131	1.750	COL1A1	0.007	2.116	0.895	0.050
hsa-miR-145-5p	0.278	8.197	MMP1	0.138	76.745	0.132	0.294
hsa-miR-30a-5p	0.014	0.669	SOX4	0.169	7.494	0.785	0.069
hsa-miR-497-5p	0.542	1.395	BIRC5	0.155	8.411	0.383	0.191
hsa-miR-1-3p	0.055	0.683	ADAM12	0.316	78.717	0.125	0.274
hsa-miR-29c-3p	0.608	0.902	ADAM12	0.316	78.717	0.649	$-$ 0.127
			COL1A1	0.007	2.116	0.937	$-$ 0.015
hsa-miR-195-5p	0.375	1.197	BIRC5	0.155	8.411	0.816	$-$ 0.081

Figure 3.

(A) The miRNA expression levels of miR-93-5p and miR-106b-5p climbed while miR-30a-5p declined in ESCCs. (B) The mRNA expression levels of KIAA0232, SASH1, VAT1 and PPL were downregulated, while the expression levels of COL1A1 were upregulated. (C–D) Pearson’s correlation analysis of miRNA-mRNA regulatory pairs in 30 paired samples. MiR-106b-5p/KIAA0232 were the two negatively regulated pairs displayed in figure. Data were presented as mean $\pm$ SEM. ${}^{*}p<$ 0.05, ${}^{**}p<$ 0.01, ${}^{***}p<$ 0.001 and ${}^{****}p<$ 0.0001 (Student’s $t$ -test). $P$ -values are listed in Table 3.

Figure 4.

Receiver operating characteristic (ROC) curves of external validation and TCGA database was separately presented in Fig. 4A and B. (A) The areas under the curve of complex model, miR-106b-5p and KIAA0232 were 0.771, 0.681 and 0.763 respectively in external validation. (B) The areas under the curve of complex model, miR106b-5p and KIAA0232 were 0.914, 0.908 and 0.733 respectively. (C–D) DCA presented satisfied diagnostic performance of miRNA-mRNA network in external validation and TCGA database.

Wilcoxon rank test conducted in TCGA-ESCC database determines the primary set of DE-mRNAs and DE-miRNAs. As is presented in Fig. 2A, totally 10 mRNA and 4 miRNA expression datasets measured by microarray met our requirements. The detailed information of 14 GEO datasets was exhibited in Table 1 After calculating with GEO2R, 64 upregulated mRNAs and 50 downregulated mRNAs appeared in all 10 datasets. After comparing to normal samples, 38 upregulated miRNAs were mentioned in 3 more GEO datasets and 31 down-regulated genes appeared in 2 or more datasets. Taking GEO datasets, TCGA, dbDEMC, HMDD and miRCancer into consideration, 13 upregulated miRNAs and 13 downregulated miRNAs were also screened out.

3.2 GO and KEGG pathway enrichment analysis

Results of GO/KEGG analysis based on DE-miRNAs are revealed in Fig. S1A and S1B. For downregulated miRNAs, Proteoglycans in cancer, Pathways in cancer, ECM-receptor interaction were the most significant KEGG pathways. For upregulated miRNAs, Endocytosis, Melanoma and Hepatitis C are the enriched pathways with smallest $P$ -value. Extracellular vesicular exosome and mitochondrial matrix were GO terms of downregulated miRNAs with the most counts. Ion binding was the GO term of upregulated miRNAs with the most counts and the smallest $P$ -value. Outcomes of GO analysis based on DE-mRNAs were separately revealed in Fig. S1C and S1D. KEGG analysis disclosed that upregulated mRNAs were enriched in DNA replication.

3.3 Determination of reverse regulation miRNA-mRNA interactions in ESCC

Via consulting MiRTarBase and TarBase, we identified 37 reverse regulation miRNA-mRNA pairs(Fig. 2B). After searching each miRNA-mRNA axis in PubMed, 14 pairs, including 6 miRNA (up)-mRNA (down) pairs and 8 miRNA (down)-mRNA (up) pairs, were verified in esophageal tissue or cell line (Fig. 2C).

3.4 Validating expression of mRNAs and miRNAs in ESCC FFPE tissue

Expression data of 14 miRNA-mRNA axes, including 10 mRNAs and 12 miRNAs, were measured in 30 paired ESCC tissues and adjacent normal tissues via RT-qPCR. Specific outcome of all miRNAs and mRNAs was listed in Table 3. Figure 3A and B separately presented miRNA and mRNA expression data. Hsa-miR-93-5p ( $p=$ 0.027) and hsa-miR-106b-5p ( $p=$ 0.034) were observed upregulated in tumor tissue. Hsa-miR-30a-5p ( $p=$ 0.014) expression of ESCC tissue was lower than that of adjacent esophageal tissue. KIAA0232 ( $p<$ 0.0001), SASH1 ( $p<$ 0.0001), VAT1 ( $p=$ 0.001) and PPL ( $p=$ 0.006) were downregulated in cancer tissue. Only COL1A1 ( $p=$ 0.007) was observed overexpressed in ESCC tissue. According to Pearson’s correlation analysis, KIAA0232 expression correlated with miR-106b-5p significantly ( $p=$ 0.046, $r=$ $-$ 0.25) (Fig. 3C and D).

3.5 Estimating diagnostic value of miRNA-mRNA pairs in ESCC

For lack of GEO datasets containing both miRNA and mRNA expression pattern, we measured predictive value of the miRNA-mRNA pair with expression profile gained from FFPE tissue by PCR and TCGA database. As is displayed in Fig. 4A–B, we drew receiver operating characteristic (ROC) curves and decision curve analysis (DCA) of hsa-miR-106b-5p, KIAA0232 and a model containing miR-106b-5p/KIAA0232 axis separately. The areas under the curve (AUC) of complex model were 0.7711 (95% CI: 0.6498 to 0.8924) in PCR and 0.9142 (95% CI: 0.8252 to 1.000) in TCGA, maximum among all 3 curves. DCA suggested satisfactory diagnosis value of miRNA-mRNA axis in TCGA-ESCC database (Fig. 4D).

3.6 Clinical subgroup analysis and survival analysis with miRNA/mRNA expression

Table 4
Analysis of CIBERSORT scores of 22 types of immune cells in esophagus cancers versus controls

	ESCC score	Control score	$P$ -value	FC
B cells naive	0.0469	0.0690	0.4487	0.6795
B cells memory	0.0022	0.0123	0.2626	0.1789
Plasma cells	0.0524	0.1084	0.0124	0.4831
T cells CD8	0.0645	0.1108	0.0081	0.5822
T cells CD4 naive	0.0003	0.0001	0.8127	2.7478
T cells CD4 memory resting	0.2229	0.2365	0.5076	0.9422
T cells CD4 memory activated	0.0366	0.0087	0.0100	4.2013
T cells follicular helper	0.0303	0.0146	0.0438	2.0685
T cells regulatory (Tregs)	0.0450	0.0388	0.2823	1.1607
T cells gamma delta	0.0007	0.0001	0.5236	7.1618
NK cells resting	0.0234	0.0196	0.4673	1.1915
NK cells activated	0.0173	0.0296	0.1976	0.5830
Monocytes	0.0070	0.0119	0.0074	0.5914
Macrophages M0	0.1502	0.0250	$<$ 0.0001	6.0197
Macrophages M1	0.0659	0.0206	0.0013	3.2049
Macrophages M2	0.1006	0.0923	0.2485	1.0896
Dendritic cells resting	0.0309	0.0241	0.4266	1.2853
Dendritic cells activated	0.0383	0.0215	0.0580	1.7820
Mast cells resting	0.0266	0.1448	$<$ 0.0001	0.1836
Mast cells activated	0.0272	0.0058	0.1137	4.6921
Eosinophils	0.0041	0.0067	0.5118	0.6055
Neutrophils	0.0089	0.0010	0.0142	9.1845

Figure 5.

Expression level of miRNA-mRNA pairs were compared in subgroups divided by clinical features of ESCC patients acquired from TCGA database. Alcohol exposure, gender, age of diagnosis and tumor status had no relationship with miR-106b-5p expression. KIAA0232 expression level was higher in female patients.

Figure 6.

(A) The bar graph of 22 immune cells differentially infiltrated fraction in ESCC versus NCs were displayed. (B) The correlation heatmap between cell fraction of 9 differently infiltrated immune cells and expression level of miRNA-mRNA pairs was constructed by Pearson’s correlation. (C) The relevance between expression level of miRNAs and targeting mRNAs and 5tumor microenvironment factors (TMB, ESTIMATE score, immune score, stromal score, and tumor).

All the clinical cases collected in TCGA-ESCC database were grouped by alcohol exposure, gender, age, and tumor stage respectively. As is shown in Fig. 5, no expression difference of miR-106b-5p was observed between all four groups. KIAA0232 expression was higher in female patients than in male ( $p=$ 0.040), while exposure to alcohol, age, and tumor status of patients were not related with KIAA0232 expression. Survival analysis was conducted based on overall survival (OS) acquired from TCGA-ESCC database. Kaplan-Meier (K-M) survival curves suggested that expression level of neither KIAA0232 nor miR-106b-5p was associated with prognosis (Fig. S2). The date of loss of follow-up for each case was also displayed in Fig. S2.

3.7 Correlation analysis between miRNA/mRNA expression and tumor-related phenotypes

There was statistically significant difference in the proportion of 9 types of immune cells between ESCC tissue and normal esophageal tissue (Fig. 6A and Table 4). Correlation heatmap between miRNA/mRNA expression and immune cell proportion was plotted in Fig. 6B. Percentage of resting mast cells was positively correlated with KIAA0232 expression ( $p<$ 0.05, $r=$ 0.15) while negatively correlated with miR-106b-3p expression ( $p<$ 0.05, $r=$ $-$ 0.36). Follicular helper T cell was only correlated with miR-106b-5p ( $p<$ 0.05, $r=$ 0.23). As is presented in Fig. 6C, miR-106b-5p expression level was relevant with tumor immune microenvironment and tumor purity. No correlation was found between KIAA0232 expression and tumor-related phenotypes.

4. Discussion

In past few decades, the essential role of miRNA dysregulation in tumor occurrence, and immune microenvironment has been gradually recognized by scholars [19, 20]. MiRNAs post-transcriptionally inhibit expression of target genes A number of miRNAs such as miR-196a-5p and miR-1-3p has been identified as candidates for ESCC diagnosis and therapy targets [21] However, miRNA-mRNA interactions and their roles in ESCC have not been fully disclosed. Firstly, we collected expression pattern in normal and esophageal squamous cell carcinoma tissue from 4 miRNA GEO datasets, 10mRNA GEO datasets and TCGA databases Expression from GEO database and TCGA-ESCC database were analyzed by GEO2R and Wilcoxon signed-rank test respectively. Taking several miRNA differential expression libraries including HMDD, miRCancer and dbDEMC into consideration, we ultimately defined 26 DE-miRNAs (13 upregulated and 13 downregulated miRNAs) and 114 DE-mRNAs (64 upregulated and 50 downregulated mRNAs). Totally 37 miRNA-mRNA pairs were confirmed using MiRTarBase and TarBase, databases collecting experimentally supported miRNA-mRNA networks, 14 of which were observed in esophagus tissue or cell line. We further detected the expression of candidate DE-miRNAs and DE-mRNAs in 30 paired FFPE esophagus samples by poly(A) RT-qPCR. After the processing of experimental outcomes, the miR-106b-5p/KIAA0232 emerged as specific miRNA/mRNA signature of ESCC ultimately.

GO/KEGG analysis based on 26 DE-miRNAs and 114 DE-mRNAs emphasized the role of DNA replication, virus infection and energy metabolism. Notably, a few well-known tumor-associated cell signaling pathways, such as WNT signaling pathway, MAPK signaling pathway, ErbB signaling pathway, are involved in ESCC progression, which is in accordance with previous studies [22, 23, 24].

In present study, has-miR-106b-5p was observed upregulated in ESCC. Previous study has disclosed that miR-106b-5p is overexpressed in ESCC tissues and cell lines and that silencing of miR-106b-5p inhibits tumor growth and pulmonary metastasis of ESCC [25]. MiR-106b-5p was covered in a disease progression risk estimation model of stage II colon cancer [26]. In addition, miR-106b-5p was reported overexpressed in esophagus adenocarcinoma [27]. MiR-106b belongs to the miR-106b-25 cluster (miR-106b, miR-25 and miR-93), which is involved in a variety of cancers such as brain cancer, breast cancer and hepatocellular carcinoma [28, 29, 30]. However, no correspondences were found between miR-106b-5p and patients’ clinical features including alcohol exposure, gender, age and tumor status. Through correlation analysis, we identified KIAA0232 as the downregulated target mRNA of miR-106b-5p. KIAA genes are uncharacterized protein-coding human genes by identified Kazusa cDNA project [31]. KIAA0232 expression was reported median in esophagus squamous epithelial cells according to The Human Protein Atlas database (https://www.proteinatlas.org/) The specific function of KIAA0232 protein has not been fully identified. According to GeneCards database (https://www.genecards.org/), KIAA0232 was predicted to enable ATP binding activity KIAA0232 expression of pancreatic beta cell is involved in insulin secretion stimulated by Glucagon-like peptide 1 [32]. However, there are few studies focusing on KIAA0232 protein and its role in tumor In this study, expression level of KIAA0232 was found higher in female patients than in male patients. The role of KIAA0232 protein and its impact on tumor occurrence, immune microenvironment, and metastasis deserve further exploration.

The diagnosis value of KIAA0232, miR-106b-5p and complex panel combining miRNA/mRNA signature were estimated by ROC and DCA curve This study indicated that the AUC of complex panel was 0.9142 (95% CI: 0.8252 to 1.000) and 0.7711 (95% CI: 0.6498 to 0.8924) in TCGA-ESCC and the external validation respectively. Previous studies recommended miR-106b as a member of 8-miRNA signature diagnosis panel in serum specimens [33]. The diagnostic value of miRNA-mRNA signature was affirmed by the decision curve analysis in TCGA-ESCC database.

The expression level of miR-106b-5p and KIAA0232 separately had negative and positive correlation with resting mast cells proportion. It is revealed in previous research that high mast cell density is associated with ESCC progression and reduce postoperative survival [34]. Mast cell expressing interleukin 17 in tumor tissues could function as prognostic marker and therapeutic target for ESCC [35]. Thus, it is reasonable to conclude that miR-106b-5p targeting at KIAA0232 may have influence on tumor immune microenvironment.

Unfortunately, the Kaplan-Meier (K-M) survival curves based on overall survival data of TCGA demonstrated that neither miR-106b-5p nor KIAA0232 could serve as prognosis predictor of ESCC. Similarly, the serum levels of miR-106b showed no correspondence with ESCC prognosis in research conducted by Isamu Hoshino et al. [36]. No previous research investigated the association between KIAA0232 and prognosis in ESCC. Due to lack of sufficient clinical survival data, detailed study on prognostic biomarkers were thwarted. Establishing novel survival prognostic model based on miRNA/mRNA axes needs to be further acknowledged in the future.

To sum up, this study investigated DE-mRNA, DE-miRNA, and their interactions in ESCC. The role of miRNA-mRNA interacting axes in tumor progression were emphasized in massive studies due to well-developed bioinformatics approaches. However, molecular mechanisms of how DE-mRNAs and DE-miRNAs affect tumor growth and immune escape deserve deeper dig into details.

5. Conclusions

In summary, this research aimed at finding the key miRNA-mRNA axes for ESCC by bioinformatic methods and verifying them in PCR. Our study may provide novel biomarkers for diagnosis and prognosis of ESCC in future.

Author contributions

Conception: Yulin Yue, Tongshan Wang, and Wei Zhu.

Interpretation or analysis of data: Shuang Peng and Shiyu Zhang.

Preparation of the manuscript: Shuang Peng, Shiyu Zhang and Xingchen Fan.

Revision for important intellectual content: Shuang Peng, Jingfeng Zhu and Cheng Liu.

Supervision: Wei Zhu.

Supplementary data

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

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

Supplemental material, sj-pdf-1-cbm-10.3233_CBM-220309.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.

Ajani

J.A.

et al., Esophageal and Esophagogastric Junction Cancers, Version 2.2019, NCCN Clinical Practice Guidelines in Oncology, J Natl Compr Canc Netw 17(7) (2019), 855–883.

Thrift

A.P.

, Global burden and epidemiology of Barrett oesophagus and oesophageal cancer, Nature Reviews Gastroenterology & Hepatology 18(6) (2021), 432–443.

Huang

F.L.

and Yu

S.J.

, Esophageal cancer: Risk factors, genetic association, and treatment, Asian J Surg 41(3) (2018), 210–215.

Rupaimoole

and Slack

F.J.

, MicroRNA therapeutics: towards a new era for the management of cancer and other diseases, Nature Reviews Drug Discovery 16(3) (2017), 203–222.

Fabian

M.R.

et al., Regulation of mRNA Translation and Stability by microRNAs, Annual Review of Biochemistry 79(1) (2010), 351–379.

Ali Syeda

et al., Regulatory Mechanism of MicroRNA Expression in Cancer, Int J Mol Sci 21(5) (2020).

Zhao

et al., A novel prognostic mRNA/miRNA signature for esophageal cancer and its immune landscape in cancer progression, Mol Oncol 15(4) (2021), 1088–1109.

Wei

et al., microRNA-messenger RNA regulatory network of esophageal squamous cell carcinoma and the identification of miR-1 as a biomarker of patient survival, J Cell Biochem 120(8) (2019), 12259–12272.

10.

Wang

et al., Bioinformatics-Based Identification of a circRNA-miRNA-mRNA Axis in Esophageal Squamous Cell Carcinomas, J Oncol 2020 (2020), 8813800.

11.

Huang

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

12.

Xie

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

13.

Yang

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

14.

Vlachos

I.S.

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

15.

Shannon

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

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.

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.

18.

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.

19.

Cho

et al., MicroRNAs and Their Targetomes in Tumor-Immune Communication, Cancers (Basel) 12(8) (2020).

20.

Iorio

M.V.

and Croce

C.M.

, MicroRNA dysregulation in cancer: diagnostics, monitoring and therapeutics. A comprehensive review, EMBO Mol Med 4(3) (2012), 143–59.

21.

Sang

et al., Identification and validation of hub microRNAs dysregulated in esophageal squamous cell carcinoma, Aging (Albany NY) 12(10) (2020), 9807–9824.

22.

et al., Demethylzelasteral inhibits proliferation and EMT via repressing Wnt/beta-catenin signaling in esophageal squamous cell carcinoma, J Cancer 12(13) (2021), 3967–3975.

23.

Lin

D.C.

et al., Genomic and molecular characterization of esophageal squamous cell carcinoma, Nat Genet 46(5) (2014), 467–73.

24.

Morgan

and Grandis

J.R.

, ErbB receptors in the biology and pathology of the aerodigestive tract, Exp Cell Res 315(4) (2009), 572–82.

25.

Yang

et al., MiR-106b-5p regulates esophageal squamous cell carcinoma progression by binding to HPGD, BMC Cancer 22(1) (2022), 308.

26.

Zhang

J.-X.

et al., Prognostic and predictive value of a microRNA signature in stage II colon cancer: a microRNA expression analysis, The Lancet Oncology 14(13) (2013), 1295–1306.

27.

and Zhang

, Epigenetic regulation on the gene expression signature in esophagus adenocarcinoma, Pathology – Research and Practice 213(2) (2017), 83–88.

28.

Gruszka

and Zakrzewska

, The Oncogenic Relevance of miR-17-92 Cluster and Its Paralogous miR-106b-25 and miR-106a-363 Clusters in Brain Tumors, Int J Mol Sci 19(3) (2018).

29.

Guarnieri

A.L.

et al., The miR-106b-25 cluster mediates breast tumor initiation through activation of NOTCH1 via direct repression of NEDD4L, Oncogene 37(28) (2018), 3879–3893.

30.

et al., Role of the miR-106b-25 microRNA cluster in hepatocellular carcinoma, Cancer Sci 100(7) (2009), 1234–42.

31.

Nagase

et al., Kazusa mammalian cDNA resources: towards functional characterization of KIAA gene products, Briefings in Functional Genomics 5(1) (2006), 4–7.

32.

Gudmundsdottir

et al., Integrative network analysis highlights biological processes underlying GLP-1 stimulated insulin secretion: A DIRECT study, PLoS One 13(1) (2018), e0189886.

33.

Miyoshi

et al., A microRNA-based liquid biopsy signature for the early detection of esophageal squamous cell carcinoma: a retrospective, prospective and multicenter study, Mol Cancer 21(1) (2022), 44.

34.

Fakhrjou

et al., Prognostic Value of Tumor-Infiltrating Mast Cells in Outcome of Patients with Esophagus Squamous Cell Carcinoma, Journal of Gastrointestinal Cancer 45(1) (2014), 48–53.

35.