We recruited histologically confirmed SCLC patients from the Shandong Cancer Hospital and Institute (SCH). All diagnoses were independently confirmed by two experienced pathologists. In addition to blood samples (2 ml), tumor tissue samples were obtained by biopsies. A strict quality inspection was carried out to remove contaminated and insufficient DNA samples. Finally, 177 patients were enrolled in our study. The overall survival (OS) time was defined as the interval between diagnosis and death, or between diagnosis and the last observation point. For surviving patients, data were censored at the last follow-up (November 26, 2020). Clinicopathological data were retrieved from the patients’ medical records. This study was approved by the Ethics Committee of the Shandong Cancer Hospital and Institute. All patients included in this study provided written informed consent.

Biopsied tumor tissues were fixed in formalin and then embedded in paraffin (FFPE). The corresponding blood samples were set as controls. Genomic DNA was extracted from each FFPE sample using the GeneRead DNA FFPE Kit (Qiagen, USA) and from the blood sample using the DNA Blood Midi/Mini Kit (Qiagen, USA).

Genomic DNA was enzymatically digested into 200 bp fragments (5× WGS Fragmentation Mix, Qiagen, USA). T-adapters were added to both ends after repairing and A tailing. For the WES library construction, purified DNA was amplified by ligation-mediated PCR. Then, final sequencing libraries were generated using the 96 rxn xGen Exome Research Panel v1.0 (Integrated DNA Technologies, USA), according to the manufacturer’s instructions. Paired-end multiplex samples were sequenced with the Illumina NovaSeq 6000 System (Illumina, USA). The sequencing depth was 200× per tissue sample and 100× per white blood cell (WBC) sample.

