Background Information

Prostate adenocarcinoma (PRAD) is a complex but common malignancy that accounts for about 1300000 new cases and 360,000 deaths every year worldwide. Notably, PRAD accounts for 15% of all new tumor-related cases, making it the second most common neoplasia in elderly males, and the fifth most frequent cause of cancer-related deaths worldwide (Bray et al., 2018; Siegel et al., 2020). In China, PRAD is fast increasing, owing to the rapid socio-economic development in the country, coupled with changes in people’s living and eating habits, as well as an increase in the aging population. Consequently, the disease has become one of the most common urogenital malignancies among elderly Chinese males (Wei et al., 2020). Current treatment options for PRAD are limited, while patient prognosis remains unsatisfactory. Therefore, prospecting for novel prognostic markers is imperative to development of effective treatment strategies and enhanced prognosis of PRAD patients. Previous studies have shown that necroptosis, which was first identified and named in 2005, is regulated by intracellular signaling pathways (Degterev et al., 2005). Notably, this phenomenon can be unmediated by caspases, thus functioning upon inhibition of apoptotic pathways. Moreover, its cellular morphology is consistent with that of conventional necrosis. Additional research evidences have shown that necroptosis is not only an important mechanism underlying cell death, but also plays a crucial role in development and progression of many immune system diseases, including chronic kidney diseases, cerebral ischemia, myocardial ischemia, acute and chronic neurodegenerative diseases, as well as tumors and many other human pathological activities (Jouan-Lanhouet et al., 2014; Barbosa et al., 2018). Although necroptosis exerts important functions in oncogenesis and anticancer processes, only a handful of studies have described its significance in PRAD. In the present study, we systematically analyzed differential expression profiles of necroptosis-related genes between normal and PRAD tissues, then developed a novel risk-score-based model for predicting prognosis of PRAD patients.

Material and Methods

Data Acquisition and Differential Gene Expression Analysis

Necroptosis-related genes were extracted from previous studies (Fan et al., 2014; Frank and Vince, 2019; Gong et al., 2019; Yuan et al., 2019; Tang et al., 2020). RNA sequence dataset belonging to 499 PRAD patients and 52 normal controls, together with corresponding clinical information were accessed and downloaded from the TCGA database. Next, we employed the “limma” package implemented in R software to identify differentially expressed genes (DEGs) between the tumor and adjacent normal tissues, based on FDR <0.05 and |log2FC| ≥ 1. Thereafter, we recruited a total of 80 PRAD patients at The First Affiliated Hospital of Zhengzhou University, as the validation cohort. All patients voluntarily signed a written informed consent prior to inclusion in the study, and ethical approval for the study protocol was obtained from The First Affiliated Hospital of Zhengzhou University.

Construction and Validation of a Prognostic Model

The prognostic value of each DEG was first assessed by univariate Cox regression analysis, then genes that were significantly correlated with OS in PRAD patients identified. To avoid overfitting of the model, we performed Lasso regression analysis to further select significant prognostic genes for OS in PRAD patients, using a penalty parameter tuning (λ) that was conducted by 10-fold cross-validation based on minimum criteria. The identified significantly expressed genes were then incorporated into a multivariate Cox regression model, and the risk score of each patient calculated using the following formula: risk score = esum (each gene’s expression × corresponding coefficient). Median risk-scores were then used to stratify patients into either high- or low-risk groups, and validation of model feasibility and accuracy conducted by generating AUC of the ROC as well as calibration plots. Next, we applied the Cox proportional hazards regression model to analyze these risk factors in PRAD, targeting risk scores, gender, age, as well as the T, N, and M stages. To validate the established model, based on TCGA, we used the median scores to divide patients in the validation cohort into high- and low-risk groups. We also validated the model by stratifying patients in the validation cohort into low- and high-risk subgroups based on the median value of risk scores using the same formula as in the TCGA cohort.

Validation of Gene Expression

