Gene expression signatures as candidate biomarkers of response to PD-1 blockade in non-small cell lung cancers.

Although anti-PD-1/PD-L1 monotherapy has achieved clinical success in non-small cell lung cancer (NSCLC), definitive predictive biomarkers remain to be elucidated. In this study, we performed whole-transcriptome sequencing of pretreatment tumor tissue samples and pretreatment and on-treatment whole blood samples (WB) samples obtained from a clinically annotated cohort of NSCLC patients (n = 40) treated with nivolumab (anti-PD-1) monotherapy. Using a single-sample gene set enrichment scoring method, we found that the tumors of responders with lung adenocarcinoma (LUAD, n = 20) are inherently immunogenic to promote antitumor immunity, whereas those with lung squamous cell carcinoma (LUSC, n = 18) have a less immunosuppressive tumor microenvironment. These findings suggested that nivolumab may function as a molecular targeted agent in LUAD and as an immunomodulating agent in LUSC. In addition, our study explains why the reliability of PD-L1 expression on tumor cells as a predictive biomarker for the response to nivolumab monotherapy is quite different between LUAD and LUSC.

Introduction

For the last several decades, there have been remarkable advances in cancer immunology and cancer immunotherapy. The clinical efficacy of immuno-oncology (IO) agents that inhibit the programmed death 1 (PD-1)-programmed death ligand 1 (PD-L1) signaling pathway in advanced NSCLC has been described. In particular, nivolumab, a fully human IgG4 anti-PD-1 monoclonal antibody, is the first approved IO agent for use in patients with previously treated advanced NSCLC. Unfortunately, nivolumab monotherapy is not effective in all patients with advanced NSCLC, with an objective response rate (ORR) of no higher than 20% [1, 2].

Preclinical and clinical studies have revealed that various factors can affect the clinical outcomes of anti-PD-1/PD-L1 monotherapy in patients with advanced NSCLC; these factors include PD-L1 expression, the presence of tumor-infiltrating lymphocytes (TILs), the tumor mutation burden (TMB), human leukocyte antigen class I (HLA-I) genotype, T-cell repertoire diversity, the gene expression profile and the gut microbiota [3-7]. Among these factors, PD-L1 expression, which is determined by immunohistochemical assays, is the only biomarker currently approved as a companion or complementary diagnostic biomarker for the response to anti-PD-1/PD-L1 monotherapy. In general, PD-L1 expression on tumor cells is positively correlated with the clinical response to anti-PD-1/PD-L1 monotherapy in patients with nonsquamous NSCLC [2, 8]. Importantly, however, some patients with a positive PD-L1 expression rate of less than 1% also benefit from anti-PD-1/PD-L1 monotherapy, with an ORR of approximately 10% [9]. Unlike in patients with nonsquamous NSCLC, PD-L1 expression has no predictive value in patients with squamous NSCLC (1). PD-L1 expression alone remains an imperfect biomarker, and none of the abovementioned factors has yet proven to be sufficiently robust for clinical use. Therefore, extensive efforts have been devoted to exploring predictive biomarkers for the response to anti-PD-1/PD-L1 monotherapy.

The recent development of next-generation sequencing (NGS) technology and computational genomic tools has enabled deep analysis of the so-called omics data, including genomic, transcriptomic, proteomic, and epigenomic data [10]. To date, biomarker studies in NSCLC have mainly focused on the histological or genomic analyses, and there have been very few transcriptomic analyses. In metastatic melanoma and other types of cancer, multiple comprehensive transcriptomic studies have identified immune-related gene expression signatures that are positively or negatively associated with the clinical response to anti-PD-1 and anti-CTLA4 monotherapy [11, 12]. We thus assumed that by combining gene expression signatures obtained from pretreatment and on-treatment transcriptomic profiles with clinical data of patients, we could identify a novel promising biomarker to predict the response to anti-PD-1/PD-L1 monotherapy in NSCLC patients.

The two main histological subtypes of NSCLC are adenocarcinoma and squamous cell carcinoma. In biomarker studies, these subtypes have often been combined and analyzed in the same way; however, each exhibits distinct mutational and genomic profiles. In this cohort, we analyzed the patient’s transcriptomic features separately to identify histology-specific gene expression signatures that are associated with the clinical response to nivolumab monotherapy. Characterization of these signatures will help us to decipher the complexity of tumor-immune interactions and deepen the understanding of the tumor microenvironment (TME) that favors a better clinical response to nivolumab monotherapy.

Materials and methods

Ethical statement

All clinical data and patient samples were collected following approval by the Sendai Kousei Hospital Institutional Review Board (IRB) (IRB number: 29–4). The study period is from 18 Jul 2017 to 31 Dec 2020. In all cases, written informed consent was obtained from the patients.

Patient characteristics and sample collection

A total of 40 patients with advanced NSCLC were enrolled in this cohort. All enrolled patients were administered at least one dose of nivolumab. Tumor tissue samples were collected before the first dose of nivolumab (pretreatment tumor tissue). Immediately after the biopsy procedure (three tumor tissue samples of approximately 1.5−3.0 mm in diameter per patient), the obtained tumor tissues were suspended in RNAlater RNA Stabilization Reagent (QIAGEN, Hilden, Germany) and stored at −80°C for batched RNA extraction. Whole blood (WB) samples were collected before the first dose and after the fourth or fifth dose of nivolumab (pretreatment WB and on-treatment WB, respectively). To obtain whole blood samples, 2.5 mL of blood was collected into PAXgene Blood RNA Tubes (BD, Franklin Lakes, NJ) and stored at −80°C for batched RNA extraction.

For 28 of the 40 patients (15 LUAD and 13 LUSC), pretreatment tumor tissue samples were available and used in subsequent analyses; the other 12 patients were excluded because of an inadequate quantity of tumor cells. Pretreatment WB was obtained from 39 patients (20 LUAD and 17 LUSC); for one patient, isolation of WB was unsuccessful due to a blood sampling error. On-treatment WB was obtained from 32 patients (15 LUAD and 15 LUSC); for the other 8 patients, blood sampling was not performed due to death, early disease progression or early treatment discontinuation. In this cohort, we employed progression-free survival (PFS) time as an outcome measure of treatment efficacy. The PFS time was defined as the length of time from the start of nivolumab monotherapy until progression. Patients with a PFS time equal to or more than 6 months were defined as responders; those with a PFS time less than 6 months, non-responders.

The PD-L1 tumor proportion score (TPS) was determined by a commercial PD-L1 IHC assay with PD-L1 IHC 22C3 pharmDx (Dako, Carpinteria, CA). During the enrollment and follow-up period of this cohort, four genetic screens of driver mutations, including EGFR and BRAF mutations and ALK and ROS1 fusions, had been approved and were commercially available for advanced NSCLC in Japan. Other genetic aberrations, such as KRAS mutations, were screened in the LC-SCRUM-Asia (formerly LC-SCRUM-Japan) consortium [13], which employed an amplicon-based next-generation sequencing (NGS) panel, Oncomine Comprehensive Assay (Thermo Fisher Scientific, Waltham, MA). Tumor responses to nivolumab monotherapy were assessed by the investigators according to RECIST (Response Evaluation Criteria in Solid Tumors) version 1.1.

RNA extraction and whole transcriptome sequencing (RNA-seq)

Total RNA was extracted from tumor tissue and WB samples using an RNeasy Mini Kit (QIAGEN) and a PAXgene Blood RNA Kit (QIAGEN), respectively, according to the manufacturer’s protocol. Homogenization of tumor tissue samples was carried out using a QIAshredder homogenizer (QIAGEN). DNase I was used during processing with the manufacturer’s protocol to minimize DNA contamination.

The RNA-seq libraries were prepared from the total RNA extracts using an Illumina TruSeq Stranded mRNA Library Kit (Illumina, San Diego, CA), following the manufacturer’s protocol. Paired-end 100 bp sequencing of these libraries was performed on an Illumina HiSeq3000 platform. We utilized HISAT2 [14] to align the RNA-seq reads to the human genome assembly GRCh38. The raw read counts were generated with featureCounts in the Rsubreads Bioconductor package (version 2.4.2) [15] and normalized by conversion to TPM (transcripts per million) [16].

Identification of differentially expressed genes (DEGs)

The RNA-seq raw read counts were employed for differential expression analysis. As a preprocessing step, filtering was applied to exclude genes with low expression, retaining genes with a raw read counts > 1 in at least n samples (where n is the number of samples in each analysis). We fitted a generalized linear model to the preprocessed count data, which is generally assumed to follow a negative binomial distribution, using the DESeq2 Bioconductor package (version 1.29.6) [17]. The statistical significance of differences between responders and nonresponders was assessed with the Wald test. DEGs were defined as genes with |Log2(fold change)| ≥ 1 and an adjusted p-value < 0.1. MA and volcano plots were generated using the ggplot2 R package (version 3.3.2). Hierarchical clustering and heatmap representation of the DEGs were implemented using the ComplexHeatmap Bioconductor package (version 2.5.3) [18]. DESeq2-normalized counts were converted to log10 values, followed by normalization to the z-score values for all genes to reduce expression variance for the genes expressed at different levels.

Fast gene set enrichment analysis (FGSEA)

