Genetic and transcriptional insights into immune checkpoint blockade response and survival: lessons from melanoma and beyond

Abstract

Background

Integration of immune checkpoint inhibitors (ICIs) with non-immune therapies relies on identifying combinatorial biomarkers, which are essential for patient stratification and personalized treatment.

Methods

We analyzed genomic and transcriptomic data from pretreatment tumor samples of 342 melanoma patients treated with ICIs to identify mutations and expression signatures associated with ICI response and survival. External validation and mechanistic exploratory analyses were conducted in two additional datasets to assess generalizability.

Results

Responders were more likely to have received anti-PD-1 therapy rather than anti-CTLA-4 and exhibited a higher tumor mutation burden (both P < 0.001). Mutations in the dynein axonemal heavy chain (DNAH) family genes, specifically DNAH2 (P = 0.03), DNAH6 (P < 0.001), and DNAH9 (P < 0.01), were enriched in responders. The combined mutational status of DNAH 2/6/9 effectively stratified patients by progression-free survival (hazard ratio [HR]: 0.69; 95% confidence interval [CI] 0.51–0.92; P = 0.013) and overall survival (HR: 0.58; 95% CI 0.43–0.78; P < 0.001), with consistent association observed in the validation cohort (HR: 0.28; 95% CI 0.12–0.61; P < 0.001). DNAH-altered melanomas exhibited upregulation of chemokine signaling, cytokine-cytokine receptor interaction, and cell cycle-related pathways, along with elevated expression of immune-related signatures in interferon signaling, cytolytic activity, T cell function, and immune checkpoints. Using LASSO logistic regression, we identified a 26-gene composite signature predictive of clinical response, achieving an area under the curve (AUC) of 0.880 (95% CI 0.825–0.936) in the training dataset and 0.725 (95% CI 0.595–0.856) in the testing dataset. High-risk patients, stratified by the expression levels of a 13-gene signature, demonstrated significantly shorter overall survival in both datasets (HR: 3.35; P < 0.001; HR: 2.93; P = 0.002).

Conclusions

This analysis identified potential molecular determinants of response and survival to ICI treatment. Insights from melanoma biomarker research hold significant promise for translation into other malignancies, guiding individualized anti-tumor immunotherapy.

Background

Melanoma is a deadly form of cancer that primarily affects the skin in various locations, with an estimated 8,390 deaths expected in the United States in 2024 [1]. Over the past 30 years, melanoma incidence has doubled, driven by factors such as increased exposure to ultraviolet light, continued use of tanning beds, and better awareness and detection methods [2]. Notably, melanoma is highly responsive to immune modulation due to its strong immunogenic properties. The standard of care for metastatic melanoma has evolved significantly from historical immune therapies, such as IFNα for high-risk primary or regional disease and high-dose IL2 for advanced/metastatic disease, to more effective treatments involving monoclonal antibodies (mAb) targeting cytotoxic T-lymphocyte associated protein 4 (CTLA-4), such as ipilimumab (approved by the Food and Drug Administration [FDA] in 2011), and/or programmed death 1 (PD-1) inhibitors, including nivolumab and pembrolizumab (approved by the FDA in 2014). These advancements have revolutionized the treatment landscape for patients with advanced melanoma, markedly improving survival outcomes compared to the pre-ipilimumab era. The KEYNOTE-001 trial reported a 5-year overall survival (OS) rate of 34% to 41% with pembrolizumab monotherapy [3]. Previous analyses of the CheckMate 067 trial demonstrated that nivolumab, either alone or in combination with ipilimumab, can achieve durable disease control in patients with advanced melanoma [4, 5]. Recent data from this trial, with a minimum follow-up of 10 years, reported a median OS of 71.9 months for nivolumab plus ipilimumab, 36.9 months for nivolumab monotherapy, and 19.9 months for ipilimumab monotherapy [6]. Notably, for patients who were alive and progression-free at three years, the 10-year melanoma-specific survival rates were 96% with nivolumab plus ipilimumab, 97% with nivolumab, and 88% with ipilimumab, underscoring the potential for long-term remission with combination immunotherapy. However, the higher incidence of treatment-related adverse events with ipilimumab-containing therapies compared to nivolumab monotherapy highlights the need for alternative strategies. Emerging treatments, such as combination therapy with anti-PD-1 and anti-lymphocyte-activation gene 3 (LAG-3) agents, are being explored, as both immune checkpoints contribute to T-cell exhaustion [7, 8]. The RELATIVITY-047 trial demonstrated that dual inhibition of LAG-3 and PD-1 with relatlimab and nivolumab expanded treatment options for metastatic melanoma, achieving a median progress-free survival (PFS) of 10.1 months without additional safety concerns compared to nivolumab monotherapy [9].

