Multi-omics Analysis Reveals Immune Features Associated with Immunotherapy Benefit in Patients with Squamous Cell Lung Cancer from Phase III Lung-MAP S1400I Trial

Immune checkpoint inhibitors (ICI) targeting programmed cell death protein 1 (PD-1; e.g., nivolumab, pembrolizumab, cemiplimab) or its ligand PD-L1 (e.g., atezolizumab) have become pillars of treatment in both frontline and salvage settings for patients with advanced non–small cell lung cancer (NSCLC; refs. 1–4). In addition, recent efforts have led to multiple approved frontline regimens incorporating chemotherapy and other ICIs with anti-PD-1/PD-L1 antibodies (5–8). However, in the salvage setting, anti-PD-1/PD-L1 monotherapy remains the treatment of choice for ICI-naïve advanced-stage NSCLC (9, 10).

Ipilimumab is an ICI targeting CTL-associated protein 4 (CTLA-4). Its dual inhibition with PD-1/PD-L1 may have synergistic effects on the anticancer immune response, given the complementary functions of these two pathways. The combination of nivolumab with ipilimumab (nivo+ipi) was demonstrated to have superior efficacy than nivolumab alone in patients with advanced melanoma (11, 12). For patients with metastatic NSCLC, ipilimumab plus nivolumab has been approved by the FDA in the frontline setting with or without concurrent chemotherapy (7, 8, 13, 14). In the salvage setting, a recent phase III study, S1400I, evaluated the efficacy of nivo+ipi versus nivolumab monotherapy (nivo) in patients without previous ICI treatment for squamous NSCLC (SqNSCLC; ref. 15). The study did not show that ipilimumab plus nivolumab improved clinical outcomes. However, progression-free survival (PFS) and overall survival (OS) curves separated during later follow-up, suggesting that a subset of patients may benefit from combination treatment with ipilimumab and nivolumab.

Understanding the mechanisms underlying response and resistance to ICIs and establishing predictive molecular and immune features to identify patients who will benefit the most from ICI therapy remain unmet clinical needs. High PD-L1 expression is associated with improved outcomes in patients receiving ICI monotherapy (1, 8). However, the geographical heterogeneity of PD-L1 expression between primary tumors and metastatic sites and even between different regions within the same tumors—as well as the potential dynamic changes in PD-L1 expression over time—have raised questions about its reliability as a predictive biomarker (16, 17). Although tumor mutational burden (TMB) has been approved as a predictive marker for anti-PD-1/PD-L1 treatment for melanoma and NSCLC, and several other cancer types (18), one study found no correlation between TMB or PD-L1 with anti-PD-1 plus anti-CTLA-4 therapy in patients with NSCLC (19). Furthermore, the predictive value of PD-L1 and TMB becomes less clear when chemotherapy is added. These findings underscore the complexity of molecular determinates of the tumor immune microenvironment and response to ICIs.

In this study, we sought to elucidate the immune and molecular mechanisms that affect benefit from ICIs in patients with advanced SqNSCLC. Toward this end, we integrated immune and multi-omics profiling platforms supported by Cancer Immune Monitoring and Analysis Centers (CIMAC) in the current study. Specifically, we performed multiplex immunofluorescence (mIF), gene expression profiling (ncounter PanCancer Immune Profiling Panel), whole-exome sequencing (WES), and Olink proteomics on tissue and blood specimens from the S1400I trial to identify molecular or immune factors associated with better prognoses in patients treated with anti-PD-1 monotherapy versus anti-PD-1/CTLA-4 dual combination.

Lung-MAP (S1400I, NCT02785952) was a multicenter, open-label, phase III randomized clinical trial. The substudy Lung-MAP-I (S1400I) was conducted from December 18, 2015, to April 23, 2018, through the National Clinical Trials Network and led by the SWOG Cancer Research Network. The study was conducted in accordance with the Declaration of Helsinki and the Lung-MAP design has been described previously (15). Briefly, the trial compared nivo+ipi with nivo in patients with chemotherapy-pretreated, immunotherapy-naïve, advanced SqNSCLC. Two hundred fifty-two patients were randomly assigned to receive nivo+ipi (n = 125) or nivo (n = 127). The clinical efficacy endpoints were OS, PFS, duration of response, and best objective response by RECIST 1.1. Each site required approval by the U.S. NCI central Institutional Review Board or approval by their local Institutional Review Board. Written, informed consent was required for all patients prior to registration.

Available tumor tissue samples and blood samples (N = 160; Supplementary Fig. S1) submitted for Lung-MAP screening were provided by the SWOG tissue bank. The clinical information for correlative studies in collaboration with the CIMAC–Cancer Immunologic Data Commons (CIDC) Network is shown in Supplementary Table S1 across the different assays.

mIF staining was performed in 82 screening tumor tissue samples (nivo+ipi = 38, and nivo = 42; Supplementary Table S1). Unstained slides from formalin-fixed, paraffin-embedded (FFPE) tissue were received from the SWOG bank and stained using methods previously described and validated (20). Briefly, 4-μm-thick FFPE tumor sections were stained using an automated staining system (Leica Microsystems) and two mIF panels with the following antibodies: Panel 1, cytokeratin (CK), CD3, CD8, PD-1/PD-L1, and CD68 and Panel 2, CK, CD3, CD8, CD45RO, granzyme B (GZB), and FOXP3. Antibody clones, dilutions, and RRIDs are included in Supplementary Table S2 and have been described previously (20). All the markers were stained in sequence using their respective fluorophore contained in the Opal 7-Color Automation IHC Kit (catalog no. NEL821001KT; Akoya Biosciences). The slides were scanned using the Vectra/Polaris 3.0.3 (Akoya Biosciences) at low magnification, 10× (1.0 μm/pixel) through the full emission spectrum and positive tonsil controls from the run staining to calibrate the spectral image scanner protocol (21). A pathologist selected representative areas inside the tumor using regions of interest for scanning in high magnification by the Phenochart Software image viewer 1.0.12 (931 × 698 μm size at resolution 20× = 0.5 μm/pixel) to capture various elements of tissue heterogeneity. Marker coexpression was employed to identify malignant cells (CK+), malignant cells expressing PD-L1 (CK+PD-L1+), and the cellular subsets of tumor-associated immune cells (TAIC) listed in Supplementary Table S3. Densities of each cell phenotype were quantified as the number of cells/mm² in the tumor compartment characterized by group or nests of malignant cells, in the stroma compartment characterized by the fibrous tissue present between the tumor nets, and in both compartments described as a total. PD-L1+ malignant cells were also expressed in percentages. The data were consolidated using R studio 3.5.3 (Phenopter 0.2.2 packet; Akoya Biosciences).