Next, we performed quantitative real-time polymerase chain reaction (qRT-PCR) analysis to quantify expression of DEGs used for model construction in the validation cohort. Briefly, total RNA was extracted from thoroughly ground (under liquid nitrogen) target tissues using the Trizol reagent (Life Technology, Grand Island, NY, United States). The RNA was reverse transcribed to complementary DNA (cDNA) using the RevertAid First Strand cDNA Synthesis Kit (Thermo Scientific, Lithuania), then subjected to qRT-PCR using the SYBR® Green qPCR mix 2.0 kit performed on the Roche LightCycler® 480 Real-Time PCR System. The primers for the genes targeted in this study were obtained from TsingKe Biological Technology (Nanjing, China), and their sequences are as follows: PYGL (forward 5′- CAG​CCT​ATG​GAT​ACG​GCA​TTC -3′, reverse 5′- CGG​TGT​TGG​TGT​GTT​CTA​CTT​T-3′), ALOX15 (forward 5′-GGG​CAA​GGA​GAC​AGA​ACT​CAA-3′, reverse 5′- CAG​CGG​TAA​CAA​GGG​AAC​CT-3′), BCL2 (forward 5′- GGT​GGG​GTC​ATG​TGT​GTG​G -3′, reverse 5′-CGG​TTC​AGG​TAC​TCA​GTC​ATC​C-3′), IFNA1 (forward 5′-GCC​TCG​CCC​TTT​GCT​TTA​CT-3′, reverse 5′-CTG​TGG​GTC​TCA​GGG​AGA​TCA -3′), TLR3 (forward 5′- TTG​CCT​TGT​ATC​TAC​TTT​TGG​GG -3′, reverse 5′-TCA​ACA​CTG​TTA​TGT​TTG​TGG​GT -3′) β-actin (Forward: 5′-CGA​GCA​CAG​AGC​CTC​GCC​TTT​GCC-3′, Reverse: 5′-TGT​CGA​CGA​CGA​GCG​CGG​CGA​TAT-3’). Relative mRNA expression was calculated using the 2-ΔΔCt method.

Functional Enrichment and Drug Sensitivity Analyses

We first used the “cluster Profiler” package in R to perform Gene Ontology (GO) functional and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses on the screened necroptosis-related DEGs (Zeng et al., 2021). Next, we applied the “org.Hs.eg.db” package to identify significantly enriched genes and classify gene clusters. All analyses were performed at a statistical significance of p < 0.05. The obtained DEG list was compared with a reference dataset from the CMap database, to obtain a correlation score based on enrichment of DEGs in the reference gene expression profile. A positive number indicated that the DEG was similar to the reference gene expression profile, thus allowed analysis of similar interventions that produce the above effects. Conversely, a negative number indicated that the DEGs had an opposite gene expression profile to that of the reference, suggesting that the drug may antagonize the DEGs’ expression. Ultimately, potentially effective drugs that can counteract drug resistance were inferred from the genetic changes in drug-resistant cell lines, and ranked according to correlation scores obtained on the reference gene expression profile.

Analysis of Immune Infiltration and Tumor Microenvironment

Considering that immune infiltration levels are correlated with survival and prognosis of cancer patients, we evaluated the relationship between risk-scores and immune infiltration levels. Specifically, we applied multiple algorithms implemented in the IOBR package in R (Zeng et al., 2021) to estimate proportions of tumor-infiltrating immune cells. Next, we calculated stromal and immune scores for each sample using the ESTIMATE algorithm implemented in the ‘limma’ and ‘estimate’ packages in R.

Results

Profiles of DEGs

A total of 499 and 80 PRAD patients from the TCGA and validation cohorts, respectively, were included in this study. Detailed clinical characteristics of these patients are shown in Tables 1, 2. A total of 16 necroptosis-related genes were differentially expressed between tumor and adjacent non-tumor tissues. Profiles of these expression patterns are presented using heatmaps in Figure 1A, while the relationships among the DEGs are shown in Figure 1B. Next, we performed LASSO regression analysis to screen covariates, then applied a 10-fold cross validation with minimum criteria to obtain an optimal λ. The final value of λ yielded a minimum cross validation error. Consequently, a total of 5 DEGs, namely ALOX15, BCL2, IFNA1, PYGL and TLR3, were identified (Figures 2A,B). Incorporation of these significant DEGs into a multivariate Cox regression model revealed coefficients of included each gene. Finally, each PRAD patient was assigned a separate risk score according to the aforementioned formula.

www.frontiersin.org

TABLE 1. Basic clinical characteristics of PRAD patients in TCGA PRAD cohort.

www.frontiersin.org

TABLE 2. Basic clinical characteristics of PRAD patients in validation cohort.

www.frontiersin.org

FIGURE 1. Identification of the candidate genes. (A). Heatmap showing differentially expressed genes between the two groups. (B). The relationship among necroptosis-related DEGs.

www.frontiersin.org

FIGURE 2. Processes of LASSO Cox model fitting. (A). Profiles of coefficients in the model at varying levels of penalization were plotted against the log(lambda) sequence. (B). Ten-fold cross-validated error (the first and second vertical lines denote the minimum and cross-validated errors, respectively, within 1 standard error of the minimum).

Construction and Validation of a Prognostic Model