In classical GSEA [19, 20], correction for multiple hypothesis testing is performed using permutation tests, where independent random gene sets are generated for each permutation and each gene set. Accordingly, standard GSEA is relatively slow because of the huge computational burden imposed by the permutation tests. FGSEA with the fgsea Bioconductor package (version 1.15.0) reuses sampling for different query gene sets to reduce the computational burden, which enables quick and accurate performance of enrichment analysis [21]. As a preprocessing step, we first generated a preranked gene list from the DESeq2 results. This preranked gene list was used for FGSEA with 7,573 gene sets in the Gene Ontology (GO) biological process ontology from MSigDB v7.2 (http://www.gsea-msigdb.org/gsea/msigdb/). The minimal and maximal thresholds of the gene sets were set to 10 and 500 genes, respectively. Statistical significance was calculated with 1,000 permutation tests. Significantly enriched gene sets were defined as those with an adjusted p-value < 0.1. Lollipop plots of the normalized enrichment scores (NESs) obtained from FGSEA were constructed using ggplot2. Running enrichment scores were plotted using the fgsea package. The enrichment maps were generated using Cytoscape App (version 3.8.2: https://cytoscape.org).

Singscore

The TPM normalized counts were used in single-sample gene set scoring with the singscore Bioconductor package (version 1.9.0) [22, 23]. Genes with low counts were filtered based on the CPM (counts per million) value in the edgeR Bioconductor package (version 3.31.4) [24] to avoid rank duplication: we retained genes with a CPM > 1 across more than 50% of the samples. Genes were ranked based on count in increasing order. For each gene set of interest, the mean rank was calculated and normalized to the theoretical minimum and maximum values, centered on zero and then summed to provide a single-sample enrichment score, which ranged between −1 and 1. As with FGSEA, single-sample enrichment scores were generated using the singscore method for gene sets in the GO biological process ontology from MSigDB.

The distributions of the single-sample enrichment scores across all the gene sets were visualized by stacked histograms. The single-sample enrichment scores of the gene sets of interest were shown by scatter plots or box plots and statistically evaluated with the Wilcoxon rank-sum test. In box plots, the lower and upper box hinges indicate the 25th and 75th percentiles, respectively; the central bold line indicates the median; and the whiskers extend to the largest and smallest scores within no more than 1.5× the interquartile range. All the statistical tests were performed using the R program (version 4.0.2; https://www.r-project.org/). Stacked histograms and box plots were generated with ggplot2.

Correlation analysis and regression model fitting using cubic regression splines

Correlation coefficients and p-values between PFS time and the single-sample enrichment scores of the gene sets of interest were calculated by Spearman rank correlation analysis using R. The correlation matrix obtained was visualized with ggplot2.

Nonparametric regression models using cubic regression splines with the R and mgcv packages (version 1.8.31) were fitted to the calculated relationship between PFS time and the single-sample enrichment scores of the gene sets of interest, from which we estimated the predicted PFS time from the scores [25]. Assuming that Log[observed PFS] is normally distributed, we employed the Gaussian family and identity link function. In addition, we used REML (restricted maximum likelihood) as a smoothing parameter estimation method and performed model selection based on AICc (second-order Akaike information criterion) with the MuMIn package (version 1.43.17) [26]. The regression splines were plotted using ggplot2.

Survival analysis

Bubble plots displaying the relationships between PFS time and two of the gene sets of interest were generated using ggplot2. All Kaplan-Meier curves were visualized with the survminer package (version 0.4.9). With the survival package (version 3.2.3), Kaplan-Meier estimates of PFS time for two independent groups were assessed using the two-sided log-rank test, and those for four independent groups were assessed using the two-sided Wald test based on the multivariable Cox proportional hazards model. When comparing the PFS times of patients with high scores to those of patients with low scores, we set the threshold as the median value of the score. Both the waterfall plots describing changes in tumor size and swimmer plots showing patient responses were generated using ggplot2. Cox proportional hazard models were built using the “coxph” function from the survival package and visualized as forest plots using ggplot2.

Results and discussion

From clinically annotated NSCLC patients treated with nivolumab monotherapy in the second- or later-line settings, we prospectively collected tumor tissue and whole blood (WB) samples before the first dose of nivolumab (pretreatment tumor tissues and WB), and WB samples after the fourth or fifth dose of nivolumab (on-treatment WB) (S1 Fig). Patient characteristics and clinical courses are summarized in S2 and S3 Figs, Table 1 and S1 Table. All tumor tissue and WB samples obtained were subjected to whole transcriptome sequencing (RNA-seq). We extracted transcriptomic datasets of LUAD and LUSC from the results, analyzed each histological subtype separately, and explored whether transcriptomic features could differentiate between responders and nonresponders to nivolumab monotherapy.

Table 1

Baseline characteristics of patients in this cohort.

Characteristic (n = 40)		no. (%)
Age	Median (Range)—yr	69 (52–92)
Sex	Male	28 (70.0)
	Female	12 (30.0)
ECOG performance status	0	21 (52.5)
	1	15 (37.5)
	2	3 (7.5)
	3	2 (2.5)
Brinkman index	0 (Never smoker)	11 (27.5)
	200–599	6 (15.0)
	600–1,199	18 (45.0)
	≥1,200	5 (12.5)
Stage of disease*	Recurrent	13 (32.5)
	III	5 (12.5)
	IV	22 (55.0)
Histological subtype	Adenocarcinoma	20 (50.0)
	Squamous cell carcinoma	18 (45.0)
	Adenosquamous carcinoma	1 (2.5)
	Large cell neuroendocrine carcinoma	1 (2.5)
PD-L1 TPS	<1%	15 (37.5)
	1–49%	17 (42.5)
	≥50%	3 (7.5)
No. of prior systemic regimens	1	26 (65.0)
	2	6 (15.0)
	3	4 (10.0)
	4	4 (10.0)
Gene alteration status	EGFR mutation	6 (15.0)
	ALK fusion	1 (2.5)
	KRAS mutation	1 (2.5)
	ERBB2 mutation	1 (2.5)
	PIK3CA mutation	1 (2.5)
	WT or Unknown	30 (75.0)
Prior immuno-oncology therapy	Anti-PD-1/PD-L1**	4 (10.0)
	No	36 (90.0)
Best response to prior systemic regimen	Complete response	0 (0.0)
	Partial response	20 (50.0)
	Stable disease	14 (35.0)
	Progressive disease	6 (15.0)
CNS metastasis	Yes	11 (27.5)
	No	29 (72.5)
Treatment for CNS metastasis	Radiosurgery	9 (81.8)
	Whole brain radiotherapy	0 (0.0)
	No	2 (18.2)
Prior radiotherapy	Yes	21 (52.5)
	No	19 (47.5)

* Clinical staging was conducted according to the 8th Edition TNM Classification for Lung Cancer.

** One patient underwent anti-PD-1 therapy plus chemotherapy; one, anti-PD-L1 therapy plus chemotherapy; two, anti-PD-L1 therapy as postoperative adjuvant therapy.

A previous study has provided a metric for immune cytolytic activity based on gene expression in tumors, where immune cytolytic activity was estimated by the average expression level of CD8A, CD8B, GZMA, GZMB and PRF [27]. Using this metric, we assessed immune cytolytic activity in tumor tissues from patients in our cohort. As the result, we found no significant association between the immune cytolytic activity and clinical outcomes (S4 Fig).

A cohort of patients with LUAD

In the LUAD cohort (n = 20), differential expression analysis between responders and nonresponders was performed using the RNA-seq datasets of pretreatment tumor tissues (n = 15), pretreatment WB (n = 20) and on-treatment WB (n = 15). A total of 15, 68 and 160 differentially expressed genes (DEGs) were identified in pretreatment tumor tissues, pretreatment WB and on-treatment WB, respectively. Of the 15 DEGs in pretreatment tumor tissues, six genes were upregulated and nine genes were downregulated in responders (n = 4) compared to nonresponders (n = 11) (S5A and S5B Fig and S2 Table). Of the 68 DEGs in pretreatment WB, 27 genes were upregulated and 41 genes were down-regulated in responders (n = 5) compared to nonresponders (n = 15) (S5C and S5D Fig and S3 Table). In addition, of the 160 DEGs in on-treatment WB, 20 genes were upregulated and 140 genes were down-regulated in responders (n = 4) compared to nonresponders (n = 11) (S5E and S5F Fig and S4 Table). Heatmaps of DEGs identified in the three datasets demonstrate the relationships between the expression levels of DEGs and clinical factors, including progression-free survival (PFS), PD-L1 tumor proportion score (TPS), Brinkman index and driver mutation status (S6 Fig).

To identify biological processes associated with clinical outcomes in LUAD patients treated with nivolumab monotherapy, we performed gene set enrichment analysis (GSEA) (Fig 1, S7 Fig and S5–S7 Tables). In pretreatment tumor tissues of responders, gene sets related to B cell function (e.g., ‘B CELL-MEDIATED IMMUNITY’ [p = 8.424 × 10−20, NES = 2.738], ‘B CELL RECEPTOR SIGNALING PATHWAY’ [p = 1.692 × 10−11, NES = 2.625], and ‘POSITIVE REGULATION OF B CELL ACTIVATION’ [p = 3.990 × 10−7, NES = 2.225]) and humoral immunity (e.g., ‘HUMORAL IMMUNE RESPONSE MEDIATED BY CIRCULATING IMMUNOGLOBULIN’ [p = 1.534 × 10−20, NES = 3.003], ‘COMPLEMENT ACTIVATION’ [p = 1.317 × 10−19, NES = 2.896], REGULATION OF HUMORAL IMMUNE RESPONSE’ [p = 1.234 × 10−15, NES = 2.853], ‘HUMORAL IMMUNE RESPONSE’ [p = 1.648 × 10−23, NES = 2.664], and ‘IMMUNOGLOBULIN PRODUCTION’ [p = 4.344 × 10−16, NES = 2.626]) were significantly enriched (Fig 1A and S5 Table). Recently, three studies demonstrated that B cells and tertiary lymphoid structures (TLSs) in the TME are associated with favorable clinical outcomes to immunotherapy in patients with melanoma, renal cell carcinoma and sarcoma [28-30]. It has also been reported that the presence of B cells and TLSs could be a potential prognostic marker for NSCLC [31, 32]. Thus, these findings indicate that the involvement of B cells and TLSs in the TME seems to facilitate a better response to immunotherapy in LUAD.

Fig 1

Classical gene set enrichment analysis (GSEA) in LUAD.

A–F, Lollipop plots depicting the GSEA results in the following samples: pretreatment tumor tissues (Pre-tissue) of responders (A) and nonresponders (B), pretreatment WB (Pre-WB) of responders (C) and nonresponders (D), and on-treatment WB (On-WB) of responders (E) and nonresponders (F). The X-axes show the normalized enrichment score (NES), and the Y-axes show gene sets ranked among the top 30 enriched gene sets with adjusted p-value < 0.10 (in descending order of NES). The dot size is proportional to the size of the corresponding gene set. The dot color indicates the adjusted p-value.

Although gene sets enriched in nonresponders were extremely similar between pretreatment tumor tissues and pretreatment WB, gene sets enriched in responders were quite different between these sources (Fig 1A–1D and S5 and S6 Tables). Notably, in pretreatment WB, gene sets significantly enriched in responders were those related to type I/II interferon (IFN) signaling (e.g., ‘RESPONSE TO TYPE I INTERFERON’ [p = 2.383 × 10−14, NES = 3.049], ‘INTERFERON-GAMMA-MEDIATED SIGNALING PATHWAY’ [p = 6.655 × 10−14, NES = 2.946], ‘RESPONSE TO INTERFERON-GAMMA’ [p = 2.968 × 10−18, NES = 2.939], ‘RESPONSE TO INTERFERON-BETA’ [p = 2.708 × 10−7, NES = 2.720], and ‘CELLULAR RESPONSE TO INTERFERON-BETA’ [p = 3.193 × 10−5, NES = 2.506]) and the process of antigen processing and presentation (e.g., ‘ANTIGEN PROCESSING AND PRESENTATION OF EXOGENOUS PEPTIDE ANTIGEN VIA MHC CLASS I’ [p = 1.604 × 10−10, NES = 2.755], ‘ANTIGEN PROCESSING AND PRESENTATION’ [p = 5.540 × 10−17, NES = 2.741], ‘ANTIGEN PROCESSING AND PRESENTATION OF PEPTIDE ANTIGEN’ [p = 1.532 × 10−15, NES = 2.728], and ‘ANTIGEN PROCESSING AND PRESENTATION OF PEPTIDE ANTIGEN VIA MHC CLASS I’ [p = 4.537 × 10−11, NES = 2.706]) (Fig 1C and S6 Table). Gene sets related to host defense against viral infection (e.g., ‘DEFENSE RESPONSE TO VIRUS’ [p = 1.799 × 10−17, NES = 2.750], ‘RESPONSE TO VIRUS’ [p = 3.428 × 10−18, NES = 2.583], ‘NEGATIVE REGULATION OF VIRAL LIFE CYCLE’ [p = 1.072 × 10−7, NES = 2.578], ‘NEGATIVE REGULATION OF VIRAL GENOME REPLICATION’ [p = 4.863 × 10−7, NES = 2.546], ‘REGULATION OF VIRAL LIFE CYCLE’ [p = 9.7864 × 10−10, NES = 2.464], and ‘NEGATIVE REGULATION OF VIRAL PROCESS’ [p = 1.642 × 10−7, NES = 2.392]), where type I IFN signaling plays a key role, were also enriched in responders (Fig 1C, S7 Fig and S6 Table). Type I IFNs, such as IFN-α and IFN-β, indirectly elicit antitumor immune responses in the TME by stimulating the maturation of dendritic cells, increasing the expression of perforin and granzymes in cytotoxic T cells, promoting the survival of memory T cells, and inactivating the suppressive function of regulatory T (Treg) cells. In addition, type I IFNs enhance antitumor immunity directly through the inhibition of tumor cell proliferation and the acceleration of senescence and apoptosis [33, 34]. Type II IFN (IFN-γ) also supports antitumor immunity by augmenting the function of tumor-infiltrating immune cells, inactivating suppressive Treg cells, and modulating stromal cell function to alter metabolism and inhibit angiogenesis [35]. Both type I and type II IFN signaling can promote the process of antigen processing and presentation, which means that gene sets related to antigen processing and presentation (hereinafter referred to as ‘APP signatures’) are closely linked to those related to IFN signaling (hereinafter referred to as ‘IFN signatures’). Moreover, type II IFN can induce the expression of PD-L1 on tumor cells and tumor-associated macrophages (TAMs) in the TME, which is an ingenious strategy that tumor cells employ to evade the host immune system [35]. In contrast to pretreatment WB of responders, pretreatment tumor tissues of responders did not show strong enrichment of IFN and APP signatures (Fig 1A and S5 Table). We presume that this difference may arise from the presence of immunosuppressive conditions in the TME mainly due to PD-1/PD-L1 signaling; that is, tumor-infiltrating immune cells are inactivated in the immunosuppressive TME of responders, while WB of responders is entirely unaffected by the TME.

In on-treatment WB, APP signatures, especially gene sets related to the process of antigen processing and presentation via MHC class I molecules (e.g., ‘ANTIGEN PROCESSING AND PRESENTATION OF EXOGENOUS PEPTIDE ANTIGEN VIA MHC CLASS I’ [p = 3.337 × 10−7, NES = 2.565], and ‘ANTIGEN PROCESSING AND PRESENTATION OF PEPTIDE ANTIGEN VIA MHC CLASS I’ [p = 3.795 × 10−7, NES = 2.453]), were significantly enriched in responders (Fig 1E and S7 Table). This enrichment of APP signatures in responders seems to reflect a durable antitumor immune response elicited by nivolumab monotherapy. Interestingly, gene sets related to mitochondrial metabolism (e.g., ‘ATP SYNTHESIS COUPLED ELECTRON TRANSPORT’ [p = 4.044 × 10−6, NES = 2.442], ‘OXIDATIVE PHOSPHORYLATION’ [p = 1.482 × 10−6, NES = 2.302], ‘RESPIRATORY ELECTRON TRANSPORT CHAIN’ [p = 3.639 × 10−5, NES = 2.219], ‘ATP SYNTHESIS COUPLED PROTON TRANSPORT’ [p = 5.748 × 10−3, NES = 2.216], and ‘MITOCHONDRIAL ELECTRON TRANSPORT, NADH TO UBIQUINONE’ [p = 4.266 × 10−3, NES = 2.042]) were significantly enriched in responders (Fig 1E and S7 Table). Mitochondrial metabolism supports proinflammatory signaling; in addition, the proinflammatory milieu can reprogram mitochondrial metabolism. The electron transport chain produces adenosine triphosphate (ATP) and reactive oxygen species (ROS) via coupling with oxidative phosphorylation (OXPHOS), which can drive the differentiation and activation of T cells [36, 37]. These findings indicate that mitochondrial metabolism may be persistently enhanced to effectively support antitumor immunity in on-treatment WB of responders.

To verify the association between PFS and the gene signatures enriched in responders irrespective of phenotypic information (e.g., responders vs. nonresponders), we used an unsupervised single-sample gene set enrichment scoring approach. Among the so-called unsupervised, nonparametric methods, ssGSEA (single-sample gene set enrichment analysis) [38] and GSVA (gene set variation analysis) [39] have been described as the most common methods. Both ssGSEA and GSVA, however, need expression data and phenotypic information for all samples—the former to normalize enrichment scores and the latter to conduct kernel density estimation to approximate the cumulative density function. Thus, we employed an alternative unsupervised method, singscore, which is a truly single-sample gene set enrichment scoring method [22].

Using the singscore method, we computed and evaluated single-sample gene set enrichment scores of the enriched gene sets identified through the supervised GSEA method above (hereinafter referred to as ‘GSEA gene sets’; Fig 2A, S8–S10 Figs and S8–S10 Tables). The single-sample gene set enrichment scores of GSEA gene sets in pretreatment tumor tissues were not significantly different between responders and nonresponders (S8 Fig). In pretreatment WB of responders, the single-sample enrichment scores showed that IFN and APP signatures were strongly enriched with high reproducibility, indicating that they were robustly enriched regardless of whether the enrichment analysis was supervised or unsupervised (Fig 2A). The single-sample enrichment scores of APP signatures were also significantly higher in on-treatment WB of responders (S10 Fig). In nonresponders, the enrichment scores of type I IFN signaling in on-treatment WB were significantly increased compared with pretreatment WB (‘IFN_I’, p = 0.0128; ‘IFNB1’, p = 0.0108; ‘IFNB2’, p = 0.0025) (Fig 2B). In responders, by contrast, the enrichment scores of type I IFN signaling in on-treatment WB showed no significant differences compared with pretreatment WB (Fig 2B). It has been reported that sustained type I IFN signaling contributes to resistance to anti-PD-1 monotherapy [34, 40]. Hence, these findings suggest that nivolumab-induced activation of type I IFN signaling may be a predictive biomarker for worse clinical outcomes in LUAD patients treated with nivolumab monotherapy.

Fig 2

Single-sample enrichment analysis (singscore) in LUAD.

A, For representative GSEA gene sets that were enriched in pretreatment WB of responders, the single-sample gene set enrichment scores and dispersions were calculated using singscore and visualized on scatter plots. Red circles denote responders (n = 5); cyan circles, nonresponders (n = 15). The enrichment scores were analyzed using the Wilcoxon rank sum test. All calculated p-values are shown on the plots (*p < 0.05, **p < 0.01, and ***p < 0.001). B, Box plots indicating the single-sample enrichment scores of IFN and APP signatures calculated from four different groups with LUAD: pretreatment WB of nonresponders (Pre-Tx/NR, n = 15), on-treatment WB of nonresponders (On-Tx/NR, n = 11), pretreatment WB of responders (Pre-Tx/R, n = 5) and on-treatment WB of responders (On-Tx/R, n = 4). Red dots denote responders; cyan dots, nonresponders. The differences in the single-sample enrichment scores between the groups (i.e., ‘Pre-Tx/NR vs. On-Tx/NR’, ‘Pre-Tx/NR vs. On-Tx/R’, ‘Pre-Tx/R vs. On-Tx/R’ and ‘On-Tx/NR vs. On-Tx/R’) were evaluated by the Wilcoxon rank sum test. Only significant p-values (p < 0.05) are shown on the plots (*p < 0.05, **p < 0.01, and ***p < 0.001).

Moreover, we noted that the enrichment scores of APP signatures exhibited no significant differences between pretreatment and on-treatment WB. As neoantigens are generated from mutations, the higher the TMB, the greater the chance that some of the neoantigens presented by MHC proteins will be immunogenic and hence enable the induction of anti-tumor immune response. In fact, accumulating evidence has indicated that high TMB is correlated with better clinical outcomes in NSCLC patients with anti-PD-1/PD-L1 therapy [41]. Additionally, several studies have reported that tumor mutation burden (TMB) in blood (or circulating tumor DNA) correlates with TMB in tumor tissue, and that high TMB in blood may serve as a potential biomarker of clinical benefit in NSCLC patients with anti-PD-1/PD-L1 therapy [4, 42]. Given that TMB in blood can be a surrogate for TMB in tumor tissue, it is rational to suggest that TMB in blood reflect the immunogenicity of tumor cells. Therefore, the comparable levels of APP signatures before and after nivolumab monotherapy indicate that nivolumab monotherapy has no significant impact on the immunogenicity of the tumor itself. Based on these findings, we presume that the tumors of responders are inherently sufficiently immunogenic to effectively elicit antigen processing and presentation for antitumor immunity. Clearly, it seems that elevated IFN signaling reveals the magnitude of antitumor immunity, which in turn upregulates PD-L1 expression on tumor cells to induce PD-1/PD-L1 signaling-mediated immune evasion. Thus, the PD-1/PD-L1 immune evasion axis can be one of the primary targets of nivolumab monotherapy in LUAD patients.

For each single-sample enrichment score, we next calculated the Spearman correlation coefficient (ρ) to estimate the correlation between the PFS time and the enrichment score (Fig 3A and S11A Fig). Consequently, we observed that both IFN and APP signatures were significantly correlated with PFS time (e.g., IFN_I, ρ = 0.590, p = 0.0061; IFNB1, ρ = 0.546, p = 0.0127; IFNG1, ρ = 0.526, p = 0.0173; APP1, ρ = 0.511, p = 0.0214). Gene sets related to host defense against viral infection (hereinafter, referred to as ‘VIRUS signatures’), which have extremely strong correlations with IFN signatures, also exhibited significant correlations with PFS time (e.g., VIRUS1, ρ = 0.651, p = 0.0019; VIRUS2, ρ = 0.624, p = 0.0033). Nonparametric regression models using cubic regression splines indicated that the enrichment scores of a single gene set in the IFN and APP signatures in pretreatment WB could be candidate biomarkers to predict PFS time (e.g., IFNB1, R-squared [R-sq] = 0.5861, AICc = 39.0470, p = 0.0007; APP1, R-sq = 0.4353, AICc = 41.9426, p = 0.0010) (Fig 3B and S11B Fig). When the enrichment scores for two gene sets in the IFN and APP signatures were combined, LUAD patients with high values for both enrichment scores tended to have longer PFS time (Fig 4A and S12 Fig). By stratifying the LUAD patients into subsets with high and low enrichment scores (stratification by the median), we found that there were quite large differences in PFS time between the patients with high enrichment of signatures and those with low enrichment of signatures (Fig 4B and S13 Fig). Multivariate Cox regression analysis identified that many combinations of two of the IFN and APP signatures remained to be independent predictive factors of PFS time (Fig 4C and S14 Fig). Collectively, these findings suggest that IFN and APP signatures in pretreatment WB may have potential predictive performance in LUAD patients treated with nivolumab monotherapy.

Fig 3

IFN and APP signatures predict the nivolumab response in LUAD.

A, Spearman correlation matrix between PFS and single-sample gene set enrichment scores of IFN and APP signatures in pretreatment WB. The upper triangular region shows the values of Spearman’s ρ correlation coefficients (significant correlations are in bold; *p < 0.05, **p < 0.01, and ***p < 0.001). In the lower triangular region, positive correlations are visualized in red and negative correlations in blue. The color intensity is proportional to the value of Spearman’s ρ, and the size of the circle to the p-value. B, Scatter plots showing the relationships between PFS and the enrichment score for a single gene set in the IFN and APP signatures in pretreatment WB, with a fitted line representing the regression model using a cubic spline and 95% confidence interval. The accuracy of the fit was assessed by calculating the adjusted R-squared (R-sq) and p-values (*p < 0.05, **p < 0.01, and ***p < 0.001).

Fig 4

Two gene sets in the IFN and APP signatures as candidate biomarkers for the nivolumab response in LUAD.

A, Bubble plots showing the relationships between PFS and the enrichment scores for two gene sets in the IFN and APP signatures in pretreatment WB. Each bubble represents a patient, and the size of each bubble is proportional to the PFS time. On a gradient color scale based on the PFS time, bubbles representing responders were assigned colors ranging from white to dark red; nonresponders, ranging from white to lavender. B, Kaplan-Meier PFS curves for patients stratified by the enrichment score (‘Low’ vs. ‘High’) for two gene sets in the IFN and APP signatures. Patients with both scores above the median were defined as ‘High’ and the others as ‘Low’. The p-values were calculated by the two-sided log-rank test (*p < 0.05, **p < 0.01, and ***p < 0.001). C, Forest plot showing multivariate Cox regression analysis for potential factors associated with prolonged or shortened PFS time. Squares represent estimated hazard ratios and whiskers represent the 95% confidence intervals. Hazard ratios less than 1 indicate improved PFS time (*p < 0.05).

A cohort of patients with LUSC

In the LUSC cohort (n = 18), differential expression analysis between responders and nonresponders was performed using the RNA-seq datasets of pretreatment tumor tissues (n = 13), pretreatment WB (n = 17) and on-treatment WB (n = 15). A total of 424 DEGs were identified in pretreatment tumor tissues. Of these 424 DEGs, 36 genes were upregulated and 388 genes were downregulated in responders (n = 3) compared to nonresponders (n = 10) (S15A and S15B Fig and S11 Table). In marked contrast to LUAD, LUSC had many more DEGs in tumor tissues than in WB. In fact, only three DEGs were identified in on-treatment WB, among which two genes were upregulated and one gene was down-regulated in responders (n = 2) compared to nonresponders (n = 13) (S15E and S15F Fig and S12 Table). We found no DEGs in pretreatment WB (S15C and S15D Fig). A heatmap of DEGs obtained from pretreatment tumor tissues shows the relations between the expression levels of DEGs and clinical factors such as PFS, PD-L1 TPS and Brinkman index (S16 Fig).

Surprisingly, the GSEA results indicated that the gene sets enriched in responders or nonresponders were fairly similar between pretreatment and on-treatment WB in LUSC (Fig 5C–5F, S17 Fig and S14 and S15 Tables). In WB of patients with LUSC, we identified no or few DEGs between responders and nonresponders and almost no changes in the enrichment pattern across the gene sets between the pre- and post-nivolumab monotherapy settings. Nivolumab monotherapy seemed to have only limited impacts on WB, indicating that a transcriptome analysis of WB may be unable to elucidate systemic effects and predict the clinical response to nivolumab monotherapy in LUSC patients. We therefore focused on pretreatment tumor tissues of LUSC for the following analysis (Fig 5A, 5B, S18 Fig and S14 Table).

Fig 5

Classical GSEA in LUSC.

A–F, Lollipop plots depicting the GSEA results in the following samples: pretreatment tumor tissues (Pre-tissue) of responders (A) and nonresponders (B), pretreatment WB (Pre-WB) of responders (C) and nonresponders (D), and on-treatment WB (On-WB) of responders (E) and nonresponders (F). The X-axes show the normalized enrichment score (NES); the Y-axes, gene sets ranked among the top 30 enriched gene sets with adjusted p-value < 0.10 (in descending order of NES). The dot size is proportional to the size of the corresponding gene set. The dot color indicates the adjusted p-value.

In pretreatment tumor tissues, gene sets related to mitochondrial metabolism (e.g., ‘MITOCHONDRIAL RESPIRATORY CHAIN COMPLEX ASSEMBLY’ [p = 4.225 × 10−14, NES = 2.919], ‘ATP SYNTHESIS COUPLED ELECTRON TRANSPORT’ [p = 1.102 × 10−12, NES = 2.852], ‘NADH DEHYDROGENASE COMPLEX ASSEMBLY’ [p = 2.407 × 10−10, NES = 2.850], ‘OXIDATIVE PHOSPHORYLATION’ [p = 2.145 × 10−15, NES = 2.838], ‘MITOCHONDRIAL GENE EXPRESSION’ [p = 5.865 × 10−17, NES = 2.803], ‘INNER MITOCHONDRIAL MEMBRANE ORGANIZATION’ [p = 1.819 × 10−9, NES = 2.792], ‘MITOCHONDRIAL TRANSLATION’ [p = 4.873 × 10−14, NES = 2.750], ‘CRISTAE FORMATION’ [p = 6.394 × 10−8, NES = 2.711], ‘RESPIRATORY ELECTRON TRANSPORT CHAIN’ [p = 5.114 × 10−11, NES = 2.694], and ‘MITOCHONDRIAL ELECTRON TRANSPORT, NADH TO UBIQUINONE’ [p = 2.189 × 10−7, NES = 2.596]) were significantly enriched in responders (Fig 5A, S18A Fig and S13 Table). Additionally, gene sets related to organization of the tumor microenvironment (e.g., ‘EXTRACELLULAR STRUCTURE ORGANIZATION’ [p = 9.555 × 10−31, NES = −2.611], ‘COLLAGEN FIBRIL ORGANIZATION’ [p = 2.039 × 10−8, NES = −2.508], ‘COLLAGEN CATABOLIC PROCESS’ [p = 4.344 × 10−6, NES = −2.291], ‘REGULATION OF CELLULAR EXTRAVASATION’ [p = 1.381 × 10−5, NES = −2.261], ‘COLLAGEN METABOLIC PROCESS’ [p = 1.235 × 10−7, NES = −2.229], ‘CELLULAR EXTRAVASATION’ [p = 2.989 × 10−6, NES = −2.226], ‘MUCOPOLYSACCHARIDE METABOLIC PROCESS’ [p = 9.872 × 10−8, NES = −2.221], ‘POSITIVE REGULATION OF CELLULAR EXTRAVASATION’ [p = 3.971 × 10−5, NES = −2.219], ‘EXTRACELLULAR MATRIX DISASSEMBLY’ [p = 1.433 × 10−6, NES = −2.199], ‘NEGATIVE REGULATION OF HOMOTYPIC CELL-CELL ADHESION’ [p = 1.000 × 10−4, NES = −2.199], and ‘TISSUE REMODELING’ [p = 8.091 × 10−8, NES = −2.114]) were significantly enriched in nonresponders (Fig 5B, S18B Fig and S13 Table) [43-45].

Among these gene sets, singscore reproducibly identified three gene sets related to mitochondrial functions, namely, ‘MITOCHONDRIAL GENE EXPRESSION’, ‘INNER MITOCHONDRIAL MEMBRANE ORGANIZATION’ and ‘CRISTAE FORMATION’ (hereinafter referred to as ‘mitochondrial signatures’), as significantly enriched in responders (S19 Fig and S16 Table). In addition, seven gene sets related to the regulation of the TME, namely, ‘COLLAGEN CATABOLIC PROCESS’, ‘COLLAGEN METABOLIC PROCESS’, ‘MUCOPOLYSACCHARIDE METABOLIC PROCESS’, ‘POSITIVE REGULATION OF CELLULAR EXTRAVASATION’, ‘EXTRACELLULAR MATRIX DISASSEMBLY’, ‘NEGATIVE REGULATION OF HOMOTYPIC CELL-CELL ADHESION’ and ‘TISSUE REMODELING’ (hereinafter referred to as ‘TME signatures’), were significantly enriched in nonresponders (Fig 6A and S16 Table). It has been well documented that the TME signatures identified contribute to the structural organization and metabolic regulation of the extracellular matrix (ECM), which primarily consists of collagens, mucopolysaccharides and proteoglycans [43-45]. ECM remodeling in the TME plays important roles in various biological processes, including proliferation, adhesion, angiogenesis and metastasis, to affect tumor progression [46]. In many solid tumors, the TME exhibits excessive deposition of collagen and other ECM components. Dense collagen and other ECM components give rise to an immunosuppressive and hypoxic microenvironment. Furthermore, the hypoxic TME can not only accelerate tumor proliferation and metastasis but also promote the development of immunosuppressive conditions favoring tumor immune evasion [47]. Thus, it is conceivable that the enrichment of TME signatures in nonresponders may reflect the more immunosuppressive TME. Conversely, the enrichment of mitochondrial signatures in responders probably means a relatively oxygen-rich and less immunosuppressive TME.

Fig 6

Singscore in LUSC.

A, For representative GSEA gene sets that were enriched in the pretreatment tumor tissues of responders or nonresponders, the single-sample gene set enrichment scores and dispersions were calculated using singscore and visualized on scatter plots. Red circles denote responders (n = 3); cyan circles, nonresponders (n = 10). The enrichment scores were analyzed using the Wilcoxon rank sum test. All calculated p-values are shown on the plots (*p < 0.05). B, Box plots indicating the single-sample enrichment scores of IFN and APP signatures calculated from four different groups with LUSC: pretreatment WB of nonresponders (Pre-Tx/NR, n = 13), on-treatment WB of nonresponders (On-Tx/NR, n = 13), pretreatment WB of responders (Pre-Tx/R, n = 4) and on-treatment WB of responders (On-Tx/R, n = 2). Red dots denote responders; cyan dots, nonresponders. The differences in the single-sample enrichment scores between the groups (i.e., ‘Pre-Tx/NR vs. On-Tx/NR’, ‘Pre-Tx/NR vs. On-Tx/R’, ‘Pre-Tx/R vs. On-Tx/R’ and ‘On-Tx/NR vs. On-Tx/R’) were evaluated by the Wilcoxon rank sum test. No significant differences were observed between the groups.

We next investigated whether the enrichment scores of the IFN and APP signatures, which were candidate predictive biomarkers in LUAD patients treated with nivolumab monotherapy, also dynamically changed between responders and nonresponders or between the pre- and post-nivolumab monotherapy settings in WB of LUSC patients (Fig 6B). We observed no significant differences in the enrichment scores of the IFN and APP signatures between responders and nonresponders in pretreatment WB, indicating that both responders and nonresponders have similar levels of preexisting antitumor immunity. The increase in type I IFN signaling after nivolumab monotherapy observed in WB of nonresponders with LUAD was not observed in those with LUSC (Fig 6B). Collectively, these findings suggest that nivolumab success in LUSC depends entirely on the extent of the immunosuppressive TME, not on the inherent immunogenicity of the tumor itself.

Spearman rank correlation analysis demonstrated significant correlations between PFS time and the enrichment scores of the mitochondrial and TME signatures (Fig 7A). Notably, the TME signatures were negatively correlated with the mitochondrial signatures, supporting the hypothesis that the oxygenation and immunomodulation status of the TME can explain the enrichment status of both the mitochondrial and TME signatures. Among the TME signature scores, the ADHES score was most strongly correlated with PFS time (ρ = −0.786, p = 0.0015). Using a cubic regression spline, we found that the ADHES score is a potential predictive biomarker for PFS in LUSC patients treated with nivolumab monotherapy (R-sq = 0.6009, AICc = 28.1764, p = 0.0009) (Fig 7B). Gene sets related to the regulation of platelet activation (PLAT1 and PLAT2; hereinafter referred to as ‘platelet signatures’) were negatively correlated with PFS time (PLAT1, ρ = −0.687, p = 0.0095; PLAT2, ρ = −0.659, p = 0.0142) and positively correlated with TME signatures (Fig 7A). Platelets can suppress T cell antitumor responses through the production and activation of immunosuppressive factors [48, 49]. Hence, platelet signatures are inferred to reflect the immunosuppressive conditions in the TME, similar to TME signatures. The fitted models of the enrichment scores of platelet signatures exhibited a relatively low but statistically significant predictive power for PFS time (PLAT1, R-sq = 0.3135, AICc = 35.2646, p = 0.0289; PLAT2, R-sq = 0.4101, AICc = 34.4125, p = 0.0116) (Fig 7B and S22 Fig). When two gene sets in the TME and platelet signatures were combined, the LUSC patients with low values for both scores tended to have longer PFS times (Fig 8A and S23 Fig). When the LUSC patients were stratified into subsets with high and low scores (stratification by the median), there were quite large differences in PFS time between the patients with high enrichment of signatures and those with low enrichment of signatures (Fig 8B and S24 Fig). In multivariate Cox regression analysis revealed that the COL1/TISSUE signature remained an independent predictive factor of PFS time (Fig 8). Other combinations of two of the TME signatures seemed to be slightly better predictive factors compared to clinical factors (age, ECOG PS, Brinkman index, PD-L1 TPS and CNS metastasis), but without statistical significance (S25 Fig). These findings indicate that TME signatures may have potential predictive performance for PFS in LUSC patients treated with nivolumab monotherapy.

Fig 7

TME signatures predict the nivolumab response in LUSC.

A, Spearman correlation matrix between PFS and single-sample gene set enrichment scores of MIT and TME signatures in pretreatment tumor tissues. The upper triangular region shows the values of Spearman’s ρ correlation coefficients (significant correlations are in bold; *p < 0.05, **p < 0.01 and ***p < 0.001). In the lower triangular region, positive correlations are visualized in red and negative correlations in blue. The color intensity is proportional to the value of Spearman’s ρ and the size of the circle to the p-value. B, Scatter plots showing the relationships between PFS and the enrichment score for single gene sets in the TME signature in pretreatment tumor tissues, with a fitted line representing the regression model using a cubic spline and 95% confidence interval. The accuracy of the fit was assessed by calculating the adjusted R-squared (R-sq) and p-values (*p < 0.05, **p < 0.01 and ***p < 0.001).

Fig 8

Two gene sets in the TME signatures as candidate biomarkers for the nivolumab response in LUSC.

A, Bubble plots showing the relationships between PFS and the enrichment scores for two gene sets in the TME signature in pretreatment tumor tissues. Each bubble represents a patient, and the size of the bubble is proportional to the PFS time. On a gradient color scale based on the PFS time, bubbles representing responders were assigned colors ranging from white to dark red; nonresponders, ranging from white to lavender. B, Kaplan-Meier PFS curves for patients stratified by the enrichment score (‘Low’ vs. ‘High’) of two gene sets in the TME signature. Patients with both scores below the median are defined as ‘Low’; the others, as ‘High’. The p-values were calculated by the two-sided log-rank test (*p < 0.05). C, Forest plot showing multivariate Cox regression analysis for potential factors associated with prolonged or shortened PFS time. Squares represent estimated hazard ratios and whiskers represent the 95% confidence intervals. Hazard ratios less than 1 indicate improved PFS time (*p < 0.05).

Discussion

PD-L1 expression on tumor cells has been widely accepted as a predictive biomarker for therapeutic decision making in NSCLC, although its accuracy is limited and it has virtually no predictive value in patients with LUSC. To the best of our knowledge, the cause of the discrepancy in the reliability of its predictive power between LUAD and LUSC remains to be clarified. In this study, we employed whole-transcriptome sequencing and a single-sample enrichment analysis method—singscore—and found that the discrepancy can be explained by the difference in the immunogenicity of the tumor itself and the immunosuppressive conditions in the TME.

In LUAD, we observed that the IFN and APP signatures, which are closely related to each other and functionally cooperate to activate the antitumor immune response, were significantly enriched in pretreatment WB of responders (Fig 1 and S6 Table). This enrichment of the IFN and APP signatures suggests that responders may have preexisting antitumor immunity prior to nivolumab monotherapy. In contrast, transcriptomic data of pretreatment tumor tissues showed no significant enrichment of either IFN or APP signatures, suggesting that antitumor immunity is locally disturbed in the TME to support tumor progression (Fig 1 and S5 Table). These findings highlighted the common features in responders as follows: First, the tumors possess sufficiently high immunogenicity to induce effective antigen presentation to host immune components (S26 Fig). Second, the host immune system activates IFN signaling to activate tumor immunosurveillance, which is disrupted in the TME so that tumor immune evasion (i.e., adaptive immune resistance) is established locally. Given that activated type II IFN signaling can upregulate PD-L1 expression on tumor cells [35], it is highly likely that the PD-1/PD-L1 axis is responsible for the establishment of adaptive immune resistance. Finally, nivolumab monotherapy restores antitumor activity by inhibiting the PD-1/PD-L1 axis. In this regard, it is possible to say that nivolumab functions just as a molecular targeted agent in LUAD patients and that PD-L1 expression on tumor cells helps to predict the efficacy of nivolumab monotherapy.

In LUSC, TME signatures were significantly enriched in nonresponders, and this enrichment indicates the immunological condition within the TME. We found that the enrichment scores of TME signatures were negatively correlated with PFS time (Figs 7 and 8), indicating that patients with tumors strongly protected by the immunosuppressive TME are unlikely to benefit from nivolumab monotherapy. More importantly, IFN and APP signatures were strongly enriched in pretreatment WB of responders with LUAD, but this relationship was not observed in those with LUSC (Fig 6B). This observation raises the possibility that there may be no difference in the level of preexisting anti-umor immunity or, alternatively, in the immunogenicity of tumors between responders and nonresponders with LUSC. In support of this hypothesis, no significant enrichment of IFN and APP signatures was observed in on-treatment WB of responders (S15 Table). In fact, we identified the similarity in the patterns of enriched gene sets between pretreatment and on-treatment WB of patients with LUSC. Thus, nivolumab monotherapy has just a local impact in LUSC patients, a striking contrast to its systemic impact in LUAD patients treated with nivolumab monotherapy. Specifically, the clinical response to nivolumab monotherapy is actually determined by the extent of the immunosuppressive TME, where immunosuppressive factors other than the PD-1/PD-L1 axis may be considered to be critical, because the preexisting IFN and APP signatures exhibited no differences between responders and nonresponders, as mentioned above. This hypothesis is in line with the well-known finding that PD-L1 expression does not correlate significantly with clinical outcomes in LUSC patients treated with anti-PD-1 monotherapy [1]. Collectively, these findings prompt us to note that nivolumab monotherapy functions just as an immunomodulating agent and cannot overcome the highly immunosuppressive TME alone ().

One limitation of our study is the lack of validation in an independent cohort. However, by using a true single-sample enrichment approach, singscore, we have devised a new workflow for identifying gene expression signatures to predict a patient’s response to immunotherapy and to gain a deeper understanding of cancer biology. For example, combination strategies to enhance the immunogenicity of the tumor itself (e.g., cancer vaccines and CAR-T therapy) [4, 41–52] can be expected to improve the clinical response to nivolumab monotherapy in patients with LUAD, whereas combinational strategies to overcome the immunosuppressive TME are needed in LUSC [53, 54]. We envision that future studies will provide a blueprint for innovating combination immunotherapy approaches and optimizing patient selection and treatment strategies to improve long-term clinical outcomes in NSCLC.

A schematic illustration of patient enrollment and sample collection.

(TIFF)

Click here for additional data file.

Visualization of tumor response.

A, Waterfall plot of the best percentage change from baseline during nivolumab monotherapy according to RECIST v1.1. B, Swimmer plot of all 40 patients treated with nivolumab monotherapy. PD, progressive disease.

(TIFF)

Click here for additional data file.

Kaplan-Meier estimates of progression-free survival (PFS) and overall survival (OS) of patients according to different stratification schemes.

A–B, the total patient cohort and C–D, LUAD or E–F, LUSC patients with PD-L1 TPS ≥ 1% versus < 1%. The p-values were calculated by the two-sided log-rank test.

(TIFF)

Click here for additional data file.

The association between the immune cytolytic activity in tumor tissues and clinical outcomes of patients in this cohort.

For each tumor sample, the immune cytolytic activity was estimated as the average expression level of the marker genes (CD8A, CD8B, GZMA, GZMB and PRF1). Patients with all expression levels above the average are defined as ‘High’; the others, as ‘Low’. In the upper panel, above-average values are shown in red. The lower panel illustrates Kaplan-Meier estimates of PFS and OS of patients stratified by the estimated immune cytolytic activity. The p-values were calculated by the two-sided log-rank test.

(TIFF)

Click here for additional data file.

Differential gene expression analysis in LUAD.

A–B, MA plot (A) and volcano plot (B) of DEGs in pretreatment tumor tissues. C–D, MA plot (C) and volcano plot (D) of DEGs in pretreatment WB. E–F, MA plot (E) and volcano plot (F) of DEGs in on-treatment WB. Red dots represent DEGs [adjusted p-value < 0.10 and |log2(fold change)| ≥ 1]. Triangles and diamonds represent genes with log2(fold change) and normalized counts, respectively, out of the plot scale. The horizontal lines in the MA plots and vertical lines in the volcano plots indicate the thresholds log2(fold change) = 1 or −1. The horizontal lines in the volcano plots indicate the threshold −log10(adjusted p-value) = 1.

(TIFF)

Click here for additional data file.

Heatmap representation of DEGs in LUAD.

Heatmaps of DEGs between responders and nonresponders with hierarchical clustering of samples: A, from pretreatment tumor tissues (n = 15), B, pretreatment WB (n = 20), and C, on-treatment WB (n = 15). The DEGs clearly differentiated between responders and nonresponders in all three datasets.

(TIFF)

Click here for additional data file.

Enrichment maps for representative gene sets significantly enriched in pretreatment WB of responders with LUAD.

Each node denotes a distinct gene set, and the size of the node is proportional to the number of genes in the set. The thickness of the edges (pale blue lines) represents the degree of overlap between the two connected gene sets.

(TIFF)

Click here for additional data file.

Scatter plots of single-sample gene set enrichment scores in pretreatment tumor tissues.

For each gene set in the top 20 GSEA gene sets that were enriched in pretreatment tumor tissues of responders (A) and nonresponders (B) with LUAD, the enrichment scores and dispersions were calculated using singscore. Red circles denote responders (n = 4); cyan circles, nonresponders (n = 11). The enrichment scores were analyzed using the Wilcoxon rank sum test.

(TIFF)

Click here for additional data file.

Scatter plots of single-sample gene set enrichment scores in pretreatment WB.

For each gene set in the top 20 GSEA gene sets that were enriched in pretreatment WB of responders (A, except for those shown in Fig 2A) and nonresponders (B) with LUAD, the enrichment scores and dispersions were calculated using singscore. Red circles denote responders (n = 5); cyan circles, nonresponders (n = 15). The enrichment scores were analyzed using the Wilcoxon rank sum test (*p < 0.05 and **p < 0.01).

(TIFF)

Click here for additional data file.

Scatter plots of single-sample gene set enrichment scores in on-treatment WB.

For each gene set in the top 20 GSEA gene sets that were enriched in on-treatment WB of responders (A) and nonresponders (B) with LUAD, the enrichment scores and dispersions were calculated using singscore. Red circles denote responders (n = 4); cyan circles, nonresponders (n = 11). The enrichment scores were analyzed using the Wilcoxon rank sum test (*p < 0.05 and **p < 0.01).

(TIFF)

Click here for additional data file.

Relationships between PFS and single-sample enrichment scores for gene sets significantly enriched in responders with LUAD.

A, Spearman correlation matrix between PFS and the enrichment scores in pretreatment WB. The upper triangular region shows the values of the Spearman’s ρ correlation coefficients (significant correlations are in bold; *p < 0.05, **p < 0.01 and ***p < 0.001). In the lower triangular region, positive correlations are visualized in red and negative correlations in blue. The color intensity is proportional to the value of Spearman’s ρ and the size of the circle to the p-value. B, Scatter plots showing the relationships between PFS and the enrichment score of the above gene sets in pretreatment WB, with a fitted line representing the regression model using a cubic spline and 95% confidence interval. The accuracy of the fit was assessed by calculating the adjusted R-squared (R-sq) and p-values (*p < 0.05 and **p < 0.01).

(TIFF)

Click here for additional data file.

Bubble plots showing the relationships between PFS and two single-sample enrichment scores in pretreatment WB of LUAD patients.

Each bubble represents a patient, and the size of the bubble is proportional to the PFS time. On a gradient color scale based on the PFS time, bubbles representing responders were assigned colors ranging from white to dark red; nonresponders, ranging from white to lavender.

(TIFF)

Click here for additional data file.

Kaplan-Meier PFS curves for patients stratified by the single-sample enrichment score (‘Low’ vs. ‘High’) for two gene sets in the IFN and APP signatures.

LUAD patients with both scores above the median are defined as ‘High’; the others, as ‘Low’. The p-values were calculated by the two-sided log-rank test (*p < 0.05, **p < 0.01, and ***p < 0.001).

(TIFF)

Click here for additional data file.

Forest plot showing multivariate Cox regression analysis for potential factors associated with prolonged or shortened PFS time in LUAD patients.

Squares represent estimated hazard ratios and whiskers represent the 95% confidence intervals. Hazard ratios less than 1 indicate improved PFS time (*p < 0.05).

(TIFF)

Click here for additional data file.

Differential gene expression analysis in LUSC.

A–B, MA plot (A) and volcano plot (B) of DEGs in pretreatment tumor tissues. C–D, MA plot (C) and volcano plot (D) of DEGs in pretreatment WB. E–F, MA plot (E) and volcano plot (F) of DEGs in on-treatment WB. Red dots represent DEGs [adjusted p-value < 0.10 and |log2(fold change)| ≥ 1]. Triangles and diamonds represent genes with log2(fold change) and normalized counts, respectively, out of the plot scale. The horizontal lines in the MA plots and vertical lines in the volcano plots indicate the thresholds log2(fold change) = 1 or −1. The horizontal lines in the volcano plots indicate the threshold −log10(adjusted p-value) = 1.

(TIFF)

Click here for additional data file.

Heatmap representation of DEGs in LUSC.

Heatmaps of DEGs between responders and nonresponders with hierarchical clustering of samples from pretreatment tumor tissues (n = 13).

(TIFF)

Click here for additional data file.

The lists of the top 25 gene sets significantly enriched in WB of responders and nonresponders with LUSC.

Top 25 enriched gene sets in responders (A) and in nonresponders (B). The left panel shows the list for pretreatment WB; the right panel, on-treatment WB. The enriched gene sets common between pretreatment and on-treatment WB are shown in bold and connected by solid lines.

(TIFF)

Click here for additional data file.

Enrichment maps for representative gene sets significantly enriched in pretreatment tumor tissues of (A) responders and (B) nonresponders with LUSC. Each node denotes a distinct gene set, and the size of the node is proportional to the number of genes in the set. The thickness of the edges (pale blue lines) represents the degree of overlap between the two connected gene sets.

(TIFF)

Click here for additional data file.

For each gene set in the top 20 GSEA gene sets that were enriched in pretreatment tumor tissues of responders (A) and nonresponders (B) (except for those shown in Fig 5A) with LUSC, the enrichment scores and dispersions were calculated using singscore. Red circles denote responders (n = 3); cyan circles, nonresponders (n = 10). The enrichment scores were analyzed using the Wilcoxon rank sum test (*p < 0.05).

(TIFF)

Click here for additional data file.

For each gene set in the top 20 GSEA gene sets that were enriched in pretreatment WB of responders (A) and nonresponders (B) with LUSC, the enrichment scores and dispersions were calculated using singscore. Red circles denote responders (n = 4); cyan circles, nonresponders (n = 13). The enrichment scores were analyzed using the Wilcoxon rank sum test (*p < 0.05).

(TIFF)

Click here for additional data file.

For each gene set in the top 20 GSEA gene sets that were enriched in on-treatment WB of responders (A) and nonresponders (B) with LUSC, the enrichment scores and dispersions were calculated using singscore. Red circles denote responders (n = 2); cyan circles, nonresponders (n = 13). The enrichment scores were analyzed using the Wilcoxon rank sum test.

(TIFF)

Click here for additional data file.

Scatter plots showing the relationships between PFS and the single-sample enrichment score for gene sets significantly enriched in pretreatment tumor tissues of responders with LUSC.

The fitted line on each scatter plot represents the regression model using a cubic spline and 95% confidence interval. The accuracy of the fit was assessed by calculating the adjusted R-squared (R-sq) and p-values (*p < 0.05 and **p < 0.01).

(TIFF)

Click here for additional data file.

Bubble plots showing the relationships between PFS and two single-sample enrichment scores in the pretreatment tumor tissues of LUSC patients.

Each bubble represents a patient, and the size of the bubble is proportional to the PFS time. On a gradient color scale based on the PFS time, bubbles representing responders were assigned colors ranging from white to dark red; nonresponders, ranging from white to lavender.

(TIFF)

Click here for additional data file.

Kaplan-Meier PFS curves for patients stratified by the single-sample enrichment score (‘Low’ vs. ‘High’) for two gene sets in the selected TME signatures.

LUSC patients with both scores below the median are defined as ‘Low’; the others, as ‘High’. The p-values were calculated by the two-sided log-rank test (*p < 0.05 and **p < 0.01).

(TIFF)

Click here for additional data file.

Forest plot showing multivariate Cox regression analysis for potential factors associated with prolonged or shortened PFS time in LUSC patients.

Squares represent estimated hazard ratios and whiskers represent the 95% confidence intervals. Hazard ratios less than 1 indicate improved PFS time.

(TIFF)

Click here for additional data file.

Schematic depicting the different mechanisms of action of nivolumab monotherapy between LUAD and LUSC.

The success of nivolumab monotherapy depends on the inherent immunogenicity of the tumor itself in LUAD and the preexisting TME favoring an antitumor immune response in LUSC.

(TIFF)

Click here for additional data file.

Assessment of tumor response to nivolumab monotherapy according to RECIST v1.1.

(DOCX)

Click here for additional data file.

The lists of DEGs identified from pretreatment tumor tissue of LUAD patients.

(XLSX)

Click here for additional data file.

The lists of DEGs identified from pretreatment whole blood of LUAD patients.

(XLSX)

Click here for additional data file.

The lists of DEGs identified from on-treatment whole blood of LUAD patients.

(XLSX)

Click here for additional data file.

The results from classical GSEA (FGSEA) of pretreatment tumor tissue of LUAD patients.

(XLSX)

Click here for additional data file.

The results from classical GSEA (FGSEA) of pretreatment whole blood of LUAD patients.

(XLSX)

Click here for additional data file.

The results from classical GSEA (FGSEA) of on-treatment whole blood of LUAD patients.

(XLSX)

Click here for additional data file.

The results from single-sample enrichment analysis (singscore) of pretreatment tumor tissue of LUAD patients.

(XLSX)

Click here for additional data file.

The results from single-sample enrichment analysis (singscore) of pretreatment whole blood of LUAD patients.

(XLSX)

Click here for additional data file.

The results from single-sample enrichment analysis (singscore) of on-treatment whole blood of LUAD patients.

(XLSX)

Click here for additional data file.

The lists of DEGs identified from pretreatment tumor tissue of LUSC patients.

(XLSX)

Click here for additional data file.

The lists of DEGs identified from on-treatment whole blood of LUSC patients.

(XLSX)

Click here for additional data file.

The results from classical GSEA (FGSEA) of pretreatment tumor tissue of LUSC patients.

(XLSX)

Click here for additional data file.

The results from classical GSEA (FGSEA) of pretreatment whole blood of LUSC patients.

(XLSX)

Click here for additional data file.

The results from classical GSEA (FGSEA) of on-treatment whole blood of LUSC patients.

(XLSX)

Click here for additional data file.

The results from single-sample enrichment analysis (singscore) of pretreatment tumor tissue of LUSC patients.

(XLSX)

Click here for additional data file.

The results from single-sample enrichment analysis (singscore) of pretreatment whole blood of LUSC patients.

(XLSX)

Click here for additional data file.

The results from single-sample enrichment analysis (singscore) of on-treatment whole blood of LUSC patients.

(XLSX)

Click here for additional data file.

The raw count data from RNS-seq.

(CSV)

Click here for additional data file.

Transfer Alert

This paper was transferred from another journal. As a result, its full editorial history (including decision letters, peer reviews and author responses) may not be present.

9 Aug 2021

PONE-D-21-19886

Gene expression signatures as candidate biomarkers of response to PD-1 blockade in non-small cell lung cancers.

PLOS ONE

Dear Dr. aiba,

Thank you for submitting your manuscript to PLOS ONE. After careful consideration, we feel that it has merit but does not fully meet PLOS ONE’s publication criteria as it currently stands. Therefore, we invite you to submit a revised version of the manuscript that addresses the points raised during the review process.

We apologize our delayed response. Reviewers raised many constructive comments about your study. Due to the limitation of the small number of patients, you need to strengthen the statistical approach to provide clinical significance, instead of extensive description of the analysis. We hope you address all the reviewers' concerns.

Please submit your revised manuscript by Sep 23 2021 11:59PM. If you will need more time than this to complete your revisions, please reply to this message or contact the journal office at plosone@plos.org. When you're ready to submit your revision, log on to https://www.editorialmanager.com/pone/ and select the 'Submissions Needing Revision' folder to locate your manuscript file.

Please include the following items when submitting your revised manuscript:

A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). You should upload this letter as a separate file labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. You should upload this as a separate file labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. You should upload this as a separate file labeled 'Manuscript'.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter. Guidelines for resubmitting your figure files are available below the reviewer comments at the end of this letter.

If applicable, we recommend that you deposit your laboratory protocols in protocols.io to enhance the reproducibility of your results. Protocols.io assigns your protocol its own identifier (DOI) so that it can be cited independently in the future. For instructions see: http://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols. Additionally, PLOS ONE offers an option for publishing peer-reviewed Lab Protocol articles, which describe protocols hosted on protocols.io. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols.

We look forward to receiving your revised manuscript.

Kind regards,

Hyun-Sung Lee, M.D., Ph.D.

Academic Editor

PLOS ONE

Journal requirements:

When submitting your revision, we need you to address these additional requirements.

1. Please ensure that your manuscript meets PLOS ONE's style requirements, including those for file naming. The PLOS ONE style templates can be found at

https://journals.plos.org/plosone/s/file?id=wjVg/PLOSOne_formatting_sample_main_body.pdf and

https://journals.plos.org/plosone/s/file?id=ba62/PLOSOne_formatting_sample_title_authors_affiliations.pdf

2. In order to ensure your study is reproducible, please provide the databases accession numbers included in your study.

3. We note that the grant information you provided in the ‘Funding Information’ and ‘Financial Disclosure’ sections do not match.

When you resubmit, please ensure that you provide the correct grant numbers for the awards you received for your study in the ‘Funding Information’ section.

4. We note that you have stated that you will provide repository information for your data at acceptance. Should your manuscript be accepted for publication, we will hold it until you provide the relevant accession numbers or DOIs necessary to access your data. If you wish to make changes to your Data Availability statement, please describe these changes in your cover letter and we will update your Data Availability statement to reflect the information you provide.

[Note: HTML markup is below. Please do not edit.]

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Partly

Reviewer #2: Partly

**********

2. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: Yes

Reviewer #2: No

**********

3. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: No

Reviewer #2: No

**********

4. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

Reviewer #2: Yes

**********

5. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: In the manuscript, the authors collected pre-treatment tumor and/samples and on-treatment blood samples of NSCLC including squamous and non-squamous lung cancers (treated by anti-PD-1 immunotherapy) to generate transcriptomic data. The gene differentially expressed between responder and non-responders was identified for all three categories of samples. GSEA and ssGSEA analyses were then applied to identify and verify geneset/pathways associated with patient response to immunotherapy. For LUAD, the IFN and APP signatures were found to be associated with patient response based on gene expression data for pretreatment whole blood (WB) but not tumor samples. However, these signatures were not predictive for LUSC samples. Results from this analysis are quite interesting and of great clinical significance. My comments are below.

Major comments:

1. The combined pdf file was not correctly created. In the main text, legends for main figures and suppl. Figure/tables were mixed with result section, which make it hard to read. This comment is non-scientific but should be addressed.

2. Lack of independent validation. The gene sets (e.g., IFN) associated with patient response to immunotherapy were first identified by using GSEA (FGSEA) analysis. Then unsupervised method was used to calculate the Singscore as a measurement sample-specific gene set scores. It should be Singscore or ssGSEA score should not be considered as a validation of the GSEA results. Those scores for any gene sets identified from GSEA by comparing responder with non-responder are likely to be associated with patient response. To show more convincing evidence, application of these signature to independent datasets should be tested. Otherwise, the predictive values of these signatures cannot be secured.

3. Confounding variables. The methods used in these methods (cubic regression spline and survival comparison by Cox method) only examined the univariate effect of signature scores on PFS. It is unclear, whether the association between signature scores and PFS remain significant after adjusting confounding variables (e.g., age, PD-L1 level). This is particularly important since the samples size in this cohort is small -- only a small number of responders.

4. Integration of different signatures. It would be interesting to examine the joint prediction values of signatures identified in this study. Multivariate models should be used to integrate these signatures.

5. The gene expression data for all samples should be provided.

Minor comments:

1. More detailed information about the data availability will be useful. The RNAseq data are quite valuable.

2. The Table S1 “Baseline characteristics of patients in this cohort”, can be shown as a main table.

Reviewer #2: This is a single-institution, retrospective cohort study that examines 40 patients with a history of advanced non-small cell lung cancer (NSCLC) receiving nivolumab designed to analyze the transcriptomic features separately to identify histology-specific gene expression signatures that are associated with the clinical response to nivolumab monotherapy. Using multiple analytic approaches, you have found that the tumors of responders with lung adenocarcinoma (LUAD, n = 20) are inherently immunogenic to promote antitumor immunity, whereas those with lung squamous cell carcinoma (LUSC, n = 18) have a less immunosuppressive tumor microenvironment. You have concluded that nivolumab may function as a molecular targeted agent in LUAD and as an immunomodulating agent in LUSC. And, you have demonstrated that your study explains why the reliability of PD-L1 expression on tumor cells as a predictive biomarker for the response to nivolumab monotherapy is quite different between LUAD and LUSC.

You have provided descriptive and informative findings through the gene set enrichment analysis (GSEA), single-sample enrichment analysis (singscore), and spearman correlation analysis between the enrichment score and progression-free survival (PFS). Despite your appreciable efforts to analyze the transcriptome data, you need to improve the connections of your findings with the small number of patients with your conclusion.

1. Due to the small number of samples, it would be better to combine all the lung cancer data. You can consider the histology as a variable. Especially, the number of responders in each histology is too small to calculate the appropriate statistical power. It is not clear if the RNAseq data of whole blood need to be separately analyzed according to the histology during immunotherapy. It looks more interesting to suggest the independent factors to predict response to nivolumab through multivariable analysis.

2. It would be better to generate a schematic illustration about patient enrollment and sample collection.

3. The comparison of pathway analyses generated from bulk RNAseq data of pretreatment tumor with those of blood samples need to be clarified. During the comparison of RNAseq data of pretreatment tumor between responders and non-responders, you can consider tumor purity, immune deconvolution, cytolytic activity, and so on.

4. Your correlative findings are jumped to the confirmatory conclusion that the PD-1/PD-L1 immune evasion axis can be the primary target of nivolumab monotherapy in LUAD patients.

1) In page 14, you have described that the comparison between the single-sample enrichment scores of pretreatment WB and those of on-treatment WB indicated that type I IFN signaling was significantly enhanced in on-treatment WB of non-responders, but type I IFN signaling showed no significant enhancement in on-treatment compared with pretreatment WB of responders. And, you have demonstrated that nivolumab-induced activation of type I IFN signaling may be a predictive biomarker for worse clinical outcomes in LUAD patients treated with nivolumab monotherapy. However, the enhanced scores of on-treatment WB in non-responders are similar as those of pretreatment and on-treatment WB in responders. Furthermore, this finding is not consistent with previous GSEA analysis.

