Maximum Somatic Allele Frequency-Adjusted Blood-Based Tumor Mutational Burden Predicts the Efficacy of Immune Checkpoint Inhibitors in Advanced Non-Small Cell Lung Cancer

Simple Summary Recent studies exhibited the unstable prediction ability of blood-based tumor mutational burden (bTMB) when predicting the response of immune checkpoint inhibitors (ICIs) therapy in patients with non-small cell lung cancer (NSCLC). Circulating tumor DNA (ctDNA) abundance, usually represented by maximum somatic allele frequency (MSAF), was one possible confounding factor influencing bTMB ability in ICIs response prediction. We herein developed a novel approach to optimize the calculation of bTMB by integrating MSAF, namely, MSAF-adjusted bTMB (Ma-bTMB), to better select beneficiaries of ICIs. Our present results showed that this novel non-invasive biomarker could reduce the confounding effect of MSAF and ITH on bTMB calculation and effectively identify beneficiaries of ICIs in patients with advanced NSCLC, warranting future clinical trials. Abstract Introduction: Recent studies exhibited the unstable prediction ability of blood-based tumor mutational burden (bTMB) when predicting the response of immune checkpoint inhibitors (ICIs) therapy in patients with non-small cell lung cancer (NSCLC). Circulating tumor DNA (ctDNA) abundance, usually represented by maximum somatic allele frequency (MSAF), was one possible confounding factor influencing bTMB ability in ICIs response prediction. Methods: MSAF-adjusted bTMB (Ma-bTMB) was established and validated in patients with advanced NSCLC among Geneplus Cancer Genome Database (GCGD, n = 1679), Zhuo (n = 35), Wang (n = 45), POPLAR (NCT01903993, n = 211) and OAK (NCT02008227, n = 642) cohorts. Results: MSAF demonstrated a modest positive correlation with bTMB and a negative one with survival benefit. Improved survival outcomes of ICIs therapy have been observed among patients with high-Ma-bTMB compared to those with low-Ma-bTMB in Zhuo and Wang cohorts. In addition, compared to low-Ma-bTMB, high-Ma-bTMB was associated with more positive clinical benefits from ICIs therapy than chemotherapy both in POPLAR and OAK cohorts. Further exploration suggested that Ma-bTMB could precisely identify more potential ICIs beneficiaries compared to bTMB and LAF-bTMB, complementary to PD-L1 expression. Conclusions: We developed Ma-bTMB, a convenient, readily available, non-invasive predictive biomarker effectively differentiates beneficiaries of ICIs therapy in advanced NSCLC, warranting future clinical trials.


