AMCO Journal

Table 1: List of differentially expressed necroptosis-related genes; CPM counts per million, FC fold change, FDR false discovery rate, HDAC Histone deacetylase 9, PLA2G4E Phospholipase A2 group IVE.

Figure 1: Network between differentially expressed NRG and NRG-related lncRNA; lncRNA long non-coding RNA, NRG necroptosisrelated gene. Lines between two nodes indicate significant correlations between lncRNA and NRG.

Construction of a prognostic prediction model for risk assessment

The complete process of model construction was conducted in the training set (N=432). Twelve NRG-related lncRNAs were significantly associated with prognosis according to univariate Cox regression analysis (Figure 2A). Seven of them were indicative of poor prognosis while the other five suggested a protective element of prognosis. Ten lncRNAs with significantly association with prognosis were acquired by LASSO analysis (Figure 2B-2C). Ultimate six prognosis-related lncRNAs were resulted from multivariate Cox analysis, in which expression of three lncRNAs [Z68871.1, SPRY4-AS1, AC011247.1, hazard ratio (HR)>1] were high-risk elements for prognosis and the other three played a protective role in prognosis (COL4A2-AS1, CCDC144NL-AS1, ADAMTS9-AS2, HR<1) (Figure 2D). Based on the results of multi-Cox regression, we obtained a risk score of each sample in the training set using the formula below:

Figure 2: Construction of a risk model with prognostic necroptosis-related lncRNAs. (a) Forest plot of univariate Cox regression analysis of prognosis-related lncRNAs; (b) 10-fold cross-validation for selection of variables in the LASSO model; (c) The LASSO coefficient profile of ten survival-associated lncRNAs; (d) Forest plot of multivariate Cox regression analysis of six prognosis-related lncRNAs; LASSO least absolute shrinkage and selection operator, lncRNA long non-coding RNA.

Risk score= expr (COL4A2-AS1) * (-2.872) + expr (CCDC144NL-AS1) * (-0.887) + expr (ADAMTS9-AS2) * (-2.694) + expr(AC011247.1) * (1.232) + expr (SPRY4-AS1) * (0.767) + expr (Z68871.1) * (0.763). (2)

Then we used this formula to calculate the risk score of samples in the testing set. In both groups, samples were stratified into highrisk and low-risk levels according to the median risk score (Supplementary Figure 2A-C). Survival analysis suggested that patients in high-risk group had a significant poorer overall survival than those of low-risk group, no matter in the training or testing set (Figure 3A3C, Supplementary Figure 2D-2F). Heatmap reaffirmed the results that lncRNAs Z68871.1, SPRY4-AS1, and AC011247.1 were highly expressed in high-risk groups while the other three lncRNAs COL4A2-AS1, CCDC144NL-AS1, and ADAMTS9-AS2 were highly expressed in low-risk groups, no matter in training, testing or the complete datasets, suggesting their prognostic meaning (Supplementary Figure 2G-2I). The association of clinicopathological factors other than risk scores with overall survival was analyzed and the results showed that risk score together with age, tumor stage, subtype, and type of surgery were independent risk factors of overall survival. Multi-Cox regression analysis demonstrated that older age, higher staging of BC, TNBC, and more invasive operation of BC were indicative of a poorer prognosis (Figure 3D-3E).

Figure 3: Prognostic analysis of risk model in different sample sets. (a-c) Kaplan–Meier survival curves of overall survival of breast cancer patients between high- and low-risk group in training, testing, and entire samples sets, respectively; (d-e) Forest plots visualizing independent clinicopathological variables of overall survival in the entire sample set.

ROC curves illustrating the area under curve (AUC) in terms of 1-, 3-, 5- , and 10-year were 0.718, 0.623, 0.670, and 0.720 (Figure 4A). Regarding 10-year ROC curves of the model, the AUC of the variable “risk score” was larger than any other clinicopathological variables (AUC=0.720), illustrating a relatively powerful predictive capability (Figure 4B). The HR of risk score in uni- and multivariate Cox model was 1.208 and 1.236, respectively (both p value<0.0001).

Figure 4: ROC curves of the risk model in different comparisons. (a) 1-, 3-, 5-, 10-year ROC curves of the model in complete sample set; (b) 10-year ROC curves for different clinicopathological variables; AUC area under curve, ROC receiver operating characteristic.

