Genetic Determinants of EGFR-Driven Lung Cancer Growth and Therapeutic Response In Vivo

During tumor evolution, cancer cells accumulate alterations in oncogenes and tumor suppressor genes, which contribute to many of the hallmarks of cancer (1). Despite their extensive genomic complexity, tumors are frequently classified based on the presence of oncogenic driver mutations, whereas the function of coincident tumor suppressor gene alterations is largely ignored. There is emerging evidence that the interplay between oncogenic drivers and tumor suppressor genes may influence tumor fitness and alter treatment response (2). However, the combinatorially vast landscape of genomic alterations makes it difficult, except in the most extreme cases, to glean information about the epistatic interactions between tumor suppressor genes and oncogenes from human cancer sequencing data alone (3). This complexity makes inferring the relationship between genotype and therapy response even more tenuous.

Recently, high-throughput, tractable systems that combine autochthonous mouse modeling and genome editing have been developed to directly uncover the functional consequences of genetic alterations during tumorigenesis in vivo (4–10). However, very few studies have investigated the biological consequences of inactivating tumor suppressor genes in the context of different oncogenic drivers in vivo, and existing knowledge about the role of specific genetic alterations in tumor suppressor genes has been primarily inferred from correlative human studies (11–13).

In lung adenocarcinoma, EGFR and KRAS are the two most frequently mutated oncogenic driver genes, and they occur within a background of diverse putative tumor suppressor gene alterations (2, 11–14). Among these, TP53 is the most commonly mutated tumor suppressor gene in both oncogenic EGFR- and KRAS-driven lung adenocarcinoma, consistent with the importance of disrupting this pathway during lung cancer development (2, 13, 15–18). Interestingly, many other putative tumor suppressor genes are mutated at different frequencies in oncogenic EGFR- and KRAS-driven lung adenocarcinomas (2, 11–13, 15). Whether these differences are due to different fitness effects that depend on the oncogenic context has never been tested experimentally. Indeed, previous studies on tumor suppressor genes in lung cancer models in vivo have been primarily conducted in the context of oncogenic KRAS-driven tumors, whereas the functional importance of different tumor suppressor genes in EGFR-driven lung tumors remains largely unstudied (Supplementary Fig. S1A).

In addition to driving growth, inactivation of tumor suppressor genes may affect therapeutic responses. In advanced EGFR-mutant lung adenocarcinomas, treatment with EGFR tyrosine kinase inhibitors (TKI) is the first line of therapy (19–22). Response rates to TKIs are high, but there is large variability in the depth and duration of response among patients, and acquired resistance inevitably occurs (14). Genomic alterations, including those in RB1 and TP53, have been found to correlate with clinical outcomes to TKI treatment (13, 17, 23). However, given the complexity and diversity of genomic alterations in these tumors, the functional contribution of tumor suppressor genotypes to TKI sensitivity and drug resistance remains poorly understood.

To quantify the functional importance of a panel of 10 diverse putative tumor suppressor genes on oncogenic EGFR-driven lung tumor growth in vivo, we coupled multiplexed CRISPR/Cas9-mediated somatic genome editing and high-throughput barcode sequencing (Tuba-seq) with a novel genetically engineered mouse model of EGFR^L858R-driven Trp53-deficient lung cancer. Through the comparison of tumor suppressor effects in oncogenic EGFR- and KRAS-driven lung cancer models, we uncovered prevalent epistasis between the oncogenic drivers and tumor suppressors, which explains the different mutational spectra of these tumor suppressor genes in oncogenic EGFR- and KRAS-driven human lung adenocarcinomas. Moreover, we established a direct causal link between tumor suppressor genotypes and differential responses to osimertinib treatment in EGFR-driven lung adenocarcinoma that is supported by correlative human mutational datasets.

The use of lentiviral vectors to initiate tumors in genetically engineered mouse models of human cancer enables control of tumor multiplicity, tumor barcoding to map clonality, and the delivery of lentivirus-encoded complementary DNAs, short hairpin RNAs (shRNA), and single guide RNAs (sgRNA) to modify neoplastic cells (6, 8, 24, 25). The simplicity of viral Cre-initiated models of oncogenic KRAS-driven lung cancer has enabled the analysis of many genes that cooperate to drive tumor growth within these autochthonous mouse models (Supplementary Fig. S1A; ref. 26). To permit the generation of virally initiated oncogenic EGFR-driven Trp53-deficient lung tumors, we bred mice with a reverse tetracycline-controlled transactivator (rtTA)/tetracycline-inducible transgene encoding the common lung adenocarcinoma–associated EGFR^L858R mutant (TetO-EGFR^L858R), a Cre-regulated rtTA transgene (Rosa26^{CAGs-LSL-rtTA3-IRES-mKate}, abbreviated R26^RIK), and homozygous Trp53 floxed alleles (p53^flox/flox; refs. 27–29). In these TetO-EGFR^L858R;Rosa26^RIK;p53^flox/flox (EGFR;p53) mice, lentiviral Cre transduction of lung epithelial cells leads to the production of rtTA and mKate, as well as the inactivation of Trp53. Coincident doxycycline treatment induces the expression of oncogenic EGFR (Fig. 1A).

We initiated tumors with a lentiviral Cre vector (30) in EGFR;p53 mice and used magnetic resonance imaging (MRI) to monitor tumor development (Fig. 1B). Tumors were first detectable in EGFR;p53 mice 8 weeks after tumor initiation. Histologic analysis of lungs 11 weeks after tumor initiation confirmed the development of multifocal lung adenocarcinomas (Fig. 1B, bottom panels). Tumors stained positively for EGFR^L858R and mKate, as well as surfactant protein C (SFTPC) and the lung lineage–defining transcription factor NKX2-1/TTF1 (Fig. 1C, top panels; Supplementary Fig. S1B; refs. 31–33). Importantly, in this model, tumors were more focal than the diffuse tumors that rapidly develop in the previous CCSP-rtTA;TetO-EGFR^L858R model (27), likely due to tumor initiation from fewer cells in the virally initiated model (Fig. 1B, top right panel). Lentivirus-induced tumors in EGFR;p53 mice were more poorly differentiated than those typically observed in the CCSP-rtTA;TetO-EGFR^L858R model (Fig. 1B, bottom panels), as shown by increased prevalence of a micropapillary pattern that is highly aggressive in human adenocarcinoma (27). Thus, this new lentiviral Cre-initiated model recapitulates the genetic and histopathologic features of human oncogenic EGFR-driven TP53-deficient lung tumors.

To enable somatic genome editing in the EGFR;p53 model, we further incorporated a conditional Cas9 allele (Rosa26^LSL-Cas9-GFP; ref. 4) to generate TetO-EGFR^L858R;R26^{RIK/LSL-Cas9-GFP};p53^flox/flox (EGFR;p53;Cas9) mice (Fig. 1A). Lentiviral Cre delivery to EGFR;p53;Cas9 mice initiated multifocal lung adenomas and adenocarcinomas that expressed EGFR^L858R, mKate, Cas9, and GFP (Fig. 1C and D). Tumors in EGFR;p53;Cas9 mice were histologically similar to those in EGFR;p53 mice (Fig. 1B and D).