Introduction
Immune checkpoint inhibitors (ICIs), targeting programmed death-ligand 1 (PD-L1) or programmed death 1 (PD-1), have been found to be associated with improved overall survival (OS) of patients with advanced non-small cell lung cancer (NSCLC) [1][2][3]. However, evidence showed that only a few patients could benefit from ICIs treatment, and identifying predictors of favorable ICIs treatment outcomes is of great importance.
Though PD-L1 expression and tissue-based tumor mutational burden (tTMB) have been proven to be useful predictive biomarkers for patients with advanced NSCLC treated with ICIs [4,5], sufficient high-quality tissue samples for standard biomarker analysis were available only for a minority of patients [6]. Hence, developing a non-invasive, convenient, and safe method is urgently needed. Blood-based TMB (bTMB) is of great interest as a candidate in place of tTMB, whereas the results showed that bTMB could not predict OS benefit for ICIs treatment [7,8]. A recently updated result of the phase 2 B-F1RST trial showed that bTMB failed to predict both OS and progression-free survival (PFS) of patients with NSCLC receiving ICIs at a median 36.5-month follow-up [9]. Moreover, Wei Nie et al. found an upside-down U-shaped curve between bTMB and survival, indicating that patients with lower or higher bTMB had longer PFS and OS compared with those with medium bTMB [10]. In short, the information above revealed the unstable prediction ability of bTMB when predicting the survival benefit of ICIs therapy.
Different from tTMB obtained through local tumor biopsy, bTMB is calculated based on circulating tumor DNA (ctDNA) derived from all tumor niches, implying intratumor heterogeneity (ITH) and tumor burden in a patient with advanced cancer [11,12]. ctDNA abundance, which is difficult to directly measure and usually represented by maximum somatic allele frequency (MSAF), has been thought to estimate tumor content, displaying a positive correlation with bTMB while negatively correlating with survival [7,13]. In addition, ITH refers to distinct subclones or subclonal mutations resulting from tumor evolution within the same tumor-bearing patient, showing a correlation with poor prognosis [14,15]. Hence, there is an inner connection between ctDNA abundance (MSAF as an alternative) and ITH as poor prognosis confounding factors in bTMB prediction for ICIs efficacy in patients with advanced NSCLC.
With the development of targeted genomic sequencing assays, a simple, steady, and propagable modified bTMB algorithm that could identify both OS and PFS beneficiaries from ICIs treatment across various panels and calling principles is needed. In this study, we identified MSAF as the core confounding factor of bTMB in survival prediction of ICIs therapy, according to which we established and validated a novel predictive biomarker MSAF-adjusted bTMB (Ma-bTMB) in patients with advanced NSCLC among Zhuo, Wang, POPLAR, and OAK cohorts. Finally, the application of Ma-bTMB was compared with that of other bTMB biomarkers and PD-L1 expression.

Patient Enrollment and Sample Collection
Eighty patients with stage III-IV NSCLC were enrolled from two centers, Peking University Cancer Hospital & Institute (Zhuo cohort, n = 35) and Cancer Hospital of Chinese Academy of Medical Sciences & Peking Union Medical College (Wang cohort, n = 45) from 25 August 2016 to 12 April 2019. Patients with MSAF = 0 were excluded due to incalculableness (Zhuo cohort, n = 4; Wang cohort, n = 5). All patients were pathologically diagnosed and underwent anti-PD-(L)-1 therapy. The baseline characteristics, mutational information, and clinical information of Wang and Zhuo cohorts are presented in Tables S1-S3. Peripheral blood samples were collected for 1021-gene-panel genomic analysis before ICI administration. The study was designed and conducted in accordance with the Declaration of Helsinki and approved by the Institutional Review Board (IRB) of Cancer Hospital of Chinese Academy of Medical Sciences & Peking Union Medical College (NCC2018-092) and Peking University Cancer Hospital & Institute (2017KT57). Written informed consent was granted prior to sample collection, genetic sequencing, and data analysis from each patient.

Datasets and Cohorts
To analyze the association between MSAF and bTMB, somatic single nucleotide variations (SNVs) and indels detected in the ctDNA by 1021-gene-panel from the plasma of 1679 Chinese NSCLC patients enrolled from 1 July 2018 to 12 March 2020, were downloaded from Geneplus Cancer Genome Database (GCGD) cohort (n = 1679, accessed in October 2020), which is an internal database of Geneplus Co., Ltd. (Bangkok, Thailand). To verify that Ma-bTMB is associated with ICIs rather than chemotherapy response, POPLAR (NCT01903993) and OAK (NCT02008227) cohorts, two large randomized clinical trials investigating the efficacy and safety of atezolizumab (MPDL3280A, anti-PD-L1) compared with docetaxel (chemotherapy) in patients with NSCLC, were obtained from a published study and included in the subsequent algorithm analysis [7]. Different from the 1021-genepanel, baseline bTMB in these two public cohorts was sequenced through FDA-approved FoundationOne (F1) CDx NGS assay (POPLAR, N = 211; OAK, n = 642). Among these patients, 106 received docetaxel therapy and 105 received atezolizumab therapy in the POPLAR cohort, with 318 vs. 324, respectively, in the OAK cohort [7]. Detailed mutational and clinical information on these two cohorts is presented in Table S4. PD-L1 expression was only available in OAK cohort; PD-L1-positive expression was defined as PD-L1 expression on ≥1% of tumor cells (TCs) or ≥1% of tumor-infiltrating immune cells (ICs).

