Tumor heterogeneity underlies resistance to tyrosine kinase inhibitors (TKI) in lung cancers harboring EGFR mutations. Previous evidence suggested that subsets of preexisting resistant cells are selected by EGFR-TKI treatment, or alternatively, that diverse acquired resistance mechanisms emerge from drug-tolerant persister (DTP) cells. Many studies have used bulk tumor specimens or subcloned resistant cell lines to identify resistance mechanism. However, intratumoral heterogeneity can result in divergent responses to therapies, requiring additional approaches to reveal the complete spectrum of resistance mechanisms. Using EGFR-TKI-resistant cell models and clinical specimens, we performed single-cell RNA-seq and single-cell ATAC-seq analyses to define the transcriptional and epigenetic landscape of parental cells, DTPs, and tumor cells in a fully resistant state. In addition to AURKA, VIM, and AXL, which are all known to induce EGFR-TKI resistance, CD74 was identified as a novel gene that plays a critical role in the drug-tolerant state. In vitro and in vivo experiments demonstrated that CD74 upregulation confers resistance to the EGFR-TKI osimertinib and blocks apoptosis, enabling tumor regrowth. Overall, this study provides new insight into the mechanisms underlying resistance to EGFR-TKIs.

Significance:

Single-cell analyses identify diverse mechanisms of resistance as well as the state of tolerant cells that give rise to resistance to EGFR tyrosine kinase inhibitors.

Recent advances in genomic sequencing have allowed identification of somatic mutations or translocations in receptor tyrosine kinases in lung cancer cells. These genetic alterations result in aberrant activation of tyrosine kinases, eventually leading to development of cancer. The most common somatic mutations, which are activating EGFR mutations [such as exon 19 deletions (Del19) or the exon 21 L858R mutation], are detectable in 10% to 30% of non–small cell lung cancers (NSCLC) and promote downstream prosurvival/antiapoptotic signals (1, 2). Dysregulation of tyrosine kinase activity by EGFR gene mutations results in “oncogene addiction,” in which cells become dependent on these mutations for survival, proliferation, invasion, and metastasis (3). Although many patients with EGFR mutations show dramatic responses to first-generation (gefitinib and erlotinib) and second-generation (afatinib and dacomitinib) EGFR-TKIs, most acquire resistance within 2 years of continued drug exposure (4). Previous studies show that resistance to first- and second-generation EGFR-TKIs occurs due to emergence of the EGFR-T790M gatekeeper mutation (5, 6), activation of alternative growth factor receptor signaling pathways such as MET (7, 8), ERBB2 amplification (9), activation of EGFR signaling due to activating mutations in downstream genes such as PIK3CA, and histologic transformation (10, 11). The third-generation EGFR TKI osimertinib was developed to target cells harboring T790M and irreversibly binds to the ATP binding site at EGFR C797 (12). Preclinical studies suggest that osimertinib inhibits not only T790M clones but also classic EGFR-sensitizing mutations, as evidenced by a phase III trial comparing osimertinib with first-generation EGFR TKIs as first-line therapy (the FLAURA trial). That trial showed that both median PFS and overall survival (OS) were significantly longer following treatment with osimertinib than with first-generation EGFR-TKIs (13, 14). However, 20% to 30% of patients are intrinsically resistant osimertinib and either do not respond or have very short PFS (less than 6 months; ref. 15). Even for patients who initially respond, acquired resistance emerges within approximately 24 months. Analysis of the AURA cohort further revealed that emergence of resistance mutations, such as C797S, or loss of the T790M mutation was not detected in 45% (10/22) of patients who developed early resistance (15). Interestingly, analysis of the FLAURA cohort indicated that resistance to first-line osimertinib treatment differs from that of second-line osimertinib in the AURA cohort (16). Mechanisms underlying this difference are unknown.

Recently, several lines of evidence suggest that drug-tolerant persister (DTP) populations of cancer cells, which survive and rapidly adapt to therapy, play an important role in emergence of resistance to targeted therapies (17–19). Such DTP cells employ diverse genetically based resistance mechanisms and contribute to spatial and temporal heterogeneity of tumor cells. Therefore, targeting DTP cells may be an efficient strategy to prevent resistance to EGFR TKIs.

Most reported resistance mechanisms have been identified by analysis of bulk tumor specimens and/or subcloned resistant cell lines. Given that intratumoral heterogeneity is a driver of drug resistance, this approach may not adequately reveal the complete set of resistance mechanisms (20). To overcome this problem, single-cell analysis has recently been used to characterize individual cell subpopulations within a bulk tumor. For example, we previously used single-cell RNA-seq analysis to show a heterogeneous transcriptome response following gefitinib treatment of EGFR-mutant cells and identified a small subpopulation of cells that expressed relatively high levels of AURKA (21), which encodes a kinase that renders EGFR-mutant cells resistant to EGFR TKIs (22).

In this study, we analyzed intratumoral heterogeneity of subpopulations of lung cancer cells sensitive to EGFR TKIs and identified key gene expression changes that occur as osimertinib resistance is acquired. We also sought to identify novel mechanisms by which DTPs emerge. Our results demonstrate transcriptomic and epigenomic responses to EGFR-TKIs at the single-cell level.

Cells