Despite tremendous progress and remarkable clinical successes in harnessing the immune system against cancer, several challenges persist. One major obstacle is understanding the immunobiology of tumors that exhibit primary or acquired resistance, as nearly 50% of melanoma patients treated with immune checkpoint inhibitors (ICIs) develop resistance [10, 11]. Additionally, elucidating the mechanisms underlying immune-related adverse events (irAEs) remains an important area of research. Furthermore, with the primary objective of ICI therapy being to eliminate the disease and/or achieve long-term durable responses, the identification of predictive and prognostic biomarkers holds great potential for tailoring treatment decisions to individual patients. Research on pretreatment tumors has uncovered potential biomarkers associated with ICI response in melanoma, including increased CD8 T-cell infiltration, interferon-gamma (IFNγ) gene signatures, tumor mutation burden (TMB), and PD-L1 expression [12]. Recent studies have explored additional biomarkers. Among these, microsatellite instability (MSI) has emerged as a predictor for immunotherapy efficacy [13]. Tumors with deficient mismatch repair (dMMR) and high MSI accumulate extensive mutations, leading to the production of neoantigens that enhance anti-tumor immune responses. Additionally, systemic inflammation markers, such as the neutrophil-to-lymphocyte ratio (NLR), have been linked to ICI outcomes in melanoma. A high NLR reflects an increased neutrophil population skewed toward the immunosuppressive N2 phenotype, correlating with poor prognosis [14]. Beyond tumor-intrinsic factors, external influences, such as concomitant medications can also modulate ICI efficacy. Cortellini et al. demonstrated the immunomodulatory effects of concomitant medications in patients receiving PD-1/PD-L1 inhibitors across various malignancies [15]. Furthermore, Mallardo et al. reported that the concurrent use of cetirizine with anti-PD-L1 therapy may promote M1 macrophage polarization via the IFNγ pathway, significantly improving PFS and OS in advanced melanoma patients [16]. However, none of these biomarkers have consistently predicted patient response to ICIs. For those who do not respond to ICIs, there is a growing need for identifying new combinatorial biomarkers, as single biomarkers often fail to capture the complex and dynamic interplay between the tumor and the host immune system.

In this study, we analyzed pretreatment tumor samples from a cohort of 342 melanoma patients treated with ICI therapies. Whole-exome sequencing (WES) data were available for all patients, while a subset of 201 patients had whole transcriptomic sequencing (WTS) data. Validation and mechanistic exploration analyses were performed using two independent external cohorts from publicly available datasets, namely the non-small cell lung cancer (NSCLC) cohort and The Cancer Genomic Atlas (TCGA) pan-cancer cohort. Given the immunogenic nature of melanoma and its suitability for identifying immunotherapy biomarkers, we believe the insights gained from melanoma may apply to other cancer types in the context of immunotherapy.

Methods

Patients and cohorts

This multicohort study initially included 356 subjects, integrating clinical, genomic, and transcriptomic data from four independent cohorts of melanoma patients treated with ICIs, obtained from publicly available datasets [17,18,19,20] (Fig. 1). Patients were excluded if they (1) lacked a documented tumor response assessment, (2) had tumor samples collected during or post-treatment instead of pretreatment, or (3) had an unknown overall survival status. A detailed overview of the patient inclusion process is provided in Figure S1. As a result, the discovery cohort comprised 342 melanoma patients treated with ICIs, all of whom had whole-exome sequencing (WES) data, while 201 had whole transcriptomic sequencing (WTS) data that passed quality control. This cohort was used to identify genes and expression signatures predictive of ICI response and survival. For validation, we combined clinical and WES data from two immunotherapeutic cohorts of NSCLC patients, referred to as the NSCLC cohort [21, 22]. Additionally, clinical, genomic, and transcriptomic data from the TCGA pan-cancer analysis [23], encompassing 942 patients with primary tumors profiled by whole-genome sequencing (WGS) and WTS, were utilized for external mechanistic exploration. This study was conducted in accordance with the principles of the Declaration of Helsinki.

Fig. 1

Overview of study design. Flow diagram outlining patient inclusion and cohorts, with subanalyses addressing specific clinical questions

Assessment of clinical outcomes

The evaluation of ICI response was based on the immune-related Response Evaluation Criteria in Solid Tumors (irRECIST) [24,25,26]. Responders were defined as patients who achieved complete or partial responses or had stable disease lasting more than 6 months, whereas non-responders were those who derived no durable benefit, exhibiting either progressive disease or stable disease lasting 6 months or less. Progression-free survival (PFS) was calculated from the date of the first application of ICI treatment to the date of documented disease progression or death. Overall survival (OS) was calculated from the start of ICI treatment to the date of death from any cause. For patients lost to follow-up, data were censored on the last date they were known to be progression-free or alive.

Differential expression analysis

