All study plans and experimental protocols were submitted to the ethics/licensing committee of Shanghai Pulmonary Hospital, and all of them had been approved. Written informed consents were obtained from all patients involved. All methods, personal training, and experiments were performed following relevant regulations and guidelines.

This study aimed to enroll lung cancer patients without EGFR [single-nucleotide variant (SNV), insertion/deletion (INDEL)] or ALK (fusion) alterations and analyze their DDR-related genomic and transcriptomic characteristics. All patients enrolled received first-line chemotherapy. All eligible patients received polymerase chain reaction (PCR) assays for EGFR and ALK before any therapy. After that, samples from 122 qualifying advanced lung cancer patients were eligible, including 86 NSCLC and 36 SCLC participants (Table S1). Tumor specimens and blood samples were collected as we previously described (21).

Sample storage, DNA extraction, and DNA sequencing were done as previously described (21).

Somatic variation was detected by Genome Analysis Toolkit (GATK, Version 4.1.7.0). We followed the workflows of Burrows–Wheeler aligner (BWA, Version 0.7.17-r1198) for aligning sequencing data to the hg19 genome (GRch37). Duplicated reads underwent subsequent marking and were removed by the GATK Picard tool. Base quality score was recalibrated via BaseRecalibrator and ApplyBQSR functions of GATK, and Mutect2 from GATK was designed for calling SNVs and INDELs from tumor-normal matched pairs. Above analyses were performed on the cloud-based genomic analysis platform: Biomedical Data Analysis Platform (BMAP, https://bmap.sjtu.edu.cn/). Significantly mutated genes mutated more frequently than expected accidentally were determined through MutSigCV (Version: 1.41) with q values <0.1 (22) in our study.

Tissue DNA or RNA was extracted based on the manufacturer’s protocol, and reversed transcript would be performed for extracted RNA for subsequent PCR amplification. EGFR mutations and EML4-ALK fusion were detected as described (21).

We evaluated seven major DDR pathways, BER, MMR, HR, NER, NHEJ, FA, and cell cycle checkpoint in our study. A total of 74 DDR genes were assembled as being associated with DDR, grouped into different functional pathways from published resources (15, 19, 23–27) (Table S2). We considered all loss-of-function alterations deleterious, including nonsense mutations, splice site, or frameshift alterations (Table S3) (19). Two diverse methods were applied to determine the functional impacts of missense mutations: (1) by in silico functional analysis, all missense mutations that were classified as “probably damaging” or “possibly damaging” in Polyphen2 (28) or “high” or “medium” in MutationAssessor (29) were recognized deleterious; (2) we manually reviewed the missense mutations in Catalogue of Somatic Mutations in Cancer (COSMIC) (30), algorithmically identified recurrent hotspot mutations (31), and annotated oncogenicity via OncoKB. DDR gene alterations were defined as deleterious DDR mutations. In our study, DDR-deficient subtype was defined as individuals with deleterious DDR mutations. DDR-proficient subtype was defined as individuals without deleterious DDR mutations. The deficient and proficient individuals in a particular pathway were identified as individuals with or without mutation in this particular pathway.

We used SignatureAnalyzer to infer mutational signatures of our samples (32, 33) (http://software.broadinstitute.org/cancer/cga/msp). SignatureAnalyzer applied a Bayesian variant of non-negative matrix factorization algorithm for signature analysis, and mutational signatures were identified by comparing with 30 COSMIC mutational signatures.

Consensus clustering was performed using the unsupervised tool named ConsensusClusterPlus (34). We set Pearson correlation distances for distance, 80% item resampling for pItem, and 10 resamplings for reps. Eventually, we determined three clusters among these patients.

We calculated TMB by dividing the total number of cancer tissue non-synonymous variations (SNV and INDEL, allele frequency >5%) by the length of the whole-exome sequencing (WES) panel. In addition, we used the MATH score as a quantitative measure for ITH, which considered the width of variant allele frequency distribution for calculation (35).

To figure out relevant expression levels of immune-related genes in lung cancer patients without EGFR/ALK mutations between DDR-deficient and DDR-proficient groups, appropriate data were obtained from databases like The Cancer Genome Atlas (TCGA) obtained from cBioportal, lung adenocarcinoma (LADC) cohort from dataset “Lung Adenocarcinoma (TCGA, Firehose Legacy)” (https://www.cbioportal.org/study/summary?id=luad_tcga), lung squamous cell carcinoma (LUSC) cohort from dataset “Lung Squamous Cell Carcinoma (TCGA, Firehose Legacy)” (https://www.cbioportal.org/study/summary?id=lusc_tcga), and small cell lung cancer (SCLC) cohort from dataset “Small Cell Lung Cancer (U Cologne, Nature 2015)” (https://www.cbioportal.org/study/summary?id=sclc_ucologne_2015). After excluding patients with EGFR/ALK genomic alterations, a total of 64, 72, and 116 patients were included in the LADC, LUSC, and SCLC cohorts, respectively.

We performed univariate and multivariate Cox regression analyses considering clinical information (including age, gender, smoking history, clinical stage, and pathological type), and TMB, ITH, and DDR gene mutations as variables. We divided patients into training and testing groups randomly (7:3). Factors selected from Cox regression were included for building a nomogram model in the training group (36). Prognostic values at 1- and 2-year survival of lung cancer patients were predicted. Then, time-dependent receiver operating characteristic (ROC) curves were performed to evaluate the performance of this nomogram in both training and testing datasets. After that, Kaplan–Meier survival analysis was further used for evaluating the clinical value of this predictive model. The cutoff value of the total point of nomogram was determined using R package “survminer.”

The findings of our data were validated using external datasets. We used datasets from cBioportal to form the external validation cohort for verifying the performance of the clinical prediction model based on the DDR mutational status. The validation cohorts we used for exploring the relationships of DDR gene alterations with TMB, ITH, and survival, and verifying model performance were obtained from cBioportal (LADC, LUSC, and SCLC cohorts mentioned above). Patients with EGFR/ALK genomic alterations were excluded.

We determined objective response rate (ORR), disease control rate (DCR), PFS, and OS based on Response Evaluation Criteria in Solid Tumors version 1.1 (37). PFS was defined as the interval from the start of first-line therapy to the date of disease progression or death. OS was defined as the interval from the start of first-line therapy to death. Wilcoxon rank-sum test was applied for comparing mutational burden between defined subgroups, and two-tailed Student’s test was performed for comparing immune-related gene expression and immune infiltration between DDR-deficient and DDR-proficient subtypes. Associations between DDR-deficient and DDR-proficient groups were analyzed by chi² or Fisher’s exact tests for qualitative data. We used maftools (38) for the detection of co-occurring or mutually exclusive sets of genes, in which pairwise Fisher’s exact test was used to explore significant gene pairs. Multiple logistic regression analysis was performed by categorizing TMB based on the median to analyze the effects of DDR mutations and smoking on TMB. The Kaplan–Meier curves with log-rank test were performed for testing survival differences between two subgroups. Cox regression was also used for the determination of clinical values of single DDR genes. We also used Benjamini–Hochberg test for p-value correction to test multiple hypotheses when appropriate. Extern mutational and clinical data were obtained from TCGA via cBioportal. All visualizations were achieved using R software. We defined “*” as statistically significant (p < 0.05), “**” as highly statistically significant (p < 0.01), and “***” as very highly statistically significant (p < 0.001). We defined p < 0.05 as statistically significant.

We identified the clinicopathological features and DDR mutation spectrum of 122 advanced lung cancer samples without EGFR or ALK mutations/translocations. In brief, 86.9% (106/122) of patients were male, and 69.7% (85/122) of patients had smoking history (Figure 1A). In our study, the proportion of pathologically determined LADC and SCLC were 34.4% (42/122) and 29.5% (36/122), respectively. LUSC accounted for 30.3% (37/122) (Figure 1A). A total of 67 patients were identified with deleterious alterations of DDR genes (67/122, 54.9%, Figure 1A; Table S3). Their clinical characteristics are summarized in Table S1. POLQ, BRCA2, ATM, ATR, PARP4, and POLD1 alterations were most commonly observed in the entire advanced lung cancer cohort (Figure 1A). Different histopathological types exhibited specific mutation characteristics in the DDR pathways (Figure 1B). For example, MMR alterations were relatively more common in LADC (8/42, 19.0%) but less frequent in SCLC (1/36, 2.8%; LADC vs. SCLC, p = 0.059). FA alterations were observed more frequently in LUSC (9/37, 24.3%), while 9.5% were observed in the LADC subtype (4/42, 9.5%, p = 0.077). However, the proportion of patients with DDR gene alterations revealed no difference either between NSCLC and SCLC or between LADC and LUSC (Figure 1C).

Then, we explored the relationship between distinct clinical characteristics and DDR mutation. No statistical difference was observed in smoking (p = 0.899) and gender (p = 0.335) between DDR-deficient and DDR-proficient subtypes by chi² or Fisher’s exact test. The same result was also observed in age by Student’s t-test (p = 0.832). For SCLC, a statistically significant result was observed between smoking and DDR mutation (p = 0.045). No significant relationship was detected between these clinical factors and DDR mutation in LADC and LUSC (all p > 0.05).

By detecting somatic interactions of genes of our data, we found that there were widespread comutations between DDR genes and significantly mutated genes. For example, significantly mutated gene KRAS alteration was observed co-occurring with alterations of DDR genes MDC1 (p = 0.026); DCAF8L2 alteration was observed co-occurring with alterations of ATM (p = 0.010), BARD1 (p = 0.022), FANCM (p = 0.022), and RAD50 (p = 0.042, Figure 1D). No mutually exclusive sets between DDR genes and significantly mutated genes were found. Unlike significantly mutated genes showing both comutations and mutual exclusions, between DDR genes, extensive comutations were extremely frequent. DDR gene coalterations occur both in the same and different DDR pathways. For instance, HR gene RAD51B was observed comutated with BAP1 (p = 0.017), BARD1 (p = 0.025), RAD50 (p = 0.033) in the same pathway, and ERCC2 (p = 0.017), ERCC4 (p = 0.017), FANCI (p = 0.017), FANCM (p = 0.025), MLH3 (p = 0.025), PARP2 (p = 0.008), and RARP3 (p = 0.017) in other DDR pathways.

In our study, the median TMB for all patients enrolled was 6.06 mutations/megabase, ranged from 0.17 to 67.95. When comparing the TMB status of groups with DDR gene alterations to those harboring no DDR-related genomic mutations, we identified that TMB was comparably higher in participants with DDR genomic alteration than DDR-proficient patients (p = 0.007 in LADC cohort and p = 0.003 in LUSC cohort; Figure 2). SCLC also showed a similar trend, although no significant difference was observed (p = 0.389). Different histological subtypes revealed substantial differences in TMB distribution in different DDR mutation states. For LADC, groups with mutations of MMR (p = 0.039), HR (p = 0.006), and cell cycle checkpoint genes (p = 0.017) showed significantly higher TMB than corresponding proficient groups, while NER (p = 0.922), FA (p = 0.075), BER (p = 0.188), and NHEJ (p = 0.171) failed to reveal the difference (Figure 2A). LUSC patients with DDR deficiencies of MMR (p = 0.012), NER (p = 0.002), HR (p = 0.044), BER (p = 0.009), and NHEJ pathways (p = 0.024) showed significantly higher TMB than corresponding proficient participants (Figure 2B). However, for SCLC, only HR (p = 0.044) and NHEJ (p = 0.013) alterations showed significant differences; MMR (p = 0.222), NER (p = 0.622), FA (p = 0.467), BER (p = 0.396), and cell cycle checkpoint alterations (p = 0.622) showed no statistical significance (Figure 2C). Then, we classified high- and low-TMB using a cutoff of median TMB value (Table S4). By binary comparison, positive associations between TMB-high group and genomic alterations were observed in the mutational status of DDR (p = 0.006), HR (p = 0.014), and cell cycle checkpoints (p = 0.013).

We also identified three mutational signatures (signatures 2, 4, and 6) compared with 30 COSMIC mutational signatures and divided the whole cohort into three subgroups (C1–C3) via unsupervised clustering analysis latter. Figure 2D shows that C2 subgroup had enriched signature 6, which might have a relationship with deficient MMR (dMMR). High proportions of DDR gene mutations of samples were detected in groups C1 and C2 (Figure 2E). By comparing the proportions of samples carrying various DDR gene mutations in each signature/cluster group, we found that C1 group had the highest proportion of MMR alterations (13.59%) in comparison with C2 and C3 groups. Samples most relevant to signature 6 (dMMR) failed to show relatively high TMB (Figure 2F). These suggested, although the signature 6 mutation feature was detected in our data, the effects of signature 6 and mutational status of DDR genes on TMB were different. The mutational status of DDR genes reflected the mutation load of lung cancer patients better. We then compared the effects of DDR mutations and smoking on TMB (Figure 2G). By Wilcoxon rank-sum test, DDR-deficient subtype was associated with significantly higher TMB among both smokers and non-smokers. No significant correlation was observed between smoking status and TMB among both DDR-deficient and DDR-proficient cohorts. We also performed multiple logistic regression analysis by categorizing TMB based on the median to analyze the effects of DDR mutation and smoking on TMB. It showed that DDR mutations and smoking were both significantly associated with increased TMB (p = 0.006 and p = 0.030, respectively). The regression coefficient of DDR mutation (1.05) was much larger than smoking (0.92), which suggested that in lung cancer patients without EGFR/ALK alterations, DDR mutation had a greater impact on TMB than smoking.

MATH score was calculated for ITH evaluation (Figure 3A). In our study, Pearson correlation was only 0.015 for TMB and ITH (p = 0.873), which negated the correlation between TMB and ITH (Figure 3B). DDR-deficient patients showed significantly higher ITH than DDR-proficient patients (p = 0.043) in LUSC, while no significance was suggested in LADC or SCLC (Figures 3C–E). For specific DDR pathways, significant difference was only observed in LADC patients with NHEJ alterations (p = 0.048; Figure 3C). In SCLC, DDR gene alterations had no significant effect on ITH (Figure 3E). Among all patients, except for NHEJ alteration (p = 0.036), mutations in the other DDR pathways were not related to MATH score. By dividing all patients into high and low MATH score group using a cutoff of median value, we found no correlation between ITH status and DDR pathway alterations (p > 0.05, Table S5).

Based on TMB and ITH, survival analyses were performed to identify their predictive and prognostic value for first-line chemotherapy in lung cancer patients without EGFR/ALK mutations. High-TMB patients suggested no survival benefits in both PFS and OS compared with TMB-low patients in all (p = 0.66, p = 0.46), NSCLC (p = 0.43, p = 0.46), LADC (p = 0.51, p = 0.92), LUSC (p = 0.26, p = 0.43), and SCLC cohorts (p = 0.2, p = 0.78; Figures S1A, S1B). We also observed no significant prognosis value of ITH in either NSCLC or SCLC (Figures S1C, S1D).

The relationship between DDR pathways and therapeutic response to initial chemotherapy in distinct populations was identified. Table S6 suggested no significant correlations among DDR gene alterations with regards to ORR and DCR.

To assess the predictive role of DDR in first-line chemotherapy, we identified PFS and OS by analyzing different gene sets and individual DDR genes. We first compared patients’ survival between DDR-deficient and DDR-proficient patients; however, no clinical association was found when analyzing PFS and OS in distinct cohorts (p > 0.05; Figure S2). We then evaluated the clinical efficacy based on levels of diverse DDR pathways. Among all patients, NER aberrations showed poor PFS (vs. NER-proficient patients, 151.000 days, 95% CI 13.569–288.431 vs. 394.000 days, 95% CI 275.953–512.047; p < 0.001) and OS (vs. NER-proficient patients, 241.000 days, 95% CI 235.868–246.132 vs. 438.000 days, 95% CI 257.813–618.187; p = 0.013). BER alterations (vs. BER-proficient patients, 231.000 days, 95% CI 118.041–343.959 vs. 397.000 days, 95% CI 163.308–630.692; p = 0.033) indicated significantly poor PFS (Figures 4A–C). BER mutation also showed pooper PFS in the LADC cohort (vs. BER-proficient patients, p = 0.026; Figure 4D). For LUSC, groups with NER aberrations showed shorter PFS (p < 0.001) and OS (p = 0.02) than NER-proficient patients (Figures 4E, F). For SCLC, we observed comparably shorter PFS in patients with alterations of HR pathways (vs. HR-proficient group, 197.000 days, 95% CI 164.656–229.344 vs. 411.000 days, 95% CI 286.087–535.913; p = 0.04; Figure 4G). Although significant survival differences were observed in MMR-deficient status (p = 0.028) for PFS and cell cycle checkpoint-deficient status (p < 0.001) for OS in SCLC, and NHEJ-deficient status for PFS in LUSC (p < 0.001), they failed to reveal the actual situation for limited mutated samples (Tables S7, S8). No significant survival differences were observed in other pathways (Tables S7, S8). At the same time, in the whole lung cancer patients, by using univariate cox analysis, somatic alterations in single genes, including MLH3 (HR = 3.311; 95% CI, 1.022–10.727; p = 0.046), ERCC2 (HR = 15.183; 95% CI, 1.828–126.123; p = 0.012), ERCC4 (HR = 15.183; 95% CI, 1.828–126.123; p = 0.012), PARP2 (HR = 15.183; 95% CI, 1.828–126.123; p = 0.012), BAP1 (HR = 5.234; 95% CI, 1.235–22.173; p = 0.025), RAD51B (HR = 15.183; 95% CI, 1.828–126.123; p = 0.012), FANCB (HR = 5.125; 95% CI, 1.197–21.951; p = 0.028), and FANCI (HR = 7.450; 95% CI, 1.711–32.448; p = 0.007; Figure 4H) showed statistically shorter OS after current therapy. The impacts of DDR-related single genes on patients’ OS were also identified in NSCLC and SCLC (Figures 4I, J).

Considering the great significance of DDR pathways and single genes in predicting the prognosis of advanced lung cancer patients, we then constructed a DDR-based predictive model to predict patients’ prognosis. Apart from the univariate analysis of DDR-related alteration (Figure 4H), we also conducted univariate analysis based on basic clinical parameters, TMB, and MATH score (ITH), which revealed that the clinical characteristics including clinical stage (p < 0.001), smoking status (p = 0.017) were significant factors for patients’ prognosis. No significant results were obtained in different groups divided by age (p = 0.257), gender (p = 0.122), histological subtype (p = 0.532), TMB (p=0.337), and ITH (p = 0.695). Factors with significant difference were included for multivariate Cox regression analysis (Table S9). The Cox regression selected three independently predictive factors. Based on the result from multivariate Cox regression, a prediction model was established using the nomogram algorithm in the training dataset, which included both clinical factors and DDR-related alterations (Figure 5A). Smoking status was also included for model construction considering its significance to patients' prognosis. This combined prediction model showed excellent performance, with area under the curve (AUC) values of 0.732 and 0.789 for 1- and 2-year receiver operating characteristic (ROC) curves in the training group, and 0.724 and 0.726 for 1- and 2-year ROC curves in the testing group (Figures 5B, C). In addition, we also evaluated its clinical value and found that patients with high total points showed worse clinical outcomes than those with low total points (third quartile of OS, high risk vs. low risk, 188.000 days vs. 726.000 days; p = 0.003; Figure 5D). We have modified Figure 5, which was uploaded in the proof website.

We verified the correlations of TMB and ITH with DDR gene alterations in patients with no EGFR or ALK alterations from cBioportal database at first. Totally, consistent with our finding, higher TMB was found in DDR-deficient groups, and HR alterations were related to higher TMB in LADC, LUSC, and SCLC (Figures S3A–C). ITH showed significant associations with DDR altered status in LUSC and SCLC but not in LADC (Figures S3D–F). After analyzing data from public databases, high-TMB also suggested no clinical benefit in OS in LUSC (p = 1) and SCLC (p = 0.84), while for LADC, the TMB-high group demonstrated poorer OS (p = 0.002; Figures S3G–I). For ITH, only SCLC with high MATH scores showed longer survival (p = 0.014; Figures S3J–L). Then, we utilized univariate and multivariate analyses considering clinical factors and DDR-related pathway and single gene alterations and constructed a nomogram model as well in public lung cancer cohorts. Interestingly, in both LADC and SCLC public cohorts, the predictive model based on DDR gene alterations showed good performance (Figure S4), which were also better performed than the model without factors of DDR gene alterations.

To elucidate underlying immune mechanisms, we first compared somatic interactions of DDR genes and immunomodulatory genes. Figure 6A demonstrated the extensive comutations between DDR and immune genes. Mutations in immune-related genes linked to specific DDR gene alterations significantly. Then, we analyzed the potential immune evasion caused by class I HLA genes (39). DDR gene alterations were observed comutated with various genes of HLA class I pathway including TAP1 and TAP2 when compared with DDR-proficient samples, although no significant difference of mutation frequency was found between the two groups (Table S10). We also analyzed type I interferon (IFN) genes that attach great importance to optimal immunosurveilance and antitumor efficacy. IFNA8 and IFNA10 mutations were observed in DDR-deficient cohort, while no relevant genes were observed in DDR-proficient cohort. However, type I IFN gene mutations were rare in our DDR-deficient cohort, with only 1 of 67 samples containing alterations in IFNA7, IFNA8, IFNA10, IFNA13, and IFNB1 (Table S10). Besides, we found no difference in the mutation frequency of immune-stimulated or inhibited genes between the DDR-deficient and DDR-proficient samples (p > 0.05; Table S10). Then, we compared the expression characteristics of immune-related genes between DDR-deficient and DDR-proficient samples without EGFR or ALK mutations from the TCGA database (Figure 6B). Table S11 summarized genes with p < 0.05. In LADC, DDR-deficient types had relatively lower expression of immune-stimulated genes (such as STING1, CD28, HLA-DRB6) compared with DDR-proficient types. Compared with DDR proficient-type, TNFRSF4 (p = 0.017) and TGFB1 expressions (p = 0.033) were lower in DDR-deficient LUSC; CD40 (p = 0.012) and TNFRSF14 expressions (p = 0.022) were relatively higher in DDR-deficient SCLC. Furthermore, we examined the immune microenvironment characteristics of lung cancer patients without EGFR mutations and ALK fusion (Figure 6C). The heatmaps suggested that immune infiltration varies with histological subtypes. Among differential immune cells (Table S12), DDR deficiency slightly decreased the infiltration of immune cells including resting mast cells (p = 0.003), memory B cells (p = 0.025), resting dendritic cells (p = 0.035) in LADC, and macrophage M0 (p = 0.049) in LUSC with a relatively lower cell-fraction than DDR-proficient type; activated mast cells in LADC (p = 0.044) and macrophage M2 in LUSC (p = 0.004) were slightly higher.