Details of the study have been described previously ( Supplementary Figure S1 ) (ClinicalTrial.gov Identifier: NCT03958097) (8). Briefly, 20 patients enrolled the trial, all the patients received Sintilimab 200mg and 3×10⁹ NK cells every 3 weeks until progressive disease, unacceptable toxicity, or withdrawal of consent. The study protocol was approved by institutional review board of first hospital of Jilin University study center, and the study was done in accordance with standards of Good Clinical Practice and the Declaration of Helsinki. All patients provided written informed consent before enrollment.

Aged ≥18 years, with pathologically confirmed locally advanced or metastatic NSCLC who had failed the first line of platinum-based chemotherapy were eligible for the study. Patients with EGFR or ALK mutation were excluded. An Eastern Cooperative Oncology Group performance status of 0 or 1 and measurable disease according to Response Evaluation Criteria in Solid Tumors version 1.1. Patients previously received immunotherapy (ICIs, Chimeric Antigen Receptor T cells, antitumor vaccine etc.) were not included.

The primary end point of the study was PFS. Secondary endpoints included ORR, OS, duration of response (DoR), disease control rate (DCR), clinical benefit rate (CBR) and safety. The exploratory predictive biomarkers were Tumor Mutational Burden (TMB), PD-1 and PD-L1 expression on NK cells.

Cell-free DNA (cfDNA) was isolated using the QIAamp Circulating Nucleic Acid Kit (Qiagen, Hilden, Germany). The germline genomic DNA from white blood cells was isolated using the QIAamp DNA Blood Mini Kit (Qiagen, Valencia, CA, USA). The DNA concentration was measured using a Qubit fluorometer and the Qubit dsDNA HS (High Sensitivity) Assay Kit (Invitrogen, Carlsbad, CA, USA). The size distribution of the cfDNA was assessed using an Agilent 2100 BioAnalyzer and a DNA HS kit (Agilent Technologies, Santa Clara, CA, USA) as previously reported (9, 10).

Sequencing libraries were prepared using the KAPA Library Preparation Kit (Kapa Biosystems, Wilmington, MA, USA). A custom-designed panel covering ~1.5Mbp genome and targeting 1021 cancer-related genes was used for hybridization enrichment with DNA libraries, and then sequenced using Gene+ Seq 2000 instrument (GenePlus-Beijing) (9, 10).

After removing adapters and low-quality reads, the clean reads were mapped to the human reference genome (hg19) using BWA18 (version 0.7.12-r1039). Reads realignment and recalibration were performed with the Picard software MarkDuplicates (v4.0.4.0; Broad Institute, Cambridge, MA, USA). Somatic single nucleotide variants (SNVs) were determined by MuTect and NChot (Geneplus-Beijing, inhouse), small insertions and deletions (indels) were called using GATK (v3.6-0-g89b7209; Broad Institute), copy number alterations (CNVs) were detected with CONTRA (2.0.8), and structural variations (SVs) were identified with the algorithm NCsv (v0.2.3 Geneplus-Beijing, inhouse). All final candidate variants were manually verified with the integrative genomics viewer browser. Somatic alterations were filtered with matched patient’s blood cell control to remove germline mutations (9).

TMB was calculated as the number of somatic nonsynonymous SNVs and small insertions/deletions per Mb in the coding region (with Variant Allele Frequency (VAF) ≥0.03 for tumor tissue, and ≥0.005 for ctDNA, respectively). TMB-high patients were identified with ≥9 mutations/Mb using the top quartile threshold of 2,000 samples from the Geneplus database (9).

Molecular tumor burden index (mTBI) was calculated using the mean allele fraction of mutations in a mutation cluster with the highest cellular prevalence of circulating tumor DNA (ctDNA) at each time point sample. mTBI reflects the percentage of ctDNA detected in cfDNA and its changes can reflect the change of tumor burden at the molecular level (11).

The ctDNA concentration was expressed in haploid genome equivalents (hGE) per mL of plasma (hGE/mL) and was calculated by multiplying the maximal VAF of detected mutations by the plasma cell-free DNA concentration and dividing by 3.3, with the assumption that each haploid genomic equivalent weighs 3.3 pg (12). The ΔctDNA was calculated as the ctDNA concentration at the on-treatment timepoint minus the baseline ctDNA concentration in our cohort. ctDNA clearance was defined as the absence of non-germline variants at the response timepoint.

