Development of a novel miR-3648-related gene signature as a prognostic biomarker in esophageal adenocarcinoma.

BACKGROUND: Esophageal adenocarcinoma (EA) is a typical immunogenic malignant tumor with a dismal 5-year survival rate lower than 20%. Although miRNA-3648 (miR-3648) is expressed abnormally in EA, its impact on the tumor immune microenvironment remains unknown. In this study, we sought to identify immune-related genes (IRGs) that are targeted by miR-3648 and develop an EA multigene signature.
METHODS: The gene expression data of 87 EA tumor samples and 67 normal tissue samples from The Cancer Genome Atlas (TCGA) database and the Genotype-Tissue Expression (GTEx) database were downloaded, respectively. Weighted gene co-expression network analysis (WGCNA), the CIBERSORT algorithm, and Cox regression analysis were applied to identify IRGs and to construct a prognostic signature and nomogram.
RESULTS: MiR-3648 was expectedly highly expressed in EA tumor tissues (P=2.6e-8), and related to the infiltration of activated natural killer cells (NK cells) and activated CD4 T lymphocytes (CD4 cells). A total of 70 miR-3648-targeted genes related to immune cell infiltration were identified. Among them, 4 genes (C10orf55, DLL4, PANX2, and NKAIN1) were closely related to overall survival (OS), and were thus selected to construct a 4-gene risk score (RS). The RS had a superior capability to predict OS [area under the curve (AUC) =0.740 for 1 year; AUC =0.717 for 3 years; AUC =0.622 for 5 years]. A higher score was indicative of a poorer prognosis than a lower score [hazard ratio (HR) =2.71; 95% confidence interval (CI): 1.45-5.09; P=0.002]. Furthermore, the nomogram formed by combining the RS and the TNM classification of malignant tumors (TNM stage) improved the accuracy of survival prediction [Harrell's concordance index (C-index) =0.698].
CONCLUSIONS: MiR-3648 may play a critical role in EA pathogenesis. The novel 4-gene signature may serve as a prognostic tool to manage patients with EA. 2021 Annals of Translational Medicine. All rights reserved.

Introduction

Esophageal adenocarcinoma (EA) is an aggressive malignancy with high mortality (1). Patients with EA are often diagnosed at an advanced stage, which results in less than half of the patients being eligible for potentially curative treatment at diagnosis (2). Treatment is largely reliant upon minimally effective cytotoxic chemotherapy, with the 5-year overall survival (OS) rate for patients with advanced EA being lower than 20% (3). There are significant prognostic discrepancies in patients with the same TNM classification of malignant tumors (TNM stage) (4). To improve risk prediction, it is necessary to identify objective biomarkers to predict response to EA treatment and to predict survival.

Tumor cells can restrain anti-tumor immunity by forming an immunosuppressive microenvironment (5). Cells around the tumor cells called the tumor microenvironment (TME) are also involved in tumor progression and therapy response (6). Different populations of immune cells in TME can also affect the anti-tumor response, which becomes a major obstacle to the treatment of tumors (7). The tumor immune microenvironment of EA remains poorly understood and an improved understanding could be crucial to improving EA prognosis and treatment.

MicroRNAs (miRNAs) regulate many important cellular processes such as cell proliferation, migration, and apoptosis (8). Several cancer-associated miRNAs have now been identified, of which miRNA-3648 (miR-3648) is one of the most established and broad-acting (9-11). MiR-3648 is abnormally expressed in many cancers, including esophageal carcinoma (11). MiR-3648 could regulate cell growth by effecting APC2 levels, which is a tumor suppressor and a negative regulator of Wnt signaling (9). Dysregulation of miR-3648 can promote cell proliferation and induce the invasion and metastasis of cancer cells (10,12). However, the impact of miR-3648 on clinical characteristics and the immune microenvironment in EA patients is less well explored. A more in-depth understanding of the immune regulatory network regulated by miR-3648 may provide better prognostic information for EA.

In the present study, we utilized weighted gene co-expression network analysis (WGCNA) to construct a co-expression network of miR-3648-targeted genes. The infiltrated immune cells were evaluated with the CIBERSORT algorithm, and immune-related genes (IRGs) in the co-expression network were determined. A Cox analysis was then performed to select genes that were closely related to patients’ survival, and a multigene prognostic signature was constructed. Finally, we developed a nomogram to improve the ability to predict patient survival. We demonstrated that miR-3648 could play a critical role in the development of EA, and constructed a novel 4-gene signature that target by miR-3648 could serve as a prognosis biomarker in EA. We present the following article in accordance with the TRIPOD reporting checklist (available at https://dx.doi.org/10.21037/atm-21-6237).

Methods

Source of data

The patients used for this study included 12 women and 75 men aged 27 to 86 years old from The Cancer Genome Atlas (TCGA) database (https://cancergenome.nih.gov/). Clinical information and expression profiles of miRNAs/mRNAs associated with EA were retrieved from TCGA database (see available online: https://cdn.amegroups.cn/static/public/atm-21-6237-1.xlsx). The outcome was defined as any death recorded in TCGA clinical data. MiRNA data were extracted from 87 EA tumor tissues and 13 normal esophageal tissues. The mRNA-seq expression profiles of 54 normal esophageal tissues as controls were downloaded from the Genotype-Tissue Expression (GTEx) database (https://www.gtexportal.org/). Finally, the expression data of 19835 genes in 87 tumor samples and 67 normal samples were summarized as mRNA data. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Differential expression analysis and target gene prediction

The R package limma was used to identify mRNAs that differed in expression between tumor and normal samples (13). For differentially expressed mRNA identification, we used a false discovery rate (FDR)-adjusted P value of <0.05 with a simultaneous absolute value of >1 for logFC.

To identify the target genes of miR-3648, 12 databases (Microt4, miRWalk, mir-bridge, miRanda, miRDB, miRMap, Pictar2, PITA, MiRNAMap, RNAhybrid, RNA22, and TargetScan) in miRWalk2.0 (http://zmf.umm.uni-heidelberg.de/apps/zmf/mirwalk2/) were employed. Genes that appeared in at least 3 databases were identified as the target genes of miR-3648.

A differentially expressed target gene (DETG) of miR-3648 was defined as any gene that was cross-linked to the target genes of miR-3648 or the differentially expressed genes (DEGs) in EA.

Evaluation of immune cell infiltration

The deconvolution algorithm CIBERSORT (https://cibersort.stanford.edu/) was applied with 100 permutations and leukocyte signature matrix 22 (LM22) to quantify the proportions of 22 infiltrating immune cell types in tumor and normal samples (14). CIBERSORT P values reflect the statistical significance of the results, and in this study unmatched samples that had P>0.05 were removed.

WGCNA

The R package “weighted gene co-expression network analysis” (15) was utilized to identify modules of highly correlated DETGs in EA tumor samples as previously described (16). Briefly, paired Pearson correlation was used to evaluate the weighted co-expression relationships among the subjects from all datasets in the adjacency matrix. A soft threshold of β=4 was selected, and the co-expression similarity matrix was changed to the adjacency matrix. Next, a topological overlap matrix (TOM) was built based on the adjacency matrix. Different modules were identified from TOM-based dissimilarity as densely interconnected gene clusters. Dynamic hybrid tree-cutting was then used to create a hierarchical clustering tree for dividing modules.

The infiltration fractions of immune cells in EA tumor tissues were employed as clinical traits. Pearson’s correlation test was then used to evaluate the module correlations with clinical traits. Among these clinical traits, the genes of the most significant module were determined to be IRGs for subsequent analysis.

Identification of candidate prognostic genes

To identify key prognostic genes, a Cox regression analysis was applied to estimate the correlation between the expression of the genes of the most significant module and the OS time of patients. Finally, four statistically significant genes for which P was <0.05 were identified as candidate prognostic genes.

Construction of the prognostic signature

A prognostic signature was constructed based on the expression levels of the candidate genes. The risk score (RS) for each tumor sample was calculated using the following formula:

In the above formula, Coef represents the coefficient of Cox regression of gene i. Exp represents the expression level of gene i, and n represents the number of genes included in the prognostic signature. To determine whether this RS could independently predict the prognosis of EA patients, univariate and multivariate Cox regression analyses were conducted. The median RS was then used as a cut-off for deciding high- or low-risk patients.

Establishment of the nomogram

A nomogram was established by incorporating all significant prognostic characteristics and the RS to estimate the probability of the 1-, 3-, and 5-year OS of EA. The calibration curve was visualized to measure the prediction performance of the nomogram. We also calculated Harrell’s concordance index (C-index) to assess the survival prediction capability of the nomogram. R package rms v6.0-2 was employed to construct the nomogram.

Statistical analysis

The area under the curve (AUC) of the receiver operating characteristic (ROC) curve and Kaplan-Meier (K-M) curve analyses were performed to evaluate the prognostic performance and diagnostic capacity. The Wilcoxon test was applied to estimate the differences between groups. All R program packages for statistical analysis were implemented on the R program (v4.0.4). P<0.05 was considered statistically significant.

Results

Upregulation of miR-3648 and its target genes in EA

MiR-3648 was expectedly more highly expressed in tumor tissues than in normal tissues (P<0.001; ) and showed excellent discrimination between tumor tissue and normal tissue (AUC =0.935; P=4.03e-15; ). We also compared the expression of miR-3648 between different clinical feature groups, but no significant difference was observed (). Based on the 12 databases of miRWalk2.0, 12,289 potential target genes of miR-3648 were identified. We identified 608 genes in at least 3 databases as the targets of miR-3648.

Figure 1

Clinical performance of miRNAs. (A) Expression level of miR-3648 in EA tumor/normal tissues; (B) ROC curve of miR-3648 in diagnosis; (C-H) expression level of miR-3648 in sex, age, the TNM classification of malignant tumors (TNM stage), and classification groups. MiRNAs, microRNAs; miR-3648, miRNA-3648; EA, esophageal adenocarcinoma; ROC, receiver operating characteristic.

A total of 5,746 DEGs were identified, including 2,817 upregulated genes and 2,929 downregulated genes (). Furthermore, we cross-linked the 608 target genes of miR-3648 and DEGs, and identified 170 DETGs ().

Figure 2

MiR-3648-targeted gene identification. (A) Volcano plot of DEGs of EA; (B) Venn diagram of miR-3648-targeted genes in DEGs. MiR-3648, miRNA-3648; DEGs, differentially expressed genes; EA, esophageal adenocarcinoma; TCGA, The Cancer Genome Atlas.

The landscape of infiltrating immune cells of EA

We applied the CIBERSORT algorithm to calculate the infiltration fraction of 22 immune cell types. Our results showed that M0 macrophages, M1 macrophages, activated dendritic cells, activated mast cells, activated memory CD4 T cells, follicular helper T cells, and regulatory T cells (Tregs) were the most common immune cell populations in EA tumor tissue. By contrast, resting dendritic cells, resting mast cells, monocytes, activated natural killer cells (NK cells), resting memory CD4 T cells, and naïve CD4 T cells were the most common immune cell populations in normal tissue (). We further evaluated the associations among the expression of miR-3648 and the infiltrating immune cells. The high expression group of miR-3648 was found to promote more infiltrating NK cells as well as memory CD4 T cells ().

Figure 3

The profiles of immune infiltration. (A) The profiles of immune infiltration of EA; (B) the relationship between the immune infiltration and the expression of miR-3648. ns, P>0.05; *, P≤0.05; ***, P≤0.001; ****, P≤0.0001. EA, esophageal adenocarcinoma; miR-3648, miRNA-3648.

Identification of IRGs

The 170 DETGs were then considered for the WGCNA. We selected β=4 to satisfy the scale-free network law. In the network analysis, 2 gene modules (denoted “turqoise” and “blue”) and 37 oligogenes were identified under the parameter settings minModuleSize =20 and mergeCutHeight =0.25 (). The turquoise module contained 70 genes, and the blue module contained 63 genes. Subsequently, the infiltration fraction of the 13 immune cells in samples were examined as a clinical trait in order to assess the associations between them and gene modules (). Turquoise exhibited the most significant negative correlation with Treg infiltration (cor =−0.48; P=3e-6). Activated mast cells and activated NK cells appeared to be negatively correlated with the turquoise module. Turquoise was also related to resting dendritic cells. The 70 genes in the turquoise module were therefore identified as IRGs.

Figure 4

WGCNA of the DETGs of miR-3648. (A) Two modules and 37 oligogenes were clustered. Each branch of the tree diagram represents genes, and genes clustered into the same module are assigned the same module color; (B) correlation analysis between modules and immune infiltration. P values <0.05 are colored red. WGCNA, weighted gene co-expression network analysis; DETGs, differentially expressed target genes; miR-3648, miRNA-3648.

Construction of the 4-gene signature

Through a Cox regression analysis, the prediction performance of 70 IRGs on patients’ OS was evaluated. Four genes (C10orf55, DLL4, PANX2, and PANX2) that were closely related to OS were identified (). Then, the RS was calculated by multiplying these 4 expression values by the Cox coefficients. The formula for calculating the RS was: RS = 0.29 × expression of C10orf55 + 0.39 × expression of DLL4 + 0.15 × expression of PANX2 + 0.16 × expression of NKAIN1.

Table 1

Cox analysis of 4 genes in the signature

Gene symbol	Gene full name	Coefficient	HR (95% CI)	P value
C10orf55	Chromosome 10 putative open reading frame 55	0.29	1.3 (1.1–1.7)	0.017
DLL4	Delta-like canonical Notch ligand 4	0.39	1.5 (1.1–2.0)	0.018
PANX2	Pannexin 2	0.15	1.2 (1.0–1.3)	0.030
NKAIN1	Na+/K+ transporting ATPase interacting 1	0.16	1.2 (1.0–1.4)	0.046

HR, hazard ratio; CI, confidence interval.

The median RS was 9.04. The 87 patients were separated into 43 high-risk patients and 44 low-risk patients. Low-risk patients had longer survival times than high-risk patients (P=0.0012; ). The AUCs of the RS for 1-, 3-, and 5-year survival prediction were 0.740, 0.717, and 0.622, respectively ().

Figure 5

Prognostic value of the RS. (A) K-M curves of high/low-RS groups with OS time; (B) ROC curves for the RS to predict 1-, 3-, and 5-year survival; (C) univariate and (D) multivariate Cox regression analysis of the RS and clinical characteristics. RS, risk score; K-M, Kaplan-Meier; OS, overall survival; ROC, receiver operating characteristic; AUC, area under the curve.

Independent prognostic analysis of the 4-gene signature

We used univariate Cox regression analysis to evaluate the prognostic performance of the RS and clinical characteristics (sex, age, TNM stage, and classification). The univariate Cox regression analysis results revealed that the RS [hazard ratio (HR) =2.71 in the high-risk group; P=0.002] and the TNM stage (stage III/IV; HR =3.63; P=0.002) were correlated with patients’ survival (). TNM stage was therefore considered for the multivariate Cox regression analysis to evaluate the prognostic performance of the RS. Multivariate Cox analysis suggested that the RS could be considered an independent prognostic marker [HR =2.29 in the high-risk group; 95% confidence interval (CI): 1.07–4.94; P=0.03] ().

Nomogram analysis

Based on the TNM stage and the RS, nomograms were constructed using the Cox regression model to predict the probability of death at 1, 3, and 5 years (). Based on the RS (C-index =0.661) and TNM stage (C-index =0.630), the nomograms proved to be good predictors. Our optimal survival prediction nomogram was built by combining the RS and TNM stage. According to the nomogram, OS had a C-index of 0.698.

Figure 6

Nomogram and calibration curves for survival prediction. (A) The nomogram for survival prediction model; (B-D) the calibration curves for survival prediction model. OS, overall survival.

Discussion

EA is a typical immunogenic malignant tumor. Because the TME in EA is highly heterogeneous, it is not predictable which patient will respond best to treatment. Therefore, it is essential to screen prognostic markers that are related to EA’s immune microenvironment. In this study, we found that miR-3648 was significantly upregulated in EA tumor tissue compared with normal tissue and identified 4 miR-3648-targeted genes related to the TME. A multigene signature for EA was then developed. This signature may be used as a potential prognostic biomarker in EA patients.

MiR-3468 is abnormally expressed and has been implicated in many kinds of cancer (10,11). A previous study suggested that miR-3648 can induce the proliferation of prostate cancer cells (12). Upregulated miR-3648 was reported to induce the invasion and metastasis of human bladder cancer (10). EA tumor tissue expresses a higher level of miR-3648 than normal tissue, indicating a positive clinical effect, which is in line with previous studies. Furthermore, the high expression of miR-3648 was related to the infiltration of activated NK cells and activated CD4 cells. This finding suggests that miR-3648 may play a crucial role in the progression of EA. Despite we demonstrated that miR-3648 was related to the outcome of EA patients, the mechanisms of miR-3648 in EA should be further explored in vitro experiments.

We applied the CIBERSORT algorithm to predict immune cell infiltration levels to better understand the complexity of the EA tumor immune microenvironment. We found that M0 macrophages, M1 macrophages, activated dendritic cells, activated mast cells, activated memory CD4 T cells, follicular helper T cells, and Tregs were the most common immune cell populations in EA tumor tissues. The IRGs targeted by miR-3648 most related to Tregs were further identified. Tregs are an immunosuppressive subset of CD4+ T cells and can suppress anticancer immunity, thus promoting tumor development and progression (17). The correlation of miR-3648 with the tumor immune microenvironment may be a prominent mechanism in the pathological process of EA.

Based on the 170 DEGs targeted by miR-3648, we identified 70 IRGs using the WGCNA algorithm. Four genes were identified as being associated with patient survival based on a Cox regression analysis. Previous reports indicate that these genes are associated with different types of cancer in various ways. For example, C10orf55 has been found to be a potential biomarker to predict the relapse of acute myeloid leukemia (18). DLL4 plays a key role in tumor angiogenesis, functioning as a negative regulator of tumor angiogenesis and is upregulated in the tumor vasculature (19,20). A previous study indicated that DLL4 could be an independent prognostic factor for predicting the OS of patients with non-small cell lung cancer and was correlated with immunocyte infiltration (21). PANX2 is a channel-forming glycoprotein and is associated with many common diseases, including cancers (22). PANX2 was reported as a growth regulator in glioma cells (23) and was associated with the survival of clear cell renal cell carcinoma patients (24). NKAIN1 was overexpressed in urinary sediments from patients with prostate cancer and could be a promising biomarker for the early diagnosis of prostate cancer (25). In our study, all these genes were suggested as potentially useful biomarkers for predicting survival outcomes in EA. Future, in vitro and in vivo experiments are needed to determine their role in the proliferation, metastasis, and invasion of tumor EA cells.

We constructed an IRG-based 4-gene (C10orf55, DLL4, PANX2, and NKAIN1) prognostic signature. By combining this signature with the TNM stage, a nomogram was developed which was extremely accurate in predicting a patient’s prognosis. The novel immune-related prognostic signature may serve as a prognostic biomarker for EA and may reflect changes in immune cell infiltration in the TME. It may also help in individualized diagnosis and treatment.

Some limitations should be noted here. First, the study included only DETGs of miR-3648, so the prognostic signature does not reflect all EA transcription. Second, the associations between the 4 genes and the biological mechanisms of EA were not clarified in this study. Finally, as a retrospective study, our results should be confirmed in prospective multicenter clinical trials.

In summary, using WGCNA and CIBERSORT, we evaluated the tumor immune microenvironment of EA versus normal tissue and identified 70 miR-3648-targeted genes involved in immune cell infiltration in EA. Subsequently, a novel 4-gene signature was identified (C10orf55, DLL4, PANX2, and NKAIN1). A nomogram based on the signature and the TNM stage showed promising predictive performance for EA survival and may enable individualized treatment.

The article’s supplementary files as