Using the point pattern distribution of the cell phenotypes relative to malignant cells, we measured the distance from malignant cells (CK+) to TAICs included in each mIF panel using R studio 3.5.3 (Phenopter 0.2.2 packet). We applied the median nearest neighbor function from malignant cells (CK+) to different cell phenotypes to determine where these TAICs were located; specifically, whether the TAICs were close to (i.e., equal to or less than the median distance) or far from (i.e., more than the median distance) the malignant cells (CK+) and associated with clinical outcomes.

Cells were subset by phenotype using the markers in the mIF panels and examined as the following: Tumor/PD-L1+ (CK+PD-L1+), Tumor (CK+), Other-Tcells (CD3+), Other-Tcells/PD-1+ (CD3+PD-1+), Macrophages (CD68+), Macrophages/PD-L1+ (CD68+PD-L1+), CTLs (CD3+CD8+), CTLs/PD-1+ (CD3+CD8+PD-1+), CTLs/GB+ (CD3+CD8+GB+), and Tregs (CD3+CD8−Foxp3+). The above phenotypes were used to visualize the spatial organization of cells by type. This analysis was carried out in R version 4.2.0 (R studio 2022.07.2).

Using the marked planar point pattern representations of each mIF image, we calculated the spatially varying probabilities for each of the phenotypes (described above). We used the spatstat toolbox (22) which provides the relrisk function to identify areas of segregation for a multitype (markers >2) marked point pattern. This function estimates for each phenotype, the spatially varying probability or the ratios of the probabilities, using kernel smoothing. The output of this function was used to plot the contour of the spatially segregated neighborhoods for each phenotype.

We identified cell clusters in each image using Euclidean distance and a hierarchical clustering method. A minimum cluster size of 10 cells and distance ≤ 20 μm was the requirement for clustering. The distance-based hierarchical clustering yielded the neighborhood information in a matrix. The cells that did not form clusters were labeled "Free_cell". The relative percentages of cells in each phenotype within a cluster were used to generate the heat map. We used the SPIAT library (SPIAT version 1.0.4) to identify cell clusters and made additions to the SPIAT functions as required for our analysis using R version 4.2.0.

DNA and RNA were coextracted from FFPE specimens received from the SWOG bank (Supplementary Table S1) and subjected to WES and gene expression. The RNA from a total of 38 FFPE samples (nivo = 23 and nivo+ipi = 15) passed the quality control (QC) and was run on the nCounter platform using the PanCancer Immune Profiling Panel (730 immune-related and 40 housekeeping genes) per the manufacturer's instructions. Briefly, samples were hybridized overnight at 65°C to probes; excess probes were washed using the automated prep station and then imaged on the digital analyzer. All runs included a Human Reference RNA control for batch correction. Data were processed and normalized with NanoString's nSolver analysis software (23). All samples passed the post-run QC metrics, and no batch effects were evident in the runs. In addition, gene expression profiles were deconvoluted by TIMER and nSolver advanced analysis tools to infer immune cells correlated to clinical outcomes.

WES analysis was conducted using the CIDC WES pipeline on tumor DNA from 50 tumors (nivo = 28 and nivo+ipi = 22; Supplementary Table S1) that passed the QC. DNA from paired peripheral blood mononuclear samples was used as germ line control. WES implements Gene Analysis Toolkit (24) best practices and identifies somatic variants using Sentieon TNScope and Haplotyper algorithms (25), respectively. Somatic variants are annotated using the Variant Effect Predictor software (26). The pipeline uses an ensemble of three callers, CNVkit (27), Sequenza (28), and Facets (29), to characterize tumor copy-number variation (CNV), and the CNV segments called by at least two callers were used to generate a high-confident consensus set. Sequenza and FACETS were used to estimate tumor purity and also PyClone-VI was utilized to infer clonal status of mutations (30). PyClone v 0.13.1 (31) was used to perform mutation clonality analysis. It is a Bayesian clustering method that enables mutations to be grouped into putative clonal clusters by integrating copy number, tumor purity (obtained from Sequenza), and variant allele frequency data.

We performed circulating serum analyte measurements using proximity extension assay (Olink) in 561 serum samples collected longitudinally from 160 patients (Supplementary Table S1). A series of 92 proteins, such as cytokines and soluble immune checkpoints included in the “immuno-oncology” panel, was measured as described previously (32). Protein levels were normalized using internal positive and negative controls and quantified as log₂ protein expressions (NPX), which were subsequently used as input for downstream analysis.