We used an improved version of our Tuba-seq approach to quantify tumor suppressor gene function in oncogenic EGFR-driven lung tumors (ref. 34; Supplementary Methods). Genomic integration of barcoded lentiviral vectors uniquely tags each transduced cell and all of the neoplastic cells within the resulting clonal tumors (8). Each barcode encodes an eight-nucleotide sgID region specific to the sgRNA followed by a random 15-nucleotide barcode (sgID-BC); thus, high-throughput sequencing of this sgID-BC region from bulk tumor–bearing lung can be used to quantify the number of cells in each tumor of each genotype (Supplementary Methods; ref. 8). The absolute number of neoplastic cells in each tumor is calculated by normalizing the number of reads of each unique barcode to the number of reads from benchmark control cells added to each sample (Fig. 2A; Supplementary Fig. S2A; Supplementary Methods). Thus, Tuba-seq enables the parallel analysis of the impact of multiple tumor suppressor gene alterations on tumor growth in vivo.

To assess the function of 10 diverse putative tumor suppressor genes that are frequently altered in human lung adenocarcinoma, we initiated tumors in EGFR;p53 and EGFR;p53;Cas9 mice with a pool of barcoded Lenti-sgRNA/Cre vectors (Lenti-sgTS^Pool/Cre; Fig. 2A). In addition to Lenti-sgRNA/Cre vectors targeting each putative tumor suppressor gene, this pool contains negative control vectors, including four Lenti-sgInert/Cre vectors and a vector with an sgRNA targeting Trp53, which is already inactivated by Cre-mediated recombination in EGFR;p53;Cas9 mice (Fig. 1A; ref. 8). Lenti-sgTS^Pool/Cre–initiated tumors were first detectable by MRI 4 weeks after tumor initiation in EGFR;p53;Cas9 mice. At 11 weeks after tumor initiation, when tumors were readily detectable in all mice, tumor-bearing lungs were collected for Tuba-seq analysis and histology (Fig. 2B; Supplementary Fig. S1E).

We used Tuba-seq to quantify the number of neoplastic cells in clonal tumors of each genotype in the EGFR;p53 model and used two summary statistics to describe the tumor size distribution: percentiles within the tumor size distributions and the log-normal (LN) mean (Fig. 2C–E; Supplementary Fig. S3A and S3B; Supplementary Methods; ref. 8). Effects of the sgRNAs on tumor growth were assessed by examining the significance of the differences in the tumor size distribution compared with controls and also the magnitude of the effects. Tumors with each sgRNA in EGFR;p53 mice, which lack the Cas9 allele, had very similar size distributions (Fig. 2D). Furthermore, tumors initiated with each Lenti-sgInert/Cre vector or Lenti-sgp53/Cre in EGFR;p53;Cas9 mice had very similar tumor size profiles (Fig. 2C and E).

Inactivation of Rb1, Apc, and Rbm10 most dramatically enhanced the growth of oncogenic EGFR-driven Trp53-deficient lung tumors (Fig. 2C and E). Although the importance of Rb1 inactivation in EGFR-driven lung adenocarcinomas has begun to be investigated (23, 35), very few studies have addressed the functional importance of the adenomatous polyposis coli (APC) pathway on EGFR-driven lung adenocarcinomas (36). Interestingly RBM10 is an RNA binding protein and splicing regulator that is poorly studied in cancer in general and has not previously been implicated as a critical regulator of EGFR-driven lung cancer growth (12, 37–39). Rbm10 is an X-linked gene; however, we did not uncover any sex-specific effects of Rbm10 inactivation on tumor size (Supplementary Fig. S2B). More broadly, both the LN mean and the tumor size at the 95th percentile relative to tumors with inert sgRNAs were similar between male and female mice across the tumor suppressor genes studied (Supplementary Fig. S2B).

Surprisingly, inactivation of either Lkb1 or Setd2, which are strong tumor suppressor genes in analogous oncogenic KRAS-driven lung tumor models, dramatically reduced tumor growth relative to sgInert tumors (8, 40, 41). These effects were consistent across multiple percentiles within the tumor size distribution and as assessed by the LN mean of tumor sizes (Fig. 2C and E). Other tumor suppressor genes (Atm, Arid1a, Cdkn2a, and Keap1) did not significantly alter tumor growth in the context of this experiment (Fig. 2C and E). Notably, the significant effects observed with sgSetd2, sgLkb1, sgSmad4, sgApc, sgRbm10, and sgRb1 were all recaptured even after we simulated a 50% reduction in cutting efficiency of sgRNAs or when we used other strategies for subsampling underscoring the robustness of our findings (Supplementary Fig. S2C and S2D; Supplementary Fig. S4A and S4B).

To assess tumor suppressor gene function at a later time point of tumor growth, we initiated tumors in EGFR;p53;Cas9 mice with 10-fold less Lenti-sgTS^Pool/Cre and performed Tuba-seq after 19 weeks of tumor growth (Supplementary Fig. S1C–S1H). At this time point, the histology of the tumors was similar to that observed after 11 weeks of tumor initiation (Supplementary Fig. S1C and S1D). Tuba-seq analysis confirmed the tumor-suppressive function of Rbm10, Apc, and Rb1 (Supplementary Fig. S3A–S3D). Because we used a 10-fold lower viral titer for this experiment, there were proportionally fewer tumors (Supplementary Fig. S1F), which reduced the resolution of the Tuba-seq analysis. Thus, although inactivation of the other genes did not significantly affect tumor growth at this time point, we cannot exclude the possibility that these genes may influence tumor growth. Interestingly, despite the decreased statistical power, inactivation of Cdkn2a or Arid1a had a positive effect on tumor growth at this 19-week time point (but not at the 11-week time point), suggesting a potential role of these tumor suppressor genes during a later phase of tumorigenesis in this model (Supplementary Fig. S3A and S3B).

We performed further experiments to confirm the role of two less-well-studied tumor suppressors, Apc and Rbm10, on the growth of EGFR-driven tumors. We initiated lung tumors in EGFR;p53 and EGFR;p53;Cas9 mice with Lenti-sgApc/Cre, Lenti-sgNeo2/Cre (sgInert), and two Lenti-sgRbm10/Cre vectors each with a unique sgRNA targeting Rbm10 (n = 3 EGFR;p53 mice/group; n = 5 EGFR;p53;Cas9 mice/group; Fig. 3A). We used two sgRNAs targeting Rbm10 to increase the rigor of our findings and because the tumor suppressive role of RBM10 is largely uncharacterized in EGFR-driven lung cancer. Inactivation of either Apc or Rbm10 in EGFR;p53;Cas9 mice gave rise to tumors that were larger than control sgNeo2 tumors in EGFR;p53;Cas9 mice (Supplementary Fig. S5A and S5B). Moreover, by quantifying the size of individual sgApc or sgRbm10 tumors in histologic sections, we observed that EGFR;p53;Cas9 tumors were larger than tumors initiated with the same vectors in EGFR;p53 mice (Fig. 3B–E). Lenti-sgApc/Cre–initiated tumors in EGFR;p53;Cas9 mice had more cancer cells with stabilization and nuclear localization of β-catenin, as well as increased levels of SOX9, consistent with Apc inactivation (Supplementary Fig. S5C and S5D; ref. 42). Furthermore, at least 50% of Lenti-sgRbm10/Cre–initiated tumors in EGFR;p53;Cas9 mice lacked or had lower RBM10 protein levels (Supplementary Fig. S5E). Tumors with either Apc or Rbm10 inactivation were histologically similar to tumors in EGFR;p53 mice at this time point and had papillary/acinar or micropapillary structures with medium or high nuclear grade (Supplementary Fig. S6A and S6B). All tumors were positive for adenocarcinoma markers (NKX2-1/TTF1 and SFTPC) and were negative for the neuroendocrine markers synaptophysin (SYP) and ubiquitin C-terminal hydrolase L1 (UCHL1) and a squamous cell carcinoma marker (p63; Supplementary Fig. S7). Cancer cells were highly proliferative, whereas apoptotic cells were rare in tumors in EGFR;p53;Cas9 mice initiated with either Lenti-sgApc/Cre or Lenti-sgRbm10/Cre (Supplementary Fig. S7). These results further confirm the importance of these tumor suppressor genes in constraining EGFR-driven tumor growth in vivo. Collectively, our findings underscore the value of coupling Tuba-seq and CRISPR/Cas9-mediated somatic genome editing with our virally induced mouse model to dissect gene function in oncogenic EGFR-driven lung cancer.