2) The enrichment scores of APP signatures exhibited no significant differences between pretreatment and on-treatment WB, indicating that nivolumab monotherapy has no impact on the immunogenicity of the tumor itself. Based on these findings, you presume that the tumors of responders are inherently sufficiently immunogenic to effectively elicit antigen processing and presentation for antitumor immunity. However, you need to provide some rationale to connect systemic immunity with local tumor immunity.

3) In page 12, in contrast to pretreatment WB of responders, pretreatment tumor tissues of responders did not show significant enrichment of IFN and APP signatures (Fig 1A and S4 Table). However, S4 Table shows GO_RESPONSE_TO_INTERFERON_GAMMA (p=0.00048, NES=1.6) and GO_INTERFERON_GAMMA_MEDIATED_SIGNALING_PATHWAY (p=0.0099, NES=1.52).

5. In Figure 3A, you have calculated the Spearman correlation coefficient (ρ) to estimate the correlation between the PFS time and the enrichment score. However, for survival comparison, it would be better to use Cox regression analysis considering time. To select the most powerful prognostic factor, it would be better to consider stepwise Cox regression analysis including hazard ratio. Your analytic approach did not support your findings that IFN and APP signatures in pretreatment WB may have robust predictive performance in LUAD patients treated with nivolumab monotherapy.

