HLA-I Evolutionary Divergence Confers Response to PD-1 Blockade plus Chemotherapy in Untreated Advanced Non–Small Cell Lung Cancer

Immune-checkpoint inhibitors (ICI) targeting PD-1 or PD-L1 (PD-1 blockade) plus cytotoxic chemotherapy have become the first-line standard of care for patients with advanced non–small cell lung cancer (NSCLC) without driver gene alterations (1). Several global or regional phase III trials reported that the objective response rate (ORR) of PD-1 blockade plus chemotherapy was approximately 45% to 60% in the intention-to-treat (ITT) population. Yet, not all of them can benefit from this regimen (2). There are three officially approved predictive biomarkers to guide the clinical application of PD-1 blockade monotherapy, including PD-L1 expression, tumor mutational burden (TMB), and microsatellite instability status (3, 4). However, they failed to show the feasibility to predict response to PD-1 blockade plus chemotherapy in advanced NSCLC. Therefore, the identification of robust predictive biomarkers for PD-1 blockade plus chemotherapy is urgently needed in current clinical practice.

Unlike PD-1 blockade monotherapy, the application of PD-1 blockade plus chemotherapy is mainly based on the undetermined hypothesis that some chemotherapeutic drugs may stimulate T cell–mediated immunity via increasing antigen release from tumor cell death as well as elimination of some immunosuppressive cells, such as regulatory T cells (Treg) and myeloid-derived suppressor cell (MDSC; refs. 5–7). As an essential component of antigen presentation, human leukocyte antigen class I (HLA-I) plays a pivotal role in T cell–mediated antitumor immunity (8). Theoretically, a more diverse HLA-I repertoire would lead to the presentation of broader antigens, increasing the odds of presenting more immunogenic antigens and increasing the likelihood of response to ICIs. In fact, several elegant studies have reported that the germline HED, a measure of sequence divergence between alleles of HLA-I genotype that is associated with the diversity of immunopeptidomes and intratumoral T-cell receptor clonality, could predict the outcomes of PD-1 blockade monotherapy in patients with various solid tumors including NSCLC, melanoma, and gastrointestinal cancer (9, 10). However, some studies did not support its predictive value for immunotherapy alone (11, 12). Under this circumstance, we supposed that a more diverse HLA-I repertoire would be more potent to predict response to PD-1 blockade plus chemotherapy in advanced NSCLC.

In this study, we integrated clinical and HLA-I genotype data of 427 patients with advanced NSCLC who received first-line PD-1 blockade plus chemotherapy from two previous phase III trials and investigated the predictive value of HED to PD-1 blockade plus chemotherapy. To delineate the immune profiles of tumors with distinct HED, we conducted single-cell RNA sequencing (scRNA-seq) of 58,977 cells collected from an independent cohort of 11 individuals with NSCLC.

Two previously published cohorts (CameL and CameL-sq) of patients from phase III trials were identified. Briefly, CameL is a randomized, open-label, multicenter, phase III trial to compare the efficacy and safety of camrelizumab plus chemotherapy with chemotherapy alone as first-line treatment for patients with non-squamous NSCLC without EGFR and ALK alteration (13). CameL-sq is a randomized, placebo-controlled, double-blind phase III trial to investigate the efficacy and safety of camrelizumab or placebo plus chemotherapy as first-line treatment for patients with advanced squamous NSCLC (14). Both studies met the primary endpoints and have predefined the biomarker analysis (13, 14). In the CameL study, patients with the following criteria were eligible: ages 18 to 70 years, histologically or cytologically confirmed stage IIIB–IV non-squamous NSCLC without EGFR or ALK alterations, Eastern Cooperative Oncology Group performance status (ECOG PS) of 0 or 1, no previous systemic chemotherapy, at least one measurable lesion per Response Evaluation Criteria in Solid Tumors version 1.1 (RECIST v1.1), and a life expectancy of ≥3 months. Patients were excluded if they had untreated central nervous system (CNS) metastases and corticosteroid use within two weeks before study treatment was recorded. In the CameL-sq study, eligible patients were: ages 18 to 75 years, pathologically confirmed stage IIIB–IV squamous NSCLC, no previously systemic therapy for metastatic disease, at least one measurable lesion per RECIST v1.1, ECOG PS of 0 or 1, archival (within 12 months from first study dose) or fresh tumor tissues for biomarker testing. Key exclusion criteria were: patients with active or symptomatic CNS metastases, EGFR or ALK alterations, history or presence of autoimmune disease or interstitial pneumonia, and use of immunosuppressants within two weeks before study treatment. Patients with adequate pretreatment fresh or archival formalin-fixed paraffin-embedded (FFPE) tissue and peripheral blood samples were eligible. This study was approved by independent ethics committees/institutional review boards at each study site and conducted according to the Declaration of Helsinki, Guidelines for Good Clinical Practice, and local laws and regulations of China.

