This was a prospective study of patients with stage III or IV NSCLC (n=15) diagnosed with adenocarcinoma, squamous cell carcinoma, or NSCLC NOS (not otherwise specified), recruited to the study before planned treatment start. Informed consent was obtained from the Regional Ethics Review Board in Gothenburg, Sweden (Permit number 953/18). Patients were recruited from April 2019 to April 2020, with a pause during the summer months of May to July ( Table 1 ). The cohort included patients with single PD-1 (n=14) or PD-L1 (n=1) blockade as 1^st (n=10) or 2^nd (n=5) line therapy. Formalin-fixed paraffin-embedded (FFPE) tumor tissues were obtained from the cohort at the time of diagnosis ( Supplementary Figure 1A ). Participation in the study did not influence the course of the treatment or clinical procedures. Peripheral blood was collected from each patient in K2EDTA tubes at five-time points, before and after each consecutive 2–4-week treatment cycle, during routine clinical assessment ( Supplementary Figure 1A ). Blood was also collected in K2EDTA tubes on one occasion from healthy controls (n=3), who were age- and sex-matched with the patient cohort. PD-L1 staining of tumor biopsies before PD-1 blockade therapy was assessed according to routine clinical testing using PD-L1 28-8 antibody (Dako). The pathologist defined PD-L1 protein expression as the percentage of tumor cells exhibiting positive membrane staining at any intensity.

The clinical response to PD-1 blockade was determined every 3 months after obtaining the results of the CT-Scan, in line with the immune-related Response Evaluation Criteria in Solid Tumors (irRECIST) algorithm but assessed by an oncologist according to clinical judgment (17). The clinical response was divided into (1) complete response, no measurable tumor; (2) partial response, shrinkage in tumor size compared with baseline; (3) stable disease, no change in tumor size compared with baseline; and (4) progressive disease with increase in tumor size compared with baseline. In this study, responders were defined as patients who maintained a complete response, partial disease, or stable disease at 9 months (3^rd assessment) after the start of therapy to certify that they were indeed clinically responding to therapy. In one patient, the assessment was made 9-10 months post-treatment. Non-responders were those patients with a progressive disease before or at 9 months ( Table 1 ).

Fresh whole blood in K2EDTA was stored at 22-24⁰C for 12–18 hours before FACS staining and analysis. Blood samples (1ml) were stained with antibody mix ( Supplementary Table 1 ) and incubated for 20 minutes at 4⁰C. After lysis of red blood cells with a lysing solution at room temperature and two rounds of centrifugation and washing of the cells using a lyse wash assistant (Becton Dickinson, BD), the cells were acquired using a BD FACSCanto analyzer. Data analysis was performed using the Diva (BD) software based on the standardized phenotyping of human immune cells (9). Lymphocytes were enumerated using BD Multitest™ 6-color TBNK reagent (BD Biosciences) according to the manufacturer’s instructions.