6. In page 18, you have described that the GSEA results indicated that the gene sets enriched in responders or non-responders were fairly similar between pretreatment and on-treatment WB in LUSC (Fig 4C–F, S14 Fig and S4 Table). However, figure legend of S14 Fig shows the lists of the top 25 gene sets significantly enriched in WB of responders and non-responders with LUAD. S14Fig shows only the non-responders’ results.

7. Please clarify the IRB number and study period.

8. All the RNAseq data reported in a submitted manuscript should be deposited in an appropriate public repository such as GEO or SRA.

9. Multiple tables are integrated in an excel file. Each supplementary table should be separately described.

**********

6. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email PLOS at figures@plos.org. Please note that Supporting Information files do not need this step.

23 Sep 2021

Response to the editor

We wish to express our appreciation to the editor for the insightful comments, which have helped us significantly improve our manuscript.

Comment 1: Please ensure that your manuscript meets PLOS ONE's style requirements, including those for file naming.

Response to Comment 1: We sincerely apologize that our original manuscript did not meet PLOS ONE's style requirements. According to the PLOS ONE style templates, we have corrected and revised the manuscript.

Comment 2: In order to ensure your study is reproducible, please provide the databases accession numbers included in your study.

Comment 4: We note that you have stated that you will provide repository information for your data at acceptance. Should your manuscript be accepted for publication, we will hold it until you provide the relevant accession numbers or DOIs necessary to access your data. If you wish to make changes to your Data Availability statement, please describe these changes in your cover letter and we will update your Data Availability statement to reflect the information you provide.

