The Prognostic and Predictive Role of Xeroderma Pigmentosum Gene Expression in Melanoma.

BACKGROUND: Assessment of immune-specific markers is a well-established approach for predicting the response to immune checkpoint inhibitors (ICIs). Promising candidates as ICI predictive biomarkers are the DNA damage response pathway genes. One of those pathways, which are mainly responsible for the repair of DNA damage caused by ultraviolet radiation, is the nucleotide excision repair (NER) pathway. Xeroderma pigmentosum (XP) is a hereditary disease caused by mutations of eight different genes of the NER pathway, or POLH, here together named the nine XP genes. Anecdotal evidence indicated that XP patients with melanoma or other skin tumors responded impressively well to anti-PD-1 ICIs. Hence, we analyzed the expression of the nine XP genes as prognostic and anti-PD-1 ICI predictive biomarkers in melanoma.
METHODS: We assessed mRNA gene expression in the TCGA-SKCM dataset (n = 445) and two pooled clinical melanoma cohorts of anti-PD-1 ICI (n = 75). In TCGA-SKCM, we applied hierarchical clustering on XP genes to reveal clusters, further utilized as XP cluster scores. In addition, out of 18 predefined genes representative of a T cell inflamed tumor microenvironment, the TIS score was calculated. Besides these scores, the XP genes, immune-specific single genes (CD8A, CXCL9, CD274, and CXCL13) and tumor mutational burden (TMB) were cross-correlated. Survival analysis in TCGA-SKCM was conducted for the selected parameters. Lastly, the XP response prediction value was calculated for the two pooled anti-PD-1 cohorts by classification models.
RESULTS: In TCGA-SKCM, expression of the XP genes was divided into two clusters, inversely correlated with immune-specific markers. A higher ERCC3 expression was associated with improved survival, particularly in younger patients. The constructed models utilizing XP genes, and the XP cluster scores outperformed the immune-specific gene-based models in predicting response to anti-PD-1 ICI in the pooled clinical cohorts. However, the best prediction was achieved by combining the immune-specific gene CD274 with three XP genes from both clusters.
CONCLUSION: Our results suggest pre-therapeutic XP gene expression as a potential marker to improve the prediction of anti-PD-1 response in melanoma.

Introduction

Immune checkpoint inhibitors (ICIs) are a standard treatment for advanced melanoma and other immunogenic tumors. For the therapy of melanoma, they include ipilimumab, a monoclonal antibody directed against the cytotoxic T-lymphocyte-associated antigen 4 (CTLA-4) receptor, and nivolumab or pembrolizumab, antibodies targeting programmed cell death-1 (PD-1) receptor (1–3). Despite the impressive and long-lasting clinical activity of ICIs in some patients, many do not respond. Furthermore, severe side effects are frequent, especially in the combined application of ipilimumab and nivolumab (4). These typically include immune-related adverse events of multiple organs and tissues, leading to inflammations such as thyreoiditis, pneumonitis, colitis or hypophysitis (1, 3, 4). Thus, predictive biomarkers of ICI response are urgently needed in order to identify those patients who achieve the greatest ICI benefit (1–3).

For the efficacy of anti-PD-1 ICIs, different predictive biomarkers have been proposed (5, 6). These can be classified as follows: tumor-intrinsic biomarkers (e.g., tumor mutational burden (TMB) or neoantigen load), which are indirect measures of tumor antigenicity generated by somatic tumor mutations, and immune-specific biomarkers (e.g., T cell-inflamed gene expression profiles (GEPs) or programmed-death-ligand 1 (PD-L1) expression), which are indicative of a T cell-inflamed tumor microenvironment (TME) (6, 7).

Particularly, many studies on immune-specific biomarkers have been conducted recently (8). For instance, Ayers et al. (9) analyzed GEPs using RNA from baseline tumor samples of patients treated with pembrolizumab and eventually defined an 18-gene GEP, hereafter referred to as the Tumor Inflammation Signature (TIS). This signature was predictive in 220 patients with nine different cancers and contained IFN-gamma–responsive genes related to antigen presentation, chemokine expression, cytotoxic activity, and adaptive immune resistance (9). In a current large-scale metanalysis of 1,008 ICI treated cases (n = 353 with melanoma), different predictive biomarkers of ICI response were compared with each other (10). In the markers of immune infiltration category, the TIS single genes CXCL9, CD8A, and TIS itself were the predictors with the strongest effect size. The gene CXCL13 was also a highly predictive gene marker in the whole tumor cohort. Intriguingly, in the three melanoma anti-PD-1 cohorts (7, 11, 12) included in this meta-analysis, CD274 (coding for PD-L1) was a further predictive marker. However, looking at each cohort individually, only in the cohort published by Cristescu et al. (7) TIS, CXCL9, and CD274 were significantly positively associated with ICI response. Finally, the authors concluded that 34 predefined biomarkers (among them the markers of immune infiltration) could only explain about 60% of the total proportion of variance in ICI response, indicating that the remaining factors determining ICI response still need to be discovered (10).