Gene-level expression values, reported as Transcripts Per Million (TPM), were normalized using the min–max normalization method. Protein-coding genes with a TPM value greater than zero in more than 50% of samples were selected for downstream analysis. Differentially expressed genes (DEGs) were identified by comparing geometric means for individual genes between subgroups using the Wilcoxon rank-sum test. To correct for multiple comparisons across all tested genes, P-values were adjusted using the Benjamini-Hochberg (BH) method to control the false discovery rate (FDR). DEGs were defined as genes with an FDR-adjusted P-value (padj) < 0.05. Pathway enrichment analysis for Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and Gene Set Enrichment Analysis (GSEA) was performed using the R package “ClusterProfiler”.

Gene expression signature analysis

Five publicly available gene expression signatures related to ICI response were analyzed in this study: the IFNγ signature [27], expanded IFNγ signature [27], cytolytic activity [28], T-cell signature [29], and immune checkpoint signature [30] (Table S1). A composite immune-related signature, consisting of 28 genes commonly reported across these five signatures, was also analyzed. Sample-wise scores for each gene signature were calculated from RNA-seq data using TPM values and logistic regression. Genes with missing expression data were excluded from the calculations of gene signature scores.

Immune cell population fraction analysis

The proportions of immune cells were estimated by analyzing normalized gene expression, measured as TPM values, using the CIBERSORT algorithm. The “Tidyestimate” function in R was employed to calculate immune, stromal, and ESTIMATE scores, providing insights into the composition of immune and stromal cells within the tumor microenvironment (TME).

Statistical analysis

All statistical analyses were performed in R (version 4.1.3). Categorical variables were compared using Chi-square tests or Fisher’s exact tests, as appropriate, while continuous variables were analyzed using Wilcoxon rank-sum tests. Non-parametric tests were chosen for their robustness in handling skewed distributions, small sample sizes, and potential outliers. The Point-Biserial test assessed the correlation between normalized gene expression and clinical response. Receiver operating characteristic (ROC) curves were generated to assess genes associated with ICI response, with the area under the curve (AUC) calculated to measure the discriminatory ability of the predictor. The optimal cut-point was determined using the Youden index to maximize diagnostic accuracy while maintaining a balance between sensitivity and specificity. Kaplan–Meier curves compared the survival across different subgroups, and the log-rank test was used to assess statistical significance. Cox proportional hazard models assessed the association between individual gene mutations and survival, estimating hazard ratios (HR) with 95% confidence intervals (CI). For predictive modeling, we applied LASSO logistic regression and LASSO Cox regression for classification and prognostication, respectively. The RNA sequencing data was randomly split into training (N = 140) and test (N = 61) datasets in a 7:3 ratio. The training dataset was used for candidate gene selection and hyperparameter tuning, while the test set was reserved for performance evaluation and generalizability assessment. The penalty parameter was optimized using tenfold cross-validation with minimal partial likelihood of deviance. All reported P-values were two-tailed, and statistical significance was defined as P < 0.05 (*P < 0.05, **P < 0.01, ***P < 0.001).

Results

Patient characteristics

The discovery cohort (N = 342) included clinical, genomic, and transcriptomic data from four independent cohorts of melanoma patients receiving ICI treatment (Fig. 1). Male patients were nearly twice as prevalent as female patients, and 72.5% (248/342) were metastatic patients with M1c disease (Table S2). For those receiving anti-CTLA-4 therapies, ipilimumab, with or without dacarbazine (a chemotherapy drug) or vemurafenib (a BRAF inhibitor), was the primary choice. Nivolumab and pembrolizumab were used in 17.8% and 32.7% of patients, respectively, to target the PD-1/PD-L1 axis. Among patients with available serum lactate dehydrogenase (LDH) data, 158 had normal levels, while 126 exhibited elevated LDH. BRAF or NRAS mutations were identified in 53.8% of patients in the discovery cohort. Prior to ICI initiation, 33.3% (114/342) of patients were treatment-naïve, whereas the majority (57.0%; 195/342) had received at least one previous therapy. Notably, prior treatment exposure and higher lines of therapy were associated with reduced ICI efficacy following ICIs, possibly due to treatment resistance, immune suppression, or the inherently lower responsiveness of patients with more aggressive diseases requiring multiple prior therapies. Responders were also more likely to have received anti-PD-1 therapy compared to anti-CTLA-4 treatment and had tumors with a higher tumor mutation load (both P < 0.001) (Figure S2 A, B). All patients in the discovery cohort had OS data, with a median OS of 22.4 months (95% confidence interval [CI] 18.1–30.8 months) and survival probabilities of 63.2% (95% CI 58.2–68.5%) and 48.5% (95% CI 43.3–54.3%) at 12 and 24 months, respectively (Figure S2 C). Among patients with available PFS data (71.6%, 245/342), the median PFS was relatively short at 3.3 months (95% CI 3.0–4.1 months), with 1-year and 2-year PFS rates of 29.2% (95% CI 24.0–35.5%) and 23.8% (95% CI 18.9–29.8%), respectively (Figure S2D).

Genomic correlates of response to ICI treatment