Nomogram construction and calibration plot assessment

Based on the results of the multivariate Cox model for survival analysis, we brought variables “risk score”, “age”, “stage”, “subtype”, and “type of surgery” into the nomogram, a visualization tool for prediction of 1-year, 3-year, 5-year, and 10-year survival rate (Figure 5A). Calibration plot demonstrates that the nomogram exhibited an excellent concordance with prediction (Figure 5B).

Figure 5: Nomogram of the model and the calibration curves for model validation. (a) Nomogram of the model for prediction of 1-, 3-, 5-, and 10-year overall survival; (b) calibration curves validating the prediction of 1-, 3-, 5-, and 10-year overall survival; Lump lumpectomy, Mast simple mastectomy, MRM modified radical mastectomy, NTNBC non-triple-negative breast cancer, OS overall survival, TNBC triple-negative breast cancer.

Biological pathway analysis

Through GSEA analysis, we summarized different KEGG biological process and pathways in high-risk (N=417) and low-risk (N=447) groups. Each five of significantly enriched biological functions and/or pathways were collected from high- and low-risk groups, respectively. Amongst, high-risk groups included DNA mismatch repair, DNA replication, nucleotide excision repair, RNA degradation, and cell cycle regulation; low-risk groups included cytokine receptor interaction, calcium signaling pathway, TGF-β signaling pathway, regulation of actin cytoskeleton, and cancer pathway (Figure 6).

Figure 6: GSEA analysis of the top five KEGG biological pathways in high- and low-risk groups, respectively. GSEA Gene Set Enrichment Analysis, KEGG, Kyoto Encyclopedia of Genes and Genomes.

Comparison of immune signature, TME between high- and low-risk groups

Through ssGSEA, significant discrepancy of immune cell abundance could be found between high- and low-risk groups, especially cytotoxic T cells, B cells, natural killer cells (NK cells), dendritic cells (DCs), and neutrophils. Low-risk group had a higher immune cell infiltration than its counterpart (Figure 7A), which is in line with the result of the activity of major immune functions (cytolytic activity, T-cell costimulation, and interferon response, etc.) between the two groups (Figure 7B).

Figure 7: ssGSEA analysis of immune landscapes in high- and low-risk group in the entire sample set. (a) Enrichment level of major immune cells in samples from high- and low-risk groups; (b) Activity level of typical immune functions or pathways between high- and low-risk groups; aDC activated dendritic cell, APC antigen presenting cell, CCR chemokine receptor, DC dendritic cell, HLA human lymphocyte antigen, iDC immature dendritic cell, IFN interferon, MHC major histocompatibility complex, NK cell natural killer cell, pDC plasmacytoid dendritic cell, ssGSEA single sample Gene Set Enrichment Analysis, Tfh follicular helper T cell, Th helper T cell, TIL tumor-infiltrating lymphocyte, Treg regulatory T cell; ***, p<0.001; **, p<0.01; *, p<0.05.

Correlation analysis between immune cell infiltration and risk score was also performed and data were derived from several algorithms (XCELL, TIMER, QUANTISEQ, MCPCOUNTER, EPIC, and CIBERSORT, Figure 8A). Take TIMER algorithm for instance, scatter plots demonstrate that T cells (CD4⁺, CD8⁺), B cells, DCs, and neutrophils were all negatively correlated with the risk score which corroborated the results of ssGSEA above (Figure 8B-8F).

Figure 8: Correlation of major innate and adaptive immune cell infiltration with different risk scores in the entire sample set. (a) Bubble plot demonstrating the correlation between immune cell infiltration with risk score via different various algorithms; (b-f) correlation of infiltration scores of CD4⁺ T cells (b), CD8⁺ T cells (c), B cells (d), myeloid Dendritic cells (e), and neutrophils (f) with risk scores.

Furthermore, we calculated immune scores and stromal scores for each BC patients and the boxplots show that patients with higher risk had lower immune/stromal/ESTIMATE scores than patients with lower risk (Figure 9A-9C).

Figure 9: Comparison of tumor microenvironment of samples in entire sample sets between groups with different risk levels. (a) Comparison of immune cell infiltration between high- and low-risk group; (b) Comparison of stromal cell infiltration between high- and low-risk group; (c) Comparison of ESTIMATE score between high- and low-risk group. ESTIMATE score is the sum of the immune and stromal score.