Recent studies revealed that mutational processes directly altering the DNA damage response (DDR) could influence response to ICI (13–16). As one mechanism, DDR defects can lead to a higher TMB, which implicates a greater abundance of immunogenic neoantigens; this is impressively illustrated by the strong clinical activity of anti-PD-1 ICI in mismatch-mediated repair (MMR) deficient tumors (17–21). Notably, besides MMR, two other pathways are responsible for the repair of DNA single-strand breaks (SSB): base excision repair (BER) and nucleotide excision repair (NER). In contrast, DNA double-strand breaks (DSB) are repaired by homologous recombination and further by more error-prone nonhomologous end joining and microhomology-mediated end joining (22–25) Another DNA repair pathway is the Fanconi Anemia/BRCA pathway that restores DNA interstrand crosslinks.

In addition to an increased TMB, other more specific mechanisms leading to altered immunogenicity have been attributed to modified DDR pathways and signaling (6, 15, 25). These mechanisms include upregulation of PD-L1 expression by enhanced DDR signaling through SSB or DSB. Expression of PD-L1 is additionally increased by depletion of BRCA2, which is involved in homologous recombination, or by depletion of Ku70/80, a critical factor of nonhomologous end joining, and by BER reduction (26, 27). Importantly, increased PD-L1 expression after DSB and SSB was associated with the activation of STAT1, STAT3 and IRF-1, which are all part of the canonical interferon (IFN)-gamma-pathway (28). Additionally, for loss of interstrand crosslink repair function in breast cancer, an increased IFN-related gene expression, namely, the two critical mediators of CD4+ and CD8+ T-cell chemotaxis, CXCL10 and CCL5, was discovered (29). Those and other cytokines are involved in T-cell inflammation, which is often a prerequisite for anti-PD-1 ICI response (30). Mechanistically, the crosstalk between immune and cancer cells within the TME, leading to PD-L1 upregulation on cancer cells, is the basis for the mode of action of anti-PD-1 ICI (31). These observations support the joint analysis of DDR pathway and immune-specific gene expression in the TME (32).

Although several case reports stated impressive anti-PD-1 ICI responses of patients with NER germline defects, and while some of them have identified a higher TMB, the further immunogenic impact by alterations of this DDR pathway is far less explored (13, 25, 33–37). Biallelic pathogenic variants in one of the seven NER genes coding for the so-called complementation groups, XPA, ERCC3, XPC, ERCC2, DDB2, ERCC4, ERCC5, the NER gene ERCC1, and the gene coding for XP variant, POLH, are the causes of the rare hereditary disease Xeroderma pigmentosum (XP) (38). They lead to an absent or inactivated protein and are hereafter referred to as the nine XP genes. The NER is mainly responsible for the repair of UV-induced DNA lesions and is divided into a global genome (GG) and transcription-coupled (TC) repair subpathway, which shares a common end section (39). XP patients under age 20 years have a 10,000-fold increased risk for non-melanoma skin cancer and a 2,000-fold increased risk for melanoma, making skin cancer the most common cause of death in this population (40). Hence, XP patients with skin tumors could benefit greatly from successful ICI treatment, requesting investigation of the role of these nine XP genes for ICI response.

A recent study correlated DDR pathway mutations irrespective of XP disease with overall survival of 1,661 ICI-treated patients and revealed that the NER pathway was predictive of ICI benefit—independent of TMB and tumor type. However, in 40,181 unique cancers, only 3.4% of melanomas possessed NER gene mutations (41). Moreover, Litchfield et al. found no predictive role of DDR pathway mutations for ICI response in seven different tumor types (10). An aspect, presumably limiting further the predictive role of DDR mutations, is that different genes in the same DDR pathway can unevenly affect the TME and the ICI response, as shown for BRCA1 and BRCA2 mutations (11, 42).

Based on the discussion above, we focused on analyzing the nine XP genes away from mutation data to gene expression data to investigate the predictive role of XP gene expression as an anti-PD-1 response marker in melanoma. Accordingly, The Cancer Genome Atlas Skin Cutaneous Melanoma project (TCGA-SKCM) (43) dataset, consisting of systemic treatment-naïve primary and metastatic melanoma samples, was used to identify two primary clusters of XP gene expression. We discovered that these were inversely correlated with the expression of TIS and single immune infiltration genes. In TCGA-SKCM, no significant negative correlations between XP genes and TMB were observed. Importantly, besides being predictive for the response to a specific treatment, biomarkers can also be prognostic by providing information about the patients overall cancer outcome, regardless of therapy (44). Because this can potentially interfere with their predictive value, we used TCGA-SKCM to analyze the prognostic role of different factors, and from the XP genes found only the expression of ERCC3 to be prognostic. In contrast, expression of XP genes and clusters thereof could better predict response to anti-PD-1 ICs than well-established immune-specific biomarkers in two pooled clinical cohorts.