Risk scores = 0.05834 expression level of ALOX15–0.0004914 expression level of BCL2+ 0.665* expression level of IFNA1+0.03093 expression level of PYGL −0.0008311 expression level of TLR3. The resulting median cut-off value was used to stratify patients into high-risk and low-risk groups, comprising 249 and 250 individuals, respectively (Figure 3A). A principal component analysis (PCA) plot showed that patients in both groups were distributed in different directions, indicating that the established model had excellent predictive ability to distinguish between high and low-risk PRAD (Figure 3B). Kaplan-Meier curves showed that patients in the low-risk group were significantly associated with higher survival rates compared to their high-risk counterparts (p < 0.05) (Figure 3C). In addition, we generated time-dependent ROC curves to estimate performance of the risk prediction model. AUC values for the prognostic model were 0.822, 0.856, and 0.795 for 1-, 3-, and 5-years survival, respectively (Figure 3D). The calibration curve was close to 45°, indicating that the model had good prognostic performance (Figure 3E).

www.frontiersin.org

FIGURE 3. Prognostic value of the 5-gene signature model in the test (TCGA) cohort. (A). Distribution and median values of the risk scores. (B). PCA plot. (C). Kaplan-Meier curves showing OS of patients in the high-risk and low-risk groups. (D). AUC values of the time-dependent ROC curves showing prognostic performance of the risk score. (E). Calibration plots for the established model.

qRT-PCR results showed that BLC2, PYGL, and TLR3 were significantly downregulated in PRAD relative to normal adjacent tissues in the validation cohort (Figure 4A). Results from risk score analysis for each patient, calculated using the same formula and critical values, indicated that patients in the high-risk group had significantly worse overall survival rates than those in the low-risk group (p < 0.05) (Figure 4B). AUC values for the 5 necroptosis-related gene signature were 0.836, 0.669, and 0 0.726 at 1, 3 and, 5 years, respectively (Figure 4C). Figure 4D shows the statistically significant results for the genes that make up the model. Multivariate cox regressionreveals that the most significant in 5 DEGs was PYGL. Based on these findings, we further explored the relationship between PYGL and multiple pathways, and found that PYGL expression was positively correlated with IFN−Gamma signature, APM signaling, Proteasome, Basal differentiation, EMT differentiation, Immune differentiation, Myofibroblasts, Interferon response, and Keratinization, but had a negative correlation with Urothelial differentiation, Luminal differentiation and Neuroendocrine differentiation (Figure 4E). Collectively, these results suggested that PYGL expression might be involved in cancer progression.

www.frontiersin.org

FIGURE 4. Validation of the 5-gene signature using a validation cohort. (A). Results of qRT-PCR analysis. (B). Kaplan-Meier curves showing OS of patients in the high-risk and low-risk groups. (C). AUCs of time-dependent ROC curves indicating the prognostic performance of the risk score. (D). Kaplan-Meier curves showing OS of patients in the high- and low-expression groups of 5 DEGs. (E) Pathway analysis targeting PYGL.

Independent Prognostic Value of the Risk Model

Next, we investigated the prognostic significance of different clinical features in PRAD patients. Results from univariate Cox proportional hazard regression analysis revealed that the risk score had prognostic significance in OS (HR = 1.67, 95% CI: 1.23–2.79, Table 3). These results remained significant (HR = 1.42, 95% CI: 1.14–1.75, Table 3) in PRAD patients. Validation of this signature’s prognostic value revealed that it was an independent prognostic factor for PRAD (Table 3). A heatmap of DEGs and distribution of patients’ age as well as survival status across both low- and high-risk subgroups are shown in Supplementary Figure S1.

www.frontiersin.org

TABLE 3. univariate cox and multivariate cox regression analyses.

High-Risk Scores Were Associated With a Hot Tumor Microenvironment

Previous studies have shown that the TME plays a crucial role in tumor extension, progression, migration and invasion (Chalmin et al., 2018). GO terms showed that the DEGs were significantly enriched in immune-related GO entries, including lymphocyte differentiation, and macrophage aggregation. Moreover, KEGG pathway enrichment results revealed that these genes were enriched in cell adhesion molecules, cytokine and cytokine receptor interactions, and necroptosis (Figure 5A). Overall, these results indicated that the DEGs were significantly enriched in immune-related activities, which implies that immune factors represent the main feature of the TME in PRAD. Furthermore, risk scores were positively correlated with the tumor microenvironment (ESTIMATEScore) (Figure 5B). Furthermore, multiple algorithms revealed that PYGL expression was strongly associated with macrophages, eosinophils and activated CD4+ T cells among others, indicating that the established prognostic model was associated with a hot TME (Figure 6A). ssGSEA results revealed significant enrichment of many types of immune cells, such as macrophages, and Tregs, among others, between the groups, while immune-related functions were also significantly different between the groups (Figure 6B).

