Analysis of the prognostic significance and potential mechanisms of lncRNAs associated with m6A methylation in papillary thyroid carcinoma

doi:10.21203/rs.3.rs-807935/v2

Download PDF

Research Article

Analysis of the prognostic significance and potential mechanisms of lncRNAs associated with m6A methylation in papillary thyroid carcinoma

https://doi.org/10.21203/rs.3.rs-807935/v2

This work is licensed under a CC BY 4.0 License

Journal Publication

published 01 Dec, 2021

Read the published version in International Immunopharmacology →

Version 2

posted

You are reading this latest preprint version

Background : m6A methylation-related long non-coding RNAs (lncRNAs) play a significant role in the progression of various tumors and can be used as prognostic markers. However, whether m6A-related lncRNAs also play the same function as prognostic markers in papillary thyroid carcinoma (PTC) remains unclear.

Methods : Consensus cluster analysis was performed to divide PTC samples obtained from The Cancer Genome Atlas database into two clusters according to the expression of m6A-related lncRNAs. Then, the least absolute shrinkage and selection operator (LASSO) regression analysis was performed to create and verify a prognostic model. Furthermore, the relationship among risk scores, clusters, programmed death-ligand 1 (PD-L1), tumor microenvironment (TME), clinicopathological characteristics, immune infiltration, immune checkpoint, and tumor mutation burden (TMB) was analyzed. In addition, a nomogram was created, and subsequently, the drug sensitivity of lncRNAs in the prognostic model was analyzed. Finally, the relationship between these lncRNAs and prognosis in pan-cancer was investigated.

Results: The prognosis, RAS, BRAF, M, and TME were found to be different in two clusters. The prognostic model included three lncRNAs: PSMG3-AS1 , BHLHE40-AS1 , and AC016747.3 . The risk score was associated with clusters, PD-L1, tumor microenvironment, clinicopathological characteristics, immune cell infiltration, immune checkpoint, and TMB, and thus, risk score was confirmed as useful prognostic indicators. Differentially expressed lncRNAs are involved in many malignancies and can be identified as cancer prognostic makers.

Conclusion : According to our research, we can regard m6A-related lncRNAs involved in the procession of PTC as a biomarker of PFS for PTC patients, and pan-cancer.

Oncology

Cancer Biology

m6A

Papillary thyroid carcinoma

lncRNA

prognostic biomarker

PD-L1

The incidence of thyroid cancer, the most common endocrine cancer, has been continuously increasing in recent decades [1, 2]. Papillary thyroid carcinoma (PTC) is the most common type of thyroid cancer with an excellent prognosis [3–5]. However, approximately 30% of patients with PTC have a poor prognosis due to local recurrence or distant metastasis [6, 7]. Individuals with poorer prognoses should be crucially identified, and the mechanism of PTC development should be deeply understood.

Post-transcription modifications have become a significant regulator of many physiological processes and disease development, which attracts researchers’ interest in cancer studies [8]. A previous study has shown that N (6)-methyladenosine (m6A) is involved in the most common internal (non-cap) alteration of the messenger RNA (mRNA) in all higher eukaryotes, which are essential for cell viability and development [9]. There are three types of m6A regulators: the methyltransferase complex (MTC) often known as “writers” that catalyzes m6A; demethylase, commonly known as “eraser” that removes the m6A; and the RNA reader protein identifies m6A, attaches to the RNA, and performs the necessary activities. In cancer etiology and progression, writers, erasers, and readers are crosslinked [10, 11].

Researchers show great interest in long non-coding RNAs(lncRNAs) after understanding their functions. LncRNAs, the non-coding transcripts, are involved in at least 200 nucleotides [12]. However, their functions are still mostly unclear. LncRNAs may function as mRNA expression control via translational regulation, histone changes, or posttranscriptional regulation [13]. They can also impact numerous tumor cell characteristics (survival, proliferation, and migration) by gene control [12]. Aberrant lncRNA expression is involved in tumor aggressiveness and procession in various malignancies, including PTC. At this point, researchers have concluded that m6A, a necessary part of the complex lncRNA process, functions on tumor development [14]. However, the mechanism and prognostic value of m6A-related lncRNAs in PTC remains unknown.

Transcriptome data of PTC were obtained from a public source (TCGA-THCA database) and divided into two groups according to the consensus cluster analysis based on the expression level of m6A-related lncRNAs. To further accurately predict the progression-free survival (PFS) of patients with PTC, the prognostic model of m6A-related lncRNAs was developed and verified. The relationship among, programmed death-ligand 1 (PD-L1), tumor microenvironment (TME), risk score, clinicopathological characteristics, immune cell infiltration, tumor mutation burden (TMB), and clusters were also analyzed. Subsequently, to explore the predicted probability of the PFS in patients with PTC, a nomogram was drawn, and the drug sensitivity of three lncRNAs was analyzed in the prognostic model. Finally, the relationship among three lncRNAs and the prognosis in pan-cancer was also analyzed.

Data collections