The extent to which different oncogenic drivers affect the landscape of tumor suppression is almost entirely unknown. We approached this question by comparing the fitness landscape of tumor suppression within the contexts of oncogenic EGFR- and KRAS-driven lung tumors. We repeated an experiment previously performed by our group in which we inactivated the same panel of tumor suppressor genes in Kras^LSL-G12D/+;p53^flox/flox;R26^LSL-Tomato;H11^LSL-Cas9 (Kras;p53;Cas9) mice (29, 30, 43) and used library preparation methods and analytic pipelines identical to those used for the lungs from EGFR;p53 and EGFR;p53;Cas9 mice (Supplementary Fig. S8A). Inactivation of Lkb1, Setd2, and Rb1 were particularly strong drivers of oncogenic KRAS-driven Trp53-deficient tumor growth, whereas inactivation of Rbm10, Apc, Cdkn2a, and Arid1a also modestly increased tumor growth (Fig. 4A and B). These results are largely consistent with our previous Tuba-seq analysis of Kras;p53;Cas9 mice, as well as other studies of these genes in oncogenic KRAS–driven lung cancer mouse models (Supplementary Fig. S1A; refs. 9, 40–42, 44–48).

Inactivation of several of the tumor suppressor genes (e.g., Rb1) had similar effects on EGFR- and KRAS-driven tumors, suggesting that these tumor suppressor genes limit lung adenocarcinoma growth regardless of the oncogenic context (Fig. 4C; Supplementary Fig. S8B). However, inactivation of either Lkb1 or Setd2 greatly increased the growth of oncogenic KRAS-driven lung tumors but decreased the growth of oncogenic EGFR-driven lung tumors (Fig. 2C and E; Fig. 4A–C; Supplementary Fig. S8B). Thus, the consequences of tumor suppressor gene inactivation in specific contexts are not limited to the magnitude of tumor-suppressive effects but can also be manifested as opposite effects (known as the sign epistasis) even when the oncogenic drivers (in this case, EGFR and KRAS) are traditionally thought to be within a linear pathway.

To compare our functional data with the spectrum of tumor suppressor gene mutations found in human lung adenocarcinomas, we queried data from the American Association for Cancer Research (AACR) Project GENIE database (49). We calculated the frequency of tumor suppressor gene mutations that co-occur with oncogenic EGFR (L858R, exon 19 deletions, L861Q, G719X) or oncogenic KRAS (at codons 12, 13, or 61) mutations and TP53 mutations. This analysis revealed that RB1, RBM10, and APC are frequently altered in EGFR/TP53-mutant lung adenocarcinomas. Interestingly, RB1 mutations are more frequent in EGFR/TP53 tumors compared with KRAS/TP53 tumors (7.5% vs. 3.1%). However, Rb1 inactivation was a major driver of tumor growth in both EGFR- and Kras-mutant tumors in mice (Fig. 2C and E; Fig. 4A–C). This difference between mouse and human may be related to the higher frequency of alterations in CDKN2A in human KRAS/TP53-mutant tumors (7.3%; Fig. 4D), which would disrupt the same cell cycle pathway as RB1 inactivation. LKB1/STK11 and SETD2 are among the most frequently mutated tumor suppressor genes in KRAS/TP53-mutant lung adenocarcinomas (Fig. 4D; refs. 50, 51). Mutations in these genes occurred at significantly higher frequencies in KRAS/TP53-mutant tumors compared with EGFR/TP53 tumors, further supporting a difference in the function of LKB1 and SETD2 in EGFR/TP53- versus KRAS/TP53-mutant lung adenocarcinomas (Fig. 4E). This asymmetry in the mutation frequency of LKB1 or SETD2 within oncogenic EGFR- or KRAS-driven lung tumors is also significant when we extend our analysis to include all tumors regardless of TP53 mutation status (Supplementary Fig. S8C). Collectively, the mouse and human data indicate that the mutational patterns can reflect the consequences of profound epistasis between tumor suppressor genes and oncogenic drivers.

Genetically engineered mouse models have provided insight into the biology of EGFR-driven lung tumors and have proved valuable for studying mechanisms of resistance to EGFR TKIs, especially on-target mechanisms of resistance (27, 52–54). The TKI osimertinib was recently approved for the first-line treatment of EGFR-driven lung adenocarcinomas. However, pathways involved in modulating the depth of response and mechanisms of resistance to osimertinib are still under investigation (55, 56). To determine how tumor suppressor genes influence the therapeutic response of lung tumors to EGFR inhibition, we treated EGFR;p53;Cas9 mice with Lenti-sgTS^Pool/Cre–initiated tumors with osimertinib for 2 weeks starting at 9 weeks after tumor initiation (Fig. 5A). Osimertinib treatment greatly reduced the overall tumor burden relative to vehicle-treated EGFR;p53;Cas9 mice (Supplementary Fig. S9A–S9D). Residual neoplastic cells were sparse, as determined by staining for EGFR^L858R, and those cells were not proliferating (Supplementary Fig. S9E and S9F). The overall response was similar when the 2-week treatment was started 17 weeks after tumor initiation (Supplementary Fig. S9G–S9J). Residual tumors in EGFR;p53;Cas9 mice retained lung adenocarcinoma features and did not have evidence of phenotypic changes, such as small cell lung cancer (SCLC) transformation, after 2 weeks of osimertinib treatment (Supplementary Fig. S10A).

To quantify the impact of inactivating each tumor suppressor gene on the response to osimertinib in vivo, we performed Tuba-seq on the lungs from osimertinib-treated EGFR;p53;Cas9 mice 11 and 19 weeks after tumor initiation and compared the results to Tuba-seq results from vehicle-treated controls (Fig. 5B; Supplementary Fig. S11A and S11B). Consistent with the imaging data and histologic analysis, osimertinib treatment greatly reduced tumor burden as assessed by Tuba-seq (compare Fig. 2C and E, Fig. 5B, Supplementary Fig. S3A and S3B, and Supplementary Fig. S11A and S11B; Supplementary Methods).

After 2 weeks of osimertinib treatment, inactivation of Apc, Rb1, or Rbm10 was still associated with larger tumors, whereas the size distribution of tumors with inactivation of other tumor suppressor genes remained similar to that of tumors with inert sgRNAs (compare Fig. 5B and Fig. 2E; Fig. 2C and Supplementary Fig. S11A; Supplementary Fig. S11B and Supplementary Fig. S3A and S3B). The striking exception was tumors with sgKeap1. In vehicle-treated mice, the size distribution of sgKeap1 tumors was almost identical to that of tumors with inert sgRNAs; however, in osimertinib-treated mice, sgKeap1 tumors were significantly larger than the tumors with inert sgRNAs. This suggests that inactivation of Keap1 limits responses to osimertinib (Fig. 5B). Osimertinib resistance conferred by Keap1 inactivation was also observed at 19 weeks after tumor initiation (Supplementary Fig. S11A and S11B).