Patients were categorized into responders and non-responders to identify genomic correlates of ICI response. As shown in Fig. 2A, the most frequently mutated genes in melanoma patients receiving ICI treatment were TTN (70%), MUC16 (67%), DNAH5 (51%), PCLO (47%), BRAF (42%), and LRP1B (41%). Figure 2B displays the top 30 mutational features with significantly different mutational frequencies between the two groups. Notably, 8 mutations, namely PCLO, LRP1B, HYDIN, CSMD2, CSMD3, XIRP2, DNAH9, and DNAH6, showed the most significant variations, and these mutations were more prevalent in responders compared to non-responders. Conversely, gene amplifications were more common in non-responders, with amplifications in MIR663 A and NCOR1P1 being exceptions, as they were more prevalent in responders (Fig. 2C). At the pathway level, alterations in the NOTCH, WNT, HIPPO, and NRF2 pathway may be associated with ICI response, while other oncogenic signaling pathways showed no significant differences (Fig. 2D). Functional analysis of genes with significantly different mutation frequencies between groups revealed that responders were more enriched with mutations in genes involved in extracellular matrix (ECM)-receptor interactions, cytoskeleton organization in muscle cells, and calcium signaling pathway. In contrast, mutations more prevalent in non-responders were likely associated with ABC transporters and nucleocytoplasmic transporter activities (Fig. 2E, F).

Fig. 2

Genomic correlates of response to ICI treatment. A Mutational landscape of melanoma patients clustered by clinical response. B Comparison of mutational frequencies for the top frequently mutated features between responders (R) and non-responders (NR). The bar plot highlights genes with a ≥ 15% difference in mutational frequency between groups. C Copy number variant (CNV) analysis depicting the proportion of patients with focal-level CNV stratified by clinical response. D Prevalence of ten oncogenic signaling pathways in patients grouped by clinical response. E, F Functional annotation of genes with differential mutational frequencies between subgroups using KEGG (E) or GO analysis (F)

TMB has emerged as a promising biomarker for ICI efficacy. Based on the median mutation load across samples, we categorized tumors into TMB-high and TMB-low groups. Higher TMB was associated with improved ICI efficacy and a favorable prognosis in melanoma patients, potentially due to the increased presence of neoantigens, which may enhance immune recognition and targeting of cancer cells (Figure S3 A-C). The most significantly different features between TMB-high and TMB-low tumors were further analyzed to identify associated biological processes, cellular components, and molecular functions (Figure S3D-F). Similar to findings reported in subgroup analysis based on clinical response, TMB-high tumors were enriched for mutations affecting cytoskeletal motor activity, as well as the PI3K-Akt and MAPK signaling pathways. Notably, genes from the dynein axonemal heavy chain (DNAH) family, which are implicated in Huntington’s disease, exhibited a significant correlation with both clinical response and elevated TMB. These results suggest that the mutational status of DNAH family genes may serve as a surrogate marker for TMB and, potentially, for predicting clinical response in melanoma patients treated with ICIs.

DNAH mutations as a potential biomarker of ICI response and survival in solid tumors

The DNAH family of genes is known for their roles in the structure and function of cilia and flagella, which are crucial for cell motility and signaling [31]. While genomic variants of DNAH genes have been reported in various malignant tumors, their role in ICI response remains poorly understood. To address this research gap, we extracted genomic data from the discovery cohort for all 13 DNAH gene members. Mutation rates for all these DNAH genes between responders and non-responders were as follows: DNAH5 (57% vs. 46%), DNAH7 (39% vs. 31%), DNAH8 (38% vs. 28%), DNAH11 (34% vs. 28%), DNAH9 (37% vs. 23%), DNAH3 (31% vs. 27%), DNAH10 (25% vs. 27%), DNAH17 (28% vs. 23%), DNAH6 (35% vs. 14%), DNAH2 (24% vs. 15%), DNAH12 (19% vs. 13%), DNAH1 (18% vs. 13%), and DNAH14 (17% vs. 14%) (Figure S4 A). However, only three DNAH genes, particularly DNAH2 (P = 0.03), DNAH6 (P < 0.001), and DNAH9 (P < 0.01) exhibited significantly different mutation rates between the two response groups (Fig. 3A–C). Tumors with alterations in any of these three genes, referred to as “DNAH 2/6/9” or “DNAH-altered,” were more frequently observed in responders (Fig. 3D). Additionally, DNAH-altered tumors exhibited significantly higher tumor mutation burden compared to DNAH-unaltered tumors (Figure S4B).

Fig. 3