Response to Comment 2 and 4: We are very sorry that we have misunderstood the meaning of “Data Availability statement.” To be frank, we thought that it means that we must deposit our data to an appropriate public repository after acceptance for publication. Originally, we are going to deposit the RNA-seq data to DDBJ. However, our manuscript has not yet been published. Once our manuscript has been accepted for publication, we will deposit all the RNA-seq data. Thus, we have changed “Data Availability statement” from “Yes” to “No.” We kindly ask for your understanding.

Comment 3: We note that the grant information you provided in the ‘Funding Information’ and ‘Financial Disclosure’ sections do not match.

Response to Comment 1: Our study was funded by Ono Pharmaceutical Co., Ltd. and Bristol-Myers Squibb Co. The companies made a contract with our hospital directly to support our study only in terms of funding, and have had no input into the conception, conduct or reporting of our study.

We don’t have any grant number for the awards: instead, we have just a written contract. This contract was made in the name of Sendai Kousei Hospital, where the principal investigator of our study is Dr. Shunichi Sugawara.

Response to reviewers

We wish to express our appreciation to the reviewers for their insightful comments, which have helped us significantly improve our manuscript.

Reviewer #1

Major comment 1: The combined pdf file was not correctly created. In the main text, legends for main figures and suppl. Figure/tables were mixed with result section, which make it hard to read. This comment is non-scientific but should be addressed.