Targeted Capture Sequencing and Genomic Data Analysis
Blood samples were collected in EDTA vacutainer tubes (BD Diagnostics, Franklin Lakes, NJ, USA) and subjected to ctDNA extraction and sequencing within 3 hours. Peripheral blood lymphocytes and plasma were separated via sequential centrifugation (2500× g, 10 min; 16,000× g, 10 min). Next, circulating free DNA (cfDNA) was extracted using a DNeasy Blood Kit (Qiagen, Hilden, Germany) and QIAamp Circulating Nucleic Acid Kit (Qiagen). White blood cells from the same patient were used as control. Sequencing libraries were constructed utilizing the KAPA DNA Library Preparation Kit (Kapa Biosystems, Wilmington, MA, USA) per the manufacturer's instructions. Barcoded libraries were subsequently hybridized to a customized panel of 1021 genes, consisting of whole exons and selected introns of 288 genes, along with selected regions of 733 genes. A detailed gene list was described in our previous study [16]. DNA sequencing was performed with the HiSeq 3000 Sequencing System (Illumina, San Diego, CA, USA) with 2 × 101 bp paired-end reads. For each mutation, somatic allele frequency was calculated by dividing the number of mapped mutational reads by the average sequencing depth. The number of somatic coding non-synonymous SNVs, small insertions, and deletions per megabase (Muts/Mb) of each genome examined was calculated as bTMB, and germline polymorphisms were removed for calculating bTMB. Then, MSAF-adjusted bTMB was calculated as Ma-bTMB = bTMB / (MSAF × 10). The number 10 is a manually selected constant, which could scale down the values of Ma-bTMB to roughly fall within a range similar to other reported bTMB values.

Subclone and Subclonal Mutation Calculation
PyClone was employed to estimate clonal and subclonal populations [17]. We defined a subclone as a collection of cells in the tumor sample that harbored the same set of genomic variants, including SNVs, insertions, and deletions. Subclonal mutations were the total mutated genes of corresponding subclone clusters.

Statistics
Efficacy evaluation was performed by using RECIST version 1.1, and PFS was defined as the time from the start of anti-PD-(L)1 treatment until either objective disease progression (assessed by an investigator using RECIST version 1.1) or death from any cause. OS was defined as the time from the start of anti-PD-(L)1 treatment until death from any cause. We censored patients who had not progressed, who had died at the time of statistical analysis, or who failed to follow up but had not died at the time of their last evaluation.
Pearson correlation tests were employed to determine the association between MSAF and OS/PFS or bTMB. Kaplan-Meier curves were plotted and, together with log-rank tests, performed to calculate the hazard ratios (HRs) and evaluate the prognostic value of various markers. The interaction p-value was obtained from an unstratified proportional Cox model, including the terms of treatment, the biomarker subgroup, and treatment by the biomarker subgroup interaction. For all analyses, p < 0.05 was considered statistically significant. IBM SPSS software (version 24.0), GraphPad Prism (version 9.0), and R 4.2.0 were used for the statistical analysis.

Study Design and Cohorts
The study design is illustrated in Figure 1. The advanced NSCLC patient is characterized by high tumor burden and ITH. cfDNA was obtained for ITH, bTMB, and MSAF analysis, according to which Ma-bTMB was established. Patients with MSAF = 0 were excluded due to incalculableness. Finally, the algorithm Ma-bTMB was explored and validated in 71 patients with advanced NSCLC receiving anti-PD-(L)-1 from Zhuo (discovery cohort, n = 31) and Wang (validation cohort, n = 40) cohorts and 851 patients with advanced NSCLC receiving anti-PD-(L)-1 or chemotherapy from POPLAR (discovery cohort, n = 209) and OAK (validation cohort, n = 642) cohorts. PD-L1 expression was only available in the OAK cohort (n = 635).