Material and Methods

Study Design

In , we outline our analysis workflow beginning with the pre-processing of our three input datasets from TCGA-SKCM (43), Hugo et al. (11), and Riaz et al. (12). In parallel, we parsed the genes to be analyzed, the nine XP genes (XPA, ERCC3, XPC, ERCC2, ERCC4, ERCC5, DDB2, POLH, ERCC1), the 18 genes of the T- cell inflamed signature (9) (here named Tumor Inflammation signature, TIS) and the predictive biomarker, CXCL13, were retrieved through literature research (6, 10). The further utilized TIS score was calculated as the weighted sum of the 18 gene expression values according to Ayers et al. and Cristescu et al. (7, 45).

Figure 1

A schematic diagram of our workflow, including TCGA-SKCM and two anti-PD-1 cohorts of melanoma patients.

Accordingly, the XP gene expression in TCGA-SKCM was clustered hierarchically to define two XP clusters consisting of the mean expression of the corresponding genes. As an additional parameter, we included the TMB for the TCGA-SKCM data in our workflow. An underlying association was investigated via Spearman correlation between the computed scores, the particular gene expressions of XP genes, and predictive biomarkers and TMB. Afterwards, we assessed in TCGA-SCKM whether XP or TIS score, single XP or immune infiltration gene expression, or TMB could be prognostic for survival. Except for the pre-processing, the specified workflow was repeated for multiple sample subgroups split by clinical parameters such as age, sample type or gender.

To evaluate the potential of the XP genes as a predictive biomarker for the ICI response, we utilized the two anti-PD-1 datasets and developed simple prediction models using Youden’s index and Xtreme gradient boosting.

All data analyses were performed using R version 3.6.3 (46). A p-value <0.05 was considered statistically significant in all analyses, and a p-value <0.005 was highly statistically significant.

Data Collection and Preprocessing