DNAH mutations as a potential biomarker for clinical response and survival. A Genomic alterations in DNAH2, DNAH6, and DNAH9 within the discovery cohort. B Pie chart showing the variant classification of DNAH 2/6/9 alterations, with data presented as percentages and counts. C Bar plots illustrate the proportion of patients with DNAH gene mutations, comparing responders to non-responders. D Comparison of the number of patients with DNAH alterations versus those without, stratified by clinical response. E, F Kaplan–Meier curves showing progression-free survival (PFS) of patients stratified by DNAH6 status (E) and DNAH 2/6/9 status (F). G–J Kaplan–Meier curves depicting overall survival (OS) in patients stratified by individual mutations in DNAH2 (G), DNAH6 (H), and DNAH9 (I), or any combinations of these mutations (J). K Cox proportional hazard models for PFS (top) and OS (bottom) in patients stratified by key clinical features. L Forest plot depicting the multivariate analysis (including DNAH alteration status and other clinicopathological factors) for PFS (top) and OS (bottom) in melanoma patients. M, N Kaplan–Meier curves showing PFS of patients stratified by DNAH6 status (M) and DNAH 2/6/9 status (N) in the NSCLC cohort

We then evaluated the prognostic impact of DNAH mutations in the clinical context of ICI treatment. Among all DNAH genes, DNAH6 emerged as the only gene significantly associated with PFS, with DNAH6-mutated patients exhibiting markedly longer PFS compared to DNAH6-wildtype patients (HR: 0.54; 95% CI 0.36–0.80; P = 0.002) (Fig. 3E). Additionally, the combined mutational status of DNAH 2/6/9 effectively stratified patients based on their PFS (HR: 0.69; 95% CI 0.51–0.92; P = 0.013) (Fig. 3F). Furthermore, five DNAH family genes showed significant associations with overall survival, including DNAH2 (HR: 0.66; 95% CI 0.44–0.98; P = 0.039), DNAH6 (HR: 0.48; 95% CI 0.32–0.73; P < 0.001), DNAH9 (HR: 0.64; 95% CI 0.46–0.90; P = 0.009), DNAH5 (HR: 0.67; 95% CI 0.51–0.90; P = 0.006), and DNAH12 (HR: 0.57; 95% CI 0.36–0.92; P = 0.019) (Fig. 3G–I; Figure S4 C, D). Patients with alterations in DNAH 2/6/9 also had significantly longer OS compared to those without such alterations (HR: 0.58; 95% CI 0.43–0.78; P < 0.001) (Fig. 3J).

PFS across key clinical subgroups favored DNAH-altered tumors over DNAH-unaltered tumors, including male patients (HR: 0.67; 95% CI 0.47–0.97; P = 0.032), those with distant metastases (HR: 0.67; 95% CI 0.49–0.91; P = 0.011), patients receiving anti-PD-1 therapies (HR: 0.67; 95% CI 0.45–1.00; P = 0.049), and those with elevated serum LDH levels (HR: 0.50; 95% CI 0.32–0.79; P = 0.003) (Fig. 3K). The prognostic impact of DNAH alterations was even more pronounced for OS (HR: 0.58; 95% CI 0.43–0.78; P < 0.001), demonstrating a significant survival benefit across subgroups. Notably, this association was independent of BRAF or NRAS mutational status. Additionally, the survival advantage was observed regardless of prior systemic therapy, though it did not reach statistical significance in treatment-naïve patients. While TMB was predictive of improved outcomes following ICI treatment, it did not retain its significance after adjusting for potential confounding factors (Fig. 3L). In contrast, multivariate analysis revealed that DNAH alteration was the most significant prognostic factor for favorable OS (HR: 0.61; 95% CI 0.43–0.87; P = 0.007) followed by anti-PD-1 treatment (HR: 0.64; 95% CI 0.46–0.89; P = 0.008). Elevated LDH levels (P < 0.001), prior therapy (P = 0.005), and distant metastases (P = 0.023) were significantly associated with worse prognosis in these patients. Collectively, these findings suggest that DNAH mutations may serve as a robust prognostic biomarker in various clinical contexts for melanoma patients treated with ICIs.

We hypothesized that the mutational status of these three DNAH genes may also have significant potential for application in other cancer types. To test this, we employed a validation cohort from publicly available datasets consisting of 91 NSCLC patients receiving ICI treatments. Worth noting that only PFS data were available for these patients, whose primary tumor samples were profiled by WES. Consistently, the presence of DNAH6 mutations, as well as alterations in DNAH 2/6/9, was associated with better prognosis in NSCLC patients, further reinforcing the favorable prognostic implications of these gene mutations in ICI-treated patients (Fig. 3M, N; Figure S4E). Nevertheless, since most of these DNAH alterations were missense mutations of unknown significance, functional analysis is needed to provide further insights into how these mutations affect protein function and their impact on tumor biology and response to ICI treatments.

Potential underlying immune response mechanisms in DNAH-altered tumors