Fresh or archival FFPE baseline samples from the CameL and CameL-sq studies were collected before the protocol-defined treatments. Fresh biopsy tissues were snap-frozen in liquid nitrogen within 30 minutes. Pretreatment blood samples (8–10 mL) were collected in ethylene diamine tetraacetic acid (EDTA)-coated tubes (RRID:SCR_013311, BD Biosciences) and centrifuged at 1,800 × g for 10 minutes within two hours of collection to separate white blood cells.

Fresh or FFPE tumor samples with ≥20% tumor cell contents were qualified and included. Genomic DNA from tumor tissues and peripheral blood cells were extracted with the QIAamp DNA FFPE Tissue Kit (RRID:SCR_008539, Qiagen) and DNeasy Blood and Tissue Kit (Qiagen) or GeneRead DNA FFPE Kit (Qiagen 180134) and Tguide S32 Magnetic Blood Genomic DNA Kit (RRID:SCR_023688, Tiangen) following the manufacturers’ protocols. DNA concentration was measured by Qubit dsDNA High-Sensitivity Assay Kit (RRID:SCR_008452, Thermo Fisher), whereas the quality of DNA was assessed by Agilent 2100 BioAnalyzer (RRID:SCR_019715, Agilent).

Following a previous report, we analyzed the LOH of HLA event of each tumor sample using the LOHHLA tool (RRID:SCR_023690, https://bitbucket.org/mcgranahanlab/lohhla/src/master/) with default parameter settings (19).

Thirty to 300 ng of genomic DNA was sheared to a length of approximately 200 base pairs (bp) by Covaris LE220 and library preparations were performed with KAPA Hyper Prep Kit (KAPA Biosystems). Libraries were quantified with quantitative PCR using the KAPA Library Quantification kit (KAPA Biosystems), and the size was determined using the Bioanalyzer 2100 (Agilent Technologies). Sequencing was performed on the Illumina HiSeq4000 platform using PE150 sequencing chemistry (Illumina).

Base calling was performed on bcl2fastq (RRID:SCR_015058) V.2.16.0.10 (20) to generate sequence reads in the FASTQ format. After removing low-quality reads by Trimmomatic (v0.36; RRID:SCR_011848; ref. 21), clean reads were aligned to the human reference genome (hg19, NCBI Build 37.5; RRID:SCR_012917) with the Burrows–Wheeler Aligner (BWA, version 0.7.17; RRID:SCR_010910). Then, the Picard toolkit (version 2.23.0, RRID:SCR_006525; ref. 22) was used to convert SAM files to compressed BAM files, which were then sorted according to chromosome coordinates. The Genome Analysis Toolkit (GATK, RRID:SCR_001876; ref. 23) was used for realignment. Single-nucleotide variants (SNV) and small insertions/deletions (InDel) were called via MuTect2 (RRID:SCR_000559; ref. 24) with tumor–normal mode. To avoid false-positive results, SNVs and InDels that appeared on the blacklist (including sequence-specific errors, repeat regions, segmental duplications, and lowly mappable regions recorded in ENCODE; RRID:SCR_006793) were removed. After annotation by ANNOVAR (RRID:SCR_012821; ref. 25), we filtered out variants in either introns or synonymous mutations. Somatic mutations with following criteria were used for analysis: (i) the sequencing depth was more than 100×; (ii) the variant allele frequency threshold of SNV was 4% and that of InDels was 5%. Furthermore, variants with a minor allele frequency ≥1% in the Exome Aggregation Consortium (ExAC, RRID:SCR_004068) and Genome Aggregation Database (gnomAD, RRID:SCR_014964) were removed. Mutations not recorded in the COSMIC (RRID:SCR_002260) database were filtered out. Somatic copy-number variations were identified using the CNVkit (v0.9.5, RRID:SCR_021917).

As previously described (26), TMB was calculated by integrating the total number of somatic, base substitutions, coding, and InDel mutations per megabase of the genome examined. We first quantified the number of somatic nonsynonymous SNVs. Then, the value was extrapolated to the whole-exome using a validated algorithm (27). Only the regions with sequencing depth >100× after deduplication were used for TMB calculation. Germline genomic alterations in the Single-Nucleotide Polymorphism database or occurring with ≥2 counts in the ExAC database were not recorded.

Similar to our previous study (26), we calculated the neoantigen burden by using the following algorithm. First, we called the HLA-I alleles from the matched normal exome sequencing data using HLA-HD (v1.2.0.1, RRID:SCR_022285; ref. 28). Neoepitope presentation was then predicted for tumor-specific peptides of length 9–11 using the eluted-ligand mode of NetMHCpan-4.0 (RRID:SCR_018182; ref. 29). The following criteria were used to select the potential neoantigens for subsequent analysis: (i) derived from tumor-specific genomic alterations (including missense, inframe indels, frameshift, and fusions); (ii) high predicted affinity to HLA-I alleles [half maximal inhibitory concentration (IC₅₀) < 500 nmol/L] with k-mer of 9–11 length; (iii) fold change >10 comparing to wild-type binding affinity. IEDB-recommended model (RRID:SCR_006604) was conducted to predict the HLA binding affinity using all variant-containing 9–11 mer for HLA-A/B/C binding estimations. HLA typing for patients was performed in silico using HLA-ATHLATES (RRID:SCR_023689; ref. 30) per the recommended algorithm.

Our previously published raw scRNA-seq data of patients with NSCLC were collected and integrated for comparison of immune cell landscape between tumors with high and low mean HED (31, 32). Their matched para-carcinoma normal tissues were used for HED calculation as previously described. Briefly, filtering of low-quality cells, data normalization, sample integration, clustering, and identification of cell types were performed using Seurat (version 4.1.1; ref. 33): (i) scRNA-seq data sets from 11 patients with NSCLC were loaded into R (version 4.1.3) using the CreateSeuratObject function of Seurat (RRID:SCR_016341) and merged together. (ii) Low-quality cells where >20% of transcripts derived from the mitochondria and either <500 expressed genes or >5,000 were filtered out. (iii) The NormalizeData function was performed to normalize raw counts. The top 2,000 highly variable genes were used for the following principal component analysis (PCA). (iv) The IntegrateData function was applied to remove batch effects among 11 patients. (v) The RunPCA and RunUMAP functions were used to reduce dimension. Subsequently, clustering cells by executing the FindNeighbors and FindClusters functions.

The Wilcoxon test was used to identify changes in the proportion of each cell type between tumors with high and low mean HED in total cells or immune cells in our cohorts. Significantly different cell types were then validated by using scRNA-seq data of a published cohort of 47 tumor samples from 36 patients with NSCLC who received PD-1–based treatments (34). A two-tailed Student t test was performed to calculate the significant differences of immune cells among pretreatment, posttreatment nonresponsive, and posttreatment responsive tumors.

DEGs between tumors with high and low mean HED in specific cell types were defined by performing the FindMarkers procedure of Seurat, remaining top 100 DEGs with both |Log₂FC| > 0.5 and adjusted P < 0.01.

GO enrichments were performed by using the enrichGO function in the clusterProfiler package (version 3.12.0; refs. 35, 36) based on the highly expressed DEGs for the HED^high arm in specific cell types (P value cutoff = 0.01, q-value cutoff = 0.05). The AddModuleScore function of Seurat was applied to calculate the activity of each GO term in each cell.

The Kruskal–Wallis test was used for comparisons of HED distributions across individual HLA-I loci. Of the CameL and CameL-sq cohorts, cutoffs for high mean HED were determined using the top quartile, and low mean HED were defined as values less than the top quartile. For TCGA cohort, the absolute cutoff value of the CameL cohort was applied to the LUAD cohort, and the absolute cutoff value of CameL-sq was used for the LUSC cohort. High TMB was defined as values greater than or equal to 10 Muts/Mb, and positive PD-L1 expression was determined using the tumor proportional score greater than or equal to 1%. When combining mean HED and PD-L1 expression levels, patients with both high mean HED and positive PD-L1 expression were defined as the “Both high” group, whereas the other patients belonged to the “Others” group. When combining mean HED and TMB, patients were defined as the “Both high” group if their mean HED and TMB were both greater than or equal to the cutoff points for mean HED and TMB, and the others were defined as the “Others” group.

Survival curves were performed using Kaplan–Meier analysis. P values were calculated using the log-rank test. Hazard ratios (HR) and corresponding 95% confidence intervals (CI) were analyzed using stratified Cox proportional hazards model. For multivariable analyses, P values, HRs, and corresponding 95% CIs were analyzed using the Cox proportional hazards model and were visualized by “function ggforest.” The Fisher exact test was introduced to analyze the difference in response rate between high and low mean HED groups. The Spearman coefficients method was used for the correlation analysis. Two-side P < 0.05 was considered statistically significant. All statistical analyses were performed with R version 4.0.5 (RRID:SCR_001905).

The data that support the findings of this study are available from the corresponding author upon reasonable request. All requests for raw data will be reviewed by the leading clinical site, Shanghai Pulmonary Hospital, Tongji University School of Medicine, and the study sponsor, Jiangsu Hengrui Pharmaceuticals, to check whether the request is subject to any intellectual property or confidentiality obligations. A signed data access agreement with the sponsor is required before accessing shared data. Source data are provided in this paper. All the relevant raw data were uploaded to the National Genomics Data Center. The accession numbers are HRA001180 and HRA004337 (https://ngdc.cncb.ac.cn/gsa-human/). We also listed the research resource identifier in Supplementary Table S4.

To investigate the predictive value of HED to PD-1 blockade plus chemotherapy, we conducted integrated analyses using the CameL study (n = 174; ref. 13), and the CameL-sq study (n = 253; ref. 14; Fig. 1A; Supplementary Fig. S1 and Table 1). The clinical data of two studies used for biomarker association analyses were based on the updated analysis of both trials. Demographic and baseline parameters in the biomarker-evaluable population were generally analogous to those in the ITT population (Supplementary Tables S1 and S2). Treatment outcomes were also representative of the ITT population in both cohorts.

We first delineated the landscape of HEDs at HLA-A, HLA-B, and HLA-C in the CameL study. Hierarchical clustering of HEDs showed that there were different clusters of high and low divergence between alleles (Supplementary Fig. S2), in line with previous reports (9, 10) and known interrelationships of the three loci. Mean HED distributions were similar between patients of the CameL and TCGA LUAD cohorts (Fig. 1B), and patients from PD-1 blockade plus chemotherapy and chemotherapy groups (Fig. 1C). Comparison of distribution patterns of HEDs for each HLA-A, HLA-B, and HLA-C heterozygous genotype showed that HLA-C pairwise divergences were significantly lower to HLA-A and HLA-B pairwise divergences (Fig. 1D), in accordance with previous studies reported HLA-C as the most recently evolved gene (37). There was no obvious difference in pairwise divergences between HLA-A and HLA-B alleles (Fig. 1D).

Similar to the previous definition (10, 38), we defined “HED^high” as the mean HED greater than or equal to the top quartile, and “HED^low” as the mean HED less than the top quartile. The demographic and baseline parameters in the HED-high group were generally analogous with those in the HED-low group in both treatment cohorts. First, we excluded the mean HED as a generally prognostic factor in both LUAD and LUSC by using TCGA data sets (Supplementary Fig. S3). Then, we surveyed the associations between mean HED and outcomes of PD-1 blockade plus chemotherapy in all included patients. The results showed that HED^high was correlated with better OS and PFS than HED^low in patients who received PD-1 blockade plus chemotherapy, but not chemotherapy, in both cohorts (Supplementary Fig. S4). Meanwhile, HED^high also associated with numerically higher ORR to PD-1 blockade plus chemotherapy in the CameL study (P = 0.113; Supplementary Fig. S4).

Given the impact of HLA-I heterozygosity on the diversity of immunopeptidomes, we hypothesized that the predictive performance of mean HED could be improved in patients with fully heterozygous HLA-I genotypes. In fully heterozygous patients of the CameL cohort (n = 143), patients with HED^high had significantly longer OS (HR = 0.40, 95%CI = 0.19–0.85; P = 0.014) and PFS (HR = 0.47, 95% CI = 0.26–0.86; P = 0.012) than those with HED^low in PD-1 blockade plus chemotherapy arm (Fig. 2A and B). In the HED^high group, patients treated with PD-1 blockade plus chemotherapy had significantly better OS (HR = 0.47, 95% CI = 0.20–1.09; P = 0.071) and PFS (HR = 0.36, 95% CI = 0.18–0.72; P = 0.003) than those treated with chemotherapy, whereas no significant difference on OS and PFS was observed in patients with HED^low (Fig. 2A and B). The associations between HED^high and treatment outcomes were absent in the chemotherapy arm (OS: HR = 1.00, 95% CI = 0.53–1.86; P = 0.987; PFS: HR = 0.80, 95% CI = 0.46–1.39; P = 0.431). Moreover, patients with HED^high had a markedly better ORR than those with HED^low in the PD-1 blockade plus chemotherapy arm (ORR: 81.8% vs. 53.2%; P = 0.032; Fig. 2C), whereas it was similar in the chemotherapy arm (ORR: 27.3% vs. 47.8%; P = 0.123; Fig. 2D). The predictive value of HED^high for PD-1 blockade plus chemotherapy held in multivariate analysis (OS: HR = 0.39, 95% CI = 0.18–0.85; P = 0.018; PFS: HR = 0.47, 95% CI = 0.26–0.88; P = 0.0.18; Fig. 2E) when adjusted for PD-L1 expression and TMB level, confirming that the effect of HED^high on survival benefit may be specific and independent of PD-L1 expression and TMB.

We next validate the predictive value of HED^high in fully heterozygous patients of the CameL-sq study (n = 196; ref. 14). Patients with HED^high had dramatically prolonged OS (HR = 0.38, 95% CI = 0.21–0.71; P = 0.002) and PFS (HR = 0.49, 95% CI = 0.29–0.81; P = 0.005) than those with HED^low in PD-1 blockade plus chemotherapy arm (Fig. 2F and G). Patients with HED^high received PD-1 blockade plus chemotherapy had markedly prolonged OS (HR = 0.40, 95% CI = 0.20–0.80; P = 0.007) and PFS (HR = 0.08, 95% CI = 0.04–0.18; P < 0.001) than those treated with chemotherapy, whereas it was similar between PD-1 blockade plus chemotherapy and chemotherapy groups in patients with HED^low (Fig. 2F and G). Superior ORR was observed in patients with HED^high than in those with HED^low in PD-1 blockade plus chemotherapy arm (89.2% vs. 62.3%; P = 0.007; Fig. 2H) whereas it was similar in the chemotherapy arm (58.6% vs. 47.1%; P = 0.378; Fig. 2I). Multivariate analysis also showed the significantly predictive performance of HED^high for PD-1 blockade plus chemotherapy (Fig. 2J), further supporting that HED^high is a reliable predictive biomarker to PD-1 blockade plus chemotherapy in untreated advanced NSCLC, especially those with fully heterozygous HLA-I genotypes.

Previously, Lu and colleagues reported that HED of HLA-A, HLA-B, or HLA-C had a distinct impact on predicting ICI-treated outcomes in gastrointestinal cancer, and HLA-B HED was correlated with improved clinical benefit after ICI treatment (9). Therefore, we evaluated the effect of HED of HLA-A, HLA-B, or HLA-C on treatment outcomes. As shown in Supplementary Figs. S5 and S6, each of HLA-A, HLA-B, or HLA-C HED cannot independently and perfectly predict the treatment response and outcomes of PD-1 blockade plus chemotherapy in both cohorts.

Considering the predictive robustness of mean HED, we then assessed the effect of distinct cutoffs of mean HED on treatment outcome prediction. Generally, a negative relationship was observed between mean HED and HR of PFS, suggesting that an increase in mean HED corresponds to improved survival benefit (Supplementary Fig. S7A and S7B). Interestingly, although several cutoffs of mean HED showed statistical significance to predict treatment outcomes, the top quartile seemed to be the optimal cutoff to stratify patients into groups with different clinical outcomes in both cohorts. In addition, previous studies revealed that HLA-I LOH is one of the key immune-evade mechanisms and negative predictors of outcomes after ICI treatment in NSCLC (19, 39). We also explored the associations between HLA-I LOH and response to PD-1 blockade plus chemotherapy. The results showed that HLA-I LOH cannot predict treatment response and outcomes of PD-1 blockade plus chemotherapy in both the cohorts (Supplementary Fig. S7C and S7F).

Taken together, these data highlight the importance of mean HED and suggest that HED^high may portend survival benefit in advanced NSCLC patients treated with first-line PD-1 blockade plus chemotherapy.

As we mentioned above, the increasing mean HED was correlated with longer PFS, suggesting a “dose-dependent manner” between mean HED and treatment benefit. To confirm this association, we utilized the STEPP approach to investigate the 2-year PFS probability along the mean HED range according to a previous publication (38). In both cohorts, we observed the nearly linear association between mean HED and improved 2-year PFS probability in the PD-1 blockade plus chemotherapy arm (Fig. 3A–B). However, higher mean HED was not associated with improved outcome in the chemotherapy arm. Notably, the threshold of these benefits was slightly different between CameL and CameL-sq cohorts. Collectively, these data suggest that HED^high could predict long-term treatment benefits of PD-1 blockade plus chemotherapy in advanced NSCLC.

To evaluate the combined effect of mean HED and other known predictive biomarkers for immunotherapy, including PD-L1 expression or TMB on response to PD-1 blockade plus chemotherapy, we first examined the correlation between mean HED and PD-L1 expression or TMB. Mean HED did not correlate with PD-L1 expression or TMB (Fig. 4A and B) in the CameL cohort, whereas mean HED showed mild correlation with PD-L1 expression but not TMB in the CameL-sq cohort (Supplementary Fig. S8A and S8B). No significant correlation was detected between mean HED and neoantigen burden (Fig. 4C). We also checked the associations between mean HED at each HLA-I locus and PD-L1 expression, TMB, or neoantigen burden, and found no correlation among them (Supplementary Fig. S9). When stratifying patients with fully heterozygous HLA-I genotypes by mean HED and TMB level, we found that the joint utility of both HED^high and high TMB showed an improved effect on the survival benefit of PD-1 blockade plus chemotherapy in the CameL cohort but not in the CameL-sq cohort, and no improved effect on treatment response in both cohorts (Supplementary Fig. S10). Of note, HED^high could predict the survival benefit in both high and low TMB groups (Supplementary Fig. S11), indicating that the predictive value of HED^high for PD-1 blockade plus chemotherapy was greatly irrespective of the TMB level.

To evaluate the impact of HED^high combined with PD-L1 expression on survival benefit and treatment response, we then stratified patients with fully heterozygous HLA-I genotypes by mean HED and PD-L1 expression. The results showed that the combination of HED^high and positive PD-L1 expression showed better predictive performance for both OS (HR = 0.29, 95% CI = 0.12–0.69; P = 0.003; Fig. 4D) and PFS (HR = 0.41, 95% CI = 0.22–0.80; P = 0.006; Fig. 4E) in the PD-1 blockade plus chemotherapy arm, as evidenced by reduced HR and corresponding P values. The combined effect was absent in the chemotherapy arm (Fig. 4D and E). Furthermore, ORR was significantly better in patients with HED^high and positive PD-L1 expression in the PD-1 blockade plus chemotherapy arm (ORR: 84.2% vs. 52.1%; P = 0.025; Fig. 4F) whereas it was similar in the chemotherapy arm (ORR: 37.5% vs. 45.7%; P = 0.771; Fig. 4G). The combined effect of mean HED and PD-L1 expression on survival benefit after PD-1 blockade plus chemotherapy was further confirmed in the CameL-sq cohort, which also showed reduced HR and corresponding P values of both OS (HR = 0.24, 95% CI = 0.10–0.62; P = 0.001; Supplementary Fig. S8C) and PFS (HR = 0.34, 95% CI = 0.18–0.65; P = 0.001; Supplementary Fig. S8D). Analogously, ORR was dramatically higher in patients with HED^high and positive PD-L1 expression in PD-1 blockade plus chemotherapy arm (ORR: 89.2% vs. 62.2%; P = 0.007; Supplementary Fig. S8E) and similar in the chemotherapy arm (ORR: 58.6% vs. 47.1%; P = 0.378; Supplementary Fig. S8F).

To determine the immune profiles of tumors with distinct HEDs, we analyzed a large scRNA-seq data of 58,977 cells collected from 11 individuals with NSCLC (Supplementary Table S3). Dimensionality reduction with t-distributed stochastic neighborhood embedding (t-SNE) and graph-based clustering identified eight major cell subsets, including T, B, natural killer like T, myeloid, fibroblast, endothelial, mast, and cancer cells (Fig. 5A; Supplementary Fig. S12A). We found an increased number of T and B cells in the HED^high group, whereas the proportion of myeloid and cancer cells was significantly lower in the HED^high group than in the HED^low group (Fig. 5B). We further identified 14 immune cell subsets using graph-based clustering (Fig. 5C; Supplementary Figs. S12B and S13A and S13B). We observed that patients with HED^high had a significantly higher proportion of CD4⁺ central memory T (Tcm), CD4⁺ Th1-like, and follicular B cells (Fig. 5D). Pathway enrichment analyses showed a strong enrichment of antigen presentation and immune response in CD4⁺ Tcm, CD4⁺ Th1-like, and follicular B cells from the HED^high group (Fig. 5E–G and I). However, several immunosuppressive cells, including classic dendritic cell type 2 (cDC2), monocyte, macrophage, neutrophil, and mast cells, was significantly lower in the HED^high group (Fig. 5D and H). These significantly different cells between HED^high and HED^low groups had distinct transcriptomic features and enriched pathways (Supplementary Fig. S13E–S13H).

Despite the huge success observed in clinical trials evaluating first-line PD-1 blockade plus chemotherapy in advanced NSCLC, not all of them can benefit from this regimen. This divergence in terms of treatment benefit highlights a critical need to identify robust predictive biomarkers. Given the potential mechanism of the synergistic antitumor effect of PD-1 blockade plus chemotherapy that builds upon the increased tumor-related antigen release by some chemotherapeutic drugs (5), this study investigated the predictive significance of HED for PD-1 blockade plus chemotherapy in untreated advanced NSCLC from two phase III trials. Our results showed that HED^high was associated with significantly better treatment response, survival and 2-year PFS probability in NSCLC patients treated with first-line PD-1 blockade plus chemotherapy, especially in those with fully heterozygous HLA-I genotypes, but not in the chemotherapy group. We also found that a combination of mean HED and PD-L1 expression showed better predictive performance. Using scRNA-seq of untreated NSCLC, we found that patients with HED^high were associated with improved antigen presentation and antitumor immunity, further supporting HED^high as a predictive biomarker in this setting.

Previously, several studies investigated the predictive value of HED for PD-1 blockade monotherapy. In 2019, Chowell and colleagues reported the first evidence that patients with HED^high received PD-1 blockade monotherapy showed a longer OS in melanoma and NSCLC (10). They further found that HED^high was associated with both improved clinical benefit and durability of response in advanced renal cell carcinoma patients treated with PD-1 blockade plus lenvatinib (38). Then, Lu and colleagues reported that high HLA-B HED was correlated with a better OS compared with the low HLA-B HED subgroup in gastrointestinal cancer (9). Nevertheless, controversial findings were reported in two studies with large samples. The first group of investigators integrated patients across 17 pembrolizumab clinical trials and reported that HLA-I diversity including HED was not associated with response to pembrolizumab (11). The second one performed a meta-analysis investigating the association between germline HLA-I heterozygosity and HED and outcomes with PD-1 blockade among various solid tumors and also reported no association with response (12). This inconsistency might be due to the different patient distribution. For example, nearly all of the patients are of European ancestry and a majority of them with PD-L1–positive expression in the pembrolizumab clinical trials (11). Moreover, most of them received pembrolizumab monotherapy in the later-line setting, which could not represent the whole population. In addition, they included a mixture of patients from distinct trials across a long-time interval, which would have an impact on the homogeneous estimation of survival benefit. Nowadays, PD-1 blockade plus chemotherapy have been established as a new standard of care in various solid tumors including NSCLC. Considering chemotherapy could enhance immunogenicity via increasing tumor antigen release, the predictive value of HED needs to be further reevaluated for PD-1 blockade plus chemotherapy in advanced NSCLC.

To our best knowledge, the current study first investigated the predictive significance of HED for PD-1 blockade plus chemotherapy in patients with untreated advanced NSCLC from two phase III trials. Our results showed that HED^high was associated with significantly better response and survival benefit in NSCLC patients treated with first-line PD-1 blockade plus chemotherapy in both trials, especially in patients with fully heterozygous HLA-I genotypes. We also found that HED^high was associated higher landmark time probability of 2-year PFS. However, it did not affect the treatment outcome in the chemotherapy group. Multivariate analysis demonstrated that HED^high was a strong determinant of efficacy for PD-1 blockade plus chemotherapy. Currently, both the approved companion diagnostics of PD-L1 expression and TMB cannot be used for predicting response to PD-1 blockade plus chemotherapy in advanced NSCLC. Our study was the first one to reveal HED as the potential predictor for first-line PD–1 blockade plus chemotherapy in advanced NSCLC, though the robustness of its predictive value remains to be investigated. More importantly, our finding highlights the necessity of identifying distinct predictive biomarkers regarding different immunotherapeutic regimens.

The effect of the combinatorial utility of mean HED and TMB on survival benefits from ICI treatment has been demonstrated to be superior to either alone in melanoma and gastrointestinal cancer (10, 40), suggesting a rational combination of different potential biomarkers could improve the predictive performance. We also explored the combined effect of mean HED and PD-L1 expression or TMB on response to PD-1 blockade plus chemotherapy. The results showed that the joint utility of HED^high and PD-L1 expression showed better performance than either alone in predicting treatment benefits from this combination. Patients with both high HED and positive PD-L1 expression are most likely to benefit from PD-1 blockade plus chemotherapy, underscoring the rationale of the combinatorial utility of mean HED and PD-L1 expression in advanced NSCLC patients treated with first-line PD-1 blockade plus chemotherapy. Herein, both HED and PD-L1 expression should be taken into consideration in NSCLC in future clinical trials.

Having noticed the potential predictive value of mean HED, we are eager to understand the immune profiles of tumors with HED^high versus HED^low. We adopted the scRNA-seq analysis of pretreated tumors with available HLA-I genotypes’ data and found that tumors with HED^high had markedly higher infiltration of CD4⁺ Tcm, CD4⁺ Th1-like, and follicular B cells, whereas cDC2, monocyte, macrophage, neutrophil, and mast cells were significantly lower in the HED^high group. CD4⁺ Tcm, CD4⁺ Th1-like, and follicular B cells were previously found to be associated with improved T cell–mediated immunity (41, 42), whereas cDC2, monocyte, macrophage, neutrophil, and mast cells were immunosuppressive in most cases (43–45). These findings indicate that tumors with HED^high would possess an immunosupportive TME phenotype. Interestingly, our further analysis revealed strong enrichment of antigen presentation and immune response in CD4⁺ Tcm, CD4⁺ Th1-like, and follicular B cells from the HED^high group. Importantly, CD4⁺ Tcm and CD4⁺ Th1-like cells were found to be significantly increased in responsive tumors after PD-1 blockade plus chemotherapy in a recent study (Supplementary Fig. S13C and S13D; ref. 34). Collectively, these data revealed that tumors with HED^high were associated with improved antigen presentation and T cell–mediated antitumor immunity, supporting it as a potential predictor for PD-1 blockade plus chemotherapy in NSCLC.

There were several limitations of this work. First, despite we analyzed the largest cohort of NSCLC patients treated with first-line PD-1 blockade plus chemotherapy to identify the potential biomarker hitherto, the current study was a retrospective analysis with hypothesis generating in nature. Hence, we should cautiously interpret these findings until they are validated in prospective, standardized biomarker studies. Second, all of the included patients were Asians. Whether HED^high was also a reliable predictor of PD-1 blockade plus chemotherapy in non-Asian ancestry patients with distinct HLA-I genotypes and immunophenotypes remained unknown. Third, we utilized the top quartile as the cutoff to define HED^high versus HED^low groups. Although the cutoff is popular in prior publications and helpful for us to clarify the relevant investigations due to its briefness, it may not be optimal. Fourth, we only focused on the genetic variables to predict treatment benefit and did not evaluate the tumor immune microenvironment features. The integration of multiple-dimensional information would be valuable for the identification of robust biomarkers. Fifth, we did not observe the robust predictive or prognostic value of TMB in this cohort by using 10 Muts/Mb as the cutoff. Although the definition and cutoff of TMB are popular, the use of 10 Muts/Mb cutoff was mostly applied to commercially targeted gene panels, and the data of this study were generated by whole-exome sequencing. Considering the great discrepancy of the test platform, process, analysis algorithm, and definitions, we still need larger populations to validate the predictive and/or prognostic value of TMB in NSCLC received first-line PD-1 blockade plus chemotherapy. Finally, although we performed mIF and scRNA-seq to depict the immune profiles of tumor with distinct HED, we did not unravel the detailed biological explanations. The scRNA-seq data were generated mainly from patients with early-stage NSCLC without driver gene status, PD-L1 expression, and TMB level, so there may be significant confounders in terms of late-stage lung cancers with HED high and immune cell landscapes. Future clinical and translational research in this field would provide mechanistic insights into how HED could improve the T cell–mediated tumor control and predict the survival benefit of NSCLC patients who received first-line PD-1 blockade plus chemotherapy.

In conclusion, HED^high represents a potential biomarker to predict the response and survival outcomes of patients with untreated advanced NSCLC who received PD-1 blockade plus chemotherapy, especially in those with fully heterozygous HLA-I genotypes. Further investigation in prospective studies would be helpful to optimize PD-1 blockade plus chemotherapy in patients with untreated advanced NSCLC who are more likely to derive benefit. Considering HED can be analyzed from normal tissue (e.g., peripheral blood cells) DNA sequencing, it would be more convenient for daily clinical practice when ample tissue samples were difficult to obtain.

C. Zhou reports honoraria as a speaker from Roche, Lily China, Boehringer Ingelheim, Merck, Hengrui, Qilu, Sanofi, Merck Sharp & Dohme, Innovent Biologics, C-Stone, Luye Pharma, TopAlliance Biosciences, and Amoy Diagnostics and advisor fees from Innovent Biologics, Hengrui, Qilu, and TopAlliance Biosciences. No disclosures were reported by the other authors.

T. Jiang: Conceptualization, resources, data curation, software, formal analysis, supervision, funding acquisition, validation, investigation, visualization, methodology, writing–original draft, project administration, writing–review and editing. Q. Jin: Resources, data curation, software, formal analysis, funding acquisition, validation, visualization, methodology. J. Wang: Resources, data curation, software, formal analysis, supervision, visualization, methodology. F. Wu: Resources, data curation, software, formal analysis, supervision, investigation, visualization. J. Chen: Resources, data curation, software, formal analysis, supervision. G. Chen: Resources, data curation, software, formal analysis. Y. Huang: Resources, data curation, software, formal analysis. J. Chen: Resources, data curation, software, formal analysis. Y. Cheng: Resources, data curation, software, formal analysis. Q. Wang: Resources, data curation, software, formal analysis. Y. Pan: Resources, data curation, software, formal analysis. J. Zhou: Resources, data curation, software, formal analysis. J. Shi: Resources, data curation, software, formal analysis. X. Xu: Resources, data curation, software, formal analysis. L. Lin: Resources, data curation, software, formal analysis. W. Zhang: Resources, data curation, software, formal analysis. Y. Zhang: Resources, data curation, software, formal analysis. Y. Liu: Resources, data curation, software, formal analysis. Y. Fang: Resources, data curation, software, formal analysis. J. Feng: Resources, data curation, software, formal analysis. Z. Wang: Resources, data curation, software, formal analysis. S. Hu: Resources, data curation, software, formal analysis. J. Fang: Resources, data curation, software, formal analysis. Y. Shu: Resources, data curation, software, formal analysis. J. Cui: Resources, data curation, software, formal analysis. Y. Hu: Resources, data curation, software, formal analysis. W. Yao: Resources, data curation, software, formal analysis. X. Li: Resources, data curation, software, formal analysis. X. Lin: Resources, data curation, software, formal analysis. R. Wang: Resources, data curation, software, formal analysis. Y. Wang: Resources, data curation, software, formal analysis. W. Shi: Resources, data curation, software, formal analysis. G. Feng: Resources, data curation, software, formal analysis. J. Ni: Resources, data curation, software, formal analysis. B. Mao: Resources, data curation, software, formal analysis. D. Ren: Resources, data curation, software, formal analysis. H. Sun: Resources, data curation, software, formal analysis. H. Zhang: Resources, data curation, software, formal analysis. L. Chen: Conceptualization, resources, data curation, software, formal analysis, validation, visualization, methodology, project administration. C. Zhou: Conceptualization, resources, data curation, software, formal analysis, funding acquisition, validation, project administration, writing–review and editing. S. Ren: Conceptualization, resources, data curation, software, formal analysis, supervision, funding acquisition, validation, investigation, visualization, writing–original draft, project administration, writing–review and editing.

We are grateful to all patients and their families and all staff at the study centers. This study was also supported in part by grants from the National Natural Science Foundation of China (No. 82102859, 12126605, T2341007, 12131020, 31930022, 12026608), JST Moonshot R&D Grant (No. JPMJMS2021), Shanghai Public Health Committee Foundation (2020CXJQ02), Shanghai Rising-Star Program (23QA1408000), Oncology development incentive program of Shanghai Pulmonary Hospital, Cultivation project for National Natural Science Foundation of China Youth Fund of Shanghai Pulmonary Hospital (fkzr2306), Innovation Group of Shanghai Pulmonary Hospital (FKCX1903), Shanghai Hospital Development Center Three-year Action Plan to Promote Clinical Skills and Clinical Innovation in Municipal Hospitals (SHDC2020CR1036B), and the Project of Shanghai Municipal Health Commission: Establishment, Promotion and Application of Multidisciplinary Collaborative Diagnosis and Treatment System for Pulmonary Non-infectious Diseases.

The publication costs of this article were defrayed in part by the payment of publication fees. Therefore, and solely to indicate this fact, this article is hereby marked “advertisement” in accordance with 18 USC section 1734.