>Gene expression data used in this manuscript were obtained from TCGA-SKCM (http://cancergenome.nih.gov/, n = 464) (43) and two datasets of anti-PD-1 ICI cohort studies in melanoma patients, Hugo et al. (n = 28, GEO: GSE78220) (11) and Riaz et al. (n = 110, GEO: GSE91061) (12). The TCGA-SKCM dataset was reduced to n = 445 samples, which are fully annotated with clinical information, such as age, gender, and survival time. Likewise, we included samples of the other datasets after filtering for the mentioned clinical data availability and exclusion of on-treatment samples from the ICI cohorts, resulting in n = 26 (11) and n = 49 (12). All analyzed RNA-seq data were formatted as FPKM and log2 transformed. For TCGA-SKCM, somatic mutations were obtained from the TCGA data portal, and the TMB was calculated as log10 of the number of non-synonymous mutations per 50 Mb (package “maftools” v.2.2.10) (47). Responder (complete response [CR] or partial response [PR]) and non-responder (stable disease [SD] or progressive disease [PD]) were defined by RECIST criteria-based radiological response (7, 10). The clinical characteristics plus the scope of the computed scores of the utilized cohorts are listed in .

Clustering

In the process of clustering method selection, multiple clustering methods and distance metrics of hierarchical clustering were tested ( ). Clusters containing only one single gene were excluded because singe genes analysis of XP genes was performed apart. Hence, as the final XP clusters we selected the best performing partition with at least two genes per cluster, which was supported by the majority of all tested clustering methods and distance metrics.

Calculation of Scores

To identify a T cell-inflamed TME, we followed Ayers et al., based on the log2 transformed FPKM values; the TIS score was calculated as the weighted sum of the expression values of the 18 genes, enumerated in , applying the predefined weights derived by Ayers et al. (7, 9, 45). Considering the generation of scores based on the sum of signature related genes, we accordingly defined two XP gene cluster scores by summing up the expression values of genes in the same cluster.

Correlation Analysis

To assess the co-expression relationship between the considered genes, we cross-correlated the specified parameters. The Spearman rank correlation with p-value adjustment using Benjamini–Hochberg was performed by R package “psych” (v.2.1.3) (48) and visualized with package “ComplexHeatmap” (v.2.2.0) (49) using complete clustering with Euclidean distance for the dendrogram displayed at the columns.

Survival Analysis

For the survival analysis of the TCGA-SKCM data, we defined the overall survival (OS) as the time between melanoma diagnosis and the death or the last follow-up of the patient. The median follow-up was 669.50 days, while the survival status was decoded by 0 (alive) and 1 (dead). The constructed univariant Cox regression model predicted the overall survival from the continuous scores and gene expression values obtained from R packages “survival” (v.3.2-11) (50, 51) and “survminer” (v.0.4.9) (52). Kaplan–Meier analysis was used to calculate the survival probability of stratified patients, and the log-rank p-values for each analysis were given.

Response Prediction Model Construction

For anti-PD-1 response analysis, expression data of the two clinical cohorts (11, 12) were downloaded and reanalyzed using the Wilcoxon test and comparing expression levels of scores and single genes between responder and non-responder samples. The Youden index with associated ROC was determined for each parameter with R package “cutpointr” (v.1.1.0) (53, 54). The analysis was extended by multivariable predictive models for classification with the machine learning algorithm XGBoost (v.1.4.1.) (55) by partitioning the samples 75%/25% to training and testing data, respectively. The performed classification into responder and non-responder used “caret” (v. 6.0-86) (56) with the “xgbTree” method (55) and 10-fold cross-validation for combinations of multiple parameters.

Results

Heterogeneity of XP Gene Expression in TCGA-SKCM

First, we explored the nine XP genes XPA, ERCC3, XPC, ERCC2, ERCC4, ERCC5, DDB2, POLH, and ERCC1 in TCGA-SKCM and observed heterogeneous expression patterns. By unsupervised clustering, we could identify two XP gene expression clusters, referred to as XP gene clusters 1 and 2 ( ). Cluster 1 comprised the genes XPA, ERCC4, and ERCC5, while cluster 2 included ERCC3, XPC, ERCC2, DDB2, POLH, and ERCC1. Remarkably, the same clustering appeared if the cohort had been priorly divided by sample type (primary or metastatic, ), age (younger or older than median age of 58, ), or gender (male or female, ). Genes of both XP clusters and their function in the NER pathway and of POLH are summed up in .

Figure 2

The heatmap of log2 transformed FPKM values of the nine XP genes for all patient samples in TCGA-SKCM. The columns are clustered by hierarchical clustering with Manhattan distance and complete linkage.

Table 1

XP genes with corresponding clusters, their encoding proteins and their functionality in the NER pathway and translesion synthesis, respectively.

Genes	Cluster Membership	Corresponding Proteins	Main Function
XPA	1	XPA	Involved in multiple NER steps, e.g., DNA damage verification; interacts with almost all other NER proteins
ERCC4	1	XPF	DNA lesion excision in a complex with ERCC1 at 5ʹ end from the lesion
ERCC5	1	XPG	DNA lesion excision at 3ʹ end from the lesion
ERCC3	2	XPB	DNA damage verification as TFIIH basal transcription factor complex DNA helicase subunit
XPC	2	XPC	DNA-damage recognition in GG-NER
ERCC2	2	XPD	DNA damage verification as TFIIH basal transcription factor complex DNA helicase subunit
DDB2	2	XPE	Auxiliary DNA-damage-recognition factor in GG-NER
POLH	2	XPV	DNA polymerase η, which is an enzyme of translesion synthesis, that bypasses unrepaired DNA damage
ERCC1	2	ERCC1	DNA lesion excision in a complex with XPF at 5ʹ end from the lesion

GG, global genome; NER, nucleotide excision repair; TFIIH, transcription initiation factor IIH.

Altogether, median XP gene expression did not vary significantly in the analysis of subgroups. However, ERCC4, XPC, and POLH were expressed substantially greater in metastatic samples, whereas DDB2 was expressed considerably higher in primary tumors. Subdividing the whole TCGA-SKCM cohort by median age, we found that in melanoma tissue from younger patients, XP cluster 1 genes and also XPC and DDB2, belonging to XP cluster 2, were expressed to a relatively higher extent ( ). In samples from female patients, all XP genes, except ERCC1 and ERCC2, were expressed to a greater extent than in males.

Correlation Analysis Between XP & Immune Infiltration Genes, TMB and Computed Scores

Next, we investigated the correlation of XP genes and associated XP expression clusters to well-established predictive biomarkers of anti-PD-1 ICI response ( ).

Figure 3

Correlation analysis between expression of the nine XP genes, of computed XP and TIS scores and of single immune infiltration genes.

The expression of XP cluster 1 score with the 18-gene immune infiltration TIS score (p = 0.00034; R = 0.1793) as well with its single genes CD274 (p = 6.858 e−07; R = 0.244), CXCL9 (p = 1.640 e−06; R = 0.237), CXCL13 (p = 2.666 e−06; R = 0.232), and CD8A (p = 0.0107; R = 0.131) showed weak but significant positive correlations. Likewise, the XP cluster 1 genes XPA and ERCC5 were significantly positively correlated with TIS score and the above-mentioned immune infiltration genes. However, expression of ERCC4 was only significantly correlated with CD274 (p = 0.004; R= 0,146). Expression of XP cluster 2 score, on the other hand, was negatively correlated with the expression of TIS score (p = 0.037; R = −0.108), CXCL9 (p = 0.013; R = −0.127), CD8A (p = 0.012; R = −0.129), and CXCL13 (p = 0.01; R = −0.133). Importantly, the XP cluster 2 score genes ERCC1 and ERCC2 were individually highly significantly negatively correlated with TIS score (p = 0.0002, R = −0.186 and p = 0.0004, R = −0.176), immune infiltration genes and the XP cluster 1 score ( ).

Tumor mutation burden (TMB) had weak positive correlations only with CD274 (p = 0.01342; R = 0.127), CXCL9 (p = 0.038; R = 0.1078), and also with the XP cluster 2 genes ERCC3 (p = 0.00196; R = 0.157) and ERCC2 (p = 0.027; R = 0.114). Importantly, no significant negative correlations between XP genes and TMB were observed.

When considering primary and metastatic samples separately, some differences were evident ( ): in primary samples only (n = 96), XP cluster 1 and cluster 2 scores had a positive correlation (p = 0.017; R = 0.279). Furthermore, the expression of the genes ERCC3 and XPC was closely correlated with the expression of the XP cluster 1 score and its genes XPA, ERCC4, and ERCC5. Now, with a few exceptions, there were no significant correlations between XP genes and TIS score or immune infiltration genes but a positive correlation of TMB with TIS score (p = 0.00993; R = 0.3), CD274 (p = 0.0048; R = 0.323), CXCL9 (p = 0.0046; R = 0.327), CD8A (p = 0.004; R = 0.325), and CXCL13 (p = 8.56 e−05; R = 0.422). Correlation of the by far larger group of metastatic samples (n = 349) revealed almost the same picture as for the whole group.

Further splitting by age and gender led to identical correlation patterns of XP gene clusters 1 and 2 with TIS score and immune infiltration genes, as we had observed for the whole TCGA-SKCM cohort ( ). Of note, the significant positive correlation of TMB with CD274 was only detected if considering just males or the younger subgroup of patients.

XP & Immune Infiltration Genes, TMB, and Computed Scores as Prognostic Biomarkers for Survival

The great majority of TCGA-SKCM samples were obtained in the pre-ICI era. Only two patients received anti-PD-1 ICI after acquiring their tumor, but the removal of these two patients did not lead to significantly different results ( ) (43). Hence, we sought to analyze if there is a linear association between the expressions of the XP genes or cluster scores with survival of patients in TCGA-SKCM, and independent of ICI. Additionally, we analyzed the prognostic value of parameters predictive of anti-PD-1 ICI response: The TIS score, selected single score genes (CD274, CXCL9, CD8A), CXCL13 TMB, and age.

demonstrated that neither single XP cluster 1 nor cluster 2 score were associated with survival in TCGA-SKCM. Considering the univariate Cox regression, out of the single XP genes, only ERCC3 expression (hazard ratio, HR = 0.66, p = 0.043) revealed a significant association with survival ( ). In contrast, TIS score (HR: 0.87), CXCL13 (HR = 0.93), TMB (HR = 0.74), and age (HR = 1.03) were all prognostic. TIS score and age also remained significant after segregation by the median. illustrate the corresponding Kaplan–Meier survival curves (log ranked p-values for age, p = 0.0014; TIS score, p = 0.0077).

Figure 4

Overview of the univariate cox regression analysis for all TCGA-SKCM patients (A). The bar indicates the reference Hazard ratio of 1. The patients split by median age into older patients (B) and younger patients (C) show different Hazard ratios for the same parameters.

The same analysis was repeated independently for age-divided sub-cohorts to decrease the influence of age as a dominant factor. Through this fractioning, the median overall survival dropped from 3,424 days to 1,927 days in the older patient group, while it increased to 4,634 days for the younger patients.

In those subgroups, we noted some differences ( and – ). Although age and TIS score remained significant prognostic parameters, ERCC3 (HR = 0.48, p = 0.0084) was the best predictor for survival and the only prognostic factor in the younger cohort. In the subgroup of older patients, CXCL13 (HR = 0.93) and, newly, XP gene cluster 1 (HR = 0.87) were additional prognostic parameters.

XP & Immune Infiltration Genes and Computed Scores as Predictive Biomarkers for Response to Anti-PD-1 ICI

Because of the remarkable anecdotal benefit of XP patients to anti-PD-1 ICI, we analyzed XP gene expression as a predictive biomarker for response to anti-PD-1 ICI in two pooled publicly available melanoma cohorts (n = 75) (11, 12).

The distribution of the responding (complete or partial response) patients differed significantly from the non-responders based on the XP cluster 1 score (p = 0.015), with a higher score indicating a greater response ( ). A similar significant difference between these two groups was also applied for the single XP cluster 1 gene ERCC5 (p = 0.026) ( ). Importantly, XP cluster 2 score, TIS score, and single genes indicative of immune infiltration, except CD27 and PSMB10, were not significantly associated with response in the pooled anti-PD-1 melanoma cohorts ( ).

Figure 5

Boxplot of ICI response data (n = 75), compared with Wilcoxon test based on the expression of (A) XP cluster 1 score and (B) ERCC5.

To assess the predictive performance of these single genes for ICI response, we computed the Youden index for each parameter (a XP gene or signature) and compared the areas under the receiver operating characteristic curves (AUC) ( ) (7). In this analysis, ERCC5 (AUC = 0.660), XP cluster 1 score (AUC = 0.654), XP cluster 2 score (AUC = 0.632), and CD8A (AUC = 0.627) were the best performing variables. Except XPA (AUC = 0.532), POLH (AUC = 0.545), and CD274 (AUC = 0.533), all other parameters outperformed TIS score (AUC = 0.586).

Expansion of these restricted single parameter analyses by combining multiple variables from two to five possible parameters led to many response classification models. Due to infinitive values for two samples, this analysis was limited to 73 patients. The combination of two parameters had the best results if one immune infiltration gene (like CD274 or CXCL13) was combined with one XP gene (CD274_ERCC4, AUC = 0.7; CXCL13_ERCC5, AUC = 0.68), or if XP gene cluster 2 gene ERCC2 was combined with ERCC5 (AUC 0.69). All these combinations outperformed any combination involving the TIS score ( and ). The three-parameter-based analysis performed better than the combination of two parameters and revealed that combining XP genes provided the best classification triplet for response (ERCC3_XPC_ERCC4, AUC = 0.8; XPC_ERCC4_ERCC1, 0.75). ( and ). The prediction of the combination of four variables had the best AUC value of 0.85, and even outperformed the combination of five parameters ( , and ). Of note, it included the combination of one immune infiltration gene (CD274), two XP cluster 1 genes (XPA, ERCC4) and one XP cluster 2 gene (ERCC2).

Figure 6

ROC curves with AUCs of top 5 combinations of (A) 2 parameters, (B) 3 parameters and (C) 4 parameters for prediction of ICI response.

Discussion

Clustering by using mRNA gene expression levels can identify higher-level structures and relationships and establish a new molecular classification of tumors (57, 58). Furthermore, hierarchical clustering based on gene expression profiles (GEPs) can be used to, e.g., reveal immune competency or sensitivity to ICI treatment (9, 59).

Because melanoma has been the “model tumor” for the development of ICI and was also the first tumor entity in which ICI was approved, we focused our exploratory analysis on melanoma. Our study analyzed XP gene expression and deduced two different expression clusters (XP gene clusters 1 and 2, see and ) with heterogeneous functions in the NER pathway and translesion synthesis. The proteins encoded by the XP cluster 1 genes ERCC4, i.e., XPF, and ERCC5, i.e., XPG, are responsible for the dual incision of DNA damage. Cluster 1 further includes XPA, the central NER coordinator, because of its interaction with almost all other NER factors (39, 60). Accordingly, it also interacts with ERCC1, an endonuclease and fulfills its function as a heterodimer with XPF (61, 62). Surprisingly, ERCC1 is part of XP cluster 2 instead of cluster 1, and there was a strong negative correlation between ERCC1 and ERCC4 expression in the TCGA-SKCM samples ( ). Besides, in NER, the ERCC1–XPF complex is involved in interstrand crosslink and DSB repair, and mutations in one of the two genes result in a very complex constellation of clinical symptoms (39, 63). In epithelial ovarian cancer, ERCC1 and ERCC4 expression correlated on mRNA and protein level with one another, however, ERCC1 mRNA was negatively correlated, and ERCC4 mRNA was unrelated with its protein expression, suggesting a posttranscriptional mode of regulation (64). Because protein expression data was unavailable for most of our nine XP genes, we could not expand our analysis to protein correlations, which might have revealed substantially different clusters due to extensive posttranslational modifications in NER (39). In global genome-NER, damage recognition is performed by XPC and involves XPE (encoded by DDB2), which are both parts of XP cluster 2. Different genes of cluster 2 include ERCC2 (codes for XPD) and ERCC3 (codes for XPB), which are DNA helicases and part of transcription initiation factor IIH complex verifying DNA damage lesions. POLH, whose defect leads to XP variant, codes for the DNA error-prone polymerase eta involved in translesion synthesis (39, 65).

Subsequently, we investigated the correlation of the XP genes and clusters with well-known, ICI predictive biomarkers of a T cell-inflamed TME: The TIS score, namely, its single genes CD8A, CXCL9, CD274, and the recently postulated biomarker CXCL13, that could be characteristic of clonal neoantigen-reactive CD8 T infiltrating lymphocytes (6, 7, 9, 10). Of note, we revealed that XP cluster 1 and its genes XPA and ERCC5 had a highly significant positive correlation with the TIS score and all other immune infiltration associated genes ( ). In line with this, Boonstra et al. (66) compared UVB suppression of ConA-induced IFN-γ production in XPA, XPC, and CSB deficient mice and demonstrated that only XPA mice showed a substantial reduction of IFN-γ production by UVB. Regarding gene correlations in TCGA-SKCM, our results are different from those of BER/SSB repair genes, which more homogeneously and almost exclusively present a negative correlation with CD274 expression (26).

A negative correlation with CD274, TIS score, or the other genes representative of immune infiltration was identified for the XP cluster 2 genes ERCC1 and ERCC2, which was true especially for samples of younger patients. However, neither the expression of ERCC1, ERCC2, nor of other XP genes was negatively correlated with TMB. Likewise, the frequency of XP gene mutations with median TMB values was not correlated in different cancers (41). Noteworthy, from the XP genes, ERCC3 expression even had a highly significant positive correlation with TMB in samples of the whole cohort and the younger subcohort ( and ). Taken together, the results from Hsiehchen et al. (41) and ourselves illustrate that the implications for tumor immunogenicity through XP gene mutations and diverse expression are presumably more complex than being based solely on the more abundant generation of neoantigens caused by somatic tumor mutations. Admittedly, we cannot precisely determine if the XP gene expression we assessed is preferentially constituted by tumor or other, e.g., immune cells in the TME (67). Circadian XP gene expression particularly affects XPA and could potentially impact the results of our analysis; however, it seems not to be relevant in actively proliferating tissues as tumors (68, 69).

Just recently, survival prognosis in melanoma was correlated with immune-related gene signatures (70–72). Of note, Danaher et al. (73) found that TIS was also highly statistically significantly prognostic in TCGA-SKCM, limiting its predictive value in melanoma patients. Accordingly, we analyzed the predictive role of all markers and found out that besides TIS, CXCL13, TMB, age, and ERCC3 were prognostic. Age and TMB have been revealed as prognostic biomarkers in different studies and cancer entities before, though results for TMB are not homogeneous and depend on the used thresholds (74–77). Correspondingly, in our univariate analysis with segregated parameters for high and low values based on medians, we found no difference between TMB groups ( ). CXCL13 was suggested as a prognostic biomarker in melanoma before, but the correlation of ERCC3 with survival needs further validation (78, 79).

Notably, after dividing the cohort by age, in younger patients, ERCC3 (HR 0.48) and TMB (0.66) revealed the lowest HR values; however, for TMB, it was not significant. This observation, together with the positive correlation between ERCC3 and TMB in younger patients, suggests that ERCC3 might have specific relevance for disease progression, especially in the young. In contrast, earlier reports have suggested that the presence of an intense immune infiltrate in older persons could have more prognostic weight (75). This assumption might explain that in our analysis, CXCL13 expression, as a specific marker of exhausted T cells, was only prognostic in melanoma tissue of older when subdividing by patient age. The same accounted for XP cluster 1, which in samples of older was highly correlated with TIS score and other immune infiltration genes ( and ).

An ongoing challenge is the identification of reliable biomarkers predictive of ICI response. Recent efforts leave single parameters and evolve into combinatorial biomarkers (6, 7, 10). Especially, the combination of tumor intrinsic factors, like TMB, and genes representative of immune infiltration in the TME, like TIS or CXCL13, show promise in exploratory studies (7, 9, 10, 80). Hence, we tested the single and combinational predictive performance of our parameters in two pooled anti-PD-1 cohorts. The existing immune-specific biomarkers were only of limited value and constructed prediction models (10). Importantly, we observed that the combination of either an immune-infiltrating gene (CD274) with three XP genes, or the combination of only three XP genes from both clusters provided the best ROC and AUC values ( ). Of note, the prediction was not improved by extension to five parameters; hence, we did not test further combinations of even more variables. Because of our small, pooled cohort and to avoid overfitting, we did conduct only a split of our data in training and testing set for our prediction models. We, therefore, admit that our constructed models lacked robustness to establish new ICI predictive biomarkers.

Table 2

Combination of the best 25 performing models, based on the AUC values, across the different number of used parameters.

Parameter for predictive model	AUC	# parameters
CD274_XPA_ERCC2_ERCC4	0.85	4
ERCC3_XPC_ERCC4	0.80	3
TIS score_ERCC2_ERCC1_cluster 1 XP score	0.78	4
TIS score_CXCL13_CXCL9_CD274_cluster 1 XP score	0.78	5
CXCL9_CD8A_CD274_ERCC2_cluster 1 XP score	0.78	5
CXCL9_CD8A_XPC_ERCC1_cluster 1 XP score	0.78	5
CD8A_CD274_ERCC3_ERCC4_cluster 1 XP score	0.78	5
XPA_ERCC2_ERCC1_cluster 1 XP score_cluster 2 XP score	0.78	5
CXCL9_ERCC3_ERCC4_cluster 1 XP score	0.77	4
CXCL9_CD8A_XPA_ERCC2_ERCC5	0.77	5
CXCL9_CD274_XPA_ERCC4_ERCC1	0.77	5
CD8A_CD274_DDB2_ERCC1_cluster 1 XP score	0.77	5
CD274_XPA_XPC_ERCC4_ERCC1	0.77	5
CD274_DDB2_ERCC4_ERCC1_cluster 1 XP score	0.77	5
CD8A_ERCC2_ERCC1_cluster 1 XP score	0.76	4
ERCC3_ERCC4_ERCC1_cluster 1 XP score	0.76	4
XPC_ERCC4_ERCC1	0.75	3
TIS score_XPA_ERCC2_cluster 1 XP score	0.75	4
CXCL9_CD8A_XPC_cluster 1 XP score	0.75	4
CXCL9_CD274_ERCC1_cluster 2 XP score	0.75	4
TIS score_CXCL13_CD8A_ERCC4_ERCC1	0.75	5
TIS score_CD8A_CD274_ERCC1_cluster 1 XP score	0.75	5
TIS score_CD8A_CD274_cluster 1 XP score_cluster 2 XP score	0.75	5
TIS score_CD8A_ERCC3_ERCC4_cluster 1 XP score	0.75	5
TIS score_CD274_ERCC2_ERCC4_cluster 1 XP score	0.75	5

Our study has several limitations. First, we restricted our analysis to melanoma, and the samples sizes of the two clinical cohorts that we analyzed were small, limiting our results’ comparability. For example, TIS score and genes indicative of immune infiltration, which, except for CD8A, performed poor as singular splitting parameters in our pooled data of two anti-PD-1 cohorts (11, 12), were significantly predictive in the study of Cristescu et al. (AUC of TIS score = 0.638) (7, 10). Second, due to the standardized pre-processing, the sample size of the two cohorts had to be further reduced. For the anti-PD-1 cohorts, we did not analyze survival data restricting our analysis to ICI response, which is not an appropriate measure of long-term treatment benefit. Third, 11 of 26 patients included from the cohort of Hugo et al. (11) and 26 of 49 patients in the cohort of Riaz et al. (12) were not treatment-naïve and had received prior MAPK inhibitor treatment or anti-CTLA-4 ICI, respectively, before sample acquisition, potentially influencing gene expression. Nevertheless, primary and acquired resistance constitutes a major problem in the systemic therapy of melanoma, suggesting the analysis of the XP gene clusters in additional patients with therapy failure (2, 3, 6, 81). Fourth, in the two clinical cohorts that we considered, mainly metastatic tissue was analyzed, and our analysis was made regardless of gender and age, which could have had an influence, especially on XP GEPs. Fifth, due to limited data, TMB and clonal TMB, which were the best performing predictive markers in the meta-analysis of Litchfield et al. (10), could not be assessed in our study. Likewise, we did also not assess the expression of genes in other DDR pathways.

Despite all these limitations, our analysis provides significant new findings that deserve attention: Firstly, XP genes are expressed in two heterogeneous clusters in melanoma. Secondly, these clusters correlate differentially with markers of a T cell-inflamed TME, and correlations depend to a certain degree on melanoma tissue origin (primary vs metastatic), age, and gender. Thirdly, a higher ERCC3 expression could be associated with a better prognosis in melanoma, especially in younger patients. Lastly, the differentiated consideration of XP gene expression in the TME and its combination with established ICI predictive biomarkers could be useful in predicting anti-PD-1 ICI response in melanoma and should be explored by further studies.

Data Availability Statement

The original contributions presented in the study are included in the article/ . Further inquiries can be directed to the corresponding author.

Author Contributions

SF and AT designed this study, conducted the analysis, and wrote the article. MH, SE, OW and GF provided conceptual advice and critically reviewed the article. All authors contributed to the article and approved the submitted version.

Funding

AT is supported by the “Hiege Foundation—the German Skin Cancer Foundation” and the “Forschungsförderung der Universitätsmedizin Rostock (FORUN)”. Support for computing equipment was provided by the European Union (EFRE, “Europäischer Fonds für regionale Entwicklung”). The funders had no role in the study design, collection, analysis, interpretation of the data, writing of the manuscript, or the decision to submit the paper for publication.

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.