To investigate the potential underlying mechanisms in tumors with DNAH 2/6/9 alterations that responded to ICI treatment, we compared gene expression profiles of DNAH-altered (n = 91) and DNAH-unaltered (n = 110) tumors using the discovery cohort. In total, we identified 476 and 570 upregulated and downregulated genes in DNAH-altered tumors, respectively (Fig. 4A). KEGG analysis revealed that genes with elevated expression in DNAH-altered melanomas were significantly enriched in immune-related pathways, including the chemokine signaling pathway, cytokine-cytokine receptor interaction, cell cycle, and p53 signaling pathways (Fig. 4B). In contrast, DNAH-unaltered tumors displayed higher expression of genes associated with glycosylphosphatidylinositol-anchor biosynthesis, propanoate metabolism, and olfactory transduction, though none of these pathways remained statistically significant after multiple testing corrections. Furthermore, we compared the expression of five publicly available immune-related gene signatures, as well as a composite immune signature comprising 28 genes identified in the discovery cohort, between DNAH-altered and DNAH-unaltered tumors (Fig. 4C; Table S1). All these immune signatures were significantly upregulated in DNAH-altered tumors, with the 28-gene immune signature showing the most pronounced difference between subgroups. Immune cell infiltration analysis also revealed that DNAH-altered tumors had higher fractions of monocytes and resting dendritic cells (P = 0.019 and P = 0.032, respectively), as well as a slightly elevated immune-related score based on gene expression data (Figure S4 F, G). However, PD-L1 expression in baseline tumors was not discriminated by DNAH alterations (Figure S4H).

Fig. 4

Transcriptional analysis of DNAH-altered versus DNAH-unaltered primary tumors. A Volcano plot showing genes with a significantly different expression between DNAH-altered and DNAH-unaltered melanomas. B KEGG pathway enrichment analysis of genes identified through differential expression analysis in the discovery cohort. C Expression of immune-related gene signatures, comparing their levels between DNAH-altered and DNAH-unaltered melanomas. D Stacked bar plot of the mutation status of DNAH genes (DNAH2, DNAH6, DNAH9) in the TCGA pan-cancer cohort. E Differentially expressed genes between DNAH-altered (n = 143) and DNAH-unaltered (n = 799) cancers within the TCGA pan-cancer cohort. F Top 10 pathways corresponding to upregulated genes in the DNAH-altered and DNAH-unaltered cancers, as determined by KEGG analysis. Pathways upregulated in DNAH-unaltered cancers are reported with adjusted P-values (padj) for statistical significance after multiple testing corrections or nominal P-values (pval) for non-significance

To assess whether the baseline gene expression profiles of DNAH-altered and DNAH-unaltered patients with other malignancies resemble those in melanoma, we analyzed transcriptional data from the TCGA pan-cancer patients stratified by mutations in DNAH2, DNAH6, and DNAH9 (Table S3). The analysis revealed moderate to high mutation rates in these genes across various cancers, with NSCLC exhibiting the highest rate exceeding 40% (Fig. 4D). In line with the discovery cohort, upregulated genes in DNAH-altered tumors were primarily enriched in pathways related to the cell cycle (padj < 0.001), DNA replication (padj < 0.001), Fanconi anemia (padj = 0.008), and homologous recombination (padj = 0.02) (Fig. 4E, F). Pathways previously enriched in DNAH-altered melanoma patients, such as the chemokine signaling pathway, p53 signaling pathway, and TNF signaling pathway, were also identified, though not statistically significant (Table S4). Collectively, these findings indicate that DNAH-altered malignancies may share a consistent mechanism involving the modulation of cell cycle-related pathways, which potentially contributes to improved survival outcomes in patients harboring these mutations.

Transcriptomic signature predicts clinical response to ICI and long-term survival

Previously reported immune-related gene expression signatures, such as PD-L1 expression and those associated with IFN signaling, cytolytic activity, T cell functions, and immune checkpoints, demonstrated significantly different expression between responders and non-responders in our discovery cohort; however, their efficacy in predicting response to ICI was limited (Figure S5 A, B). Analysis of immune cell infiltration patterns revealed that responsive tumors had a significantly higher immune score (P = 0.007), characterized by an increased fraction of monocytes and a decreased fraction of eosinophils (P = 0.035 and P < 0.001, respectively) (Figure S5 C, D). Despite these, none of the immune-related features yielded sufficient predictive power for ICI response (Figure S5E). This underscores the need for more robust expression biomarkers to reliably predict clinical response and overall survival in patients receiving ICI therapies.

To address this, we compared the gene expression profiles of melanoma patients stratified by clinical response, which identified 516 and 934 genes with significantly higher expression in responders and non-responders, respectively (Fig. 5A). As expected, tumors with a better response to ICI treatment exhibited elevated expression of genes related to immune response processes, including cytokine-cytokine receptor interaction, T cell receptor signaling pathway, B cell receptor signaling pathway, natural killer cell-mediated pathway (Fig. 5B). Conversely, non-responsive tumors showed higher expression of genes primarily associated with ribosomal activity and thermogenesis. GO and GSEA analyses further supported these findings, revealing greater immunogenicity in primary tumors from responders, while non-responsive tumors exhibited enhanced DNA repair and ribosomal function (Figure S5 F, G).

Fig. 5