Comparison of immune checkpoint molecules and potential immunotherapy response between different risk levels

47 immune checkpoint molecules were collected and the expression of these genes were extracted from the original FPKM file of female BC samples. We compared the expression level of these genes between high- and low-risk group, which is demonstrated in Figure 10.

Figure 10: Expression of typical types of immune checkpoint molecules between high- and low-risk groups in the entire sample set; ***, p <0.001; **, p<0.01; *, p<0.05.

Results showed that expression of genes in low-risk group were higher than that of high-risk group except CD80. IPS of patients in high/low-risk groups were obtained. Score 2 (IPS-PD-1/PD-L1/PD-L2 blocker), score 3 (IPS-CTLA-4 blocker), and score 4 (IPS- PD-1/PD-L1/PD-L2 and CTLA-4 blocker) in low-risk group were significantly higher than those in high-risk group (Figure 11A-11D), suggesting a better expected response of ICIs (anti-PD-1/anti-CTLA-4 monoclonal antibodies) in the former group. The results were consistent with the counterpart of immune landscape and TME mentioned in the previous section.

Figure 11: Comparison of IPS of anti-PD-1 and anti-CTLA-4 monoclonal antibody between high- and low-risk group in the entire sample set. (a) pure IPS-1 score between two groups; (b) IPS-1/PD-1 blocker score between two groups; (c) IPS-1/CTLA-4 blocker score between two groups; (d) IPS-1/PD-1/CTLA-4 score between two groups; CTLA-4 cytotoxic T lymphocyte-associated antigen-4, IPS immunophenoscore, neg negative, pos positive, PD-1 programmed death-1.

Discussion