We applied an analytic approach that we previously developed and validated to quantify the genotype-specific responses (Fig. 5C; Supplementary Fig. S11C–S11G; Supplementary Methods; ref. 34). By comparing the LN mean of the observed tumor size distributions in osimertinib-treated mice with the expected tumor size distribution based on the overall drug effects, we can estimate genotype-specific drug responses (ScoreRLM; Fig. 5C; Supplementary Fig. S11C). At 11 weeks after tumor initiation, following 2 weeks of osimertinib treatment, sgRb1 tumors were 25% smaller than expected (P = 0.04). Conversely, tumors with Keap1 inactivation were 48% larger than expected (P = 0.07; Fig. 5C). The effect of Keap1 inactivation was even greater at 19 weeks after tumor initiation, when tumors were 274% larger than expected (P = 0.13; Supplementary Fig. S11C). Given the magnitude of the ScoreRLM for sgKeap1 at both time points (ScoreRLM = 0.57 and 1.90 after 2 weeks of treatment at 11 and 19 weeks after tumor initiation, respectively), we combined the two independent P values and confirmed that Keap1 inactivation significantly reduced the therapeutic response to osimertinib (Fisher method, P = 0.05). Other statistical measures of genotype-specific responses, including relative tumor number (ScoreRTN) and relative geometric mean (ScoreRGM), did not significantly differ between the treated and untreated groups (Supplementary Fig. S11D; Supplementary Methods). Our analytic methods allow us to uncover when effects are more pronounced on larger tumors (ScoreRTN and ScoreRGM have much lower sensitivity when the effects are greater on larger tumors; Supplementary Fig. S11E–S11G; Supplementary Methods). Thus, our data are consistent with resistance conferred by Keap1 inactivation being more pronounced in larger tumors.

To further investigate these findings, we initiated tumors with Lenti-sgKeap1/Cre in EGFR;p53 and EGFR;p53;Cas9 mice followed by treatment with osimertinib or vehicle (Fig. 5D and E). Osimertinib treatment reduced the size and number of tumors in EGFR;p53 mice compared with vehicle-treated EGFR;p53 mice (Fig. 5D; Supplementary Fig. S12A and S12B). Conversely, osimertinib treatment of Lenti-sgKeap1/Cre–initiated tumors in EGFR;p53;Cas9 mice did not decrease tumor size or number compared with vehicle-treated EGFR;p53;Cas9 mice (Fig. 5D and E; Supplementary Fig. S12A and S12B). Consistent with the inefficiency of CRISPR/Cas9-mediated genome editing in somatic cells, some tumors initiated with Lenti-sgKeap1/Cre in EGFR;p53;Cas9 mice retained KEAP1 protein (Fig. 5E, top panels). However, the tumors that remained in osimertinib-treated EGFR;p53;Cas9 mice all had medium to low or negative KEAP1 levels (Fig. 5E). Together, these data indicate that, although Keap1 inactivation is not positively selected for during oncogenic EGFR-driven tumor growth, osimertinib treatment selects for cancer cells with low or absent KEAP1, thus reducing the therapeutic response to the drug.

To correlate our findings with clinical data, we analyzed the effects of KEAP1 pathway alterations on patient outcomes to EGFR inhibition in EGFR/TP53-mutant lung adenocarcinomas. Oncogenic EGFR–driven tumors with mutations in genes in the KEAP1 pathway have been suggested to be less responsive to TKIs, and we confirmed this association in EGFR/TP53 tumors (Fig. 6A and B; ref. 57). Mutations in KEAP1/NFE2L2/CUL3 were associated with a significantly shorter time on EGFR TKI therapy in contrast to KEAP1/NFE2L2/CUL3 wild-type tumors (5.8 vs. 14.3 months; P = 0.01, log-rank test; Fig. 6A; Supplementary Tables S1 and S2). This remained significant even after adjustment for potential confounders such as age, sex, race, and smoking status (Supplementary Table S3). Among several other tumor suppressor genotypes, KEAP1 pathway alterations were the most significant driver of limited sensitivity after correction for multiple hypothesis testing (Fig. 6B).

We also analyzed a dataset of oncogenic EGFR lung adenocarcinoma patient samples collected through the Yale Lung Rebiopsy Program prior to first-line TKI treatment and after the development of resistance to TKIs (Fig. 6C; Supplementary Table S4). Among 18 patients with EGFR and TP53 mutant tumors, two had KEAP1 alterations at relapse (TKI-R samples). One case had an acquired missense mutation (Y206N) that lies in a domain of KEAP1 involved in forming the complex with Cullin3 to mediate ubiquitination and degradation of NRF2, which is encoded by NFE2L2 (58). The other tumor, which was analyzed only at resistance, had heterozygous loss of KEAP1 (Fig. 6C). Moreover, we also observed two cases of NFE2L2 copy-number gain (one detected prior to treatment and maintained at resistance, the other detected only at resistance), and one case of CUL3 heterozygous loss at resistance. Taken together, our in vivo functional results are consistent with human data and support a role for alterations in the KEAP1 pathway in reducing EGFR TKI sensitivity.

The genomic landscape of lung adenocarcinoma is combinatorially complex, with a large number of genomic alterations in oncogenic drivers and tumor suppressor genes occurring in various combinations (2, 3). Although critical for precision medicine approaches, dissecting the functional contribution of these co-occurring genomic alterations to tumor growth and response to therapy is extremely challenging. Using somatic CRISPR/Cas9-mediated gene editing and Tuba-seq (8, 9), we evaluated the fitness effects of inactivating 10 tumor suppressor genes that are commonly altered in lung adenocarcinoma. We uncovered distinct roles for specific tumor suppressor genes during tumor growth in the different oncogenic contexts and on sensitivity to the TKI osimertinib in vivo.

EGFR and KRAS are the most frequently mutated oncogenic drivers in lung adenocarcinoma, collectively occurring in ∼40% to 50% of all tumors (2). Oncogenic alterations in these two genes are mutually exclusive, consistent with them being in the same canonical receptor tyrosine kinase EGFR/RAS/MAPK signaling pathway (59, 60). For this reason, one might have anticipated that inactivating tumor suppressor pathways in EGFR-driven and KRAS-driven lung cancers would have similar biological consequences. It is tempting to make a simple assumption that tumor suppressor genes have a constant marginal effect independent of other genomic alterations, or, in other words, that the fitness landscape of tumor suppression is smooth and that epistatic interactions are infrequent and/or small in magnitude. Here, we tested this assumption by comparing the fitness landscape of tumor suppression in EGFR-driven Trp53-deficient and KRAS-driven Trp53-deficient lung adenocarcinoma and found distinct patterns of epistasis between tumor suppressor genes and oncogenes. Inactivation of Rb1, Rbm10, and Apc had a similar effect on EGFR- and KRAS-driven lung tumors (Fig. 2C and E; Fig. 4A–C; Supplementary Fig. S8B). In contrast, although Lkb1 and Setd2 were among the most potent tumor suppressor genes in Kras/p53 tumors, inactivation of Lkb1 or Setd2 led to reduced tumor growth of EGFR/p53 tumors (Fig. 2C and E; Fig. 4A and B). This observation is consistent with human data, where alterations in either LKB1 or SETD2 are much less common in lung tumors with oncogenic EGFR than in tumors with oncogenic KRAS (Fig. 4D and E; Supplementary Fig. S8C). Thus, we reveal a highly context-dependent fitness landscape of tumor suppression that depends on the nature of the oncogenic driver. Investigation of the mechanisms that underlie the sign epistasis of LKB1 and SETD2 may uncover new biological insights and vulnerabilities of EGFR-mutant tumors. More broadly, sign epistasis within seemingly similar cancer contexts could help identify genetic interactions for further functional investigation and should be considered when interpreting cancer genomic data.

