Introduction
Lung adenocarcinoma (LUAD), the most common cause of cancer death worldwide, exhibits significant heterogeneity in tumor cell identity and overall differentiation state (Travis et al., 2011). The state of LUAD differentiation correlates closely with prognosis, intrinsic therapeutic sensitivity, and drug resistance (Campos-Parra et al., 2014; Rotow and Bivona, 2017; Russell et al., 2011). Recent work by our lab and others has shown that the pulmonary lineage specifier NKX2-1/TTF1 is a central regulator of LUAD growth and identity (Maeda et al., 2012; Snyder et al., 2013). NKX2-1 is expressed in ~75% of human LUAD, and NKX2-1 negative tumors confer a worse prognosis than NKX2-1-positive tumors (Barletta et al., 2009). The specific role of NKX2-1 in LUAD depends, in part, on the driver oncogene (Maeda et al., 2012; Skoulidis and Heymach, 2019; Snyder et al., 2013). The majority of LUADs harbor mutually exclusive mutations in driver oncogenes that signal through the mitogen-activated protein kinase (MAPK) pathway, including
Treatment-naive
MAPK signaling is generally considered to drive proliferation and survival in LUAD, and has become a therapeutic target in some contexts (Ferrara et al., 2020). For example, small molecule inhibitors of BRAF and its downstream kinase MEK have recently shown clinical efficacy in
We have sought to dissect the role of NKX2-1 in
Results
Activation of oncogenic BRAF in the absence of NKX2-1 from the pulmonary alveolar epithelium leads to the development of invasive mucinous adenocarcinoma
To dissect the role of NKX2-1 in mutant BRAF-induced lung adenocarcinoma, we utilized two established mouse strains bearing recombinase-activatable alleles of
Figure 1.
(A) Representative photomicrographs of lung neoplasia arising 8 weeks after initiation with PGK-Cre lentivirus (5 × 103 pfu/mouse). BP mice are
Figure 1—figure supplement 1.
(A, B) Histologic phenotype of lung neoplasia arising 8 weeks after initiation with PGK-Cre lentivirus (5 × 103 pfu/mouse) in
Figure 1—figure supplement 2.
Percentage of incomplete recombinant tumors in
Graph represents mean ± S.D.
Despite similar morphologic phenotypes in the BRAFV600E and KRASG12D models,
In the KRASG12D model, we previously showed that inducing loss of lineage specifiers in established neoplasms can have distinct consequences when compared to lineage specifier loss at the time of tumor initiation (Camolotto et al., 2018). Models that enable gene manipulation in established tumors may be more physiologically relevant to tumor progression than models in which all genetic perturbations are present during tumor initiation. We therefore used a dual recombinase system in the present study to assess the role of NKX2-1 in the established BRAFV600E lung adenocarcinomas. Mice harboring FlpO-recombinase-sensitive alleles of
Figure 2.
(A–C) All mice are
Figure 2—figure supplement 1.
(A) Immunostaining for indicated mucin/gastric markers in NKX2-1-positive and -negative LUADs. Lung tumor sections were obtained from corn oil or tamoxifen treated
Figure 2—figure supplement 2.
Percentage of incomplete recombinant tumors in
Graph represents mean ± S.D.
We next evaluated the effect of NKX2-1 loss on the growth of established BRAFV600E -driven tumors. At early timepoints (3 weeks after tamoxifen treatment),
Loss of NKX2-1 stimulates MAPK signaling downstream of mutant BRAF
We have previously shown that
Figure 3.
NKX2-1 regulates MAPK pathway in BRAFV600E-driven LUAD.
(A, B) IHC for indicated proteins in lung neoplasia present at 8 weeks post-initiation in
Figure 3—figure supplement 1.
NKX2-1 regulates MAPK pathway in BRAFV600E-driven LUAD.
(A) IHC for indicated phospho-proteins in NKX2-1-positive and -negative LUADs. Lung tumor sections were obtained from corn oil or tamoxifen-treated
Finally, we performed IHC for pERK on a panel of mucinous (n = 17, see Materials and methods for inclusion criteria) and NKX2-1-positive primary human lung adenocarcinomas (n = 51). Strong pERK staining was detectable in the majority of mucinous tumors (14/17). Among positive cases, seven tumors exhibited diffuse staining (>90% of tumor cells positive), and another seven exhibited partial positivity (20–90% of tumor cells positive). In contrast, NKX2-1-positive LUAD showed a trend toward lower pERK staining than the mucinous adenocarcinomas. Only 20% of NKX2-1-positive tumors were diffusely positive (vs. 41% of mucinous tumors), and 35% of NKX2-1-positive tumors were negative for pERK (vs. 18% of mucinous tumors). Quantitation and representative pictures are shown in Figure 3D, Figure 3—figure supplement 1D. The ERK phosphorylation site is more labile than epitopes recognized by antibodies that detect total protein levels. Therefore, staining patterns likely reflect both biologic heterogeneity and variable processing of clinical specimens. Nevertheless, these data show that the high levels of pERK in IMA mouse models can be observed in many cases of human IMA.
Response of NKX2-1-positive and NKX2-1-negative tumors to RAF/MEK inhibition
The increase in MAPK activity we observed in both in vivo and in vitro models led us to ask whether NKX2-1 status influenced response to MAPK pathway inhibition. To investigate this, we directly compared the effect of RAF/MEK inhibitors on the growth and proliferation of autochthonous BP and BPN lung tumors. We chose dual inhibition of MEK and BRAF because this combination is the standard of care for BRAFV600E-mutant human lung adenocarcinoma (Planchard et al., 2016) and because in other genetically engineered mouse models driven by BRAFV600E example for thyroid cancer (McFadden et al., 2014), the combination of the two inhibitors is more effective than either one alone at inhibiting the MAPK pathway in vivo. We found that administration of BRAFV600E inhibitor (PLX4720) in combination with MEK inhibitor (PD0325901) effectively suppressed ERK phosphorylation in lung tumors of both genotypes (Figure 4A). Combined BRAF/MEK inhibition led to a dramatically lower tumor burden in both BP and BPN mice when assessed after 2 weeks (Figure 4B) and 4 weeks (Figure 4C) of MAPKi treatment. The decrease in tumor burden in BPN mice was greater by two-fold compared to that in BP mice, suggesting that BPN tumors were overall more sensitive to MAPK inhibition. In fact, we noted that residual MAPK-inhibited BPN tumors resembled hyperplasia induced by
Figure 4.
NKX2-1 status modulates response to MAPK pathway inhibitors.
(A) Representative H and E and phospho-ERK1/2 immunostaining photomicrographs of paraffin-embedded lung sections from
Figure 4—figure supplement 1.
NKX2-1 status modulates response to MAPK pathway inhibitors.
(A) Representative histology of BP and BPN tumors after 4 weeks of combined BRAFV600E and MEK inhibitor chow treatment, as well as at different timepoints following drug withdrawal demonstrating rapid tumor relapse. (B) Graphs comparing absolute expression levels of direct NKX2-1 targets (
To gain insights into the mechanisms of drug response in each model, we performed transcriptome profiling by sequencing RNA from whole tumors and single-cell RNA sequencing. Using the Cre-activated
We first performed Principal Component Analysis (PCA) on RNA sequencing data from whole tumors using the top 500 most variable genes. This analysis revealed that replicate samples cluster together but that the four groups (control and MAPKi-treated BP and BPN tumors) could be distinguished from each other at the transcriptomic level. As represented in Figure 4D, the greatest source of transcriptome diversity was related to deletion of
Despite these differences, NKX2-1 targets (e.g.
We characterized the transcriptome of 5065 control and 5563 MAPKi drug-treated single BPN cells (Figure 4—figure supplement 1E). In our initial analysis, rare non-tumor cells clustered separately from tumor cells, and their identity was further validated by the expression of well-known stromal markers including
Figure 5.
Differential impact of MAPK inhibitors on the cell cycle in NKX2-1-positive and NKX2-1-negative tumors.
(A, B) Immunostaining for and quantitation of the proliferation marker MCM2 in tumors from BP and BPN mice at 8 weeks post-initiation. Control or MAPK-inhibitor infused chow feeding started at 6 weeks and was maintained for 2 weeks. BP C (25 tumors from four mice), BP Tx (15 tumors from three mice), BPN C (20 tumors from three mice), BPN Tx (20 tumors from three mice). Scale bar: 100 µm. Graphs represent mean ± S.D. ****p<0.0001 by Student’s
Figure 5—figure supplement 1.
Differential impact of MAPK inhibitors on the cell cycle in NKX2-1-positive and NKX2-1-negative tumors.
(A) Detection of proliferation markers, BrdU incorporation and phospho-histone H3 Ser10, by IHCs in neoplasia arising 8 weeks after tumor initiation with PGK-Cre lentivirus (5 × 103 pfu/mouse). At the 6 week timepoint, BP and BPN mice were given control or PLX4720 plus PD0325901 inhibitor chow for 2 weeks. Scale bar: 100 µm. (B) Quantitation of BrdU-positive tumor cells from mice treated as in (A). Graphs represent mean ± S.D. BP C (34 tumors from three mice), BP Tx (23 tumors from three mice), BPN C (39 tumors from three mice), BPN Tx (40 tumors from three mice) were analyzed. ****p<0.0001 by Student’s
In contrast, the remaining cells (30%) exhibited frequent and abundant
NKX2-1 status influences cell cycle response to RAF/MEK inhibition
The MAPK pathway regulates proliferation by multiple mechanisms, including activating expression of the D-type cyclins, which can drive cells out of quiescence and into the cell division cycle (Lavoie et al., 1996; Tuveson et al., 2004). We therefore evaluated the impact of MAPKi drug treatment on cell cycle status of BP and BPN lung tumors using IHC for cell cycle markers as well as analysis of RNA-seq datasets. We first evaluated drug response by IHC for MCM2, a helicase detectable throughout the cell cycle but not in quiescence (G0). In BP tumors, 2 weeks of MAPKi treatment led to a significant decline in the percentage of MCM2-positive cells (Figure 5A,B). Additional IHC analysis showed a decline in the percentage of BP cells positive for BrdU incorporation and phospho-histone H3 (pHH3, M phase marker) (Figure 5—figure supplement 1A–C). Taken together, these data show that MAPKi treatment of BP tumors blocks proliferation and induces cell cycle exit, thereby increasing the proportion of cells in quiescence.
In contrast to BP tumors, MAPKi treatment of BPN tumors led to a paradoxical increase in the percentage of MCM2-positive cells (Figure 5A,B) despite decreased tumor burden (Figure 4B). In contrast, the percentage of BPN cells positive for BrdU and pHH3 declined to the same extent as BP cells after MAPKi treatment (Figure 5—figure supplement 1A–C). These data suggest that even though BPN cell proliferation is impaired by BRAF/MEK treatment, most residual BPN cells fail to exit the cell cycle and enter quiescence in the absence of MAPK activity.
Analysis of RNA-seq data from whole tumors further demonstrated that BP and BPN tumor cells exhibit a differential cell cycle response to MAPK inhibition. Drug-induced transcriptomic changes were highly distinct between genotypes. Pathway analysis (via Illumina Correlation Engine) demonstrated key differences between BP and BPN MAPK-inhibitor treated lung tumors in the most significant Gene ontology (GO) terms. Specifically, multiple cell cycle-related pathways decline in treated BP tumors, but not in treated BPN tumors, relative to their respective untreated controls (Figure 5C; Supplementary file 5, 6). Instead, the top GO pathway terms in MAPK-inhibited BPN tumors relative to vehicle controls included pathways related extracellular reorganization and cell motility (Figure 5D; Supplementary file 6).
We calculated cell cycle scores of control and MAPKi-treated tumor cells in scRNA-seq data using the Seurat package and phased cell cycle gene signatures previously defined in mouse cells (Mizuno et al., 2009). The results of this analysis are consistent with the general conclusion that BP and BPN cells have distinct cell cycle responses to RAF/MEK inhibition (Figure 5E). However, the Seurat cell cycle score approach is not well adapted to distinguish G0 from G1 phase cells. We therefore developed a novel methodology for analyzing quiescence in scRNA-seq data using diffusion mapping (Haghverdi et al., 2015) and the complete set of fine-scale cell cycle phase signatures from Mizuno et al., 2009, which includes both quiescent and cycling cells (Figure 5—figure supplement 1D). Using this novel approach, cells with high expression of S-phase or G2/M-phase signatures mapped to successive extremes at the right of the cyclic graph (Figure 5F), while the majority of cells (>80%) mapped to the left-hand portion corresponding to G0 and G1 signature enrichments (Figure 5F). Remarkably, these scRNA-seq analyses also revealed NKX2-1-dependent effects of drug treatment on cell cycle status. Treatment of BP cells caused a marked redistribution toward the left hand extreme of the graph, which is most distal from cycling cells and is enriched for the G0 signature, while the distribution of BPN cells was minimally affected by MAPK-inhibition (Figure 5G).
Given the ability of the MAPK pathway to drive cell cycle entry by activation of D-type cyclin expression, we evaluated the expression patterns of
These data raise the possibility that increased Cyclin D2 levels in BPN tumors upon MAPK inhibition might maintain CDK4/6 activity, thereby preventing cells from entering quiescence despite their MAPK-low state. To test this idea further, we used IHC to evaluate RB phosphorylation at S807/S811, two sites that can be phosphorylated by CDK4/6 (Figure 5I). Two weeks of combined BRAF/MEK inhibition significantly decreased the percent of phospho-RB-positive cells in BP tumors (Figure 5I). In contrast, MAPK-inhibition caused no change in phospho-RB S807/S811 levels in BPN tumors.
To evaluate specifically the potential role of Cyclin D2 in the cell cycle response to MAPK inhibitors, we knocked down Cyclin D2 in BPN organoids using two different shRNAs (Figure 5—figure supplement 1H). We then transplanted these lines into NSG mice subcutaneously, allowed tumors to form and grow for 8 weeks, and treated mice with RAF/MEK inhibitor for 4 weeks. Histopathologic analysis revealed that organoid-derived tumors exhibited a triphasic morphology in vivo, including large cystic structures, glandular adenocarcinoma, and small foci of poorly differentiated adenocarcinoma. Cystic structures appeared to accumulate fluid at a variable rate, precluding the use of tumor size as a metric of drug response. We quantitated MCM2 levels in the glandular adenocarcinomas, the component with the greatest similarity to the autochthonous model. We found that MAPK inhibition failed to reduce the MCM2 rate in control adenocarcinomas, whereas the MCM2 rate declined in tumors with Cyclin D2 knockdown upon RAF/MEK inhibitor treatment (Figure 5J). Although organoid-derived subcutaneous tumors do not perfectly mimic the autochthonous BPN model, these data provide additional orthogonal evidence that Cyclin D2 plays a role in the response of NKX2-1-negative lung adenocarcinoma to MAPK inhibition.
Taken together, these data show that NKX2-1 status has a profound effect on the cell cycle response of BRAFV600E-driven lung tumors to targeted therapy. We propose that Cyclin D2 helps maintain CDK4/6 activity and RB phosphorylation in MAPKi-treated BPN tumors, thus preventing cell cycle exit. The inability of residual BPN cells to exit the cell cycle suggests that distinct therapeutic strategies may be needed to eliminate BP versus BPN tumor cells after MAPK inhibition.
MAPK pathway regulates identity of NKX2-1-negative mucinous adenocarcinoma within the gastric lineage
Further analysis of gene expression changes in BPN tumors revealed surprising switch-like changes in the expression of lineage markers associated with specific gastric cell types. BPN tumor cells normally express high levels of transcripts that mark surface mucous cells in the stomach, such as
Figure 6.
MAPK activity regulates cellular identity in NKX2-1-negative tumors.
(A) Transcriptome analysis comparing control BPN cells with MAPKi-treated BPN tumor cells. Heatmap depicts gastric surface mucous cell markers whose expression diminished versus gastric chief cell or tuft cell markers whose expression increased in MAPK-inhibitor-treated BPN cells relative to untreated.
Figure 6—figure supplement 1.
MAPK activity regulates cellular identity in NKX2-1-negative tumors.
(A) IHC for detection of the surface neck mucous marker, GKN1, and Alcian Blue staining in control BPNs and BPN tumors that were treated with drug chow for 2 weeks. Scale bar: 100 µm. (B) IHC for the tuft cell marker, POU2F3, in human non-mucinous lung adenocarcinoma and IMA tumors. Scale bar: 100 µm. (C) Percentage of NKX2-1-positive and NKX2-1-negative human lung adenocarcinomas that contain POU2F3-positive tumor cells. (D) Immunostaining for tdTomato on BPN and KPN organoid FFPE sections indicating that primary spheroid cultures contain only tumor-derived epithelial cells. (E-G) qRT-PCR determination of the expression levels of select stomach epithelium lineage markers in multiple lung tumor organoid lines. Spheroid cultures were established from
Figure 6—figure supplement 2.
Durability of BRAF/MEK inhibitor-induced cell identity changes in NKX2-1-negative tumors.
IHC for indicated markers, plus Alcian Blue staining, on BPN tumors from mice treated with BRAF/MEK inhibitor chow for 1 month. Tumors were harvested 2 weeks after drug cessation. Scale bar: 100 µm.
A shift in expression of gastric lineage markers was also evident at the single-cell level. Our single-cell transcriptomics confirmed that gastric surface mucous markers like
Given this apparent lineage switch, we used IHC to evaluate a subset of lineage markers at the protein level in vivo. GKN1 protein (surface mucous marker) was abundant in control BPN tumors, but was depleted in MAPK inhibitor-treated BPN tumors (Figure 6—figure supplement 1A). Alcian Blue staining for mucin production showed the same pattern as GKN1 protein. Control BPN tumors were entirely negative for Pepsinogen C, whereas a subset of drug-treated BPN cells were Pepsinogen C-positive (Figure 6D, upper row). In contrast, a small subset of control BPN tumor cells were POU2F3 positive, and the fraction of positive cells was higher upon MAPK inhibition (Figure 6D, lower row). These data are concordant with our single-cell analysis of these two markers in terms of changes in proportion of positive cells vs. induction of absolute levels on a per-cell basis. Both Pepsinogen C and POU2F3-positive cells were found in alveolar hyperplasia induced by
We also examined POU2F3 by IHC in a panel of primary human LUAD tissues (Figure 6—figure supplement 1B). We found that almost all IMA tumors (15/16) harbor a minority population of POU2F3-positive cells (Figure 6—figure supplement 1B). Like in BPN tumors, these POU2F3-positive cells were relatively rare (~5% or less of tumor cells overall). In contrast, most NKX2-1-positive human lung adenocarcinomas were entirely POU2F3-negative (44/51) (Figure 6—figure supplement 1B,C). Interestingly, we detected a minor population of POU2F3-positive tumor cells in seven NKX2-1-positive cases (Figure 6—figure supplement 1B). This expands upon other recent observations that POU2F3 can be upregulated in specific subtypes of lung cancer (Huang et al., 2018). Thus, a rare population of POU2F3-positive tumors cells is readily detectable in both human IMA and murine models.
We next asked whether regulation of gastric lineage by the MAPK pathway is a general feature of IMA, or limited to BRAFV600E-driven tumors. To address this question, we generated organoids from autochthonous IMA models driven by KRASG12D (KPN) or BRAFV600E (BPN). IHC for the tdTomato reporter on sections of primary organoid cultures (Figure 6—figure supplement 1D) showed that all cultures examined are ~90–100% positive for tdTomato. Rare tdTomato-negative organoids have an exclusively epithelial morphology, and thus are unlikely to represent stromal cell contamination in these cultures. We treated BPN and KPN organoids with a MEK inhibitor (Cobimetinib) for 3 days and found that this was sufficient to stimulate the expression of chief (
MAPK inhibition leads to an increase in canonical WNT activity in NKX2-1 negative lung adenocarcinoma
LGR5 marks a subset of normal murine chief cells that can function as facultative stem cells (Leushacke et al., 2017).
We examined a set of WNT pathway genes derived from the ‘WNT signaling’ ontology category on AmiGO (Carbon et al., 2009) in our whole-tumor RNA sequencing of BP and BPN tumors. Several of these genes have previously been identified as direct transcriptional targets of canonical WNT signaling, including
Figure 7.
MAPK inhibition activates WNT signaling in NKX2-1-negative tumors.
(A) Transcriptome analysis of genes comprising the canonical WNT pathway gene ontology (AmiGO) in RNA purified from FACS-sorted BPN tumor cells 1 week post- treatment with PLX4720+PD0325901 or control chow.
Figure 7—figure supplement 1.
MAPK inhibition activates WNT signaling in NKX2-1-negative tumors.
(A) List of the top-scoring genes in IPA Upstream Regulator analysis when comparing differentially expressed genes between control- and BRAFi/MEKi-chow treated BPN tumors in whole-tumor RNA sequencing data. (B) Regulator Motif analysis by Illumina Correlation Engine identifies significant differences in TCF1 binding site genesets between control and treated BPN tumors. (C, D) Plots of normalized raw and imputed expression values for
Figure 7—figure supplement 2.
Multiple sources for WNT ligands in BRAFV600E-driven lung adenocarcinoma.
Relative expression values for all
Intriguingly, the pattern of induction was distinct for these two WNT target genes (Figure 7—figure supplement 1C,D). Cells with the highest levels of
To further investigate this question, we performed in vitro time course experiments, which showed that MEK inhibition induces
Previous studies in other cancer types have documented negative feedback loops between the MAPK and WNT signaling pathways. Mechanistically, these have been reported to be mediated by changes in AXIN1 (Zhan et al., 2019), TCF4 isoform levels (Heuberger et al., 2014) or phosphorylation of FAK (Chen et al., 2018) after MAPK inhibition. However, we have not been able to detect consistent changes in any of these parameters after treatment of our organoids with MEK inhibitors (Figure 7—figure supplement 1E).
To identify potential sources of WNT production, we analyzed the levels of
WNT signaling and FoxA1/2 contribute to lineage switching and cell cycle response in NKX2-1-negative lung adenocarcinoma treated with RAF/MEK inhibitor
WNT/β-catenin signaling plays an essential role in gastric epithelial patterning during embryogenesis and promotes chief cell differentiation (McCracken et al., 2017). We therefore asked whether the WNT pathway played a role in the lineage switch induced by MAPK inhibition in IMA. We used two different approaches to modulate WNT signaling in organoid cultures. First, we cultured organoids in 5% conditioned media (L-WRN), corresponding to a 10 fold reduction in the amount of exogenous WNT3A and R-spondin three in the culture media relative to standard, 50% L-WRN, conditions. Second, we used the small molecule LGK-974, which blocks Porcupine-mediated posttranslational modification of WNT ligands that is required for their secretion and paracrine signaling. Levels of
Figure 8.
WNT signaling and the transcription factors FoxA1/FoxA2 are partially required for lineage switching induced by MAPK inhibition.
(A) Analyses of indicated gene expression levels in BPN (1243 and 988) and KPN (1268) tumor organoid lines by qRT-PCR at 24 hr and under different treatment conditions. Organoids were cultured in 50% L-WRN (a, b) or reduced 5% L-WRN media (c, d, e, f) and treated with DMSO (a, c), single agent Cobimetinib (b, d) or the Porcupine inhibitor, LGK-974 (e), and both inhibitors (f). Graphs indicate mean ± S.D. [p values are for
Figure 8—figure supplement 1.
WNT signaling and the transcription factors FoxA1/FoxA2 are partially required for lineage switching induced by MAPK inhibition.
(A) Analyses of indicated gene expression levels in BPN (1243 and 988) and KPN (1268) tumor organoid lines by qRT-PCR at 24 hr and under different treatment conditions. Organoids were cultured in 50% L-WRN (a, b) or reduced 5% L-WRN media (c, d, e, f) and treated with DMSO (a, c), single agent Cobimetinib (b, d) or the Porcupine inhibitor, LGK-974 (e), and both inhibitors (f). Graphs indicate mean±S.D. ****p<0.0001, ***p<0.001, **p<0.01, *p<0.05, ns = not significant by Student’s
Figure 8—figure supplement 2.
Characterization of MAPK inhibitor-driven signaling changes.
(A, B) Active (unphosphorylated on T41, S37, S33 residues) and total β-catenin staining in tumors under indicated treatment conditions. For BPN tumors, see Figure 8C for experimental details. BP tumors are 8-week autochthonous lung tumors from mice that received control or MAPK-inhibitor chow for 2 weeks (See Figure 1A for more experimental details). Small intestine: positive control for detection of active β-catenin. Scale bar: 100 µm.
Previously, our lab reported that FoxA1 and FoxA2 are required for gastric differentiation in mouse models of IMA, including surface mucous cell marker expression (Camolotto et al., 2018). Here, using a novel NKX2-1-negative organoid line that harbors conditional alleles of
Finally, we asked whether WNT signaling might play a role in failure of BPN tumors to exit the cell cycle after RAF/MEK inhibition. WNT signaling contributes to oncogenesis in a variety of settings (Zhan et al., 2017), including in NKX2-1-positive LUAD GEMMs driven by BRAFV600E (Juan et al., 2014) or KRASG12D (Tammela et al., 2017). Further, β-catenin activation upon
To investigate this possibility, we first correlated known signatures of MEK (Dry et al., 2010) and WNT activity (‘positive regulation of canonical WNT signaling’ category on AmiGO, Carbon et al., 2009) with levels of
Given the relationship between WNT activity and
We also assessed the status of β-catenin in drug-treated BPN tumors in vivo. Using an antibody that recognizes non-phosphorylated (active) β-catenin or a second antibody for total β-catenin (Figure 8—figure supplement 2A,B), we found that vehicle and LGK-974-treated tumors exhibit predominantly membranous staining by IHC, with no evidence of nuclear β-catenin. MAPK inhibition alone elicited accumulation of β-catenin throughout the cell, including the nucleus. Levels of nuclear/cytoplasmic β-catenin in tumors diminished in the presence of combined MAPK and LGK-974. These findings are consistent with the possibility that β-catenin mediates MAPK-inhibition-induced WNT pathway activation and lineage switching in BPN tumors. In contrast, control BP tumors had lower levels of active β-catenin than BPN tumors, and there was no appreciable increase in staining with MAPK inhibition.
Finally, we performed a larger experiment in which we combined RAF/MEK inhibition with either WNT inhibition (LGK-974) or CDK4/6 inhibition (Palbociclib) for 4 weeks. We analyzed a subset of mice at the end of 4 weeks of drug treatment and performed a survival study on the remaining mice. In mice analyzed immediately after the final dose of drug, we found that adding either LGK-974 or Palbociclib to RAF/MEK inhibitors increased the number of cells in quiescence compared to RAF/MEK inhibitor alone (Figure 8D). Addition of Palbociclib to the RAF/MEK inhibitor also reduced the amount of phospho-RB (S807/811) in the tumors (Figure 8E). These data suggest a model in which elevated WNT signaling in MAPKi-treated BPN tumors prevents cell cycle exit by maintaining higher levels of CDK4/6 activity than in BP tumors treated with the same drug. Our gene expression data suggest that Cyclin D2 is the CDK4/6 partner that is most likely to maintain its activity in MAPKi-treated BPN tumors.
Despite the short term effects of adding a third drug, we did not see a significant increase in survival when either LGK-974 or Palbociclib was added to RAF/MEK inhibitor (Figure 8—figure supplement 1E). This indicates that driving residual BPN cells out of the cell cycle is not sufficient to prevent tumor rebound after drug cessation.
Discussion
LUAD progression is driven not only by activation of MAPK signaling but also by changes in cellular identity. However, the direct impact of dysregulated MAPK activity on LUAD identity itself remains an unexplored problem in lung cancer. This problem is important not only because of the correlation between cellular identity and intrinsic malignant potential, but also because lineage switching has become an increasingly recognized, common mechanism of resistance to therapies targeting the MAPK pathway. A better understanding of this phenomenon would provide new insights into the natural history of LUAD progression as well as the changes in cellular identity that occur as a result of targeted therapy.
Analyses of the histologic spectrum of LUAD in patients and mouse models have identified two cellular differentiation programs that impose barriers to lung tumor progression (Camolotto et al., 2018; Snyder et al., 2013; Winslow et al., 2011). Initially, NKX2-1/FOXA1/2 transcriptional networks maintain a well-differentiated pulmonary identity. We have previously shown that downregulation of NKX2-1 and subsequent relocalization of FOXA1 and FOXA2 to the regulatory elements of gastrointestinal lineage genes allows for loss of surfactant proteins, aberrant upregulation of HNF4A and conversion to gastric cell identity. Loss of both differentiation programs can lead to poorly differentiated tumors containing high levels of HMGA2.
Human IMAs are associated with adverse clinical outcomes and harbor a distinct spectrum of driver mutations. Compared to LUAD overall,
GEMMs have been valuable in explaining the distinct role of NKX2-1 in LUAD driven by
We have previously shown that the gastric differentiation state of IMA is driven by transcription factors such as FoxA1/2 and HNF4α. Here, we show that MAPK and WNT signaling provide an additional layer of regulation, modulating the specific cell type that IMA cells most closely resemble within the gastric lineage. Our data reveal that high MAPK activity drives a gene expression program characteristic of the surface mucous cells of the stomach. In contrast, low ERK levels, in concert with high WNT activity, activate gene expression signatures associated with other gastric cell types, including LGR5-positive chief cells and tuft cells. POU2F3-positive tuft-like cells have also been identified in pancreatic neoplasia, where they restrain tumorigenesis (DelGiorno et al., 2020). Interestingly, hyperactivation of MAPK signaling (via TGF-α mediated activation of EGFR) is observed in Menetrier’s disease, a hypertrophic gastropathy characterized by hyperproliferation of isthmus progenitors and the preferential differentiation of these progenitors into surface mucous cells at the expense of parietal and chief cells (Fiske et al., 2009). Blockade of EGFR results in regression of overgrown surface epithelium and restoration of lineage fidelity in progenitors.
Importantly, our results identify crosstalk between the MAPK and WNT pathways in IMA tumors, as suppression of MAPK pathway rapidly activates WNT target gene expression. In the mammalian intestine, reciprocal gradients of WNT and phosphorylated ERK1/2 control the balance between proliferation and differentiation of intestinal stem cells (Kabiri et al., 2018). Some molecular mechanisms driving the crosstalk between WNT and ERK signaling have been reported. In one study, reduction of MEK activity via gut-specific ablation of
Our results also have translational implications for LUAD treatment.
Additional work is clearly needed to identify specific vulnerabilities of residual IMA cells after MAPK inhibition. The lineage switch we observe in both KRAS and BRAF-driven IMA models is likely to be relevant to the identification of such vulnerabilities. MAPK inhibition causes IMA cells to undergo an identity shift from a surface mucous cell-like phenotype toward a more gastric stem cell-like state or tuft cell identity. These observations support the broader notion that tumor cell plasticity includes the capacity to explore myriad cell states from the developmental repertoire. Lineage switching may be particularly beneficial for tumor cells under the selective pressure of targeted therapy, as diverse lineage programs have distinct dependencies on specific growth and survival mechanisms. We speculate that cell identity changes within the gastric lineage may render residual IMA cells susceptible to specific therapeutic interventions. This would be reminiscent of prior studies of basal cell carcinoma of the skin showing that an identity switch confers resistance to inhibition of the Hedgehog pathway while rendering the cells susceptible to blockade of WNT signaling (Biehs et al., 2018; Sánchez-Danés et al., 2018).
In summary, this study characterizes the context-dependent role of NKX2-1 in LUAD and identifies novel mechanisms of cell identity regulation and therapeutic response in this disease (Figure 8F). The data provide new insights into the complex interplay between lineage specifiers, oncogenic signaling pathways and the susceptibility of lung cancer cells to targeted therapy. Thus, in the interest of developing fully effective therapies, these studies call for deep cataloguing of cell states sampled by tumor cells as well as their state-specific vulnerabilities.
Materials and methods
Key resources table
| Reagent type | Designation | Source or reference | Identifiers | Additional information |
|---|---|---|---|---|
| Genetic reagent ( |
| Dankort et al., 2007 | MGI:3711771 | Dr. Martin McMahon (HCI, Salt Lake City, Utah); mixed C57BL/6J × 129SvJ background |
| Genetic reagent ( |
| Shai et al., 2015 | Dr. Martin McMahon (HCI, Salt Lake City, Utah); mixed C57BL/6J × 129SvJ background | |
| Genetic reagent ( |
| Jonkers et al., 2001 | MGI:1931011 | Dr. Anton Berns, University of Amsterdam; Jackson Laboratories (Bar Harbor, Maine); mixed C57BL/6J × 129SvJ background |
| Genetic reagent ( |
| Lee et al., 2012 | MGI:5306612 | Dr. David G Kirsch (Duke University Medical Center, Durham, North Carolina); mixed C57BL/6J × 129SvJ background |
| Genetic reagent ( |
| Jackson et al., 2001 | MGI:2429948 | Dr. Tyler Jacks (MIT, |
| Genetic reagent ( |
| Young et al., 2011 | MGI:5007794 | Dr. Tyler Jacks (MIT, |
| Genetic reagent ( |
| Kusakabe et al., 2006 | MGI: 3653706 | Dr. Shioko Kimura (NCI, NIH, Bethseda, Maryland); mixed C57BL/6J × 129SvJ background |
| Genetic reagent ( | Madisen et al., 2010 | MGI: 4436847 | Jackson Laboratories (Bar Harbor, Maine); mixed C57BL/6J × 129SvJ background | |
| Genetic reagent ( |
| Schönhuber et al., 2014 | MGI: 5616874 | Dr. Dieter Saur (Technische Universitat Munchen, Munchen, Germany); mixed C57BL/6J × 129SvJ background |
| Genetic reagent ( |
| Gao et al., 2008 | MGI: 3831163 | Dr. Klaus H. Kaestner (Univ. of Pennsylvania School of Medicine, Philadelphia, PA); mixed C57BL/6J × 129SvJ background |
| Genetic reagent ( |
| Sund et al., 2000 | MGI: 2177357 | Dr. Klaus H. Kaestner (Univ. of Pennsylvania School of Medicine, Philadelphia, PA); mixed C57BL/6J × 129SvJ background |
| Genetic reagent ( | NOD/SCID-gamma chain deficient (NSG) | The Jackson Laboratory | 005557 | |
| Cell line | 293T | DuPage et al., 2009 | ||
| Cell line ( | L-WRN | ATCC | CRL-3276 | |
| Recombinant DNA reagent | d8.9 (plasmid) | DuPage et al., 2009 | ||
| Recombinant DNA reagent | VSV-G (plasmid) | DuPage et al., 2009 | ||
| Recombinant DNA reagent | 7TGP (plasmid) | Addgene | 24305 | |
| Chemical compound, drug | Tamoxifen | Sigma Aldrich | T5648 | |
| Chemical compound, drug | Tamoxifen supplemented chow | Envigo | TD.130858 | |
| Chemical compound, drug | PLX4720 supplemented chow | Plexxikon/ | Tsai et al., 2008 | |
| Chemical compound, drug | PD0325901 supplemented chow | Plexxikon/ | Trejo et al., 2012 | |
| Chemical compound, drug | PLX4720 | Selleckchem | S1152 | |
| Chemical compound, drug | PD0325901 | Selleckchem | S1036 | |
| Chemical compound, drug | GDC-0994 | Genentech | Blake et al., 2016 | |
| Chemical compound, drug | Cobimetinib | Genentech | ||
| Chemical compound, drug | Palbociclib | LC Laboratories | P-7744 | |
| Chemical compound, drug | LGK-974 | Selleckchem | S7143 | |
| Antibody | Anti-NKX2-1 (Rabbit monoclonal) | Abcam | Cat# ab76013 | (1:2000) |
| Antibody | Anti-DUSP6 (Rabbit monoclonal) | Abcam | Cat# ab76310 | (1:500) |
| Antibody | Anti-SPRY2 (Rabbit monoclonal) | Abcam | Cat# ab180527 | (1:1000) |
| Antibody | Anti-pERK (Rabbit monoclonal) | Cell Signaling Technology | Cat# 4370 | (1:500) |
| Antibody | Anti-ERK (Mouse monoclonal) | Cell Signaling Technology | Cat# 4696 | (1:2000) |
| Antibody | Anti-pRSK S380 (Rabbit monoclonal) | Cell Signaling Technology | Cat# 11989 | (1:300) |
| Antibody | Anti-p4EBP1 (Rabbit monoclonal) | Cell Signaling Technology | Cat# | (1:400) |
| Antibody | Anti-pS6 S235/236 (Rabbit monoclonal) | Cell Signaling Technology | Cat# 4858 | (1:400) |
| Antibody | Anti-pS6 S240/244 (Rabbit monoclonal) | Cell Signaling Technology | Cat# 5364 | (1:1000) |
| Antibody | Anti-pMEK S221 (Rabbit monoclonal) | Cell Signaling Technology | Cat# 2338 | (1:100) |
| Antibody | Anti-PGC (Rabbit polyclonal) | Cell Signaling Technology | Cat# HPA031718 | (1:100) |
| Antibody | Anti-POU2F3 (Rabbit polyclonal) | Cell Signaling Technology | Cat# HPA019652 | (1:200) |
| Antibody | Anti-pRB S807/S811 | Cell Signaling Technology | Cat# 8516 | (1:800) |
| Antibody | Anti-CD36 (Rat monoclonal) | R and D Systems | Cat# MAB25191 | (1:300) |
| Antibody | Anti-proSPB (Rabbit polyclonal) | Millipore | Cat# AB3430 | (1:3000) |
| Antibody | Anti-proSPC (Rabbit polyclonal) | Millipore | Cat# AB3786 | (1:4000) |
| Antibody | Anti-Gastrokine 1 (Mouse monoclonal) | Abnova | Cat# H00056287-M01 | (1:50) |
| Antibody | Anti-Muc5AC (Mouse monoclonal) | Abnova | Cat# MAB13117 | (1:100) |
| Antibody | Anti-Histone H3 phospho-Ser10 (Rabbit polyclonal) | Abcam | Cat# ab5176 | (1:200) |
| Antibody | Anti-HNF4A (Rabbit monoclonal) | Cell Signaling Technology | Cat# 3113 | (1:400) |
| Antibody | Anti-MCM2 (Rabbit monoclonal) | Abcam | Cat# ab108935 | (1:5000) |
| Antibody | Anti-PDX1 (Mouse monoclonal) | DSHB | Cat# F109-D12 | (1:20) |
| Antibody | Anti-RFP (Rabbit polyclonal) | Rockland Immunochemicals | Cat# | (1:1200) |
| Antibody | Anti-β-catenin (Mouse monoclonal) | BD Biosciences | Cat# 610153 | (1:200) |
| Antibody | Anti-Non-phospho β-catenin (Rabbit monoclonal) | Cell Signaling Technology | Cat# 8814 | (1:1600) |
| Antibody | Anti-Vinculin (Rabbit monoclonal) | Abcam | Cat# ab129002 | (1:20,000) |
| Antibody | Anti-β-tubulin (Mouse monoclonal) | DSHB | Cat# E7 | (1:15,000) |
| Antibody | Anti-Axin1 (Rabbit monoclonal) | Cell Signaling Technology | Cat# 2087 | (1:1000) |
| Antibody | Anti-TCF4 (Rabbit monoclonal) | Cell Signaling Technology | Cat# | (1:1000) |
| Antibody | Anti-pFAK Y397 (Rabbit polyclonal) | Sigma Aldrich | Cat# SAB4504403 | (1:1000) |
| Antibody | Anti-Galectin 4 (Goat polyclonal) | R and D Systems | Cat# AF2128 | (1:200) |
| Antibody | Anti-Cathepsin E (Rabbit polyclonal) | Lifespan Biosciences | Cat# LSB523 | (1:12,000) |
| Antibody | Anti-BrdU (Rat monoclonal) | Abcam | Cat# ab6326 | (1:400) |
| Antibody | Anti-CyclinD1 (Rabbit monoclonal) | Abcam | Cat# ab137875 | (1:100) |
| Antibody | Anti-CyclinD2 (Mouse monoclonal) | NeoMarkers | Cat# MS-213-P1ABX | (1:500) |
| Other | PGK-Cre | DuPage et al., 2009 | Lentivirus | |
| Other | Ad5CMVFlpo | Gene Transfer | VVC-U of | Adenovirus |
| Other | Ad5CMVCre | Gene Transfer | VVC-U of | Adenovirus |
| Other | Ad5SpcCre | Gene Transfer | VVC-Berns- | Adenovirus |
| Other |
| ThermoFisher Scientific | Mm00518185_m1 | TaqMan Gene Expression Assay (FAM) |
| Other |
| ThermoFisher Scientific | Mm00482488_m1 | TaqMan Gene Expression Assay (FAM) |
| Other |
| ThermoFisher Scientific | Mm01204823_m1 | TaqMan Gene Expression Assay (FAM) |
| Other |
| ThermoFisher Scientific | Mm00443610_m1 | TaqMan Gene Expression Assay (FAM) |
| Other |
| ThermoFisher Scientific | Mm00438890_m1 | TaqMan Gene Expression Assay (FAM) |
| Other |
| ThermoFisher Scientific | Mm00433596_m1 | TaqMan Gene Expression Assay (FAM) |
| Other |
| ThermoFisher Scientific | Mm02342429_g1 | TaqMan Gene Expression Assay (FAM) |
| Other |
| ThermoFisher Scientific | Mm00492318_m1 | TaqMan Gene Expression Assay (FAM) |
| Other |
| ThermoFisher Scientific | Mm00478293_m1 | TaqMan Gene Expression Assay (FAM) |
| Other |
| ThermoFisher Scientific | Mm00432359_m1 | TaqMan Gene Expression Assay (FAM) |
| Other |
| ThermoFisher Scientific | Mm00438070_m1 | TaqMan Gene Expression Assay (FAM) |
| Other |
| ThermoFisher Scientific | Mm00484713_m1 | TaqMan Gene Expression Assay (FAM) |
| Other |
| ThermoFisher Scientific | Mm01976556_s1 | TaqMan Gene Expression Assay (FAM) |
| Other | RBC Lysis Buffer | eBioscience | 00-4333-57 | |
| Other | Collagenase type I | Thermofisher Scientific | 17100017 | Enzyme |
| Other | Dispase | Corning | 354235 | Enzyme |
| Other | Deoxyribonuclease I | Sigma Aldrich | DN25 | Enzyme |
Mice, tumor initiation, and drug treatment in vivo
All animal work was done in accordance with a protocol approved by the University of Utah Institutional Animal Care and Use Committee. The following mouse strains were used.
Rodent Lab Diet (AIN-76A) was formulated with the vemurafenib-related compound PLX4720 at 200 mg/kg (Tsai et al., 2008) and PD0325901 7 mg/kg (Trejo et al., 2012). Drug formulation was by Plexxikon and chow manufacture was by Research Diets. AIN-76A was used as control chow. Mice were maintained on the indicated drug pellets for the indicated time periods. BrdU incorporation was performed by injecting mice at 50 mg/kg (Sigma) intraperitoneally 1 hr prior to necropsy. Mice in survival studies were monitored for lethargy and respiratory distress, at which time animals were euthanized.
LGK-974 and Palbociclib (Selleck Chemicals) were delivered to mice via oral gavage. LGK-974 was formulated in 0.5% (w/v) methylcellulose/0.5% (v/v) Tween-80 solution and given at 7.5 mg/kg dose with 5 days ON/2 days OFF schedule for 2 weeks (Figure 8C, S8D) or 4 weeks (Figure 8D and E, S8E). Weights were monitored once or twice weekly and no toxicity was observed in mice. Palbociclib was formulated in 50 mM Sodium Lactate and initially given at 120 mg/kg dose with 5 days ON/2 days OFF schedule. Weights were monitored once or twice weekly. Due to weight loss in female mice receiving Palbociclib + MAPK-inhibitor chow combination, the following modifications were made for all mice receiving Palbociclib as single or dual agent. On week 2: 120 mg/kg dose with 4 days ON/3 days OFF schedule; week 3: Palbociclib was skipped altogether; week 4: 100 mg/kg dose with 4 days ON/3 days OFF schedule.
Tamoxifen administration
Tumor-specific activation of CreERT2 nuclear activity was achieved by intraperitoneal injection of tamoxifen (Sigma) dissolved in corn oil at a dose of 120 mg/kg. Mice received a total of six injections over the course of 9 days. For survival experiments, mice were additionally given pellets supplemented with 500 mg/kg tamoxifen (Envigo).
Cell lines and cell culture
HEK-293T cells were cultured in DMEM/High Glucose medium (Gibco). L-WRN cells (ATCC) were cultured in DMEM/High Glucose medium containing G418 (Sigma) and Hygromycin (InvivoGen) media. All media contained 10% FBS (Sigma). To eliminate mycoplasma contamination, cell lines were treated with 25 µg/mL Plasmocin (InvivoGen) for 2–3 weeks. To maintain cultures mycoplasma free, media were supplemented with 2.5 µg/mL Plasmocin. Cell line identity was authenticated using STR analysis at the University of Utah DNA Sequencing Core.
Establishing primary lung organoids
Eight to 10 weeks after tumor initiation with PGK-Cre, tumor bearing mice were euthanized and lungs were isolated. Lungs were then minced under sterile conditions and digested at 37°C for 30 mins with continuous agitation in a solution containing the enzymes, Collagen Type I (Thermo Fisher Scientific, at 450 U/ml final); Dispase (Corning, at 5 U/ml); DNaseI (Sigma, at 0.25 mg/ml) and Antibiotic-Antimycotic solution (Gibco). The digested tissue was passed through an 18 or 20-gauge syringe needle. Enzyme reactions were then stopped by addition of cold DMEM/F-12 with 10% FBS and the cell suspension was dispersed through 100 µm, 70 µm, and 40 µm series of cell strainers to obtain single-cell suspension. Subsequently, cell pellets were treated with erythrocyte lysis buffer (eBioscience). Finally, cell pellets were reconstituted in Advanced DMEM/F-12 (Gibco) supplemented with L-glutamine, 10 mM HEPES, and Antibiotic-Antimycotic. Thereafter, 100,000 tumor cells were mixed with 50 µl of Matrigel (Corning) at 1 to 10 ratio and plated in 24-well plates. For the first week of organoid initiation, Matrigel droplets were overlaid with recombinant organoid medium (Advanced DMEM/F-12 supplemented with 1X B27 (Gibco), 1X N2 (Gibco), 1.25 mM nAcetylcysteine (Sigma), 10 mM Nicotinamide (Sigma), 10 nM Gastrin (Sigma) and containing growth factors (100 ng/ml EGF [Peprotech], 100 ng/ml R-spondin1 [Peprotech], 100 ng/ml Noggin [Peprotech], 100 ng/ml FGF10 [Peprotech], as well as the ROCK inhibitor Y27632 [R and D Systems] and the TGF-β type I receptor-inhibitor SB431542 [R and D Systems]) as described in Barker et al., 2010; Sato et al., 2009).
After the initial spheroid culture propagation, organoids were switched to 50% L-WRN conditioned media generated as detailed in Miyoshi and Stappenbeck, 2013. Briefly, L-WRN cells (ATCC) seeded and maintained in 10% FBS, DMEM high glucose containing 500 µg/ml G418 (Sigma) and 500 µg/ml Hygromycin (InvivoGen) media until confluency. Once confluent, L-WRN cells were switched to Advanced DMEM/F-12 containing 20% FBS, 2 mM L-glutamine, 100 U/ml penicillin, 0.1 mg/ml streptomycin. Daily, conditioned supernatant medium was collected from L-WRN cultures, filtered through 0.2 µm vacuum filters and saved at −20°C as stock (100%) L-WRN medium. To use for culturing spheroids, stock L-WRN was diluted with equal volume Advanced DMEM/F-12 (final concentration 50%). Where specified, stock L-WRN media was diluted with Advanced DMEM/F-12 to a final concentration of 5% while still keeping FBS concentration at 10% as described in VanDussen et al., 2015. For drug studies, organoid media was supplemented with Cobimetinib (GDC-0973), GDC-0994 (Genentech), LGK-974, PD0325901, PLX-4720 (Selleck Chemicals) at indicated concentrations. Drug containing media was refreshed every 48–72 hr.
Lentiviral production and transduction
HEK293T cells were transfected with CRE-encoding lentiviral vector, d8.9 packaging vector and VSV-G envelope vector mixed with TransIT-293 (Mirus Bio). Virus containing supernatant was collected at 36, 48, 60, and 72 hr after transfection. Ultracentrifugation at 25,000 r.p.m. for 2 hr was necessary to concentrate virus for in vivo infection (DuPage et al., 2009).
For stable transduction of organoids, organoid cultures were first prepared into single cell suspensions by subjecting them to successive incubations with Cell Recovery Solution (Corning) and TrypLE (Gibco). Cells were then resuspended in a 1:1 mixture, by volume, of 50% L-WRN and lentivirus containing supernatant. After supplementation with 8 µg/ml polybrene, cells were incubated for 2 hr. Cells were then pelleted, mixed back with Matrigel, and seeded. 72 hr later, Puromycin selection for 1 week was performed to achieve stable lines.
Flank tumor transplantation
For subcutaneous allograft experiments, 3 × 105 single-cell suspension of 988 BPN organoid cells stably expressing either the Mission pLKO.1-puro Non Target shRNA control (shscr) or one of the validated
Real-time PCR
Total RNA was isolated using TRIzol followed by PureLink RNA Mini Kit (Thermo Fisher Scientific). RNA was treated with RNase-Free DNAse I (Invitrogen) on-column. For low cell numbers (e.g. after FACS), Animal Tissue RNA Purification Kit (Norgen Biotek) was used instead. RNA was converted to cDNA using iScript Reverse Transcription Supermix (BioRad). cDNA was analyzed either by SYBR green real-time PCR with SsoAdvanced Universal SYBR Green Supermix (BioRad) or by Taqman real-time PCR with SsoAdvanced Universal Probes Supermix (BioRad) using a CFX384 Touch Real-Time PCR Detection System (BioRad). Gene expression was calculated relative to
Cell viability and growth assay
For organoids, cells were seeded at a density of 2000 cells per 15 µl Matrigel dome, four domes per line per time point, each dome in a single well of a 96-well plate. After overnight culture, organoids were treated with different concentrations of the indicated inhibitors and incubated at 37°C for indicated times. Subsequently, the relative number of viable cells was measured by by CellTiter-Glo 3D Cell Viability Assay (Promega), according to the manufacturer’s instructions. Luminescence was then measured by a microplate reader (EnVision 2105 Multimode Plate Reader, PerkinElmer). Replicate values for each experimental group were averaged and all values were normalized to control treatment group for each line.
Immunoblot analysis
Cells were lysed on ice using RIPA buffer (50 mM HEPES, 140 mM NaCl, 1 mM EDTA, 1% triton X-100, 0.1% sodium deoxycholate and 0.1% SDS) supplemented with protease and phosphatase inhibitor cocktails (Roche). Protein extracts were clarified and concentrations were measured with Pierce Coomassie Plus Protein Assay reagent (Thermo Fisher Scientific). Lysates were then resolved on SDS-PAGE gels (BioRad), and transferred to Nitrocellulose blots. Membranes were probed with primary antibodies against AXIN1 (2087, CST), TCF4 (2569, CST) pFAK Y397 (SAB4504403, Sigma), DUSP6 (EPR129Y, Abcam), SPRY2 (EPR4318(2)(B), Abcam), pERK (4370, CST), tERK (4696, CST), Vinculin (EPR8185, Abcam), and β-tubulin (E7, DSHB). Blots were subsequently incubated with HRP conjugated secondary antibodies. For signal development, Supersignal West Femto Substrate kit (Thermo Fisher Scientific) was used, followed by image acquisition using darkroom development. ImageJ was used for band intensity quantitation.
Histology and immunohistochemistry
Lung tissues were fixed overnight in 10% neutral buffered formalin, processed through 70% ethanol, embedded in paraffin, and sectioned at 4 µm thickness. Staining of hematoxylin and eosin, as well as detection of mucin by Alcian Blue were carried out at the HCI Research Histology Shared Resource. Immunohistochemistry was performed manually as detailed in Camolotto et al., 2018. The following primary antibodies were used. BrdU (BU1/75, Abcam), Cathepsin E (LS-B523, Lifespan Biosciences), Galectin 4 (AF2128, R and D Systems), Gastrokine 1 (2E5, Abnova), Histone H3 pSer10 (ab5176, Abcam), HNF4A (C11F12, CST), MCM2 (ab31159, Abcam), Muc5AC (SPM488, Abnova), NKX2-1 (EP1584Y, Abcam), PDX1 (F109-D12, DSHB), RFP (600-401-379, Rockland Immunochemicals), pERK1/2 (D13.14.4E, CST), pRSK S380 (D3H11, CST), p4EBP1 T37/46 (2855, CST), pS6 S235/236 (D57.2.2E, CST), pS6 S240/244 (D68F8, CST), pMEK1/2 S221 (166F8, CST), Pepsinogen C (HPA031718, Sigma), POU2F3 (HPA019652, Sigma), pRB S807/811 (8516, CST), Cyclin D1 (SP4, Abcam), Cyclin D2 (DCS-3.1, NeoMarkers), CD36 (324205, R and D Systems), pro-SPB (AB3430, Millipore), pro-SPC (AB3786, Millipore), β-catenin (610153, BD Biosciences), active β-catenin (D13A1, CST). To visualize Cyclin D2 signal, IHC was performed using the CSA II, Biotin-Free Catalyzed Amplification System (Dako) following manufacturer’s instructions. Pictures were taken on a Nikon Eclipse Ni-U microscope with a DS-Ri2 camera and NIS-Elements software. Tumor quantitation was performed on hematoxylin and eosin-stained or IHC-stained slides using NIS-Elements software. All histopathologic analyses were performed by a board-certified anatomic pathologist (E.L.S.).
Primary human tumors
De-identified formalin fixed, paraffin-embedded (FFPE) tumors were obtained from the Intermountain Biorepository, which collects samples in accordance with protocols approved by the Intermountain Healthcare Institutional Review Board. Mucinous tumors were included in the study based on the following criteria: 1. Diagnosis of invasive mucinous adenocarcinoma; OR 2. Diagnosis of mucinous bronchioloalveolar carcinoma (older diagnostic term for IMA); OR 3. Mucinous lung adenocarcinoma with strong morphologic features of IMA, such as abundant intracellular mucin. Colloid adenocarcinomas were excluded, as these are considered to be a discrete subtype of lung cancer.
Fluorescence-activated cell sorting (FACS) for in vivo tumors and in vitro organoid cultures
Tumor bearing lungs were isolated in ice-cold PBS. Lungs were enzymatically digested as described above in the procedure for establishing primary organoids. Once single-cell suspensions were obtained, cells were reconstituted in phenol red-free DMEM/F-12 with HEPES containing 2% FBS, 2% BSA and DAPI. Cells were sorted using BD FACS Aria for tdTomato-positive and DAPI-negative cells into DMEM/F-12% and 30% serum. After sorting, cells were pelleted and flash frozen or resuspended TRIzol then frozen at −80°C.
To study the dynamics of WNT signaling in vitro, BPN organoids were stably transduced with the β-catenin/TCF-dependent GFP reporter lentiviral construct 7TGP (Addgene plasmid #24305). Single cells were isolated from Matrigel using Cell Recovery Solution and TrypLE incubations. Cells were analyzed on BD LSRFortessa or BD FACSAria and sorted on BD FACSAria. Events first passed through a routine light-scatter and doublet discrimination gates, followed by exclusion of DAPI-positive dead cells. Viable tumor organoids were further identified as tdTomato-high cells. 7TGP organoids grown in 5% L-WRN media for 2–3 days were used as control for setting the GFP-negative gate. This gating strategy was applied to 7TGP organoids that were cultured in standard 50% L-WRN media in order to identify and sort for the WNT-reporter-low fraction. This fraction was collected in 50% L-WRN and expanded for ~1 week. After expansion, sorted 7TGP organoids were seeded and treated with 50% L-WRN (± Cobimetinib) or 5% L-WRN (± Cobimetinib) for 24 hr, 48 hr, and 72 hr. Then, the fraction of WNT-reporter-high population under the above four conditions was analyzed. Flow cytometric data analyses was performed using FlowJo software.
RNA sequencing from whole tumors
Library preparation was performed using the TruSeq Stranded mRNA Library Preparation Kit with poly(A) selection (Illumina). Purified libraries were qualified on an Agilent Technologies 2200 TapeStation using a D1000 ScreenTape assay. The molarity of adapter-modified molecules was defined by quantitative PCR using the Kapa Biosystems Kapa Library Quant Kit. Individual libraries were normalized to 5 nM and equal volumes were pooled in preparation for Illumina sequence analysis. Libraries were sequenced on Illumina HiSeq 2500 (50 cycle single-read sequencing v4).
Whole-tumor RNA-seq analysis
Adapters in raw FASTQ files containing 50 bp single-end reads were trimmed with cutadapt. QC metrics were generated for each sample with FastQC and Picard’s CollectRnaSeqMetrics after aligning reads to genome with STAR (Dobin et al., 2013). QC metrics were summarized by MultiQC. From this QC analysis, one sample (id 14489 × 13) did not group with other samples in the MultiQC report. It was later confirmed that very limited starting material was a contributing factor to its difference from the others. This sample was removed from downstream analysis because the noise in the sample was likely a larger negative trade-off than the gains from increasing biological replicates for that condition.
Next, Salmon v0.9.1 quantified RNA transcript expression with two steps (Patro et al., 2017); first, the mouse transcriptome downloaded from Gencode (release 26) was indexed using k-mers of length 19 with type ‘quasi’; second, Salmon quasi-mapped reads and corrected for sequence-specific bias (see option –seqBias). Salmon-based transcript expression estimates were converted to gene expression estimates with R package tximport (Soneson et al., 2015).
Differential gene expression modeling with DESeq2 (Love et al., 2014) and EdgeR (Robinson et al., 2010) evaluated four distinct contrasts testing (i) genotype effects within control groups [BP C vs. BPN C]; (ii) genotype effects within treatment groups [BP Tx vs. BPN Tx]; (iii) treatment effects within BP genotype [BP C vs. BP Tx]; (iv) treatment effects within BPN genotype [BPN C vs. BPN Tx]. Prior to fitting a single-factor DESeq2 regression (four factor levels: BP C, BPN C, BP Tx, BPN Tx), genes were filtered if there were less than five total counts across all samples. PCA was done on the top 500 most variable genes subject to a regularized log transform to confirm within-group variation was similar. A Wald test was applied to the fitted model for each of the four contrasts specified above, where the null hypothesis was gene expression differences less than or equal to log2(2) fold change in absolute value and the alternative hypothesis was that differences exceed log2(2). Significance of gene-wise differences was controlled by a false discovery rate of 10% (Benjamini and Hochberg, 1995).
Differential gene expression modeling with EdgeR for volcano plots used package edgeR_3.24.3 with ggplot2_3.1.1. Significance of gene-wise differences was calculated by tag-wise exact test. WNT signaling genes were delineated in the whole-tumor sequencing matrix using AmiGO pathway annotations (http://amigo.geneontology.org) for ‘Wnt Signaling’ (Carbon et al., 2009). WNT genes were restricted to those annotated for
Gene set enrichment analysis (GSEA) was carried out on RNA-seq data from whole tumors (and single cell data (below)) by comparing gene expression profiles with archived gene sets from Hallmarks [ftp.broadinstitute.org://pub/gsea/gene_sets/h.all.v7.0.symbols.gmt] and Oncogenic signatures [ftp.broadinstitute.org://pub/gsea/gene_sets/c6.all.v7.0.symbols.gmt] as well as with the cell type specific genes sets described in Haber et al., 2017; Leushacke et al., 2017; Montoro et al., 2018; Zhang et al., 2019.
Single-cell RNA sequencing
All protocols to generate scRNA-seq data on 10x Genomics Chromium platform including library prep, instrument and sequencing setting can be found at: https://support.10xgenomics.com/single-cell-gene-expression.
The Chromium Single Cell Gene Expression Solution with 3’ chemistry, version 3 (PN-1000075) was used to barcode individual cells with 16 bp 10x Barcode and to tag cell specific transcript molecules with 10 bp Unique Molecular Identifier (UMI) according to the manufacturer’s instructions. The following protocol was performed at High-Throughput Genomics Shared Resources at Huntsman Cancer Institute, University of Utah. First, FACS sorted single-cell suspensions of tdTomato-positive lung tumors were resuspended in phosphate buffered saline with 0.04% bovine serum albumin. The cell suspension was filtered through 40 micron cell strainer. Viability and cell count were assessed on Countess I (Thermo Scientific). Equilibrium to targeted cell recovery of 6,000 cells, along with 10x Gel Beads and reverse transcription reagents were loaded to Chromium Single Cell A (PN-120236) to form Gel Beads-in emulsions (GEMs), the nano-droplets. Within individual GEMs, cDNA generated from captured and barcoded mRNA was synthesized by reverse transcription at the setting of 53°C for 45 min followed by 85°C for 5 min. Subsequent A tailing, end repair, adaptor ligation and sample indexing were performed in bulk according to the manufacturer’s instructions. he resulting barcoding libraries were qualified on Agilent D1000 ScreenTape on Agilent Technology 2200 TapeStation system and quantified by quantification PCR using KAPA Biosystems Library Quantification Kit for Illumine Platforms (KK4842). Multiple libraries were then normalized and sequenced on NovaSeq 6000 with 2 × 150 PE mode.
Single-cell RNA-seq analysis
Demultiplexing and data alignment
Single-cell RNA-seq data were demultiplexed with 10x cellranger mkfastq version 3.1.0 to create fastq files with the I1 sample index, R1 cell barcode+UMI and R2 sequence. Reads were aligned to the mouse genome (mm10 with custom tdTomato reference) and feature counts were generated using cellranger count 3.1.0 with expected-cells set to 6000 per library. QC reporting, clustering, dimension reduction, and differential gene expression analysis were performed in 10x Genomics’ Cell Loupe Browser (v3.1.1). For the BPN control sample, we captured 5065 cells and obtained 46,894 mean reads per cell; 3544 median genes per cell. For the BPN treated sample, we captured 5563 cells and obtained 37,488 mean reads per cell; 2771 median genes per cell. For further details of the primary Cell Ranger data processing, see https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/algorithms/.
Visualization and expression analysis
Single-cell expression data was analyzed using Seurat (3.1.0) (Butler et al., 2018; Stuart et al., 2019) and visualized in Loupe browser (3.1.1). Alternatively, single-cell data was processed mostly in R (3.5.1). Cells with unique feature counts over 7500 or less than 200 and more than 20% mitochondrial counts have been filtered out. For visualization, the umap-learn package, from Seurat (McInnes et al., 2018), was used to reduce the dimensionality of the data.
Clustering
The Seurat package clustering method was used to designate cluster membership of single cell transcriptomics profiles. Seurat clustering is heavily inspired by Levine et al., 2015 and Xu and Su, 2015, which has two steps in the first step, the findNeighbors function is employed to embed single-cell profiles in a K-nearest neighbor (KNN) graph based euclidean distance in a (20 PCs) PCA space, followed by refining the edge weights based on jaccard similarity (https://satijalab.org/seurat/v3.0/pbmc3k_tutorial.html). In the second step, the FindClusters function employs Louvain algorithm modularity optimization techniques to iteratively group cells together (https://satijalab.org/seurat/v3.0/pbmc3k_tutorial.html).
Stromal cell designations
Clusters of cells sharing highly similar positions in dimensionality reduced UMAP-space were examined for differential expression using Loupe browser’s built-in differential expression calculator, and genes showing significant cluster specific expression were screened manually against the literature for markers of known cell types. Those clusters representing minority populations that were highly distinct in Umap-space from the bulk of cells in the analysis and which bear distinctive expression of stromal marker genes were considered as stromal contaminants and removed from subsequent, tumor-specific Umap data visualization and cluster specific gene expression level analysis.
Differential gene expression using seurat
For the DEG test the ‘FindMarkers’ function from the Seurat package with default setting was used, including Wilcoxon Rank Sum test, and Bonferroni multi-testing correction based on the total number of genes in the dataset.
Cell cycle phase modeling
We calculated cell cycle scores of control and treated tumor cells in scRNA-seq data using the Seurat package and the cell cycle genes defined by Mizuno et al., 2009. The Seurat package uses the scoring strategy described in Tirosh et al., 2016 and assigns scores for cells based on the expression of G2/M and S phase marker genes, with low scoring cells defaulting to G1. In an alternative approach, we modeled cell cycle phase by determining mean expression per cell of each phased cell cycle signature described in Mizuno et al., 2009. A matrix of these values was then graphed as a diffusion map using the Destiny package (Angerer et al., 2016) in R. A coherent cell cycle loop through three dimensional space using the first three diffusion components was projected onto two dimensions using an heuristic arithmetic, diffusion component combination. To this graph a principal curve was fit using Princurve (Hastie and Stuetzle, 1989) and tangents for each point (lambdas) were used as cellular positions along a trajectory designated as Q-depth based on the observed position of G0-signature enriched cells (an indicator related to quiescence) at the extreme distal to cells enriched in G1 signatures or the S and G2/M signatures.
MEK activation signature and WNT signature expression
Average expression of an 18 gene signature from Dry et al., 2010 was determined per cell across the single cell data set to model MEK activation state. Enrichment for the WNT13 signature described above was similarly examined across the single-cell dataset.
Data imputation
We used a novel data imputation approach to model the expression of low-detection genes in single-cell data (https://github.com/TheSpikeLab/FIESTA; Mehrabad et al., 2021a; Mehrabad et al., 2021b). Briefly, a weighted non-negative matrix factorization was carried out on the single-cell data matrix using machine-learning-based factor optimization. Subsequently, optimized factors were used to reconstruct an idealized completed matrix with values recommended for many gene nodes that were originally below the detection threshold.
Recombination fidelity
Cluster-specific BAM files were parsed from the concatenated scRNA-Seq BAM files using subset-bam tool (https://github.com/10XGenomics/subset-bam, 10xGenomics, 2020). The presence of select splice junctions was visualized using IGV (version 2.7.2) (Robinson et al., 2011) against mouse genome version (mm10).
Statistics p-values were calculated using the unpaired two-tailed t-test. RNA-seq statistics are described above.
© 2021, Zewdu et al. This work is published under https://creativecommons.org/licenses/by/4.0/ (the “License”). Notwithstanding the ProQuest Terms and Conditions, you may use this content in accordance with the terms of the License.