Identification of Potential Prognostic and Predictive Biomarkers for Immune-Checkpoint Inhibitor Response in Small Cell Lung Cancer

Background

Small cell lung cancer (SCLC) is one of the deadliest and most aggressive lung cancers due to its rapid growth, early metastasis, and acquired therapeutic resistance [1]. Unlike non-small cell lung cancer (NSCLC), which can be treated with a plethora of drugs, few treatment options are available for patients with SCLC. The current standard of care, which was defined several decades ago [2], relies on platinum-based chemotherapy and has a median survival time of only 8–13 months [3].

However, the treatment landscape has recently changed with the introduction of immune-checkpoint inhibitors (ICIs). Incorporation of the anti-PDL1 antibody atezolizumab with platinum-doublet chemotherapy has become the first-line therapy for extensive SCLC, with an improvement in median overall survival (OS) from 10.3 to 12 months for patients with metastatic disease [4]. Ipilimumab (an anti-CTLA4 monoclonal antibody) also showed increased survival in a phase II trial [5].

Unfortunately, only a subset of patients with SCLC benefits from immunotherapy [4]. More recent studies have reported the use of drugs such as CheckMate-032 and CheckMate-331 for maintenance and second-line therapy, but these fail to improve survival for patients with SCLC [6,7]. Recent research suggest that high tumor mutational burden (TMB) may be correlated with efficacy of immunotherapy in SCLC [8]. Moreover, no study used genetic marker to select patients more likely to benefit from these treatments, which may explain their clinical failure. Therefore, there is an urgent need to identify biomarkers that can predict the responsiveness to immune-checkpoint blockade treatment.

To address this, we analyzed the expression of different genes in healthy and tumor tissues from patients with SCLC in 2 independent cohorts. We then identified genes associated with immunotherapeutic response. Using subclass mapping and survival analysis, we found 6 prognostic genes predictive for response to immune-checkpoint inhibitors. These results were validated by meta-analysis and a third cohort (a summary of the statistical analysis is presented in Figure 1).

Material and Methods

DATA AND RESOURCES:

We used 3 different cohorts for this study. Expression data for the first 2 cohorts were downloaded from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo) [9]. The first cohort was obtained from the Affymetrix U133plus2 GEO SCLC Cohort (accession number GSE43346) and comprised 42 normal samples and 23 clinical SCLC samples. The second cohort downloaded from GEO (accession number GSE60052) consisted of 7 normal lung samples and 79 SCLC samples, of which only 48 had associated clinical data. Gene expression profiles were obtained using the Illumina HiSeq2000 platform. The third cohort was downloaded from the Cancer Cell Line Encyclopedia (CCLE) database (https://portals.broadinstitute.org/ccle) [10]. The CCLE datasets contain information from cell lines derived from different tumor types; we used data from SCLC tumors, resulting in 51 samples for this cohort.

Additionally, we analyzed a comprehensive list of 5826 immunotherapy-related genes available from the Tumor-Immune System Interactions Database (TISIDB) (http://cis.hku.hk/TISIDB/) [11]. These genes are known to differ in expression between patients that are responders vs those who are non-responders to immunotherapy.

IDENTIFICATION OF DIFFERENTIALLY EXPRESSED GENES RELATED TO IMMUNOTHERAPY:

We identified differentially expressed genes (DEGs) by comparing the normalized expression data of SCLC samples to that of normal samples from GSE43346 and GSE60052. Only DEGs with |log2FC| >2 and FDR <0.05 were extracted using the Limma package in R [12]. We calculated the difference between the 2 groups using the t test and adjusted the P value using the Benjamini-Hochberg (BH) method. To identify immunotherapy-related genes among the SCLC DEGs, Venn analysis was utilized to determine overlapping DEGs among the 2 GEO datasets and the TISIDB dataset above. To narrow the scope of candidate genes, we selected the top 300 DEGs, in ascending order of FDR, from each GEO dataset.

FUNCTIONAL ENRICHMENT ANALYSIS AND PPI NETWORK CONSTRUCTION:

We conducted function and protein-protein interaction (PPI) enrichment analyses to predict the biological function of the aforementioned DEGs using Metascape online tools (http://metascape.org/gp/index.html#/main/step1) [13]. The DEGs were annotated according to gene ontology biological processes. The P value was set at <0.01, indicating a statistically significant enrichment score of the GO or pathway terms. The PPI network was also visualized by Metascape.

SURVIVAL ANALYSIS:

Fifty-five patients (48 SCLC samples and 7 normal lung samples) from GSE60052 were divided into high- and low-expression groups based on optimal cut-off values defined by R package “Survminer” [14] for candidate gene expression. We used Kaplan-Meier survival curves to show the differences in patients’ OS between the high and low expression groups. Genes with P values <0.05 were identified as prognosis-related genes.

IMMUNOTHERAPEUTIC/CHEMOTHERAPEUTIC RESPONSE PREDICTION:

We divided the expression profiles from GSE60052 (n=86) into high- and low-expression groups based on the median expression values of our key genes. We further utilized an unsupervised subclass mapping (https://cloud.genepattern.org/gp/) to compare these expression profiles with those of another dataset containing 47 patients with melanoma who responded to immunotherapy [15,16] to predict the response to ICI therapy in high- and low-subgroup. R package “pheatmap” was used for visualizing the results. We also analyzed the expression profiles of SCLC cell lines (n=51) from the CCLE database using the same method to validate the immunotherapeutic response prediction. Groups with Bonferroni-corrected P values <0.05 were identified as immunotherapy responders.

In addition to the immunotherapeutic response prediction, the chemotherapeutic response for 86 samples was predicted using the Genomics of Drug Sensitivity in Cancer database (GDSC) (https://www.cancerrxgene.org/) [17]. We predicted the responses to cisplatin and etoposide using the pRRophetic package in R [18]. Each sample’s half-maximal inhibitory concentration (IC50) was assessed by ridge regression, and the prediction accuracy was evaluated by 10-fold cross-validation relying on the GDSC training model. We summarized duplicate gene expression as the mean expression, and set the tissue type to “allSolidTumours” and the batch effect to “combat”. Meanwhile, the values of other parameters were set as default.

SURVIVAL META-ANALYSIS OF KEY HENES:

Using survival meta-analysis on data from the Lung Cancer Explorer (LCE) database (http://lce.biohpc.swmed.edu/lungcancer/) [19], we were able to validate whether our chosen key genes were prognosis-related in SCLC based on the meta-analysis module in LCE database. Survival meta-analysis was based on 56 studies, including over 6700 patients with lung cancer. Cox hazard models were fitted to evaluate the association between key genes expression and patient survival outcome in each dataset. After heterogeneity testing, the HR from individual datasets were then calculated to summary HR utilizing a random effects model.

STATISTICAL AND BIOINFORMATIC ANALYSIS:

The t test was used to compare gene expression between tumor tissues and adjacent nontumorous tissues. The OS between different groups was compared by Kaplan-Meier analysis with the log-rank test. Bonferroni correction was used for adjust P value in subclass mapping. The Wilcox test was used to compare the IC50 between different subgroups. R software (Version 4.0.0, http://www.r-project.org), Metascape online tool, subclass mapping and LCE database were used to analyze data and plot graphs. If not specified above, a P value less than 0.05 was considered statistically significant.

Results

IDENTIFICATION OF CANDIDATE GENES RELATED TO IMMUNOTHERAPY:

To select genes differentially expressed between SCLC and normal tissues, we conducted Limma analysis of data from GSE43346 and GSE60052. We identified 2951 and 3365 DEGs from each dataset, respectively, and the 2 datasets had 777 DEGs in common (co-DEGs). We then identified common genes between these co-DEGS and the immunotherapy-related genes from TISIDB, which included DEGs between immunotherapy responders and non-responders, resulting in 443 candidate genes retained for further analysis (Figure 2A).

To obtain further insight into their function, the candidate genes were uploaded to the Metascape database. GO and KEGG (Gene Ontology and Kyoto Encyclopedia of Genes and Genomes) analysis results showed that the candidate genes were mainly enriched through mitotic nuclear division, chromosome segregation, nuclear division, organelle fission, nuclear chromosome segregation, and regulation of chromosome segregation. Moreover, candidate genes for organelle fission were significantly enriched (Figure 2B).

Using the Metascape database, we also constructed PPI networks selecting samples that had a PPI minimum network size above 3. From these constructed PPI networks, we identified closely related proteins using molecular complex detection (MCODE). The candidate genes in the PPI networks were stratified into 11 subclusters through the MCODE algorithm. Figure 2C shows MCODE components.

KEY GENES AS PROGNOSTIC BIOMARKERS:

The prognostic importance of the 11 selected genes modulated in SCLC was evaluated using the GSE60052 clinical data, from which we retrieved 55 SCLC patient samples with clinical data. The gene expression values were categorized as high or low according to the optimal cut-off values set by the Survminer package for R. As a result, 6 genes were found to be independent prognostic factors for OS. Low expression of enhancer of zeste homolog 2 (EZH2) was associated with poor OS (P=0.049), whereas high expression of minichromosome maintenance complex component 2 (MCM2) (P=0.031), cyclin-dependent kinase inhibitor 2A (CDKN2A) (p=.008), centromere protein K (CENPK) (P<0.001), checkpoint kinase 1 (CHEK1) (P=0.001), and exosome component 2 (EXOSC2) (P=0.031) were associated with poor OS (Figure 3A). The heatmap showed that the expression level of 6 genes is lower in the normal tissue (Figure 3B). Moreover, the correlation heatmap showed a significant positive correlation among the 6 genes (Figure 3C).

KEY GENES AS PREDICTIVE BIOMARKERS OF IMMUNOTHERAPEUTIC/CHEMOTHERAPEUTIC RESPONSE:

We subsequently performed a SubClass Mapping (SubMap) analysis to determine whether groups with high and low expression of these key genes were similar to other groups that clearly responded to immunotherapy [15]. SubMap results indeed showed that low expression of the 6 key genes was correlated with favorable response to anti-PD-1 treatment CDKN2A (Bonferroni-corrected P=0.008), CENPK (Bonferroni-corrected P=0.007), CHEK1 (Bonferroni-corrected P=0.023), EXOSC2 (Bonferroni-corrected P=0.001), EZH2 (Bonferroni-corrected P=0.031), and MCM2 (Bonferroni-corrected P=0.008). Moreover, lower expression of EXOSC2 was correlated with more favorable response to anti-CTLA4 treatment (Bonferroni-corrected p=0.039) (Figure 4).

Because chemotherapy is the standard of care for SCLC, we assessed the GDSC data for chemotherapeutic responses to cisplatin and etoposide based on a predictive model. The estimated IC50 of cisplatin, but not of etoposide, significantly differed between the high- and low-expression groups. Low expression of MCM2 (P=0.019), EZH2 (p=0.0008), CENPK (P=0.008), and CHEK1 (P=0.0013) was correlated with favorable cisplatin response, whereas the expression of CDKN2A (p=0.598) and EXOSC2 (P=0.096) was not predictive of chemotherapeutic response. There were no statistically significant differences between the high and low expression groups in terms of predicting response to etoposide (MCM2, P=0.824; EZH2, P=0.473; CENPK, P=0.775; CHEK1, P=0.832; CDKN2A, P=0.993; and EXOSC2, P=0.924). Overall, our results suggest that MCM2, EZH2, CENPK, and CHEK1 are predictive biomarkers for chemotherapeutic response to cisplatin (Figure 5).

VALIDATION OF THE PROGNOSTIC AND PREDICTIVE VALUE OF THE 6 KEY GENES:

To validate the expression of the 6 potential biomarkers for the prognosis of SCLC, we examined their impact on the OS of patients with SCLC. In multiple selected lung cancer studies, we found that the survival meta-analyses for MCM2, EZH2, CENPK, CHEK1, CDKN2A, and EXOSC2 were significant in patients with SCLC (P<0.01). Our results indicate that these 6 genes can predict the prognosis of SCLC (Figure 6).

Finally, we selected 51 SCLC cell lines to validate the ability of our 6 key genes to predict immunotherapeutic response. Using SubMap analysis, we found that low expression of 4 of the 6 genes predicted response to anti-PD-1 therapy (CENPK (Bonferroni-corrected P=0.020), CHEK1 (Bonferroni-corrected P=0.021), and EZH2 (Bonferroni-corrected P=0.048)). Low expression of MCM2, CDKN2A, and EXOSC2 did not significantly predict response (MCM2 Bonferroni-corrected P=0.479, CDKN2A Bonferroni-corrected P=1.000, and EXOSC2 Bonferroni-corrected P= 0.315). Interestingly, the high expression of CENPK predicted response to anti-CTLA4 therapy (CENPK Bonferroni-corrected P=0.030) (Figure 7).

Discussion

The addition of ICIs to chemotherapy regimens for first-line management is a significant advancement in the treatment of patients with SCLC. Predictive biomarkers are crucial to identify patients with SCLC who could benefit from ICIs. Recent studies have shown that TMB could serve as a predictive factor for immunotherapy in SCLC [20]. However, genetic markers that predict ICI response have not been well studied. To the best of our knowledge, this is the first study to assess potential genes that predict response to ICIs in SCLC. In this study, we also assessed chemotherapeutic response in patients with SCLC. Our results revealed that lower expression of 6 key genes (MCM2, EZH2, CDKN2A, CENPK, CHEK1, and EXOSC2) in patients with SCLC predicted responses to both anti-PD-1 therapy and cisplatin. We expect that low expression of these 6 key genes will be a significant predictor of clinical response to the combination of chemotherapy and ICIs in patients with SCLC, and the expression levels of 6 key genes were lower in the adjacent normal tissues. However, due to the lack of detailed clinical information, the expression level of the key genes was unclear in different stages. Thus, further studies are urgently needed to explore the correlation of these 6 key genes and tumor progression.

A limitation of this study is that we did not explore the potential mechanisms of our 6 key genes underlying the response of patients with SCLC to immunotherapy, as all the data analyzed in our study were retrieved from online databases. The prognostic significance of the 6 genes is validated in lung cancers, not specifically in small cell lung cancer. Due to the lack of valid verifying datasets of SCLC, further studies with larger sample sizes are required to validate our findings and to explore the detailed mechanisms of the 6 key genes in SCLC.

MCM2, which is essential for DNA replication and cell cycle, has been recently identified as a prognostic biomarker in oral, gastric, colon, and breast cancers [21–24]. Fujii et al reported that MCM2 could be a strong predictor of poor outcome in SCLC [25]. These studies suggest that high expression of MCM2 may play a role in the poor prognosis of cancers, which is consistent with our results [26]. However, the specific immunotherapeutic function of MCM2 in SCLC has not yet been explored. In our study, low MCM2 expression was correlated with increased response to anti-PD-1 treatment in the GEO datasets, but this was not statistically significant in the CCLE dataset. Further studies are needed to determine whether MCM2 plays a vital role in immunotherapeutic response of SCLC.

EZH2 is highly expressed in various cancer types and is more frequently overexpressed in SCLC than in NSCLC [27]. Hubaux et al demonstrated that EZH2 inhibits apoptosis and promotes the cell cycle in SCLC by promoting the oncogenic RB1/E2F pathway [28]. Gardner reported that inhibition of EZH2 enhances the effectiveness of current standard chemotherapy [29]. Toyokawa et al found that the expression of EZH2 was not significantly associated with postoperative survival [30]. In SCLC, EZH2, and other PRC2 components contribute to gene repression, which may be associated with worse prognosis in SCLC [31]. Intriguingly, we found a correlation between low EZH2 expression and poor OS in our SCLC cohorts.

EZH2 is further reported to inhibit immunity by repressing the production of Th1-type chemokines CXCL9 and CXCL10 in colon and ovarian carcinomas [32,33]. Additionally, Bugide et al reported that inhibition of EZH2 enhances hepatocellular carcinoma eradication by NK cells [34]. In our study, low EZH2 expression was correlated with a favorable response to anti-PD-1 treatment, validated by 2 independent SCLC cohorts. Nevertheless, the potential mechanism of EZH2 in the regulation of SCLC immune responses remains unclear.

Hans et al reported that CDKN2A is a major driver of pancreatic carcinogenesis in patients with pancreatic ductal adenocarcinoma (PDAC) [35]. However, the prognostic role of CDKN2A in SCLC is yet to be investigated. In this study, we demonstrated that lower CDKN2A expression was correlated with poor OS and response to anti-PD-1 therapy, but was not correlated with response to cisplatin.

To date, little is known about the expression and roles of CENPK, CHEK1, or EXOSC2 in SCLC. Lee et al reported that overexpression of CENPK in ovarian cancer is correlated with poor patient survival [36]. Liu et al discovered that high expression of CHEK1 in NSCLC is associated with poor OS [37]. Other studies have found that CHEK1 inhibitors may affect sensitivity to radiotherapy and chemotherapy [38–40]. In our study, low CENPK, CHEK1, and EXOSC2 expression was correlated with both poor OS and favorable response to anti-PD-1 therapy. Low CENPK and CHEK1 expression was correlated with response to cisplatin, but EXOSC2 expression was not significantly correlated with this response.

Conclusions

In conclusion, our results showed that expression of the 6 key genes was significantly associated with OS of patients with SCLC, and that lower expression of these 6 genes was correlated with response to anti-PD-1 therapy. We also showed that low expression of MCM2, EZH2, CENPK, and CHEK1 could predict response to cisplatin, but not etoposide, indicating that MCM2, EZH2, CENPK, CHEK1, CDKN2A, and EXOSC2 could be prognostic and predictive biomarkers for response to ICIs in patients with SCLC.