A novel necroptosis-related genes signature for predicting the prognosis of thyroid cancer

doi:10.21203/rs.3.rs-1548457/v1

Download PDF

Research Article

A novel necroptosis-related genes signature for predicting the prognosis of thyroid cancer

https://doi.org/10.21203/rs.3.rs-1548457/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Background: Necroptosis, a type of programmed cell death, has been implicated in a variety of cancer-related biological processes. However, the roles of necroptosis-related genes in thyroid cancer yet remain unknown.

Methods: A necroptosis-related gene signature was constructed using the least absolute shrinkage and selection operator (LASSO) regression analysis and Cox regression analysis. The predictive value of the prognostic signature was validated in an internal cohort. Additionally, the single-sample gene set enrichment analysis (ssGSEA) was used to examine the relationships between necroptosis and immune cells, immunological functions, and immune checkpoints. Finally, the modeled gene expression was validated in clinical tumor and normal tissue samples.

Results: In this study, the risk signature of seven necroptosis-related genes was created to predict the prognosis of papillary thyroid cancer (PTC) patients, and all patients were divided into high- and low-risk groups. Patients in the high-risk group fared worse in terms of overall survival than those in the low-risk group. The area under the curve (AUC) of the receiver operator characteristic (ROC) curves proved the predictive capability of created signature. The risk score was found to be an independent risk factor for prognosis in multivariate Cox analysis. The low-risk group showed increased immune cell infiltration and immunological activity, implying that they might respond better to immune checkpoint inhibitor medication. Finally, GEO database and qRT-PCR were used to validate the expression of the seven modeled genes in PTCs, and the results were compatible with TCGA database.

Conclusions: The created necroptosis-related risk signature may have an effective predictive in the and prognosis of PTC. The findings of this study could pave the way for further research into the link between necroptosis and tumor immunotherapy.

Papillary thyroid cancer

Necroptosis

Prognosis

Gene signature

TCGA

Thyroid cancer is becoming more common in the world. One key explanation is the improved detection of small papillary thyroid cancers (PTC) using imaging technologies such as ultrasonography (US) has advanced and become more widely used (1). Aside from PTC, there are three more kinds of thyroid cancer: follicular thyroid cancer (FTC), medullary thyroid cancer (MTC), and anaplastic thyroid cancer (ATC). PTC is the most frequent kind, accounting for more than 80% of all pathological types (2–4). Although PTC has an overall excellent prognosis, aggressive PTC variants such as the tall cell (TC) and diffuse sclerosing (DS) variants are on the rise and are associated with more aggressive pathologic characteristics that can lead to a poor prognosis or even recurrence in certain patients (5–7). Therefore, searching for new prognostic markers is crucial for the diagnosis, prognosis, and treatment of PTC patients.

Necroptosis, discovered and described in 2005, is a programmed form of necrotic cell death that is distinct from apoptosis, ferroptosis, and pyroptosis (8). AS a caspase-independent and regulated necrotic cell death mechanism, necroptosis that is mechanistically similar to apoptosis and morphologically similar to necrosis (9). Necroptosis is regulated by molecules that are also known to govern apoptosis, but it is dependent on the creation of the necrosome, which consists of receptor-interacting serine/threonine-protein kinase 1 (RIPK1) and RIPK3, which activates the pseudokinase mixed lineage kinase (MLKL). MLKL then mediates the release of intracellular materials, which results in the execution of the necroptotic program (10, 11). Necroptosis is not only crucial in the development and progression of numerous immune system disorders and viral infections, but it also plays an important role in cancer biology regulation, including oncogenesis, cancer metastasis, cancer immunity, and cancer subtypes (12, 13). Some dual effects of necroptosis on cancer have been demonstrated as a combination of apoptosis and necrosis. Several articles have reported that necroptosis can trigger inflammatory responses and reportedly promotes cancer metastasis and immunosuppression as a necrotic cell death modality (14, 15). Necroptosis is viewed as a barrier that can limit tumor formation when apoptosis is disrupted. As a result, necroptosis has been described as both a friend and a foe of cancer. The expression of necroptosis factors in cancer, as well as their impact on cancer prognosis, is likewise complex. Multiple studies have discovered that numerous important components in necroptotic signaling pathways are downregulated in a variety of malignancies, implying that cancer cells may avoid necroptosis to survive. RIPK3 is downregulated in a variety of malignancies, including breast cancer, colorectal cancer, and melanoma(16–18). Furthermore, decreased RIPK3 expression has been observed to independently predict lower overall survival in colorectal cancer and breast cancer (16, 18). It has also been observed that RIPK1 expression is downregulated in head and neck squamous cell carcinoma, which is associated with disease progression (19). These findings imply that the necrosis pathway has an anti-cancer effect in cancer. However, it does not appear that necroptotic factors are downregulated in all malignancies. These variables have been identified to be elevated in various malignancies, and their overexpression is positively connected with tumor growth. In glioblastoma, for example, the RIPK1 expression is overexpressed in around 30% of cases (grade IV) which leads to a worse prognosis (20). A recent study has attempted to discover new necroptosis-related signatures in some cancers. Based on necroptosis-related genes, a predictive signature for pancreatic adenocarcinoma was created (21). Necroptosis-related genes also play an essential role in tumor immunity in prostate cancer and can be used to predict the prognosis of prostate adenocarcinoma (22). In stomach adenocarcinoma, a necroptosis-related prognostic signature and a lncRNA SNHG1/miR-21-5p/TLR4 regulatory axis were discovered. However, the prognostic significance of necroptosis-related genes in PTC remains unknown.

Therefore, this study aimed to construct a novel prognostic signature for PTC by detecting the expression and prognostic value of necroptosis-related genes (NRGs) and evaluation of their prognostic value from many aspects. Additionally, a thorough analysis of the link between the signature and the tumor immune microenvironment (TIME) using an external validation via the Gene Expression Omnibus (GEO) database was also conducted.

A workflow diagram is shown in Figure S1.

Datasets and Preprocessing