Response to Major comment 1: We sincerely apologize for the inconvenience caused this time. We regret that our original manuscript did not meet PLOS ONE's style requirements. According to the PLOS ONE style templates, we have corrected and revised the manuscript.

Major comment 2: Lack of independent validation. The gene sets (e.g., IFN) associated with patient response to immunotherapy were first identified by using GSEA (FGSEA) analysis. Then unsupervised method was used to calculate the Singscore as a measurement sample-specific gene set scores. It should be Singscore or ssGSEA score should not be considered as a validation of the GSEA results. Those scores for any gene sets identified from GSEA by comparing responder with non-responder are likely to be associated with patient response. To show more convincing evidence, application of these signature to independent datasets should be tested. Otherwise, the predictive values of these signatures cannot be secured.

Response to Major comment 2: We appreciate the reviewer’s comment. However, we have not considered Singscore as a validation of the classical GSEA results. The classical GSEA excels at identifying which predefined gene sets are behaving differently between the two groups defined. Due to the lack of phenotypic information (i.e., responders vs. nonresponders), we cannot employ the classical GSEA to predict whether a new and independent patient is a responder or not. In contrast, Singscore can calculate enrichment scores of predefined gene sets even in a new and independent patient. Moreover, the unsupervised enrichment scoring of interested gene sets enables regression modeling to predict PFS time. We screened positively or negatively enriched gene sets using the classical GSEA first, and then performed scoring of the positively or negatively enriched gene sets to apply the obtained enrichment scores to further analyses including correlation analysis, regression analysis and survival analysis.