Relevant genetic alterations and pathways that influence the tumor genesis, progression, differentiation, and metastasis [15] can be investigated because of a tremendous amount of genomic and clinical data available in the Cancer Genome Atlas (TCGA). The TCGA-THCA dataset contains 58 adjacent nontumorous tissues (N) and 510 patients with PTC (T) and 487 with simple nucleotide variation data (VarScan), which are used to calculate the TMB of patients with PTC. The University of California Santa Cruz (UCSC) genome database (https://xenabrowser.net; accessed 12 July 2021) provided the clinical data of 501 patients with PTC [16, 17]. Patients with incomplete follow-up information were excluded.

To broaden the scope of our study, 23 traditional regulators were included based on previous studies: writers (METTL3, METTL14, METTL16, WTAP, VIRMA [KIA1499], ZC3H13, RBM15, and RBM15B), readers (YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, HNRNPC, FMR1, LRPPRC, HNRNPA2B1, IGFBP1, IGFBP2, IGFBP3, and RBMX), and erasers (FTO and ALKBH5) [18-20].

Screening out lncRNAs associated with m6A methylation

After extracting expressions of 23 m6A regulators from the TCGA-THCA dataset, co-expression analysis between m6A regulators and m6A-related lncRNAs is performed using the “limma” package. Based on the “|correlation coefficient| of ≥0.5” and “p-value of < 0.001,” co-expression analysis was conducted and regulators of m6A methylation-related lncRNAs were selected using the Pearson correlation test. The expression information of the co-expression network for m6A RNA methylation-related lncRNA was also achieved in the “igraph” R package. m6A-related lncRNAs associated with the PFS of patients with PTC were screened using univariate COX progression analysis, and lncRNAs with p < 0.01 were included in the follow-up study.

Consensus Clustering Analysis

In the TCGA-THCA cohort, patients with PTC were divided according to clusters and examined the biological features of m6A-related lncRNAs using the R package “ConsensusClusterPlus.” Then, the PFS of clusters were estimated using the Kaplan-Meier analysis. The probable functionalities and downstream access of two clusters were observed using the Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression data (ESTIMATE) to evaluate the composition of the TME including the estimate score (stromal-immune comprehensive score), immune score (immune content), stromal score (stromal content). The R package “corrplot” was used to explore the correlation between the expression levels of PD-L1 and m6A-related lncRNAs.

Construction and validation of the prognostic model

After evenly dividing samples collected from TCGA-THCA into training and testing datasets, the least absolute shrinkage and selection operator (LASSO) regression analysis method was performed using the “glmnet” package [21] in the R4.1.0 to develop the prognostic model of m6A-related lncRNAs for PTC samples in the training dataset. The risk score of the LASSO regression model can be calculated as follows [22]:

where β_i represented the lncRNA coefficient and E_i represented the lncRNA expression. The risk score was calculated using this method in the training, testing, and TCGA-THCA datasets. Three datasets were divided into the high- and low-risk groups according to the median risk score. The area under the curve (AUC) was used to analyze the predictive power of signatures using the receiver operating characteristic (ROC) curves. The “stats” R package with the “prcomp” function was used to show principal component analysis (PCA) for the lncRNA expression. t-Distributed Stochastic Neighbor Embedding (t-SNE) was also performed using the “Rtsne” R package to explore variation from different groups in terms of their distribution. Heatmaps were drawn to display the distribution of clinicopathological characteristics between the high- and low-risk groups.

Analysis of the value of prognostic models

To assess the usefulness of risk score as independent prognostic variables, Cox regression models were used with other clinical characteristics (age [≥55 vs. <55 years], gender [male vs. female], American Joint Committee on Cancer [AJCC] stage [III+IV vs. I+II], AJCC M [M1 vs. M0], AJCC N [N1 vs. N0], AJCC T [T3+4 vs. T1+2], BRAF [mutated vs. wild-type), RAS [mutated vs. wild-type], and risk score (high vs. low]). A stratification study was conducted to determine whether risk score retained their capacity in analyzing PFS among different subgroups according to the order stratified analysis based on clinicopathological characteristics, which further evaluates the prognostic capability. A heatmap was created using the R package “pheatmap” that depicted the correlation between clinicopathological variables and groups. Given the significance of checkpoint inhibitor-based immunotherapies, differences in 47 immune checkpoint expressions were explored further between the two subgroups, and the connection among these lncRNAs was shown (47 immune checkpoints are provided in Supplementary List 1). Subsequently, we investigated the correlation among risk scores, immune cells, or TMB using the Pearson correlation analysis.

Sensitivity Analysis of Chemotherapy

In the cancer research community, a database, known as Cellminer, was developed to determine the integration and analysis of molecular and pharmacological data for the NCI-60 malignant cell lines collected by the United States Developmental Therapeutics Program. More than 100,000 chemical compounds and natural items will be screened by the National Cancer Institute (since 1990) (https://discover.nci.nih.gov/cellminer/home.do). The Pearson correlation analysis was used to evaluate the connection between prognostic lncRNAs expression and drug sensitivity. The correlation study was used to explore the effectiveness of 263 medicines that have been authorized by the Food and Drug Administration or are currently undergoing clinical studies.

Prognosis-related lncRNA’s pan-cancer value promotion

The Gene Expression Profiling Interactive Analysis 2 (GEPIA2) (http://gepia2.cancer-pku.cn/#index), which can integrate the TCGA and Genotype-Tissue Expression Project (GTEx) database, was used to analyze the pan-cancer value about prognostic-related lncRNAs [23]. Expression levels of three prognostic-related lncRNAs were first investigated in 33 tumors and normal tissues (p-value <0.05, |logFC| >0.1, differential methods: limma). We also explored the relationship between lncRNA and the prognosis of 33 tumors.

Statistical analysis

The Kaplan-Meier survival analysis (PFS considered as the end-point) was performed using the “survival” package in R. The chi-square test was performed to explore differences in clinical parameters of the high- and low-risk groups in the training, testing, and overall TCGA-THCA datasets. The Pearson’s correlation analysis was performed to assess related variables. A pair of independent t-test samples was used to compare the results in the Mann-Whitney test to explore differences between high- and low - risk samples. A log-rank test was used to calculate the p-values for PFS. Statistical analyses were performed and visualized in R version 4.1.0, IBM SPSS Statistics for Windows version 25.0 (IBM Corp., Armonk, NY, USA), and GraphPad Prism version 8.0 (GraphPad Software Inc., La Jolla, CA, USA).

Analyzing the consensus cluster based on median expression levels of m6A-related lncRNAs in patients with PTC

We firstly found that 630 lncRNAs co-expressed with 15 m6A methylation regulators and then used them to visualize a co-expression network (Supplementary Figure 1). Red nodes show the m6A RNA methylation regulators, and blue nodes display m6A-related lncRNAs.

To further investigate the relationship between lncRNAs and the PFS of patients with PTC, we found prognostic-related lncRNAs (Figure 2A, B) using the univariate COX analysis and collected 21 lncRNAs (AC010168.2, CAPN10-DT, PSMG3-AS1, AL662844.4, AC010761.3, AC093752.3, LMNTD2-AS1, AL031705.1, CD44-AS1, U73166.1, BHLHE40-AS1, MZF1-AS1, AL360181.2, AC138956.1, AL359878.1, AC002128.1, AC084125.2, AC139795.2, AC016747.3, AL355987.4, and AC138956.2) for subsequent studies.

Based on the expression of 21 m6A-related lncRNAs, the TCGA-THCA cohort was divided into distinct clusters using the “ConsensusClusterPlus” R package (Figure 2C, D, E). Regarding 2 as the consensus matrix’s k-value, the least crossover was obtained from all PTC samples. The difference in PFS was explored between two clusters using the Kaplan-Meier technique. PFS was significantly longer in cluster 1 (PFS, p = 0.005) as compared to cluster 2 (Figure 2F). After comparing the clinicopathological characteristics of two clusters and analyzing their connection, RAS, BRAF, and AJCC M stage were all significantly related to our cluster analysis (Figure 2G). The percentage of the 22 immune cell types in two clusters was analyzed using the CIBERSORT. CD8 T cells, CD4 memory actived T cells, follicular helper T cells, macrophages M1, and activated dendritic cells were significantly different between the two clusters (Figure 2H). The TME was examined using the ESTIMATE, stromal, and immune scores and found significant differences in the two groups (Figure 2I, J, K). Considering the important function of PD-L1 in tumor treatment, a correlation was observed among the PD-L1 and 21 m6A-related lncRNAs (Figure 2L) and the expression of PD-L1 between PTC samples and normal tissues in clusters 1 and 2 (Supplementary Figure 2A, B). The expression of PD-L1 was negatively correlated with 18 lncRNAs. Statistically significant changes in PD-L1 expression were observed between tumor and healthy control. However, compared to cluster 2, expression levels of PD-L1 expression were higher in cluster 1.

Identification and verification of the prognostic model of m6A-related lncRNAs

Next, the effectiveness of the prognostic model was evaluated. After randomly dividing the entire TCGA-THCA dataset into two groups (training dataset [n = 252] and testing dataset [n = 249]), a prognostic model was obtained in the training dataset using the LASSO regression analysis based on the expression of 21 lncRNAs (Supplementary Figure 3A, B). A total of three lncRNAs (PMSG3-AS1, BHLHE40-AS1, and ACO16747.3) were included in the prognostic model. The risk score of patients in the training, testing, and entire TCGA-THCA dataset were calculated by the prognostic model. The following prognostic model formula was used: risk score = (0.264992826432626*PSMG3-AS1 expression level) + (-0.777201837590316*BHLHE40-AS1 expression level) + (1.34421059388075*AC016747.3 expression level). According to the median risk score, patients in all three datasets were divided into high- and low-risk groups. The associations among the risk score, prognostic status, and expressive signatures of three m6A-related lncRNAs in the training, testing, and entire TCGA-THCA cohorts are shown (Supplementary Figure 4A-I). The PFS between the two groups was further analyzed, with significant differences, whether in the training dataset (p = 0.004, Figure 3A), the testing dataset (p = 0.001, Figure 3B), or the entire TCGA-THCA dataset (p < 0.001, Figure 3C). As shown in the time-dependent ROC curves, significant differences were found between the high- and low-risk groups. AUC analyses for 1-, 2-, 3-year PFS prediction according to the lncRNA signature were 0.747, 0.626, and 0.722 in the training dataset (Figure 3D); 0.709, 0.651, and 0.771 in the testing dataset (Figure 3E); 0.733, 0.634, and 0.739 in the entire TCGA-THCA dataset (Figure 3F). PCA and t-SNE analyses displayed that patient in the three datasets were dispersed in two methods (Figure 3G-L). Thus, the risk score is relatively accurate in determining the PFS of patients with PTC because it is greatly distinct among patients.

Univariate and multivariate Cox analyses were performed to examine the association between the risk score and PFS. A risk score was found as a risk factor for PFS in all three datasets, as demonstrated by the univariate cox analysis results (training dataset: hazard ratio [HR] = 2.767, 95% confidence interval [CI] = 1.891-4.050, p<0.001, Figure 3M; testing dataset: HR = 6.817, 95%CI = 2.375-19.561, p<0.001, Figure 3N; entire TCGA-TCHA dataset: HR = 3.147, 95%CI = 2.260-4.383, p<0.001, Figure 3O). In the multivariate Cox regression analysis, the risk score remained a significant predictor of PFS even considering the effects of the other factors (training dataset: HR = 2.849, 95%CI = 1.702-4.771, p<0.001, Figure 3P; testing dataset: HR = 40.147, 95%CI = 3.020-533.628, p=0.005, Figure 3Q; entire TCGA-THCA dataset: HR = 2.458, 95%CI = 1.593-3.794, p<0.001, Figure 3R).

Stratification analysis based on clinicopathological characteristics

The prognostic profile associated with m6A-mediated lncRNAs might be used to predict PFS in both high- and low-risk groups. A stratified analysis on clinicopathological variables showed appropriate results, including age (≤55 years vs. >55 years), gender (female vs. male), AJCC stage (stage Ⅰ + Ⅱ vs. Ⅲ + Ⅳ), AJCC T stage (T Ⅰ + Ⅱ vs. T Ⅲ + Ⅳ), AJCC M stage (M0 vs. M1), AJCC N stage (N0 vs. N1), BRAF (wild-type vs. mutated), and RAS (wild-type vs. mutated) (p < 0.05, Figure 4). The Kaplan-Meier method of survival analysis has been performed, and results have shown that the high-risk group has a shorter PFS than the low-risk group in the number of clinical strata, except for patients aged <55 years, and patients with RAS mutation (Figure 4A). This additional proof confirmed our hypothesis on its significance.

The correlation between risk score and clinicopathology characteristics, cluster analysis results, or immune scores

When comparing high- with low-risk samples from the TCGA-THCA dataset with clinicopathological factors, immune scores, and cluster analysis results, expression levels of PSMG3-AS1 and AC016747.3 were found to be higher in the high-risk group (Figure 5A), indicating high-risk lncRNAs. Conversely, the expression level of BHLHE40-AS1 in low-risk samples was higher than in high-risk samples, indicating a low-risk lncRNA. A higher-risk value can be observed in the RAS-mutated (Figure 5B), BRAF wild-type (Figure 5C), N0 (Figure 5D), age ≥55 years (Figure 5E), ImmuneScore low (Figure 5F), and cluster 2 (Figure 5G). Then, we found differences in the expression levels of 30 immune checkpoints between high- and low-risk samples (Figure 5H).

Correlation between risk score with immunocytes and tumor mutation burden and prognostic nomogram

After performing the CIBERSORT analysis, Pearson correlation was used to analyze the correlation between risk scores and 22 immune cells. A positive correlation can be observed between regulating T cells (Tregs) and risk score (Figure 6C), whereas the correlation between resting dendritic cell (Figure 6A) or resting mast cell (Figure 6B) and risk scores was negative. A positive correlation also be observed between TMB and risk score (Figure 6D). To quantitatively predict the patient’s prognosis, a nomogram was created based on clinicopathological characteristics and risk scores to accurately predict the patient’s PFS (Figure 6E).

Sensitivity analysis of chemotherapy

Two prognosis-related lncRNAs (PSMG3-AS1, BHLHE40-AS1) were found to correspond to 11 drugs (Figure 7). As the PSMG3-AS1 expression increases, the sensitivity of cancer cells to these drugs (mitomycin C, mitoxantrone, tegafur, idarubicin, gemcitabine, fulvestrant, and teniposide) increases. As the PSMG3-AS1 expression increases, the resistance of cancer cells to the drug (bosutinib) increases. The sensitivity of cancer cells to epirubicin, mitoxantrone, and doxorubicin and their resistance to allopurinol increase while increasing the BHLHE40-AS1 expression levels. Remarkably, as the expression of these two lncRNAs increases, the sensitivity of cancer cells to mitoxantrone increases (Figure 7C, H).

Prognosis-related lncRNA’s pan-cancer value promotion

We found that these three lncRNAs are differentially expressed in most tumors (Figure 8A, C, E). The results show that the overall survival rate is significantly different with the increasing PSMG3-AS1 and BHLHE40-AS1 expression levels (Figure 8B, D). A significant difference in disease-free survival rate was observed between patients with low- and high expressions of AC016747.3 (Figure 8F).

The role of post-transcriptional modification in life activities has increasingly received attention [8]. Among the many RNA modifications, m6A is the most common type of modification and is involved in the occurrence and development of different tumors [14, 24]. Some studies have reported that METTL3 inhibits the progression of thyroid papillary carcinoma by affecting neutrophil infiltration through m6A /c-Rel/IL-8 [25] and regulating m6A methylation on TCF1 [26]. Another study has found that the protein content of FTO correlated with lymph node metastasis and tumor grade is an independent prognostic factor of thyroid cancer grade [27]. Some studies have shown that m6A-methylation-related regulators can predict the prognosis of patients with PTC using bioinformatics methods [28–30]. Therefore, the formation and prognosis of patients with PTC can be evaluated with m6A methylation. The function of lncRNAs in different cancers has greatly increased the interest of researchers. Thus, studies have claimed that m6A is a key step in the lncRNA process that plays a role in the pathological process of tumors [31], and we can infer the prognosis of patients through the expression of m6A-methylation-related lncRNAs, such as pancreatic adenocarcinoma [32], ovarian cancer [33], hepatocellular carcinoma [34], triple-negative breast cancer [35], primary head and neck squamous cell carcinoma [14], kidney renal clear cell carcinoma [36], metastatic skin cutaneous melanoma [37], and lower-grade glioma [38]. No studies have shown the possibility of predicting the prognosis of patients with PTC by m6A-related lncRNAs. Therefore, we conducted this study using Pearson’s method to obtain 630 lncRNAs co-expressed with m6A genes based on the expression levels of 15 m6A genes.

We performed the univariate COX analysis to screen out 21 lncRNAs associated with prognosis and draw forest maps and heat maps. The results showed that 20 lncRNAs are risk factors for prognosis, and one lncRNA is a protective factor for prognosis. According to the expression of 21 lncRNAs, TCGA-THCA cohert were divided into two cluster by consensus clustering method. The results show that patients in cluster2 have a worse prognosis, with higher expression of 21 lncRNAs.

To further study the function of these lncRNAs, we analyzed fractions of infiltrating immune cells between two clusters by CIBERSORT, and the results showed that the content of CD8 T cells, CD4 memory active T cells, follicular helper T cells, macrophages M1, and dendritic cells activated were statistically different, indicating that the immune infiltration between the two clusters was different. Therefore, the ESTIMATE analysis method was used to further explore the tumor microenvironment between the two clusters. Cluster 2 was also found to show a lower estimate score, matrix score, and immune score, which emphasized the role of these lncRNAs in the TME. Considering the importance of PD-L1 during the tumor immunity process [39–42], a correlation between PD-L1 and 21 lncRNAs was observed. The results showed that PD-L1 was negatively correlated with several lncRNAs, and most of the lncRNAs are positively related.

Two datasets (a training group and a testing group) were created by a random method and obtained a prognostic model containing three m6A-related lncRNAs (PSMG3-AS1, BHLHE40-AS1, and AC016747.3) by the LASSO regression analysis in the training group, all of which can be used to further increase the feasibility of our prediction.

Previous reports have indicated that PSMG3-AS1 sponging miR-613 to upregulate SphK1 increases the growth of non-small cell lung cancer cells [43]; and sponging miR‑143‑3p to increase the proliferation and migration of breast cancer cells [44]; and can function on the migration, proliferation, and invasion of cervical squamous cell carcinoma [45], lung adenocarcinoma [46], hepatocellular carcinoma [47], renal cell carcinoma [48], and breast cancer cells; can distinguish glioblastoma and sarcoidosis [49]; and can potentially predict breast cancer [50] and ovarian cancer [51]. Some studies have shown that early breast cancer progresses by increasingly expressed BHLHE40-AS1 functioning on IL-6/STAT3 signaling [52]. AC016747.3 has not been reported, and these three lncRNAs have not been investigated in PTC; therefore, follow-up in vivo and in vitro experiments on these lncRNAs in PTC is necessary.

Furthermore, the prognostic model was verified in the entire TCGA-THCA dataset and the testing group. Based on the prognostic model, the risk score of patients was calculated. According to the median risk value, we divided patients with PTC into high- and low-risk groups. Our analytical results show that patients from high-risk groups have worse PFS than low-risk groups (training dataset: p = 0.004; testing dataset: p = 0.001; TCGA-THCA: p < 0.001); ROC curve, PCA, and t-SNE analysis displayed that the prognostic model can be used to predict the patient’s PFS to a certain extent. We got an analytical result according to the univariate and multivariate COX analyses: risk score is a prognostic factor independent of other clinical traits for patients with PTC. In the stratified analysis based on clinicopathological characteristics, the prognosis of patients in different groups is also statistically different. The analytical results also showed that some significant difference was observed in the risk scores among several clinicopathological characteristics (RAS, BRAF, N, age, immune score, and cluster) and is correlated with three types of immune cells (resting dendritic cell, resting mast cell, and regulatory T cells). These results increase our understanding of the function of m6A-related lncRNAs in PTC and provide directions for further experiments. To further expand the clinical value of our study, we developed a nomogram to predict the PFS of patients with PTC, performed a drug sensitivity analysis among these three lncRNAs, and found that with the expression of two lncRNAs increased sensitivity or resistance to certain chemotherapeutic drugs. Notably, as the expression level of these two lncRNAs increases, the sensitivity of tumor cells to mitoxantrone increases. Moreover, these three lncRNAs also have a certain prognostic value in pan-cancer.

This study has some limitations that should be considered. First, it is the data source from the TCGA-THCA cohort used to analyze and verify the prognostic significance of m6A-related lncRNAs. To further verify the prognostic model, we still lack more central data. Second, only the correlation between the expression of these three lncRNAs and clinicopathological characteristics, immune checkpoints, drug sensitivity, or pan-cancer was considered; however, we did not analyze pathways and mechanisms of their effects in in vivo and in vitro experiments.

In summary, m6A-related lncRNAs of patients with PTC were obtained based on public databases and constructed and verified a prognostic model. The three m6A-related lncRNAs in the prognostic model can be used as prognostic markers for patients with PTC. In addition, the functional analysis also provides a new perspective in future studies on the role and mechanism of m6A-related lncRNAs in patients with PTC.

LncRNAs: long non-coding RNAs; PTC: papillary thyroid carcinoma; LASSO: least absolute shrinkage and selection operator; TCGA: The Cancer Genome Atlas database; PD-L1:programmed death-ligand 1, TME: tumor microenvironment; TMB: tumor mutation burden; m6A:N (6)-methyladenosine; mRNA: messenger RNA; UCSC: University of California Santa Cruz; PFS: progression-free survival; ESTIMATE: Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression data; AUC: area under the curve; ROC: receiver operating characteristic; PCA: principal component analysis; t-SNE: t-Distributed Stochastic Neighbor Embedding; AJCC: American Joint Committee on Cancer; GEPIA2: Gene Expression Profiling Interactive Analysis 2; GTEx:Genotype-Tissue Expression Project database; HR: hazard ratio; CI: confidence interval.

Acknowledgements

We appreciate the TCGA database for allowing access to their data.

Authors' contributions

YDH created the analytical strategies and drafted the paper. XL verified the data analysis method. XYL, YL, XL, YZH, SW, SYW, YCH, and HJ performed data analyses. JZ designed and wrote the paper. Each author equally contributed to the article, and all of them approved the final draft that was eventually published.

Funding

The National Natural Science Foundation of China sponsored this study (81600602).

Availability of data and materials

The data of current study are available from TCGA dataset.

Declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

The authors agree for publication.

Competing interests

The writers indicated that their work was not influenced by competing financial interests or personal connections.

Pellegriti G, Frasca F, Regalbuto C, Squatrito S, Vigneri R (2013) Worldwide increasing incidence of thyroid cancer: update on epidemiology and risk factors. J Cancer Epidemiol 2013:965212
Rajabi S, Dehghan MH, Dastmalchi R, Jalali Mashayekhi F, Salami S, Hedayati M (2019) The roles and role-players in thyroid cancer angiogenesis. Endocr J 66:277–293
Schneider DF, Chen H (2013) New developments in the diagnosis and treatment of thyroid cancer. CA Cancer J Clin 63:374–394
Fagin JA, Wells SA (2016) Jr.: Biologic and Clinical Perspectives on Thyroid Cancer. N Engl J Med 375:1054–1067
Xing M (2013) Molecular pathogenesis and mechanisms of thyroid cancer. Nat Rev Cancer 13:184–199
Lamartina L, Grani G, Durante C, Filetti S (2018) Recent advances in managing differentiated thyroid cancer. F1000Res 7:86
Davies L, Morris LG, Haymart M, Chen AY, Goldenberg D, Morris J, Ogilvie JB, Terris DJ, Netterville J, Wong RJ, Randolph G (2015) AMERICAN ASSOCIATION OF CLINICAL ENDOCRINOLOGISTS AND AMERICAN COLLEGE OF ENDOCRINOLOGY DISEASE STATE CLINICAL REVIEW: THE INCREASING INCIDENCE OF THYROID CANCER. Endocr Pract 21:686–696
He L, Li H, Wu A, Peng Y, Shu G, Yin G (2019) Functions of N6-methyladenosine and its role in cancer. Mol Cancer 18:176
Wang X, Lu Z, Gomez A, Hon GC, Yue Y, Han D, Fu Y, Parisien M, Dai Q, Jia G et al (2014) N6-methyladenosine-dependent regulation of messenger RNA stability. Nature 505:117–120
Deng X, Su R, Weng H, Huang H, Li Z, Chen J (2018) RNA N(6)-methyladenosine modification in cancers: current status and perspectives. Cell Res 28:507–517
Panneerdoss S, Eedunuri VK, Yadav P, Timilsina S, Rajamanickam S, Viswanadhapalli S, Abdelfattah N, Onyeagucha BC, Cui X, Lai Z et al (2018) Cross-talk among writers, readers, and erasers of m(6)A regulates cancer growth and progression. Sci Adv 4:eaar8263
Ramilowski JA, Yip CW, Agrawal S, Chang JC, Ciani Y, Kulakovskiy IV, Mendez M, Ooi JLC, Ouyang JF, Parkinson N et al (2020) Functional annotation of human long noncoding RNAs via molecular phenotyping. Genome Res 30:1060–1072
Huang Q, Yan J, Agami R (2018) Long non-coding RNAs in metastasis. Cancer Metastasis Rev 37:75–81
Feng ZY, Gao HY, Feng TD (2021) Immune Infiltrates of m(6)A RNA Methylation-Related lncRNAs and Identification of PD-L1 in Patients With Primary Head and Neck Squamous Cell Carcinoma. Front Cell Dev Biol 9:672248
Cancer Genome Atlas Research Network, Weinstein JN, Collisson EA, Mills GB, Shaw KR, Ozenberger BA, Ellrott K, Shmulevich I, Sander C, Stuart JM. The Cancer Genome Atlas Pan-Cancer analysis project. Nat Genet. 2013 Oct;45(10):1113–1120. doi:10.1038/ng.2764. PMID: 24071849; PMCID: PMC3919969
Goldman MJ, Craft B, Hastie M, Repečka K, McDade F, Kamath A, Banerjee A, Luo Y, Rogers D, Brooks AN et al (2020) Visualizing and interpreting cancer genomics data via the Xena platform. Nat Biotechnol 38:675–678
Haeussler M, Zweig AS, Tyner C, Speir ML, Rosenbloom KR, Raney BJ, Lee CM, Lee BT, Hinrichs AS, Gonzalez JN et al (2019) The UCSC Genome Browser database: 2019 update. Nucleic Acids Res 47:D853–Dd858
Tang J, Wan Q, Lu J (2020) The prognostic values of m6A RNA methylation regulators in uveal melanoma. BMC Cancer 20:674
Tian S, Lai J, Yu T, Li Q, Chen Q (2020) Regulation of Gene Expression Associated With the N6-Methyladenosine (m6A) Enzyme System and Its Significance in Cancer. Front Oncol 10:623634
He PC, He C (2021) m(6) A RNA methylation: from mechanisms to therapeutic potential. Embo j 40:e105977
Friedman J, Hastie T, Tibshirani R (2010) Regularization Paths for Generalized Linear Models via Coordinate Descent. J Stat Softw 33:1–22
Huang R, Li G, Wang Z, Hu H, Zeng F, Zhang K, Wang K, Wu F (2020) Identification of an ATP metabolism-related signature associated with prognosis and immune microenvironment in gliomas. Cancer Sci 111:2325–2335
Tang Z, Kang B, Li C, Chen T, Zhang Z (2019) GEPIA2: an enhanced web server for large-scale expression profiling and interactive analysis. Nucleic Acids Res 47:W556–Ww560
Liu X, Wang P, Teng X, Zhang Z, Song S (2021) Comprehensive Analysis of Expression Regulation for RNA m6A Regulators With Clinical Significance in Human Cancers. Front Oncol 11:624395
He J, Zhou M, Yin J, Wan J, Chu J, Jia J, Sheng J, Wang C, Yin H, He F (2021) METTL3 restrains papillary thyroid cancer progression via m(6)A/c-Rel/IL-8-mediated neutrophil infiltration. Mol Ther 29:1821–1837
Wang K, Jiang L, Zhang Y, Chen C (2020) Progression of Thyroid Carcinoma Is Promoted by the m6A Methyltransferase METTL3 Through Regulating m(6)A Methylation on TCF1. Onco Targets Ther 13:1605–1612
Tian R, Zhang S, Sun D, Bei C, Li D, Zheng C, Song X, Chen M, Tan S, Zhu X, Zhang H: M6A Demethylase FTO Plays a Tumor Suppressor Role in Thyroid Cancer. DNA Cell Biol 2020
Hou J, Shan H, Zhang Y, Fan Y, Wu B (2020) m(6)A RNA methylation regulators have prognostic value in papillary thyroid carcinoma. Am J Otolaryngol 41:102547
Wang X, Fu X, Zhang J, Xiong C, Zhang S, Lv Y (2020) Identification and validation of m(6)A RNA methylation regulators with clinical prognostic value in Papillary thyroid cancer. Cancer Cell Int 20:203
Xu N, Chen J, He G, Gao L, Zhang D (2020) Prognostic values of m6A RNA methylation regulators in differentiated Thyroid Carcinoma. J Cancer 11:5187–5197
Lan Y, Liu B, Guo H (2021) The role of M(6)A modification in the regulation of tumor-related lncRNAs. Mol Ther Nucleic Acids 24:768–779
Yuan Q, Ren J, Li L, Li S, Xiang K, Shang D (2021) Development and validation of a novel N6-methyladenosine (m6A)-related multi- long non-coding RNA (lncRNA) prognostic signature in pancreatic adenocarcinoma. Bioengineered 12:2432–2448
Zheng J, Guo J, Cao B, Zhou Y, Tong J (2021) Identification and validation of lncRNAs involved in m6A regulation for patients with ovarian cancer. Cancer Cell Int 21:363
Li L, Xie R, Lu G: Identification of m6A methyltransferase-related lncRNA signature for predicting immunotherapy and prognosis in patients with hepatocellular carcinoma. Biosci Rep 2021, 41
Wu J, Cai Y, Zhao G, Li M (2021) A ten N6-methyladenosine-related long non-coding RNAs signature predicts prognosis of triple-negative breast cancer. J Clin Lab Anal 35:e23779
Yu J, Mao W, Sun S, Hu Q, Wang C, Xu Z, Liu R, Chen S, Xu B, Chen M (2021) Identification of an m6A-Related lncRNA Signature for Predicting the Prognosis in Patients With Kidney Renal Clear Cell Carcinoma. Front Oncol 11:663263
Huang S, Lyu S, Gao Z, Zha W, Wang P, Shan Y, He J, Li Y (2021) m6A-Related lncRNAs Are Potential Biomarkers for the Prognosis of Metastatic Skin Cutaneous Melanoma. Front Mol Biosci 8:687760
Tu Z, Wu L, Wang P, Hu Q, Tao C, Li K, Huang K, Zhu X (2020) N6-Methylandenosine-Related lncRNAs Are Potential Biomarkers for Predicting the Overall Survival of Lower-Grade Glioma Patients. Front Cell Dev Biol 8:642
Ishida Y, Agata Y, Shibahara K, Honjo T (1992) Induced expression of PD-1, a novel member of the immunoglobulin gene superfamily, upon programmed cell death. Embo j 11:3887–3895
Francisco LM, Sage PT, Sharpe AH (2010) The PD-1 pathway in tolerance and autoimmunity. Immunol Rev 236:219–242
Asano T, Meguri Y, Yoshioka T, Kishi Y, Iwamoto M, Nakamura M, Sando Y, Yagita H, Koreth J, Kim HT et al (2017) PD-1 modulates regulatory T-cell homeostasis during low-dose interleukin-2 therapy. Blood 129:2186–2197
Sun Y, Tan J, Miao Y, Zhang Q (2021) The role of PD-L1 in the immune dysfunction that mediates hypoxia-induced multiple organ injury. Cell Commun Signal 19:76
Shen W, Hong X, Jin C, Xi Y (2021) LncRNA PSMG3AS1 promotes proliferation of non-small cell lung cancer cells by sponging miR-613 to upregulate SphK1. Cell Cycle 20:829–838
Cui Y, Fan Y, Zhao G, Zhang Q, Bao Y, Cui Y, Ye Z, Chen G, Piao X, Guo F et al (2020) Novel lncRNA PSMG3–AS1 functions as a miR–143–3p sponge to increase the proliferation and migration of breast cancer cells. Oncol Rep 43:229–239
Man S, Li X, Zhu W (2021) miR-4417 targets lncRNA PSMG3-AS1 to suppress cell invasion and migration in cervical squamous cell carcinoma. Oncol Lett 22:502
Yue N, Ye M, Zhang R, Guo Y (2020) MiR-449b-5p targets lncRNA PSMG3-AS1 to suppress cancer cell proliferation in lung adenocarcinoma. BMC Pulm Med 20:152
Zhang J, Huang J, Chen W, Hu Z, Wang X (2020) miR-143-3p Targets lncRNA PSMG3-AS1 to Inhibit the Proliferation of Hepatocellular Carcinoma Cells. Cancer Manag Res 12:6303–6309
Cheng G, Liu D, Liang H, Yang H, Chen K, Zhang X (2019) A cluster of long non-coding RNAs exhibit diagnostic and prognostic values in renal cell carcinoma. Aging 11:9597–9615
Chen L, Wang G, Xu Z, Lin K, Mu S, Pan Y, Shan M (2020) Overexpression of LncRNA PSMG3-AS1 Distinguishes Glioblastomas from Sarcoidosis. J Mol Neurosci 70:2015–2019
Gong PJ, Shao YC, Huang SR, Zeng YF, Yuan XN, Xu JJ, Yin WN, Wei L, Zhang JW (2020) Hypoxia-Associated Prognostic Markers and Competing Endogenous RNA Co-Expression Networks in Breast Cancer. Front Oncol 10:579868
Li N, Zhan X (2019) Identification of clinical trait-related lncRNA and mRNA biomarkers with weighted gene co-expression network analysis as useful tool for personalized medicine in ovarian cancer. Epma j 10:273–290
DeVaux RS, Ropri AS, Grimm SL, Hall PA, Herrera EO, Chittur SV, Smith WP, Coarfa C, Behbod F, Herschkowitz JI (2020) Long noncoding RNA BHLHE40-AS1 promotes early breast cancer progression through modulating IL-6/STAT3 signaling. J Cell Biochem 121:3465–3478

SupplementaryFigure1.tif
Supplementary Figure 1. Interaction network between m6A-methylated genes and related lncRNAs.
SupplementaryFigure2.tif
Supplementary Figure 2. PD-L1 expression between tumor and normal tissues in clusters 1/2.
SupplementaryFigure3.tif
Supplementary Figure3. The process of constructing a prognostic model through LASSO regression analysis.
SupplementaryFigure4.tif
Supplementary Figure4. The median value of risk score, prognostic status, and expression of m6A-related lncRNAs are in three datasets
SupplementaryList1.txt
Supplementary List1. Forty-seven immune checkpoints

Download PDF

Journal Publication

published 01 Dec, 2021

Read the published version in International Immunopharmacology →

Version 2

posted

You are reading this latest preprint version

Analysis of the prognostic significance and potential mechanisms of lncRNAs associated with m6A methylation in papillary thyroid carcinoma

Status:

Journal Publication

Version 2

Abstract

Figures

Introduction

Materials And Methods

Results

Discussion

Conclusion

Abbreviations

Declarations

References

Supplementary Files

Status:

Journal Publication

Version 2