Raw sequencing data were preprocessed by Fastp to trim adaptor sequences (12). Then, clean reads in FastQ format were aligned to the reference human genome (hg19/GRCh37) by Burrows-Wheeler Aligner (BWA, v0.7.15) (13). SAM tools (14) and Picard (2.12.1) (http://picard.sourceforge.net/) were used to sort mapped BAM files and process PCR duplicates. To compute the sequencing coverage and depth, the final BAM files were generated by GATK (Genome Analysis Toolkit 3.8) for local realignment and base quality recalibration (15).

Somatic mutational signatures were de-novo analyzed from the clean WES data by the “Somatic Signatures” R package (v2.20.0) (16), according to the non-negative matrix factorization (NMF) method. Four highly confident mutational signatures were derived from the SCH cohort. Then, they were compared with the consensus signatures in the COSMIC dataset (https://cancer.sanger.ac.uk/cosmic/), based on the cosine similarity analysis to nominate each derived signature with the highest COSMIC dataset similarity [i.e., SBS4 (S4), SBS2 (S2), SBS6 (S6), and SBS5 (S5), respectively] for the SCH cohort.

To further determine the distribution of mutational signatures and the frequencies of each patient, deconstructSigs (v1.9.0) was used as previously described (17). Patients harboring the S2, S4, S5, and S6 mutations, as well as S2, S4, S5, and S6 weights, were compared using the Wilcoxon test among the two groups.

Somatic single nucleotide variations (SNVs) were identified from the clean sequencing data by MuTect (18), and somatic small insertions and deletions (InDels) were detected by the GATK Somatic Indel Detector. The ANNOVAR software was used for the annotation of variants based on multiple databases (19), including variant (HGVS), population frequency (1000 Genomes Project, dbSNP, ExAC), variant functional prediction (PolyPhen-2 and SIFT), and phenotype or disease (OMIM, COSMIC, ClinVar) databases. After the annotation, the retained non-synonymous SNVs were screened from disease databases for further analysis with variant allele frequency (VAF) (cutoff ≥ 3%) or VAF for cancer hotspots (cutoff ≥ 1%). Tumor mutation burden (TMB) was calculated with the total numbers of non-synonymous SNVs and indel variants per megabase of coding regions. Dominant tumor neoantigens were predicted using OptiType to infer the individual HLA type (20). Tumor neoantigen burden (TNB) was calculated with the total numbers of neoantigens per megabase of coding regions. Significant driver genes were identified by combining MutsigCV and dNdScv, as previously described (21, 22), with a false discovery rate (FDR) cutoff <5%. Genes with significantly different mutation frequencies among the two groups were determined based on the gene mutation rates in each group using a two-sided Fisher’s exact test with a p-value of 0.05.

Copy number variations (CNVs) for all patients in the SCH cohort were first identified using the Genome Identification of Significant Targets in Cancer (GISTIC) 2.0 algorithm (23). At the chromosomal arm level, significant amplifications or deletions were screened with FDR (cutoff < 10%) for further analyses. At a focal CNV level, significant amplification was screened with FDR (cutoff < 5%) and G-score (cutoff > 0.3). Significant deletion was screened with FDR (cutoff < 5%) and G-score (cutoff < −0.2) for further analyses.

Focal CNV-related gene analysis was performed for each patient based on paired tumor-normal WES data using the GATK depth of coverage with parameters (–minBaseQuality 0 –minMappingQuality 20 –start 1 –stop 500 –nBins 200 –includeRefNSites –countType COUNT_FRAGMENTS). The amplified genes were defined by a copy number ratio of tumor vs. normal >4, while deleted genes were defined by a copy number ratio of tumor vs. normal <0.5. Then, focal CNV-related genes were filtered according to the COSMIC Cancer Gene Census database (https://cancer.sanger.ac.uk/cosmic/) to obtain a cancer-related focal CNV gene list. Genes with significantly different CNV frequencies among the two groups were determined based on the gene alteration rates in each group using a two-sided Fisher’s exact test with a p-value of 0.05.

Somatic mutation and focal CNV-related genes with enriched biological functions and involved pathways were analyzed using the online tool Metascape (https://metascape.org/gp/index.html#/main/step1), based on the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases (https://www.kegg.jp/kegg/kegg1.html).

To investigate intratumor heterogeneity (ITH), mutant allele tumor heterogeneity (MATH) values for each tumor sample were calculated from the median absolute deviation (MAD) and the median of its mutant–allele fractions at tumor-specific mutated loci: MATH = 100 × MAD/median. These analyses were performed in R with default parameters as previously reported (24). Cancer cell fraction (CCF) and clonal and subclonal mutations in each tumor specimen were calculated based on the proportion of mutated reads (VAF) as previously reported (25).

Regarding genome instability analyses, cellular purity, ploidy, and the segmented allele-specific copy number profiles of each specimen’s tumor cells were estimated using Sequenza (26). The fraction of genome altered (FGA) was defined as the percentage of a tumor genome harboring copy number variations against the whole genome. Loss-of-heterozygosity (LOH) segments or mutations were defined by the minor allele copy number or mutation ratio <0.25 (27). Whole-genome doubling (WGD) events were defined as the major allele ploidy >1.5 on at least 70% of at least 11 autosomes as the duplicated autosome number per sample (28). The Wilcoxon rank-sum test was used to compare the median values of the variables between the two groups. A p-value of 0.05 was considered significant.

Immunohistochemical staining was conducted using the Enhance Labelled Polymer System (ELPS). First, the tumor specimen sections were incubated with anti-PD-L1 (CST, 13684, 1:100) and anti-CD8⁺ (CST, 85336, 1:100) at 4°C overnight. Then, they were washed three times with PBS (5 min per wash). Next, the slides were incubated with the corresponding secondary antibodies at 37°C for 30 min, and they were washed three times with PBS (5 min per wash). Furthermore, the slides reacted with 3,3-diaminobenzidine (DAB) and then washed with distilled water. Next, dehydration was conducted, followed by clearing and mounting with neutral gums. Finally, the stained tissue images were captured by the Digital Pathology Slide Scanner (KF-PRO-120, KF-BIO).

To evaluate programmed cell death ligand 1 (PD-L1) expression, a tumor proportion score (TPS) was defined as the number of PD-L1-staining tumor cells divided by the total number of viable tumor cells multiplied by 100. A combined positive score (CPS) was defined as the number of PD-L1-staining cells divided by the total number of viable tumor cells multiplied by 100 (29). Tonsil PD-L1 staining was adopted to ensure the eligibility of the enrolled specimens. Qualified staining was defined as strong positivity for PD-L1 in the intratonsillar cleft epithelium, whereas negative staining was for PD-L1 in lymphocytes (mantle zone and germinal center B cells) and superficial epithelial cells.

We also evaluated whether CD8⁺ T cells were uniformly distributed in the tumor stroma at lower magnification. If CD8⁺ T cells were equally distributed, they were measured in three randomly chosen areas (0.1 mm²) at a 200-fold magnification. If unequally distributed, the corresponding areas were selected at a 200-fold magnification according to CD8⁺ T-cell percentages in areas of different densities (0.1 mm²), as referred from the PD-L1 expression evaluation criteria. The calculation was defined as follows: CD8⁺ T-cell count/0.1 mm² × 10 or CD8⁺ T-cell count/mm².

Raw sequencing data were deposited in the Genome Sequencing Archive (GSA) of the China National Center for Bioinformation (CNCB) (https://ngdc.cncb.ac.cn/gsa/) under accession number subHRA001430.

The R Foundation for Statistics Computing Package (R package, version 3.3.3) was used to perform the statistical analyses. The Fisher’s exact test (for categorical variables) and the Wilcoxon rank-sum test (for continuous variables) were used to analyze the relationship between the two groups. The Kaplan–Meier method was used to estimate the effects on OS and PFS time based on the log-rank tests. A p-value <0.05 was defined as statistically significant. Hazard ratios of multiple factors on OS and PFS time were obtained from the Cox proportional hazards model.

A further expedition of SCLC genomic landscape features was presented through the WES of the large SCH SCLC cohort. First, the mutation spectrum analysis showed that the two most frequent nucleic acid base substitutions of somatic mutations were transversions (C>A/G>T) followed by transitions (C>T/G>A) ( Supplementary Figure 2 ), consistent with previous SCLC studies. Then, mutational signatures were de-novo calculated and characterized from all 177 specimens based on the 96 possible mutation types, according to a previously published method.

Two signatures showed differences between the two groups, compared to the COSMIC mutational signature database: smoking-related S4 and unknown S5 ( Figure 1 ). Generally, SCLC is associated with heavy tobacco smoking. Chemosensitive patients tended to harbor smoking-related S4 mutations more than chemorefractory patients. However, smoking history did not significantly differ between the two groups. The situation had previously been reported that C>A/G>T transversions, typically prevalent in S4, have no significant correlation with the SCLC smoking status (30, 31). In turn, S5 mutations were more frequently found in the chemorefractory group. This indicated that therapeutic vulnerabilities for different SCLC subtypes may be related to molecular change during SCLC tumor development.

SCLC is considered to be a disease of genomic alterations and relatively higher mutation rate per megabase (Mb) (1, 2). We used MutSigCV and dNdScv (FDR q < 0.1) to identify the differential somatic mutation genes between the two groups, and the genes with a mutation frequency greater than 5% were indicated in the graph ( Figure 2A ). Our results showed that KMT2D (30%), LRP2 (30%), OR1N2 (17.5%), KIAA1109 (15.83%), LAMA4 (15%), ZNF469 (14.17%), GPR158 (11.67%), NRK (10.83%), APBA2 (10.83%), RNF213 (10.83%), ABCC1 (9.17%), GLI2 (9.17%), RP1 (9.17%), ADAMTS13 (8.33%), and IQSEC2 (8.33%) were more frequently predicted in the chemosensitive group, while FNDC1 (22.81%), FAT2 (21.05%), SPATA31E1 (17.54%), AOC1 (17.54%), SYNE2 (17.54%), THSD7A (17.54%), TRIM58 (17.54%), OGDHL (15.79%), NTRK3 (14.04%), OR4C6 (14.04%), PTPN13 (14.04%), COL9A1 (12.28%), MYO18A (12.28%), and KDM3B (10.53%) more commonly appeared in the chemorefractory group. This suggested the relativity between the high-frequency mutations of somatic mutation genes and disease recurrence.

Further functional analysis showed that the differential genes more frequently predicted in the chemorefractory group were most significantly enriched in tight junction, lysosome, Hippo signaling pathway, and olfactory transduction. Moreover, the differential genes more frequently predicted in the chemosensitive group were most significantly enriched in neuron projection morphogenesis and localization within the membrane ( Figure 2B ). The results suggested the potential mechanisms of chemoresistance on somatic mutation levels, which warrants further study.

In addition, among the differential genes found between the chemorefractory group and the chemosensitive group, we applied survival correlation analysis and identified eight somatic mutation genes mentioned earlier that significantly reduced PFS time, compared to wild-type genotypes ( Figure 2C ; Supplementary Figures 3 , 4 ). Meanwhile, three somatic mutation genes that were more frequently predicted in the chemosensitive group significantly increased PFS time. This further strengthens the ( Figure 2D ) fact that they have pivotal roles in SCLC relapse and chemotherapy resistance. There was no significant difference, however, in OS time apart from LRP2 ( Supplementary Figure 5 ).

The VAF analysis reflected that VAF in the Re group was higher than that in the Se group in all gene mutations (p = 6.2e−11) and clonal gene mutations (p = 0.022), but there was no difference concerning driver mutations (p = 0.41) or LOH (p = 0.27). No difference was observed between the two groups with respect to TMB (p = 0.12) and MATH score (p = 0.28) ( Supplementary Figure 6 ).

CNVs, novel structural variations in human chromosomes, are extremely common in SCLC, and our results had come to similar conclusions ( Figure 3 ). We identified a high frequency of significant mutations in the two groups: C8orf82CEP72 (42.11%), EXOC3 (28.07%), PLEKHG4B (24.56%), CEP72 (22.81%), ING1 (21.05%), SYNGR3 (19.3%), HAGH (19.3%), MAP7D1 (17.54%), VPS28 (17.54%), STX10 (17.54%), HSD11B1L (17.54%), FAM195A (15.79%), BTNL3 (15.79%), RPL8 (15.79%), ZNF414 (15.79%), TNFRSF12A (15.79%), LRRC24 (15.79%), KIFC2 (15.79%), FBXL16 (15.79%), and GLI4 (14.04%) in the chemorefractory group and CTRB2 (22.5%), LOC101928018 (18.33%), PQLC1 (13.33%), CLEC18C (11.67%), CTIF (10.83%), and SIVA1 (10.83%) in the chemosensitive group. The above somatic CNV analyses showed that the chemorefractory group seemingly experienced more mutations of the CNV site compared with the chemosensitive group ( Figure 3A ). We have noted that these changes involved certain oncogenic signaling pathways ( Figure 3B ) whose alterations significantly impacted patients’ PFS time ( Figure 3C ; Supplementary Figures 7 - 11 ). Similar to previous results, these changes had very little impact on OS time ( Supplementary Figures 7 - 12 ).

The subsequent analysis revealed differences between the two groups. The segmented copy numbers were visualized in a heatmap and the significance of chromosome alterations was determined by GISTIC analysis ( Figures 4A, B ). The differing clones of mutations were prevalent in the two types of specimens with a few similar parts. We found more amplified mutations in the chemorefractory group, which are related to certain functions, such as transcription factor binding, regulation of hemopoiesis, leukocyte differentiation, peptidyl-tyrosine phosphorylation, positive regulation of cell death, transcription regulator complex, regulation of cellular response to stress, chromatin binding, histone modification, regulation of kinase activity, damaged DNA binding, response to radiation, homeostasis of a number of cells, positive regulation of endothelial cell proliferation, ubiquitin protein ligase binding, rhythmic process, regulation of cellular response to growth factor stimulus, response to hypoxia, ankyrin binding, and negative regulation of catabolic process ( Figures 4C, D ). Amplified mutations in the chemosensitive group inhibited the regulation of the cell cycle process, transcription factor binding, and B-cell differentiation. Analysis of deletion mutation also suggested that the detection of the chromosomal-level variation might influence our treatment strategies ( Supplementary Figure 13 ).

Immune checkpoint inhibitors (ICIs) have altered the treatment of SCLC (32–34). Multiple biomarkers, such as PD-L1 and CD8⁺ tumor-infiltrating lymphocytes (TILs), have been identified to help tailor immunotherapy. According to our current study, higher PD-L1 expressions (based on TPS methods) were mostly present among chemosensitive patients (p = 0.026; Figure 5 ), while there were no differences between the two groups in terms of PD-L1 expressions (based on CPS methods) and CD8⁺ TILs. Significantly, hyperprogressive disease (HPD), a pattern of progression in which a flare-up of tumor growth occurred during immunotherapy, was similar between the two groups. This suggests that chemosensitive patients might more likely benefit from immunotherapy.

It has been shown that some crucial genetic mutations could influence the efficiency of ICI treatments (35). However, our study suggests that there were no differences in genomic instability between the two groups. Analyses that take more factors into account are needed in our future studies ( Supplementary Figure 14 ).

There was a significant difference in the survival rates of different resistance levels in SCLC patients ( Figure 6A ). In this study, we attempted to screen optimal subsets of features between the two groups and set up a predictive model of drug resistance of SCLC. The statistical analysis of data from three different angles offered us some clues ( Table 1 ; Supplementary Table 1 ). The first part included the clinical characteristics, age, gender, stage, family history, smoking, drinking, metastasis, PD-L1 expressions, and CD8⁺ TILs. The factors with p-value ≤0.2 were selected by univariate logistic regression analysis, and age, stage, family history, and CD8⁺ TILs were chosen. Similarly, differences at the molecular level, such as ABCC1, APBA2, GPR158, KMT2D, NTRK3, TRIM58, FNDC1, FAT2, OR1N2, LRP2, and KIAA1109, were obtained. The same goes for features of chromosome variation that included C8orf82, CTRB2, EXOC3, PQLC1, and BTNL3.

There were 20 eigenvalues in total as discussed previously. Step function was used to determine multivariate logistic regression through a stepwise regression process with resistance as the target variable. Eventually, 16 eigenvalues were selected ( Figures 6B, C ).

According to the model determined by logistic regression, each sample was endowed with a predictive probability value (PV). Then, the samples were divided into a high-risk group (>0.55) and a low-risk group (≤0.55), and the survival analysis showed obvious differences between the groups. This suggests that the predictive model, to a certain extent, could predict the drug resistance of SCLC.