Concurrent TP53 and RB1 mutations have been linked to the transformation of lung adenocarcinomas to SCLC that is observed in some EGFR-mutant tumors upon the emergence of TKI resistance (23, 61, 62). Additional histologic changes have also been described in TKI-resistant tumors, including squamous cell carcinoma differentiation (63). In our screen and validation studies, we found that the tumors had lung adenocarcinoma features and did not express neuroendocrine or squamous cell carcinoma markers (Supplementary Figs. S7 and S10A). We also examined a small cohort of EGFR;p53;Cas9 mice with tumors initiated using Lenti-sgRb1/Cre, and these tumors were positive for adenocarcinoma markers and negative for neuroendocrine and squamous cell carcinoma markers at baseline (Supplementary Fig. S10B). Long-term TKI treatment studies in these new models will provide us with the opportunity to determine how tumors on different tumor suppressor gene backgrounds evolve through TKI treatment and whether specific tumor suppressor genes contribute to altering the histologic properties of TKI resistant tumors.

A major advantage of autochthonous genetically engineered mouse models of human cancer is that they can be used to study the impact of inactivation of putative tumor suppressor genes on therapy responses in vivo (10, 34, 64). We found that inactivation of Keap1 decreases sensitivity to osimertinib in vivo (Fig. 5). This finding was supported in a clinical cohort of EGFR/TP53-mutated lung adenocarcinomas, in which patients with tumors harboring mutations in KEAP1/NFE2L2/CUL3 had a significantly shorter time to treatment discontinuation with EGFR TKIs compared with matched controls (Fig. 6A and B; Supplementary Fig. S13A–S13D; refs. 17, 57). These results are consistent with inactivation of this tumor suppressor gene limiting the response to therapy and explain the occasional presence of alterations in KEAP1 pathway genes in TKI-resistant human tumors.

The reduced sensitivity of tumors initiated with Lenti-sgKeap1 may be associated with the accumulation and transactivation of genes regulated by NRF2 (65–67). Indeed, we found that KEAP1-deficient tumors had increased nuclear NRF2 and NQO1 expression relative to KEAP1-proficient tumors, suggesting the presence of enhanced NRF2 transcriptional activity (Supplementary Fig. S12C; refs. 65–67). Activation of both the MAPK and AKT/mTOR pathways downstream of EGFR were strongly decreased in both KEAP1-deficient and KEAP1-proficient tumors after 2 weeks of treatment with osimertinib, thus excluding the possibility that the reduced sensitivity could be due to sustained activation of these pathways (Supplementary Fig. S12D). NRF2 is a pleiotropic transcription factor that can lead to increased production of antioxidants, regulate ferroptosis, and alter cellular metabolism (68). Which of the functional outputs of activation of the NRF2 pathway are responsible for the decreased sensitivity to osimertinib in EGFR-mutant tumors remains to be determined. Experiments in lung adenocarcinoma cell lines treated with targeted therapies against the EGFR/MAPK pathway suggest that KEAP1 loss leads to decreased levels of drug-induced reactive oxygen species and increased glutathione synthesis (66). A systematic analysis of the KEAP1 pathway in EGFR-mutant genetically engineered mouse models and patient-derived tumors will be necessary to fully elucidate how NRF2 mediates TKI resistance and reveal potential therapeutic vulnerabilities. Indeed, there are potential therapeutic avenues for treating KEAP1-deficient tumors that should be explored in EGFR-mutant tumors. In KRAS-driven lung tumors, KEAP1/NRF2-dependent metabolic reprogramming leads to a dependency on glutaminolysis, which can be targeted via glutaminase inhibition (69–73). Future studies will be important to establish whether similar metabolic rewiring occurs in KEAP1-deficient EGFR-mutant tumors that can be leveraged therapeutically.

Genetic alterations in KEAP1 occur in less than 2% of EGFR-mutant tumors at baseline and are found in less than 10% of TKI-resistant EGFR-mutant lung adenocarcinomas (Fig. 6C; Supplementary Fig. S13E; ref. 17). It is possible that non-genetic alterations that increase NRF2-dependent gene expression also occur in EGFR-mutant tumors and that examination of the mutational status of the pathway underestimates the frequency with which it is disrupted (17, 74–78). Collectively, these results raise the possibility that NRF2 should be investigated as a potential target to reduce or delay the onset of drug resistance in EGFR-driven lung adenocarcinomas. More broadly, our study demonstrates that our approach can be used to identify clinically relevant pathways that modulate response to therapy in vivo. By uncovering the driving forces of the heterogeneity of responses to therapy observed in patients, these types of studies could help define high-risk versus low-risk patient populations and guide therapeutic interventions.

This study provides insight into the complex interplay between tumor suppressor genes and other co-occurring mutations in EGFR-driven lung adenocarcinoma tumorigenesis and thus has significant clinical implications. Across many cancers, including lung adenocarcinoma, most tumor suppressors are mutated in a low percentage of tumors. Several of the tumor suppressor genes that we investigated were among the most frequently mutated genes in EGFR-driven TP53-mutant lung cancer, such as RBM10 and RB1 (Fig. 4D; Supplementary Fig. S14A; refs. 79, 80). However, even the tumor suppressor genes that we found to have an effect on in vivo tumor growth are not as frequently altered in human tumors as might be expected from the experimental findings. This could be related to the overall higher mutation burden in human tumors versus mouse models (81), with human tumors harboring alterations in genetically redundant pathways (cell-cycle regulation, WNT pathway) or inactivation of tumor suppressor genes by epigenetic mechanisms (11, 36, 82, 83). Finally, the type of mutational process that gives rise to the mutations could also affect the frequency with which certain tumor suppressor genes are mutated (84, 85). Tuba-seq allows the analysis of genetic effects without being confounded by the frequency with which alterations in that gene are generated by natural processes in human tumors.

The lack of genomic complexity in the genetically engineered mouse models (81) allows us to evaluate interactions between co-occurring alterations by avoiding confounding factors pervasive in human genomic data (e.g., tumor mutation load, mutation frequency, additional mutations) and environmental factors such as smoking. Our data provide clear quantitative data on mutual exclusivity and synergistic biological effects of genetic alterations that mirror observations in patient specimens, highlighting the strength of these interactions. Additional studies in models with a higher tumor mutation burden, more reflective of that present in human tumors, will be important to uncover the extent to which additional genetic alterations further modulate tumor growth and drug sensitivity in specific oncogenic driver/tumor suppressor gene contexts. Notably, tumors with other oncogenic drivers (e.g., ALK rearrangements) also have unique spectra of co-occurring tumor suppressor gene alterations, further suggesting widespread interactions between tumor suppressor gene pathways and oncogenic drivers (Supplementary Fig. S14B). Future in vivo Tuba-seq studies should investigate tumor models driven by other oncogenes to uncover a broader understanding of the genetic interactions between diverse oncogenes and large panels of tumor suppressor genes. Precise mapping of the fitness consequences of combinations of genetic alterations during tumor evolution will help uncover the biological and clinical relevance of specific alterations during carcinogenesis and identify pathways that can be exploited as therapeutic targets to prevent or overcome resistance to targeted therapies.