Genomic DNA was extracted from paired blood samples and FFPE tumor biopsies collected from the same patient. DNA isolated from the blood was used as a genomic DNA control for tumor tissue. Tumor tissue biopsies from 14 patients ( Table 2 ) were included in the genetic analysis. The NGS panel used included 597 genes (Oncopanel All in One v2.8, Eurofins Genomics (Europe Sequencing GmbH, Germany), and the design covered 10 bp flanking regions of all exons ( Supplementary Table 2 ). All steps from extraction, quantification, library preparation (Agilent Technologies Santa Clara, CA, USA), and sequencing were performed at Eurofins Genomics using optimized in-house protocols. Sequencing was performed on the Illumina NovaSeq 6000 platform (Illumina, San Diego, CA, USA) using 2 × 150 bp paired-end reads.

Data quality assessment, mapping, and variant calling were performed using Eurofins Genomics in-house pipeline (Europe Sequencing GMB, Germany). Variant filtration was performed using Alissa Interpret software (Agilent Technologies, Santa Clara, CA, USA). A cutoff of 5% presence of the mutational allele was used for all filtrations. Alamut Visual (version 2.15; Sophia Genetics, Lausanne, Switzerland) and cBioportal (18, 19) were used to interpret variants. Mutations were further classified as a benign, likely benign, variant of unknown significance (VUS), likely pathogenic, or pathogenic, using the model described by Froyen et al. (20) and in accordance with ACMG and AMP guidelines. For mutational signature analysis, FastQC (version 0.11.2) was used to assess the quality of the data, and samtools (version 1.9) were used to sort, index, and assess mapping statistics. Paired-end reads were aligned to the human reference genome (hg19) using Burrows-Wheeler Aligner (BWA mem version, BWA_0.7.13) (21). Picard (version 2.2.4) was used to remove duplicates. The Genome Analysis ToolKit (GATK, version 4.1.3.0) (22) was used for base quality score recalibration, Mutect2 was used for calling and filtering somatic variants, and SigProfiler Extractor was used to extract mutational signatures (23). An analysis of variants in a selection of genes classified as pathogenic, likely pathogenic, or VUS in TP53, KRAS, serine/threonine kinase 11 (STK11), KEAP1, and NOTCH1 was combined into a gene score for each patient, where 0 indicates no variants, gene score of 1 indicates a single pathogenic/likely pathogenic variant or VUS in one gene, gene score of 2 indicates two pathogenic/likely pathogenic variants or VUS in one gene or a pathogenic/likely pathogenic variant or VUS in two genes, and a gene score of 3 indicates three pathogenic/likely pathogenic variants or VUS in one gene or a pathogenic/likely pathogenic variant or VUS in three different genes ( Table 2 ).

TMB was calculated using Eurofins Genomics by dividing the number of mutations by the size of the targeted coding region in megabases (Mb). Only non-synonymous missense variants were included in the calculation. The calculation was performed using the following exclusion criteria: non-coding mutations, mutations listed as known somatic mutations (according to cosmic v71), known germline mutations (in dbSNP), mutations with depth below 50x and allele frequency below 0.05, germline mutations with more than two counts in genome AD mutations (https://gnomad.broadinstitute.org/) in tumor suppressor genes (24, 25). Patients with a TMB of 10 mutations/Mb or higher are referred to as TMB-H ( Table 2 ).

Principal component analysis (PCA) was performed using the pca3d package in R. To examine the relationship between clinical responder or nonresponder (Y-variable) and the immune and tumor genetic signatures of the patients studied (X-variables), orthogonal partial least squares discriminant analysis (OPLS-DA) was performed using the SIMCA-P+ software (Sartorius GmbH, Göttingen, Germany). The quality of the OPLS-DA was based on the parameter R2Y, that is, the model’s goodness of fit (values ≥0.5, which define good discrimination and best possible fit, R2Y=1), and Q2, the goodness of prediction of the model (26). A Q2 value >0.4 is considered satisfactory with biological variables. Furthermore, the difference between the Q2 and R2Y values should not exceed 0.4. A combination of variable influence on projection (VIP) and VIPcvSE was used to exclude variables that were less likely to contribute to building the model. VIPcvSE is the confidence interval of the VIP. The significance of the separation between the groups in the OPLS-DA was calculated using CV-ANOVA (26).

Non-parametric Mann–Whitney U-tests or Wilcoxon matched-pair signed rank tests were performed for unpaired and paired analyses, respectively, using the GraphPad Prism software (GraphPad Software, San Diego, USA). Progression-free survival (PFS) was estimated using the Kaplan–Meier method. The log-rank test was used to assess the differences in overall survival and PFS between the groups. Statistical significance was set at p <0.05, and no adjustments were made for multiple comparisons. Data analysis was performed using IBM SPSS Statistics version 27 (IBM, New York, USA) and GraphPad Prism software.

In this patient cohort, the clinical assessment at the 9-10-month follow-up was based on iRECIST criteria. Responders (n=10) included one patient with a complete response, six patients with partial response and three patients with stable disease ( Table 1 and Figure 1A ). Non-responders included patients with progressive disease observed at 3 (n=2), 6 (n=2), or 9 (n=1) months post-treatment ( Figure 1A ).

The events leading to clinical response or progression after PD-1 blockade therapy are multi-factorial, and the immune status before treatment initiation could impact the subsequent clinical response. We first investigated the baseline frequencies and fold-change after treatment of NK and B cells in our patient cohort ( Supplementary Figure 1B ). Although we could not quantify any significant change in the frequencies of circulating B cells (CD19⁺) or NK cells (CD3^-CD56⁺) in the blood before and after treatment (data not shown), we found a trend for a higher baseline frequency of B cells (p=0.15) and NK cells (p=0.11) in responders than in non-responders ( Figure 1B ). The analysis of the frequencies of CD3⁺ T cells at baseline revealed a trend for higher frequencies of T cells in non-responders (p=0.07) than in responders ( Figure 1B ). However, no significant change was observed in the frequencies of CD4⁺ and CD8⁺ T cells in the blood before and after the 1^st cycle of treatment (data not shown) or in the baseline frequencies of CD4⁺ T cells and CD8⁺ T cells between responders and non-responders ( Figure 1C ). Finally, to test the hypothesis that increased regulatory T cell (Treg) frequencies correlate with poor prognosis, we analyzed the frequencies of CD4⁺CD25⁺CD127^-CD45RO⁺CCR4⁺ Tregs in circulation ( Supplementary Figure 1B ). We did not observe any difference in baseline levels of Tregs between responders and non-responders ( Figure 1C ), nor any change that could be associated with clinical response after the 1^st cycle of treatment (data not shown).

In summary, while frequencies of the bulk CD4⁺ and CD8⁺ populations could not predict treatment response, a trend for higher baseline frequencies of B cells and NK cells in responders and CD3⁺ T cells in non-responders was observed after PD-1 blockade in our cohort of patients.

Since the frequencies of bulk T cell populations did not change after treatment, we continued to analyze specific subsets of effector or memory T cells using a clinical assay based on the expression of CD45RA and CCR7 markers and activation markers CD38 and HLA-DR (9) ( Supplementary Figure 1B ). To determine the time at which changes in the frequencies of the immune cell subsets were observed after treatment, we analyzed pre-treatment and four cycles post-treatment (A-E) fresh whole blood samples. Our analysis revealed that for all patients except one with progressive disease, there was an increase in the frequencies of activated effector memory CD4⁺ and CD8⁺ T cells after the 1^st or 2^nd treatment cycle, followed by a decrease at the 3^rd or 4^th treatment cycle, reaching either pre-treatment levels or somewhat higher levels ( Supplementary Figure 1C ). In the responder group, an increase (p<0.01) in the frequency of both CD4⁺ and CD8⁺ activated effector memory (CD45RA^-CCR7^-CD38⁺HLA-DR⁺) T cells was observed after the 1^st cycle of treatment compared with the pre-treatment levels ( Figures 2A ). When analyzing the fold-change between pre- and post treatment values, a significant difference between responders and non-responders (p=0.048) was only observed for the activated effector memory CD8⁺ T cell subset ( Supplementary Figure 2A ).

Analysis of central memory (CD45RA^-CCR7⁺) CD4⁺ and CD8⁺ T cells for all treatment cycles revealed a decrease in frequency after the 1^st treatment cycle in 6/10 patients in the responder group and 2/5 patients in the non-responder group ( Supplementary Figure 1D ). The frequencies of central memory cells increased over time in individual patients and reached pre-treatment frequencies at the 4^th treatment cycle for most patients, independent of clinical response ( Supplementary Figure 1D ). The trend for a decrease in the frequency of central memory CD4⁺ and CD8⁺ T cells at 1^st cycle after treatment compared to pre-treatment was however not significant ( Figure 2B ). Further, we observed no fold-change differences in the frequencies of central memory CD4+ or CD8+ T cells between the responders and non-responders ( Supplementary Figure 2B ).

Finally, the analysis of effector (CD45RA⁺CCR7^neg) CD4⁺ and CD8⁺ T cells did not reveal any changes in the pre- and post-treatment levels over time in either the responder or the nonresponder group ( Supplementary Figures 1E ). Interestingly, a high frequency of CD4⁺ effector T cells was observed in one patient with a clinically durable response (>2 years) ( Figure 2C ). The fold-change difference in the frequencies of effector T cells pre- and post-treatment could not distinguish the reponders from non-responders ( Supplementary Figure 2C ).

In summary, in our cohort of patients, a post-treatment significant increase in the frequency of CD8⁺ activated effector memory cells could indicate clinical response to PD-1 blockade.

We performed NGS of isolated DNA from pre-treatment tumor biopsies to explore the potential of a combined genetic and immune cell signature as an early biomarker of clinical response to PD-1 blockade. We analyzed gene variants in tumor-specific genes, including TP53, STK11, KRAS, and KEAP1 ( Supplementary Table 2 ), and an in silico panel of 59 immune-related genes. The analysis revealed that TP53 mutations were frequent, with missense, nonsense, in-frame, and splice mutations detected in 10/14 patients, independent of the clinical response. Our analysis also showed no significant difference in the time to progression at 9 months for patients with wild-type or mutated TP53 ( Supplementary Figure 3A ). Furthermore, in our patient cohort, KRAS mutations (5/10 patients) or KRAS and TP53 mutations were detected only in the responders ( Table 2 ). Although published data (27) suggest that patients with mutations in KRAS^mut or TP53 treated with PD-1 blockade experience a longer PFS than patients with wildtype (WT) tumors, we observed only a trend in patients with KRAS ^mut tumors (p=0.072), possibly because of the low number of patients in this study cohort ( Figure 3A ).

Analysis of the 59-gene immune panel revealed truncating and missense variants in several genes. However, only one nonsense mutation in the tumor suppressor gene NOTCH1 in responders was predicted to be likely pathogenic. Variants in the additional genes were detected but predicted to be VUS ( Supplementary Table 3 ). The number of pathogenic, likely pathogenic, and VUS mutations detected in individual patients was the basis for a gene score of ≥ 2 in 7/10 patients in the responder group and ≤ 1 in 5/5 patients in the nonresponder group. This indicates that pathogenic and likely pathogenic mutations were more common in the responders ( Table 2 and Supplementary Table 3 ). We conclude that both the number of variants in genes and the pathogenicity of these variants could contribute to the clinical response. Thus, interpreting the pathogenicity of gene variants could be important to consider for future studies.

Finally, we investigated the relationship between TMB score, PDL-1 status, and clinical response in our cohort. The patients with TMB-H (>10 mutations/Mb) were distributed among both responders and non-responders ( Table 2 ). Neither TMB-H nor PD-L1 expression (>1%) could be associated with differences in the time to PFS between responders and non-responders ( Supplementary Figure 3B ).

In summary, we found that driver mutations, particularly in KRAS, was related to clinical response after PD-1 blockade but not TMB score or PD-L1 expression.

Next, we performed multivariate analyses to test our hypothesis by analyzing the relationship between immune and genetic variables and their relationship to clinical response after PD-1 blockade. The baseline frequencies of B cells, NK cells, CD3 ⁺ T cells, CD4 ⁺ and CD8 ⁺ T cells, and Tregs were included in the analysis. In addition, we also added pre- and post- 1^st cycle of treatment, fold-change activated effector memory T cells, central memory T cells, effector T cells, and Tregs. The clinical and genetic parameters included the gene score, PD-L1 status, and TMB score. PCA, including all these variables, indicated a separation between responders and non-responders, although only about 45% of the variance was explained by the two principle components ( Figure 3B ). To further define the variables most important for discrimination between the groups, OPLS-DA was performed using a VIP cutoff >0.8 and a VIPcvSE cutoff <1.5. The analysis showed that the two clinical response groups, responders and non-responders, were separated with high discrimination and predictability (R2Y:0.95, Q2:0.847, P=0.01) ( Figure 3C ) based on the combination of immune and genetic variables. OPLS-DA analysis, including only immune or genetic variables alone, resulted in poor discrimination and predictability values (immune parameters: R2Y: 0.756; Q2: 0.389 and genetic parameters R2Y: 0.52; Q2: 0.119). The variables best defining the clinical responder group were baseline frequencies of B cells, fold-change activated effector memory CD8⁺ T cells, and gene score >1 ( Figure 3D ). In summary, we report that a combination of genetic and immune parameters analyzed before or within 3–4 weeks after the start of PD-1 blockade therapy could differentiate responders from non-responders.