Single-Cell Evolutionary Analysis Reveals Drivers of Plasticity and Mediators of Chemoresistance in Small Cell Lung Cancer

This article is featured in Selected Articles from This Issue, p. 865

Small cell lung cancer (SCLC) is a devastating disease characterized by a 5-year survival rate of only 6% (1). The majority of patients (75%) present with extensive stage disease at diagnosis, and survival for these patients is generally under 1 year (2). Most patients are treated with first-line chemotherapy as the disease is commonly observed as extensive stage disease at diagnosis, and therefore surgical resection is rare (3). While initial response to chemotherapy generally appears favorable, the majority of patients have tumors that will rapidly acquire chemoresistance and relapse within months (2). Despite the therapeutic advances made in other tumor types, the grim outlook for patients with SCLC has been largely unchanged in the last 60 years (2, 4, 5). Many clinical trials have investigated the feasibility to utilize targeted therapies in SCLC, such as inhibition of PARP and other checkpoint inhibitors, inhibition of apoptotic pathways, Notch pathway inhibitors and others have shown either minimal responses or negative toxicities (6, 7). Recently, a PD-L1 inhibitor was approved for patients with extensive stage disease, but with an only two-month increase in overall survival time (8), further advances are critically needed.

SCLC is highly correlated with smoking, and almost universal loss of the tumor suppressors RB (RB1) and p53 (TP53) is observed (9). SCLC historically has been thought of as a primarily neuroendocrine disease where patients can be classified in to “classic” or “variant” subtypes, with most patients being the neuroendocrine-high “classic” subtype (10). As our knowledge of SCLC genetics evolves, it has become increasingly clear that patients can be stratified in to one of four molecular subtypes based on key transcription factor expression: NEUROD1, ASCL1, YAP1, or POU2F3 (11). These four molecular subtypes (SCLC-N, SCLC-A, SCLC-Y, SCLC-P, respectively) seem to have a differential effect on tumor behavior and patient outcome, but the regulation and behavior of subtypes is still largely unknown (11–16). Also unknown are the factors that lead to determining which subtype will be prominent, how dynamic an individual tumor's subtype remains during tumor development, and how these factors relate to chemosensitivity/chemoresistance.

Cells in a single tumor are often phenotypically different from one another and therefore contribute differentially to tumor dynamics, a phenomenon known as intratumoral heterogeneity (ITH; refs. 17–19). ITH has been shown to occur in SCLC, particularly postchemotherapy and metastatic events (13–15, 20–26). As tumors progress, the predominant molecular subtype can change, which generally progresses from an ASCL1-high lineage to a NEUROD1-high lineage (14, 15). A number of potential regulators of SCLC identity, including SOX2, MYC, and NOTCH have emerged; however, a definitive lineage has yet to be constructed (13–16, 20). A greater understanding of SCLC ITH is needed, knowing that ITH in SCLC poses a formidable challenge in providing successful treatments. To uncover or design better therapeutics for SCLC, it is first critical that we understand the contribution of different clonal populations to ITH, and how those populations change in response to therapy.

SCLC is almost always driven by loss of tumor suppressors Rb and p53; however, the very early stages of tumor development are not always captured in animal models that are ages out to 3 or even 6 months after tumor initiation. We have yet to understand what happens in the very early stages after tumor initiation. Here we use genetic barcoding and subclonal analysis combined with single-cell RNA sequencing (scRNA-seq) in a mouse model of SCLC and xenografts of both the SCLC-A and SCLC-N subtypes to understand tumor heterogeneity and evolution. By this approach we aim to (i) investigate pathways involved with early tumor development and chemoresistance and (ii) any markers that might be specific to chemoresistant cells. Indeed, we identify a population of tumor cells with high expression of AP-1 network genes that is maintained throughout the course of disease and is required for colony formation. In barcoded SCLC xenografts, the cancer testis antigens (CTA) PAGE5 and GAGE2A were significantly upregulated in chemoresistant tumors and were validated to be mediators of chemoresistance in SCLC. Other members of the CTA family have been successfully targeted for immunotherapy in other tumor types and could represent a promising therapeutic target for SCLC.

Mice were maintained according to the guidelines set forth by the NIH and were housed in the Sanford Research Animal Research Center, accredited by Association for Assessment and Accreditation of Laboratory Animal Care International using protocols reviewed and approved by our local Institutional Animal Care and Use Committee (IACUC). The staining and scoring, as well as storage of the data for all human specimens were approved by the Sanford Health Institutional Review Board.

SCLC cell lines NJH29 (H29), NCI-H82 (H82, CVCL_1591), NCI-H1836 (H1836, CVCL_1498), and NCI-H209 (H209, CVCL_1525) were used in this study. The cells were maintained in suspension and cultured in RPMI with 10% bovine growth serum and penicillin/streptomycin. All cell lines regularly tested negative for Mycoplasma contamination by PCR (quarterly or when thawing a new vial) and validated by IDEXX BioAnalytics (annually). The NHJ29 cell line was a kind gift from Julien Sage.