4–5 μm thickness sections were cut from tumor FFPE tissues. The specific steps are as follows: baking at 63°C for 1h, dewaxing at room temperature, antigen repair is completed in EDTA solution (pH=9.0) and cleaned by TBST buffer for three times. Then the markers (CD4, CD8, PD-1, PD-L1, CD56) was labeled and stained. Upon the completion of the staining, the slides were scanned using Vectra 3.0.5 (PerkinElmer, Massachusetts, USA) in accordance with the manufacturer’s instructions. Multispectral images were unmixed using spectral library built from images of single stained tissues using the inform Advanced Image Analysis software (inForm 2.3.0; PerkinElmer, Massachusetts, USA). 12–15 representative original multispectral images were selected to train the inForm software (tissue segmentation, cell segmentation, phenotyping tool, and positivity score). All the settings applied to the training images were saved within an algorithm for batch analysis of multiple original multispectral images of the same tissue. The positive percentage of CD4, CD8, PD-1, PD-L1 or CD56 were calculated by the number of CD4, CD8, PD-1, PD-L1 or CD56 positive cell divided by the total number of nucleated cells.

In this exploratory analysis of OS, mOS and OS rates were estimated using the Kaplan-Meier method. The Hazard Ratios (HRs) and confidence intervals (CIs) were estimated using stratified Cox proportional hazard models. Patients without events were censored at the date they were last known to be alive. Analyses were performed on the data in the overall trial population. In the current extended follow-up (data cutoff, January 31, 2024), a sufficient number of events was reported to provide an estimate of mOS.

The study enrolled a total of 20 patients between May 2019 and October 2020. As of the analysis cutoff date (January 31, 2024), six patients remained alive, with four maintaining progression-free disease status. Among those who experienced disease progression, progression-free survival was 29.6 and 25.6 months, respectively. Of the patients who remained progression-free, one discontinued treatment owing to adverse events, while another withdrew for personal reasons. The remaining two patients continued to participate in the study treatment. The median follow-up duration for the trial, from the time of consent to the cutoff date, was 49.5 months, with a minimum follow-up of 39.9 months. Patients received a median of 12 cycles of NK cells treatment (range: 2–36 cycles) and 11.5 cycles of Sintilimab (range: 2–64 cycles). The discrepancy in the maximum number of cycles between NK cell infusions and Sintilimab administration was mainly because disruptions during the COVID-19 pandemic sometimes prevented timely NK cell infusions while Sintilimab treatment continued as planned.

The primary endpoint PFS, secondary endpoints included ORR, DoR, DCR, CBR have been reported previously (8). PFS was 11.7 months at the extended follow-up analysis, almost the same as previously reported ( Supplementary Figure S2 ).

The median OS was 27.3 months (95% CI, 0.76 to 53.8) in the intention-to-treat (ITT) population ( Figure 1 ). OS was assessed across various subgroups, including age, smoking status, histology, and TMB level. No statistically significant differences were observed among these subgroups ( Supplementary Table S1 ).

With this extended follow-up analysis, no new safety events were observed ( Supplementary Table S2 ). In conclusion, the adverse events (AEs) encountered by the ITT population corresponded with those delineated in the preceding report of this study. The most common grade 3 AEs observed were hyperglycemia, hypertriglyceridemia, creatine kinase increased, white blood cell decreased, and neutrophil count decreased.

Liquid biopsy samples were collected at baseline and at each subsequent efficacy evaluation.

Among the 20 participants, liquid biopsy samples were obtained from 19 individuals. Of these, 18 patients had ctDNA measurements available for both baseline and post-baseline dynamics analyses. The top three identified mutated genes were TP53, LRP1 and NOTCH1.Notably, 7 out of 9 patients who achieved the best response of partial response (PR), or complete response (CR) exhibited decreased ctDNA levels during therapy. In contrast, 3 out of 4 patients who had the best response of progressive disease (PD) showed increased ctDNA levels at the first efficacy evaluation timepoint ( Figure 2 ).