In our study, seven DEGs were found related to necroptosis between NTNBC group and TNBC group. Amongst, PLA2G4E and HDAC9 were upregulated in TNBC samples while H2AC 4/12/14/20 were downregulated in TNBC samples. PLA2G4E belongs to the cytosolic phospholipase A2 group IV family. This group of genes regulate membrane tubule-mediated transport. Proteins produced by this gene get involved in recreating tubules to move specific cargo molecules back to cell surface (see details in https://www.ncbi.nlm.nih.gov/gene/123745). Previous studies have found that PLA2G4E was related to necroptosis. Inhibition of PLA2G4E by microRNAMir504-5p reduced necroptosis of skin and prolonged survival of skin cells [50]. Previous studies have found that inactivated function of RIPK1, RIPK3, and MLKL, especially by the necroptosis inhibitor Nec‐1, protected neurons from being damaged, degenerated, and keeping the viability of neural cells in these neurodegenerative diseases [51-55]. In an animal model of Alzheimer ́s disease, PLA2G4E served as a protective role for the survival of cerebral tissue. Overexpression of PLA2G4E could resume memory and cognition, indicating a potential role in treatment of neural degenerative disease [56]. It is worthwhile to further explore relationship between this gene and necroptosis in neurodegenerative diseases. HDAC9 is a member of the class IIa family of histone deacetylases. This gene family get involved

in gene silencing by modifying chromatin into a transcriptionally repressed conformation. HDAC9 are expressed in many human tissues and organs, especially in brain and skeletal muscles [57,58]. HDAC9 negatively regulates muscle differentiation by interacting with myocyte enhancer factor (MEF2) and inhibiting its transcriptional activity [59]. Similar to PLA2G4E, HDAC9 also plays a protective role in the nervous system. One study has found that reduced expression of HDAC9 caused death of neural cells [60]. HDAC9 possesses several functions in carcinogenesis as well. It has been found to be upregulated in various types of tumors including BC, hepatocellular carcinoma, pancreatic cancer, osteosarcoma, and hematological malignancies [4, 61-65]. In BC, researchers found that overexpression of HDAC9 was associated with metastasis and poor prognosis [62]. HDAC9 upregulated NF-κB expression and its downstream inflammatory pathways, stabilizing TME and leading to recurrence of BC [61]. Another study led by Salgado et.al found that HDAC9 was overexpressed in TNBC, which was consistent with our results [66]. HDAC contributed to increased invasion and angiogenesis in TNBC tissues via VEGF and MAPK3 signaling pathway, which is negatively affected by miR-206. A bench study from China using genome-wide CRISPR screening found that HDAC9 was a pivotal molecule in regulating paclitaxel-resistance of TNBC. Inhibition of HDAC9 caused tubulin instability and cell cycle arrest, leading to cancer cell death and increasing the sensitivity to chemotherapy [67]. Other related studies also found that HDAC9 suppressed the expression of Erα and endowed antiestrogen-resistance to BC cells [68]. These results indicated a potential HDAC-targeting therapy in the treatment of TNBC to suppress its aggression and improve the sensitivity of TNBC to systemic therapy. Different from the upregulated genes, the significantly downregulated NRGs in TNBC samples were all H2A clustered genes. Proteins encoded by these genes were all belong to replication-dependent histone H2A family (see https://www.ncbi.nlm.nih.gov/gene/8338). A study led by Su found that H2AC was upregulated in ER-positive BC. Suppression of this replication-dependent histone H2A isotype caused abnormal estrogen signaling and cell cycle arrest of BC [69]. In another study, the authors also furtherly found that H2AC was related to telomeres in human cells and interacted with telomeres protein 1 (POT1) as well as telomere repeat factor 2 (TRF2). Depletion of H2AC would dysfunction the telomere, destabilize chromosome structure, and cause DNA damage [69].

According to the correlation analysis, 153 lncRNAs were significantly associated with six of the seven necroptosisrelated DEGs (PLA2G4E did not appear in the network). After univariate Cox regression, LASSO regression, and multivariate Cox regression analysis, final 6 prognosis-related lncRNAs were brought into the risk prediction model, including Z68871.1, SPRY4-AS1, AC011247.1, COL4A2-AS1, CCDC144NL-AS1, and ADAMTS9-AS2. Several bioinformatic studies found that high expression of Z68871.1 and SPRY4-AS1 were correlated with worse OS in solid tumors including BC, HCC, and squamous cell carcinoma [69-73]. High expression of ADAMTS9-AS2 was related to a better treatment response and prognosis in BC [74,75], which was coherent with our results. However, high expression of COL4A2-AS1, CCDC144NL-AS1, and ADAMTS9-AS2 were related to proliferation and worse prognosis of breast, gastric, colorectal cancer, HCC, and lung cancer [76-81]. Given that details about functions and involved pathways are unclear, more intense research should be conducted to explore the role of these lncRNAs in BC.

According to the expression and coefficient of lncRNAs, we constructed the model and measured the risk scores for each sample in the training set and testing set. Samples in both sets were stratified into two subgroups: high-risk group and low-risk group. Survival analysis in both sets were similar that patients in high-risk group had a worse prognosis compared with low-risk counterparts. In addition to risk score, age, tumor stage, subtype, and type of surgery were all independent variables for prognosis. Nomogram including these variables was created to predict 1-, 3-, 5-, and 10-year survival rates for BC patients. Predictive efficacy was testified by ROC curves and calibration plots. Results showed that risk score was the best independent predictor of overall survival.

Based on the results of GSEA, we found that significantly enriched biological pathways were focused on cell cycle, cellular metabolism, and cellular proliferation as well as other special pathways including TGF-β and Ca2⁺ signaling pathways. These are important pathways that necroptosis gets involved in. For instance, tumor-induced endothelial necroptosis triggered tumor cell extravasation and metastasis. TGF-β-activated kinase 1 (TAK1) could inhibit necroptosis of endothelial cells [82,83]. Studies also found that dysfunction or deficiency of TAK1 enhanced the necroptosis in which RIPK3 was upregulated [84]. In another study, scientists found that thymoquinone, an apoptosis-inducing agent with anticancer activity, could also induce necroptosis through elevating cytosolic calcium concentration by endoplasmic reticulum calcium depletion and activation of store-operated calcium entry (SOCE) in lymphoma cell [85].

We compared the immune signatures of BC samples between high-risk and low-risk group. Results showed that low-risk group had higher level of immune cell infiltration, higher immune, and stromal scores, and higher expression level of immune checkpoint molecules. However, expression of CD80 was lower in low-risk group than that in high-risk group. Also named as B7-1, CD80 could not only bind to immune-stimulating signaling molecule CD28 but also combine with immunosuppressive molecule CTLA4 [86]. As reported, the affinity of CD80 to CTLA-4 was higher than that to CD28, providing its immune-inhibiting feature [87]. A study led by Arutha Kulasinghe found that patients with refractory TNBC had a higher level of CD80 [88], corroborating with our results. More investigation should be performed to validate this conclusion. Correlation analysis illustrates that the abundance of immune cells (DCs, innate, and adaptive immune cells) was all inversely proportional to the risk score. The above results were consistent with the immunotherapy response prediction that lowrisk group had a potentially better response to PD-1 and CTLA-4 monoclonal antibodies. Necroptosis not only directly eliminates tumor cells by promoting cell death, but also indirectly suppress tumor growth and proliferation by inducing immune response in TME [89]. During necroptosis, rupture of the plasma membrane led to outpouring damage associated molecular patterns (DAMPs) into intercellular space, then necroptotic cells recruited immune and inflammatory cells to clear the cell debris and helped tissue repair [89,90]. Animal studies showed that CD8⁺ T cells and NKT cells played major roles in necroptosis-inducing immune response [91,92]. Necroptotic cells offered tumor antigens and inflammatory stimuli to antigen-presenting cells (APC), which then activated CD8 ⁺ T cells by antigen cross-priming [93]. CD8⁺ T cells exerted antitumor cytotoxic immunity by producing and secreting multiple cytokines. The above biological process was highly dependent on RIP1 and NF-κB signaling [91]. Different from CD8⁺ T cells, NKT-induced immune responses were RIPK3dependent, of which downstream signaling mediated by pGAM5 and Drp1 [92]. Counterintuitively, some researchers suggested that inflammation derived from necroptosis was a pro-cancer factor in TME [94]. In the process of carcinogenesis, inflammation could promote tumor cell growth, proliferation, aggressiveness, invasion, and metastasis by providing bioactive molecules [90,95,96]. ROS and other mutagenic chemicals also facilitated evolution of tumor cells by genetic alterations [90]. Given that the above conclusions were controversial, the contribution of immune and inflammatory responses induced by necroptosis in tumor prevention and development needs to be further elucidated.

Despite a comprehensive analysis of necroptosis in TNBC and its relationship with prognosis and immune landscape, some limitations still existed in our study. Firstly, the results were not further validated in basic studies (in vitro and in vivo). The biological functions of prognosis-related lncRNAs were needed to be clarified. Secondly, different extents of biases during data processing and case inclusion/exclusion were inevitable because of a retrospective feature of this study. Thirdly, the expression of PD-1 (CD279) and PD-L1 (CD274) were not included in the original transcriptome data file of TCGA-BRCA. Hence, we failed to directly compare the expression level of these two biomarkers between high- and low-risk groups though data of anti- PD-1 and PD-L1 response could be retrieved from TCIA database.

To conclude, seven differentially expressed NRGs were selected between NTNBC and TNBC groups. Six NRG-related lncRNAs were significantly associated with prognosis and were brought into risk model construction. Risk score is an excellent predictor for overall survival of BC and the model was of high predictive quality. Compared with low-risk group, high-risk group has a gloomier prognosis, less abundance immune cell infiltration, and less responsive ICI treatment. The results pave way for further exploration of precision therapy for TNBC.

Acknowledgements

The authors of this manuscript sincerely appreciate the efforts of all researchers who have contributed the data to the public databases of TCGA, cbioportal, TCIA, and TIMER; GSEA software and algorithm of ESTIMATE. The interpretation and reporting of these data are the sole responsibility of the authors.

Statements & Declarations

Funding

This work was supported by CAMS Innovation Fund for Medical Sciences (CIFMS) under Grant [2021-I2M-1-061 and 2021-1-I2M-003 and 2018-I2M-3-001]; CAMS Clinical and Translational Medicine Research Funds under grant [2019XK320006].

Competing Interests

The authors have no relevant financial or non-financial interests to disclose.

Data Availability Statement

The authors confirm that the data supporting the findings of this study are available within the article and its supplementary materials.

Author Contributions

All authors contributed to the study conception and design. Material preparation, data collection, and analysis were performed by Zijun Zhao, Xiayao Diao, and Qing Cao. The first draft of the manuscript was written by Zijun Zhao and Xiayao Diao and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.