www.frontiersin.org

FIGURE 5. (A). Profiles of GO functional terms and KEGG pathway enrichment for the identified DEGs (B). An overview of the ESTIMATE algorithm.

www.frontiersin.org

FIGURE 6. (A) The relationship between PYGL expression and immune cell infiltration. (B). Comparison of ssGSEA scores between the two risk groups in the test (TCGA) cohort.

DEG-Based Tumor Classification

Next, we performed K-means cluster analysis to stratify PRAD patients into subtypes according necroptosis-related DEGs, and obtained the highest intra- and low inter-group correlations when k = 2 (Figure 7A). Based on results of the 16 DEGs, the 499 PRAD patients were divided into 2 clusters. We generated a heatmap to illustrate the resulting gene expression profiles and found that they were matched to the clinical data (Figure 7B). A comparison of OS revealed significant differences between the clusters (Figure 7C).

www.frontiersin.org

FIGURE 7. Tumor classification based on the identified pyroptosis-related DEGs (A). Patients were grouped into 2 clusters according to the consensus clustering matrix (k = 2). (B). Kaplan–Meier curves showing OS of patients in the 2 clusters. (C). A heatmap and clinicopathologic characteristics of patients in the 2 clusters based on the identified DEGs.

Drug Sensitivity Analysis

The highest negative correlation score was obtained in Thioridazine (-0.703), a drug used for treatment of acute schizophrenia, mania and depression. This suggests that this drug has potential therapeutic effects in PRAD. On the other hand, the highest scores were obtained in trichostatin A (an antitumor inhibitor), LY-294002 (the first synthetic protein kinase inhibitor), Sirolimus (an immunosuppressant), Tanespimycin (an antitumor agent), Monorden (an antibiotic), Sirolimus (an immunosuppressant), Tanespimycin (an antitumor agent), and Monorden (an antibiotic). Estradiol (estradiol) is a transdermal estrogenic therapeutic agent that has been used to treat advanced PRAD. The above-mentioned drugs exhibited a strong negative correlation suggesting that they may have potential form treatment of PRAD (Figure 8).

www.frontiersin.org

FIGURE 8. Drug sensitivity analysis.

Discussion

Tumor growth occurs due to an imbalance between tumor cell death and growth (Meier et al., 2000). Notably, inhibition of excessive cell proliferation or normal cell death in the body markedly exacerbates incidence of malignant tumors. Therefore, some researchers believe that unlimited cell proliferation and death inhibition represent the two distinguishing features of malignant tumors. Necroptosis is a newly discovered form of cell death that morphologically manifests with similar characteristics to those observed in necrosis. However, the 2 phenomena differ in that while necrosis is a passive death caused by external physicochemical stresses, such as infection or inflammation, and is not regulated by signaling pathways, necroptosis is a regulated by programmed death. Numerous studies have demonstrated that necroptosis plays an important role in cancer initiation, progression, and metastasis (Jouan-Lanhouet et al., 2014; Barbosa et al., 2018). In the present study, we evaluated the relationship between necroptosis-related genes and prognosis of PRAD patients. Summarily, we screened for key genes that can independently predict prognosis of PRAD patients, then used them to construct a prognostic prediction model. Next, we verified the predictive power of the established model, then applied a multifactorial Cox regression model to assess the effect of other clinicopathological parameters on the signature’s prognostic value in PRAD patients.

A total of 5 DEGs were included in the current model, with ROC and calibration curves in both training and validation cohorts showing that the model had excellent power in predicting PRAD patients. In addition, we selected PYGL, a gene located on chromosome 14q22.1 with a total of 20 exons that has been widely used as a building block for predictive models (Luo et al., 2020), and investigated its role in PRAD. Results showed that high PYGL expression was an independent predictor of poor prognosis in PRAD patients, consistent with a pervious study that reported similar findings in glioma patients (Luo et al., 2020). Our results further demonstrated that high PYGL expression was closely associated with infiltration of immune cells in tumors, with this expression pattern also positively correlated with cancer associated fibroblasts (CAFs). Previous studies have shown that CAFs inhibit T cell infiltration by secreting peritumoral collagen and TGF-Beta/PD-L1 specific antibody YM101, while M7824 could effectively suppress CAFs activity thereby promoting T cell infiltration (Lan et al., 2018; Yi et al., 2021a; Yi et al., 2021b). Additional research evidences have shown that the extracellular matrix produced by CAFs can also limit efficacy of tumor therapy (Ziani et al., 2018; Dou et al., 2020). Therefore, CAFs have become a new target for tumor therapy. For example, inhibition of tumor progression by targeting and regulating CAFs or overcoming their barrier effect has become a new approach for improving efficacy of tumor therapy (Ziani et al., 2018). Since only a handful of reports have described the specific function of PYGL in tumorigenesis of PRAD, results of the present study provide new insights into the relationship between necroptosis and PRAD based on PYGL functions.