NCI-H1975 (H1975) cells, Phoenix-Ampho cells, and HEK293T cells were purchased from the ATCC, and PC9 cells were kindly provided by Dr. Pasi Jänne. All cell lines were authenticated by short-tandem repeat DNA profiling and were routinely tested for Mycoplasma using the Mycoalert Mycoplasma Detection Kit (Lonza). PC9-ER cells were established as described (23). To obtain resistant lines, we cultured each cell line with increasing concentrations of osimertinib (LC laboratories, #O7200; 10 nmol/L to 2,000 nmol/L). All lines were maintained in RPMI1640 (Sigma Aldrich, #R8758-500ML) plus heat-inactivated 10% FBS (Sigma, #F9423) and antibiotics-antimycotics (Gibco, #15240062). Cells were cultured in 5% CO2 at 37°C.

Western blotting

Cell cytoplasmic extracts were prepared using protein lysis buffer plus protease and phosphatase inhibitors (Sigma Aldrich, #8340). Whole cell extracts were prepared as follows: cells were rinsed once with ice-cold 0.9% NaCl, and protein was precipitated with 100% trichloroacetic acid. Precipitates were lysed in cell lysis buffer (0.05 M Tris pH 6.8, 3% SDS, 15% sucrose, 0.01% bromophenol blue, and 10% 2-mercaptoethanol) and incubated at 95% for 3 minutes. SDS-PAGE was performed and transferred to PVDF membranes (catalog no. 456-1094). Membranes were blocked in 5% BSA in TBS-T buffer. Anti-CD74 (#77274), anti-β-actin (#4970), and anti-BCL-XL (#2764) antibodies were purchased from Cell Signaling Technology (Danvers). After incubation with HRP-conjugated secondary antibodies, we detected signals using ECL Prime Western blotting detection reagent (GE Healthcare). Both primary and secondary antibodies were diluted 1:1,000. Images were acquired and recorded using ImageQuant LAS4000 (GE Healthcare).

siRNA studies

ON-TARGETplus SMARTpool siRNAs were obtained from Thermo Fisher Scientific and resuspended in 1× Dharmacon siRNA buffer prior to use. Indicated cell lines were transfected with CD74 siRNA (Thermo Fisher Scientific, #L-012667-00-0005) or control siRNA (Qiagen, #102207), according to the manufacturer's instructions. Knockdown efficacy was validated by Western blotting.

Cell viability assays

H1975 and PC9 cells were plated into 96-well plates (2,000 cells per well) and treated with osimertinib, erlotinib (Santa Cruz Biotechnology, #sc202154a), alisertib (MLN8237, #1133, Selleckchem), or DMSO, as indicated. Cell viability was assessed using CellTiter 96 Aqueous One Solution Cell Proliferation Assay (Promega, #G3581). Combination Index or IC50 values were calculated using Excel software (Microsoft). Experiments were independently repeated at least three times.

Clinical specimens

All patients provided written informed consent before sampling, in accordance with the Declaration of Helsinki. This study was performed in a blinded manner and was approved by the National Cancer Center Ethics Committee.

Single-cell RNA-seq analysis

To obtain single-cell RNA-seq data, we used 10× genomics chromium single-cell 3′ v2 and a V(D)J Kit (10x Genomics) according to the manufacturer's instructions. Briefly, cultured cells were diluted and mixed with single-cell master mix, gel beads, and partition oil into the bottom of the wells on a chip. Libraries were sequenced at 150 paired-end cycles using the HiSeq3000 system (Illumina). Sequencing bcl2 files were converted into FASTQ files using the Cell Ranger pipeline provided by 10x Genomics. Sequence was mapped to the human genome (UCSC hg38) using STAR, which was included in Cell Ranger. Using Seurat ver2.3.4 (24), low-quality reads and PCR “sister” duplicates were removed. To filter out low-quality cell data, the following threshold was determined for each sample. For cell lines, (i) >7,500 UMI per cell, (ii) >1,000 genes per cell, and (iii) <10% mitochondrial gene expression; for Pt-1 and -3; (i) >5,000 UMI per cell, (ii) >1,000 genes per cell, and (iii) <20% mitochondrial gene expression; for Pt-2; (i) >5,000 UMI per cell, (ii) >1,000 genes per cell, and (iii) <10% mitochondrial gene expression. Stats of the scRNA-seq analysis are shown in Supplementary Table S1.

Single-cell ATAC sequencing analysis

To obtain single-cell assay for transposase-accessible chromatin using sequencing (scATAC-seq) data, we used 10x Genomics Chromium single-cell ATAC Reagent Kit (10x Genomics) according to the manufacturer's instructions. We set the number of target cells as 3,000. We processed scATAC-seq data using Cell Ranger and Signac (ver1.0.0). bcl2 files were converted into FASTQ files using the Cell Ranger pipeline provided by 10x Genomics. Sequence was mapped to the human genome (UCSC hg38) using STAR, which was included in Cell Ranger. Using Signac, low-quality reads were removed. Stats of the scATAC-seq analysis are shown in Supplementary Table S2.

Target sequence using the NCC oncopanel

A total of 1.0 × 106 H1975 and PC9 cells for each phase (a total of seven samples) were used to extract genomic DNA using the DNeasy Blood & Tissue Kit (Qiagen). Genomic DNA was processed, and sequencing libraries were subjected to the SureSelect NCC Oncopanel (Agilent).

CD74 expression constructs

pcDNA3.1+/C-(K)DYK-CD74, which encodes human CD74, was purchased from GenScript (#OHu19995). CD74 coding sequence was PCR-amplified using Pfu Turbo DNA polymerase (Agilent). PCR products were verified by sequencing and then ligated into the MIGR1 retroviral vector (a gift of Dr. Warren Pear; Addgene plasmid #27490; http://n2t.net/addgene:27490; RRID:Addgene_27490) between the XhoI and HpaI sites.

Stable cell lines

Phoenix-Ampho cells were transfected 72 hours with MIGR1 vector encoding CD74 or with empty vector using the TransIT-X2 Dynamics Delivery System (Mirus). Culture supernatants were collected thereafter at 24, 48, and 72 hours through a 0.45-μm filter, incubated overnight at 4°C with a 25% poly ethylene glycol solution (Sigma-Aldrich, #81260), and centrifuged at 1,500 × g for 30 minutes at 4°C to pellet retrovirus. Target cells were infected with retrovirus in the presence of 8 μg/mL polybrene (Thermo Fisher Scientific, NC9200896) and then subjected to GFP sorting to obtain either control or CD74-overexpressing cells.

Generation of CD74 knockout cell lines

CD74 knockout was generated in H1975 cells using the lentiCRISPR V2 system. Briefly, HEK293T cells were cotransfected with lentiviral vectors encoding guide RNA targeting human CD74 (Genscript; ATGCAGAATGCCACCAAGTA) or lentiCRISPRv2 backbone vectors (Addgene, #52961) and pCMV-dR8.91 (a gift of Dr. Bob Weinberg; Addgene plasmid #8455; http://n2t.net/addgene:8455; RRID:Addgene_8455) and pMD2.G (a gift of Dr. Didier Trono; Addgene plasmid #12259; http://n2t.net/addgene:12259; RRID:Addgene_12259) using the TransIT-X2 Dynamics Delivery System. H1975 cells were then infected as described above, followed by selection in 3 μg/mL puromycin (Thermo Fisher Scientific, #AAJ672368EQ) to obtain targeted gene knockout cells. Cells were diluted and grown until visible colonies could be individually selected for Western blot analysis to confirm loss of CD74 expression. Each clone was verified by sequencing.

qPCR

Total RNA was isolated from cells using an RNA Isolation Kit (Qiagen, RNeasy Mini Kit) and subjected to reverse transcription (RT) using Invitrogen Superscript II Reverse Transcriptase (Thermo Fisher Scientific). mRNA abundance was normalized to that of GAPDH. Real-time PCR analysis was performed in triplicate with using iTaq Universal SYBR Green Supermix (Bio-Rad). Respective forward and reverse PCR primers were: CD74 (5′-GGAGAAGCAGGAGCTGTCGG-3′ and 5′-CCAGCATGGGCAGTTGCTCA-3′) and GAPDH (5′-CCACATCGCTCAGACACCAT-3′ and 5′-CCAGGCGCCCAATACG-3′).

Flow cytometric analysis

Cells were stained with PE-labeled mouse monoclonal antibodies recognizing human CD74 (BioLegend, #326808) or PE-labeled mouse IgG1, κ immunoglobulin isotype control (BioLegend, #400114) for flow cytometric analysis using a CytoFLEX LX Flow Cytometer (Beckman Coulter).

Soft agar colony formation assay

A total of 1.5 × 104 cells in 1 mL RPMI containing 0.3% agar was plated in each well of 6-well plates coated with 2 mL RPMI containing 0.5% agar. RPMI medium (200 μL) was added as supernatant twice weekly. Supernatants contained either 30 nmol/L Osimertinib, 50 ng/mL macrophage migration inhibitory factor (MIF) human recombinant (MIF; ProSpec, #CYT-575), or DMSO vehicle. After 4 to 5 weeks, colonies were stained with 0.005% crystal violet solution and quantified using ImageJ software.

Apoptosis assay

Apoptosis was assessed by measuring caspase activity using the Caspase-Glo 3/7 Assay System (Promega, #G8091). Cells (2,000 cells/well) were plated in 96-well plates, incubated for 24 hours and then treated 24 more hours with osimertinib (100 nmol/L) or DMSO vehicle. Caspase activity was measured according to the manufacturer's instructions. Results are presented as luminescence units obtained after subtracting the luminescence value of a blank reaction and dividing by cell number.

Transient transfection

H1975-KO and control cells were plated at 70 to 80% confluence in 6-well plates or at 8,000 cells/well in 96-well plates and incubated for 24 hours. Cells were then transiently transfected with pcDNA3.1-empty or pcDNA3.1-CD74 vectors using a TransIT-X2 Dynamics Delivery System. After overnight incubation, the medium was changed and osimertinib (100 nmol/L), MIF (50 ng/mL), or DMSO was added. Cells in 6-well plates were collected for Western blot analysis, and cells in 96-well plates were assayed for caspase activity.

Murine xenograft models

All animal studies were approved by the Institutional Animal Care and Use Committee at Beth Israel Deaconess Medical Center. Female Foxn1 (nu/nu) mice (JAX stock, #002019) were obtained from The Jackson Laboratory. To establish xenografts, H1975-OE, -KO, or control cells (3–10 × 106) were subcutaneously injected. After an optimal tumor volume was reached, mice were randomized into groups (5–6 mice per group) and administered osimertinib (5 mg/kg/day), or vehicle by oral gavage once daily during the treatment period (11 days for KO cells and 28 days for OE cells) followed by a 14-day withdrawal period. Tumor size was measured three times a week, with caliper and body weight once a week. Tumor volume was calculated according to the following formula; (tumor volume) = (major axis) × (minor axis)2 × 0.5. Tumor volume change was calculated as follows: [(tumor volume) – (baseline volume)]/maximum volume of vehicle group. Data are presented as the mean % change in tumor volume ± SD. Asterisks in figures indicate P value as follows: *, P < 0.05 and **, P < 0.005.

Generation of cell line models representing resistance to first- or second-line osimertinib

To investigate cellular responses that occur during acquisition of resistance to osimertinib, we generated polyclonal osimertinib-resistant cells using three EGFR-mutant osimertinib-sensitive cell lines: H1975, PC9, and PC9-ER cells. H1975 cells harbor L858R (exon 21) + T790M (exon 20) mutations and are resistant to first- and second-generation EGFR TKIs due to the T790M mutation. PC9 cells harbor the EGFR Del19 (E746-A750) and are sensitive to all generations of EGFR TKIs. Erlotinib-resistant PC9-ER cells were established as described (23) and harbor the EGFR exon 19 deletion and the T790M mutation. Thus, H1975 and PC9 cells resistant to osimertinib (H1975-OR and PC9-OR, respectively) represent EGFR-mutant tumor cells treated with first-line osimertinib, while PC9-ER cells rendered resistant to osimertinib (PC9-EROR) represent those treated with second-line osimertinib after failure of first- or second-generation TKIs (Fig. 1A). To create OR lines, cells were exposed to increasing concentrations of osimertinib up to 2 μmol/L (Supplementary Fig. S1A) and then stored when their growth reached a steady-state at various concentrations. In all lines we performed cell viability assays to confirm osimertinib resistance (Supplementary Fig. S1B). Paclitaxel sensitivity remained similar in OR cells cultured in the absence of osimertinib for at least 2 weeks compared with corresponding parental cells, suggesting that multidrug resistance mechanisms do not underlie osimertinib resistance in OR cells (Supplementary Fig. S1C; ref. 25). PC9-OR and PC9-EROR cells were morphologically similar to corresponding parental cells. However, relative to parental H1975 cells, H1975-OR cells were larger and grew in discrete groups or patches of cells rather than in an even distribution (Supplementary Fig. S1D).

Figure 1.

scRNA-seq analysis of osimertinib-resistant lung cancer line models. A, Schematic showing establishment of osimertinib-resistant lines used in this study. B, t-Distribution stochastic neighbor embedding (tSNE) plot based on scRNA-seq in merged H975 parental and H975-OR2,000 (left) cells and four PC9 datasets (right; PC9 parental, PC9-ER parental, PC9-OR2,000, and PC9-EROR2,000). C–E, t-SNE plots based on scRNA-seq of merged datasets for each line: H1975 (C), PC9-ER (D), and PC9 (E) datasets. F–H, Cell ordering-based pseudotime analyzed by Monocle for each line: H1975 (F), PC9-ER (G), and PC9 (H) datasets. (I–K) Heatmaps showing differentially expressed genes for each cluster in H975 (I), PC9-ER (J), and PC9 (K) datasets. Normalized gene expression is shown.

Figure 1.

scRNA-seq analysis of osimertinib-resistant lung cancer line models. A, Schematic showing establishment of osimertinib-resistant lines used in this study. B, t-Distribution stochastic neighbor embedding (tSNE) plot based on scRNA-seq in merged H975 parental and H975-OR2,000 (left) cells and four PC9 datasets (right; PC9 parental, PC9-ER parental, PC9-OR2,000, and PC9-EROR2,000). C–E, t-SNE plots based on scRNA-seq of merged datasets for each line: H1975 (C), PC9-ER (D), and PC9 (E) datasets. F–H, Cell ordering-based pseudotime analyzed by Monocle for each line: H1975 (F), PC9-ER (G), and PC9 (H) datasets. (I–K) Heatmaps showing differentially expressed genes for each cluster in H975 (I), PC9-ER (J), and PC9 (K) datasets. Normalized gene expression is shown.

Close modal

To determine whether osimertinib treatment induced acquired mutations that underlie potential changes in gene expression, we performed deep sequencing at a mean depth of approximately 3000x. As expected, H1975 parental and -OR cells exhibited L858R and T790M mutations, and PC9 parental and PC9-ER cells showed Del19. PC9-ER cells also carried the T790M mutation, which was not detected in PC9-parental or PC9-OR cells (Supplementary Fig. S1E). No additional mutations were detected in EGFR or relevant genes such as PIK3CA or RAS, which reportedly confer EGFR-TKI resistance. We also checked T790M allele frequency in PC9-ER cells at each osimertinib concentration. Mutation frequency was stable at 20%, suggesting that T790M loss does not underlie acquisition of resistance to osimertinib in PC9-EROR cells (Supplementary Fig. S1F).

Single cell RNA-seq analysis of osimertinib- and/or erlotinib-resistant cells

To detect transcriptome changes occurring during osimertinib exposure at the single-cell level, we subjected parental cells and cells resistant to 2,000 nmol/L (OR2,000) osimertinib to scRNA-seq analysis. After noise reduction procedures (see Materials and Methods for details), we obtained scRNA-seq data from a total of 7,464, 7,757, and 7,787 cells in the H1975-OR, PC9-OR, and PC9-EROR series, respectively (Supplementary Table S1). Among these cells, approximately 3,000 were analyzed for each cell line. t-SNE plot analysis of four end-point data clusters showed that each group of parental or OR cells exhibited a distinct gene expression profile (Fig 1B). Interestingly, expression patterns seen in PC9-OR and PC9-EROR cells were distinct, although both cells showed resistance to osimertinib (Supplementary Fig. S1B).

Tumor cells reportedly maintain a transient DTP state following treatment with targeted therapies, allowing cell subpopulations to adapt to drugs by acquiring genetic or epigenetic alterations (18, 26). We hypothesized that we would detect transcriptomic or epigenetic changes in cells in this tolerant state, given that it is a reversible phenotype characterized by reduced drug sensitivity (18). Thus, in addition to parental and OR2,000 cells, we examined transcriptomic changes in cultured in 30 nmol/L osimertinib (OR30 cells) as a model of the drug-tolerant state. Osimertinib IC50 of were approximately 30 nmol/L in all OR30 cells (Supplementary Fig. S1B). t-SNE plots using merged datasets indicated clusters of H1975-OR30 (Fig. 1C) and PC9-EROR30 (Fig. 1D) cells, partially overlapping with parental and PC9-EROR2,000 clusters, respectively, whereas three clusters of parental, OR30, and OR2,000 cells were independent in PC9-OR cells (Fig. 1E). Accordingly, pseudotime analysis of cellular trajectory revealed that H1975-OR30 and PC9-EROR30 were ordered closely with parental H1975 and PC9-EROR2,000, respectively, whereas each osimertinib-treated cluster of PC9 cells showed distinct transcriptome patterns (Fig. 1F,H). Hierarchal clustering analysis showed that H1975, PC9-ER, and PC9 cells comprised 7, 8, and 11 clusters, respectively (Fig. 1I,K). Cellular heterogeneity was observed at every stage in all cell lines. Taken together, these results suggest that osimertinib exposure induces significant and complex transcriptome changes as cells acquire resistance and that that process is distinct among modeled cell lines.

AURKA and TPX2 upregulation in response to first-line osimertinib treatment is a mechanism of intrinsic and acquired resistance

Next, we asked whether our scRNA-seq strategy could identify genes mediating osimertinib-resistance. aurora kinase A is a serine/threonine kinase functioning in microtubule formation and/or stabilization at the spindle pole during chromosome segregation, and its activation by its coactivator TPX2 promotes osimertinib resistance (22). Violin plots (Fig. 2A,C) showed that AURKA expression was downregulated by osimertinib treatment in PO9-EROR cells (Fig. 2B), suggesting that AURKA does not participate in mechanisms of resistance to “second-line” osimertinib in this context. In contrast, we detected relatively high AURKA expression in several cell clusters in tolerant (OR30) and resistant (OR2,000) H1975-OR cells (Fig. 2A), whereas AURKA expression was upregulated in the tolerant state and remained high in a cluster in highly resistant (OR2,000) PC9-OR cells (Fig. 2C). Interestingly, TPX2 expression overlapped with that of AURKA in all cases tested (Fig. 2D,F). AURKA expression was limited to several clusters in all states of both H1975 and PC9 cell types (Fig. 2G and H). To confirm Aurora kinase A function in osimertinib resistance, we treated the H1975 series with a combination of osimertinib and alisertib, an Aurora kinase A specific inhibitor. Alisertinib cotreatment enhanced cell growth suppression induced by osimertinib in all states of H1975 cells (Fig. 2I,K). These results suggest that cells with higher AURKA and TPX2 expression show intrinsic (H1975) or acquired (PC9-OR) resistance to first-line osimertinib.

Figure 2.

Identification of Aurora kinase A as functioning in intrinsic and/or acquired resistance to first-line osimertinib. A–C, t-SNE plots (left) and violin plots (right) showing AURKA expression in H1975 (A), PC9-ER (B), and PC9 (C) datasets. D–F, t-SNE plots (left) and violin plots (right) showing TPX2 expression in H1975 (D), PC9-ER (E), and PC9 (F) datasets. G and H, t-SNE plots highlighting clusters showing high AURKA expression in H1975 (G) and PC9 (H) datasets. I–K, Cell viability assay using H1975, H1975-OR30, and H1975-OR2,000 cells treated with osimertinib alone or in combination with alisertib. *, P < 0.05.

Figure 2.

Identification of Aurora kinase A as functioning in intrinsic and/or acquired resistance to first-line osimertinib. A–C, t-SNE plots (left) and violin plots (right) showing AURKA expression in H1975 (A), PC9-ER (B), and PC9 (C) datasets. D–F, t-SNE plots (left) and violin plots (right) showing TPX2 expression in H1975 (D), PC9-ER (E), and PC9 (F) datasets. G and H, t-SNE plots highlighting clusters showing high AURKA expression in H1975 (G) and PC9 (H) datasets. I–K, Cell viability assay using H1975, H1975-OR30, and H1975-OR2,000 cells treated with osimertinib alone or in combination with alisertib. *, P < 0.05.

Close modal

scRNA-seq reveals that the epithelial–mesenchymal transition is associated with the DTP state

The epithelial–mesenchymal transition (EMT) and expression of related genes reportedly play a pivotal role in acquisition of drug resistance in numerous cancer types (27). The EMT is also reportedly involved in an EGFR-TKI-induced drug-tolerant state, and some mesenchymal phenotypes are maintained after acquisition of additional mutations (17, 28). Violin plots of the EMT-related gene vimentin (VIM) in H1975 (Fig. 3A), PC9-ER (Fig. 3B), and PC9 cells (Fig. 3C) showed that its expression increased in OR30 cells but returned to basal levels in OR2,000 cells in H1975 and PC9 cells (Fig. 3A and C). In contrast, it remained high in PC9-EROR2,000 cells (Fig. 3B). CDH1 expression was also examined in H1975 (Supplementary Fig. S2A), PC9-ER (Supplementary Fig. S2B), and PC9 cells (Supplementary Fig. S2C). Of note, it decreased inversely with VIM expression in some H1975 and PC9 cell clusters (Supplementary Figs. S2A and S2C: red arrows), supporting the idea that osimertinib induces the EMT in some, if not all, cell populations. In contrast, high VIM expression and relatively low CDH1 expression were observed at all states of PC9-EROR cells (Supplementary Fig. S2B), indicating that PC9-ER cells acquire EMT characteristics and maintain them during osimertinib treatment, consistent with a previous report (17).

Figure 3.

The EMT is associated with a DTP state. A–C, t-SNE plots (left) and violin plots (right) showing VIM expression in H1975 (A), PC9-ER (B), and PC9 (C) datasets. D–F, t-SNE plots (left) and violin plots (right) showing AXL expression in H1975 (D), PC9-ER (E), and PC9 (F) datasets. G, scATAC-seq data showing peaks in genomic regions of VIM (left) and AXL (right) in indicated parental and H1975-OR series. Bottom, relevant ChIP-Seq data based on ENCODE3. H, scATAC-seq data showing peaks in genomic regions of VIM (left) and AXL (right) in 10 clusters from H1975 datasets. I, t-SNE diagrams based on scATAC-seq embedding scRNA-seq (left) and scRNA-seq (right) in H1975 datasets. J, Violin plots showing indicated gene expression in each H1975-OR30 cluster. *, P < 0.05.

Figure 3.

The EMT is associated with a DTP state. A–C, t-SNE plots (left) and violin plots (right) showing VIM expression in H1975 (A), PC9-ER (B), and PC9 (C) datasets. D–F, t-SNE plots (left) and violin plots (right) showing AXL expression in H1975 (D), PC9-ER (E), and PC9 (F) datasets. G, scATAC-seq data showing peaks in genomic regions of VIM (left) and AXL (right) in indicated parental and H1975-OR series. Bottom, relevant ChIP-Seq data based on ENCODE3. H, scATAC-seq data showing peaks in genomic regions of VIM (left) and AXL (right) in 10 clusters from H1975 datasets. I, t-SNE diagrams based on scATAC-seq embedding scRNA-seq (left) and scRNA-seq (right) in H1975 datasets. J, Violin plots showing indicated gene expression in each H1975-OR30 cluster. *, P < 0.05.

Close modal

Previous reports suggest that the tyrosine kinase receptor AXL is upregulated in tumors resistant to EGFR-TKIs, including osimertinib, and that AXL functions in intrinsic osimertinib resistance (29, 30). Although AXL was expressed in a pattern similar to VIM in all cell datasets (Fig. 3D,F), AXL and VIM expression did not overlap in all cell lines (Supplementary Figs. S2DS2F), suggesting that AXL is not the sole mechanism underlying EMT induction following EGFR TKI treatment.

To assess a potential epigenetic basis for underlying transcriptomic changes associated with altered VIM and AXL expression, we performed a scATAC-seq analysis (31) using the H1975 cell series. Overall, we generated a total of 10,904 cells for scATAC datasets: 1,315, 5,089 and 4,500 cells for H1975-parental (Supplementary Fig. S3A), H1975-OR30 (Supplementary Fig. S3B) and H1975-OR2,000 (Supplementary Fig. S3C), respectively. Parental (Supplementary Fig. S3A) and OR2,000 (Supplementary Fig. S3C) cells showed no clear separation of individual cells, forming relatively loose overall cluster structures. However, scATAC-seq showed that H1975-OR30 cells formed a series of tight clusters (Supplementary Fig S3B). Using the H1975 dataset, we asked whether changes in VIM expression were controlled at the epigenetic level. Consistent with scRNA-seq results, the VIM promoter region showed an open chromatin structure in H1975-OR30 cells as indicated by aggregated read analysis (Fig. 3G). When separated into different clusters, the “peak” of ATAC-seq tags was most significant for cluster 7 in the OR30 state (indicated by red arrow in Fig. 3H). On the basis of ENCODE3 ChIP-seq data, common peaks in other clusters are predicted to be ZEB1 and ZEB2 binding sites, and the specific peak in OR30 is predicted to be the CREB1 binding site (Fig. 3G). All of these transcription factors are known to function in the EMT (32). To further characterize cluster 7 cells, we associated scATAC-seq data with scRNA-seq data using Signac. Cluster 7 was associated with a particular scRNA-seq cluster in OR30 cells (cluster B; Fig. 3I). That cluster showed the highest VIM expression levels (Fig. 3J).

CD74 is induced by osimertinib at a drug-tolerant state, based on scRNA-seq

Hierarchal clustering analysis indicated several genes upregulated in OR30 cell clusters that may underlie the transition from parental cells to a DTP state. Analysis of differentially expressed genes (DEG) in the H1975 data set showed that CD74 was induced in some H1975-OR30 cell clusters and its expression remained only in a small fraction of H1975-OR2,000 cells (Fig. 4A and B). In PC9-EROR cells, CD74 was expressed in parental cells and that expression increased in PC9-EROR30 cells and remained in some clusters of PC9-EROR2,000 cells (Fig. 4B). We detected no overlap between VIM-high and CD74-high clusters (Supplementary Fig. S4A), suggesting that CD74 upregulation occurs independently of the EMT. scATAC-seq revealed a specific peak of open chromatin only in H1975-OR30 cells (Fig. 4C). Of note, we observed other peaks at a neighboring 5′-proximal intronic region, which may represent a distinct enhancer region (Fig. 4C). These peaks remained unchanged during state progression, possibly providing a background chromatin context for induction of gene expression. We observed similar patterns in PC9-EROR but not PC9-OR cells (Supplementary Fig. S4B). We next used ENCODE ChIP-seq datasets to identify transcription factors that bind to promoter and enhancer regions identified above. As expected, we observed potential binding of RFX5, CREB1, and NFYB, all of which form a complex with CIITA, a positive regulator of class II major histocompatibility complex gene transcription (Fig. 4C).

Figure 4.

scRNA-seq identifies CD74 as a novel gene expressed in DTPs. A, t-SNE diagram showing clusters in merged H1975 datasets (left) and clusters that show higher CD74 expression in indicated parental and H1975 H1975-OR lines (right). B, t-SNE plots (top) and violin plots (bottom) showing CD74 expression in H1975 datasets. C, Top, scATAC-seq data showing peaks in the CD74 genomic region in indicated H1975 datasets. Bottom, RFX5, CREB1, and NFYB binding sites in indicated lines based on the ENCODE database. D, Gene coexpression networks involving CD74 based on our previous study (36). E, Transitions in gene expression of CD74-related genes: genes upstream of CD74 genes (top and middle rows) and downstream genes (bottom). F, Predicted signaling events downstream of CD74 following osimertinib treatment.

Figure 4.

scRNA-seq identifies CD74 as a novel gene expressed in DTPs. A, t-SNE diagram showing clusters in merged H1975 datasets (left) and clusters that show higher CD74 expression in indicated parental and H1975 H1975-OR lines (right). B, t-SNE plots (top) and violin plots (bottom) showing CD74 expression in H1975 datasets. C, Top, scATAC-seq data showing peaks in the CD74 genomic region in indicated H1975 datasets. Bottom, RFX5, CREB1, and NFYB binding sites in indicated lines based on the ENCODE database. D, Gene coexpression networks involving CD74 based on our previous study (36). E, Transitions in gene expression of CD74-related genes: genes upstream of CD74 genes (top and middle rows) and downstream genes (bottom). F, Predicted signaling events downstream of CD74 following osimertinib treatment.

Close modal

The transmembrane protein CD74 plays an essential role in MHCII antigen presentation and regulates endosomal trafficking, cell migration, and cellular signaling as a surface receptor for macrophage MIF (33, 34). Upon stimulation, the CD74 intracellular domain is cut by gamma secretase, processed, and eventually becomes nuclear and serves as a transcriptional regulator by interacting with factors including RUNX and NF-κB (35). Previous studies show that CD74 is downstream of CIITA and regulates proinflammatory reactions stimulated by cytokines, although its role in lung cancer cells is unknown. To assess CD74 function in osimertinib resistance in lung cancer, we referred to a catalog of gene coexpression networks that we constructed from a panel of lung cancer lines in our previous study (36). CD74 was present in one of these networks, which included a total of 135 genes (Supplementary Fig. S4C). Those genes comprised several categories, including STAT1/IRF-associated interferon signaling, MHC class I/II-associated antigen processing and presentation, and TNF-related apoptotic signaling (Fig. 4D). Thus, we hypothesized that STAT1/IRF-associated defense responses or cell survival mechanisms may be upregulated in the drug-tolerant state as STAT1/IRF signaling activates CIITA and several TFs that regulate HLA class II-associated genes, leading to upregulation of HLA genes and CD74. Analysis of H1975-OR scRNA-seq datasets detected expression of STAT1, CIITA, RFX5, and IRF1 in parental and OR30 cells (Fig. 4E, top) and upregulation of the downstream genes HLA-B and CD74 in H1975-OR30 cells (Fig. 4E, left middle,). Furthermore, cytokine genes (e.g., IL6 and IL11) and genes associated with resistance to apoptosis (e.g., BCL2L1, CREB3, TRAF1, and TRAF2), all downstream of RUNX/NF-kB signaling, were upregulated simultaneously or after CD74 upregulation (Fig. 4E). Importantly, BCL2L1 encodes BCL-XL, an anti-apoptotic BCL-2 family member that plays an important role in emergence and maintenance of DTP cells (17). Taken together, we hypothesized that osimertinib-induced CD74 would lead to inhibition of apoptosis (Fig. 4F).

CD74 upregulation suppresses apoptosis and contributes to osimertinib resistance

Next, we used H1975 cells to confirm scRNA-seq data showing CD74 upregulation by osimertinib (Fig. 4A and B). Osimertinib treatment of H1975 cells induced CD74 mRNA (Supplementary Fig. S5A) and protein (Supplementary Fig. S5B) expression. Moreover, flow cytometry analysis indicated upregulation of CD74 expression on the surface of osimertinib-treated H1975 cells (Supplementary Fig. S5C). To determine whether these changes contribute to osimertinib resistance, we used siRNA to knock down CD74 expression to non-detectable levels in osimertinib-resistant H1975-OR30 cells (Supplementary Fig. S5D). We then treated these and siControl H1975-OR30 cells with increasing concentrations of osimertinib and assayed their viability using an MTS assay. Relative to corresponding siControl cells, CD74 knockdown H1975-OR30 cells showed decreased viability following osimertinib treatment (Supplementary Fig. S5E).

To further assess CD74 function in osimertinib resistance, we established H1975 cells that overexpressed CD74 (OE) as well as H1975 cells deficient in CD74 (KO). Immunoblot analysis confirmed CD74 overexpression in OE cells (Supplementary Fig. S5F, left), and showed that CD74 was undetectable in two different KO clones (KO#1 and KO#2; Supplementary Fig. S5F, right). As indicated in our scRNA-seq dataset (Fig. 4E), osimertinib treatment upregulated BCL-XL expression in H1975 control cells (Fig. 5A; lanes 1 vs. 2, 3), an effect enhanced in H1975-OE cells (Fig. 5A; lanes 1, 2, 3 vs. 4, 5, 6). In contrast, we observed decreased expression of BCL-XL protein in H1975-KO cells after osimertinib treatment, particularly at the 48-hour time point (Fig. 5B; lanes 3 vs. 6, 9). We next asked whether CD74 expression suppresses apoptosis induced by osimertinib. To do so, we treated H1975-OE, H1975-KO, or control cells with osimertinib and assayed caspase-3/7 activity as an apoptotic marker. Relative to control cells, H1975-OE cells showed significantly decreased caspase-3/7 activity following osimertinib treatment (Fig. 5C). In contrast, caspase-3/7 activity was significantly increased in osimertinib-treated H1975-KO cells compared with control cells (Fig. 5D). To confirm that effects were CD74-specific, we repeated this analysis using H1975 KO cells transiently transfected with a CD74 overexpression construct. Upregulation of caspase-3/7 activity in response to osimertinib treatment seen in untransfected cells was attenuated in H1975 KO cells transiently overexpressing CD74 cells (Fig. 5E; Supplementary Fig. S5G).

Figure 5.

CD74 upregulation by osimertinib contributes to a drug-tolerant state in EGFR-mutant tumor cells. A, Western blotting showing upregulation of CD74 and BCL-XL following osimertinib treatment of H1975 cells overexpressing CD74 (OE) or control H1975 cells (Cntl). Both groups were incubated for 0, 24, and 48 hours with osimertinib at 100 nmol/L. Long exposure increased CD74 expression in control cells (Cntl). B, Western blotting showing CD74 knockout efficiency and BCL-XL downregulation in H1975 CD74 KO cells. Two different KO clones (KO#1 and KO#2) and control cells were treated with osimertinib at 100 nmol/L for 0, 24, and 48 hours in the presence of macrophage MIF (50 ng/mL). C and D, Caspase-3/7 activity indicative of inhibition of apoptosis induced by osimertinib treatment of H1975-OE cells (C) and increased apoptosis induced by osimertinib treatment in H1975-KO#1 cells (D). Cells were incubated 24 hours with osimertinib (100 nmol/L) prior to analysis of caspase-3/7 activity, which was normalized to cell number. E, Changes in caspase-3/7 activity elicited by osimertinib treatment of H1975-KO#1 transiently overexpressing CD74 or empty vector. Cells were treated with osimertinib (100 nmol/L) in the presence of MIF (50 ng/mL) for 24 hours prior to determination of caspase activities, which were normalized to cell number. F, Comparison of time course to acquire resistance to osimertinib in H1975-KO and control cells. Cells were treated with chronic exposure of osimertinib at gradually increasing doses. G, Analysis of osimertinib resistance in mouse xenograft tumor models harboring CD74-knockout (KO) or -overexpressing (OE) H1975 cells. Left, nude mice bearing H1975-KO or control cells (n = 6 per group) were treated with vehicle or osimertinib (5 mg/kg orally once daily) for 11 days, followed by osimertinib withdrawal over the next 14 days. Mice were monitored for changes in tumor volume. Right, nude mice bearing H1975-OE or control cells (n = 5 per group) were comparably analyzed, except that drug was administered for 28 days, followed by withdrawal over the next 14 days. In both graphs, data are presented as the mean percent change in tumor volume ± SD. *, P < 0.05; **, P < 0.005.

Figure 5.

CD74 upregulation by osimertinib contributes to a drug-tolerant state in EGFR-mutant tumor cells. A, Western blotting showing upregulation of CD74 and BCL-XL following osimertinib treatment of H1975 cells overexpressing CD74 (OE) or control H1975 cells (Cntl). Both groups were incubated for 0, 24, and 48 hours with osimertinib at 100 nmol/L. Long exposure increased CD74 expression in control cells (Cntl). B, Western blotting showing CD74 knockout efficiency and BCL-XL downregulation in H1975 CD74 KO cells. Two different KO clones (KO#1 and KO#2) and control cells were treated with osimertinib at 100 nmol/L for 0, 24, and 48 hours in the presence of macrophage MIF (50 ng/mL). C and D, Caspase-3/7 activity indicative of inhibition of apoptosis induced by osimertinib treatment of H1975-OE cells (C) and increased apoptosis induced by osimertinib treatment in H1975-KO#1 cells (D). Cells were incubated 24 hours with osimertinib (100 nmol/L) prior to analysis of caspase-3/7 activity, which was normalized to cell number. E, Changes in caspase-3/7 activity elicited by osimertinib treatment of H1975-KO#1 transiently overexpressing CD74 or empty vector. Cells were treated with osimertinib (100 nmol/L) in the presence of MIF (50 ng/mL) for 24 hours prior to determination of caspase activities, which were normalized to cell number. F, Comparison of time course to acquire resistance to osimertinib in H1975-KO and control cells. Cells were treated with chronic exposure of osimertinib at gradually increasing doses. G, Analysis of osimertinib resistance in mouse xenograft tumor models harboring CD74-knockout (KO) or -overexpressing (OE) H1975 cells. Left, nude mice bearing H1975-KO or control cells (n = 6 per group) were treated with vehicle or osimertinib (5 mg/kg orally once daily) for 11 days, followed by osimertinib withdrawal over the next 14 days. Mice were monitored for changes in tumor volume. Right, nude mice bearing H1975-OE or control cells (n = 5 per group) were comparably analyzed, except that drug was administered for 28 days, followed by withdrawal over the next 14 days. In both graphs, data are presented as the mean percent change in tumor volume ± SD. *, P < 0.05; **, P < 0.005.

Close modal

These results suggest that CD74 upregulation by osimertinib treatment suppresses apoptotic activity in EGFR-mutant tumor cells, and that effect is due, at least in part, to upregulation of BCL-XL.

CD74 upregulation is associated with acquisition of osimertinib resistance

Next, to determine whether CD74 expression contributes to emergence of resistance to osimertinib, we treated H1975-KO or -OE cells with increasing concentrations of osimertinib. We observed that a longer time period was required for H1975-KO cells to acquire osimertinib resistance than that seen in control cells (Fig. 5F), whereas H1975-OE cells required a relatively shorter period of time to become resistant (Supplementary Fig. S5H). Next, we examined the effects of osimertinib on anchorage-independent growth of H1975-OE or -KO cells using a soft agar colony formation assay. Although we observed increased colony numbers of osimertinib-treated H1975-OE cells relative to control cells (Supplementary Fig. S5I, left), and comparably decreased colony numbers in osimertinib-treated H1975-KO cells (Supplementary Fig. S5I, right), those differences were not statistically significant.

Next, we established mouse xenograft models by injecting H1975-KO, H1975-OE, or control cells SC into nude mice. After tumors became detectable, we administered osimertinib and monitored tumor growth based on volume. Tumor growth was analyzed in two phases, starting with an initial period of daily osimertinib treatment followed by several days of drug withdrawal. In the case of H1975-KO xenografts, mice were treated 11 days with osimertinib followed by a 14-day withdrawal period. During the treatment period, tumors in mice implanted with H1975-KO cells significantly shrank relative to those seen in animals implanted with control cells (Fig. 5G, left). Moreover, during the 14-day withdrawal period, tumors regrew more rapidly in control compared with H1975-KO xenografted mice, suggesting that a greater number of residual tumor cells persisted in mice bearing tumors comprised of H1975 control cells. We then performed comparable analysis in mice bearing CD74-overexpressing xenografts, although animals were drug-treated for 28 days prior to withdrawal. Overall, we observed a comparable decrease in tumor volume in both groups during the treatment period, whereas H1975-OE tumors regrew significantly more rapidly than did tumors in mice bearing control tumors during the 14-day withdrawal period (Fig. 5G, right). Taken together, our results suggest that CD74 upregulation contributes to osimertinib resistance.

scRNA-seq analysis of clinical samples identifies clusters showing AURKA, VIM, or CD74 upregulation

Finally, we evaluated the clinical relevance of results obtained from cell line models by analyzing specimens from patients with EGFR mutant lung cancer whose tumors had become EGFR-TKI resistant. Characteristics and treatment of patients are summarized in Table 1. Tumor cells were isolated from liver biopsy (liver metastasis) samples, pleural effusion, and transbronchial biopsy samples and then subjected to scRNA-seq analysis. The data were obtained for 3,763, 2,651, 747, and 936 cells from Pt-1, Pt-2, Pt-3, and Pt-4 specimens, respectively (Table 1). After clustering analysis, clusters were annotated on the basis of representative marker gene expression. EPCAM, PTPRC, VCAM1, and VWF served as markers of epithelial cells, CD45+ immune cells, cancer-associated fibroblasts (CAF), and endothelial cells, respectively (Supplementary Fig. S6). Although fewer tumor cells were detected in pleural effusion than in biopsy samples, we were able to analyze the transcriptome by scRNA-seq in all samples (Fig. 6).

Table 1.

Patient characteristics.

CharacteristicPt-1Pt-2Pt-3Pt-4
Histology Adenocarcinoma Adenocarcinoma Adenocarcinoma Adenocarcinoma 
EGFR Mutations EGFR exon19del EGFR exon19del EGFR L858R+T790M EGFR L858R+T790M 
Sex 
Age 76 42 80 57 (43 at diagnosis) 
Smoking 20 pack-years 5 pack-years Never Never 
Treatments TKI failure (Osimertinib failure) TKI failure (Afatinib failure) TKI failure (Gefitinib failure) TKI failure (Osimertinib failure) 
Sample types Liver biopsy Pleural effusion Bronchoscopy biopsy Bronchoscopy biopsy 
Raw reads 354,246,362 312,086,257 243,573,120 444,150,695 
Total cells after cut-off 3,763 2,651 747 936 
Included tumor cells 3,697 118 237 166 
CharacteristicPt-1Pt-2Pt-3Pt-4
Histology Adenocarcinoma Adenocarcinoma Adenocarcinoma Adenocarcinoma 
EGFR Mutations EGFR exon19del EGFR exon19del EGFR L858R+T790M EGFR L858R+T790M 
Sex 
Age 76 42 80 57 (43 at diagnosis) 
Smoking 20 pack-years 5 pack-years Never Never 
Treatments TKI failure (Osimertinib failure) TKI failure (Afatinib failure) TKI failure (Gefitinib failure) TKI failure (Osimertinib failure) 
Sample types Liver biopsy Pleural effusion Bronchoscopy biopsy Bronchoscopy biopsy 
Raw reads 354,246,362 312,086,257 243,573,120 444,150,695 
Total cells after cut-off 3,763 2,651 747 936 
Included tumor cells 3,697 118 237 166 
Figure 6.

Analysis of clinical specimens shows tumor heterogeneity and confirms gene expression changes identified in lung cancer lines. A, C, E, and G, t-SNE plots showing all isolated single cells (left; purple) and tumor cells (right) from indicated patients. B, D, F, and H, violin plots showing expression of AURKA, TPX2, VIM, AXL, and CD74 in each tumor cluster isolated from Pt-1 (B), Pt-2 (D), Pt-3 (F), and Pt-4 (H).

Figure 6.

Analysis of clinical specimens shows tumor heterogeneity and confirms gene expression changes identified in lung cancer lines. A, C, E, and G, t-SNE plots showing all isolated single cells (left; purple) and tumor cells (right) from indicated patients. B, D, F, and H, violin plots showing expression of AURKA, TPX2, VIM, AXL, and CD74 in each tumor cluster isolated from Pt-1 (B), Pt-2 (D), Pt-3 (F), and Pt-4 (H).

Close modal

Tumor cells from Pt-1 or Pt-4 were isolated from osimertinib-resistant tumors, and are thus considered OR cells. Cells from Pt-2 or Pt-3 were isolated from tumors resistant to afatinib or gefitinib, respectively, and are therefore considered ER-type cells. EGFR-T790M was detected in Pt-3 and Pt-4 but neither Pt-1 nor Pt-2 showed acquired resistance mutations in the EGFR gene. Tumor cells from Pt-1 were clustered into nine subgroups (Fig. 6A). Violin plots showed high AURKA and TPX2 expression in clusters 3 and 7, indicating that aurora kinase A may function in osimertinib resistance. VIM was expressed in cluster 8, suggesting that these cells were EMT-like. CD74 was also expressed in clusters 6 and 8 (Fig. 6B). Tumor cells from Pt-2 were clustered into two subgroups (Fig. 6C), and violin plots indicated high VIM and CD74 expression in both (Fig. 6D). Tumor cells from Pt-3 or Pt-4 were clustered into four subgroups (Fig. 6E and G). Interestingly, in both cases, VIM and CD74 upregulation was reciprocal, with VIM upregulated in cluster 2 and CD74 in clusters 0, 1, and 3 (Fig. 6F and H). These suggest that CD74 contributes to drug resistance independent of the EMT. Interestingly, AURKA or TPX2 upregulation was not detected in Pt-2, Pt-3, or Pt-4 (Fig. 6D, F, and H), suggesting that AURKA-related mechanisms underlie osimertinib-specific resistance.

Tumor heterogeneity and activities that confer drug resistance are major hurdles for successful cancer treatment. Given that essentially all patients who suffer advanced NSCLC with EGFR mutations exhibit EGFR-TKI resistance, it is essential to elucidate relevant mechanisms to provide strategies to prevent them. In this study, we sought to assess intratumoral heterogeneity that leads to intrinsic and/or acquired resistance and define mechanisms by which DTPs emerge using single-cell sequencing technology. By performing scRNA-seq and scATAC-seq analysis using cell line models and clinical specimens, we identified CD74 as a novel key factor functioning in the transition from a naïve to a drug-tolerant state of tumor cells. In addition, we also assessed known resistance mechanisms induced by aurora kinase A or EMT. Finally, we successfully identified cell clusters that express AURKA, VIM, or CD74 in clinical specimens.

Osimertinib is widely used to treat patients with EGFR-sensitizing mutations, with or without T790M mutations. A subset of patients whose tumors harbor EGFR-sensitizing mutations remain intrinsically resistant to TKIs (37). Although previous reports show that co-occurrence of ctnnb1 mutations with EGFR mutations predicts intrinsic resistance (38, 39) to EGFR TKIs, little information is available to predict intrinsic resistance to EGFR-TKIs of any generation. A previous report demonstrates that activation of aurora kinase A and its co-activator TPX2 leads to acquired resistance to osimertinib and rociletinib, two third-generation EGFR-TKIs (22). Accordingly, we demonstrated upregulation of AURKA and TPX2 expression in H1975 and PC9 cells during osimertinib treatment. Of note, some cell clusters among H1975 parental cells showed upregulation of both genes, indicating that AURKA and TPX2 expression may contribute to both intrinsic and acquired resistance to osimertinib. Interestingly, neither AURKA nor TPX2 was upregulated in any clusters of PC9-EROR cells or in clinical samples resistant to gefitinib or afatinib, supporting the idea that resistance induced by aurora kinase A may be specific to first-line osimertinib. Moreover, there were several distinct clusters that did not show upregulation of AURKA nor TPX2 among OR30 or OR2,000 cells, each of which may exhibit a unique osimertinib-resistance mechanism. Efforts to characterize each cluster are ongoing.

Several studies suggest mechanisms underlying DTP cell emergence, including the EMT (28), the NOTCH3–β-catenin axis (40), ERK/YAP pathways (41) and miR-147b activity (42). Of these, EMT is the most studied. Our scATAC-seq identified CREB1, but not ZEB1, as potentially responsible for increased VIM expression. Interestingly, PC9-ER cells already exhibit EMT-like properties together with EGFR-T790M; thus, it is unlikely that the EMT contributes to emergence of osimertinib resistance.

We identified CD74 as a potential inducer of a drug-tolerant state; thus, CD74 inhibition may be a novel strategy to suppress DTP emergence. In normal tissue, CD74 is expressed on HLA class II-positive cells. In cancer cells, it is expressed in more than 90% of B-cell malignancies and in some solid tumors, including NSCLC (34). CD74 is also a known partner of ROS1 (43) or NRG1 (44) fusions in NSCLC. We report here CD74 upregulation in the drug-tolerant state in some clusters of H1975 and PC9-EROR cells. Cells with high CD74 expression did not overlap with those with high VIM expression, suggesting that CD74 independently contributes to DTP emergence. Although CD74 may indirectly induce the EMT via IL6 expression, based on our catalog of gene coexpression network study, we predict that CD74 may promote DTP emergence due to its ability to inhibit apoptosis by inducing antiapoptotic proteins such as BCL-XL. Indeed, we showed that CD74 upregulation in H1975 cells induced BCL-XL expression and suppressed apoptosis. Furthermore, CD74 deletion significantly delayed tumor regrowth after osimertinib withdrawal in vivo, whereas CD74 overexpression enhanced (Fig. 5G). Further investigation is required to define these mechanisms.

In summary, we used single-cell sequencing technology to identify intratumor heterogeneity in resistant cell line models and clinical specimens. Our analysis opens a new avenue for identification of target molecules in naïve and drug-tolerant states, which may eventually prevent emergence or resistance to targeted therapies.

Y. Kashima reports grants from JSPS KAKENHI during the conduct of the study. H. Udagawa reports personal fees from AstraZeneca during the conduct of the study. H. Izumi reports grants from JSPS KAKENHI (JP20K17215) outside the submitted work and research support from Amgen and personal fees (honoraria) from Ono. Y. Shibata reports grants and personal fees from Pfizer Inc., AstraZeneca K.K., and Ono Pharmaceutical Co., Ltd.; personal fees from Bristol-Myers Squibb K.K. and Taiho Pharmaceutical Co. outside the submitted work. K. Tanaka reports grants from JSPS outside the submitted work. A. Ohashi reports nonfinancial support from Takeda Pharmaceutical Company Ltd. and Arjuna Therapeutics, personal fees from Ono Pharmaceutical Company Ltd. and Craif Inc., grants from Astellas Pharma Global Development Inc., and Daiichi Sankyo Company, Ltd. outside the submitted work; and reports employment with Takeda Pharmaceutical Company Ltd. K. Goto reports grants and personal fees from Amgen Astellas BioPharma K.K., Amgen Inc., AstraZeneca K.K., Boehringer Ingelheim Japan, Inc., Bristol-Myers Squibb K.K., Chugai Pharmaceutical Co., Ltd., Daiichi Sankyo Co., Ltd., Eisai Co., Ltd., Eli Lilly Japan K.K., Janssen Pharmaceutical K.K., Kyowa Kirin Co., Ltd., Novartis Pharma K.K., MSD K.K., Ono Pharmaceutical Co., Ltd., Pfizer Japan Inc., Taiho Pharmaceutical Co., Ltd., and Takeda Pharmaceutical Co., Ltd.; personal fees from Amoy Diagnosties Co., Ltd., Guardant Health Inc., Life Technologies Japan Ltd., Otsuka Pharmaceutical Co., Ltd.; grants from Astellas Pharma Inc., Ignyta, Inc., Kissei Pharmaceutical Co., Ltd., Loxo Oncology, Inc., Medical & Biological Laboratories Co., Ltd., Merck Biopharma Co., Ltd., Merus N.V., NEC Corporation, Sumitomo Dainippon Pharma Co., Ltd., Spectrum Pharmaceuticals, Inc., Sysmex Corporation, Haihe Biopharma Co., Ltd. outside the submitted work. K. Tuchihara reports grants from National Cancer Center Research and Development Fund during the conduct of the study; personal fees from Takeda, personal fees from Chugai, personal fees from Taiho, personal fees from Novartis, personal fees from Bristol-Myers Squibb, and personal fees from AstraZeneca outside the submitted work. S.S. Kobayashi reports grants from Japan Society for the Promotion of Science, NCI, and National Cancer Center Research and Development Fund during the conduct of the study, grants and personal fees from Boehringer Ingelheim, grants from MiNA Therapeutics, Taiho Therapeutics, personal fees from Bristol-Myers Squibb, Takeda Pharmaceuticals, AstraZeneca, and Chugai Pharmaceutical outside the submitted work. No disclosures were reported by the other authors.

Y. Kashima: Data curation, formal analysis, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. D. Shibahara: Data curation, formal analysis, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. A. Suzuki: Data curation, formal analysis, validation, investigation, writing–review and editing. K. Muto: Data curation, software, investigation, methodology, writing–review and editing. I.S. Kobayashi: Formal analysis, investigation, methodology. D. Plotnick: Investigation. H. Udagawa: Data curation, investigation, writing–review and editing. H. Izumi: Validation, investigation, methodology. Y. Shibata: Validation, investigation, methodology. K. Tanaka: Validation, investigation. M. Fujii: Data curation, investigation, writing–review and editing. A. Ohashi: Formal analysis, writing–review and editing. M. Seki: Data curation, formal analysis, writing–review and editing. K. Goto: Resources, supervision, writing–review and editing. K. Tuchihara: Conceptualization, supervision, writing–review and editing. Y. Suzuki: Conceptualization, resources, supervision, writing–review and editing. S.S. Kobayashi: Conceptualization, resources, data curation, formal analysis, supervision, funding acquisition, validation, investigation, writing–original draft, project administration, writing–review and editing.

The authors thank all members of the Kobayashi laboratory for helpful discussions. They also thank Dr. Pasi Jänne for PC9 cells. This work was supported by MEXT KAKENHI grant no. 19K16879 (to Y. Kashima), JSPS KAKENHI grant no. 16K21746 (to S.S. Kobayashi), the National Cancer Center Research and Development Fund 31-A-6 (to S.S. Kobayashi), and NIH 1R01CA240257 (to S.S. Kobayashi). The authors thank the University of Tokyo for divided fabrication and hgc for computer.

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.

1.
Citri
A
,
Yarden
Y
. 
EGF-ERBB signalling: towards the systems level
.
Nat Rev Mol Cell Biol
2006
;
7
:
505
16
.
2.
Nguyen
KS
,
Kobayashi
S
,
Costa
DB
. 
Acquired resistance to epidermal growth factor receptor tyrosine kinase inhibitors in non-small-cell lung cancers dependent on the epidermal growth factor receptor pathway
.
Clin Lung Cancer
2009
;
10
:
281
9
.
3.
Pagliarini
R
,
Shao
W
,
Sellers
WR
. 
Oncogene addiction: pathways of therapeutic response, resistance, and road maps toward a cure
.
EMBO Rep
2015
;
16
:
280
96
.
4.
Gazdar
AF
. 
Activating and resistance mutations of EGFR in non-small-cell lung cancer: role in clinical response to EGFR tyrosine kinase inhibitors
.
Oncogene
2009
;
28
:
S24
31
.
5.
Kobayashi
S
,
Boggon
TJ
,
Dayaram
T
,
Janne
PA
,
Kocher
O
,
Meyerson
M
, et al
EGFR mutation and resistance of non-small-cell lung cancer to gefitinib
.
N Engl J Med
2005
;
352
:
786
92
.
6.
Pao
W
,
Miller
VA
,
Politi
KA
,
Riely
GJ
,
Somwar
R
,
Zakowski
MF
, et al
Acquired resistance of lung adenocarcinomas to gefitinib or erlotinib is associated with a second mutation in the EGFR kinase domain
.
PLoS Med
2005
;
2
:
e73
.
7.
Bean
J
,
Brennan
C
,
Shih
JY
,
Riely
G
,
Viale
A
,
Wang
L
, et al
MET amplification occurs with or without T790M mutations in EGFR mutant lung tumors with acquired resistance to gefitinib or erlotinib
.
PNAS
2007
;
104
:
20932
7
.
8.
Engelman
JA
,
Zejnullahu
K
,
Mitsudomi
T
,
Song
Y
,
Hyland
C
,
Park
JO
, et al
MET amplification leads to gefitinib resistance in lung cancer by activating ERBB3 signaling
.
Science
2007
;
316
:
1039
43
.
9.
Blakely
CM
,
Bivona
TG
. 
Resiliency of lung cancers to EGFR inhibitor treatment unveiled, offering opportunities to divide and conquer EGFR inhibitor resistance
.
Cancer Discov
2012
;
2
:
872
5
.
10.
Oser
MG
,
Niederst
MJ
,
Sequist
LV
,
Engelman
JA
. 
Transformation from non-small-cell lung cancer to small-cell lung cancer: molecular drivers and cells of origin
.
Lancet Oncol
2015
;
16
:
e165
72
.
11.
Suda
K
,
Tomizawa
K
,
Fujii
M
,
Murakami
H
,
Osada
H
,
Maehara
Y
, et al
Epithelial to mesenchymal transition in an epidermal growth factor receptor-mutant lung cancer cell line with acquired resistance to erlotinib
.
J Thorac Oncol
2011
;
6
:
1152
61
.
12.
Cross
DA
,
Ashton
SE
,
Ghiorghiu
S
,
Eberlein
C
,
Nebhan
CA
,
Spitzler
PJ
, et al
AZD9291, an irreversible EGFR TKI, overcomes T790M-mediated resistance to EGFR inhibitors in lung cancer
.
Cancer Discov
2014
;
4
:
1046
61
.
13.
Ramalingam
SS
,
Vansteenkiste
J
,
Planchard
D
,
Cho
BC
,
Gray
JE
,
Ohe
Y
, et al
Overall survival with osimertinib in untreated, EGFR-mutated advanced NSCLC
.
N Engl J Med
2020
;
382
:
41
50
.
14.
Soria
JC
,
Ohe
Y
,
Vansteenkiste
J
,
Reungwetwattana
T
,
Chewaskulyong
B
,
Lee
KH
, et al
Osimertinib in untreated EGFR-mutated advanced non-small-cell lung cancer
.
N Engl J Med
2018
;
378
:
113
25
.
15.
Oxnard
GR
,
Hu
Y
,
Mileham
KF
,
Husain
H
,
Costa
DB
,
Tracy
P
, et al
Assessment of resistance mechanisms and clinical implications in patients with EGFR T790M-positive lung cancer and acquired resistance to osimertinib
.
JAMA Oncol
2018
;
4
:
1527
34
.
16.
Leonetti
A
,
Sharma
S
,
Minari
R
,
Perego
P
,
Giovannetti
E
,
Tiseo
M
. 
Resistance mechanisms to osimertinib in EGFR-mutated non-small cell lung cancer
.
Br J Cancer
2019
;
121
:
725
37
.
17.
Hata
AN
,
Niederst
MJ
,
Archibald
HL
,
Gomez-Caraballo
M
,
Siddiqui
FM
,
Mulvey
HE
, et al
Tumor cells can follow distinct evolutionary paths to become resistant to epidermal growth factor receptor inhibition
.
Nat Med
2016
;
22
:
262
9
.
18.
Sharma
SV
,
Lee
DY
,
Li
B
,
Quinlan
MP
,
Takahashi
F
,
Maheswaran
S
, et al
A chromatin-mediated reversible drug-tolerant state in cancer cell subpopulations
.
Cell
2010
;
141
:
69
80
.
19.
Terai
H
,
Kitajima
S
,
Potter
DS
,
Matsui
Y
,
Quiceno
LG
,
Chen
T
, et al
ER stress signaling promotes the survival of cancer "persister cells" tolerant to EGFR tyrosine kinase inhibitors
.
Cancer Res
2018
;
78
:
1044
57
.
20.
Dagogo-Jack
I
,
Shaw
AT
. 
Tumour heterogeneity and resistance to cancer therapies
.
Nat Rev Clin Oncol
2018
;
15
:
81
94
.
21.
Kashima
Y
,
Suzuki
A
,
Liu
Y
,
Hosokawa
M
,
Matsunaga
H
,
Shirai
M
, et al
Combinatory use of distinct single-cell RNA-seq analytical platforms reveals the heterogeneous transcriptome response
.
Sci Rep
2018
;
8
:
3482
.
22.
Shah
KN
,
Bhatt
R
,
Rotow
J
,
Rohrberg
J
,
Olivas
V
,
Wang
VE
, et al
Aurora kinase A drives the evolution of resistance to third-generation EGFR inhibitors in lung cancer
.
Nat Med
2019
;
25
:
111
8
.
23.
Ogino
A
,
Kitao
H
,
Hirano
S
,
Uchida
A
,
Ishiai
M
,
Kozuki
T
, et al
Emergence of epidermal growth factor receptor T790M mutation during chronic exposure to gefitinib in a non small cell lung cancer cell line
.
Cancer Res
2007
;
67
:
7807
14
.
24.
Butler
A
,
Hoffman
P
,
Smibert
P
,
Papalexi
E
,
Satija
R
. 
Integrating single-cell transcriptomic data across different conditions, technologies, and species
.
Nat Biotechnol
2018
;
36
:
411
20
.
25.
Wani
MC
,
Taylor
HL
,
Wall
ME
,
Coggon
P
,
McPhail
AT
. 
Plant antitumor agents. VI. The isolation and structure of taxol, a novel antileukemic and antitumor agent from Taxus brevifolia
.
J Am Chem Soc
1971
;
93
:
2325
7
.
26.
Rusan
M
,
Li
K
,
Li
Y
,
Christensen
CL
,
Abraham
BJ
,
Kwiatkowski
N
, et al
Suppression of adaptive responses to targeted cancer therapy by transcriptional repression
.
Cancer Discov
2018
;
8
:
59
73
.
27.
Thomson
S
,
Buck
E
,
Petti
F
,
Griffin
G
,
Brown
E
,
Ramnarine
N
, et al
Epithelial to mesenchymal transition is a determinant of sensitivity of non-small-cell lung carcinoma cell lines and xenografts to epidermal growth factor receptor inhibition
.
Cancer Res
2005
;
65
:
9455
62
.
28.
Raoof
S
,
Mulford
IJ
,
Frisco-Cabanos
H
,
Nangia
V
,
Timonina
D
,
Labrot
E
, et al
Targeting FGFR overcomes EMT-mediated resistance in EGFR mutant non-small cell lung cancer
.
Oncogene
2019
;
38
:
6399
413
.
29.
Taniguchi
H
,
Yamada
T
,
Wang
R
,
Tanimura
K
,
Adachi
Y
,
Nishiyama
A
, et al
AXL confers intrinsic resistance to osimertinib and advances the emergence of tolerant cells
.
Nat Commun
2019
;
10
:
259
.
30.
Zhang
Z
,
Lee
JC
,
Lin
L
,
Olivas
V
,
Au
V
,
LaFramboise
T
, et al
Activation of the AXL kinase causes resistance to EGFR-targeted therapy in lung cancer
.
Nat Genet
2012
;
44
:
852
60
.
31.
Buenrostro
JD
,
Wu
B
,
Litzenburger
UM
,
Ruff
D
,
Gonzales
ML
,
Snyder
MP
, et al
Single-cell chromatin accessibility reveals principles of regulatory variation
.
Nature
2015
;
523
:
486
90
.
32.
Xu
X
,
Zhu
Y
,
Liang
Z
,
Li
S
,
Xu
X
,
Wang
X
, et al
c-Met and CREB1 are involved in miR-433-mediated inhibition of the epithelial-mesenchymal transition in bladder cancer by regulating Akt/GSK-3beta/Snail signaling
.
Cell Death Dis
2016
;
7
:
e2088
.
33.
Schroder
B
. 
The multifaceted roles of the invariant chain CD74–More than just a chaperone
.
Biochim Biophys Acta
2016
;
1863
:
1269
81
.
34.
Stein
R
,
Mattes
MJ
,
Cardillo
TM
,
Hansen
HJ
,
Chang
CH
,
Burton
J
, et al
CD74: a new candidate target for the immunotherapy of B-cell neoplasms
.
Clin Cancer Res
2007
;
13
:
5556s
63s
.
35.
Gil-Yarom
N
,
Radomir
L
,
Sever
L
,
Kramer
MP
,
Lewinsky
H
,
Bornstein
C
, et al
CD74 is a novel transcription regulator
.
PNAS
2017
;
114
:
562
7
.
36.
Suzuki
A
,
Onodera
K
,
Matsui
K
,
Seki
M
,
Esumi
H
,
Soga
T
, et al
Characterization of cancer omics and drug perturbations in panels of lung cancer cells
.
Sci Rep
2019
;
9
:
19529
.
37.
Maemondo
M
,
Inoue
A
,
Kobayashi
K
,
Sugawara
S
,
Oizumi
S
,
Isobe
H
, et al
Gefitinib or chemotherapy for non-small-cell lung cancer with mutated EGFR
.
N Engl J Med
2010
;
362
:
2380
8
.
38.
Blakely
CM
,
Watkins
TBK
,
Wu
W
,
Gini
B
,
Chabon
JJ
,
McCoach
CE
, et al
Evolution and clinical impact of co-occurring genetic alterations in advanced-stage EGFR-mutant lung cancers
.
Nat Genet
2017
;
49
:
1693
704
.
39.
T de Montpreville
V
,
Lacroix
L
,
Rouleau
E
,
Mamodaly
M
,
Leclerc
J
,
Tutuianu
L
, et al
Non-small cell lung carcinomas with CTNNB1 (beta-catenin) mutations: a clinicopathological study of 26 cases
.
Ann Diagn Pathol
2020
;
46
:
151522
.
40.
Arasada
RR
,
Shilo
K
,
Yamada
T
,
Zhang
J
,
Yano
S
,
Ghanem
R
, et al
Notch3-dependent beta-catenin signaling mediates EGFR TKI drug persistence in EGFR mutant NSCLC
.
Nat Commun
2018
;
9
:
3198
.
41.
Kurppa
KJ
,
Liu
Y
,
To
C
,
Zhang
T
,
Fan
M
,
Vajdi
A
, et al
Treatment-induced tumor dormancy through YAP-mediated transcriptional reprogramming of the apoptotic pathway
.
Cancer Cell
2020
;
37
:
104
22
.
42.
Zhang
WC
,
Wells
JM
,
Chow
KH
,
Huang
H
,
Yuan
M
,
Saxena
T
, et al
miR-147b-mediated TCA cycle dysfunction and pseudohypoxia initiate drug tolerance to EGFR inhibitors in lung adenocarcinoma
.
Nature metabolism
2019
;
1
:
460
74
.
43.
Takeuchi
K
,
Soda
M
,
Togashi
Y
,
Suzuki
R
,
Sakata
S
,
Hatano
S
, et al
RET, ROS1 and ALK fusions in lung cancer
.
Nat Med
2012
;
18
:
378
81
.
44.
Fernandez-Cuesta
L
,
Plenker
D
,
Osada
H
,
Sun
R
,
Menon
R
,
Leenders
F
, et al
CD74-NRG1 fusions in lung adenocarcinoma
.
Cancer Discov
2014
;
4
:
415
22
.