As the reviewer pointed out, our study lacks validation in an independent cohort. We totally agree with this point, as we have emphasized in the Discussion part [page 25, line 13] that “One limitation of our study is the lack of validation in an independent cohort.” Unfortunately, we could not find any independent dataset, which contain both RNA-seq data from tumor tissue and whole blood and the associated clinical information of NSCLC patients treated with nivolumab monotherapy.

Major comment 3: Confounding variables. The methods used in these methods (cubic regression spline and survival comparison by Cox method) only examined the univariate effect of signature scores on PFS. It is unclear, whether the association between signature scores and PFS remain significant after adjusting confounding variables (e.g., age, PD-L1 level). This is particularly important since the samples size in this cohort is small -- only a small number of responders.

Major comment 4: Integration of different signatures. It would be interesting to examine the joint prediction values of signatures identified in this study. Multivariate models should be used to integrate these signatures.

Response to Major comment 3 and 4: We appreciate the reviewer’s comment. It is absolutely right, and we agree that it is unclear whether the association between signature scores and PFS remain significant after adjusting confounding variables. In response to the comment, we performed multivariate Cox regression analysis. As the results of the multivariate Cox regression analysis contain quite important data, we have included them in the main figure rather than in a supplementary figure. Also, we believe that the results would be the answer to “Major comment 4.”

As an aside, to predict PFS time, generalized additive models (GAM) using the 'gam' function in R (just as shown in Fig 3B and S10B Fig) were fit on the enrichment scores for two gene sets in the IFN and APP signatures. However, this attempt resulted in model overfitting. In our opinion, this is surely because the sample size in this cohort is too small to deal with two independent variables. We therefore removed the GAM data from the original manuscript.

Changes:

• We have divided the original “Figs 3A–D” into two figures (“Figs 3A–B” and “Figs 4A–B”) and have included the result from the multivariate Cox regression analysis in “Fig 4C”; the remaining figures have been renumbered accordingly.

• We have added the following text in the revised manuscript:

“Multivariate Cox regression analysis identified that many combinations of two of the IFN and APP signatures remained to be independent predictive factors of PFS time (Fig 4C and S14 Fig).” [page 17, lines 3–5]

• We have divided the original “Figs 6A–D” into two figures (“Figs 7A–B” and “Figs 8A–B”) and have included the result from the multivariate Cox regression analysis in “Fig 8C.”

• We have added the following text in the revised manuscript:

“In multivariate Cox regression analysis revealed that the COL1/TISSUE signature remained an independent predictive factor of PFS time (Fig 8). Other combinations of two of the TME signatures seemed to be slightly better predictive factors compared to clinical factors (age, ECOG PS, Brinkman index, PD-L1 TPS and CNS metastasis), but without statistical significance (S25 Fig).” [page 22, lines 21–25]

• We have added two new supplementary figures, “S14 Fig” and “S25 Fig.”

• We have added the following text in “Materials and methods - Survival analysis” part of the revised manuscript:

“Cox proportional hazard models were built using the “coxph” function from the survival package and visualized as forest plots using ggplot2.” [page 9, lines 17–19]

Major comment 5: The gene expression data for all samples should be provided.

Response to Major comment 5: We agree with the reviewer’s advice and have provided the raw count data of RNA-seq as “S19 Table.”

Minor comments 1: More detailed information about the data availability will be useful. The RNAseq data are quite valuable.

Response to Minor comment 1: We appreciate and greatly agree with the reviewer’s comment. That’s the plan, and we are going to deposit the RNA-seq data to DDBJ. However, our manuscript has not yet been published. Once our manuscript has been accepted for publication, we will deposit all the RNA-seq data. We kindly ask for the reviewer’s understanding.

Minor comments 2: The Table S1 “Baseline characteristics of patients in this cohort”, can be shown as a main table.

Response to Minor comment 2: We appreciate the reviewer's suggestion. As the reviewer suggested, we have included the original “Table S1” in the main manuscript as “Table 1.”

Changes:

• The original “S1 Table” has been renumbered as “Table 1”; the remaining supplementary tables have been renumbered accordingly.

Reviewer #2

Major comment 1: Due to the small number of samples, it would be better to combine all the lung cancer data. You can consider the histology as a variable. Especially, the number of responders in each histology is too small to calculate the appropriate statistical power. It is not clear if the RNAseq data of whole blood need to be separately analyzed according to the histology during immunotherapy. It looks more interesting to suggest the independent factors to predict response to nivolumab through multivariable analysis.

Response to Major comment 1: We appreciate the reviewer’s comment. We understand that the reviewer wonders why we further divided the relatively small number of samples into two smaller groups based on histology. “Due to the small number of samples, it would be better to combine all the lung cancer data” is exactly what we were thinking. In fact, we first performed classical GSEA on all the NSCLC data, which identified that only a few gene sets were significantly enriched in responders or nonresponders. For this reason, we analyzed LUAD and LUSC separately and found that baseline transcriptomic features in responders or nonresponders are clearly different between LUAD and LUSC.

Major comment 2: It would be better to generate a schematic illustration about patient enrollment and sample collection.

Response to Major comment 2: We thank the reviewer for this suggestion. In respond to the suggestion, we have included a new “S1 Fig.”

Changes:

• We have added a new “S1 Fig”; the remaining supplementary figures have been renumbered accordingly.

Major comment 3: The comparison of pathway analyses generated from bulk RNAseq data of pretreatment tumor with those of blood samples need to be clarified. During the comparison of RNAseq data of pretreatment tumor between responders and non-responders, you can consider tumor purity, immune deconvolution, cytolytic activity, and so on.

Response to Major comment 3: We appreciate the reviewer for raising this very interesting point. Only tumor tissue samples with a tumor content of more than 30% were included. We believe that this reflects the reality in daily clinical practice.

Additionally, a previous study has provided a metric for immune cytolytic activity based on gene expression in tumors, where immune cytolytic activity was estimated by the average expression level of CD8A, CD8B, GZMA, GZMB and PRFref.1. Using this metric, we assessed immune cytolytic activity in tumor tissues from patients in our cohort. As the result, we found no significant difference in PFS and OS between high and low immune cytolytic activity. We have included the result as a new “S4 Fig”.

ref.1 Rooney MS, et al. Molecular and genetic properties of tumors associated with local immune cytolytic activity. Cell 2015;160: 48–61.

Changes:

• We have included a new “S4 Fig”; the remaining supplementary figures have been renumbered accordingly.

• We have added the following text in the revised manuscript:

“A previous study has provided a metric for immune cytolytic activity based on gene expression in tumors, where immune cytolytic activity was estimated by the average expression level of CD8A, CD8B, GZMA, GZMB and PRF [27]. Using this metric, we assessed immune cytolytic activity in tumor tissues from patients in our cohort. As the result, we found no significant association between the immune cytolytic activity and clinical outcomes (S4 Fig).” [page 10, lines 13–18]

• We have included a new reference [27]; the remaining references have been renumbered accordingly.

Major comment 4: Your correlative findings are jumped to the confirmatory conclusion that the PD-1/PD-L1 immune evasion axis can be the primary target of nivolumab monotherapy in LUAD patients.

1) In page 14, you have described that the comparison between the single-sample enrichment scores of pretreatment WB and those of on-treatment WB indicated that type I IFN signaling was significantly enhanced in on-treatment WB of non-responders, but type I IFN signaling showed no significant enhancement in on-treatment compared with pretreatment WB of responders. And, you have demonstrated that nivolumab-induced activation of type I IFN signaling may be a predictive biomarker for worse clinical outcomes in LUAD patients treated with nivolumab monotherapy. However, the enhanced scores of on-treatment WB in non-responders are similar as those of pretreatment and on-treatment WB in responders. Furthermore, this finding is not consistent with previous GSEA analysis.

2) The enrichment scores of APP signatures exhibited no significant differences between pretreatment and on-treatment WB, indicating that nivolumab monotherapy has no impact on the immunogenicity of the tumor itself. Based on these findings, you presume that the tumors of responders are inherently sufficiently immunogenic to effectively elicit antigen processing and presentation for antitumor immunity. However, you need to provide some rationale to connect systemic immunity with local tumor immunity.

3) In page 12, in contrast to pretreatment WB of responders, pretreatment tumor tissues of responders did not show significant enrichment of IFN and APP signatures (Fig 1A and S4 Table). However, S4 Table shows GO_RESPONSE_TO_INTERFERON_　GAMMA (p=0.00048, NES=1.6) and GO_INTERFERON_GAMMA_MEDIATED_　SIGNALING_PATHWAY (p=0.0099, NES=1.52).

