N6-Methyladenosine-Related Long Noncoding RNAs as Potential Prognosis Biomarkers for Endometrial Cancer

Introduction

Endometrial cancer (EC) is one of the most common diseases of the female reproductive system, ranking the sixth most common gynaecologic malignancy globally, with nearly 382,000 new diagnoses and 89,900 deaths in 2018.¹ Epidemiological analysis of late age at menopause, postmenopausal hormone replacement therapy, and early age at menarche may serve as risk factors for the development of EC.² Most cases are diagnosed due to abnormal introitus bleeding and the prognosis is good in the early diagnosed cases, while the prognosis is poor in progressive cases, with a 5-year survival rate of only 15%.³ Therefore, new treatment strategies and molecular indicators are crucial for the prognosis of endometrial cancer.

N6-methyladenosine modification (m6A) is the most abundant and prevalent RNA modification present in eukaryotic cells and contains three categories of m6A methyltransferases: writer, reader, and and eraser.⁴ M6A is implicated in multiple aspects of cellular metabolism, including RNA splicing, transcription, translation, and degradation.⁵ Studies have shown that m6A participates in the development of a multitude of tumours, including EC. IGF2BP1 acts as an m6a reader to recognize and stabilize PEG10 mRNA and enhance its expression, thereby accelerating cell proliferation and promoting the progression of EC.⁶ FTO removes YTHDF2-mediated degradation of HOXB13 and increases its expression, thereby triggering the wnt pathway and facilitating invasion and metastasis of EC cells.⁷ YTHDF3 is involved in multiple aspects of breast cancer brain metastasis by binding to its mRNA, promoting its translation, and regulating downstream genes.⁸

lncRNAs are a class of non-coding RNAs with molecular mass greater than 200 and less than 2000, which were previously shown to be unable to encode proteins, however some have been found to have coding functions.⁹ lncRNAs participate in the pathogenesis of many diseases, such as Mendelian disorders, cardiovascular diseases, and tumours.¹⁰ lncRNA H19 promotes hepatocellular carcinoma bone metastasis by mediating PPP1CA/p38MAPK-dependent reduction in OPG expression.¹¹ LINC01969 is predominantly located in the cytoplasm and promotes the malignant biological behaviour of ovarian cancer by acting as a sponge of miR-144-5p adsorbed to regulate the expression of the downstream target gene LARP1.¹² Therefore, lncRNAs have essential roles in several aspects of tumour initiation and development.

The study of the tumour microenvironment (TME) has become increasingly prevalent. The tumour microenvironment consists of soluble molecules, tumour stromal cells, and extracellular matrix, of which the most investigated are immune cells and stromal cells.^13–15 Several studies have indicated that the tumour microenvironment may be involved in multiple aspects of tumour cell invasion, metastasis, drug resistance, and immune escape.^16,17 Furthermore, the study of immune-related lncRNAs has attracted increasing attention. lncRNA Xist modulated the expression of C/EBPa and KLF6 through competitive binding of miR-101 to mediate macrophage polarisation and thus affect the proliferation and migration of breast and ovarian cancer cells.¹⁸ lncRNA LUCAT1 inhibits the transcription of interferon genes by interacting with STAT1 in the nucleus, thereby inhibiting the immune response of human cells.¹⁹ lncRNA AFAP1-AS1 acts on KYSE410 cells through exosomal delivery of M2 type macrophages, down-regulating miR-26a, and promoting the expression of ATF2, thereby promoting cancer cell migration, invasion, and lung metastasis.²⁰

Currently, lncRNA studies are more common in EC, whereas m6A studies are scarce. In this study, we extracted three m6A-related lncRNAs based on data from EC patients in The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/) using bioinformatics and statistical analysis, and constructed prognostic models for m6A-associated lncRNAs to analyse the prognostic-related risk factors of EC. Endometrial patients were divided into high-and low-risk groups, and survival analysis and immune microenvironment analysis were performed. Finally, an accurate nomogram was constructed to validate and perform KEGG and GO enrichment analyses. These analyses suggest that these lncRNAs may serve as biomarkers for the diagnosis and prognosis of patients with EC.

Materials and Methods

Databases and m6A-Related Genes