TetO-EGFR^L858R, p53^flox/flox, Rosa26^{CAGs-LSL-rtTA3-IRES-mKate}, Rosa26^{CAGs-LSL-Cas9-GFP}, Kras^LSL-G12D, Rosa26^LSL-Tomato, and H11^LSL-Cas9 mice have been described (4, 25, 27–30, 43, 86, 87). EGFR;p53 and EGFR;p53;Cas9 were on a mixed BL6/129/FVB background and Kras;p53;Cas9 mice were on a mixed BL6/129 background. Approximately equal numbers of males and females were used for each experiment, and the number of mice used for each experiment is listed in each figure legend. Lung tumors were initiated by intratracheal administration of lentiviral Cre vectors as previously described (26). Tumor burden was assessed by MRI, fluorescence microscopy, lung weight, and histology, as indicated. Doxycycline was administered by feeding mice with doxycycline-impregnated food pellets (625 ppm; Harlan-Teklad). Osimertinib (AstraZeneca) was resuspended in 0.5% (w/v) methylcellulose (vehicle) and was administered orally (25 mg/kg 5 days a week). All animals were kept in pathogen-free housing under guidelines approved by either the Yale University Institutional Animal Care or the Stanford University Institutional Animal Care and Use Committee (IACUC) guidelines.

The barcoded vectors in the Lenti-sgTS^Pool/Cre have been previously described (Supplementary Table S5; ref. 8). The second Lenti-sgRbm10/Cre vector used in the validation experiments was generated by site-directed mutagenesis (Supplementary Table S6). Briefly, lentiviral U6-sgRNA/Cre vectors contain an eight-nucleotide defined sequence (sgID) that identifies the sgRNA followed by a 15-nucleotide random barcode (BC) to uniquely tag each tumor (8). To avoid BC–sgRNA uncoupling driven by lentiviral template switching during reverse transcription of the pseudo-diploid viral genome, each barcoded Lenti-sgRNA/Cre vector was generated separately (88, 89). We cultured HEK293T cells in Dulbecco's Modified Eagle Medium with 10% fetal bovine serum and transfected them with individual barcoded Lenti-sgRNA/Cre plasmids (sgLkb1, sgp53, sgApc, sgAtm, sgArid1a, sgCdkn2a, sgKeap1, sgNeo1, sgNeo2, sgNeo3, sgNT1, sgRb1, sgRbm10, sgRbm10#2 unbarcoded, sgSetd2, or sgSmad4) along with pCMV-VSV-G (8454; Addgene) envelope plasmid and pCMV-dR8.2 dvpr (8455; Addgene) packaging plasmid using polyethylenimine. We treated the cells with 20 mmol/L sodium butyrate 8 hours after transfection, changed the culture medium 24 hours after transfection, and collected supernatants 36 and 48 hours after transfection. Subsequently, we removed the cell debris with a 0.45 μm syringe filter unit (SLHP033RB; Millipore Sigma), concentrated each lentiviral vector by ultracentrifugation (25,000 g for 1.5 hours at 4°C), resuspended the virus in phosphate-buffered saline, and stored the virus at −80°C. To determine the titer of each vector, we transduced Rosa26^LSL-YFP mouse embryonic fibroblasts (a gift from Dr. Alejandro Sweet-Cordero, University of California, San Francisco), determined the percentage of yellow fluorescent protein–positive cells by flow cytometry, and normalized the titer to a lentiviral preparation of known titer. The cell line was authenticated at the University of California, San Francisco, and tested negative for Mycoplasma contamination. Lentiviral vectors were thawed and pooled immediately prior to delivery to mice. All of these plasmids are available at https://www.addgene.org/Monte_Winslow/.

Anticipated growth rates were determined by monitoring tumor development through MRI in pilot experiments, and the analysis time points were selected to ensure that tumors were detectable by MRI such that their response to treatment could be evaluated. Viral titers used in the experiments were chosen to balance the total tumor burden across mice at the time of analysis after tumor initiation. For analysis of tumor growth 11 weeks after tumor initiation, the Lenti-sgTS^Pool/Cre titer administered to EGFR;p53 mice was 2 × 10⁶ infectious units/mouse, whereas for EGFR;p53;Cas9 mice we used 1 × 10⁶ infectious units/mouse. We reasoned that using a higher viral titer in the control EGFR;p53 mice would increase our confidence that any differences observed between EGFR;p53 and EGFR;p53;Cas9 mice were due to inactivation of tumor suppressor genes in the latter model. At the 19-week time point, in EGFR;p53;Cas9 mice we initiated tumors with 1 × 10⁵ infectious units/mouse. Two weeks before collection, mice were treated with either vehicle or osimertinib. For the validation experiments in which we used a single vector to initiate tumors—Lenti-sgApc/Cre, Lenti-sgRbm10/Cre, Lenti-sgRbm10#2/Cre, or Lenti-sgNeo2/Cre (sgInert)—we used 1 × 10⁵ infectious units/mouse and analyzed the mice after 14 weeks of tumor growth (Supplementary Table S6). For the validation with Lenti-sgKeap1/Cre virus (1 × 10⁵ infectious units/mouse), both EGFR;p53 and EGFR;p53;Cas9 mice were treated 15 weeks after tumor initiation, and lungs were collected after 2 weeks of treatment with either vehicle or osimertinib. Kras;p53;Cas9 were analyzed 14 weeks after tumor initiation with 2.2 × 10⁴ infectious units/mouse.

All procedures were performed in accordance with protocols approved by the Yale University IACUC and in agreement with the National Institutes of Health Guide for the Care and Use of Laboratory Animals. Respiratory gated, gradient-echo MRI images of mice were collected with a 4T (31-cm bore) small-animal AVANCE horizontal-bore spectrometer (Bruker). All data were collected as previously described (52). Tumor volume was quantified by calculating the area of visible lung opacities present in each image sequence per mouse using BioImage Suite 3.01 (90). For the drug treatment studies, when tumors were detected using MRI, mice were assigned to different groups to ensure a comparable tumor volume distribution prior to treatment between the vehicle-treated and the osimertinib-treated group.

Genomic DNA was isolated from bulk tumor–bearing lung tissue from each mouse as previously described (8). Briefly, three benchmark control cell lines (∼5 × 10⁵ cells per cell line) carrying unique sgID-BCs were added (“spiked-in”) to each sample prior to lysis to enable the calculation of the absolute number of neoplastic cells in each tumor from the number of sgID-BC reads. Following homogenization and overnight protease K digestion, genomic DNA was extracted from the lung lysates using standard phenol–chloroform and ethanol precipitation.

sgID-BC sequencing libraries were prepared by polymerase chain reaction (PCR) amplifying the sgID-BC region from total genomic DNA. To enable the identification and subsequent computational elimination of index-hopped reads after high-throughput sequencing, the sgID-BC region of the integrated Lenti-sgRNA-BC/Cre vectors was PCR amplified using unique dual-indexing primer pairs. To increase the sequence diversity, we added 6 to 9 Ns before the sequence-specific primers (34). For each mouse, we performed eight 100-μL PCR reactions per sample (4 μg DNA per reaction, 32 μg per mouse) using the Q5 High-Fidelity 2X Master Mix (M0494X; New England Biolabs). The PCR products were purified with Agencourt AMPure XP beads (A63881; Beckman Coulter) using a double-size selection protocol. The concentration and quality of the purified libraries were determined using the Agilent High Sensitivity DNA Kit (5067-4626; Agilent Technologies) on the Agilent 2100 Bioanalyzer (G2939BA; Agilent Technologies). The libraries were pooled based on lung weight (to have sequencing depth more evenly distributed across samples), cleaned up, size-selected using AMPure XP beads, and sequenced on the Illumina HiSeq 2500 platform to generate paired-end 150-bp reads (Admera Health). To increase the diversity at the beginning of sequencing reads and improve sequencing quality, 5% to 15% PhiX control DNA was added to the library.

