Complicated associations of the 13 TNFRSF genes with the type I anti-tumor response, oncogenes, and ICB-related genes
Using the TCGA/TARGET/GTEx RNA-seq datasets that include 18,299 samples across 23 cancer types and corresponding normal tissues, we investigated the expression of the 13 TNFRSF genes and their associations with markers of major histocompatibility complex class II (HLA-DRA), dendritic cells (BATF3), B cells (PAX5), macrophages (CD68), granulocytes (mainly neutrophils, CD177), T-cells (CD4 and CD8A), type-I anti-tumor responses (IFNG and GZMB), cell proliferation (MKI67), cell death (BAX), cell survival (BCL2), type-II pro-tumor responses (IL5), several known oncogenes (ERBB2, HRAS, KRAS, NRAS, MYC, SRC, TERT, and EGFR), and ICB-related genes (PDCD1 (also known as PD-1), CD274 (also known as PD-L1), CTLA4, and CD86) in a pan-cancer analysis (Fig. 1a).
Using a threshold of false discovery rate (FDR) < 0.05 and |r| ≥ 0.3, the expression of 11 TNFRSF genes was positively correlated with at least one immune cell maker (BATF3, PAX5, CD68, CD177, CD4, and CD8A) in both normal (n = 8,163) and tumor tissues (n = 10,136) (Fig. 1b). Notably, their correlations differed between normal tissues and tumors. For example, in normal tissues, LTBR showed a positive correlation with the macrophage marker (CD68), but not with the neutrophil marker (CD177), whereas the opposite was true in tumors. Another example is TNFRSF11A, which was found to be positively correlated with the neutrophil maker (CD117) in tumors but not in normal tissues. These differences in correlation status may be partly explained by the significant role of these TNFRSF genes in immune cell turnover and immune homeostasis in normal tissues, as evidenced by the observation in normal tissues that cell proliferation (MKI67) and proapoptotic makers (BAX) were tightly linked and both were positively correlated with nine TNFRSF genes. Furthermore, eight TNFRSF genes were positively correlated with IFNG, GZMB, HLA-DRA, CD4, or CD8A, implying that their associations with type-I T-cell and APC responses in the tumor microenvironment. In tumors, CD8A, IFNG, and GZMB were determined to be positively correlated with the APC markers HLA-DRA (MHC class II molecule), BATF3 (dendritic cells), PAX5 (B cells), and CD68 (macrophages). These findings indicate that these TNFRSF genes are linked to cytotoxic T-cell responses (CD8A, IFNG, and GZMB), although they may be provided by APCs in a context-dependent way. In both normal and malignant tissues, there was no significant correlation between the 13 TNFRSF genes and type-II responses (IL5).
In addition, LTBR, TNFRSF9, TNFRSF12A, TNFRSF13C, and RELT were substantially correlated with at least one oncogene in tumors, suggesting the possible associations of these genes with oncogenic pathways. We also discovered that 12 TNFRSF genes were positively correlated with at least one ICB-related gene in tumors. Collectively, our results demonstrate that the 13 TNFRSF genes have complex associations with cytotoxic T-cell responses, oncogenic signaling pathways, and ICB therapy.
Potential prognostic value of the 13 TNFRSF genes in cancers
To evaluate whether the 13 TNFRSF genes can be utilized to improve clinical outcomes, we investigated the relationship between the expression of these TNFRSF genes and cancer patient overall survival (OS) using a pan-cancer analysis (n = 9,346). High levels of TNFRSF14, CD27, and TNFRSF13B expression in tumors were found to be associated with improved OS (all FDR < 0.05; Fig. 1c). This could be due to beneficial cytotoxic T-cell responses (CD8A, IFNG, and GZMB). By contrast, high expression patterns of TNFRSF1B, TNFRSF9, TNFRSF11A, TNFRSF8, TNFRSF13C, TNFRSF12A, LTBR, and RELT, which included five TNFRSF genes that were positively correlated with oncogenes in tumors, were associated with poorer OS (all FDR < 0.05; Fig. 1c and d). Interestingly, although TNFRSF9 and RELT were positively linked to cytotoxic T-cell responses, patients nevertheless had a poorer prognosis, which may be underlined by their positive correlation with the oncogene NRAS. This suggests that a subset of TNFRSF genes influencing cancer patient OS is not only dependent on T-cell-mediated cytotoxicity but may also be associated with oncogenic signaling pathways. As such, the development of prognostic biomarkers based on these TNFRSF genes should take this complexity into account.
Cancer-type-specific regulation of the 13 TNFRSF genes in tumorigenesis
Next, expression alterations between tumors and normal tissues in 23 cancer types were investigated separately to establish the role of the 13 TNFRSF genes in the tumorigenesis of each cancer type. These genes showed differential expression patterns between tumors and normal tissues in multiple cancers with an FDR threshold of 0.05 (Fig. 2a). Notably, these genes exhibited cancer-type-specific expression changes in tumors when compared to normal tissues. RELT, for example, exhibited higher expression levels in tumors of the bile duct, bladder, brain, breast, cervix, colon, endometrium, esophagus, head and neck, kidney, ovary, pancreas, rectum, stomach, and uterine compared to corresponding normal tissues, whereas it had lower expression levels in malignancies of the adrenal, blood, liver, lung, prostate, testis, and thyroid. Furthermore, except for pancreatic cancer, the direction of expression changes for these genes was inconsistent when comparing tumors and normal tissues for each cancer type, suggesting that these TNFRSF genes are more likely to play distinct roles in the tumorigenesis of each cancer type.
Given that melanocytes make up only 10% of the normal skin population [6], a comparison of normal skin and melanoma tissues may be biased due to the predominance of keratinocytes in the epidermal basal layer. Therefore, we validated the expression of these TNFRSF genes in nevus samples and compared them to expression in melanoma tissues using the RNA-seq dataset collected from GSE112509. TNFRSF9, TNFRSF17, and RELT showed higher expression levels in NRASwt melanoma (n = 18) than in NRASwt nevi (n = 3) (all FDR < 0.05; Fig. 2b). This matches well with our above results in the TCGA-melanoma dataset, where TNFRSF9 and TNFRSF17 expression levels were higher in melanoma tissues (n = 102) than in normal skin tissues (n = 557).
The NRAS mutation correlates with reduced expression of the 13 TNFRSF genes in melanoma
Through a pan-cancer analysis, our study establishes significant associations between the expression of a subset of the 13 TNFRSF genes and the expression of oncogenes, including NRAS. To better understand the functional relationships between them, the expression of these genes in melanoma tissues was explored under NRAS mutations, which are observed in 10–25% of melanomas [7]. These TNFRSF genes were demonstrated to be commonly downregulated by the NRAS mutation. With an FDR threshold of 0.05, the expression of TNFRSF1B, TNFRSF8, TNFRSF11A, TNFRSF12A, TNFRSF13B, TNFRSF18, and RELT in NRASmut melanoma tissues (n = 19) was significantly lower than in NRASwt melanoma tissues (n = 18) (Fig. 2b). Four additional genes, CD27, TNFRSF9, TNFRSF13C, and TNFRSF17, were identified when the FDR cutoff was set to 0.1.
Because bulk tumor samples contain a wide range of cell types, cellular heterogeneity among samples may be a substantial confounding factor in the analysis. Therefore, we used the single-cell RNA-seq dataset retrieved from GSE81383 to further investigate the expression of the 13 TNFRSF genes under NRAS mutations in melanoma cells. Eight of them (LTBR, TNFRSF9, TNFRSF11A, TNFRSF12A, TNFRSF13C, TNFRSF14, TNFRSF17, and RELT) were available to be tested in this dataset. TNFRSF9, TNFRSF14, and RELT expression were significantly downregulated in NRASmut/BRAFwt melanoma cells (n = 77) compared to NRASwt/BRAFwt melanoma cells (n = 92) with an FDR threshold of 0.05 (Supplementary Fig. 1). In addition, reduced expression levels of the three genes were detected when comparing NRASwt/BRAFmut melanoma cells (n = 54) to NRASwt/BRAFwt melanoma cells. Mutations in the oncogenes BRAF and NRAS have previously been shown to trigger constitutive activation of the mitogen-activated protein kinase (MAPK) pathway, which promotes tumor development and illness progression [8]. It is plausible that at least a subset of TNFRSF genes is linked to RAS-mediated oncogenic signaling pathways.
To test this, pathway enrichment analysis was performed on the 13 TNFRSF genes, and it was discovered that they were significantly enriched in RAS downstream signaling pathways, such as the nuclear factor kappa B (NF-κB) signaling pathway, regulation of the extracellular signal-regulated kinase 1 (ERK1) and ERK2 cascade, and regulation of the stress-activated MAPK cascade (all FDR < 0.05) (Supplementary Table 1). This link between these TNFRSF genes and RAS-mediated oncogenic signaling pathways is also supported by the functional interaction of these TNFRSF receptors and their intracellular binding partners, TRAFs, which are adaptor molecules regulating various RAS downstream signaling pathways such as the NF-κB, MAPK, and the phosphatidylinositol 3‑kinase (PI3K)/protein kinase B (AKT) [9–11]. Furthermore, it has been proposed that CD40 expression and the RAS/RAF/PI3K pathway have a compensatory role in activating RAS-downstream signaling activity [3, 4]. The proposed mechanism may also function for the 13 TNFRSF genes, as it does for CD40. Further biological experiments in cancer cells are required to corroborate the proposed mechanism.
The expression of the 13 TNFRSF genes has the potential to be exploited as biomarkers for immunotherapy
Our above results elucidate the importance of these TNFRSF genes in tumorigenesis and determine the association between them and ICB-related genes in cancers. As a result, we assessed the potential predictive value of these TNFRSF genes in first-line ICB therapy, such as PD-1/PD-ligand 1 (PD-L1) and CTLA-4 blockade. To accomplish this, we used the Tumor Immune Dysfunction and Exclusion (TIDE) web platform (tide.dfci.harvard.edu), which is a computational tool developed to evaluate the potential for tumor immune escape using gene expression profiles from cancer samples [12]. The area under the receiver operating characteristic curve (AUC) for the 13 TNFRSF genes was measured and compared to four existing immunotherapy biomarkers: CD40, tumor mutational burden (TMB), CD274, and CD8 [4, 13]. Using an AUC threshold of 0.7, we discovered that these TNFRSF genes performed well in at least one cancer type. Compared to four existing biomarkers, these TNFRSF genes were outperformed in at least one ICB cohort, highlighting the potential benefit of using them to predict ICB response (Fig. 2c). Notably, these TNFRSF genes were identified to exhibit cancer-type-specific predictive performance. For instance, among the 13 TNFRSF genes and four existing biomarkers, LTBR is the only biomarker that predicts ICB response in kidney cancer with an AUC ≥ 0.7. Furthermore, in contrast to the four existing biomarkers that only perform effectively in post-treatment tumors, TNFRSF8, TNFRSF9, TNFRSF12A, TNFRSF13C, and TNFRSF18 displayed good predictive power in both post-treatment and pre-treatment/treatment-naïve cancers. This suggests the predictive ability of the five TNFRSF genes not only relies on treatment-induced responses but also on immune escape mechanisms in pre-treatment or treatment-naïve biopsies. Moreover, although most of the 13 TNFRSF genes performed well in predicting the response to both anti-CTLA-4 and anti-PD-1 immunotherapy, TNFRSF8, TNFRSF13C, and TNFRSF17 had proper prediction ability only for anti-PD-1 immunotherapy. This suggests that these three TNFRSF genes may be more predictive of PD-1 blockade-related signals since PD-1 and CTLA-4 blockade show distinct impacts on the melanoma-reactive CD8 + T cell response [14]. Overall, these findings reveal that the 13 TNFRSF genes have context-dependent predictive potential in immunotherapy responses. Supplementary Fig. 2 illustrates their context-dependence in a systematic manner.
The next question is how well a custom biomarker that integrates expression patterns of the 13 TNFRSF genes can predict response outcomes. Although the context-dependent predictive potential of these genes in immunotherapy response has been proven above, we employed the average expression level of these genes to simplify our study. The custom biomarker provided an AUC greater than or equal to 0.7 in a greater number of ICB cohorts (n = 8) than individual TNFRSF genes (n ranging from 1 to 7) (Fig. 2c). Furthermore, we discovered that this custom biomarker outperformed individual genes in two cohorts: the Ruppin study in non-small cell lung cancer treated with PD-1 (AUC = 0.87) and the Nathanson study in melanoma treated with CTLA-4 (AUC = 0.82). These findings imply that these TNFRSF genes, at least to some extent, play a complementary role in prediction and that integrating their expression could facilitate the development of a robust immunotherapy predictive biomarker. Moreover, this custom biomarker was associated with prolonging OS in four cohorts of melanoma patients treated with anti-PD-1 or anti-CTLA4 antibodies (all P < 0.05; Supplementary Fig. 3). This is in line with previous observations that anti-PD-1 and/or anti-CTLA-4 therapy combined with the engagement of TNFRSF members improves cancer immunotherapy success rates [15]. Altogether, our findings show that these TNFRSF genes have the potential to be exploited as immunotherapy biomarkers. Based on purpose, choosing a single gene or effective integration of these genes will help boost prediction power. It is noteworthy that our results are based on analyses of publicly available datasets. Further biological and clinical research is needed to confirm the predictive abilities of these TNFRSF genes before they are used in clinical practice.