Transcriptional signature predicting clinical response and overall survival. A Differences in mean gene expression between responders (R) and non-responders (NR). B Top 10 pathways enriched in responder (left) and non-responders (right) determined by KEGG analysis. Significant pathways with an adjusted P-value (padj) < 0.05 are labeled in red for the NR group. C Estimate effect of gene expression levels on ICI response and prognosis. Seven genes were identified as both predictive and prognostic (blue), while those identified solely as predictive or prognostic are marked in red and green, respectively. D The receiver operating characteristic (ROC) curve shows the predictive power of the 26-gene signature established using LASSO logistic regression in the training dataset (N = 140). E Prediction score distribution in patients stratified by ICI response. F Stacked bar plot shows the proportion of patients classified according to prediction score cutoff thresholds, as defined by the Youden index in the training dataset. G ROC curve for the predictive gene signature in the testing dataset (N = 41). H Prediction score distribution of patients in the testing dataset. I Proportion of patients in the testing dataset classified using the same cutoff threshold as in the training dataset. J, K Kaplan–Meier curves illustrating the overall survival of patients stratified by the median risk score, with the cutoff built from the training dataset using LASSO Cox regression modeling, shown for both the training (J) and testing datasets (K), respectively

To narrow down the scope of candidate gene expression features potentially correlated with clinical response and survival, we conducted a correlation analysis on the 1,450 genes that displayed significantly different expressions between responders and non-responders. From this, 84 genes with a correlation coefficient of at least ± 0.2 were further selected for LASSO logistic regression modeling. To develop and validate the prediction model, we split the RNA data into training and testing datasets. No significant differences in clinical features were observed between the two cohorts (Table S5). The LASSO model on the training set identified a 26-gene composite signature capable of predicting clinical response with an AUC of 0.880 (95% CI 0.825–0.936) (Fig. 5C, D; Table S6). The prediction score was significantly higher in non-responders compared to responders (Fig. 5E). Using the Youden index (0.5508) as the cutoff, patients with higher scores were more prevalent in non-responders (Fig. 5F). In the testing dataset, the composite signature retained an AUC of 0.725 (95% CI 0.595–0.856), with non-responders consistently showing significantly higher scores (Fig. 5G–I). Furthermore, multivariate LASSO Cox regression analysis identified 13 out of the 84 genes as being significantly associated with overall survival (Fig. 5C; Table S6). Patients were divided into high-risk and low-risk groups using the median risk score calculated from the training dataset as the cutoff. High-risk patients, as defined by this gene signature, exhibited notably shorter overall survival compared to low-risk patients in both the training (HR: 3.35; 95% CI 2.08–5.41; P < 0.001) and testing datasets (HR: 2.93; 95% CI 1.46–5.89; P = 0.002) (Fig. 5J, K).

Discussion

Although immune checkpoint blockade therapies are routinely employed in advanced/metastatic melanoma, tumor response rates remain suboptimal, and disease progression is sometimes inevitable. A critical focus of research is the identification of reliable predictive and prognostic biomarkers to enable personalized treatment, tailoring therapy to each patient’s genetic profile. In this study, we report potential associations of DNAH gene mutations, and gene expression signatures, to ICI response and prognosis in melanoma patients. These findings may have broader implications across various malignancies in the context of immunotherapy.

Two major categories of ICIs were included in our study. Here, we showed that melanoma patients treated with anti-PD-1 mAbs, such as nivolumab and pembrolizumab, exhibited better clinical responses compared to those receiving anti-CTLA-4 therapies, primarily with ipilimumab. This result was indeed aligned with expectations, as single-agent trials and randomized studies have demonstrated superior efficacy and lower rates of irAEs with anti-PD-1 therapies compared to ipilimumab [11, 32], leading to the approval of PD-1 inhibition as the standard first-line therapy for metastatic melanoma in 2014.

DNAH genes encode axonemal dynein heavy chain proteins, which are important components of microtubules, a key target in cancer treatments. Functionally, DNAH proteins affect ATPase activity [33], drive microtubule motor functions [34], and participate in several biological processes, such as cilium assembly, cilium movement, and inner/outer dynein arm assembly [35]. Genomic alterations in DNAH genes have been reported in various malignant tumors, such as clear cell renal cell carcinoma, breast cancer, colorectal cancer, and gastric adenocarcinoma [31, 36,37,38]. However, no prior studies have explored their potential association with ICI response. In this regard, we are the first study to present evidence that DNAH mutations, particularly DNAH2, DNAH6, and DNAH9, are associated with both initial clinical response and long-term overall survival in melanoma patients receiving ICI therapies.

