Using m5C-related long non-coding RNAs to predict outcomes and treatment responses in esophageal squamous cell cancer

We have included patients diagnosed with ESCC in The Cancer Genome Atlas (TCGA) program (https://portal.gdc.cancer.gov/repository?facetTab=cases↗). Patients without complete clinical information and transcriptome profiling, or diagnosed with esophageal adenocarcinomas were excluded. Finally, the transcriptome, clinicopathologic, and somatic mutation data of 80 ESCC and 11 adjacent normal tissues were downloaded. Additionally, RNA microarray profiles and corresponding clinical information of 60 ESCC patients were downloaded from the Gene Expression Omnibus database (GSE53622↗, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=gse53622↗). The clinicopathological parameters of TCGA-ESCC and GSE53622↗ cohorts were summarized in Table S1. Immunohistochemical staining images of normal esophageal tissues were obtained from the Human Protein Atlas (HPA) (https://www.proteinatlas.org/↗). And this study started in February 2022 and finished in July 2022.

Paraffin embedded sections of 54 ESCC patients were obtained from the First Affiliated Hospital of Xi’an JiaoTong University. Tissues were collected during surgery and were used for Fluorescence in situ hybridization (FISH) examination. And 14 ESCC and corresponding normal tissue samples were collected for the detection of lncRNA expression.

We identified 16 m⁵C regulators from previous literature and extracted their expression from RNA-seq profiles of ESCC and adjacent normal tissues. Then, the differential expression of 16 m⁵C regulators was determined in the ESCC tissues versus adjacent normal tissues. The differential expression of these regulators was analyzed, and their interrelationships were visualized using the ‘corrplot’ R package. A protein–protein interaction (PPI) network of m⁵C regulators was constructed using the STRING database (https://cn.string-db.org/↗) with the gene interaction score ≥ 0.5²³. We selected the lncRNAs existed in both TCGA-ESCC and GSE53622↗ cohorts for widespread use of m⁵C-LPS. Pearson correlation coefficient was calculated between the expression of 16 m⁵C regulators and lncRNAs using the built-in function ‘cor.test’ in R. We identified 4279 m⁵C-related lncRNAs with |correlation coefficient| > 0.35 and the p-value < 0.01 for further analysis.

We used the ‘survival’ R package to perform univariate Cox regression analysis on the candidate m⁵C-related lncRNAs, filtering out those with significant prognostic value (p < 0.05). Then, we used the least absolute shrinkage and selection operator (LASSO) regression analysis with the ‘glmnet’ R package to establish a prognostic signature and calculate the coefficients for each lncRNA^24–26. These coefficients were used to generate a risk score formula: RiskScore=∑icoefficientofm5CrelatedlncRNAi×lncRNAiexpressionlevel. Patients were stratified into high- and low-risk groups based on their calculated risk scores. Kaplan–Meier (K–M) analysis was performed to assess the overall survival (OS) of different groups, and time-dependent receiver operating characteristic (ROC) curve analysis was used to evaluate the predictive value of the risk score with the ‘survivalROC’ R package.

A Sankey diagram was used to illustrate the one-to-one match between the m⁵C genes, m⁵C-related lncRNAs, and the corresponding risk types. Furthermore, a correlation circle graph was generated using the ‘corrplot’ and ‘circlize’ R package to visualize the co-expression status of the 9 identified lncRNAs. We also investigated the association between m⁵C-LPS and clinicopathological parameters. Both univariate and multivariate Cox regression analyses were conducted to investigate the independent value of the m⁵C-LPS and other parameters. Based on the significant prognostic variables, we constructed a nomogram to predict 1-, 2-, and 3-year survival rates using the ‘rms’ R package. And calibration curves were used to verify the agreement between nomogram-predicted survival and actual survival probabilities. Additionally, we evaluated the prognostic value of clinicopathological features by using ROC curves and calculating the area under the curve (AUC).

We conducted functional enrichment analyses based on gene ontology (GO)²⁷, Kyoto Encyclopedia of Genes and Genomes (KEGG)²⁸, and Reactome²⁹ databases to explore the biological functions and pathways associated with m⁵C-LPS through ‘clusterProfiler’³⁰ and ‘ReactomePA’ R package.

Estimate³¹, single-sample gene set enrichment analysis (ssGSEA)³², Cibersort³³, and xCell³⁴ algorithms were utilized to estimate the relative abundance of immune and stromal cells in the tumor microenvironment. We also calculated the Pearson coefficients between risk scores and immune checkpoint genes and immunomodulators, such as chemokines, receptors, MHC, immunoinhibitors, and immunostimulators, which were obtained from the TISIDB database³⁵.

The ‘maftools’ R package was utilized to identify the top 20 mutated genes based on the mutation rate across low- and high-risk groups³⁶. Subsequently, we further investigated the fraction of affected samples and pathways based on alterations in 10 canonical oncogenic signaling pathways for different risk groups³⁷.

We employed the ‘pRRophetic’ R package³⁸ to predict the half-maximal inhibitory concentration (IC50) for each patient using three publicly available drug sensitivity databases (Cancer Genome Project (CGP)³⁹, Cancer Therapeutics Response Portal (CTRP)⁴⁰, and Genomics of Drug Sensitivity in Cancer (GDSC)⁴¹). Additionally, we utilized the genomic-adjusted radiation dose (GARD) model⁴² to predict the radiotherapy response of each patient, with higher GARD values indicating increased sensitivity to radiotherapy.

The human normal esophageal cell line HET-1A was purchased from American Type Culture Collection (ATCC, Virginia, USA), and the human ESCC cell lines TE-1 and KYSE150 were purchased from the Cell Bank of the Chinese Academy of Sciences Typical Culture Preservation Committee (Shanghai, China). HET-1A was cultured in Dulbecco’s Modified Eagle’s Medium (DMEM, Gibico, USA) supplemented with 10% fetal bovine serum (FBS, Gibco, USA), while TE-1, and KYSE150 were cultured in Roswell Park Memorial Institute (RPMI) 1640 medium (Gibico, USA) supplemented with 10% FBS. All cells were cultured in a 5% CO₂ incubator at 37 °C.

Total RNA was extracted using the RNAfast200 kit (Fastagen, China) according to the manufacturer’s instructions. RNA concentration was quantified using NanoDrop 3000 (ThermoFisher, USA). Then, 1.0 μg of total RNA in a 20 μl reaction system was reversely transcribed into cDNAs using Evo M-MLV RT Kit with gDNA Clean for polymerase chain reaction (PCR, Accurate Biotechnology, China). Quantitative real-time PCR (qRT–PCR) was performed using 2× RealStar Green Fast Mixture (GeneStar Technology, China). GAPDH expression was used as an internal reference. The relative expression level of lncRNAs was calculated using the 2^−ΔΔCT method. Each experiment was performed in triplicate. The primer sequences used in this study are listed in Table S2.

The FISH probe of lncRNA AC002091.2↗ was synthesized by Servicebio (Wuhan, China). Paraffin embedded sections were dewaxed, rehydrated, digested, and dehydrated with dimethylbenzene, graded ethanol, protease K. Then the FISH probe was added to the hybridization mixture and incubated overnight. Next, the section was washed in the dark with washing buffer containing saline sodium citrate and PBS. Sections were stained with DAPI for 10 min and then visualized by fluorescence microscope.

All statistical analyses were performed using R software (version 4.1.1) and GraphPad Prism (version 8.0, USA). The differences between the two groups were compared using student’s t-test, while one-way analysis of variance (ANOVA) was used for multiple groups. Fisher’s exact test was used to compare categorical variables. The correlation between two continuous variables was analyzed using Pearson’s test. p-value < 0.05 was considered statistically significant.

This study was approved by the Ethics Committee of The First Affiliated Hospital of Xi’an Jiaotong University (Approval Number: 2017-146).

We extracted the expression profiles of 16 m⁵C regulators in the TCGA-ESCC cohort and subsequently compared their expression levels between 80 ESCC tumor samples and 11 normal adjacent samples. Our analysis revealed that the expression of most genes, including DNMT3B, NOP2, DNMT1, ALYREF, NSUN2, NSUN5, TET2, TET3, DNMT3A, TET1, and YBX1 were significantly higher in ESCC tissues than in normal adjacent tissues (Fig. 1A). Moreover, the immunohistochemical staining images of normal esophageal tissues from HPA showed that 10 of 15 m⁵C regulators were not more than medium expression, while NSUN3 and NSUN5 were not detected (Fig. S2). To investigate the interrelationships among these 16 m⁵C regulators, we obtained a PPI network using the STRING database. After setting the minimum interaction score as 0.5, we identified the PPI network contains all m⁵C genes and 90 edges (Fig. 1B). And TRDMT1 was found to be the hub gene of the network with 11 edges (Table S3). In addition, the correlation analysis revealed significant positive correlations between TRDMT1 and other 8 m⁵C genes. Interestingly, all m⁵C regulators showed a general positive correlation, with DNMT1 exhibiting the highest correlation with ALYREF (r = 0.67) (Fig. 1C).

Subsequently, we performed Pearson analysis based on the lncRNAs and m⁵C regulators in TCGA-ESCC profiles, and a total of 4279 lncRNAs were significantly correlated with m⁵C regulators (|Pearson coefficient| > 0.35 and p < 0.01). After filtering lncRNAs with the sum expression < 0.01, univariate Cox regression analysis was conducted to further explore the m⁵C-related lncRNAs associated with prognosis. Finally, we identified 41 lncRNAs that were significantly associated with the OS of ESCC patients (Table S4).

To eliminate the collinearity of variables and minimize estimation variance, LASSO regression analysis was applied to establish a prognostic signature using the 41 aforementioned lncRNAs. Subsequently, an m⁵C-LPS comprising 9 lncRNAs was identified based on the optimal λ value (Fig. 2A,B). Subsequently, the risk score was calculated based on the coefficients of the nine identified lncRNAs and their corresponding expression levels, yielding a concordance index (C-index) of 0.83, indicating strong discriminatory power (Fig. 2C,D). Besides, the model exhibited a sensitivity of 0.880, specificity of 0.643, positive likelihood ratio of 2.464, negative likelihood ratio of 0.187, positive predictive value of 0.524, and negative predictive value of 0.923. The m⁵C-LPS formula was calculated as follows: RiskScore=-1.45488×ATP2B1-AS1+0.78504×LINC2057+-3.09357×UBAC2-AS1+0.09339×CAHM+-1.19951×AC64807.1+-2.22323×AC37459.3+0.87974×AC2091.2+-0.60497×AC6329.1+0.2869×AC9275.1.

Subsequently, we categorized the 80 ESCC patients into low- and high-risk groups based on the median risk score. And the vital status and expression levels of the corresponding 9 lncRNAs in the cohort from TCGA-ESCC have presented in Fig. 2E. K–M analysis revealed that the patients in the high-risk group had relatively poorer OS and disease-free survival (DFS) compared with the low-risk group (OS: p < 0.0001, DFS: p = 0.064, Fig. 2F,G). Moreover, time-dependent ROC curves implied that m⁵C-LPS exhibited a promising ability to predict prognosis in the TCGA-ESCC cohort (1-year AUC = 0.839, 2-year AUC = 0.919, 3-year AUC = 0.898; Fig. 2H).

To validate the prognostic value of m⁵C-LPS, we calculated risk scores for another 60 ESCC patients from the GSE53622↗ cohort using the same formula. ESCC patients were divided into low- and high-risk groups according to the medium value. The distribution of the risk score, survival status, and lncRNAs expression showed that patients with higher risk scores had shorter OS and higher mortality status (Fig. 2I). Consistent with the findings in the TCGA-ESCC cohort, patients in the high-risk group presented significantly poorer prognoses (p = 0.012, Fig. 2J). And the AUC of the m⁵C-LPS was 0.7 at 2 years, 0.715 at 3 years, and 0.79 at 4 years (Fig. 2K).

We examined the co-expression status and differential expression of the 9 m⁵C-related lncRNAs. The Sankey plot showed one-to-one matches between the 7 m⁵C genes (5 writers: DNMT1, NSUN3, NSUN5-7; 2 erasers: TET1-2) and the 9 lncRNAs used in constructing the m⁵C-LPS. Additionally, the Sankey plot also depicted the risk type of each lncRNA (risk lncRNAs: AC002091.2↗, AC009275.1↗, CAHM, and LINC02057.1; protect lncRNAs: AC0006329.1, AC037459.3↗, AC064807.1↗, ATP2B1-AS1, and UBAC2-AS1, Fig. 3A). Moreover, the correlation circle plot revealed a general positive correlation among these m⁵C-related lncRNAs, except for CHAM and AC037459.3↗ had negative relationship with AC009275.1↗ and LINC02057, and UBAC2-AS1 showed negative correlation with AC002091.2↗ (Fig. 3B). Then, we compared the expression levels of these lncRNAs in normal esophageal and ESCC samples and observed that 7 lncRNAs were upregulated and 1 lncRNA was downregulated in ESCC samples (Fig. 3C).

Subsequently, we performed qRT-PCR using normal esophageal cell line HET-1A and ESCC cell lines TE-1 and KYSE150. The boxplot revealed that the upregulation of LINC02057, UBAC2-AS1, CAHM, AC002091.2↗, AC006329.1↗, and AC009275.1↗ in ESCC cells, while AC037459.3↗ was downregulated in ESCC cells (Fig. S3A). And we also detected the expression of these lncRNAs in ESCC and adjacent normal tissues and found that most of lncRNAs were upregulated in ESCC tissues (Fig. S3B). These results indicated that the expression patterns of m⁵C-related lncRNAs are consistent with the findings from the TCGA database. Since AC002091.2↗ was upregulated in ESCC cell lines and tissues and was of great prognostic value for ESCC patients, we subsequently investigated the relationship between the expression of AC002091.2↗ and patients’ survival. The FISH results showed that AC002091.2↗ was located in the cytoplasm (Fig. 3D). And K–M plot revealed that patients with higher AC002091.2↗ expression had relatively poor prognosis (Fig. 3E, p = 0.0028).

To evaluate the clinical significance of m⁵C-LPS, we assessed its association with various clinicopathological parameters of ESCC. Subgroup analysis stratified by T stage revealed a significantly higher risk score in T4 ESCC patients compared to T3 ESCC patients (p = 0.038, Fig. 4A). Stratification by M stage indicated an increased risk score in M1 patients, although the difference did not reach statistical significance (p = 0.15, Fig. 4B). No significant differences were observed between age, gender, race, tumor location, histologic grade, N stage, stage, reflux history and risk score (p > 0.05, Fig. S4A–G, Fig. 4C). Besides, ESCC patients with alcohol history exhibited a significantly elevated risk score than those without alcohol history (p = 0.022, Fig. 4D), while ESCC patients with or without smoking history had similar risk score (p = 0.71, Fig. S4H). In subgroup analysis stratified by adjuvant postoperative therapy, there was a trend towards a higher risk score in the pharmaceutical therapy and radiotherapy subgroup, although statistical significance was not achieved (p = 0.069 and 0.19, Fig. 4E,F). And ESCC patients with or without complete response after radiotherapy exhibited comparable risk scores (p = 0.47, Fig. S4I). Furthermore, ESCC patients with tumor presence, recurred/progressed, and deceased status had significantly increased risk scores (p < 0.05, Fig. 4G–I), consistent with previous results highlighting the value of m⁵C-LPS as a valuable prognostic marker.

Univariate and multivariate Cox regression analyses were conducted to determine the independent prognostic value of m⁵C-LPS and other clinicopathological parameters for ESCC patients. The forest plots showed that the N stage and risk score were independent factors for the poor prognosis (p < 0.05, Fig. 5A,B). Subsequently, we used time-dependent ROC curves to evaluate the prognostic potential of the risk score, age, gender, grade, stage, and TNM stage. The AUC values of the risk score were higher than those of other clinicopathological factors for 1-, 2-, and 3-year survival (Fig. 5C). These findings highlight the significant value of the risk score in predicting patient prognosis. Meanwhile, a nomogram was constructed based on the risk score and N stage of each ESCC patient, which could be a quantitative tool to predict 1-, 2-, and 3-year survival probability (Fig. 5D). Moreover, the calibration curves showed partial agreement between the predicted and observed survival probabilities (Fig. 5E).

We further investigated the relationship between the immune microenvironment and the risk score obtained from m⁵C-LPS. The relative abundance of immune and stromal cells of each sample was estimated using Estimate, Cibersort, ssGSEA, and xCELL algorithms. The heatmap revealed the different distribution patterns of various cell types between the low- and high-risk groups (Fig. 6A). Comparison of the Cibersort results revealed significant enrichments of CD8⁺ T cells, memory activated CD4⁺ T cells, and T follicular helper cells in the low-risk group, while M2 macrophages were found to be enriched in the high-risk group (Fig. 6B). The ssGSEA results showed that central memory CD8⁺ T cell, gamma delta T cell, macrophage, NK cell, plasmacytoid dendritic cell, Tregs, and T follicular helper cell were significantly enriched in the high-risk group (p < 0.05, Fig. S5A). However, there were no significant differences in stromal cells between the low- and high-risk groups (Fig. S5B). The correlation heatmap identified three main clustering modules: function immune cells, resting immune cells, and stromal cells (Fig. S6A). Furthermore, the correlation coefficient indicated a negative association between the risk scores and multiple well-known immune checkpoint molecules, except for IDO1 (Fig. 6C). The histogram and heatmap revealed inverse relationships between the risk score and most immunomodulators, including chemokines, receptors, MHC, immunoinhibitors, and immunostimulators (Fig. 6D, Fig. S6B). These findings indicated that the activation of immune components in the tumor microenvironment may contribute to better outcomes for patients in the low-risk group.

To elucidate the biological functions of the differentially expressed genes associated with m⁵C-LPS, we performed GO, KEGG, and Reactome enrichment analyses. The prominent GO terms in molecular function (MF), cellular component (CC), and biological process (BP) were catalytic activity acting on RNA, nuclear speck, and ncRNA metabolic process, respectively (Fig. S7A). Furthermore, the top five enriched KEGG terms included spliceosome, cell cycle, ribosome biogenesis in eukaryotes, RNA degradation, and Homologous recombination (Fig. S7B). The three main key modules identified in the Reactome analysis were rRNA processing, mRNA splicing and processing, and DNA damage repair response (Fig. S7C).

By analyzing the MuTect2 mutation annotation files, we identified the top 20 most frequently mutated genes in the low- and high-risk groups, as illustrated in Fig. S8A,B, respectively. The waterfall plots revealed that TP53, TTN, and KMT2D were most frequently mutated in both groups. However, the ranking of mutated genes showed slight changes between the two groups. For example, the mutation frequency of MUC16 was ranked third in the high-risk group (20%), but it dropped out of the top 20 mutated genes in the low-risk group. Furthermore, the mutation rates of seven oncogenic pathways (NOTCH, WNT, PI3K, MYC, TP53, TGF-beta, Cell-Cycle) were higher in the high-risk group compared to the low-risk group (Fig. S8C,D). These findings suggested that ESCC patients in the low- and high-risk groups may have different mutation driver genes and pathways.

Given that chemotherapy and radiotherapy are crucial in ESCC treatments, and DNA damage repair response plays a pivotal role in regulating chemoradiotherapy response, we attempted to evaluate the therapeutic response of the low- and high-risk groups. We estimated the IC50 levels of several commonly used chemotherapeutic drugs in each patient using CGP, CTRP, and GDSC-derived drug response data. The heatmap showed that the estimated IC50 levels of these drugs were reduced in the low-risk group, indicating that patients in low-risk group were more sensitive to chemotherapy (Fig. 7A). Boxplots further demonstrated that patients in the low-risk group exhibited greater sensitivity to five CGP-derived compounds (5-fluorouracil, cisplatin, docetaxel, vinorelbine, and etoposide), two CTRP-derived compounds (docetaxel and gemcitabine), and four GDSC-derived compounds (docetaxel, paclitaxel, oxaliplatin, and vinorelbine). And significant differences in the IC50 level of docetaxel were observed among three database-derived results (Fig. 7B–D). Besides, the radiation-sensitivity index (RSI) increased in the high-risk group, suggesting that patients in the high-risk group might require a higher radiotherapy dose, although there was no statistical significance (p = 0.23, Fig. 7E).

The data presented in this study can be found in TCGA (https://portal.gdc.cancer.gov/repository?facetTab=cases↗) and GSE53622↗ (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=gse53622↗) databases.