Through the analysis of PFS and OS in relation to ctDNA dynamics, we determined that changes in ctDNA levels (ΔctDNA) serve as a robust prognostic indicator. Median PFS for patients exhibiting decreased ctDNA levels was 26.2 months versus 3.2 months for patients with increased ctDNA (HR=0.19 [95%CI: 0.06-0.67]; P=0.01) ( Figure 3A ; Supplementary Figure S3A ). Median OS for patients exhibiting decreased ctDNA levels was not achieved (54.5% of patients in this group were still alive) versus 13.6 months for patients with increased ctDNA (HR=0.23 [95%CI: 0.07-0.80]; P=0.021) ( Figures 3A, B ). Additionally, patients who achieved ctDNA clearance had significantly longer PFS (median PFS: 29.6 months vs 3.2 months; HR=0.08 [95%CI: 0.02-0.41]; P =0.002) and OS (median OS: NA vs 14.65 months; HR=0.11 [95%CI:0.03-0.46]; P=0.002) compared to those who never achieved ctDNA clearance ( Figure 3C ; Supplementary Figure S3B ). Of particular note, patients who achieved ctDNA clearance demonstrated significantly prolonged survival outcomes, regardless of the timing of clearance or their best radiographic response, demonstrated significantly prolonged survival outcomes ( Figure 3A ). According to the latest ctDNA analysis, 3 out of 4 progression-free patients continue to show ctDNA clearance. The Response Evaluation Criteria in Solid Tumors (RECIST) were less effective in distinguishing patients with PR from those with stable disease (SD) in terms of PFS and OS, underscoring the heterogeneity of radiographic response assessments. In contrast, ctDNA more accurately reflects tumor burden dynamics. ΔctDNA demonstrated significant value in capturing survival outcomes. Furthermore, ctDNA analysis has revealed alterations in the mutational landscape in patients who develop resistance to combined NK cells and PD-1 monoclonal antibody (mAb) therapy ( Supplementary Figure S4B ).

We previously reported that mTBI closely matches radiographic evaluations (8). Long-term monitoring of mTBI in patients has demonstrated a strong correlation between mTBI dynamics and radiographic changes. Specifically, a decline in mTBI levels is observed when patients achieve a radiographic response. Conversely, an increase in mTBI levels typically precedes tumor progression ( Supplementary Figure S4 ).

We evaluated the expression of PD-1 and PD-L1 on NK cells. By comparing the long-survival group (>24months) to the short-survival group (≤24 months) (6), we found no significant differences in baseline PD-1 and PD-L1 expression on NK cells ( Supplementary Figures S5A, B ). However, patients who exhibited an increased percentage of PD-L1+ NK cells after treatment showed significantly longer OS compared to those with a decreased percentage (median OS: 43.5 months vs 6.0 months; HR=0.07 [95%CI: 0.01-0.39]; P=0.002). In contrast, no significant difference in prognosis was noted between patients with increased versus decreased percentages of PD-1+ NK cells (median OS: 14.8 months vs 32.2 months; HR=1.38 [95%CI: 0.40-4.73]; P=0.613) ( Figures 4A, B ).

Multiplex immunofluorescence analysis was successfully conducted on adequate tissue specimens obtained from eight patients: two with adenocarcinoma and six with squamous carcinoma. In terms of treatment outcomes, two patients exhibited best response of PD, three had SD, two achieved PR, and one attained CR. The median PFS for these patients was 15.8 months, and the median OS was 22.5 months.

The analysis, which involved comparing the proportions of different cell types to patients’ survival, clearly facilitates the differentiation of patients based on their cellular phenotypes ( Figure 5 ). To stratify patients into high and low groups for comparison, cutoff values are established at decile intervals. This approach yields bubble charts for various markers as depicted below, with HR and p-values organized into a table. Folders are created according to different deciles, and for each marker at each cutoff value, Kaplan-Meier survival curves were plotted and annotated with HR and p-values. The HR were calculated using the Cox proportional hazards model, while between-group comparisons were evaluated using the log-rank test ( Figure 6 ). The analysis revealed that the presence of CD56+, CD8+, PD-1+, and PD-L1+ cells within the tumor microenvironment serves as a harbinger of a favorable therapeutic response, demonstrating a significant association with extended patient survival ( Supplementary Figure S6 ). Notably, the CD56+PD-L1+ cellular phenotype emerged as a particularly potent predictor of beneficial treatment outcomes, consistently correlating with increased survival durations ( Supplementary Figures S6 , S7 ).