Previous studies have demonstrated the role of the tumor microenvironment, especially the immune microenvironment, in tumor biology (Sumida et al., 2011; Xu et al., 2020; Li et al., 2021). Results of the present study showed that the identified DEGs were mainly enriched in immune-related pathways. Consequently, we performed an immune cell infiltration analysis to assess the relationship between risk scores and overall survival of patients in the context of immune system response. Our results indicated that low-risk patients recruited more immune cells and triggered higher activation of immune pathways than their high-risk counterparts. Notably, the former group exhibited marked enrichment of NK cells. Previous studies have showed that NK cells have a strong cytocidal activity and do not require MHC activation to kill tumor cells (Shu and Cheng, 2020). In addition, tumors have been found to induce production of type II NKT cells, which in turn secrets IL-13, thereby not only causing aggregation of MDSC in the tumor microenvironment but also activating the STAT6 signaling pathway and suppressing CD8+ T cell function (Terabe et al., 2000; Terabe et al., 2005). Consequently, the use of NK cells for tumor immunotherapy has gained numerous attention from all over the world. Some of the current NK cell-based clinical approaches for tumor therapy include cytokine-mediated NK cell activation, autologous or allogeneic NK cell transfer, CAR gene-modified NK cells and memory cells. In addition to hematological tumors, modified NK and memory-like NK cells have shown great potential for treatment of liver, non-small cell lung, colorectal, ovarian and breast cancers (Li et al., 2015). Results from our drug sensitivity analysis revealed several potential drugs that might modulate the necroptosis-related genes. Notably, expression of these necroptosis-related genes was negatively correlated with thioridazine, trichostatin A, LY-294002, sirolimus, tanespimycin and monorden, suggesting that these could be potential new options as therapeutic drugs.

Conclusion

In summary, we identified differentially expressed necroptosis-related genes, between PRAD and normal adjacent tissues, and used them to establish a model for predicting prognosis of PRAD patients. Moreover, we revealed the correlation between risk scores and immune activities. However, further studies are needed to elucidate the mechanisms underlying necroptosis in tumor immunology.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Ethics Statement

The studies involving human participants were reviewed and approved by the All patients signed a written informed consent prior to inclusion in the study, and ethical approval for the study protocol was obtained from The First Affiliated Hospital of Zhengzhou University. The patients/participants provided their written informed consent to participate in this study.

Author Contributions

X-TY, J-XY, L-XS designed experiments; L-YZ, X-YL, carried out experiments, and wrote the manuscript, X-TY, L-XS performed manuscript review. All authors read and approved the final manuscript.

Funding

This study received Fundamental research program funding of Ninth People’s Hospital affiliated to Shanghai Jiao Tong university School of Medicine (No. JYZZ076), Clinical Research Program of Ninth People’s Hospital, Shanghai Jiao Tong University School of Medicine (No. JYLJ201801, JYLJ201911), the China Postdoctoral Science Foundation (No. 2017M611585) and the National Natural Science Foundation of China (No. 81871458).

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.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

Home for Researchers editorial team (www.home-for-researchers.com).

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fbioe.2021.814813/full#supplementary-material

Abbreviations

AUC, Area under the curve; APC, Antigen presenting cell; CCR, Cytokine-cytokine receptor; CI, Confidence interval; DEGs, Differentially expressed genes; FDR, False discovery rate; GO, Gene Ontology; HLA, Human leukocyte antigen; HR, Hazard ratio; KEGG, Kyoto Encyclopedia of Genes and Genomes; OS, Overall survival; PRAD, Prostate adenocarcinoma; ROC, Receiver operating characteristic; ssGSEA, Single-sample gene set enrichment analysis; PCA, Principal component analysis; aDC, Activated dendritic cell; iDC, Immature dendritic cell; pDC, Plasmacytoid dendritic cell; TCGA, The Cancer Genome Atlas; Tcm, T central memory; Tfh, T follicular helper cell; TIL, Tumor Infiltrating Lymphocyte; Tgd, T gamma delta.

References


Disclaimer:

This article is autogenerated using RSS feeds and has not been created or edited by OA JF.

Click here for Source link (https://www.frontiersin.org/)