The auxiliary lobe of the right lung was collected for each experimental mouse from the mice transduced with the Lenti-sgTS^Pool/Cre pool, whereas both the auxiliary and left lobes were collected for the validation experiments with individual sgRNAs (sgApc, sgRbm10, sgRbm10#2, and sgNeo2). Right and auxiliary lobes were collected for the experiment with Lenti-sgKeap1/Cre virus. Lung lobes were fixed in 4% paraformaldehyde overnight at room temperature, placed in 70% ethanol, and paraffin-embedded and sectioned (Histology@Yale). Four-micrometer sections were used for hematoxylin and eosin (H&E) staining and immunohistochemistry. Tumor sizes were determined by measuring the longest diameter for each tumor in H&E-stained sections. Tumor size and tumor area were quantified using ImageJ (National Institutes of Health). The limited tissue collected for histologic analysis reduced the overall number of tumors that could be measured. P values were calculated from the Mann–Whitney U test. The following antibodies were used for immunohistochemistry: anti-EGFR^L858R (1:200, 3197; Cell Signaling Technology), anti-SP-C (SFTPC; 1:200, ab40876; Abcam), anti-TTF1/NKX2–1 (1:200, ab76013; Abcam), anti–phospho-histone H3 (1:200, 9701; Cell Signaling Technology), anti–cleaved caspase-3 (1:300, 9664; Cell Signaling Technology), anti-mKate (1:2500, AB233; Evrogen), anti-Cas9 (1:500, 7A9-3A3 and NBP2-36440; Novus Biologicals), anti-β-catenin (1:500, 8814; Cell Signaling Technology), anti-RBM10 (1:200, ab224149; Abcam), anti-KEAP1 (1:500, ab227828; Abcam), anti-SOX9 (1:500, AB5535; Merck Millipore), anti-NRF2 (1:500, ab137550; Abcam), anti-NQO1 (1:500, AB28947; Abcam), anti-phospho-ERK (1:100, 4370; Cell Signaling Technology), anti-phospho-AKT (1:100, 4060; Cell Signaling Technology), anti-Phospho-S6 Ribosomal Protein (P-S6; 1:1000, 5364; Cell Signaling Technology), anti-UCHL1 (1:500, HPA005993; Sigma Millipore), anti-Synaptophysin (1:400, ab32127; Abcam), and anti-p63 (1:2500, ab124762; Abcam). Nuclear protein expression in tumors was binned as high (over 75% of positive nuclei in a tumor area), medium (between 25% and 75% of positive nuclei), or low/negative (positive nuclear staining below 25%). Cytosolic protein levels were binned as high (more than 75% of positive cells), medium (between 25% and 75% of positive cells), or low/negative (positive staining below 25%). Statistical analyses were performed with GraphPad Prism 8 and QuPath.

The AACR Project GENIE is a registry that contains Clinical Laboratory Improvement Amendments– and International Organization for Standardization–certified genomic data collected from the records of more than 9,000 patients who were treated at each of the consortium's participating institutions (49). Data from GENIE 4.1-public were accessed through the Synapse Platform. We generated a comprehensive list of all missense, nonsense, and frameshift mutations for all screened genes across all participating centers in the GENIE project, documenting these mutations at the gene and amino acid levels. Based on the sets of genes included in the different screening panels that contribute to Project GENIE, we annotated all lung adenocarcinomas within the database as being wild-type, mutated, or not screened for each gene. From this complete catalog of mutations in each tumor, we determined the rates of co-occurrence of known oncogenic KRAS mutations (G12X, G13X, Q61X) or EGFR mutations (L858R, exon 19 deletions, L861Q, G719X) with missense, nonsense, or frameshift mutations in our set of 10 tumor suppressor genes. Furthermore, we determined the frequency of coincident tumor suppressor gene mutations in tumors with oncogenic KRAS and inactivating TP53 mutations (KRAS/TP53) and in tumors with oncogenic EGFR and inactivating TP53 mutations (EGFR/TP53) using GENIE 8.1 (Supplementary Fig. S14A). A similar analysis was also performed on tumors with ALK rearrangements and TP53 mutations using GENIE 8.1 (Supplementary Fig. S14B). The pipelines for these analyses are publicly available at www.github.com/dgmaghini/GENIE and https://github.com/lichuan199010/Tuba-seq-analysis-and-summary-statistics, and they allow users to input OncoTree codes of interest, generate mutational profiles for the corresponding GENIE tumors, and identify the co-occurrence of up to two background mutations (e.g., EGFR and TP53) with mutations in other genes of interest. Analyses were performed using Phyton 3.6 and R 3.6.

Patients with lung cancer who were evaluated at the Stanford Cancer Center and had their tumors analyzed using the Stanford Solid Tumor Actionable Mutation Panel (STAMP; ref. 91) were included in the analysis after obtaining written informed consent. This retrospective study was conducted under a molecular analysis protocol approved by the Stanford University Institutional Review Board. All STAMP cases performed between 2015 and 2019 were included; during this time, there were two different assays used, with one covering 198 genes (302 kb) and the other covering 130 genes (232 kb). STAMP was performed as the standard of care and thus at the discretion of the treating physician and could have occurred at the time of diagnosis or at the time of progression.

From the available STAMP cases, patients were selected if they had stage IV non–small cell lung cancer (NSCLC) and their tumor had a pathogenic EGFR and TP53 mutation. Patients were excluded if they had incomplete data or were lost to follow-up prior to analysis of the primary endpoint, if they elected not to receive treatment, or if they received adjuvant TKI therapy for early-stage disease. Within this cohort, patients were further selected for the presence of at least one of the tumor suppressor genes investigated in the preclinical setting (KEAP1, LKB1/STK11, SETD2, SMAD4, RB1, APC, ARID1A, and/or CDKN2A). To ensure that there was enough power to conduct an analysis, genes were analyzed individually if there were at least eight patients with tumors with a mutation of that gene who met the inclusion criteria. Due to the small numbers, patients with tumors with LKB1 (n = 1), SETD2 (n = 0), and SMAD4 (n = 1) mutations were not analyzed separately. However, as two tumors harboring LKB1 and SMAD4 mutations were also EGFR/TP53 mutant, they were included as part of the control arm. Univariate analysis of time to treatment failure (TTF) was completed on cases with stage IV, EGFR/TP53-mutant NSCLC that were treated with a TKI (13) and had mutations in a KEAP1 pathway component (KEAP1, NFE2L2, CUL3), RB1, APC, ARID1A, or CDKN2A. Mutations that were significant on univariate analysis were identified and multivariate analyses accounting for (i) cofounding variables and (ii) comutations were run. Demographic data including sex, age at diagnosis, smoking history, and ethnicity were extracted (Supplementary Tables S1–S3). For each tumor suppressor, patients were matched 1:3 with a wild-type cohort on the basis of sex, smoking history, ethnicity, age, and treatment type. TTF was determined by subtracting the date of discontinuation of TKI due to progression, toxicity, or death from the date of initiation of TKI and is reported in months. Statistical analysis was performed using GraphPad Prism 8 and R studio. The Kaplan–Meier method was used to estimate TTF. Comparison of survival curves was made using the log-rank test. Significance was defined as P < 0.05. Hazard ratios (HR) were generated from multivariate regression analysis performed in R and reported with a 95% CI.

Patients with EGFR-mutant lung adenocarcinoma were enrolled to a Yale University institutional review board–approved protocol after obtaining written informed consent allowing the collection and analysis of clinical data, use of archival and fresh tissue and blood, and the generation of patient-derived xenografts. This study was conducted in accordance with recognized ethical guidelines (Declaration of Helsinki, CIOMS, Belmont Report, U.S. Common Rule). Patients who received targeted therapy (erlotinib, gefitinib, or afatinib) as first-line therapy, either alone or in combination with other therapies, such as chemotherapy or cetuximab, were included (Supplementary Table S4). For genomic studies, formalin-fixed paraffin-embedded tissue was macrodissected to enrich for tumor material. All of the tumor samples (n = 29) analyzed had matched normal tissue and cancer cell purity > 20% as assessed by the loss of heterozygosity (LOH) events from the whole-exome sequencing data.

DNA was extracted and analyzed as previously described (92). Briefly, genomic DNA was captured on the Roche NimbleGen 2.1M human exome array and subjected to 74-bp paired-end reads on the Illumina HiSeq2000 instrument. The mean coverage for normal was 109.1×, and the mean coverage for tumor was 189×, with 92.78% and 95.56% for the bases covered by at least 20 independent sequence reads, respectively. Sequence reads were mapped to the human reference genome (GRCh37) using the Burrows–Wheeler Aligner-MEM (BWA-MEM) program. Sequence reads outside the targeted sequences were discarded, and the statistics on coverage were collected from the remaining reads using in-house per scripts.

For all matched tumor–normal pairs, somatic point mutations and indels were called by MuTect2 using Bayesian classifiers. For all somatic mutations called, we extracted base coverage information in all samples and considered the mutations that were supported by at least two independent sequence reads covering non-reference alleles and present in more than 5% of all sequencing reads. Identified variants were further filtered based on their presence in repositories of common variations (1000 Genomes, National Heart, Lung, and Blood Institute exome variant server, and 2,577 noncancer exomes sequenced at Yale) and annotated using the ANNOVAR program (93). All somatic indels were visually inspected to remove the false-positive calls.

Copy-number analysis was performed as previously described (92). Briefly, copy-number variants were identified from the whole-exome sequencing data using EXCAVATOR software, which normalizes the nonuniform whole-exome sequencing data taking guanine-cytosine content, mappability, and exon size into account (94). The hidden Markov model was utilized to classify each copy number variant segment into five copy number states (homozygous deletion, heterozygous deletion, normal copy number, homozygous copy gain, or multiple copy gain). Tumor purity was estimated from LOH using in-house scripts.

The datasets generated during and/or analyzed in the current study are available in the NCBI Gene Expression Omnibus database (GSE146550) and dbGAP (phs002334.v1.p1; https://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/study.cgi?study_id=phs002334.v1.p1). The GENIE genomic data analysis is publicly available at www.github.com/dgmaghini/GENIE. More GENIE analyses and other data supporting the findings are available at https://github.com/lichuan199010/Tuba-seq-analysis-and-summary-statistics or upon request.

G. Foggetti reports grants from the American-Italian Cancer Foundation and non-financial support from AstraZeneca during the conduct of the study. H. Cai reports grants from the University of California during the conduct of the study. S.N. Gettinger reports grants from the National Cancer Institute–National Institutes of Health during the conduct of the study. M. Diehn reports grants and personal fees from AstraZeneca; personal fees from BioNTech, Genentech, Gritstone Oncology, Novartis, RefleXion, and Roche Sequencing Solutions; personal fees and non-financial support from Illumina; other support from CiberMed and Foresight Diagnostics outside the submitted work; and a patent for cancer biomarkers pending, issued, and with royalties paid from Roche. H.A. Wakelee reports grants and personal fees from AstraZeneca, Helsinn, and Xcovery; grants from ACEA Biosciences, Arrys Therapeutics, Bristol Myers Squibb, Celgene, Clovis Oncology, Genentech/Roche, Gilead, Merck, Novartis, Pharmacyclics, and Seattle Genetics; and personal fees from Blueprint, Daiichi Sankyo, Janssen, and Mirati outside the submitted work, as well as unpaid consultancy work with Genentech/Roche and Merck. D.A. Petrov reports other from D2G Oncology outside the submitted work and a patent for compositions and methods for multiplexed quantitative analysis of cell lineages issued, licensed, and with royalties paid from D2G Oncology. M.M. Winslow reports grants from the National Institutes of Health during the conduct of the study and is the founder of, and an equity holder in, D2G Oncology. K. Politi reports grants from the National Cancer Institute–National Institutes of Health during the conduct of the study; grants and personal fees from AstraZeneca; grants from Kolltan, Roche/Genentech, Boehringer Ingelheim, and Symphogen; and personal fees from Dynamo Therapeutics, Halda, Maverick Therapeutics, and Tocagen outside the submitted work; and a patent for EGFR^T790M mutation testing issued, licensed, and with royalties paid from Molecular Diagnostics/Memorial Sloan Kettering Cancer Center. No disclosures were reported by the other authors.

G. Foggetti: Data curation, formal analysis, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. C. Li: Data curation, formal analysis, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. H. Cai: Data curation, validation, investigation, methodology, writing–original draft, writing–review and editing. J.A. Hellyer: Data curation, formal analysis, investigation, methodology, writing–review and editing. W. Lin: Data curation, investigation, methodology, writing–review and editing. D. Ayeni: Investigation, methodology, writing–review and editing. K. Hastings: Investigation, methodology, writing–review and editing. J. Choi: Data curation, formal analysis, methodology, writing–review and editing. A. Wurtz: Data curation, methodology, writing–review and editing. L. Andrejka: Data curation, validation, methodology, writing–review and editing. D.G. Maghini: Data curation, formal analysis, investigation, methodology, writing–review and editing. N. Rashleigh: Methodology, writing–review and editing. S. Levy: Methodology, writing–review and editing. R. Homer: Formal analysis, validation, writing–review and editing. S.N. Gettinger: Data curation, writing–review and editing. M. Diehn: Data curation, writing–review and editing. H.A. Wakelee: Data curation, writing–review and editing. D.A. Petrov: Conceptualization, resources, supervision, funding acquisition, writing–original draft, writing–review and editing. M.M. Winslow: Conceptualization, resources, supervision, funding acquisition, visualization, writing–original draft, writing–review and editing. K. Politi: Conceptualization, resources, supervision, funding acquisition, writing–original draft, writing–review and editing.

The authors would like to acknowledge the American Association for Cancer Research and its financial and material support in the development of the AACR Project GENIE registry, as well as members of the consortium for their commitment to data sharing. Interpretations are the responsibility of study authors. We thank members of the Politi, Winslow, and Petrov laboratories for helpful feedback and support, particularly Jacqueline Starrett and Emily Forrest for editing and Fernando de Miguel for creating resources to help with the whole-exome sequencing analysis. We also thank Henning Stehr for providing the sequencing data from the STAMP cohort and the Bosenberg Laboratory for sharing their Leica M205 FA stereomicroscope. G. Foggetti was supported by an American-Italian Cancer Foundation postdoctoral fellowship. C. Li is a Connie and Bob Lurie Fellow of the Damon Runyon Cancer Research Foundation (DRG-2331). H. Cai was supported by a Tobacco-Related Disease Research Program Postdoctoral Fellowship (28FT-0019). D. Ayeni was supported by National Institutes of Health (NIH) F31 CA203488. K. Hastings was supported by NIH F32 CA210516. This work was supported by NIH R01 CA231253 and R01 CA234349 (M.M. Winslow and D.A. Petrov), NIH R01 CA120247 and P50 CA196530 (K. Politi), the Ginny and Kenneth Grunley Fund for Lung Cancer Research (K. Politi), and the Yale Cancer Center Class of ’61 Research Award (K. Politi). Gilead Sciences, Inc., supported the sequencing of a subset of the Yale Cancer Center specimens. Osimertinib is an agent developed by AstraZeneca.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.