bTMB, MSAF, and ITH Analysis
In the Zhuo cohort as the discovery cohort, we found that bTMB could not identify either PFS or OS benefits from ICIs ( Figure S1); similar results were observed in Wang cohort. Notably, MSAF was found a significant correlation with poor OS (R = −0.42, p < 0.001) and PFS (R = −0.28, p = 0.017; Figure S2). Moreover, there was a positive correlation between MSAF and bTMB (R = 0.41, p < 0.001) as well as between MSAF and subclones (R = 0.37, p < 0.001) in GCGD cohort (Figures 2A and S3), indicating that MSAF might be the core confounding factor in the bTMB prediction of ICIs effect. Therefore, we developed the modified biomarker Ma-bTMB by dividing a patient's bTMB value with the corresponding MSAF and bTMB ( Figure 1).
To explore whether Ma-bTMB balances the poor prognosis effect of ITH, we compared subclones and subclonal mutations between high and low groups divided by bTMB and Ma-bTMB. A 25% increase in the bTMB value with a cutoff score of 14 was defined as the cutoff point; the Ma-bTMB cutoff score of 10 was selected as the HR and p-value reaching a minimum in PFS prediction compared with other cutoffs ( Figure 2B). Even though no significant differences in subclones numbers were observed between high-bTMB and low-bTMB groups (Mean: 25.3 vs. 26.4, p = 0.877), there was a higher level of subclonal mutations in the high-bTMB group than in the low-bTMB group (Mean: 234.2 vs. 150.7, p = 0.032) in Wang and Zhuo cohorts (hereafter referred to as Wang & Zhuo cohort; Figure 2C). As for Ma-bTMB, no difference was observed between high-and low-Ma-bTMB groups regardless of subclones (Mean: 11.6 vs. 13.0, p = 0.710) or subclonal mutation level (Mean: 43.9 vs. 70.1, p = 0.113), suggesting that Ma-bTMB could adjust bTMB by balancing the effect of ITH on stratification.

bTMB, MSAF, and ITH Analysis
In the Zhuo cohort as the discovery cohort, we found that bTMB could not identify either PFS or OS benefits from ICIs ( Figure S1); similar results were observed in Wang cohort. Notably, MSAF was found a significant correlation with poor OS (R = −0.42, p < 0.001) and PFS (R = −0.28, p = 0.017; Figure S2). Moreover, there was a positive correlation between MSAF and bTMB (R = 0.41, p < 0.001) as well as between MSAF and subclones (R  = 0.37, p < 0.001) in GCGD cohort (Figure 2A, Figure S3), indicating that MSAF might be the core confounding factor in the bTMB prediction of ICIs effect. Therefore, we developed the modified biomarker Ma-bTMB by dividing a patient's bTMB value with the corresponding MSAF and bTMB ( Figure 1). To explore whether Ma-bTMB balances the poor prognosis effect of ITH, we compared subclones and subclonal mutations between high and low groups divided by bTMB and Ma-bTMB. A 25% increase in the bTMB value with a cutoff score of 14 was defined as the cutoff point; the Ma-bTMB cutoff score of 10 was selected as the HR and p-value reaching a minimum in PFS prediction compared with other cutoffs ( Figure 2B). Even though no significant differences in subclones numbers were observed between high-bTMB and low-bTMB groups (Mean: 25.3 vs. 26.4, p = 0.877), there was a higher level of subclonal mutations in the high-bTMB group than in the low-bTMB group (Mean: 234.2 vs. 150.7, p = 0.032) in Wang and Zhuo cohorts (hereafter referred to as Wang & Zhuo cohort; Figure  2C). As for Ma-bTMB, no difference was observed between high-and low-Ma-bTMB