To evaluate whether the baseline biomarkers are prognostically associated with survival, we dichotomized biomarker data by the median and performed univariate survival analysis with the log-rank test. OS and PFS were evaluated. The Cox proportional hazards regression model was used for multivariate survival analysis (R package Survival, https://CRAN.R-project.org/package=survival; ref. 33). We included TMB (≥10 or <10 mutations per Mb), PD-L1 (≥5 or <5%), and other statistically significant biomarkers identified from univariate analysis in Cox models. Thresholds for TMB and PD-L1 were determined from previous clinical studies (18). To assess whether continuous biomarker data were associated with response and other clinical variables, we used nonparametric tests: Spearman rank correlation for continuous clinical variables, Mann–Whitney U test for categorical clinical variables with two groups, and Kruskal–Wallis test for categorical variables with more than two groups. In parallel, we also dichotomized biomarker data and used the χ² test for a robust assessment with responders. The Benjamini-Hochberg method (34) was used for multiple testing adjustment of P values. The analysis was performed on all samples and on samples in two treatment arms separately.

To explore the association of each baseline protein level with clinical outcomes from the Olink data, we used logistic regression models for best objective response and Cox proportional hazards models for PFS and OS (R package Survival, https://CRAN.R-project.org/package=survival, RRID: SCR_021137; ref. 33). In separate regression models, univariate analyses included only the protein expression values, while the multiple variable analyses adjusted for additional covariates (i.e., treatment, age, sex, race, smoking). Then, to investigate the longitudinal changes in serum protein associated with treatment, we used mixed linear models (R package Dream and lme4; refs. 35, 36) and the timepoints baseline, cycle 2 week 3, cycle 4 week 7, and cycle 5 week 9 to quantify the effect of these variables and additional relevant clinical parameters. These were analyzed with the treatment arms nivo+ipi and nivo. In our models, each protein NPX was considered an independent variable. In contrast, phase, timepoints, and treatments were considered dependent variables and other covariates as random effects. This approach allowed us to quantify the variance across proteins and approximate degrees of freedom of the hypothesis test for each protein, thereby minimizing false-positive results. We used F-tests for multiple coefficient comparisons and moderate t tests for single coefficient comparisons.

To identify the differences between responders and nonresponders at each timepoint and longitudinally, we used time as a dependent variable. We jointly modeled survival with cytokine expression [R packages lme4, rstanarm: Bayesian applied regression modeling via Stan (RRID:SCR_024605), bayestestR, bayesplot: Plotting for Bayesian Models (https://mc-stan.org/bayesplot/, RRID: SCR_024588); refs. 36, 37] to investigate the association of longitudinal protein levels with survival outcomes. The model used Cox proportional hazards and liner mixed regression and assessed the association of dynamic biomarker changes with survival outcomes. In the random intercept, the independent variable was the number of months from baseline to biomarker collection, set as a natural spline with three knots (at most three changing timepoints between baseline and progression/death). The dependent variable was the Olink analyte NPX value. In the survival analysis component, the independent variable includes the treatment arms. The convergence of the Markov chain Monte Carlo samples was assessed using several diagnostics: potential scale reduction factor, autocorrelation and trace plots, adequate sample size, and Monte Carlo standard error (32, 38–40). Finally, we used the FDR as the preferred method to correct for multiple hypothesis testing. The thresholds for significance in the mixed linear models for differential expression tests were a log₂ fold change of at least 0.5 and an FDR < 0.05. The joint model's threshold for significance was at least 1 unit increase in log₂ NPX expression and FDR < 0.05.

For integrative analysis, we applied recursive partitioning tree analysis (RPART; rpart library in R, https://cran.r-project.org/web/packages/rpart/vignettes/longintro.pdf) and random forest (refs. 41, 42; RF, randomForest and randomForestSRC libraries in R) on Olink (N = 159) and mIF (N = 82) data. We fitted RPART tree using responder status as the dependent variable, 92 baseline level Olink proteins and 17 mIF markers as predictors. We also created decision tree survival prediction model. Separate RPART trees were fitted for mIF markers from different compartments along with Olink proteins. The minimum number of observations in a node for a split was set to be 15; 10-fold cross-validation was carried out and results used for tree-pruning. For RF, we used 81 samples with both Olink and mIF data. Bootstrap the data to create bootstrap samples; grow a survival tree for each bootstrap sample with split criteria based on the log-rank statistics; continue the recursive partition; and calculate importance of each predictor by averaging over the forest.

In conjunction with the clinical study principal investigator/chair, the NCI-sponsored network and CIDC, deidentified data will be made publicly available by request under the dbGaP PHS accession number: phs003412.v1.p1 (https://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/study.cgi?study_id=phs003412.v1.p1). Questions and requests for additional data can be directed to the corresponding authors.

There were 31 responders (19.4%), including complete responses, partial responses, and unconfirmed partial/complete response; 63 patients with stable disease (39.4%); and 61 patients (38.1%) with progressive disease (i.e., increasing disease and symptomatic deterioration). Overall, 31 patients (19.4%) were alive at the end of the study, and the median OS was 10.02 months (range: 0.3–40.3). One hundred and forty-eight patients (92.5%) had disease progression, with a median PFS of 3.4 months (range: 0.3–36.6; Supplementary Table S1).

We first analyzed the mIF and gene expression profiling data from baseline tissue samples taken before ICI treatment to identify immune features associated with clinical benefit. mIF data revealed higher densities of various immune cells in the stroma compartment compared with the tumor compartment across the whole cohort (nivo+ipi and nivo arms), with no significant differences between the nivo+ipi and nivo arms (Supplementary Table S4). The overall immune cell densities were higher in the responders across both arms, although the difference did not reach statistical significance. In the whole cohort, higher median densities of PD-1+ cytotoxic T cells (CTL; CD3+CD8+PD-1+) in the stroma (>4.1 cells/mm², P = 0.042), presence of GZB+ CTLs in the tumor compartment (>0 cells/mm², P = 0.011), and higher median densities of memory T cells (CD3+CD45RO+; >23.4 cells/mm², P = 0.041) and PD-1+ T cells (CD3+PD-1+; >16.0 cells/mm², P = 0.023) in the total compartment (tumor plus stroma) were associated with longer PFS (Supplementary Table S5). Similarly, transcriptomic analysis demonstrated that patients having tumors with a higher expression of genes associated with myeloid infiltration, immune recruitment, and inflammation had superior clinical outcomes in the whole cohort (Table 1; Supplementary Table S6). The associations between higher expression of CD163, BLNK, IRF1, and FCGR2A with better OS (P < 0.05) and higher expression of MAPK11 with worse OS remained significant in subsequent multivariate analyses after adjustments for known predictive biomarkers, including TMB and PD-L1.

In the nivo arm, higher densities of memory T cells (CD3+CD45RO+) in the total compartment (median > 24.6 cells/mm², P = 0.028) and memory/regulatory T cells (Treg; CD3+CD8-CD45RO+FOXP3+) in the total compartment (median > 4.6 cells/mm², P < 0.001) and the stroma compartment (median > 12.0 cells/mm², P = 0.049) were associated with longer PFS (Table 2; Fig. 1A and B). Higher densities of memory/Tregs (CD3+CD8−CD45RO+FOXP3+) in the total compartment (median > 4.6 cells/mm², P = 0.026) were associated with better OS (Table 2; Fig. 1C). In the nivo+ipi arm, higher densities of PD-1+ T cells (CD3+PD-1+) in the total compartment (median > 16.0 cells/mm², P = 0.0347) and the presence of GZM+ CTLs (CD3+CD8+GZB+) in the tumor compartment (>0 cells/mm², P = 0.0154) were associated with longer PFS (Table 2; Fig. 1D and E). Conversely, higher densities of Tregs (CD3+CD8−FOXP3+) in the total compartment (median > 12.4 cells/mm², P = 0.0418) were associated with worse OS (Table 2; Fig. 1F). In the nivo+ipi arm, deconvolution of transcriptomic profiling data by TIMER and nSolver demonstrated significantly higher total immune cells (CD45+), a higher exhausted CD8+ T-cell score, and a higher neutrophil score in responders versus nonresponders (P < 0.05, Fig. 1G–I), further supporting that overall higher immune infiltration is associated with superior clinical benefit from ICI treatment.

Next, we specifically investigated exceptional responders, defined as patients who had no progression for at least 18 months and were still alive by 24 months, versus early progressors, who survived more than 1 month but had progressive disease and died within 6 months after initiating ICI treatment. By these definitions, there were 11 exceptional responders and 44 early progressors across the total trial cohort (Fig. 2A). There were more exceptional responders in the nivo+ipi arm than in the nivo arm (7 of 73 vs. 4 of 87, P = 0.35). Among these patients, 8 exceptional responders and 21 early progressors had tissues available for mIF, and 6 exceptional responders and 8 early progressors had tissues available for gene expression analysis.

By mIF, we observed higher densities of CTLs (CD3+CD8+) and memory CTLs (CD3+CD8+CD45RO+) in the total compartment in exceptional responders than in early progressors (CTLs: median, 152.1 vs. 27.3 cells/mm²; P = 0.032; memory CTLs: median, 31.2 vs. 2.1 cells/mm²; P = 0.040; Supplementary Table S7). Representative images from an exceptional responder are shown in Fig. 2B and C. Moreover, in the tumor compartment, we observed higher densities of GZB+ CTLs (CD3+CD8+GZB+) in the exceptional responders than the early progressors (median, 3.6 vs. 0 cells/mm², P = 0.027; Fig. 2D). Representative images from an early progressor showing a lower density of immune infiltration are shown in Fig. 2E and F.

Furthermore, distinctive spatial neighborhoods and cell organization in tumor microenvironment (TME) were observed in exceptional responders (n = 6) relative to early progressors (n = 6; Fig. 2G–J; Supplementary Figs. S2 and S3). Shown in Fig. 2G and Supplementary Figs. 2A and 3A is the distribution of different cell subsets relative to each other, with higher immune infiltration and higher CTL densities in the TME of exceptional responders versus early progressors. We then used spatially varying probabilities of different cell phenotypes to determine the segregation among immune subsets and malignant cells, and a contour plot to represent the neighborhoods within the TMEs. These analyses revealed a higher spatial segregation of immune cell subsets relative to malignant cells in the early progressors as compared with the exceptional responders (Fig. 2H; Supplementary Figs. S2B and S3B) in line with above observation that higher densities of CTLs in the tumor region positively associated with survival. Through distance-based hierarchical clustering, we identified local cell clusters within the TME and observed distinct compositions in exceptional responders versus early progressors. The clusters (cells ≥ 10 within interacting distance 20 μm) in exceptional responders often consisted of CTLs and other T-cell populations in the proximity to tumor cells (Fig. 2I and J; Supplementary Figs. S2C, S2D, S3C, and S3D). Finally, infiltration of neutrophils inferred from gene expression profiling data was significantly higher in exceptional responders than early progressors (P = 0.029; Supplementary Fig. S4).

Given the recognized importance of distance between different cells and the clustering of CTLs and tumor cells observed in exceptional responders, we expanded our analysis to understand the spatial relationship between the cell types associated with clinical outcome described above with other cells within the TME (43, 44). In the whole cohort, shorter distances from malignant cells (CK+) as well as PD-L1+ malignant cells (CK+PD-L1+) to CTLs (CD3+CD8+; median, 139 and 148 μm) were associated with better PFS (P = 0.045 and P = 0.027, respectively); shorter distances from malignant cells (CK+) to GZB+ CTLs (CD3+CD8+GZB+) were associated with significantly longer PFS (P = 0.035) and a trend toward longer OS (P = 0.054; Fig. 3A–E). In the nivo+ipi arm, shorter distances of CTLs (CD3+CD8+) as well as GZB+ CTLs (CD3+CD8+GZB+) from malignant cells (CK+; P = 0.045 and, P = 0.026, respectively) and shorter distances between CTLs (CD3+CD8+) from PD-L1+ malignant cells (CK+PD-L1+; P = 0.033) were associated with longer PFS (Fig. 3F–H). In addition, shorter distances of GZB+ CTLs (CD3+CD8+GZB+; median, 245 μm) from malignant cells (CK+) were associated with longer OS (P = 0.045; Fig. 3I). Taken together, these results suggest that the immune cells’ density and spatial distribution may impact response to ICI therapy.

We next performed WES (n = 50) with the intent of identifying the genomic basis underlying the immune features associated with benefit from nivo+ipi versus nivo in these metastatic squamous cell carcinomas. A total of 30,081 nonsilent mutations were detected with transversions, particularly C>A, as the predominant substitutes, which was expected because most patients were smokers (Supplementary Fig. S5A). The commonly mutated cancer genes included TP53, LRP1B, CDKN2A, AR1D1A, and PIK3CA (Supplementary Fig. S5B). High CNV burden was associated with a colder tumor immune microenvironment, as evidenced by lower infiltration levels of overall CD3⁺ T cells from mIF and lower levels of various immune signatures derived from gene expression profiling (Supplementary Fig. S6A and S6B). Importantly, CNV burden was not associated with estimated tumor purity, suggesting the correlation between high CNV burden and cold tumor immune microenvironment was not due to relative high tumor cell density leading to dilution of immune cells. Taken together, these results indicate that chromosomal instability may be an underlying genomic feature associated with immune evasion in metastatic SqNSCLC. Among the commonly mutated cancer genes, mutations in LRP1B, a recently recognized potential regulator of the inflammatory response, was associated with less infiltration of GZB+ CTLs (CD3+CD8+GZB+; Supplementary Fig. S7A; refs. 45, 46). Interestingly, LRP1B mutations were enriched in nonresponders but not in responders (18/19 vs. 2/11, P = 0.049). Furthermore, patients with LRP1B mutations had significantly (P = 0.008) shorter PFS (Supplementary Fig. S7B) and numerically shorter OS in the overall cohort (Supplementary Fig. S7C). LRP1B-mutant tumors were not associated with short PFS in the nivo+ipi arm (Supplementary Fig. S7D) but were in the nivo arm (P = 0.033; Supplementary Fig. S7E).

Blood-based biomarkers are attractive because they are noninvasive, dynamic, and less impacted by intratumor heterogeneity than tissue-based markers (47). We performed Olink proximity extension assay using the immuno-oncology panel assaying a series of 92 proteins in 561 serum samples collected longitudinally from 160 patients. Using mixed models to account for demographic and relevant clinical covariates with multiple testing adjustments, several serum chemokines (CXCL9, CXCL10, CXCL13, CCL19) and activated T-cell markers (PD-1, IFNγ, IL12, IL10) were found durably increased from baseline with either nivo or nivo+ipi (Fig. 4A and B), indicating the ICIs’ immune regulating effect. Multiple markers of immune activation and priming (ICOS-L, LAMP3/DC-LAMP, IL4, IL13, NRC1, CD5) were found increased at baseline or early during treatment in responders, regardless of treatment type (P < 0.05; Fig. 4C and D), in line with associations of these important immune processes with clinical response to ICI. Conversely, macrophage-derived and hyperinflammation markers, such as IL6, IL8, CXCL13, CSF-1, TNFSF14/LIGHT, and CCL23, as well as likely stromal or tumor-derived markers, such as VEGFA, HGF, and HO-1, were significantly upregulated in nonresponders at baseline or after ICI preceding radiologic progression, with some differences based on treatment received for CXCL13 and CSF-1 (P < 0.05; Fig. 4C and D). Joint modeling of survival with Olink analytes showed an increased risk of death (HR > 1) with higher longitudinal serum levels of CXCL13, MMP12, CSF-1, and IL8, which was confirmed with independent Kaplan–Meier analyses based on median protein levels at baseline (Fig. 4E). Similar results were generally observed in the subset of patients with extreme outcomes (exceptional responders and early progressors), where LAMP3/DC-LAMP was higher while CXCL13, CCL23, and TNFSF14 were lower in exceptional responders at baseline compared with nonresponders (Fig. 4F; P < 0.05). Together with the above-described data, and considering only baseline markers, these results suggest that an activated T-cell signature (cytotoxic effector T cells and DC-LAMP) was important for responsiveness to treatment with either nivo or nivo+ipi, while a hyperinflammatory milieu (IL6, IL8, CXCL13, CCL23, TNFSF14/LIGHT, CSF-1, MMP12) had an adverse impact on response and OS.

The antitumor immunity and response to ICIs is often determined at different molecular levels. The multiomics profiling in this study provided a unique opportunity for integrative analysis to understand the molecular and immune features associated with ICI benefit. We first performed recursive partitioning on Olink, mIF, NanoString, and WES data for classification of responders (Supplementary Fig. S8A and S8B). We identified that proteins from Olink provide good prediction on response. However, mIF markers did not contribute significantly in the decision tree, which might be due to relatively small sample size for mIF (n = 159 for Olink and n = 82 for mIF). We next created a decision tree survival prediction model and observed that Cytotoxic.T.cells.antigen.experienced (CD3+CD8+PD-1+) together with IL6, LAG3 and MICA.B separate patients into sub-populations with different survival. Furthermore, we applied random forest classifier, which identified Cytotoxic.T.cells.antigen.experienced (CD3+CD8+PD-1+) and IL6 as important variables. NanoString and WES did not contribute to the association between omics markers and outcomes likely due to insufficient samples with data from these platforms.

Identifying novel biomarkers for ICI response is challenging because the molecular determination of TME and host immune response is complex and heterogeneous across different patients. A large sample size to control interpatient heterogeneity and multi-omics to identify the determinates at different molecular levels are ideal but challenging. Therefore, maximizing the use of clinical, pathologic, molecular data and learning from each patient, particularly from clinical trials and careful analysis is key to pave the way to advance our understanding and ultimately the efficacy of ICI.

In this study, we performed mIF, gene expression profiling, WES, and OLINK on the previous samples from S1400I and identified known and novel molecular features associated with nivo versus nivo+ipi combination. Responders demonstrated higher densities of multiple immune cell types defined by mIF. Analysis of CTL populations revealed that GZB+ CTLs (CD3+CD8+GZB+) located in the tumor compartment were associated with better PFS. This was corroborated by analyzing the spatial organization of cell phenotypes, whereas higher immune cell population in the tumor region was seen in the TME of exceptional responders than that of early progressors. On the other hand, higher densities of Tregs (CD3+CD8−FOXP3+) in the total compartment correlated with worse OS in the nivo+ipi arm. This highlights the emerging dichotomy regarding the impact of ICI therapies on Treg subsets and function (48), as these combinations may not modulate some Treg subsets and dominance of CTLs is needed to overcome local immune suppression. Conversely, higher densities of memory T cells (CD3+CD45RO+) and Treg/memory T cells (CD3+CD8−CD45RO+FOXP3+) were associated with better PFS in the nivo arm. As PD-1 targeting has been shown to result in the reactivation of T cells already present within the tumor immune microenvironment, the presence of Treg/memory T cells at baseline may be an essential biomarker to delineate the need for inclusion of ipilimumab as opposed to a nivolumab alone approach. Estimating the immune subsets using TIMER and nSolver software demonstrated that a higher immune presence was associated with improved outcome. These results emphasize that an active immune response within the TME is required for a favorable clinical outcome in this setting, which is supported by multiple findings identifying mechanisms of response to ICI-based therapeutic strategies across various cancer types (49–51).

The TME is composed of various immune cells and stroma cells entangled with cancer cells. In addition to the densities of different cells, the spatial distribution and proximity among various cell types are also essential features with important impact on the functional status of the tumor immune microenvironment (52, 53). mIF data from this study provided an opportunity to assess the spatial relationship of different cellular components within the tumor immune microenvironment and their association with clinical outcomes from ICI treatment. Using the spatial point metrics through the nearest neighbor analysis, we observed that tumors with higher densities of GZB+ CTLs close to malignant cells were associated with better PFS and OS in the nivo+ipi arm, suggesting that cell-to-cell distribution and especially CTLs play an important role in response to ICIs as shown by others studies in NSCLC (54). The organization of cells into clusters based on distance also demonstrated that the CTLs and malignant cells cluster more frequently in exceptional responders than early progressors suggesting a preformed antitumor response that is aided by the ICI.

We used WES to identify genomic features underlying particular immune features and found that a higher CNV burden was associated with a lower level of immune cell infiltration overall. Similarly, CNV burden was negatively associated with immune scores derived from immune gene expression profiling. These findings are in line with previous findings in different cancer types suggesting that chromosomal instability may be a common genomic alteration underlying immune evasion across human malignancies (53, 55–58). Interestingly, we also found that patients with LRP1B-mutant tumors had a reduced survival compared with patients without LRP1B mutations. LRP1B has been identified as a putative tumor suppressor and is frequently inactivated in NSCLCs (45). Recently, LRP1B mutation was reported to be associated with better prognosis in melanoma and NSCLC after anti-PD-1 therapy (46). However, in our cohort, we observed that LRP1B mutation was associated with a worse OS and PFS in both the nivo arm and nivo+ipi combination therapy arms. It is still unclear whether the difference was due to different histology (predominantly adenocarcinoma in the previous study vs. exclusively squamous cell carcinoma in the current study) or different ICI (anti-PD-1 vs. anti-PD-1 with/without anti-CTLA-4) or low sample size. Of note, the impact of LRPB1 mutations on cancer biology and response to ICIs has not been clearly defined in different cancers. For example, a study on renal clear-cell carcinoma reported worse prognosis and suppressive antitumor immunity when LRPB1 was overexpressed (59), and another found that LRPB1 mutations were associated with inferior clinical outcomes after ICI treatment in patients with hepatocellular carcinoma (60).

Although tissue-based assays remain the gold standard for molecular profiling for oncology practice, liquid biopsy, particularly peripheral blood–based assays, have gained more attention for molecular profiling and disease monitoring across various cancers because they are noninvasive, “real-time,” and less affected by intratumor heterogeneity (61, 62). In the era of immune-oncology, the Olink soluble protein detection platform has emerged as a promising tool to assess and monitor host immune response. Using Olink, we identified a high level of protumorigenic factors, such as VEGFA and CCL23, and inflammatory markers, such as IL6, IL8, and MMP12, that were associated with inferior survival in this cohort of patients. These findings suggest that general inflammation is detrimental in the context of cancer and ICI therapy. In contrast, proteins involved in T-cell and natural killer cell activation, such as LAMP3/DC-LAMP, IFNγ/IL4/IL13, and NRC-1, were associated with improved outcomes after ICI therapy. It was unexpected that a high level of CXCL13 was associated with poor response to ICI therapy and shorter survival, given the recent studies reporting this chemokine working together with DC-LAMP and playing essential roles in the establishment of tertiary lymphoid structures in NSCLC (63). It is possible that the relatively high levels of circulating CXCL13 in the serum do not reflect relatively rare CD4⁺ T cell–derived tumor tissue–specific expression of CXCL13, and this emphasizes the limitations of soluble analytes as a surrogate for local tumor events. Some analytes, such as CXCL9/10 and soluble PD-1, were dynamically increased with treatment and marginally associated with outcomes in exceptional responders, in line with previous reports (64). Of particular interest, other markers showed the strongest association with objective response during ICI treatment, for example, lower CSF-1 or IL6 or higher IL13 at cycle 2. This is reminiscent of findings from melanoma studies in which the on-treatment biopsy was more informative for long-term benefit from ICIs than the baseline biopsy, as the actual changes after treatment reflect the host immune system's response to ICIs more accurately (65, 66). Although soluble analytes are not ideal predictive biomarkers to select an optimal initial treatment regimen, if validated, these markers will be extremely helpful in switching ineffective therapy to effective alternatives to save time and potential toxicity, which is critically important for patients with stage IV lung cancers, for whom time and quality of life are essential attributes. In addition, these on-treatment markers are also valuable in distinguishing pseudoprogression from real progression—another critical clinical dilemma that the oncologists face in the era of immuno-oncology.

As a post hoc profiling of samples from a completed clinical trial, our study has several inherent limitations, including inadequate tumor specimen availability, which precluded us from generating comprehensive data integration from all platforms; imbalanced distribution of samples from the nivo versus nivo+ipi arms or responders versus nonresponders; inadequate tissues for multiomic analysis and cross-platform integrative analyses; and lack of detailed information regarding the time and anatomic sites of tumor specimens, which limited our ability to perform in-depth, organ-specific analysis. In spite of these challenges, integration of peripheral cytokine profiling and cellular profiling within the TME confirmed our single platform findings highlighting the negative association with hyperinflammation with reduced PFS and the presence of PD-1+ CTLs in the TME with increased survival. Finally, while we presented several candidates in this study, we recognize the need for additional validation and replication of our findings. Specifically, several circulating serum proteins, such as IL6, IL8, CSF-1, MMP12, and CXCL13 are promising candidates for future prospective or post hoc confirmatory studies due to their ease of collection and quantification from blood. In addition, investigation of tissue composition using spatial profiling technologies to better understand the complex interplay between tumor tissue and immune infiltrating cells may shed light on the mechanisms of immune-tumor cell–cell interactions and identify key biomarkers that can identify patients who will have the most benefit from ICIs. As a proof-of-principle study, using the S1400I trial as an example, we showcased that multi-omics, multi-institutional analyses of patient samples are feasible and can provide valuable insights for future trial development, which is one of the major goals of the CIMAC-CIDC Network.

D.Y. Duose reports grants from NIH during the conduct of the study, as well as other support from Chrysalis Biomedical outside the submitted work; in addition, D.Y. Duose has been issued patent application # 62/866,130 entitled “Methods for Duplex Sequencing of Cell-Free DNA and Applications Thereof.” F. Michor is a co-founder of and has equity in Harbinger Health and equity in Zephyr AI, and serves as a consultant for both companies; in addition, F. Michor is on the board of directors of Exscientia Plc. J.E. Gray reports personal fees from AbbVie, Axiom HC Strategies, Blueprint Medicines, Celgene, Daiichi Sankyo, Inc, EMD Serono - Merck KGaA, Inivata, Janssen Scientific, Jazz Pharmaceuticals, Loxo Oncology Inc, OncoCyte Biotechnology, Sanofi Pharmaceuticals, and Takeda Pharmaceuticals; grants from Boehringer Ingelheim, G1 Therapeutics, Ludwig Institute of Cancer Research, and Pfizer; and grants and personal fees from AstraZeneca, Bristol Myers Squibb, Genentech, Merck & Co, Inc, and Novartis outside the submitted work. R.S. Herbst reports other support from AstraZeneca, personal fees and other support from Eli Lilly and Company, and personal fees from Genentech/Roche and Merck and Company during the conduct of the study. R.S. Herbst also reports personal fees from AbbVie Pharmaceuticals, Immunocore, Junshi Pharmaceuticals, Bristol Myers Squibb, Candel Therapeutics, Cybrexa Therapeutics, DynamiCure Biotechnology, eFFECTOR Therapeutics, EMD Serono, Gilead, HiberCell Inc, I-Mab Biopharma, Immune-Onc Therapeutics, Inc, Janssen, Johnson & Johnson, Loxo Oncology, Mirati Therapeutics, NextCure, Novartis, Ocean Biomedical, Inc, Oncocyte Corp, Oncternal Therapeutics, Pfizer, Regeneron Pharmaceuticals, Revelar Biotherapeutics, Inc, Ribbon Therapeutics, Sanofi, Seattle Genetics, and Xencor, Inc; other support from Bolt Biotherapeutics; and personal fees and other support from Checkpoint Therapeutics and Normunity outside the submitted work. In addition, R.S. Herbst is a board member for American Association for Cancer Research and International Association for the Study of Lung Cancer, as well as committee member for Society for Immunotherapy of Cancer and Southwest Oncology Group. I.I. Wistuba reports grants and personal fees from Genentech/Roche, Bayer, Bristol Myers Squibb, AstraZeneca, Pfizer, Merck, Guardant Health, Sanofi, and Novartis; personal fees from GlaxoSmithKline, Flame, Daiichi Sankyo, Janssen, Merus, Regeneron, Oncocyte, and MSD; and grants from Adaptive, Adaptimmune, EMD Serono, Karus, Johnson & Johnson, Iovance, 4D, and Akoya outside the submitted work. S. Gettinger reports grants from Bristol Myers Squibb outside the submitted work. K. Kelly reports other support from BMS during the conduct of the study. L. Bazhenova reports personal fees from Neuvogen, BioAtla, Summit Therapeutics, Anheart, Gilead, Merck, AbbVie, Regeneron, Janssen, Genentech, Novocure, Intervenn, Elevation Oncology, ORIC, Mirati, Turning Point, Daiichi, Bayer, Sanofi, Novartis, BMS, and BI outside the submitted work. S. Gnjatic reports grants from Regeneron, Boehringer-Ingelheim, Genentech, Takeda, BMS, Celgene, and Janssen R&D outside the submitted work. J. Zhang reports grants from Merck; grants and personal fees from Johnson & Johnson and Novartis; and personal fees from Bristol Myers Squibb, AstraZeneca, GenePlus, Innovent, Hengrui, and Varian outside the submitted work. C. Haymaker reports grants from NCI during the conduct of the study; C. Haymaker also reports grants from Sanofi, Dragonfly, BTG, Trisalus, Iovance, Avenge, and Obsidian, as well as personal fees from Briacell, Nanobiotix, SWOG, and SITC outside the submitted work. No disclosures were reported by the other authors.

E.R. Parra: Conceptualization, investigation, methodology, writing–original draft, writing–review and editing. J. Zhang: Formal analysis, methodology, writing–original draft, writing–review and editing. D.Y. Duose: Investigation, methodology, writing–original draft, writing–review and editing. E. Gonzalez-Kozlova: Investigation, visualization, methodology, writing–original draft, writing–review and editing. M.W. Redman: Formal analysis, investigation, writing–original draft, writing–review and editing. H. Chen: Formal analysis, investigation, methodology, writing–review and editing. G.C. Manyam: Formal analysis, investigation, methodology, writing–review and editing. G. Kumar: Visualization, methodology, writing–review and editing. J. Zhang: Investigation, methodology, writing–review and editing. X. Song: Investigation, methodology, writing–review and editing. R. Lazcano: Investigation, writing–review and editing. M.L. Marques-Piubelli: Investigation, writing–review and editing. C. Laberiano-Fernandez: Investigation, writing–review and editing. F. Rojas: Investigation, writing–review and editing. B. Zhang: Investigation, writing–review and editing. L. Taing: Investigation, methodology, writing–review and editing. A. Jhaveri: Investigation, methodology, writing–review and editing. J. Geisberg: Investigation, methodology, writing–review and editing. J. Altreuter: Resources, investigation, methodology, writing–review and editing. F. Michor: Resources, investigation, methodology, writing–review and editing. J. Provencher: Resources, investigation, methodology, writing–review and editing. J. Yu: Resources, investigation, methodology, writing–review and editing. E. Cerami: Resources, investigation, methodology, writing–review and editing. R. Moravec: Resources, writing–review and editing. K. Kannan: Investigation, methodology, writing–original draft, writing–review and editing. R. Luthra: Resources, investigation, writing–review and editing. G. Alatrash: Funding acquisition, writing–original draft, writing–review and editing. H.-H. Huang: Investigation, methodology, writing–review and editing. H. Xie: Investigation, methodology, writing–review and editing. M. Patel: Investigation, methodology, writing–review and editing. K. Nie: Investigation, methodology, writing–review and editing. J. Harris: Investigation, methodology, writing–review and editing. K. Argueta: Investigation, methodology, writing–review and editing. J. Lindsay: Resources, investigation, methodology, writing–review and editing. R. Biswas: Resources, investigation, methodology, writing–review and editing. S. Van Nostrand: Investigation, methodology, writing–review and editing. S. Kim-Schulze: Resources, investigation, methodology, writing–review and editing. J.E. Gray: Investigation, writing–review and editing. R.S. Herbst: Investigation, writing–review and editing. I.I. Wistuba: Supervision, funding acquisition, writing–original draft, writing–review and editing. S. Gettinger: Investigation, writing–original draft, writing–review and editing. K. Kelly: Investigation, writing–original draft, writing–review and editing. L. Bazhenova: Investigation, writing–original draft, writing–review and editing. S. Gnjatic: Supervision, writing–original draft, writing–review and editing. J.J. Lee: Data curation, formal analysis, supervision, validation, methodology, writing–original draft, writing–review and editing. J. Zhang: Conceptualization, supervision, investigation, writing–original draft, writing–review and editing. C. Haymaker: Conceptualization, supervision, funding acquisition, investigation, writing–original draft, writing–review and editing.

Scientific and financial support for the CIMAC-CIDC Network is provided through the NCI Cooperative Agreements U24CA224319 (to the Icahn School of Medicine at Mount Sinai CIMAC), U24CA224285 (to MD Anderson Cancer Center CIMAC) and U24CA224316 (to the CIDC at Dana-Farber Cancer Institute). SWOG Study S1400I was supported in part by NIH/NCI grants U10CA180888 and U10CA180819 and in part by Bristol Myers Squibb Company, through the Foundation for the NIH, in partnership with Friends of Cancer Research. Scientific and financial support for the Partnership for Accelerating Cancer Therapies (PACT) public-private partnership (PPP) are made possible through funding support provided to the FNIH by: AbbVie Inc., Amgen Inc., Boehringer-Ingelheim Pharma GmbH & Co. KG., Bristol-Myers Squibb, Celgene Corporation, Genentech Inc, Gilead, GlaxoSmithKline plc, Janssen Pharmaceutical Companies of Johnson & Johnson, Novartis Institutes for Biomedical Research, Pfizer Inc., and Sanofi. S. Gnjatic was partially supported by NIH grants CA224319, DK124165, CA263705 and CA196521. Computational and data resources at Mount Sinai were partially supported by the Clinical and Translational Science Award (CTSA) grant UL1TR004419.

We thank the members of the Translational Molecular Pathology Immune-Profiling (TMP-IL) MoonShots Platform at MD Anderson for their support. We thank Beatriz Sanchez-Espiridion, Julia Mendoza Perez, Shani Wijeratne, and Celia Garcia-Prieto for their assistance with sample procurement and inventory; Mei Jiang, Auriole Tamegnon, Heladio Ibarguen, Salome McAllen, Saxon Rodriguez, Jianling Zhou, Mauricio E. Vasquez Juarez, and Renganayaki Krishna Pandurengan for their technical assistance; Dawen Sui and Juan Ruiz Posadas for data transfer; Curtis Gumbs and Latasha Little for genomics support; and Seunghee Kim-Schulze, for his help in running Olink assays. We thank Diane Del Valle for project management from the Mt. Sinai CIMAC; Magdalena Thurin, Helen Chen, Minkyung Song from NCI; and Rebecca Enos from Emmes. Editorial support was provided by Bryan Tutt, Scientific Editor, Research Medical Library at MD Anderson Cancer Center. We thank the patients and their families for participating in the study.