To investigate the role of DNAH alterations in immune modulation in the context of ICI response, we first analyzed the expression of several immune-related signatures, including the IFN-γ signature [27]. Although these signatures showed a slight increase in DNAH-altered tumors, IFN-γ alone demonstrated only a modest predictive value for ICI response in these tumors. However, DNAH-altered tumors exhibited higher TMB than DNAH-unaltered tumors, which is often associated with an increased number of neoantigens, enhancing immune recognition and activation. This heightened immune activity may be partially facilitated by IFN-γ signaling, which promotes immune cell recruitment and antigen presentation. Supporting this notion, immune cell infiltration analysis revealed higher fractions of monocytes and resting dendritic cells in DNAH-altered tumors. Monocytes, as key players in the innate immune response, can differentiate into macrophages and antigen-presenting dendritic cells, contributing to immune surveillance and inflammatory signaling [39, 40]. Meanwhile, dendritic cells, even in their resting state, serve as critical initiators of adaptive immunity by processing tumor-derived neoantigens and priming T cells [41]. The enrichment of these cells in DNAH-altered tumors, alongside elevated TMB, suggests a more immunogenic tumor microenvironment (TME), potentially contributing to improved ICI response. Furthermore, KEGG pathway analysis identified significant enrichment in immune-related pathways in DNAH-altered melanomas, including those associated with chemokine signaling, cytokine-cytokine receptor interaction, cell cycle, and p53 signaling. These findings were corroborated by a pan-cancer TCGA cohort, which demonstrated a consistent association between DNAH mutations and increased immune activity. Additionally, Zhu et al., reported a protective role of DNAH mutations in gastric cancer, correlating them with favorable prognosis and improved chemotherapy sensitivity, potentially due to higher TMB and increased genetic instability in DNAH-mutated tumors [31]. The association between upregulated cell cycle pathways and improved ICI response may be explained by two key mechanisms. First, increased cell cycle activity often leads to greater genomic instability and more frequent DNA replication errors, generating a higher TMB and increasing neoantigen presentation, thereby enhancing tumor recognition by the immune system. Second, tumors with upregulated cell cycle pathways tend to exhibit accelerated cell proliferation, fostering a more inflamed TME that promotes immune cell infiltration and proliferation-associated immunogenicity. Collectively, DNAH mutations emerge as a promising biomarker for predicting improved ICI response and prognosis, though further mechanistic studies are needed to fully elucidate their role in tumor progression and immune modulation.

At present, PD-L1 expression remains the most widely used predictive biomarker for ICI treatment. Although correlative, PD-L1 expression alone is insufficient to reliably predict a response following anti-PD-1 therapy. In our study, responders did exhibit higher PD-L1 expression compared to non-responders. However, neither PD-L1 expression nor the five publicly available immune-related signatures were effective in predicting clinical response, underscoring the need for additional transcriptional signatures to improve predictions of ICI response and prognosis. As expected, responsive tumors exhibited upregulated expression of genes related to active immune response, including lymphocyte proliferation and differentiation, T-cell activation, and cell–cell interactions. The candidate gene list was further narrowed down using correlation analysis to estimate each gene’s predictive and prognostic value following ICI treatment. This revealed a total of 84 genes with correlation coefficients of ≥  ± 0.2, of which 34 were associated with a favorable response and 50 with an unfavorable response. LASSO regression analyses further refined these to composite predictive (26-gene) and prognostic (13-gene) features, with 7 genes overlapping between the two groups. The composite predictive signature distinguished responders from those lacking response, while the prognostic feature effectively stratified overall survival, providing insights for personalized treatment and surveillance strategies for high-risk patients.

This study has several limitations. One notable limitation is the level of validation, which could be strengthened through immunohistochemistry or comparative analysis with external datasets involving patients treated with immunotherapy, whose pretreatment tumors were profiled by WES and WTS. Continued efforts are needed to translate these biomarkers into clinical practice to optimize immunotherapy-based treatment strategies across various cancer types. Additionally, subgroup analyses for special populations, such as those with mucosal and ocular melanomas who are known to be less responsive to ICIs, are needed. While we identified novel biomarkers potentially associated with ICI response and survival benefits, the features at different molecular levels were not well integrated, which limited the depth of our findings. Previous studies have highlighted the value of combining multi-omics data in uncovering immune response mechanisms and improving treatment prediction. Specifically, Zhou et al. demonstrated that SIRPA plays a complex and cell-type-dependent role in immunotherapy, exerting antagonistic effects in tumors and macrophages through both in vivo and in vitro analyses [42]. On the other hand, Mallardo et al. identified a combined protein and gene signature significantly associated with anti-PD-1 therapy outcomes, showing superior predictive value compared to single-omics approaches [43]. To address these limitations, future research incorporating multi-omics approaches, such as proteomics and single-cell RNA sequencing, could bridge the gap between DNA and RNA features and provide a more comprehensive understanding of immune response mechanisms.

Conclusion

In conclusion, the identification of combinatorial biomarkers through integrative genomic and transcriptomic profiling may hold significant promise for broader application across various malignancies, which may guide clinical trial design, advance monitoring ICI efficacy and prognosis surveillance, ultimately guiding the development of more personalized and effective anti-tumor immunotherapies.