Discussion
In this study, we discovered that high-Ma-bTMB could predict the PFS and OS benefit of anti-PD-(L)1 therapy while serving as an effective predictive biomarker for selecting patients who could benefit from atezolizumab rather than chemotherapy. Patients with high-Ma-bTMB showed approximately 11 months longer OS in Zhuo (17.6 vs. 6.6 months) and Wang (14.4 vs. 3.6 months) cohorts after receiving ICIs, as well as about 8 months improved OS in POPLAR (15.6 vs. 7.8 months) and OAK (16.2 vs. 8.3 months) cohorts after receiving ICIs compared to chemotherapy. Even in patients with PD-L1 negative expression, those with high-Ma-bTMB receiving atezolizumab still had longer OS than chemotherapy. Our study showed that Ma-bTMB could steadily select ICIs beneficial populations in advanced NSCLC.
The establishment of Ma-bTMB is needed. With the revelation of subsequent clinical trials with long follow-ups, such as Keynote-001 and CheckMate-067, etc., it becomes clearer that ICIs therapy provides long-term benefits over conventional treatments, and identifying OS benefit is of great importance [16][17][18]. Although bTMB, whether tested by whole-exome sequencing or targeted gene panels, is deemed as a non-invasive biomarker candidate predicting the clinical benefit of ICIs therapy [7,8,19], recent evidence showed that bTMB failed to stably predict the survival benefit of ICIs therapy, particularly OS [9,10,13,20].
The establishment of Ma-bTMB is reasonable. Distinct from tTMB, the analysis of bTMB had been found to be impacted by several factors in the blood, including ctDNA, AF, MSAF grouping, ITH, etc. [13,[20][21][22][23][24]. Considering that a patient with advanced NSCLC is characterized by high tumor burden and ITH, the relationship among bTMB, MSAF, and ITH has been analyzed. Consistent with previous studies [7,13,22], MSAF displayed a positive correlation with bTMB while negatively correlating with the survival of patients with advanced NSCLC receiving ICIs. Moreover, we found that MSAF was positively correlated with ITH. Hence, MSAF was hypothesized as the core confounding factor of bTMB in the prediction ability of ICIs therapy. In addition, previous evidence also indicated that MSAF alone was not an effective predictive biomarker of ICIs response [13]; hence, Ma-bTMB was established.
Genetically distinct cellular populations, resulting from clonal expansions, spatial segregation, and incomplete selective sweeps, were manifested as ITH [25]. Moreover, nearly all tumors display subclonal expansions, even at limited read depth [26]. Schmitt et al. evaluated the clonality of actionable driver mutations, reasoning that targeting subclonal mutations will likely result in treatment failure [27]. Considering the different numbers of subclones and subclone mutations, which were manifested as ITH, may also affect our evaluation of the subsequent treatment outcome; we analyzed their difference in different levels of Ma-bTMB. In our study, results showed that Ma-bTMB could decrease the effect of ITH on the stratification of beneficiaries and non-beneficiaries compared to bTMB. Different from the high-bTMB group presenting higher levels of subclonal mutations over the low-bTMB group, there were no differences between the high-Ma-bTMB and low-Ma-bTMB groups. In addition, this result is consistent with the evidence of the negative impact of ITH on PD-(L)1 inhibition therapy [28][29][30]. In other words, Ma-bTMB could decrease the effect of ITH on the stratification of beneficiaries and non-beneficiaries compared to bTMB, making it useful for predicting ICI-related benefits.
Although PD-L1 expression, tested via immunohistochemistry assays, is currently considered the most widely accepted biomarker guiding the selection of patients to receive anti-PD-L1/PD-1 treatment [31], patients with low PD-L1 expression could also benefit from anti-PD-L1/PD-1 treatment [4,32]. Indeed, the interaction P for OS between Ma-bTMB and treatment subgroups was significant in patients with PD-L1 negative expression and not significant in patients with PD-L1 positive expression. Patients with PD-L1 positive expression and high-Ma-bTMB showed the longest OS benefit (17.3 months) among all, and those with high-Ma-bTMB had 9 months longer OS after ICI therapy over chemotherapy (15.9 vs. 6.3 months) in patients with PD-L1 negative expression. Intriguingly, regardless of PD-L1 expression, patients with low-Ma-bTMB receiving ICIs therapy showed neither PFS nor OS benefit compared to chemotherapy. In this scenario, a combination biomarker of PD-L1 expression and Ma-bTMB analysis could guide the treatment selection between ICIs and chemotherapy.
Previously, we developed an adjusted bTMB algorithm, LAF-bTMB, by removing high AF mutations to predict the clinical benefit of PD-(L)1 inhibitors in patients with advanced NSCLC [13]. In this study, we compared Ma-bTMB, bTMB, and LAF-bTMB in POPLAR&OAK and Zhuo & Wang cohorts. Unfortunately, bTMB could not effectively identify the OS benefit of ICIs therapy in either cohort. Compared to LAF-bTMB, Ma-bTMB could identify more potential OS and PFS beneficiaries in both cohorts. The value distributions of three biomarkers showed that LAF-bTMB had the highest and narrowest "slim" curve while Ma-bTMB had the lowest and widest "flat" curve, indicating that Ma-bTMB could better disperse the population and reduce the accumulation of patients according to the characteristic of bTMB. Given that patients with higher Ma-bTMB/bTMB/LAF-bTMB have better ICIs response, the "flat" curve may decrease the selection bias caused by the definition of the cutoff value. Hence, Ma-bTMB could steadily identify more beneficiaries to minimize potential patients who miss the chance of ICIs benefits.
There are several limitations to this study. First, the algorithm was developed retrospectively, and the sample size of the Wang and Zhuo cohorts was relatively small, which could bias the results. In addition, the detection platforms and methodologies differed in POPLAR&OAK and Wang & Zhuo cohorts. A larger prospective study is warranted for further validation. Nevertheless, the stable prediction ability of Ma-bTMB between these two platforms (1021-gene-panel and F1CDx NGS assay) and four cohorts (Zhuo, Wang, POPLAR, and OAK) indicated its feasibility, to some extent, as a predictive biomarker of ICI benefit in patients with advanced NSCLC.

Conclusions
In summary, although a larger prospective study is necessary for further verification, we developed and validated a convenient and readily available novel algorithm, Ma-bTMB, that can predict clinical benefits to guide ICI therapy in patients with advanced NSCLC.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/cancers14225649/s1, Figure S1:The forest plots for HRs and p-values of OS and PFS comparing bTMB-H and bTMB-L subgroups with corresponding cutoff points in Zhuo cohort and Wang cohort; Figure S2: The correlation between MSAF and OS/PFS in Wang & Zhuo cohort; Figure S3: The distribution of subclones between high MSAF and low MSAF changed with MSAF cutoff points in GCGD cohort. Figure S4: The forest plots for HRs and p values of PFS between atezolizumab and docetaxel comparing Ma-bTMB-H and Ma-bTMB-L subgroups with corresponding cutoff points in POPLAR cohort. Figure S5: The forest plots for HRs and p-values of OS comparing LAF-bTMB-H and LAF-bTMB-L changed with corresponding cutoff points in Wang & Zhuo cohort. Table S1: Basic characteristic of patients in Wang cohort; Table S2: Basic characteristic of patients in Zhuo cohort; Table S3: Clinical data and somatic mutations of Wang & Zhuo cohort; Table S4: Clinical data and somatic mutations of OAK & POPLAR cohort.
Author Contributions: Z.W. and J.W. developed the concept and designed the study. All authors acquired, analyzed, or interpreted of data. Y.Z., Y.Y., Y.G. and S.H. performed the statistical analysis. Y.D., Y.X., J.Y., M.Z. and Z.Y. wrote the manuscript. Z.W., M.Z., X.X. and X.Y. performed a critical revision of the manuscript for the intellectual content. All authors have read and agreed to the published version of the manuscript.