Fragments per kilobase of transcript per million mapped reads (FPKM) data for EC were downloaded from the TCGA public database (https://portal.gdc.cancer.gov/),²¹ which comprised 23 normal samples and 552 tumour samples. The clinical database of EC patients was downloaded from UCSC Xena (https://xenabrowser.net/).^21–23 Then, the Perl software was used for data processing to generate a matrix of all differential genes in EC, which was then divided into lncRNA matrix and mRNA matrix. Additionally, 21 m6A-related genes were identified based on the current study, including readers (YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, IGF2BP1, IGF2BP2, IGF2BP3, RBMX, HNRNPA2B1, and HNRNPC), writers (METTL3, METTL14, METTL16, RBM15, RBM15B, WTAP, VIRMA (KIA1499), and ZC3H13), and erasers (FTO and ALKBH5). By analysing the mRNA expression matrix and m6A-related genes, the expression matrix of m6A-related genes was extracted. The m6A-related lncRNA expression matrix and m6A-related lncRNA analysis were performed using the limma package of R software and BiocManage (corFilter=0.6, pvalueFilter=0.001).

Survival Analysis

Using the Perl programming language, the expression matrix of m6A-related lncRNAs was organised and analysed with reference to the clinical data of EC. Using the survival package, univariate and multivariate Cox regression analyses were performed on the clinical data of m6A-associated lncRNAs to obtain a series of m6A-associated lncRNAs with P-value, HR values, and risk scores. The risk score calculating formula is: Risk score = , where means the coefficients, is the FPKM value of each m6A-related lncRNAs.²⁴ Patients were divided into high and low-risk groups based on the risk scores, and survival analysis and heat mapping were conducted for the high and low-risk groups using the survivor, “survmine” and “pheatmap” packages in R software. Furthermore, survival analysis was performed individually on the m6A-associated lncRNAs associated with the construction model, and survival curves were plotted. To verify the independence of the high- and low-risk group model construction, all risk-related genes and m6A-associated lncRNAs were subjected to principal component analysis using the Limma and scatterplot3d packages in R software.

Independent Prognostic Analysis of m6A-Associated lncRNA and Validation of Clinical Characteristics

To assess the impact of certain clinical characteristics on the prognosis of m6A-related lncRNA, univariate independent prognostic analyses were performed in terms of age, pathological grade, and FIGO staging by using the “survival” packages in R software, respectively. To validate the accuracy of the model, the multivariate independent prognostic analysis was performed and ROC curves were plotted for the above clinical characteristics with the “survROC” package. To authenticate the role of the clinical characteristics on prognosis from the different perspectives, we used multivariate COX regression to plot nomogram by using the “rms”, “foreign”, and “survival” packages in R software.

Tumour Immune Microenvironment Analysis

The immune cells and stromal cells in the tumour microenvironment (TME) were scored according to the ESTIMATE algorithm using the limma package and the ESTIMATE package of R language. The tumour microenvironment scores and clinical data were collated and combined for analysis using Perl programs to obtain immune and stromal scores in different risk states. By utilising the R software, box plots were drafted based on the immune score and stromal scores of patients in the high and low-risk group. The ESTIMATE algorithm in the R software ESTIMATE package was used to estimate the ratio of the immune matrix component in the TME for each sample, presented in the form of three scores: ESTIMATE score, stromal score, and immune score. This means that the higher the corresponding score, the greater the ratio of the corresponding component in the TME. Analysing and collating the risk scores of m6A-related genes and immune cells, correlation analysis was performed using the scales, “ggplot2” and “ggtext” packages in R software to plot bubble plots. The limma package, corrplot package, and BioManager package in R software were used for co-expression analysis of all differential gene matrices of endometrial EC with prognosis-associated m6A-related lncRNA matrices, where the asterisks (*) represents statistically significant differences.

Functional Annotation

The expression matrix of all genes and clinical data of different risk states were collated, followed by gene enrichment analysis (GO analysis) using the GSEA (4.0.2) software.

Samples and Quantitative Real-Time Polymerase Chain Reaction (qRT-PCR)

We collected 20 EC tissue and 11 normal endometrial tissue specimens. These specimens were obtained from surgical patients at Wuhan Union Hospital and were performed in strict accordance with the relevant regulations of medical ethics (Ethics Committee of Tongji Medical College of Huazhong University of Science and Technology), and informed consent was obtained from the participating patients.

To further verify the expression levels of m6A-related lncRNAs in the tissue samples, we extracted total RNA using TRIzol (RNAiso Plus, Takara, Japan) and cDNA was obtained using a reverse transcription kit (ABclonal, Wuhan, China). Expression of three m6A-related lncRNAs was calculated using 2-CT values, and GAPDH mRNA expression was used as a reference. The primer sequences used in this study were as follows: lncSCARNA9 forward 5ʹ-AGTCTTTCCAGTCTACCTGATGC-3ʹ and reverse 5ʹ-CATTGCCCAGAAATGATTAGGCT-3ʹ; lncTRAF3IP2 -AS1 forward 5ʹ-ACTGGTTGACAGAGCACCAAC-3ʹ, and reverse 5ʹ-AAATCCCATCCGTCCTTGCCT; lncAL133243.2 forward 5ʹ- AGTCCACCATTGCTCAACCGA-3ʹ; and reverse 5ʹ-ATCGGCCTTACATCTCCTGGC-3ʹ.

Bioinformatic Analysis

Perl programming language (Strawberry-Perl-5.30.0.1) was used to perform data collation and consolidation, and the databases were analysed using different packages in R software (R-3.6.1), such as limma, survival, ggplot2, and rms. The Mann–Whitney U-test was employed to explore the mRNA levels of m6A methylation-associated lncRNAs. Differences between the two subgroups were assessed with Student’s t-test. Survival curves were plotted using the Kaplan-Meier method. Univariate and multivariate analyses were performed with Cox regression models to explore the independent prognostic value of risk scores in combination with other clinical characteristics. The Q-PCR data were collated and analysed using GraphPad software. A value of p < 0.05 was considered statistically significant.

Results

Univariate and Multifactorial Cox Regression Analysis of Prognosis and Co-Expression Analysis of m6A-Related Genes with lncRNA

To determine whether m6A-associated lncRNAs are associated with the prognosis of EC patients, we identified six m6A-associated lncRNAs (TRAF3IP2-AS1, Z83843.1, AL133243.2, AL163051.2, Z98884.2, and SCARNA9) using univariate Cox regression analysis (Figure 1A and B). And the ROC curves (Supplementary Figure 1A) demonstrated that differential genes of EC harbored a promising ability to predict OS in the TCGA cohort (1-year AUC = 0.585, 3-year AUC = 0.639, 5-year OS = 0.644). SCARNA9 was negatively correlated with the prognosis of EC, while the remaining lncRNAs were positively correlated; subsequently, we performed co-expression analysis of m6A-related genes with lncRNAs which were differentially expressed in EC. Figure 1C demonstrates that the variably expressed lncRNAs in EC were mainly concentrated in three m6A methylation enzymes, RBM15, YTHDF3, and WTAP. Multifactorial Cox regression analysis identified three m6A-associated lncRNAs (TRAF3IP2-AS1, AL133243.2, and SCARNA9) associated with the prognosis of EC (Table 1). Finally, the Sankey diagram was plotted to visualize the co-expression relationship between the three m6A-related lncRNAs and m6A-related genes based on risk scores (Figure 1D).

Table 1 Multivariate Cox Regression Analysis of the m6A-Related lncRNAs

Figure 1 Univariate and multifactorial Cox regression analysis of prognosis and co-expression analysis of m6A-related genes with lncRNA. (A) Forest plot of m6A-related lncRNAs. HR>1(red)represents risk-associated lncRNAs, HR<1 (green) indicates protective lncRNA. (B) Heatmap of m6A-related prognostic lncRNAs. *p < 0.05, ***p < 0.001. (C) Network of the correlations between m6A-related genes and m6A-related prognostic lncRNAs. (D) Co-expression Sankey diagram. The left column represents m6A-related genes, the middle column represents m6A-related lncRNAs, and the right column represents risk types.

Survival Analysis of the Three m6A-Related Prognostic lncRNA

We classified patients with EC into high-risk and low-risk groups based on the median risk scores. The risk score formula was calculated as follows: risk score = (expression level of TRAF3IP2-AS1 multiplied by 0.695) + (expression level of AL133243.2 multiplied by 0.673) + (expression level of SCARNA9 multiplied by −0.346). To compare the overall survival (OS) rates of patients with different risk score statuses and to analyse the effect of different lncRNAs on the prognosis of patients, survival curves were plotted (Figure 2). As illustrated in Figure 2A, EC patients with lower risk scores had higher overall survival rates; however, the OS rates of the patients in the high-risk group were diametrically opposite. Figure 2B and C show the risk scores and the distribution of patient survival status associated with the different risk statuses and it shows that patients in the high-risk status group had higher risk scores and dead status. The heat map in Figure 2D illustrated the expression of three m6A-related (SCARNA9, TRAF3IP2-AS1, and AL133243.2) prognostic lncRNAs at different risk states. The expression of SCARNA9 was higher in the low-risk group than in the high-risk group, while the opposite trend was observed for TRAF3IP2-AS1 and AL133243.2. Then, we performed survival analysis of these three lncRNAs separately and found that low expression of TRAF3IP2-AS1 and AL133243.2, and high expression of SCARNA9 was associated with excellent prognosis (Figure 2E–G). To verify the difference in expression of m6A-related lncRNAs between EC tissue and normal endometrium, q-PCR was performed. The results (Figure 2H–J) suggested that the expression level of lncSCARNA9 was higher in normal endometrial tissues than in EC tissues, while the opposite trend was true for lncTRAF3IP2-AS1 and lncAL133243.2. The lncRNA expression levels were associated with the same prognostic trend.

Figure 2 Survival analysis of the three m6A-related prognostic lncRNA. (A) A survival curve based on risk scores of EC. (B and C) Risk score distribution and survival status of patients with EC. (D) Hierarchical clustering of three m6A-related lncRNAs expression levels based on risk scores. (E–G) KM curves illustrating the relationship between different expression levels of three m6A-related lncRNAs and overall survival. (H–J) The expression level of lncRNAs (SCARNA9, TRAF3IP2-AS1, and AL133243.2) were demonstrated by q-PCR. **p < 0.01.

Hierarchical Analysis of m6A-Related lncRNA

To verify the influence of some clinical characteristics such as age, pathological grade, and FIGO stage on the prognosis of patients with EC, we examined univariate and multivariate COX analyses and plotted multi-index ROC curves to verify their reliability. We demonstrated that age, pathological stage, and FIGO classification were all risk factors for the patient’s prognosis as shown in Figure 3A–C. M6A-related lncRNAs were significantly correlated with OS (HR,1.859; 95% CI:1.472–2.346, p value<0.001) by univariate Cox analysis and verified by multivariate Cox analyses. Moreover, the ROC curves suggested that the AUC was greater than 0.5, indicating a high degree of authenticity. Furthermore, we analysed the correlation of each m6A-related lncRNA with age, pathological stage, and FIGO classification separately; patients with a high expression level of AL133243.2, were typically accompanied by age >65 years, high pathological grade, and late FIGO staging. SCARNA9 expression was unrelated to age and FIGO stage, and negatively correlated with pathological grade. The expression of TRAF3IP2-AS1 increased in patients with more advanced FIGO staging and pathological grade (Figure 3D–F).

Figure 3 Hierarchical Analysis of m6A-related lncRNA. (A and B) Univariate and multivariate analyses revealed that risk score was an independent prognostic predictor. (C) Multi-index ROC curve for the risk score, age, grade and stage, accompanied by AUCs of 0.585,0.614, 0.652, and 0.792, respectively. (D–F) Patients with different clinicopathological features (including age, pathologic grade and FIGO stage) had different expression levels of the m6A-related lncRNAs. **p < 0.01, ***p < 0.001, ns, p > 0.05.

Principal Component Analysis and Construction of the m6A-Related Based Nomogram

The principal component analysis (PCA) at different levels was used to validate whether the risk status was distinguished among the different influencing factors. As shown in Figure 4A–D, the results suggested that the distributions of the high-risk and low-risk groups were relatively independent. Then, to establish a clinical quantitative tool applicable to the assessment of OS (one-year survival, two-year survival, and three-year survival) in patients with EC, we constructed a nomogram based on the m6A-related risk scores, FIGO staging, pathological grading, and age in the TCGA database (Figure 4E). The EC patients were randomly divided into the training (70%) and validation groups (30%). We then plotted ROC curves to verify the accuracy of the model. The outcomes revealed that the AUC values of the 3-year and 5-year survival rates were 0.793 and 0.816 for the training group (Figure 4F) and 0.721 and 0.769 for the validation group (Figure 4E), respectively, indicating the high stability of the nomogram model.

Figure 4 Principal Component Analysis and construction of the m6A-related based Nomogram. (A–D) Principal component analysis of low and high-risk groups based on all differential genes, m6A-related genes, m6A-related lncRNAs and risk-associated genes. (E) Nomogram based on risk score, age, pathologic grade and FIGO stage. (F and G) Time-dependent receiver operating characteristic (ROC) curves of training group and validation group for the nomogram (for predicting 3, and 5-year OS).

Tumour Immune Microenvironment Score and Immune Cell Infiltration Analysis

A correlation between m6A and tumour-infiltrating immune cells has recently been demonstrated. For instance, METTL3 activated T cells and enhanced the activation of the TLR4/NF-κB signalling pathway by mediating the methylation of TRL4, CD80, and CD40 signalling ligand Tirap transcripts, thus enhancing their translation in dendritic cells, which in turn stimulates cytokine production.¹¹ Persistent neoantigen-specific immunity is modulated by mRNA m6A methylation via the m6A-binding protein YTHDF15.¹² This research evaluated the correlation of m6A-related genes with tumour-infiltrating immune cells in EC. First, we assessed the stromal and immune scores of patients with EC based on risk status, and the results revealed that both the stromal and immune scores were lower in the high-risk group than in the low-risk group, which was statistically different (Figure 5A). Therefore, we assessed the relationship between different risk scores and tumour-infiltrating immune cells. We collected and collated the risk status of different immune cells from different algorithms (including XCELL, TIMER, QUANTISEQ, MCPCOUNTER, EPIC, CIBERSORT-ABS, and CIBERSORT), and the risk status of lncRNAs associated with m6A was collated and merged to obtain the bubble graph (Figure 5B). The results suggested that CD8+ T cells, NK cells, and M0 macrophages were negatively correlated in different databases, while M2 macrophages were positively correlated. To further clarify the correlation between prognosis-related lncRNAs and immune genes, we verified the co-expression of PD-L1 with six lncRNAs. The results indicated that PD-L1 was positively correlated with Z83843.1, AL133243.2, AL163051.2, Z98884.2, and SCARNA9, but not with TRAF3IP2-AS1 (Figure 5C). Additional immunoassay sites, such as PD-1, cytotoxic T lymphocyte associated antigen (CTLA-4), lymphocyte-activation gene 3 (LAG-3), T cell immunoglobulin domain and mucin domain-3 (TIM-3), T-cell immune receptor with immunoglobulin and ITIM domain (TIGIT), and killer cell immunoglobulin like receptor, three Ig domains and long cytoplasmic tail 1 (KIR3DL1), were verified and their results were consistent with the PD-L1 results described above (Supplementary Figure 1B).

Figure 5 Tumor immune microenvironment score and immune cell infiltration analysis. (A) ESTIMATE score, stromal score, and immune score of different groups were evaluated by ESTIMATE. (B) Immune cell infiltration analysis based on m6A-related lncRNAs risk scores. (C) Co-expression analysis of m6A-related prognostic lncRNAs and immune gene PD-L1.

Pathway and Process Enrichment Analysis

To investigate the differences in signalling pathways, molecular mechanisms, and biological processes of m6A-related lncRNAs between the high-and low-risk groups, we performed the KEGG enrichment analysis. The enrichment analysis (Figure 6) indicated that the high-risk group was mainly enriched in metabolic and RNA-related pathways, such as tyrosine metabolism, galactose metabolism, RNA polymerase, and DNA replication. Related genes enriched in the above pathways were provided in the Supplementary Table 1. These results may inform us about the biological pathways associated with the m6A-related lncRNAs.

Figure 6 Functional annotation. Gene set enrichment analysis (GSEA) indicated significant enrichment of m6A‐related phenotype in the high-risk patients.

Discussion

In the present study, 552 cases of EC from the TCGA database were analysed in order to investigate the relevance of m6A-related lncRNAs in the prognosis and immune microenvironment of EC. Univariate Cox regression analysis identified six m6A-associated lncRNAs affecting EC prognosis, and further multivariate Cox regression analysis yielded three lncRNAs, SCARNA9, TRAF3IP2-AS1, and AL133243.2. Moreover, the correlation of these three lncRNAs with clinical characteristics and expression levels was performed separately. Based on the median risk scores, we categorised patients with EC into high-risk and low-risk groups, and the overall survival rate of the high-risk group was lower than that of the low-risk group. Combining the age, pathological grading, FIGO stage, and risk score of the cases from the TCGA databases, we constructed a nomogram model and used the ROC curve to test its excellent accuracy and stability to predict the overall survival rate of EC patients. Immunisation scores of high-and low-risk groups were compared according to risk scores as variables, and the association of m6A-related lncRNAs with immune cell types and immune gene mutation loci was analysed using different algorithms.

M6A affects disease development by influencing the expression of target genes, and is involved in almost all processes of mRNA, including shearing, translation, degradation, export, and folding.^25,26 An increasing number of studies have shown that long non-coding RNAs are involved in tumorigenesis and progression as primary regulators at the level of the m6A epigenetic modifications.²⁷ Knockdown of ALKBH5 inhibits the malignant biological behaviour of colon cancer cells through NEAT1 demethylation.²⁸ Linc00460 binds to the 3ʹUTR region of HMGA1 through interaction with IGF2BP2 and DHX9 to enhance the stability of HMGA1 mRNA and induce proliferation, invasion, and migration of colon cancer cells.²⁹ In cervical cancer cells, IGF2BP3 enhances the stability of lnRNA KCNMB2-AS1 by binding to it, thereby promoting tumour cell proliferation.³⁰ M6A-modified NEAT1-1 can facilitate the development of bone metastasis in prostate cancer cells. Furthermore, the m6A level of NEAT1-1 can be used as a predictor of death in patients with prostate cancer bone metastases.³¹ ETTL3 upregulates the expression of LINC00958 by promoting its transcriptional stability, while LINC00958 acts as miR-378a-3p endogenous competing RNA that promotes YY1 and thus breast cancer progression.³² Collectively, these findings suggest that m6A, as a regulator of lncRNA, plays an essential role in the initiation and development of cancers. Hence, we should focus on the interactions between m6A modifications and lncRNAs to provide new insights into the diagnosis and treatment of tumours.

We analysed the prognostic association of six m6A-related lncRNAs with 552 cases of EC, of which three were identified to be highly relevant by multivariate Cox regression analysis. Microarray analysis in cervical cancer revealed that SCARNA9 expression is down-regulated in cancer tissues.³³ An analysis demonstrated that SCARNA9, an autophagy gene-related lncRNA, is associated with the prognosis of EC.³⁴ In glioma cells, TRAF3IP2-AS1 was identified using the lncRNA-mediated ceRNA network as one of the pivotal genes that may influence glioma initiation and development through signalling pathways such as cytokine-cytokine interactions and mTOR.³⁵ In bladder cancer, autophagy gene-associated lncRNA Z83843.1 is tightly connected with the overall survival of patients and can be considered as an independent prognostic predictor.³⁶ Moreover, a study showed that in lung adenocarcinoma, Z83843.1 was intimately associated with the metastasis of cancer cells.³⁷ Among the six lncRNAs analysed in this study, some have been reported to be associated with tumours, however, few have been reported in EC, and almost none were associated with m6A. Therefore, our bioinformatics analysis helps us to identify which m6A-related lncRNAs can be used as prognostic, immune microenvironmental targets in EC, and provide new insights for subsequent studies of tumour progression in EC.

However, in this study, we used the information of EC cases from the TCGA database and clinical data to analyse and validate them using different analytical methods. Nevertheless, some limitations exist in this study, as this work is based solely on bioinformatics analysis, and the interaction between these lncRNAs and m6A was not validated in vivo and in vitro.

Conclusion

The above analyses revealed that these lncRNAs were involved in endometriogenesis development through the regulation of m6A-related enzymes. Therefore, these lncRNAs may play a key role in the diagnosis and treatment of endometrial cancer as targets of prognosis and immune microenvironment.