To clone the retroviral barcoding library, a CAG-GFP retroviral plasmid was used as a backbone (gift from Fred Gage, Addgene plasmid #16664; ref. 27). A poly-A sequence was subcloned from TetO-FUW-sox2, a gift from Rudolf Jaenisch (Addgene plasmid #20326; ref. 28), to the 3′ end of the GFP at the PmeI site using InFusion Cloning (TaKaRa) and screened with PCR and restriction digests (Supplementary Table S8). DNA oligos containing the barcode sequence were ordered from Eurofins and annealed by heating to 95°C for 5 minutes and allowed to cool to room temperature over the course of several hours. The CMV-GFP-polyA plasmid was digested at PmeI and HindIII (added in the polyA cloning step), and the annealed barcode was inserted by InFusion cloning (TaKaRa) 3′ of the GFP and 5′ of the polyA sequence. Forty of the initial colonies were screened for presence of the barcode via PCR, and all four of the four sequences were validated to contain the barcode using Sanger sequencing. Once presence of the barcode was confirmed in a number of colonies, the cloning product was transformed and plated to grow on five 15-cm plates of Luria Broth (LB) agar. Following overnight growth, all colonies were collected and pooled by flushing the plates with prewarmed LB broth. Plasmid libraries were purified, and the product was pooled. Gamma-Retrovirus was produced by cotransfection of 293Ts with the barcoding retrovirus and retroviral packaging plasmid pCL-Amph. Retroviral supernatant was collected 48 and 72 hours after transfection and concentrated using the TaKaRa Retro-X concentrator.

To clone the AAV-LBC plasmid, a single-guide RNA (sgRNA) targeted to Rosa26 was subcloned from pU6-sgRosa26-1_CBh-Cas9-T2A-BFP, a gift from Ralf Kuehn (Addgene plasmid #64216; refs. 29, 30) and inserted in to the AAV-KPL backbone (AAV:ITR-U6-sgRNA(Kras)-U6-sgRNA(p53)-U6-sgRNA(Lkb1)-pEFS-Rluc-2A-Cre-shortPA-KrasG12D_HDRdonor-ITR (AAV-KPL), a gift from Feng Zhang (Addgene plasmid #60224; ref. 31) at SacI and MulI using InFusion cloning (TaKaRa), and screened using PCR and Sanger sequencing. Left and right homology arms to Rosa26, as well as a polyA signal were subcloned from pR26 CAG/GFP Dest, a gift from Ralf Kuehn (Addgene plasmid #74281; refs. 29, 30) in to the AAV9-r26 plasmid at PmlI. The barcode library, attached to a GFP was ordered from Thermo Fisher Scientific GeneArt program, and was inserted in to the AAV9-r26 plasmid at PmlI and BamHI via InFusion cloning. After the first 30 colonies were screened for insertion of the barcode by restriction digest, PCR, and Sanger sequencing, chemically competent cells were transformed with the plasmid pool and plated on ten 10-cm LB plates. After overnight growth, all plates were washed with prewarmed LB and the pooled colonies were purified. The AAV-LBC plasmid was used to generate AAV9 at the University of Michigan Viral Vector Core.

To profile the diversity in the barcode pools, targeted amplicon sequencing was performed on the barcode region. PCR was used to amplify the barcode and add partial adaptors for Illumina sequencing. The minimal number of cycles needed to amplify the barcode region was used to minimize the risk of introducing variants or oversaturate the sample with a limited number of barcodes that had been disproportionately amplified. Samples were sequenced on an Illumina platform at Genewiz. The targeted amplicon sequencing of the barcodes was trimmed, and QC performed via CutAdapt. GREP (bash) was used to identify barcodes and export them to R Studio for analysis. The true number of barcodes was determined by finding the optimal number of clusters in a similarity-based hierarchical clustering model that minimizes the distance in the PCR error rates between the random regions of the barcode and the constant region of the barcode. The total number of unique barcodes was estimated using Chao2 asymptotic richness estimator in R (32).

To ensure sufficient overlap in barcodes between the “pregrowth” sample and xenografts, we sought to determine the optimal number of doublings before sufficient overlap in two independent samples was observed. A total of 1.6 million cells were seeded and barcoded using the CAG-GFP-BC retrovirus and a spinfection at 940 × g for 2 hours. At each doubling for five doublings following spinfection, one well of cells was harvested and split in two. The barcode region was amplified off the cDNA, and partial adapters for Illumina-based sequencing were added. Amplicon sequencing was performed at Genewiz using an Illumina miSeq platform. The barcodes were analyzed using a custom R script after trimming and QC via CutAdapt. To verify the results of the barcode sequencing, computational modeling was used. We simulated the same doubling time experiment 1,000 times using a custom R package. On the basis of the results of the sequencing, modeling, and previously published work (33), three doublings after barcoding gives sufficient overlap between two independent samples and will be used to generate the barcoded xenografts.

The well-characterized Rb^lox/lox, p53^lox/lox, p130^lox/lox SCLC mouse model (34) was bred to the H11^{lox-stop-lox-Cas9} mouse (35) model (Jax #026816) to generate the RPR-Cas9 mouse used in this work. Tumors were initiated by intratracheal injection with Ad-CMV-Cre (Baylor Viral Vector Core) to delete the Rb, p53, and p130 loci, and induce expression of Cas9. At 1-month intervals for 5 months after tumor initiation, the AAV-LBC virus was delivered to separate cohorts of mice via intratracheal instillation (36) to barcode the forming tumors at the Rosa26 locus using the CRISPR-Cas9 system. Mice were euthanized at 1-month intervals after their Ad-CMV-Cre injection, up to 6 months, or when they became moribund according to institutional IACUC guidelines. Two mice received chemotherapy at 5 mg/kg cisplatin and 10 mg/kg etoposide on day 1, and 10 mg/kg etoposide on days 2 and 3 with subcutaneous saline, repeated for 3 weeks, and were allowed to progress until they were moribund (22).

NCI-H209 and NCI-H82 SCLC cell lines were validated as pathogen free (IDEXX BioAnalytics) and then barcoded with the rCMV-GFP-BC retrovirus by spinfection at 940 × g for 2 hours. Each barcoded line was allowed to double three times, to ensure overlap in barcodes in the cells sampled and cells used for xenografting. Immediately prior to xenografting, a portion of the cells were removed to generate the scRNA-seq library (“pregrowth” sample). To make the xenografts, 2,500 cells were mixed in a 1:1 ratio with matrigel (Corning Life Sciences) and injected into the hind flank of NOD-SCID mice. Xenografts were measured daily after the tumors were palpable by hand in the hind flank. After xenografts reached 3 cm³ total volume, mice were euthanized, and the xenografts harvested. Tumors were dissected, dissociated, and a portion of the cells were used for scRNA-seq (“prechemotherapy” sample). The remainder of cells were mixed in a 1:1 ratio with matrigel and injected into a new NOD-SCID mouse. The mice that received these serial xenografts received chemotherapy at 5 mg/kg cisplatin and 10 mg/kg etoposide on day 1, and 10 mg/kg etoposide on days 2 and 3, repeated for 3 weeks after tumors were palpable to generate chemoresistant xenografts (22). When these mice reached a total tumor burden of 3 cm³, they were euthanized, and the tumors were dissociated and used to generate scRNA-seq libraries (“postchemotherapy” sample). Xenograft growth was plotted by determining the best fit between an exponential growth model or a cubic polynomial model by Akaike's information criteria. Significance was determined by the extra sum-of-squares F test.

The “pregrowth” from the xenograft model were prepared for scRNA-seq according to the 10X Genomics protocols. After xenografts reached a cumulative volume of 2 cm³, mice were euthanized, and tumors dissected. The MACS (Miltenyi Biotec) human tumor dissociation kit was used to digest the tumors, and cells were prepared for scRNA-seq using the 10X Genomics protocols. After the RPR-Cas9 reached their endpoint, they were euthanized via cervical dislocation and lungs were harvested. Tissues were prepared for flow cytometry according to the 10X Genomics protocols, using the MACS mouse tumor dissociation kit. After dissection, the lungs were sorted for GFP+ cells, which are the barcoded population, and cells were prepared for scRNA-seq according to the 10X Genomics protocol. Tumors were microdissected from three mice by cutting out one tumor lesion, which was sequenced without undergoing FACS to capture the stromal and microenvironment cells. A 10X Genomics Chromium Controller was used for the library preparation of all tumors.

One lobe of the lung and one lobe of the liver from each of the RPR-Cas9 mice were taken for histology to verify the presence of tumors in this sample. Samples were fixed in 4% paraformaldehyde for 15 minutes and transferred to 30% sucrose for 24 to 48 hours. The samples were then embedded and cryosectioned before staining. Tumors were stained for GFP and Cas9 to confirm that the barcoding system was successfully induced by the Adenovirus induction of Cas9 expression and AAV9 induction of GFP expression.

The scRNA-seq data were processed using 10X Genomics CellRanger software (10X Genomics, version 6.0.1). Initially, CellRanger count was run on each sample followed by CellRanger aggregate to generate one expression matrix, a vector of genes and a vector of barcodes. CellRanger was run with the –expect-cells parameter set to the expected number of cells. The expression matrix and the vectors of genes and barcodes were imported to R and analyzed using the package Seurat (ref. 37; version 4.1.1) to find cell clusters and their markers. Prior to clustering cells were further filtered on the basis of percent mitochondrial RNA (autochthonous model percent mitochondrial reads < 20%, NCI-H209 and NCI-H82 xenografts < 15%) and number of genes expressed and cell-cycle differences were regressed from the data. Dimensions for Uniform Manifold Approximation and Projection (UMAP) plots were determined by minimal elbow plot inflections (autochthonous data, 20 dims, NCI-H209 xenografts 19 dims, combined NCI-H209 and NCI-H82 xenografts 20 dims). Monocle3 (refs. 37–40; version 1.2.9) was used similarly to both cluster cells and to run psuedotime analysis. Unbiased cell-type determination was performed using Garnett (ref. 41; version 0.2.8) using the provided mmLung_markers.txt classifier.

For the analysis of clonal evolution of cells and retrieval of lineage barcode (LBC) from cells, sample-level scRNA fastq data were partitioned into per-cell fastq. First cell barcodes from cells that were deemed to be viable in the CellRanger were retrieved using Seurat. These cell barcodes were used to simplify the demultiplexing step and to assign reads to cells using cutadapt (ref. 42; version 2.0). Each cell-specific fastq was mapped to human and to mouse genome using STAR (ref. 43; version 2.7.10a) and the percentage mapping rate was used to identify and filter mouse cells. Presence of LBC in the reads of each cell was checked by using cutadapt to capture reads that possess the constant flanking regions of the LBC. Per cell SNPs and Indels were called using GATK best practices (44). To speed up the analysis and avoid computational limitations, haplotypeCaller was run in the Genomic Variant Call Format (GVCF) mode followed by GenomeicDBImport and GenotypeGVCFs for each gene that was expressed in the scRNA data. Per gene vcf files were finally merged using bcftools (ref. 45; version 1.9) to produce a joint vcf for all cells and genes. Cell relatedness and clonal evolution was investigated using Dendro package in R (46). Muller plots were designed using Fishplot (ref. 47; version 0.5). Gene set enrichment analysis was performed (GSEA) using GSEA (ref. 48; version 4.2.3) by exporting the normalized gene counts from Monocle3 and running a weighted analysis against the c2.all.v7.5.1.symbols_smith.gmt genesets downloaded from the MSigDB without additional normalization from GSEA. Epithelial-to-mesenchymal transition (EMT) scores defined by imogimap (ref. 49; version 0.0.0.9000)

An alamar blue assay was used to determine the IC₅₀ value of cisplatin and etoposide in NCI-H82 and NCI-H209 cell lines. Cells were seeded in 96-well plates and treated with either drug, with concentrations spanning 3 orders of magnitude. Cellular viability was assessed daily via Alamar Blue. The IC₅₀ for cisplatin was determined to be 2.876 μmol/L. The IC₅₀ for etoposide was determined to be 0.110 μmol/L. These values were used for the resulting experiments. SCLC-N lines NJH29 and H82, and SCLC-A lines NCI-H1836 and NCI-H209 were treated with cisplatin and etoposide at the IC₅₀ values for 3 days at cycles resembling the in vivo chemotherapy treatment (cisplatin and etoposide day 1, etoposide only days 2 and 3). Cells were harvested on days 2 and 3 and RNA was extracted following the TRIzol (Ambion Biosciences) protocol, and RNA was converted to cDNA using the NEB ProtoScript Reverse Transcriptase Kit. Expression levels of PAGE5 and GAGE2A were quantified via qPCR.

Short hairpin RNAs (shRNA) targeting PAGE5 and GAGE2A (Supplementary Table S3) were designed using pSicoligoMaker3 (Ventura lab). shRNA oligos were cloned in to the lentiviral backbone pSicoR, a gift from Tyler Jacks (Addgene plasmid #11579; ref. 50). The pSicoR-shPAGE5 and pSicoR-shGAGE2A were made into a second-generation lentivirus using pMD2.G and psPAX2 (Addgene plasmid #12260) as packaging plasmids and concentrated overnight using the TaKaRa Retro-X retroviral concentrator. SCLC-N lines NJH29 and H82, and SCLC-A lines H209 and H1836 were infected with pSicoR-shPAGE5 or pSicoR-shGAGE2A and sorted by GFP expression using the BD FACS Jazz. Knockdown of PAGE5 and GAGE2A expression was validated in the sorted cells with qPCR. To generate a double knockdown, the sorted cells were infected with the reciprocal virus and expression of both PAGE5 and GAGE2A was assessed via qPCR. The single- and double-knockdown cells were used to generate xenografts to investigate the dependence of the chemoresistance phenotype on PAGE5 or GAGE2A expression. A total of 150,000 cells were mixed in a 1:1 ratio with GelTrex (Gibco) and injected in to the hindflank of NOD-SCID mice. Because of supply chain disruptions, a switch from Matrigel to GelTrex was necessary; however, they both function the same way in providing some extracellular matrix to aid in xenograft injection. Tumors were allowed to grow until a cumulative volume of 3 mm³ was reached, at which point mice were euthanized and the tumors kept for histology to validate the knockdown of PAGE5 and GAGE2A. Half of the mice were treated with chemotherapy at 5 mg/kg cisplatin and 10 mg/kg etoposide on day 1, and 10 mg/kg etoposide on days 2 and 3, repeated for 3 weeks after tumors were palpable. In culture, the shPAGE5, shGAGE2A, and double-knockdown cells were treated with the IC₅₀ value of cisplatin for 2 days and viability was assessed via Annexin V and propidium iodide staining by flow cytometry (BioLegend APC Annexin V Apoptosis Detection Kit).

GAGE2A and PAGE5 overexpression retroviruses were generated by amplifying the transgenes from SCLC cell lines and cloning them into the CAG-GFP retroviral backbone. NJH29, NCI-H82, NCI-H1836, and NCI-H209 cells were transfected with either the rCAG-PAGE5-GFP or rCAG-GAGE2A-GFP vectors by spinfection with concentrated virus at 940 × g for 2 hours. Transduced cells were treated with the IC₅₀ value of cisplatin, etoposide, or cisplatin and etoposide for 2 days and response to chemotherapy was assessed by Annexin V and propodium iodide staining, and efficiency of overexpression assessed by qPCR.

To inhibit the AP-1 complex, cJun was knocked down by transfection with a dominant-negative cJun construct. This is a common method for inhibiting the formation of the Jun/Fos AP-1 complex (51, 52). pMIEG3-JunDN was a kind gift from Alexander Dent (Addgene plasmid #40350; ref. 51). NCI-H82, NJH29, NCI-H1836, and NCI-H209 SCLC cell lines were transfected with pMIEG3-JunDN using Lipofectamine 3000. Upon visual GFP detection, cells were sorted using FACS for GFP expressing cells. Cells were seeded in to 6-well plates for a soft agar colony formation assay. Briefly, 0.8% Seaplaque agar (Lonza) was used as a bottom layer and 10,000 cells per well were seeded in 1.2% agar in the top layer. Plates were fed with full RPMI as needed to prevent drying out. Ten days after seeding, colonies were observed by eye and the plates were stained with 0.001% crystal violet for 1 hour, and plates were photographed. The number of crystal violet colonies stained was quantified with a custom CellProfiler script.

Human SCLC biopsies were obtained from the Sanford Health Biobank. Slides were stained with anti-PAGE5 (Invitrogen PA5-50470) or anti-GAGE2A (Aviva Systems ARP64957-P050). Stained slides were scanned using an Apereo AT2 slide scanner. Three independent researchers viewed and scored the scanned slides based on positivity, distribution, and intensity of staining. Positivity was a binary score, with the sample earning a positive score from any singular positive cell. Distribution was scored on a 0 to 3 scale: 0 for 0% to 5% of the sample staining positively, 1 was assigned to samples 6% to 30% positive, 2 for samples 31% to 60%, and 3 for samples more than 60% positive for PAGE5 or GAGE2A. For staining intensity, a score 0 to 3 was assigned. 0 for samples with no PAGE5 or GAGE2A staining, 1 for samples with light staining, 2 for samples with moderate staining, and 3 for samples with intense staining.

Data from this study may be accessed at the Gene Expression Omnibus, accession no. GSE216182.

To barcode tumors as they form in their native environment, we used the well-characterized SCLC mouse model with Cre-mediated deletion of tumor suppressors Rb, p53, and p130 (RPR2; ref. 34), bred with an H11^{lox-stop-lox-Cas9} mouse (35) to generate the RPR2-Cas9 mouse model (Fig. 1A). Upon intratracheal Cre injection, the mice rapidly develop tumors and concurrently express Cas9 (Supplementary Fig. S1). Mice were barcoded by intratracheal AAV injection at 1-month intervals starting at time of Cre injection, through 5 months after tumor initiation, and were euthanized at 1-month intervals following barcoding (Fig. 1A; Supplementary Table S1). Designing the barcoding of lungs in this way allows for the barcoding of tumors at early, middle, and late stages of tumor development. After dissection, barcoded tumor cells were isolated via FACS, and used for scRNA-seq.

SCLC tumor cells were readily identifiable from other cells of the tumor microenvironment by automatic annotation using known lung markers (Supplementary Fig. S2) and were characterized by classic SCLC neuroendocrine markers, Ascl1, and Mycl (Supplementary Fig. S3). Other cell types that are identified by this automatic annotation are myeloid, club, alveolar, ciliated cells, and T and B immune cells (Fig. 1B; Supplementary Fig. S2). Upon analysis of the scRNA-seq data, no barcodes were detectable in the tumors (Supplementary Fig. S1). Despite having tumors that are immunoreactive to antibodies against GFP, and the detection of GFP⁺ cells via flow cytometry, the depth of scRNA-seq in this case was not sufficient to pick up reads from the low-expressed GFP and barcode (Supplementary Fig. S4). Despite not detecting barcodes in these samples, due to the timewise design of the animal studies we then set out to determine what developmental pathways might be altered during SCLC development. When evaluating the cell-cycle composition of the cells within the tumor populations, the proportion of cells in G₁ decreases with each subsequent month that tumors are allowed to form, and the percentage of cells in G₂ or metaphase increases nearly 6-fold after 5 months of tumor development (Fig. 1C; Supplementary Fig. S5). When plotting the distribution of tumor cells by developmental time, we observe an “early” tumor signature that arises in the tumors isolated after 2 months of development and persists through later stages of tumor development (Fig. 1D and E). In the later months of tumor development (month 4 and following), a “late” signature emerges, which eventually is responsible for the majority of the tumor (62.5% ± 7.2%) in the longest developed tumors (Fig. 1D–F). While the differences between the early and late clusters are subtle, indicated by their relative proximity in the UMAP projection, we set out to investigate the transcriptomic differences that separated these cells in the unbiased clustering. The early and late tumor populations can be distinguished by the gene expression of several distinct expression modules (Fig. 1G; Supplementary Table S3). We observe a decrease in neuroendocrine markers (including Ascl1, Chga, Cgrp, and Mycl) in the late tumor cells (modules 1 and 6; refs. 14, 53), and regulation of many cell-cycle genes which could drive the decrease of cells in G1 at later stages (module 4; Fig. 1C). Noteworthy however, was module 5 which is comprised of members of the AP-1 complex, and most strongly (P = 4.16 × 10⁻⁷) distinguishes between the early and late populations (Fig. 1G).

In evaluating the transcriptomic differences between the early and late tumor clusters, members of the AP-1 transcriptional network were not only upregulated in the late population but are also expressed in a majority of the cells indicating it is a potential driver of this group and not simply a rare subtype (Fig. 2A and B). Furthermore, an AP-1–specific geneset was found to be enriched in the late cell cluster [reactome activation of the AP-1 family of transcription factors, normalized enrichment score (NES) = −0.9091, q-value = 0.0000; Supplementary Fig. S6B]. Members of the AP-1 network including JUN, FOS, and ATF proteins are indeed implicated in tumor development (52, 54–63) although their role has not yet been characterized in SCLC. We therefore sought to investigate the impact of disruption of the AP-1 complex on tumor development in SCLC. Using a dominant-negative cJun construct (JunDN; ref. 51), AP-1 signaling was disrupted in the human SCLC cell lines NJH29 and NCI-H82, and cells were sorted by GFP expression from the JunDN vector to capture transduced cells. The cells were seeded into soft agar to assess colony formation abilities. SCLC cells with disruption of AP-1 signaling due to JunDN formed significantly fewer colonies (98.6% reduction, P = 5.0 × 10⁻⁵) than wild-type cells (Fig. 2C and D; Supplementary Fig. S6C) while cell viability was not affected (Fig. 2E). Given that resistance to chemotherapy is a formidable challenge in SCLC, we had set aside two mice to be treated with cisplatin and etoposide at 5 months of tumor progression and performed scRNA-seq on the chemotherapy-treated cells after recurrence at both 6 and 7 months of tumor progression. Using these samples, we queried whether the AP-1 expression was correlated with chemotherapy resistance. The majority (64.2%) of the chemoresistant cells from the tumors treated with chemotherapy had gene expression signatures consistent with the late cell cluster, which has enrichment of AP-1 family members (Fig. 2F and G).

While the autochthonous model was able to give insights into the early mechanisms of SCLC tumor development in the native microenvironment, we set out to investigate mechanisms of human SCLC development. To that end, we again used single-cell clonal analysis, but this time using human SCLC cell lines which were labeled prior to their growth as xenografts in immunocompromised mice. Like the autochthonous barcoding model (Fig. 1A), a barcoding retrovirus was designed that contains GFP, followed by the lineage barcode, and finally a polyA sequence (Fig. 3A). To be able to follow barcode populations over time, SCLC cell lines NCI-H209 (SCLC-A subtype) and NCI-H82 (SCLC-N subtype) were barcoded in culture and the cells grown as xenografts (Fig. 3B). A portion of the cells that were barcoded in culture were assessed by scRNA-seq (“preinjection” sample), and a portion are used to generate xenografts. The initial xenografts that formed served as “prechemotherapy” samples, and upon dissection, a portion are used for scRNA-seq, with the remainder serially transplanted into a new mouse. The mouse harboring this secondary xenograft was then treated with cisplatin and etoposide and tumors harvested and assessed by scRNA-seq to observe how chemotherapy treatment changes the transcriptional profile of these cells (“postchemotherapy” sample). By this system, “preinjection,” “prechemotherapy,” and “postchemotherapy” samples would be used to gain insight into the clonal dynamics in these populations over time and treatment.

Prior to the generation of xenografts, we sought to thoroughly profile the diversity of barcodes in the generated barcode retroviral pool. As this methodology starts with a split of the cells into the preinjection and xenograft samples, it is important to assess the doubling rates of these cell lines as clonal growth of singly labeled cells would be required for a significant overlap of the same barcodes between these samples. The doubling time of the two cell lines used to generate xenografts was therefore determined to be 2.02 days (NCI-H209) or 2.39 days (NCI-H82; Fig. 3C). To determine the even distribution of unique barcoded cellular clones across two evenly divided samples, we used an Illumina targeted amplicon sequencing platform, to sequence the barcode region of the retroviral plasmid pool and the retroviral cDNA (Fig. 3D and E). To estimate the diversity of barcodes most accurately in the two samples, the PCR error rate in the nonvariable region of the LBC was used to differentiate between unique barcodes and artifactual barcodes due to a PCR error (Fig. 3F). Finally, using the Chao2 asymptotic richness estimator (32, 64), to account for subsampling and PCR amplification errors (Supplementary Fig. S7), we determined the diversity of the barcode pool to be around 5,946 unique barcodes, well in excess of the 2,500 cells used to seed each xenograft (Fig. 3E). To ensure the barcodes captured in the preinjection sample and the xenograft have sufficient overlap, prior to xenografting we profiled the barcodes in two independent samples taken from a population of barcoded cells at every doubling (Fig. 3C). As the number of cellular doublings increases, the fractional overlap of barcodes in two halves of a barcoded cellular population increases (48.3% overlap at three doublings), until a point at which the overlap is reduced (32.6% overlap at five doublings), presumably due to unequal growth of subclones or subclone death during culture (Fig. 3G). By assessing the barcode overlap over five doublings, the optimal number of doublings for sufficient overlap is identified to be three doublings (Fig. 3H).

Xenografts were generated by injecting barcoded NCI-H209 or NCI-H82 cell lines into the hind flank of immunocompromised NOD-SCID mice, and growth was monitored over time (Supplementary Fig. S8; Supplementary Table S2). After dissection, tumors were sequenced via scRNA-seq and the remainder of cells were serially xenografted into a new mouse, which was treated with chemotherapy after palpable xenografts were formed (∼36 days for NCI-H209 and ∼19.5 days for NCI-H82) and the subsequent chemotherapy-treated tumors were sequenced via scRNA-seq (Fig. 3A). While the bulk of these cells were of human origin, some cells from an immune or fibroblast mouse origin, as determined by mapping statistics between the human and mouse genome, were collected yet removed for subsequent analyses (Supplementary Fig. S9). The majority of the NCI-H82 (SCLC-N) xenograft cells, regardless of chemotherapy status cluster together, and there are very few changes in the global transcriptomes before and after chemotherapy (Fig. 4A; Supplementary Table S4). In contrast, the NCI-H209 (SCLC-A) postchemotherapy cells display a robust transcriptomic shift from the preinjection and prechemotherapy samplings, and cluster more closely with the SCLC-N xenografts, indicating a shift to a more SCLC-N–like phenotype after chemotherapy (Fig. 4A; Supplementary Table S4). Both the tumor initiation and chemotherapy processes serve as cellular selection events, as the percentage of uniquely barcoded cells in each subsequent sample decreases nearly 10-fold at each stage, indicating only a subset of cells are able to initiate xenografts, and a smaller subset are able to survive chemotherapy and re-establish a tumor (Fig. 4B). However, similar to the autochthonous model, the overall percentage of sequenced LBCs linked to GFP expression was generally low (Supplementary Fig. S8). Therefore, we identified clonal populations in these cells by the identification of single-nucleotide variants (SNV) in the scRNA-seq data, focusing on the NCI-H209 cells as they showed greater dynamics upon chemotherapy exposure. Over time, the clonal distribution within the tumors shifts (Fig. 4C), indicating that this approach is amenable to identifying clonal subpopulations with phenotypic roles during SCLC tumor development. This xenograft study was designed with three unique serial transplantations from xenografts to postchemotherapy xenografts which might each present with their own unique spontaneous SNVs, therefore we sought to further investigate the clonal dynamics over time in these three independent study arms (Fig. 4D; Supplementary Fig. S10; Supplementary Table S5). The three lineage arms show distinct modes of clonal evolution, but ultimately display increasing clonal diversity in the postchemotherapy xenografts. This is reflected in the distribution of more disparate subclones in the chemotherapy-treated xenograft populations (Fig. 4E), yet all arising from the same subpopulations in the original cell line (Supplementary Fig. S10). To assess the relative plasticity of these subclonal populations, we divided the cells into multiple discreet clusters and determined the degree of cluster residence at different stages of tumor development, that is, from cell line to xenograft, or from xenograft to chemotherapy-treated xenograft (Fig. 4F). We observe that subclonal plasticity increases after exposure to chemotherapy (Fig. 4G).

Further exploring the earliest subclones in original NCI-H209 cells, we clearly identify one subclone that does not contribute to xenograft formation (subclones 1a, 2a, and 3a), while the remaining subclones are found to populate the xenografts (Fig. 4D). We therefore identify these subclones (1b, 1c, 2b, 2c, and 3b) as tumor-initiating clone (TIC) as it appears the cells capable of forming a tumor reside in these subclones. Differential gene expression between the cell line subclones readily distinguish TICs from non-TICs (Fig. 5A; Supplementary Fig. S11; Supplementary Table S6). As EMTs are often associated with tumor invasiveness and growth we tested the relative EMT status of these clones and observe that EMT scores derived from the expression of known markers clearly separates TICs from non-TICs with the non-TICs possessing more of an epithelial state and the TICs displaying more of a mesenchymal state (Fig. 5B). This appears to closely relate to MYC activity as the non-TICs are enriched for genes associated with genetic amplification of MYCL and lack of amplification of MYC and MYCN (NES = 2.16, q-value = 2.4 × 10⁻⁴, rank = 5), while the TICs show marked upregulation of genes associated with MYC and MYCN amplifications (NES = −3.47, q-value = <1 × 10⁻⁵, rank = 16; Fig. 5C; ref. 65). We believe that this reflects higher MYCL/MYC activity in the non-TICs/TICs, respectively, and not amplification of these genes (65). SCLC cells that are high in both EpCAM and CD24 and low for CD44 are enriched for long-term tumor propagating cells (66), and consistent with that discovery we observe that the TICs display similar characteristics in the starting cell line; however, upon xenograft formation, the levels of EpCAM are markedly reduced consistent with the acquisition of a more mesenchymal state (Fig. 5B and D). Furthermore, as we observed in the mouse SCLC tumors that the AP-1 complex is highly expressed in tumor development and associated with colony formation (Fig. 2), in the TICs we also observe significant upregulation of JUN family members (JUN, JUND) and ATF5 in growing xenografts, although there is a slight repression in the chemotherapy-treated xenografts which differs from the mouse tumors (Fig. 5E).

Intriguingly, one of the gene families most upregulated upon chemotherapy treatment (five of the top 20 genes) were the CTAs (Supplementary Table S4). CTAs are a large class of proteins almost exclusively expressed in male germ cells and many cancers where they act as oncogenes with key effects on proliferation, genomic stability, invasion, colony formation, and resistance to apoptosis (67–69). Of the broad class of CTA members of the PAGE and GAGE families appear most upregulated upon chemotherapy treatment (Fig. 5F; Supplementary Table S7). We identified GAGE2A and PAGE5 (CT16) with the highest fold increase; however, there is marked upregulation also of PAGE1 (Fig. 5E–J). We therefore sought to determine the role of CTAs in mediating response to chemotherapy. We queried a well-characterized dataset of 120 SCLC patient samples (9) for expression of GAGE2A and PAGE5. While GAGE2A was not annotated in this dataset, the related family member, GAGE2B was included. Consistent with our observations in the xenograft model, both PAGE5 and GAGE2B were only expressed in chemotherapy-treated samples (Supplementary Fig. S12A). To better query expression of GAGE2A specifically, we screened patient biopsies and observed immunopositivity for both GAGE2A and PAGE5, with a more disperse staining of PAGE5 (Fig. 6A). Xenografts made with human SCLC cell lines additionally showed expression of GAGE2A and PAGE5, but only after treatment with cisplatin and etoposide (Supplementary Fig. S12B). Furthermore, SCLC cell lines treated with cisplatin and etoposide in culture also showed significant upregulation of GAGE2A and PAGE5 (Supplementary Fig. S12C), indicating that their expression is most likely linked to chemotherapy treatment, rather than serial transplantation. To understand the requirement for CTA expression on chemoresistance, both GAGE2A and PAGE5 were knocked down by retroviral expression of shRNAs in four SCLC cell lines (NCI-H209 and NCI-H1836, from the SCLC-A subtype; NCI-H82 and NJH29 from the SCLC-N subtype). Cells were treated with chemotherapy in culture (Supplementary Fig. S13) and the percentage of live cells was assessed via flow cytometry. SCLC-A type cell lines with both GAGE2A and PAGE5 knocked down are nearly 2-fold more sensitive to cell death (P = 0.013) caused by cisplatin (Fig. 6B; Supplementary Fig. S13); however, knockdown in SCLC-N type cells showed a trend toward higher sensitivity yet this was not significant (P = 0.087) perhaps due to the inefficiency of the shRNAs to provide a large enough knockdown in the SCLC-N type cells which already possess high CTA expression (Fig. 5I and J). Subsequently, we tested whether overexpression of GAGE2A or PAGE5 via retroviral transduction altered response to cisplatin, etoposide, or combination cisplatin and etoposide. Cells with overexpression of either GAGE2A or PAGE5 were significantly more resistant to cell death by chemotherapy in both SCLC-A cell lines and SCLC-N cell lines (Fig. 6C and D; Supplementary Fig. S14). Overexpression of GAGE2A or PAGE5 had no impact on cellular survival in the absence of chemotherapy. To investigate the requirement for expression of GAGE2A and/or PAGE5 for chemoresistance in vivo, xenografts were generated from the shPAGE5, shGAGE2A, or double-knockdown SCLC cells, and the resulting xenografts treated with chemotherapy. Knockdown of GAGE2A or PAGE5 had a modest impact on reducing growth of SCLC-A tumors not treated by chemotherapy (Fig. 6E). In SCLC-A xenografts treated with chemotherapy, knockdown of GAGE2A, PAGE5, or both had a robust increase in the duration of treatment response (Fig. 6F). In SCLC-N xenografts, which are generally more inherently chemoresistant (70), knockdown of GAGE2A, PAGE5, or both had an impact on growth in the absence of chemotherapy (Fig. 6G). In addition, PAGE5 knockdown or the double GAGE2A and PAGE5 knockdown increased sensitivity to chemotherapy for the SCLC-N tumors, albeit to a lesser degree than the SCLC-A type cells due to incomplete knockdown of these CTAs (Fig. 6H). Overall, CTA knockdown significantly increased treatment response in human SCLC xenografts (Fig. 6F and G).

SCLC is a devastating disease for which tumor heterogeneity has functional impacts for the treatment and course of the disease (20, 71, 72). Tumor heterogeneity over time and under chemotherapeutic pressure is not well understood, but single-cell technologies have emerged in recent years as a powerful tool for the understanding of tumor evolution. Here we use genetic barcode lineage tracing and subclone analysis to understand tumor dynamics in a mouse model of SCLC and human SCLC xenografts.

We performed scRNA-seq analysis on a cohort of RPR2 mice to investigate the genetic drivers of tumor development. We found the initial emergence of a population of neuroendocrine-high “early” tumor cells, which are maintained through the later stages of tumor development but make up proportionately less of the bulk of the tumor at later stages of tumor development. In their place, a “late” population of tumor cells arises and over later timepoints begins to be responsible for the majority of the tumor mass. Cells from chemoresistant tumors cluster most closely with this “late” tumor population. Members of the AP-1 network, especially Jun and Fos were highly upregulated in the late tumor cluster, and human SCLC cells in culture rely on the activation of the AP-1 network to maintain colony-forming abilities. We observe that chemotherapy-relapsed autochthonous SCLC tumors also favor this AP-1–positive phenotype. Whether this reflects a positive selection for AP-1–expressing cells by chemotherapy, or an induction of AP-1 activity posttreatment is unclear, but due our observation of the requirement for AP-1 in colony formation, these observations may favor the former hypothesis. In our SCLC xenograft model, members of the AP-1 network are expressed in tumor-initiating clones that also possess a more mesenchymal phenotype. JUN is capable of upregulating Vimentin and driving the EMT (73). Furthermore, this EMT is linked to metastasis (20, 74), which would be expected for a tumor-initiating population of cells and a characteristic feature of human SCLC. Finally, there are many reports of coregulation of the AP-1 network and MYC either directly or mediated through the JNK (75–77). AP-1 has been described as a regulator of cellular plasticity (78–80), and therefore it is likely that its activity can be unveiling other genetic networks. Therefore, AP-1 may be responsible for the maintenance of the tumor-initiating state and required for the formation of recurrent tumors after chemotherapy, perhaps in conjunction with MYC activity.

Using a human xenograft model from two SCLC subtypes, SCLC-A and SCLC-N, we were able to investigate mechanisms of tumor development and treatment resistance. The SCLC-A cells, with high expression of ASCL1 in a chemotherapy-naïve environment, displays a profound transcriptomic shift after chemotherapy to a more SCLC-N–like profile (70) and upregulation of GAGE2A and PAGE5. In contrast, the SCLC-N cell line, NCI-H82, was derived from a patient who has already exposed to chemotherapy (81); therefore, this cell line already displays a postchemotherapy phenotype complete with upregulation of GAGE2A and to a lesser extent, PAGE5. Further investigating the SCLC-A cells to better query the switch from a chemotherapy naïve to a postchemotherapy phenotype, we observe two bottleneck events based on barcode diversity—the first in the generation of xenografts from cell lines, and the second in generating a chemoresistant xenograft, indicating that there are indeed discreet subpopulations of cells in an SCLC tumor capable of driving tumor initiation and chemoresistance. Subclonal analysis revealed that these two steps are not equally associated with cellular plasticity. While plasticity is increased slightly upon xenograft formation, the primary stage of enhanced plasticity is in response to chemotherapy (20). Furthermore, it is noteworthy that we observe both expansion of clones found in the original xenograft after chemotherapy treatment, and novel subclones that arise after chemotherapy. This indicates that while enhanced plasticity of tumor cells can lead to chemoresistance, many chemoresistant clones may reside in the tumor prior to treatment and then become unmasked upon the selective pressure of chemotherapy. Therefore, the mechanisms of evolution in SCLC tumors may not be uniform.

We, like Gardner and colleagues (82), have observed the upregulation of CTAs in chemo-treated SCLC. In particular, GAGE2A and PAGE5 are increased after chemoresistance in SCLC-A tumors and have relatively high expression in the inherently resistant SCLC-N tumors. To investigate the impact of GAGE2A or PAGE5 on mediating chemoresistance in SCLC, we knocked down GAGE2A and PAGE5 in SCLC-A and SCLC-N cells and found increased sensitivity to cisplatin in culture. Knockdown of either GAGE2A, PAGE5, or both GAGE2A and PAGE5 led to increased sensitivity in vivo in SCLC-A and SCLC-N xenografts. Conversely, overexpression of GAGE2A or PAGE5 in SCLC cell lines led to increased resistance to chemotherapy in vitro. Finally, human SCLC biopsies show expression of GAGE2A and PAGE5 via immunostaining. PAGE5 has been identified to be expressed in some cancers, and in melanoma was elevated as an antiapoptotic gene in response to platinum-based chemotherapy, where it was shown to be promote cell survival pathways (83). GAGE2A is another antiapoptotic CTA, that seems to be related to treatment resistance in medulloblastoma (84). Another CTA, CAGE, has been implicated in the regulation of upregulation of cyclin D and cyclin E, resulting in G₁ to S progression mediated by AP-1 (85). This cell-cycle regulation was mostly mediated by JUN and JUND, in line with our observation that these are the AP-1 members most highly expressed in SCLC TIC subpopulations; however, additional investigation will be required to determine whether there is a direct link between AP-1 activity and CTA expression. CTAs have shown promise as potentially targetable, unique cancer antigens. For this reason, they make excellent candidates for immunotherapy such as chimeric antigen receptor-T therapy and cancer vaccines (67–69, 86–89). Taken together, this indicates that the CTAs GAGE2A and PAGE5 are mediators of chemoresistance in SCLC and are expressed in human SCLC.

Initially, we set out to perform single-cell DNA barcoding to trace ITH during SCLC development in xenografts and an autochthonous mouse model. While our cellular labeling was robust in the xenograft model, it was greatly reduced in the Cas9-mediated autochthonous model due to the low efficiency of CRISPR recombination. Even still, our LBC detection levels were at best approximately 12% in the xenograft model and undetectable in the autochthonous model. We surmise that this is a result of the relatively low expression on GFP in the transduced cells or the autochthonous model which is driven by the Rosa26 locus expression. Future investigations may benefit from utilizing either a stronger promoter, or an inducible promoter to drive LBC expression. Furthermore, sequence capture techniques, such as the Feature Barcoding technology from 10X Genomics could be useful to enrich for LBC sequences. We did, however, demonstrate the PCR error is a considerable issue when performing genetic lineage barcoding as error introduced by PCR processing can artificially inflate the number of identified LBCs, and therefore confound downstream analyses of cellular heterogeneity by potentially subdividing unique clones. We therefore describe a framework where PCR error determined by the constant regions of an LBC can serve as threshold for any barcodes that appear at a lower frequency.

SCLC is a difficult disease to investigate with meager progress toward generating truly targeted therapeutics. The expansion of the single-cell clonal analysis has allowed us to examine ITH in SCLC with unprecedented resolution. Identifying the AP-1 network as being involved for SCLC tumor dynamics has provided knowledge of the critical early days of tumor formation in SCLC. Two potentially targetable antigens, GAGE2A and PAGE5 have been identified as direct mediators of chemoresistance in SCLC for the first time and may represent a new avenue for overcoming therapy resistance in SCLC in the future.

H. Wollenzien reports grants from NIH - NCI and National Science Foundation during the conduct of the study. Y. Afeworki Tecleab reports grants from NIH NIGMS Center for Pediatric Research grant #5P20GM103620 during the conduct of the study. M.S. Kareta reports grants from NIH and NSF during the conduct of the study. No disclosures were reported by the other authors.

H. Wollenzien: Conceptualization, resources, formal analysis, funding acquisition, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. Y. Afeworki Tecleab: Conceptualization, resources, software, formal analysis, investigation, methodology, writing–original draft, writing–review and editing. R. Szczepaniak-Sloane: Formal analysis, validation, investigation, visualization, writing–original draft, writing–review and editing. A. Restaino: Formal analysis, validation, investigation, visualization, writing–original draft, writing–review and editing. M.S. Kareta: Conceptualization, resources, data curation, software, formal analysis, supervision, funding acquisition, validation, investigation, visualization, methodology, writing–original draft, project administration, writing–review and editing.

We would like to acknowledge the NIH NIGMS Center for Cancer Research P20GM103548 for pilot award support and the NIH NCI/NIGMS grant R01CA233661 for research support (to M.S. Kareta). The NIH NCI grant 5F31CA243149 for training and research support (to H. Wollenzien). The Sanford Research Functional Genomics and Bioinformatics Core is supported by the NIGMS Center for Pediatric Research 5P20GM103620, and the Pathology and Flow Cytometry Cores are supported by the NIGMS Center for Cancer Research P20GM103548. H. Wollenzien is thankful to the University of South Dakota-Neuroscience, Nanotechnology and Networks (USD-N3) Program, supported by a grant from the National Science Foundation Research Traineeship program DGE-1633213.

The publication costs of this article were defrayed in part by the payment of publication fees. Therefore, and solely to indicate this fact, this article is hereby marked “advertisement” in accordance with 18 USC section 1734.