The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/) was used to obtain the gene expression profiles of thyroid cancer (THCA) (510 tumor and 58 normal samples). Clinical data from TCGA were also retrieved, including age, gender, clinical stage, and survival. Additionally, gene expression profiles were extracted (N = 197) from seven GEO series (GSE6004 (N = 14), GSE33630 (N = 49), GSE3467 (N = 9), GSE35570 (N = 65), GSE29265 (N = 20), GSE60542 (N = 33), and GSE3678 (N = 7)), processed, and normalized in batches using the “sva” and “limma” packages in R language (version 4.1.2).

To reduce the statistical bias, the clinical data with missing follow-up or fewer than 30 days of follow-up were eliminated. Finally, 498 tumor samples were received with RNA expression data and clinical data, which were randomly divided into the training and test sets with no significant variations in clinical characteristics (Table 1). Furthermore, 67 NRGs were obtained from previous reviews (Additional Files 1).

Differential gene expression analysis

The “limma” package in R was employed to identify differentially expressed genes (DEGs) in all tumor and adjacent normal tissues, based on FDR < 0.05.

Construction of the protein-protein interaction (PPI) network and Functional enrichment of DEGs

A PPI network was constructed to investigate the DEG interactions using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database (http://www.string-db.org/). For PPI study, the minimum required interaction score was set at 0.9 (The highest confidence). Then, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis using Metascape (https://metascape.org), a comprehensive tool for gene functional analysis was conducted.

Consensus clustering based on DEGs

Tumor samples were clustered into distinct subgroups based on necroptosis-related DEGs using R software “ConsensusClusterPlus.” To investigate the overall survival of subgroups, the “survival” and “survminer” packages were employed. The “limma” and “heatmap” packages were used to create a heatmap illustrating the differential expression and association between NRGs and clinical characteristics within each subgroup.

Construction of necroptosis-related prognostic signature

A training set to identify prognostic DEGs was used to create a risk signature for prognosis. The test set and total set were utilized to validate the predictive potential of created signature. Using the “glmnet” package, the prognostic significance of each DEG was determined using univariate Cox regression analysis. To avoid omissions, we established a cut-off value p < 0.2. To avoid model overfitting, the least absolute shrinkage and selection operator (LASSO) penalized Cox proportional hazards regression to identify additional important prognostic genes for overall survival (OS) was used. Then, the detected genes were combined into a multivariate Cox regression model, and the risk scores for each sample were generated using the following formula: risk score = esum (each necroptosis-related gene expression level × corresponding coefficient). Samples were categorized into high- and low-risk categories based on the computed median risk score. To assess predictive capacity, Kaplan–Meier survival curves were examined using the “survival” and “survminer” packages, and receiver operating characteristic (ROC) curves were created using the “timeROC” package. The “stats” package was used to perform principal component analysis (PCA).

Finally, univariate and multivariate Cox regression analysis was done to determine whether the risk classification derived from the risk signature is an independent prognostic factor.

Functional enrichment and GSEA enrichment analyses

DEGs were found between high- and low-risk groups using FDR < 0.05 and |log₂FC| ≥ 0.585. According to DEGs, the “clusterProfiler” package in R to conduct GO and KEGG enrichment studies and the “ggplot2” tool to create pictures was used, while GSEA 4.2.1 was used to discover related functions and pathway variations in the Hallmark gene set “h.all.v7.4.symbols.gmt.”

Analysis of tumor microenvironment and immune checkpoints

We started by calculating the scores for 16 immune infiltrating cells and 13 immune-related activities in each sample in both high- and low-risk groups using single-sample Gene Set Enrichment Analysis (ssGSEA) and R package “GSVA.” Then, using the Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data (ESTIMATE) program in conjunction with the “estimate” package, the immune-score, stromal-score, and ESTIMATES-score of each sample were determined. Finally, the differences in the scores mentioned above between the two groups were compared. Using the “ggpubr” package, the immune checkpoint activation between high- and low-risk groups was also compared.

The same formula to divide the normalized GEO data into high and low-risk groups based on the median risk score was used. The methods described above are used to externally validate the link between the signature and TIME.

Quantitative PCR

The First Affiliated Hospital of China Medical University provided 48 matched tumorous and non-tumorous tissue specimens of PTC. The following table summarizes the clinicopathological characteristics of 48 PTC patients treated at our hospital (Table 2). Total RNA was extracted from tissue samples using RNAiso (Takara, Dalian, China), and then reverse transcribed into cDNA using the QuantiTect Reverse Transcription Kit (Takara, Shiga, Japan). To validate gene expression, quantitative real-time PCR (qRT-PCR) analysis was done using TB-Green (Takara, Shiga, Japan), and the level of GAPDH served as an internal reference. The comparative Ct (2^−ΔΔCt) method was used to calculate the relative expression. The following primer sequences are available (Table 3).

Identification of necroptosis-related DEGs

A total of 50 differentially expressed genes between 510 tumors and 58 adjacent nontumor tissues were identified. Out of 67 necroptosis-related genes, 21 were up-regulated and 29 were downregulated in tumors. A heatmap was drawn for showing the 50 genes expression levels (Figure 1A). A PPI network was constructed to show the relationships among all necroptosis-related genes (Figure 1B). GO and KEGG analysis were used to investigate the potential roles of these DEGs, where GO enrichment analysis revealed that these necroptosis-related DEGs are mostly engaged in programmed necrotic cell death, transcription factor activity regulation, and cytokine response. According to KEGG analysis, these genes were primarily involved in cancer and apoptotic pathways (Figure 1C).

Tumor classification based on DEGs

To investigate the association between necroptosis-related DEGs and PTC subtypes, consensus clustering analysis on the 50 DEGs using the “ConsensusClusterPlus” program was performed. The letter “k” denoted the number of clusters. When k was adjusted to 2, intragroup correlations were highest while the intergroup correlations were lowest, showing that the 498 PTC patients were classified into two clusters (Figure 2A), C1 (N = 361) and C2 (N = 137). Additionally, the principal component analysis demonstrated that these patients can be classified into two distinct groups (Figure 2B). However, survival analysis revealed no statistically significant difference in OS between the two clusters (Figure 2C). Following that, a heatmap was created to depict the association between clusters, gene expression profiles, and clinical characteristics (Figure 2D), however, it was discovered that there are minimal clinical distinctions between the two clusters.

Development of a prognostic gene model

The models were constructed using the train set. To begin, a univariate Cox regression analysis was used to find survival-related genes. The results indicated that 16 genes in the training set satisfying a threshold with a p < 0.2 were associated with survival (Figure 2E). More precisely, 12 genes (MAPK8, IPMK, CYLD, AXL, ID1, CDKN2A, BCL2, ATRX, SPATA2, KLF9, BNIP3, and BCL2L11) were detrimental with HR > 1, whereas four genes (FADD, FAS, TNFRSF1B, and MPG) were protective with HR < 1. The overfitting was reduced by using LASSO Cox regression on the optimum λ value and identified 13 genes (FADD, FAS, IPMK, TNFRSF1B, CYLD, AXL, ID1, CDKN2A, MPG, SPATA2, KLF9, BNIP3, and BCL2L11) that were substantially linked with prognosis (Figure 2F & G). Further study of these genes using multivariate Cox regression analysis resulted in the construction of a prognostic signature consisting of seven necroptosis-related genes. For each sample, a risk score was calculated using the following formula: Risk score = FAS*(-1.1721) + IPMK*(2.5404) +TNFRSF1B*(0.9750) + AXL*(1.5358) + CDKN2A*(1.6495) + SPATA2*(4.1484) + KLF9*(2.0030). The patients were classified in the training set into high- and low-risk groups based on their risk score median value and then ranked and assessed their risk score distributions (Figure 3A). Additionally, the survival status of each sample in the training set was visualized using dot plots and discovered that all deceased patients were concentrated in the high-risk group (Figure 3B). A heatmap was created to depict the differential expression of the seven prognostic genes between the two groups (Figure 3C). Survival analysis revealed that patients with high-risk scores had a significantly worse outcome than those with low-risk scores (p = 0.002) (Figure 3J). ROC analysis was utilized to determine the sensitivity and specificity of created prognostic signature, and the areas under the curve (AUC) were found to be 0.834 at one year, 0.938 at three years, and 0.960 at five years respectively (Figure 3M). The PCA revealed that patients with varying risks were dispersed in two directions (Figure 3P).

Following that, the predictive ability of prognostic signature in the test and total sets was validated. The test and total sets were separated into high- and low-risk groups, respectively, based on the median risk score in the training set. The risk score distribution, survival status, and expression of seven necroptosis-related genes are provided below for each patient in the test set and total set (Figure 3D,E,F; Figure 3G, H, I). In both the test and total sets, Kaplan-Meier survival curves indicated that high-risk patients had a shorter OS than low-risk patients (Figure 3K; 3L). In the test set, the one-year AUC was found to be 0.703, the three-year AUC was 0.656, and the five-year AUC was 0.897 in the test set respectively (Figure 3N). In the total set, the AUC for the one-year was found to be 0.795, the AUC for the three years was 0.822, and the AUC for the five years was 0.902. (Figure 3O). PCA verified the separation of patients into two clusters according to their risk status (Figure 3Q; 3R).

Independent prognostic value of the risk signature

The univariate Cox regression analysis revealed that the factors affecting the prognosis were age, clinical stage, and risk grouping (Figure 4A). After adjusting for other confounding factors, it was discovered that age and risk grouping were independent predictors of OS in patients with PTC by multivariate Cox regression analysis (Figure 4B). Furthermore, the ROC curves of risk grouping and other clinicopathological variables showed that risk grouping had a higher prediction level, and its prediction ability increased with time (Figure 4C-E). Furthermore, survival analysis revealed that individuals with high-risk ratings had lower overall survival in all clinical subgroups (Figure S2A-L).

Functional analyses based on the risk signature

A heatmap was created to depict the link between predictive gene expression and clinical features (Figure 5A). Then, using FDR < 0.05 and |log₂FC| ≥ 0.585, 268 differentially expressed genes were identified in high- and low-risk groups. Based on DEGs, GO and KEGG analysis was performed to further elucidate the putative biological activities and pathways associated with the risk signature. According to findings, these DEGs were primarily enriched in cell adhesion, immune-related activities, and several infections or immune-related pathways (Figure 5B, C). The GSEA was used to examine the transcript messages of PTC patients who were grouped by risk score into high- and low-risk categories. According to findings, most NRGs prognostic signatures regulated immunological and tumor-related pathways. In the low-risk group, biological pathways such as interferon-alpha and gamma response, allograft rejection, P53 pathway, inflammatory response, and IL6-JAK-STAT3 signaling were enriched. (Figure 5D).

Comparison of the immune activity among subgroups

The ssGSEA was used to compare the enrichment scores of 16 different types of immune cells and the activity of 13 immune-related functions in low- and high-risk groups. The findings revealed that, except for CD8⁺T cells, which did not differ substantially between the two groups, other immune cells infiltrate at a higher rate in the low-risk subgroup (Figure 6A). Furthermore, the great majority of immune-related pathways were activated more in the low-risk group than in the high-risk group (Figure 6B). Each sample's immune-score, stromal-score, and ESTIMATES-score were calculated to determine the estimated amount of stromal and immune cells in the tumor sample. It was observed that the low-risk group contained more immune cells than the high-risk group, but the difference in stromal cell composition between the two groups was not statistically significant. Furthermore, ESTIMATES-score was greater in the low-risk group, indicating reduced tumor purity (Figure 6C). The potential differences in immune checkpoint expression between the two groups were investigated, and the results revealed that most immune checkpoints were more expressed in the low-risk group (Figure 6D). Similar outcomes were reached when analyzing the immunological state of GEO cohort. The low-risk group had more immune infiltration and immune checkpoint expression than the high-risk group (Figure 7).

The expression levels of seven prognostic genes

Following that, bioinformatics analysis was performed of mRNA expression data from GEO database and qRT-PCR experiments to validate mRNA expression profiles using 48 paired tumor and normal tissues from PTC patients at China Medical University, in addition to TCGA database, to validate the expression profile of seven prognostic genes in both tumor and normal tissues. The qRT-PCR results revealed that AXL (p = 0.0308), FAS (p = 0.0134), TNFRSF1B (p = 0.0250), and CDKN2A (p = 0.0224) mRNA expression were considerably greater in tumor tissues. In contrast, tumor samples had lower levels of IPMK (p = 0.0010), SPATA2 (p = 0.0071), and KLF9 (p = 0.0002). The expression of these genes was consistent with TCGA and GEO datasets, which adds to the authenticity and confidence of our developed signature (Figure 8).

Thyroid cancer is the most frequent endocrine malignancy, with PTC accounting for more than 85% of all follicular-derived well-differentiated thyroid tumors (23). Even though most PTCs are well differentiated, with a low risk of local invasion, recurrence, or metastasis (24), some patients with more advanced or aggressive variations frequently exhibit heterogeneity with distinct clinical, pathological, and molecular features (25, 26). According to the most recent American Thyroid Association (ATA) guidelines, these pathological subtypes confer an intermediate risk of recurrence (27), and these variants are associated with higher rates of recurrence and metastasis, as well as, in some cases, minimal efficacy of radioiodine therapy and may have inferior survival (28). Due to a lack of understanding of the natural history of these more aggressive variants, treatment for these patients is frequently inadequate or suboptimal. Therefore, carrying out risk stratification management is critical for further optimizing the diagnosis and treatment of PTC patients.

Necroptosis is a type of programmed cell death that can address apoptosis resistance while also activating and enhancing antitumor immunity in cancer treatment (29), which can be triggered by cytokines, danger signals, or pathogen infection (30). After being induced, complex upstream signals eventually converge on the activation of RIPK3, which binds and directly phosphorylates MLKL, a crucial effector of necroptosis and the pathway's most downstream component known to be necessary (31). In contrast to apoptosis, this type of cell death is distinguished by plasma membrane permeabilization and the production of damage-associated molecular patterns (DAMPs), which have significant immunological effects (11). Necroptotic cells may supply antigens and inflammatory cytokines to dendritic cells (DCs) for antigen cross-priming, which activates cytotoxic CD8⁺ T lymphocytes, resulting in tumor cell eradication through the release of DAMPs into the tissue milieu (10). However, necroptotic cells may recruit immune inflammatory cells and induce inflammation, which can promote tumor development by boosting angiogenesis, cancer cell proliferation, and metastasis. A study showed that human and murine tumor cells stimulate endothelial cells to undergo programmed necrosis (necroptosis), which facilitates tumor cell extravasation and metastasis. RIPK1 inhibitor necrostatin-1 or endothelial-cell-specific RIPK3 deletion inhibited tumor cell-induced endothelial necroptosis, tumor cell extravasation, and metastasis in mice (15). As can be observed, the role of necrotic apoptosis in tumors is extremely complex. The potential significance of necroptosis-related genes in PTC, on the other hand, is uncertain. As a result, an attempt was made to identify potential diagnostic markers of necroptosis using targeted immunotherapy to increase the survival of PTC patients and to investigate the prognostic and diagnostic relevance of necroptosis. This study’s findings imply that inducing necroptosis with immunotherapy may be a viable treatment strategy for improving patient outcomes.

In this work, we looked at the expression of 67 necroptosis-related genes in PTC samples and adjacent samples and discovered that 50 of them were expressed differently, where 21 of these genes were up-regulated, and 29 were down-regulated. A consistent cluster analysis on PTC samples was conducted based on NRG expression and discovered that PTC samples can be separated into two subgroups, implying that NRGs may be involved in PTC subclassification. Following that, seven genes related to prognosis (IPMK, AXL, CDKN2A, SPATA2, KLF9, FAS, and TNFRSF1B) were obtained and constructed a seven genes prognostic risk signature from NRGs using a series of processes such as univariate Cox regression analysis, LASSO regression analysis, and multivariate Cox regression analysis. Inositol polyphosphate multikinase (IPMK) is a conserved enzyme that belongs to the inositol phosphokinase 6-kinase family and plays a key role in the phosphorylation of inositol phosphates (IPs) (32, 33). During necroptosis, IPMK is also required for activated phospho-MLKL to oligomerize and relocate to the plasma membrane. IPMK, in particular, may generate highly phosphorylated IPs in cells from lowly phosphorylated precursors, and highly phosphorylated IPs can influence necroptosis by directly binding MLKL and controlling its function (34). It was discovered that IPMK expression was lower in PTC tumor tissues than in normal tissues, but that high expression was linked with a poor prognosis in PTC. AXL is a tyrosine kinase receptor that belongs to TAM (TYRO3, AXL, and MER) receptor family. TAM family regulates cell growth, survival, and proliferation (35, 36). AXL is frequently overexpressed in cancer, which can increase the invasiveness and motility of the normally non-invasive MCF-7 cell line (37). AXL has been demonstrated to boost the motility, invasion, proliferation, and survival of breast cancer cells (38). Importantly, a study evaluating transcriptome analysis of 941 cancer cell lines found that high AXL mRNA levels indicate resistance to necroptosis while low RIPK3 mRNA levels predict resistance to necroptosis (39). AXL behaved as a cancer-promoting gene in our sample validation due to its overexpression in tumor tissues and its negative connection with survival time. Although cyclin-dependent kinase inhibitor 2A (CDKN2A) is frequently altered or deleted in a wide range of cancers and is known to be an important tumor suppressor gene, a study suggested that CDKN2A plays a key role in the formation and progression of larynx squamous cell carcinoma (40). Furthermore, CDKN2A hypermethylation may be a risk factor for a poor prognosis of pancreatic cancer (41). CDKN2A expression differed considerably between tumor and normal tissues in our investigation. It was more abundant in tumors, and increased expression in malignancies is associated with a poor prognosis. A genome-wide siRNA screen identified Spermatogenesis-associated 2 (SPATA2) as a gene implicated in necroptosis mediation (42). SPATA2 modulated RIPK1 activation via altering its M1 ubiquitination. More particular, SPATA2 deficiency can hasten the M1 ubiquitination of RIPK1, resulting in necroptosis resistance (43). According to research, increased SPATA2 expression is related to a poor prognosis in various cancers, such as ovarian or cervical cancer (44, 45). We found high SPATA2 expression in tumor tissues from patients with a relatively poor prognosis. Kruppel-like factor 9 (KLF9), also known as basic transcription element-binding protein 1 (BTE-B1), is a transcriptional regulator involved in cellular adhesion, differentiation, and proliferation. It is reported to be downregulated in various cancers, including endometrial carcinoma and colorectal cancer (46, 47). Sadia et al. discovered that KLF9 expression was significantly lower in cervical cancer patients compared to healthy controls and that it was significantly lower in the advanced tumor stage and distant metastatic groups compared to the lower tumor stage and non-metastatic groups (48). However, Zizhen Chen et al. discovered that KLF9 expression was positively associated with acute myeloid leukemia (49). Consequently, KLF9 may have a distinct role in various cancer types. KLF9 appeared to be a cancer suppressor gene in our investigation, as it was downregulated threefold in tumor tissues; yet, it also contributed to patient survival shortening since it was enriched in the high-risk group. Fas cell surface death receptor (FAS) has been found to play a critical role in the physiological regulation of programmed cell death as a member of TNF-receptor superfamily and has been implicated in the pathogenesis of several malignancies and immune system diseases (50). FAS mRNA and protein expression levels were considerably lowered in breast carcinoma, whereas high FAS expression implies a better prognosis in breast cancer patients (51). FAS was demonstrated to be increased in tumor samples in our investigation, and its high expression in PTC patients predicted better survival, implying that FAS may be a tumor suppressor gene. Tumor necrosis factor receptor-2 (TNFR2/TNFRSF1B) is a cell-surface receptor that regulates cell survival and proliferation and is abundantly expressed on the surface of many human tumors (52). Recently, the possibility of targeting this receptor as a next-generation cancer therapy method has emerged (53, 54). It was discovered that the gene was more abundant in tumor tissues and that increased expression was associated with a better prognosis.

The risk score for each PTC patient was estimated using the risk signature and classified into high- and low-risk groups based on the median value of all patient risk scores. Kaplan–Meier survival curves revealed that high-risk patients had a shorter OS than low-risk individuals. Additionally, the AUC of the receiver operating characteristic curve demonstrated that the risk signature was effective at predicting survival prognosis. Cox regression analysis, both univariate and multivariate, established that the created risk signature was an independent risk factor for prognosis. Through the above analysis, it was proved that the signature has a good ability to predict prognosis. Functional enrichment analysis was also performed on differentially expressed genes across high- and low-risk groups and discovered that these genes were primarily involved in cell adhesion, immunological function, and pathogen infection. GSEA revealed a high enrichment of immune-related pathways in the low-risk categories. As a result, it is hypothesized that necroptosis may be associated with the immunological microenvironment of PTC. The calculated scores for immune cell infiltration and immunological function in high- and low-risk groups using ssGSEA revealed that the patients in the low-risk group had a higher total immunological activity. Additionally, ESTIMATE algorithm was employed to estimate the tumor purity of each sample. To be more precise, the greater the immune- and stromal-scores, the less pure the tumor. This indicated that the high-risk group has a higher tumor purity, which is frequently associated with a poor prognosis (55). Cancer immunotherapy has experienced remarkable advances in recent years and more and more potential immune checkpoints have been found (56, 57). As a result, the degree of expression of common immune checkpoint proteins in high- and low-risk groups was also investigated, which revealed that the great majority of immune checkpoint proteins were expressed at a higher level in the low-risk group, implying that immunotherapy may be more successful for low-risk patients.

Since there are few studies on the involvement of necroptosis in thyroid cancer at the moment, this study gives some theoretical foundations and research ideas for the way forward. However, this research encountered certain limitations. To begin, due to a lack of sample survival data in GEO database, the prediction ability of the signature to prognosis is only internally validated in TCGA database, and hence were unable to obtain sufficient data from other sources to externally validate the signature-prognosis association. Furthermore, additional fundamental investigations are required to study the relationship between the signature and the tumor microenvironment. As a result, additional research is required to substantiate the above-mentioned conclusions.

In general, this study revealed the expression and prognostic value of NRGs in PTC and developed a risk prediction signature of necroptosis-related genes. The created signature contributes to prognostic prediction as well as immunotherapy evaluation in PTC patients. These results may provide ideas for further investigation of the role of necroptosis in PTC and guide clinicians towards individualized treatments for PTC patients with different clinical characteristics.

LASSO: least absolute shrinkage and selection operator

ssGSEA: single-sample gene set enrichment analysis

PTC: papillary thyroid cancer

AUC: area under the curve

ROC: receiver operator characteristic

US: ultrasonography

NRGs: necroptosis-related genes

TIME: tumor immune microenvironment

GEO: Gene Expression Omnibus

TCGA: The Cancer Genome Atlas

THCA: thyroid cancer

DEGs: differentially expressed genes

PPI: protein-protein interaction

STRING: the Search Tool for the Retrieval of Interacting Genes/Proteins database

GO: Gene Ontology

KEGG: Kyoto Encyclopedia of Genes and Genomes

OS: overall survival

PCA: principal component analysis

Acknowledgements

We would like to acknowledge the TCGA and the GEO (GSE6004, GSE33630, GSE3467, GSE35570, GSE29265, GSE60542 and GSE3678) network for providing data.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Authors’ contributions

Z.Wang performed the data analysis and wrote the manuscript. P.Wu and J.Shi performed the experiments and interpretated results. X.Ji helped to interpretated results and write the manuscript. H.Zhang and W.Sun reviewed and revised the manuscript. All authors read and approved the final manuscript.

Ethics approval and consent to participate

The research was approved by the Institutional Research Ethics Committees of the First Affiliated Hospital of China Medical University. Informed consent for publication was obtained from all patients for collection of tissue samples prior to the surgery.

Funding

This work was supported by the National Natural Science Foundation of China (81902726), the China Postdoctoral Science Foundation (2018M641739).

Availability of data and materials

Gene expression profiles and clinical information of THCA and gene expression profiles of GEO in this study are available from the public database (TCGA, https://portal.gdc. cancer.gov/; GEO, https://www.ncbi.nlm.nih.gov/geo/).

Sugitani I, Ito Y, Takeuchi D, Nakayama H, Masaki C, Shindo H, et al. Indications and Strategy for Active Surveillance of Adult Low-Risk Papillary Thyroid Microcarcinoma: Consensus Statements from the Japan Association of Endocrine Surgery Task Force on Management for Papillary Thyroid Microcarcinoma. Thyroid. 2021;31(2):183-92.
Ito Y, Miyauchi A, Kihara M, Fukushima M, Higashiyama T, Miya A. Overall Survival of Papillary Thyroid Carcinoma Patients: A Single-Institution Long-Term Follow-Up of 5897 Patients. World J Surg. 2018;42(3):615-22.
Kitahara CM, Sosa JA. The changing incidence of thyroid cancer. Nat Rev Endocrinol. 2016;12(11):646-53.
Lim H, Devesa SS, Sosa JA, Check D, Kitahara CM. Trends in Thyroid Cancer Incidence and Mortality in the United States, 1974-2013. JAMA. 2017;317(13):1338-48.
Albareda M, Puig-Domingo M, Wengrowicz S, Soldevila J, Matias-Guiu X, Caballero A, et al. Clinical forms of presentation and evolution of diffuse sclerosing variant of papillary carcinoma and insular variant of follicular carcinoma of the thyroid. Thyroid. 1998;8(5):385-91.
Carling T, Ocal IT, Udelsman R. Special variants of differentiated thyroid cancer: does it alter the extent of surgery versus well-differentiated thyroid cancer? World J Surg. 2007;31(5):916-23.
Feng J, Shen F, Cai W, Gan X, Deng X, Xu B. Survival of aggressive variants of papillary thyroid carcinoma in patients under 55 years old: a SEER population-based retrospective analysis. Endocrine. 2018;61(3):499-505.
Degterev A, Huang Z, Boyce M, Li Y, Jagtap P, Mizushima N, et al. Chemical inhibitor of nonapoptotic cell death with therapeutic potential for ischemic brain injury. Nat Chem Biol. 2005;1(2):112-9.
Christofferson DE, Yuan J. Necroptosis as an alternative form of programmed cell death. Curr Opin Cell Biol. 2010;22(2):263-8.
Gong Y, Fan Z, Luo G, Yang C, Huang Q, Fan K, et al. The role of necroptosis in cancer biology and therapy. Molecular cancer. 2019;18(1):100.
Pasparakis M, Vandenabeele P. Necroptosis and its role in inflammation. Nature. 2015;517(7534):311-20.
Seehawer M, Heinzmann F, D'Artista L, Harbig J, Roux P-F, Hoenicke L, et al. Necroptosis microenvironment directs lineage commitment in liver cancer. Nature. 2018;562(7725):69-75.
Stoll G, Ma Y, Yang H, Kepp O, Zitvogel L, Kroemer G. Pro-necrotic molecules impact local immunosurveillance in human breast cancer. Oncoimmunology. 2017;6(4):e1299302.
Seifert L, Werba G, Tiwari S, Giao Ly NN, Alothman S, Alqunaibit D, et al. The necrosome promotes pancreatic oncogenesis via CXCL1 and Mincle-induced immune suppression. Nature. 2016;532(7598):245-9.
Strilic B, Yang L, Albarrán-Juárez J, Wachsmuth L, Han K, Müller UC, et al. Tumour-cell-induced endothelial cell necroptosis via death receptor 6 promotes metastasis. Nature. 2016;536(7615):215-8.
Feng X, Song Q, Yu A, Tang H, Peng Z, Wang X. Receptor-interacting protein kinase 3 is apredictor of survival and plays atumor suppressive role in colorectal cancer. Neoplasma. 2015;62(4):592-601.
Geserick P, Wang J, Schilling R, Horn S, Harris PA, Bertin J, et al. Absence of RIPK3 predicts necroptosis resistance in malignant melanoma. Cell death & disease. 2015;6:e1884.
Koo G-B, Morgan MJ, Lee D-G, Kim W-J, Yoon J-H, Koo JS, et al. Methylation-dependent loss of RIP3 expression in cancer represses programmed necrosis in response to chemotherapeutics. Cell Res. 2015;25(6):707-25.
McCormick KD, Ghosh A, Trivedi S, Wang L, Coyne CB, Ferris RL, et al. Innate immune signaling through differential RIPK1 expression promote tumor progression in head and neck squamous cell carcinoma. Carcinogenesis. 2016;37(5):522-9.
Park S, Hatanpaa KJ, Xie Y, Mickey BE, Madden CJ, Raisanen JM, et al. The receptor interacting protein 1 inhibits p53 induction through NF-kappaB activation and confers a worse prognosis in glioblastoma. Cancer research. 2009;69(7):2809-16.
Wu Z, Huang X, Cai M, Huang P, Guan Z. Novel necroptosis-related gene signature for predicting the prognosis of pancreatic adenocarcinoma. Aging (Albany NY). 2022;14(2):869-91.
Li X-Y, You J-X, Zhang L-Y, Su L-X, Yang X-T. A Novel Model Based on Necroptosis-Related Genes for Predicting Prognosis of Patients With Prostate Adenocarcinoma. Frontiers in bioengineering and biotechnology. 2021;9:814813.
Wiltshire JJ, Drake TM, Uttley L, Balasubramanian SP. Systematic Review of Trends in the Incidence Rates of Thyroid Cancer. Thyroid. 2016;26(11):1541-52.
Cabanillas ME, McFadden DG, Durante C. Thyroid cancer. Lancet. 2016;388(10061):2783-95.
Chen B, Shi Y, Xu Y, Zhang J. The predictive value of coexisting BRAFV600E and TERT promoter mutations on poor outcomes and high tumour aggressiveness in papillary thyroid carcinoma: A systematic review and meta-analysis. Clin Endocrinol (Oxf). 2021;94(5):731-42.
Póvoa AA, Teixeira E, Bella-Cueto MR, Melo M, Oliveira MJ, Sobrinho-Simões M, et al. Clinicopathological Features as Prognostic Predictors of Poor Outcome in Papillary Thyroid Carcinoma. Cancers (Basel). 2020;12(11).
Haugen BR, Alexander EK, Bible KC, Doherty GM, Mandel SJ, Nikiforov YE, et al. 2015 American Thyroid Association Management Guidelines for Adult Patients with Thyroid Nodules and Differentiated Thyroid Cancer: The American Thyroid Association Guidelines Task Force on Thyroid Nodules and Differentiated Thyroid Cancer. Thyroid. 2016;26(1).
Lam KY, MD, FRCPA, Lo CY, MS, FRCS, et al. Papillary carcinoma of thyroid: A 30-yr clinicopathological review of the histological variants. Endocrine Pathology. 2005;16(4):323-30.
Kowalski S, Hać S, Wyrzykowski D, Zauszkiewicz-Pawlak A, Inkielewicz-Stępniak I. Selective cytotoxicity of vanadium complexes on human pancreatic ductal adenocarcinoma cell line by inducing necroptosis, apoptosis and mitotic catastrophe process. Oncotarget. 2017;8(36):60324-41.
Linkermann A, Green DR. Necroptosis. N Engl J Med. 2014;370(5):455-65.
Zhao J, Jitkaew S, Cai Z, Choksi S, Li Q, Luo J, et al. Mixed lineage kinase domain-like is a key receptor interacting protein 3 downstream component of TNF-induced necrosis. Proc Natl Acad Sci U S A. 2012;109(14):5322-7.
Irvine RF, Schell MJ. Back in the water: the return of the inositol phosphates. Nature reviews Molecular cell biology. 2001;2(5):327-38.
Otto JC, Kelly P, Chiou S-T, York JD. Alterations in an inositol phosphate code through synergistic activation of a G protein and inositol phosphate kinases. Proc Natl Acad Sci U S A. 2007;104(40):15653-8.
Dovey CM, Diep J, Clarke BP, Hale AT, McNamara DE, Guo H, et al. MLKL Requires the Inositol Phosphate Code to Execute Necroptosis. Molecular cell. 2018;70(5).
Lemke G. Biology of the TAM receptors. Cold Spring Harb Perspect Biol. 2013;5(11):a009076.
Lemke G, Rothlin CV. Immunobiology of the TAM receptors. Nat Rev Immunol. 2008;8(5):327-36.
Zhang Y-X, Knyazev PG, Cheburkin YV, Sharma K, Knyazev YP, Orfi L, et al. AXL is a potential target for therapeutic intervention in breast cancer progression. Cancer research. 2008;68(6):1905-15.
Gjerdrum C, Tiron C, Høiby T, Stefansson I, Haugen H, Sandal T, et al. Axl is an essential epithelial-to-mesenchymal transition-induced regulator of breast cancer metastasis and patient survival. Proc Natl Acad Sci U S A. 2010;107(3):1124-9.
Najafov A, Zervantonakis IK, Mookhtiar AK, Greninger P, March RJ, Egan RK, et al. BRAF and AXL oncogenes drive RIPK3 expression loss in cancer. PLoS biology. 2018;16(8):e2005756.
Sasiadek MM, Stembalska-Kozlowska A, Smigiel R, Ramsey D, Kayademir T, Blin N. Impairment of MLH1 and CDKN2A in oncogenesis of laryngeal cancer. Br J Cancer. 2004;90(8):1594-9.
Xing X, Cai W, Shi H, Wang Y, Li M, Jiao J, et al. The prognostic value of CDKN2A hypermethylation in colorectal cancer: a meta-analysis. Br J Cancer. 2013;108(12):2542-8.
Hitomi J, Christofferson DE, Ng A, Yao J, Degterev A, Xavier RJ, et al. Identification of a molecular signaling network that regulates a cellular necrotic cell death pathway. Cell. 2008;135(7):1311-23.
Wei R, Xu LW, Liu J, Li Y, Zhang P, Shan B, et al. SPATA2 regulates the activation of RIPK1 by modulating linear ubiquitination. Genes & development. 2017;31(11):1162-76.
Wieser V, Abdel Azim S, Sprung S, Knoll K, Kögl J, Hackl H, et al. TNFα signalling predicts poor prognosis of patients with endometrial cancer. Carcinogenesis. 2020;41(8):1065-73.
Wieser V, Tsibulak I, Degasper C, Welponer H, Leitner K, Parson W, et al. Tumor necrosis factor receptor modulator spermatogenesis-associated protein 2 is a novel predictor of outcome in ovarian cancer. Cancer science. 2019;110(3):1117-26.
Kang L, Lü B, Xu J, Hu H, Lai M. Downregulation of Krüppel-like factor 9 in human colorectal cancer. Pathol Int. 2008;58(6):334-8.
Simmen FA, Su Y, Xiao R, Zeng Z, Simmen RCM. The Krüppel-like factor 9 (KLF9) network in HEC-1-A endometrial carcinoma cells suggests the carcinogenic potential of dys-regulated KLF9 expression. Reprod Biol Endocrinol. 2008;6:41.
Safi S, Badshah Y, Shabbir M, Zahra K, Khan K, Dilshad E, et al. Predicting 3D Structure, Cross Talks, and Prognostic Significance of in Cervical Cancer. Front Oncol. 2021;11:797007.
Chen Z, Song J, Wang W, Bai J, Zhang Y, Shi J, et al. A novel 4-mRNA signature predicts the overall survival in acute myeloid leukemia. Am J Hematol. 2021;96(11):1385-95.
Vandenabeele P, Galluzzi L, Vanden Berghe T, Kroemer G. Molecular mechanisms of necroptosis: an ordered cellular explosion. Nature reviews Molecular cell biology. 2010;11(10):700-14.
Zhang Y, Shao X, Gao C, Xu D, Wu J, Zhu X, et al. High FAS expression correlates with a better prognosis and efficacy of taxanes and target regents in breast cancer. Cancer biomarkers : section A of Disease markers. 2021;32(2):207-19.
Chen X, Bäumel M, Männel DN, Howard OMZ, Oppenheim JJ. Interaction of TNF with TNF receptor type 2 promotes expansion and function of mouse CD4+CD25+ T regulatory cells. J Immunol. 2007;179(1):154-61.
Chen X, Oppenheim JJ. Targeting TNFR2, an immune checkpoint stimulator and oncoprotein, is a promising treatment for cancer. Sci Signal. 2017;10(462).
Vanamee ÉS, Faustman DL. TNFR2: A Novel Target for Cancer Immunotherapy. Trends in molecular medicine. 2017;23(11):1037-46.
Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nature communications. 2013;4:2612.
Charoentong P, Finotello F, Angelova M, Mayer C, Efremova M, Rieder D, et al. Pan-cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade. Cell Rep. 2017;18(1):248-62.
van den Bulk J, Verdegaal EM, de Miranda NF. Cancer immunotherapy: broadening the scope of targetable tumours. Open Biol. 2018;8(6).

Table 1

The clinical characteristics in training, test and total sets

Variables	Group	Training set(N=249)	Test set(N=249)	Total set(N=498)	P value
Age(year)	≤60	196 (78.71%)	190 (76.31%)	386 (77.51%)	0.520
	＞60	53 (21.29%)	59 (23.69%)	112 (22.49%)
Gender, n (%)	Male	62 (24.90%)	73 (29.32%)	135 (27.11%)	0.267
	Female	187 (75.10%)	176 (70.68%)	363 (72.89%)
Stage, n (%)	Stage I-II	164 (65.86%)	167 (67.07%)	331 (66.47%)	0.362
	Stage III-IV	83 (33.33%)	82 (32.93%)	165 (33.13%)

	Unknow	2 (0.80%)	0 (0.00%)	2 (0.40%)
T, n (%)	T1-2	143 (57.43%)	161 (64.66%)	304 (61.04%)	0.252
	T3-4	105 (42.17%)	87 (34.94%)	192 (38.55%)
	Tx/unknow	1 (0.40%)	1 (0.40%)	2 (0.40%)
M, n (%)	M0	152 (61.04%)	130 (52.21%)	282 (56.63%)	0.112
	M1	5 (2.01%)	4 (1.61%)	9 (1.81%)
	Mx/unknow	92 (36.95%)	115 (46.18%)	207 (41.57%)
N, n (%)	N0	105 (42.17%)	123 (49.40%)	228 (45.78%)	0.206
	N1	115 (46.18%)	105 (42.17%)	220 (44.18%)
	Nx/unknow	29 (11.65%)	21 (8.43%)	50 (10.04%)

Table 2

The clinicopathological features of PTC (N=48)

Clinical variables	Group	Sample (N=48)	Percentage (%)
Age(year), n (%)	<=60	44	91.67%
	>60	4	8.33%
Gender, n (%)	Female	38	79.17%
	Male	10	20.83%
Stage, n (%) T, n (%) M, n (%)	Stage I-II Stage III-IV Unknow T1-2 T3-4 Tx/unknow M0 M1 Mx/unknow	40 8 0 33 15 0 47 1 0	83.33% 16.67% 0.00% 68.75% 31.25% 0.00% 97.92% 2.08% 0.00%
N, n (%)	N0	14	29.17%
	N1	34	70.83%
	Nx/unknow	0	0.00%

Table 3

Premier sequences for qRT-PCR analysis

Premier	Sequences (5′–3′)
IPMK-F	GTGCTTGGCATGAGGGTTTATC
IPMK-R	TGGCAGCAACAGCATCTTTTC
CDKN2A-F	GCGGAAGGTCCCTCAGAAATG
CDKN2A-R	GCCAGCTTGCGATAACCAAA
SPATA2-F	GACTTATTTCGGAAGTACGTGC
SPATA2-R	GATCAGCCGGAATCGATAAAAG
KLF9-F	GTGTCTGGTTTCCATTTCGAAC
KLF9-R	GATCCCATATCCTCATCTGGAC
TNFRSF1B-F	CGGCTCAGAGAATACTATGACC
TNFRSF1B-R	ACAGAAGACTTTTGCATGTTGG
FAS-F	GTACACAGACAAAGCCCATTTT
FAS-R	TTTGGTTTACATCTGCACTTGG
AXL-F	AGATTTATGACTATCTGCGCCA
AXL-R	TGACATAGAGGATTTCGTCAGG
GAPDH-F	GTCTCCTCTGACTTCAACAGCG
GAPDH-R	ACCACCCTGTTGCTGTAGCCAA

No competing interests reported.

67NRGs.docx
workflowdiagram.png
Figure S1. The flowchart of the study.
FIGURES2.tif
Figure S2. The stratification analysis of clinical factors based on age, gender, TNM stage and clinical stage.

Download PDF

Version 1

posted

You are reading this latest preprint version

A novel necroptosis-related genes signature for predicting the prognosis of thyroid cancer

Status:

Version 1

Abstract

Figures

Background

Methods

Datasets and Preprocessing

Differential gene expression analysis

Construction of the protein-protein interaction (PPI) network and Functional enrichment of DEGs

Consensus clustering based on DEGs

Construction of necroptosis-related prognostic signature

Functional enrichment and GSEA enrichment analyses

Analysis of tumor microenvironment and immune checkpoints

Quantitative PCR

Results

Identification of necroptosis-related DEGs

Tumor classification based on DEGs

Development of a prognostic gene model

Independent prognostic value of the risk signature

Functional analyses based on the risk signature

Comparison of the immune activity among subgroups

The expression levels of seven prognostic genes

Discussion

Conclusion

Abbreviations

Declarations

References

Tables

Additional Declarations

Supplementary Files

Status:

Version 1