Response to Major comment 4: We appreciate the reviewer’s insightful questions and constructive suggestions in comment 1) – 3).

Response to comment 1): We first apologize for our confusing description. The reviewer has pointed out that “the enhanced scores of on-treatment WB in non-responders are similar as those of pretreatment and on-treatment WB in responders focuses on the enrichment scores.” However, it is our intention to describe the following:

i. In nonresponders (NR), the enrichment scores of type I IFN signaling in on-treatment WB were significantly increased compared with pretreatment WB.

ii. In responders (R), the enrichment scores of type I IFN signaling in on-treatment WB showed no significant differences compared with pretreatment WB.

iii. In pretreatment WB, the enrichment scores of type I IFN signaling in R were significantly increased compared with NR.

Here, we have focused on the difference in the effect of nivolumab on type I IFN signaling between R and NR. The activation of type I IFN signaling after nivolumab monotherapy was observed only in WB from NR, not in WB from R. We believe that these findings are consistent with the results of our classical GSEA. We hope that the reviewer will agree with us on this point.

Response to comment 2): We agree with the reviewer that “we need to provide some rationale to connect systemic immunity with local tumor immunity.” We should have provided a more in-depth description of the point. As neoantigens are generated from mutations, the higher the TMB, the greater the chance that some of the neoantigens presented by MHC proteins will be immunogenic and hence enable the induction of anti-tumor immune response. In fact, accumulating evidence has indicated that high TMB is correlated with better clinical outcomes in NSCLC patients with anti-PD-1/PD-L1 therapyref.2. Additionally, several studies have reported that tumor mutation burden (TMB) in blood (or circulating tumor DNA) correlates with TMB in tumor tissue, and that high TMB in blood may serve as a potential biomarker of clinical benefit in NSCLC patients with anti-PD-1/PD-L1 therapyref.3, 4. Given that TMB in blood can be a surrogate for TMB in tumor tissue, it is rational to suggest that TMB in blood reflect the immunogenicity of tumor cells. We have added this description and relevant references to the manuscript.

ref.2 Sha D, et al. Tumor mutational burden as a predictive biomarker in solid tumors. Cancer Discov 2020;10: 1808–1825.

ref.3 Gandara DR, et al. Blood-based tumor mutational burden as a predictor of clinical benefit in non-small-cell lung cancer patients treated with atezolizumab. Nat Med 2018;24: 1441–1448.

ref.4 Wang Z, et al. Assessment of blood tumor mutational burden as a potential biomarker for Immunotherapy in patients with non-small cell lung cancer with use of a next-generation sequencing cancer gene panel. JAMA Oncol 201;5: 696–702.

Response to 3) As mentioned in “Materials and methods,” we defined significantly enriched gene sets as those with an adjusted p-value < 0.1. The adjusted p-values are 0.01558 and 0.12613 for “GO_RESPONSE_TO_INTERFERON_GAMMA” and “GO_ INTERFERON_GAMMA_MEDIATED_SIGNALING_PATHWAY” respectively. Due to the discrepancy between the IFN-γ signaling gene sets, we concluded that “pretreatment tumor tissues of responders did not show significant enrichment of IFN and APP signatures.” However, we acknowledge that it seems misleading. Our attention here is directed on the magnitude of the enrichment of IFN and APP signatures, not just statistical significance. Hence, we have revised the sentence that the reviewer pointed out.

Changes:

• To better clarify what we meant, we have changed the sentences that the reviewer pointed out in comment 1) from:

“Comparison between the single-sample enrichment scores of pretreatment WB and those of on-treatment WB indicated that type I IFN signaling was significantly enhanced in on-treatment WB of nonresponders (‘IFN_I’, p = 0.0128; ‘IFNB1’, p = 0.0108; ‘IFNB2’, p = 0.0025) (Fig 2B). In contrast, type I IFN signaling showed no significant enhancement in on-treatment compared with pretreatment WB of responders (Fig 2B).”

to

“In nonresponders, the enrichment scores of type I IFN signaling in on-treatment WB were significantly increased compared with pretreatment WB (‘IFN_I’, p = 0.0128; ‘IFNB1’, p = 0.0108; ‘IFNB2’, p = 0.0025) (Fig 2B). In responders, by contrast, the enrichment scores of type I IFN signaling in on-treatment WB showed no significant differences compared with pretreatment WB (Fig 2B).” [page 15, lines 5–6]

• In response to the reviewer’s comment 2), we have changed the following text from:

“Moreover, the enrichment scores of APP signatures exhibited no significant differences between pretreatment and on-treatment WB, indicating that nivolumab monotherapy has no impact on the immunogenicity of the tumor itself.”

to

“Moreover, we noted that the enrichment scores of APP signatures exhibited no significant differences between pretreatment and on-treatment WB. As neoantigens are generated from mutations, the higher the TMB, the greater the chance that some of the neoantigens presented by MHC proteins will be immunogenic and hence enable the induction of anti-tumor immune response. In fact, accumulating evidence has indicated that high TMB is correlated with better clinical outcomes in NSCLC patients with anti-PD-1/PD-L1 therapy [41]. Additionally, several studies have reported that tumor mutation burden (TMB) in blood (or circulating tumor DNA) correlates with TMB in tumor tissue, and that high TMB in blood may serve as a potential biomarker of clinical benefit in NSCLC patients with anti-PD-1/PD-L1 therapy [42, 43]. Given that TMB in blood can be a surrogate for TMB in tumor tissue, it is rational to suggest that TMB in blood reflect the immunogenicity of tumor cells.

Therefore, the comparable levels of APP signatures before and after nivolumab monotherapy indicate that nivolumab monotherapy has no significant impact on the immunogenicity of the tumor itself.” [page 15, lines 11–25]

• We have included new references [41]–[43] to support the above description; the remaining references have been renumbered accordingly:

• In response to the reviewer’s comment 3), we have changed the following sentence from:

“In contrast to pretreatment WB of responders, pretreatment tumor tissues of responders did not show significant enrichment of IFN and APP signatures (Fig 1A and S5 Table).”

to

“In contrast to pretreatment WB of responders, pretreatment tumor tissues of responders did not show strong enrichment of IFN and APP signatures (Fig 1A and S5 Table).” [page 13, lines 18–20]

• Finally, to avoid exaggeration and overstatement, we have changed the following text from:

“Thus, the PD-1/PD-L1 immune evasion axis can be the primary target of nivolumab monotherapy in LUAD patients.”

to

“Thus, the PD-1/PD-L1 immune evasion axis can be one of the primary targets of nivolumab monotherapy in LUAD patients.” [page 15, lines 30–31]

Major comment 5: In Figure 3A, you have calculated the Spearman correlation coefficient (ρ) to estimate the correlation between the PFS time and the enrichment score. However, for survival comparison, it would be better to use Cox regression analysis considering time. To select the most powerful prognostic factor, it would be better to consider stepwise Cox regression analysis including hazard ratio. Your analytic approach did not support your findings that IFN and APP signatures in pretreatment WB may have robust predictive performance in LUAD patients treated with nivolumab monotherapy.

Response to Major comment 5: We appreciate and completely agree with the reviewer’s insightful comment. The reviewer #1 has also raised this point. Please see response to “Major comment 3 and 4” from the reviewer #1

Major comment 6: In page 18, you have described that the GSEA results indicated that the gene sets enriched in responders or non-responders were fairly similar between pretreatment and on-treatment WB in LUSC (Fig 4C–F, S14 Fig and S4 Table). However, figure legend of S14 Fig shows the lists of the top 25 gene sets significantly enriched in WB of responders and non-responders with LUAD. S14Fig shows only the non-responders’ results.

Response to Major comment 6: We appreciate the reviewer for pointing out the mistake. We sincerely apologize for our careless mistake and have corrected it in the original “S14 Fig”. The corrected figure has been included in a new “S17 Fig.”

Major comment 7: Please clarify the IRB number and study period.

Response to Major comment 7: The IRB number 29-4. The study period is from 18 Jul 2017 to 31 Dec 2020. We have added the above information to the revised manuscript. [page 5, lines 5–6]

Major comment 8: All the RNAseq data reported in a submitted manuscript should be deposited in an appropriate public repository such as GEO or SRA.

Response to Major comment 8: We appreciate the reviewer’s comment. The reviewer #1 has also provided the same advice. Please see response to “Minor comment 1” from the reviewer #1.

Major comment 9: Multiple tables are integrated in an excel file. Each supplementary table should be separately described.

Response to Major comment 9: We appreciate the reviewer’s advice. In response to the advice, we have separated the supplementary tables.

Changes:

• The original “S2 Table” has been separated into the following new tables:

“S2 Table” (The lists of DEGs identified from pretreatment tumor tissue of LUAD patients), “S3 Table” (The lists of DEGs identified from pretreatment whole blood of LUAD patients), “S4 Table” (The lists of DEGs identified from on-treatment whole blood of LUAD patients), “S11 Table” (The lists of DEGs identified from pretreatment tumor tissue of LUSC patients), and “S12 Table” (The lists of DEGs identified from on-treatment whole blood of LUSC patients).

• The original “S3 Table” has been separated into the following new tables:

“S5 Table” (The results from classical GSEA of pretreatment tumor tissue of LUAD patients), “S6 Table” (The results from classical GSEA of pretreatment whole blood of LUAD patients), “S7 Table2” (The results from classical GSEA of on-treatment whole blood of LUAD patients), “S13 Table” (The results from classical GSEA of pretreatment tumor tissue of LUSC patients), “S14 Table” (The results from classical GSEA of pretreatment whole blood of LUSC patients), and “S15 Table” (The results from classical GSEA of on-treatment whole blood of LUSC patients).

• The original “S4 Table” has been separated into the following new tables:

“S8 Table” (The results from single-sample enrichment analysis of pretreatment tumor tissue of LUAD patients), “S9 Table” (The results from single-sample enrichment analysis of pretreatment whole blood of LUAD patients), and “S10 Table” (The results from single-sample enrichment analysis of on-treatment whole blood of LUAD patients).

Submitted filename: Response to Reviewers.docx

Click here for additional data file.

11 Nov 2021

Gene expression signatures as candidate biomarkers of response to PD-1 blockade in non-small cell lung cancers.

PONE-D-21-19886R1

Dear Dr. aiba,

We’re pleased to inform you that your manuscript has been judged scientifically suitable for publication and will be formally accepted for publication once it meets all outstanding technical requirements.

Within one week, you’ll receive an e-mail detailing the required amendments. When these have been addressed, you’ll receive a formal acceptance letter and your manuscript will be scheduled for publication.

An invoice for payment will follow shortly after the formal acceptance. To ensure an efficient process, please log into Editorial Manager at http://www.editorialmanager.com/pone/, click the 'Update My Information' link at the top of the page, and double check that your user information is up-to-date. If you have any billing related questions, please contact our Author Billing department directly at authorbilling@plos.org.

If your institution or institutions have a press office, please notify them about your upcoming paper to help maximize its impact. If they’ll be preparing press materials, please inform our press team as soon as possible -- no later than 48 hours after receiving the formal acceptance. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information, please contact onepress@plos.org.

Kind regards,

Hyun-Sung Lee, M.D., Ph.D.

Academic Editor

PLOS ONE

Additional Editor Comments (optional):

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. If the authors have adequately addressed your comments raised in a previous round of review and you feel that this manuscript is now acceptable for publication, you may indicate that here to bypass the “Comments to the Author” section, enter your conflict of interest statement in the “Confidential to Editor” section, and submit your "Accept" recommendation.

Reviewer #1: All comments have been addressed

**********

2. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Yes

**********

3. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: Yes

**********

4. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

**********

5. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

**********

6. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: All my comments to the previous version of the manuscript have been fully addressed in the revision. I don't have additional comments.

**********

7. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

16 Nov 2021

PONE-D-21-19886R1

Gene expression signatures as candidate biomarkers of response to PD-1 blockade in non-small cell lung cancers.

Dear Dr. aiba:

I'm pleased to inform you that your manuscript has been deemed suitable for publication in PLOS ONE. Congratulations! Your manuscript is now with our production department.

If your institution or institutions have a press office, please let them know about your upcoming paper now to help maximize its impact. If they'll be preparing press materials, please inform our press team within the next 48 hours. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information please contact onepress@plos.org.

If we can help with anything else, please email us at plosone@plos.org.

Thank you for submitting your work to PLOS ONE and supporting open access.

Kind regards,

PLOS ONE Editorial Office Staff

on behalf of

Dr. Hyun-Sung Lee

Academic Editor

PLOS ONE