In addition to gene expression control through epigenetic mechanisms, transcription factors are master regulators of important cellular phenotypes, including stemness, proliferation, and migration23,24. Forkhead box protein P2 (FOXP2) is a member of the FOX transcription factor family that share conserved forkhead or “winged-helix” DNA-binding domains25. In addition to its role in the development of speech, FOXP2 is also known to be expressed in developing embryonic lung and can be used to mark induced pluripotent stem cell-derived alveolar lung progenitors26,27. Knockout of Foxp2 led to defective postnatal lung alveolarization in murine models28. In the adult, FOXP2 is expressed in the distal lung epithelium and inhibits lung cell differentiation genes29,30. Despite clues suggesting the importance of FOXP2 in maintaining stemness in lung epithelium, the roles it plays in lung adenocarcinoma have yet to be explored.
Results
Ezh2 haplo- and full Insufficiency drive distinct phenotypes in KRAS+/Trp53-null lung cancer
We generated the LSL:KrasG12D/+; Trp53flox/flox (LSL: lox-stop-lox) genetically engineered mouse model with either zero, one, or two floxed alleles of Ezh2 and induced autochthonous lung tumors by intranasal adeno-Cre (AdCre) virus administration (Fig. 1a). We analyzed the resulting KRAS+/Trp53-null lung tumors at a variety of timepoints. Mice with one floxed allele of Ezh2 (Ezh2 heterozygous; ∆/+) had significantly lower tumor burden when compared with Ezh2 wild-type (WT) mice (Fig. 1b, c and Supplementary Fig. 1a). In contrast, mice with two floxed alleles of Ezh2 (Ezh2 null, ∆/∆) had significantly higher tumor burden, significantly shorter survival, and significantly higher nuclear grades when compared to mice with Ezh2 WT or heterozygous tumors (Fig. 1d, e). Strikingly, Ezh2 null tumor cells had a significantly higher propensity to metastasize to areas including lymph nodes, pleural space, chest wall, diaphragm, and liver than Ezh2 WT or heterozygous tumor cells (Fig. 1f). Immunohistochemical (IHC) staining showed that Ezh2 null tumor cells had significantly lower EZH2 expression than Ezh2 WT, and that Ezh2 heterozygous cells had intermediate levels of EZH2 expression (Fig. 1g and Supplementary Fig. 1b, c). Confirming that PRC2 activity in Ezh2 null tumors was deficient, H3K27me3 staining was significantly lower in Ezh2 null tumors than in the other two genotypes (Supplementary Fig. 1d). Unlike Eed null tumors, Ezh2 null tumors had rare areas of alcian blue (mucin) stain, but did sometimes show more collagen staining, consistent with an epithelial-to-mesenchymal (EMT) phenotype (Supplementary Fig. 1e).
a Schematic for LoxP-mediated deletion of Ezh2 to generate EZH2 WT, heterozygous, and homozygous null KRAS/Trp53-null GEMMs. b Representative H&E-stained cross sections from lungs of closely related mice that were induced and sacrificed at the same timepoint, scale bar = 5 mm. c Tumor burden relative to the whole lung at the different timepoints post Adeno-Cre graphed as mean values +/− SEM. * indicates P = 0.0461 Ezh2 WT vs. null Bin 2, P = 0.028 Bin 3, P = 0.0437 Bin 4, ** indicates P = 0.005, *** indicates P = 0.0001 Ezh2 heterozygous vs. null by one-way ANOVA with multiple comparisons. Bin 1 n = 4, 9, 6, Bin 2 n = 10, 10, 9, Bin 3 n = 11, 14, 6, Bin 4 n = 9, 11, 9 for Ezh2 WT, heterozygous, null individual mice, respectively. d Percentage survival of mice of indicated genotypes. * indicates P = 0.0376 Ezh2 WT vs. Ezh2 heterozygous, *** indicates P = 0.0003 Ezh2 heterozygous vs. Ezh2 null, and * indicates P = 0.0397 Ezh2 WT vs. Ezh2 null by log-rank test. n = 8, 18, 18 for Ezh2 WT, heterozygous, null individual mice, respectively. e Nuclear grade was scored on sectioned tissue by a blinded pathologist and graphed as mean values +/− SEM. * indicates P < 0.042 by one-way ANOVA with multiple comparisons, n = 6, 9, 8 for Ezh2 WT, heterozygous, null individual mice, respectively. f Mice of indicated genotypes were scored as having metastasis (Met) or not by gross examination and histology at days 91–133 post tumor induction and percentage of each are graphed. * indicates P = 0.012 and ** indicates P = 0.0024 by Fisher’s exact test. n = 30, 36, 25 for Ezh2 WT, heterozygous, null individual mice, respectively. g Representative images of tumor sections stained with H&E, EZH2 or H3K27me3, scale bar = 100 µm. n = 26 or greater individual mouse tumor sections were stained and data are summarized in Supplementary Fig. 1a, b, d. Genotypes are signified as: WT, +/+; heterozygous, ∆/+; null, ∆/∆. Source Data are provided for this figure.
To further understand the mechanisms at play in the immune microenvironment, we performed RNA-sequencing (RNA-seq) on bulk tumor tissue containing immune cells. We first utilized TIMER2.031,32 and specifically used the CIBERSORT algorithm33 to estimate the immune cell populations. Ezh2 heterozygous tumors contained significantly more tumor-infiltrating CD8+ T cells than the other two Ezh2 genotypes, and significantly more monocytes than Ezh2 WT tumors (Supplementary Fig. 1f). Gene set enrichment analysis (GSEA)34 demonstrated that Ezh2 heterozygous tumors were significantly less enriched for cold tumor gene signatures, but had enrichment for lymphocyte signatures, indicating an anti-tumor immune microenvironment (Supplementary Fig. 1g and Supplementary Table 1). We next used HALO® artificial intelligence software to analyze the immune cell infiltration (Supplementary Fig. 1h). Analysis showed that Ezh2 heterozygous tumors had significantly more tumor-associated lymphocytes/monocytes when compared with Ezh2 null tumors (Supplementary Fig. 1i). We confirmed by IHC that Ezh2 heterozygous tumors had significantly higher infiltration of CD3+ T cells when compared with the other two genotypes (Supplementary Fig. 1j, k). Collectively, these data suggest that full and partial insufficiency for Ezh2 drive distinct phenotypic outcomes in KRAS+/Trp53-null lung adenocarcinomas, with full insufficiency driving more aggressive and metastatic tumors, while partial deficiency generates an anti-tumor immune microenvironment.
Model systems differ in cell cycle and gene expression
Next, we aimed to understand the epigenetic and transcriptional differences that underlie the in vivo tumor phenotypes, but prior to doing this, we characterized four model systems that can be used for this purpose: (1) tumor-derived 2D cell lines that are easy to grow and manipulate; (2) 3D tumoroids in Matrigel transwells that can be expanded, but less easily than 2D cells; (3) freshly sorted tumor cells that represent tumor cells in vivo, but are difficult to obtain in large numbers (Supplementary Fig. 2a); and (4) total (bulk) tumors that have the caveat of containing numerous stromal cell types. Consistent with previous reports10, 2D Ezh2 null cultures exhibited more elongated mesenchymal-like cells than 2D Ezh2 WT cultures, whereas in cross sections of 3D tumoroids and in vivo tumors, the mesenchymal phenotype was less obvious (Fig. 2a). To confirm a functional EMT in 2D Ezh2 null cells, we performed transwell migration assays, and observed significantly more migration of Ezh2 null cells than Ezh2 WT cells (Supplementary Fig. 2b). As expected, we observed that in both 2D and 3D states, Ezh2 null cells expressed lower levels of Ezh2 than the other two genotypes (Supplementary Fig. 2c). To confirm excision of the Ezh2 gene, we used endpoint PCR of the Ezh2 genomic DNA to verify the WT, floxed intact and floxed excised Ezh2 alleles in the tumoroids with the corresponding genotypes (Supplementary Fig. 2d).
a Representative images of in vivo tumor histology, in vitro 3D tumor spheroid histology and 2D cell cultures of the indicated Ezh2 genotypes, scale bar = 100 µm. n = 3 or greater individual mouse tumor-derived cultures were examined for each genotype and culture system. b Heatmap of sample-to-sample correlations using highest loading genes from PC3 analysis, where the Pearson correlations are weighted by sample-sample information and sample-sample consistency (see “Methods”). The 1 − correlation value was then used as a pseudo distance for hierarchical clustering and ordered using default argument of dendsort. n = 17, 13, 20, 21 for total tumor, sorted tumor, 2D cell cultures, 3D tumoroid individual mouse tumor-derived samples, respectively. c Rank-ordered gene lists were queried against the MSigDB databases and enrichment scores of selected gene signatures for 2D, 3D, or in vivo sorted tumors were plotted. Dot size reflect FDR estimates. Data are shown in Supplementary Table 3. d Percentages of cells in each cell cycle phase were measured by 7AAD cell cycle flow cytometry in 2D cells and 3D tumoroids of the indicated Ezh2 genotypes plotted as mean values +/− SEM. * indicates P < 0.035, ** indicates P < 0.009, *** indicates P < 0.0009 by one-way ANOVA with multiple comparisons and Holm–Šídák’s post hoc test for comparing genotypes within culture system, and by two-tailed t test for comparing the same genotype in the two culture systems. 2D: n = 4, 3, 4; 3D: n = 3, 4, 3 for Ezh2 WT, heterozygous, null individual mouse tumor-derived cultures, respectively. e 2D cells and 3D tumoroids were examined for expression of the indicated proteins by immunoblot. f Normalized protein abundance of the indicated proteins in 2D cells and 3D cultures by immunoblot. **** indicates P = 0.0000086 for pH3 in 2D vs. 3D and *P = 0.024 for PCNA in 2D vs. 3D for all Ezh2 genotypes combined by two-tailed t test. n = 2 individual mouse tumor-derived cultures per genotype/culture system. Genotypes are signified as: WT, +/+; heterozygous, ∆/+; null, ∆/∆. Source Data are provided for this figure.
To understand the transcriptional landscapes in each model systems, we compared the four models (independent of Ezh2 genotype) by RNA-sequencing. Unsupervised hierarchical clustering of samples using distances derived from weighted correlations showed that 2D cell lines form their own cluster that is distinct from all other samples, while in vivo samples clustered closely together and tumoroids were clustered between in vivo and 2D samples, regardless of Ezh2 phenotype (Fig. 2b). Principal component (PC) analysis performed on all samples determined that PC3 showed a stratification of the 2D samples away from the other three groups (Supplementary Fig. 2e). This PC is enriched for genes associated with epithelial cell proliferation, lung development, wound healing, cell–cell interactions, and immune response (Supplementary Table 2). Similarly, GSEA showed that when compared with 2D cultures, tumoroids and sorted tumor cells had enrichment for transcriptional signatures related to immunity and lung lineage determination, while 2D cells were enriched for cell cycle gene signatures (Fig. 2c and Supplementary Table 3).
To functionally verify differences in cell cycle, we next compared the cell cycle profiles of 2D cell lines and 3D tumoroids and found that they were dramatically different. 2D cultures showed much higher percentages of cells in S-phase and much lower percentages of cells in G0/G1 than 3D tumoroids for every Ezh2 genotype (Fig. 2d). Intriguingly, Ezh2 heterozygous tumoroids had increased G0/G1, while Ezh2 null tumoroids had more G2/M when compared to both Ezh2 WT and Ezh2 heterozygous, and less G0/G1 compared with mono-allelic deletion. These results may partially explain the opposing phenotypes driven by Ezh2 heterozygous and homozygous deletion in vivo. Consistent with the cell cycle results, 2D cultures had significantly higher levels of proliferating cell nuclear antigen (PCNA) and phospho-H3 than tumoroids (Fig. 2e, f). Intriguingly, Ezh2 null tumoroids expressed significantly higher levels of PCNA than Ezh2 WT and heterozygous tumoroids, while all 2D cells expressed comparable levels of PCNA. These data demonstrate critical differences in transcriptional states and growth patterns among culture systems, which should be considered when designing and testing novel therapeutic strategies.
Ezh2 deficiency leads to decreased PRC2-chromatin occupancy and diverse gene expression in different models
To further study the molecular consequences of Ezh2 deletion in vitro and in vivo, we used RNA-sequencing to investigate gene expression differences between the Ezh2 genotypes in each model system (Supplementary Fig. 2f). GSEA comparing Ezh2 null to Ezh2 WT cells, or Ezh2 null to Ezh2 heterozygous cells demonstrates the level of complexity of PRC2-mediated gene regulation in the different systems (Fig. 3a and Supplementary Table 4). When comparing Ezh2 null to Ezh2 WT cells, Polycomb targets genes, HOX genes, EMT, and immunity-related signatures were all enriched in Ezh2 null 2D cultures. In contrast, many of these same signatures were enriched in Ezh2 WT (depleted in Ezh2 null) in both sorted cells and tumoroids. These results suggest that compared to 2D cells, the tumoroids more faithfully recapitulate gene expression patterns of in vivo samples. Interestingly, Polycomb targets genes were more uniformly enriched (expressed) in both the Ezh2 null and Ezh2 WT cells when compared to Ezh2 heterozygous cells, including in the sorted cells and tumoroids (Supplementary Fig. 3a). This result suggests that Ezh2 heterozygous tumors may have stabilization of canonical PRC2 activity rather than a partial loss, and this is consistent with the observation that genes involved in histone methylation, including Polycomb complex genes, are more highly enriched in Ezh2 heterozygous cells when compared with the other two genotypes. Lastly, relative to Ezh2 null tumoroids, Ezh2 heterozygous tumoroids had enrichment of immunity gene signatures, predominately involving type-I interferon response pathways, suggesting that Ezh2 heterozygous cells may have hyperactivation of the damage-associated molecular pattern (DAMPs) and STING signaling35.
a Rank-ordered gene lists were queried against the MSigDB databases and enrichment scores of selected gene signatures enriched in the Ezh2 null or Ezh2 heterozygous compared to Ezh2 WT for the indicated models were plotted. Dots size estimates FDR. Data are summarized in Supplementary Table 4. b Venn diagram and heatmaps depict significant differentially expressed genes in Ezh2 null vs. WT cells for the indicated culture systems, with log-fold change LFC > 0.5 or LFC < −0.5, FDR < 0.05. Genes higher in Ezh2 null are in red, genes lower in Ezh2 null are in blue. Heatmaps depict LFC of selected genes identified by the indicated contrasts. c Heatmap representation of EZH2 and H3K27me3-bound chromatin peaks centered across a ± 5-kb window that showed decreased occupancy of Ezh2 heterozygous, Ezh2 null, and EPZ6438-treated Ezh2 WT mouse tumor spheroids compared to Ezh2 WT. d Venn diagram showing EZH2 and H3K27me3 peaks significantly called by SICER at FDR < E-10 for genes in the indicated Ezh2 genotypes or EZH2 inhibitor-treated organoids. e GREAT analysis heatmaps depicting the H3K27me3 peaks of mouse 3D organoid ChIP-seq samples enriched in the indicated gene sets. Data are summarized in Supplementary Table 6. f Venn diagram showing overlap of shared DEGs from (b) with H3K27me3 ChIP peaks found in Ezh2 WT tumoroids that were decreased by more than twofold or absent in Ezh2 null tumoroids, with Foxp2 indicated as a shared gene. Genotypes are signified as: WT, +/+; heterozygous, ∆/+; null, ∆/∆. Source Data are provided for this figure.
We next queried significantly differentially expressed genes (DEGs) between the Ezh2 null and Ezh2 WT cells in the different systems (Fig. 3b). Consistent with the GSEA, Hox genes were upregulated in Ezh2 null 2D cells, but Hox genes were not expressed in tumoroids and sorted tumors. In tumoroids, genes relating to lung cell differentiation such as Spdef, Muc5ac, Lamp3, and Trp63 were differentially regulated in Ezh2 null cells, whereas many of these same genes were too lowly expressed to assess in 2D cell lines. We closely examined the seven genes that were significantly upregulated in both tumoroids and sorted tumor Ezh2 null cells. Of these genes, Forkhead Box Protein 2 (Foxp2) encodes for a transcription factor vital in embryonic lung development and was significantly higher in Ezh2 null tumoroids and in vivo sorted tumors than in Ezh2 WT counterparts.
In order to understand the chromatin occupancy in PRC2-deficient cells, we performed chromatin immunoprecipitation and sequencing (ChIP-seq) for EZH2 in Ezh2 WT and heterozygous tumoroids, and for H3K27me3 in Ezh2 WT, heterozygous, null, and EPZ6438-treated WT tumoroids. Analysis of chromatin peaks showed EZH2 occupancy was comparable between the Ezh2 heterozygous and Ezh2 WT tumoroids, while the H3K27me3 occupancy on the Ezh2 heterozygous tumoroids was 75% of that in the Ezh2 WT tumoroids (Fig. 3c and Supplementary Fig. 3b–d). In Ezh2 null tumoroids, H3K27me3 peaks were 13% of the Ezh2 WT sample and nearly completely depleted in EPZ6438-treated Ezh2 WT tumoroids. Interestingly, despite similar numbers of EZH2 peaks in Ezh2 WT and heterozygous cells, only 45% of the peaks were shared between the genotypes (Fig. 3d). Similarly, there were a large number of H3K27me3 peaks that were unique to either Ezh2 WT and Ezh2 heterozygous cells. However, the majority of H3K27me3 peaks retained in Ezh2 null cells were also found in the other two genotypes and corresponded to genes encoding master transcriptional factor gene families, including the Hox, Sox and Nkx families (Supplementary Table 5 and Supplementary Fig. 3e). Further analysis of H3K27me3-bound regions showed that the majority of peaks retained in the Ezh2 null and EZH2 inhibitor-treated tumoroids were in gene promoter regions, while peaks in other genomic locations were lost, suggesting that the small amount of H3K27me3 remaining in cells is preferentially retained at promoter regions to best control gene expression (Supplementary Fig. 3f). GREAT analysis revealed that although Ezh2 heterozygous cells had lower enrichment of EZH2 at Polycomb targets genes, they retained H3K27me3 at Polycomb targets (Fig. 3e and Supplementary Table 6). Relative to Ezh2 WT, Ezh2 heterozygous tumoroids showed increased enrichment of EZH2 at genes associated with DNA damage repair, stemness and cell cycle, although these same gene signatures had only modest increases in H3K27me3 enrichment. Lastly, we examined H3K27me3 peaks in the Ezh2 WT and null tumoroids and discovered that 32% of the DEGs we previously identified had peaks that were reduced by at least twofold in Ezh2 null cells, including all seven of the DEGs shared between sorted and 3D cells, and notably, Foxp2 (Fig. 3f and Supplementary Fig. 3g).
Loss of PRC2 activity drives increased FOXP2 expression in normal and malignant lung cells
Given its critical roles in embryonic lung stem cells26,27, we next sought to characterize FOXP2 expression in the murine samples. First, we confirmed the expression of FOXP2 in 2D and 3D cells at the mRNA and protein levels. Foxp2 was elevated in Ezh2 null cells relative to Ezh2 WT cells in both 2D and 3D cultures (Fig. 4a). Furthermore, Foxp2 mRNA levels were significantly higher in Ezh2 null 3D tumoroids than in Ezh2 null 2D cells. FOXP2 protein expression in Ezh2 null tumoroids was also much higher than Ezh2 WT and heterozygous tumoroids and all the 2D cells (Fig. 4b). To confirm the correct band for FOXP2 using this antibody, we used small interfering RNAs to inhibit the transcription of Foxp2, and observed a decrease only in the band between 75 and 100 kDa (Supplementary Fig. 4a, b). To test whether EZH2/PRC2 regulates FOXP2 expression through its enzymatic activity, we re-expressed EZH2 WT and enzymatic-incompetent F667I mutant EZH2 in the Ezh2 null mouse tumor cells. We observed a significant reduction in FOXP2 protein expression only in EZH2 WT re-expression cells, but not in the F667I-mutant EZH2 re-expression cells (Fig. 4c, d). To support these in vitro findings, we verified that the FOXP2 expression in the mouse autochthonous tumors was significantly higher in the Ezh2 null tumors than in the other two genotypes by immunostaining (Fig. 4e). ChIP-seq showed that H3K27me3 and EZH2 were both enriched at the 5’ TSS of Foxp2 in the Ezh2 WT cells, but these peaks were substantially reduced in the Ezh2 null and EPZ6438-treated tumoroids (Fig. 4f and Supplementary Fig. 4c). ChIP-qPCR also showed reductions in H3K27me3 enrichment at the Foxp2 5’ transcription start site (TSS), regulatory element, and 3’ TSS, and significantly reduced enrichment at the sequence corresponding to the Foxp2 translation start site (ATG) in Ezh2 null cells compared with Ezh2 WT (Supplementary Fig. 4d).
a Relative expression of Foxp2 mRNA in mouse primary 2D tumor cells and 3D tumor spheroids with the indicated Ezh2 genotypes graphed as mean values +/− SEM. * indicates P < 0.043 by one-way ANOVA with multiple comparisons and Holm–Šídák’s post hoc test between genotypes and two-tailed t test between conditions. n = 3 individual mouse tumor-derived cultures per genotype/culture system. b FOXP2 protein expression in 2D cells and 3D tumoroids of the indicated Ezh2 genotypes was examined by immunoblot. Total histone H3 is the same as Fig. 2e. ǂ indicates the non-specific band seen with FOXP2 antibody from Cell Signaling Technologies. c Immunoblot and d the normalized abundance of the indicated proteins in Ezh2 WT, Ezh2 null, Ezh2 null with WT Ezh2 re-expression (WT), and Ezh2 null with F667I-mutant Ezh2 re-expression (F6771). Quantification data are plotted as mean values +/− SEM. * indicates P < 0.0261 by one-way ANOVA with multiple comparisons and Holm–Šídák’s post hoc test. n = 3 individual blotting experiments from two cell cultures. e Percentages of positively stained tumor nuclei for FOXP2 in the tumors of the indicated genotypes plotted as box-and-whisker plots, error bars are min–max, box bounds are 25th and 75th percentiles and center line is median. * indicates P = 0.0233 by one-way ANOVA with multiple comparisons and Holm–Šídák’s post hoc test for both comparisons. n = 29, 38, 23 for Ezh2 WT, heterozygous, null individual mice, respectively. f Visualization of H3K27me3 ChIP-seq peaks enriched in the enhancer and promoter regions of Foxp2 in the indicated 3D tumor spheroid samples. g Expression of FOXP2 mRNA in human lines HBEC3KT, BEAS2B, BEAS-KP, H460, H2030, A549, and H2009 with or without 5 µM EPZ6438 treatment for 10 days plotted as mean values +/− SEM. * indicates P < 0.043, ** indicates P < 0.009, *** indicates P = 0.0003 by two-tailed t test. n = 3 individual cell cultures per cell line/condition. h Immunoblot and i the normalized abundance of the indicated proteins in human lines HBEC3KT, BEAS2B, BEAS-KP, H460, H2030, A549, and H2009 with or without 5 µM EPZ6438 treatment for 10 days. Top FOXP2 panel is SIGMA antibody, bottom FOXP2 panel is Cell Signaling Technologies antibody. Quantification data are plotted as mean values +/− SEM. * indicates P < 0.032, *** indicates P = 0.000778 with two-tailed t test. n = 4 with two individual blotting experiments derived from one cell culture per cell line/condition each probed with two distinct FOXP2 antibodies. Genotypes are signified as: WT, +/+; heterozygous, ∆/+; null, ∆/∆. Source Data are provided for this figure.
Next, to test if EZH2 inhibition leads to derepression of FOXP2 in human samples, we measured FOXP2 expression in human cell lines treated with EPZ6438. Inhibition of EZH2 significantly elevated FOXP2 mRNA expression in human cancer lines H460, H2030, A549, and H2009, and immortalized pulmonary epithelial cells BEAS2B and BEAS2B with additional KRASG12V and TP53R175H mutations (BEAS-KP) (Fig. 4g). EZH2 inhibition also upregulated the protein expression of FOXP2 in human cancer lines H460, H2030, A549, and BEAS-KP, but not in immortalized epithelial cell lines HBEC3KT and BEAS2B (Fig. 4h, i). ChIP-qPCR showed a loss of H3K27me3 enrichment at a putative regulatory site 5’ of the FOXP2 gene in H460 and BEAS-KP cells treated with EPZ6438 (Supplementary Fig. 4e, f). To understand if the 2D vs. 3D expression differences observed in the murine cells could also be observed in human cells, H460 and H2030 were grown in both conditions and treated with EPZ6438. While EPZ6438 was able to drive increased FOXP2 expression in both systems, 3D human cells exhibited much higher FOXP2 mRNA levels than their 2D counterparts (Supplementary Fig. 4g). Together, these data suggest that EZH2/PRC2 directly regulates FOXP2 expression mainly through its EZH2 enzymatic activity, and that loss or inhibition of EZH2 increases FOXP2 expression in malignant lung cells of both mouse and human origin.
FOXP2 dictates stemness and migratory capacity of lung epithelial cells
Having established FOXP2 as a direct target of PRC2 in lung epithelial cells, we sought to understand if FOXP2 could drive the phenotypic differences observed in Ezh2 null tumor cells. We overexpressed FOXP2 in two immortalized pulmonary epithelial lines: HBEC3KT and BEAS2B, and two KRAS-driven lines: H460 and H2030 (Fig. 5a and Supplementary Fig. 5a). When grown in Matrigel cultures, FOXP2 overexpressing BEAS2B cultures yielded more, and significantly larger, organoids than the control cultures (Fig. 5b). Given the increased metastatic potential of Ezh2 null tumor cells, we tested if FOXP2 overexpression changed migratory capability by transwell migration assays. We observed that FOXP2 overexpressing cultures had significantly increased migratory capacity when compared to control cultures for all the cell lines tested (Fig. 5c, d). Furthermore, FOXP2 overexpression conferred BEAS2B culture’s resistance to etoposide and carboplatin (Supplementary Fig. 5b, c), consistent with the theory that cells in more stem-like states are relatively resistant to chemotherapy36.
a FOXP2 protein expression in HBEC3KT, BEAS2B, H460, and H2030 cells transduced with FOXP2 lentivirus was examined by immunoblot. EV is the empty vector control line. b Average diameters and organoid counts of BEAS2B 3D organoids with or without FOXP2 overexpression plotted as mean values +/− SEM. * indicates P = 0.034 by two-tailed t test. n = 3 individual cell cultures per group. c Representative images, scale bar = 100 µm, and d average numbers of migrated cells normalized to growth rates in HBEC3KT, BEAS2B, H460, and H2030 cells with or without FOXP2 overexpression plotted as mean +/− SEM. * indicates P < 0.047, ** indicates P = 0.0041 by two-tailed t test. n = 3, 4, 5, 5 for HBEC3KT, BEAS2B, H460, H2030 individual cell cultures per group, respectively. e Rank-ordered gene lists were queried against the MSigDB databases and enrichment scores of selected gene signatures enriched in FOXP2 overexpressing cell lines relative to EV lines were plotted. Dot size estimates FDR. Data are summarized in Supplementary Table 7. f FOXP2 protein expression in 2D H460 cells transduced with the indicated small hairpins was examined by immunoblot and representative of n = 4 blotting experiments from one cell culture per group. shGFP is the control cell line. g Crystal violet growth assays were performed on 2D H460 cells with or without FOXP2 knockdown at the indicated days of culture and mean values +/− SEM were graphed. ** indicates P < 0.0024, *** indicates P < 0.001 shFOXP2 compared to shGFP by two-tailed t test. n = 6 individual cell cultures per group. h Average diameters of 3D H460 spheroids with or without FOXP2 knockdown plotted as mean values +/− SEM. *** indicates P = 0.000017 shGFP vs shFOXP2.1 and P = 0.00000038 shGFP vs. shFOXP2.2 tumoroids by two-tailed t test. n = 6 individual cell cultures per group. i Average percentages of migrated cells normalized to shGFP and relative cell growth rates graphed as mean values +/− SEM. *** indicates P = 0.000014 compared to shGFP cells by two-tailed t test. n = 4 individual cell cultures per group. j Representative images of migrated cells in H460 with or without FOXP2 knockdown. Source Data are provided for this figure.
To discover the transcriptional programs changed by FOXP2 overexpression, we performed RNA-sequencing on FOXP2 overexpressing and control cultures. GSEA revealed that FOXP2 overexpression drove transcriptional amplification of many genes associated with EMT, cell adhesion, TNFα, RAS signaling, and lung development (Fig. 5e and Supplementary Table 7). The finding of RAS signaling enrichment was not unique to the H460 and H2030 KRAS mutated lines, but also observed in the HBEC3KT and BEAS2B lines. Lastly, to link FOXP2 overexpression in human cell lines to gene expression changes seen in Ezh2 null spheroids, we examined NDRG4 and MYCN in FOXP2 overexpressing or EZH2 inhibited H460 and H2030 tumoroids (Supplementary Fig. 5d, e). These two genes were highly upregulated in Ezh2 null tumoroids, and are potential PRC2/FOXP2 downstream genes based on ENCODE data37. In these cultures, the combination of FOXP2 overexpression and EZH2 inhibition led to significant increases in NDRG4 and MYCN compared to control cultures, confirming that these genes may be regulated by both EZH2 and FOXP2. Next, we tested whether the reduction of FOXP2 could reduce aggressive phenotypes of lung cancers. Knocking down FOXP2 with two different short hairpins in the human cell lines H460, A549, and H2009 remarkably impeded cell growth in 2D cultures (Fig. 5f, g and Supplementary Fig. 5f–i). FOXP2 knockdown also led to significantly smaller 3D tumoroids compared with the control shGFP in H460 cultures (Fig. 5h). Moreover, the H460 cells with FOXP2 knockdown migrated significantly less than the control, indicating a key role of FOXP2 in cell migration (Fig. 5i, j). Altogether, these discoveries identify a previously un-appreciated role for FOXP2 in the stemness and migration of normal and malignant adult lung epithelial cells.
Polycomb deficiency drives BET and histone demethylase inhibitor sensitivity
We next aimed to identify epigenetic therapeutic vulnerabilities of PRC2-deficient KRAS+/Trp53-null lung adenocarcinomas. We hypothesized that because Ezh2 null tumors were the most aggressive, methods to stabilize H3K27me3 or prevent the reading of H3K27ac could re-normalize the aberrant epigenetic state and reverse the aggressive phenotypes. Therefore, we focused our efforts on testing the histone demethylase inhibitor GSK-J4, and the BET inhibitor JQ1 (Fig. 6a). In 2D cultures, there were no differences in half-maximal inhibitory (IC50) values among the Ezh2 WT, heterozygous, and null cultures treated with GSK-J4 or JQ1 (Fig. 6b). However, Ezh2 null tumoroids were significantly more susceptible to GSK-J4 or JQ1 when compared with Ezh2 WT tumoroids (Fig. 6c and Supplementary Fig. 6a). Ezh2 heterozygous tumoroids were similarly sensitive to JQ1 as Ezh2 null, but only modestly sensitive to GSK-J4. The lack of consensus between 2D and 3D cultures prompted us to investigate drug vulnerabilities in vivo. Tumoroids were grafted subcutaneously onto immunocompromised mice. Intriguingly, consistent with the 3D culture system, the growth of Ezh2 null grafts was inhibited by GSK-J4, and significantly inhibited by JQ1 treatments, when compared to Ezh2 WT grafts, while the therapeutic effect of JQ1 on Ezh2 heterozygous grafts was between the other two genotypes (Fig. 6d and Supplementary Fig. 6b).
a Schematic depicting the mechanisms to restore epigenetic states lost in PRC2-deficient tumors. b, c GSK-J4 (upper) and JQ1 (lower) percentage survival of (b) 2D cultures after 4 days culture with drugs plotted as mean values +/− SEM. GSK-J4: n = 7, 7, 7; JQ1; n = 5, 6, 5 for Ezh2 WT, heterozygous, null individual mouse tumor-derived cultures, respectively. c 3D tumoroids after 12–14 days culture with drugs plotted as mean values +/− SEM. GSK-J4: * indicates P = 0.0158 Ezh2 WT vs. Ezh2 null, ***P = 0.0006 Ezh2 WT or heterozygous vs. Ezh2 null; JQ1: * indicates P < 0.039 and **P = 0.0084 for Ezh2 WT vs. Ezh2 null (lower) or Ezh2 WT vs. heterozygous (upper) by one-way ANOVA with multiple comparisons and Holm-Šídák’s post hoc. For both drugs, n = 7, 5, 5 for Ezh2 WT, heterozygous, null individual mouse tumor-derived cultures, respectively. d Changes in volumes of grafts in mice treated with GSK-J4 and JQ1 relative to changes in placebo-treated mice plotted over time as mean values +/− SEM. *** indicates P < 0.0008, **P < 0.0091, *P < 0.0472 for Ezh2 WT vs. Ezh2 null (lower) or Ezh2 heterozygous vs. Ezh2 null by one-way ANOVA with repeated measures and multiple comparisons. GSK-J4: n = 5, 5, 4; JQ1: n = 5, 5, 6 for Ezh2 WT, heterozygous, null individual mouse tumor-derived grafts, respectively. e Quantification of relative tumor volume change after 1 week and 2 weeks of placebo or JQ1 treatment of Ezh2 null or Ezh2 heterozygous tumor-bearing mice, ** indicates P = 0.0021 for week 1, P = 0.0003 for weak 2 by two-sided t test on log2-transformed values. f Representative magnetic resonance images of mouse lungs at baseline (week 0) and week 1 of JQ1 treatment for an Ezh2 null mouse. g Heatmap of Bliss synergy scores for EPZ6438 combined with JQ1 in H2009 and BEAS-KP cultures. The most synergistic areas had Bliss scores of 10.15 (H2009) and 10.99 (BEAS-KP), indicating synergy. n = 3, 4 for H2009 and BEAS-KP individual cell cultures, respectively. Genotypes are signified as: WT, +/+; heterozygous, ∆/+; null, ∆/∆. Source Data are provided for this figure.
To test if JQ1 could target tumors growing in the native lung microenvironment, we treated autochthonous tumor-bearing mice with JQ1 and measured changes in tumor volume by magnetic resonance imaging (MRI). Seven days of JQ1 treatment reduced the volumes of lung tumor nodules in Ezh2 null and heterozygous GEMMs, while tumors continued to grow in the placebo group (Fig. 6e, f and Supplementary Fig. 6c). An additional seven days of JQ1 administration was sufficient to maintain tumors in their regressed state, though no additional tumor burden reduction was observed. To further explore the clinical significance of the drug vulnerability discoveries from murine tumors, we investigated several KRAS mutated human non-small cell lung cancer lines. EZH2 knockdown sensitized the human line H2030 to JQ1 (Supplementary Fig. 6d, e), and such sensitization was more observable in 3D spheroid culture than 2D culture in H460 (Supplementary Fig. 6f). Additionally, drug matrix experiments demonstrated synergy when combining JQ1 with EPZ6438 in H2009 and BEAS-KP cells (Fig. 6g).
The differential response to BET inhibitor inspired us to interrogate whether FOXP2 was involved in this therapeutic effect in EZH2-deficient cells. First, we performed RNA-sequencing on the 3D tumoroids and subcutaneous grafts treated with JQ1, which revealed that BET inhibition led to significantly decreased Foxp2 expression in Ezh2 null cells (Supplementary Fig. 6g). c-Myc, a gene often changed by JQ1 treatment38, was also significantly decreased in expression in the JQ1-treated graft, confirming the drug’s efficacy. RT-qPCR on mouse 3D tumoroids showed that Foxp2 was significantly inhibited by 100 nM JQ1 in Ezh2 heterozygous and null tumoroids, but not in Ezh2 WT tumoroids (Supplementary Fig. 6h), which may explain the difference in drug sensitivity among the genotypes. In the JQ1-treated H2030 cultures, FOXP2 levels were significantly decreased in the shEZH2 lines, but not in the shGFP line, again suggesting that regulation of FOXP2 correlated to the increased sensitivity of cells to JQ1 (Supplementary Fig. 6i). Similarly, in EZH2 WT rescued Ezh2 null tumoroids, which we demonstrated had lowered FOXP2 levels, JQ1 was significantly less effective than in Ezh2 null tumoroids re-constituted with the F677I EZH2 mutant that retain high levels of FOXP2 (Supplementary Fig. 6j). Lastly, we observed H460 overexpressing FOXP2 were more sensitive to JQ1 than the control cells, and again the increase in sensitivity was more apparent in 3D cultures than in 2D (Supplementary Fig. 6k, l). Taken together, these data indicate reduction in PRC2 activity confers sensitivity to BET inhibition and demethylase inhibition in lung adenocarcinoma, and FOXP2 could be a promising target of BET inhibition.
Decoupling of EZH2 from PRC2 activity may explain co-expression of EZH2 and FOXP2 in patient samples
To extend our findings to human lung cancer patient samples, we performed immunohistochemical staining for EZH2, H3K27me3, and FOXP2 on lung cancer tissue microarray. Consistent with previous research6,9, the H3K27me3 stain was significantly higher in ADC than SCC, while EZH2 had the opposite pattern (Fig. 7a). FOXP2 expression was variable in ADC samples, but found to be significantly higher in poorly differentiated tumors than in the other three groups. Because the mouse model we used was adenocarcinoma, and became more poorly differentiated when Ezh2 was deleted, we next focused only on ADCs and poorly differentiated lung tumors. In this cohort, we observed that both well and moderately differentiated tumors had significantly lower EZH2, and significantly higher H3K27me3 marks, than poorly differentiated tumors, which had lower H3K27me3, but paradoxically had higher EZH2 levels (Fig. 7b, c). Intriguingly, EZH2 expression showed a strong positive correlation to FOXP2 expression (Fig. 7d). In tumors that were FOXP2 positive (>10% of nuclei) that mirror the high expression found in some Ezh2 heterozygous and most Ezh2 null tumors, there was a significant inverse relationship between EZH2 and H3K27me3 (Supplementary Fig. 7a). To study the correlation of EZH2/FOXP2 expression with patient prognosis, we queried published RNA-sequencing data from The Cancer Genome Atlas using KM-Plotter39 on 366 lung adenocarcinoma samples. Consistent with previous reports, higher EZH2 expression correlated with significantly worse relapse-free survival in lung adenocarcinoma (Supplementary Fig. 7b). In our tissue microarray, EZH2 expression also correlated with a worse overall survival, but it is important to note that in our dataset, SCC patients had a significantly shorter overall survival than did lung ADC patients, which coincides with the fact that SCC tumors usually have higher EZH2 expression (Supplementary Fig. 7c, d). Importantly, higher FOXP2 expression correlated with significantly shorter relapse-free survival up to 250 months post diagnosis (Fig. 7e).
a Positively stained tumor nuclei in ADC, ADSCC, poorly differentiated tumor and SCC specimens for the markers EZH2, H3K27me3 and FOXP2 plotted as box-and-whisker plots, error bars are min–max, box bounds are 25th and 75th percentiles and center line is median. * indicates P = 0.0164, ** indicates P < 0.0084, *** indicates P < 0.0006 by one-way ANOVA with multiple comparisons and Holm–Šídák’s post hoc test. n = 92, 24, 17, 104 for ADC, ADSCC, poorly differentiated and SCC individual tumors, respectively. b The percentage of EZH2 and H3K27me3 positive nuclei in lung adenocarcinoma and poorly differentiated cancer specimens. * indicates P < 0.0425, ** indicates P < 0.0086 by one-way ANOVA with multiple comparisons and Holm–Šídák’s post hoc test. n = 13, 51, 97 for well, moderately, poorly differentiated individual tumors, respectively. c Representative images of lung adenocarcinoma specimens defined as well, moderately, and poorly differentiated stained for the indicated markers, scale bar = 100 µm. The number of tumors stained are summarized in (b). d Correlation between EZH2 and FOXP2 based on the percentage of positively stained tumor nuclei in primary lung cancer specimens. Pearson’s correlation coefficients and p values are shown. n = 161 individual tumors. e Kaplan–Meier relapse-free survival curves for FOXP2-high and FOXP2-low tumors as measured by RNA-sequencing 250 months post diagnosis; positive and negative groups were determined by best cut-off. The P value was calculated by log-rank test. n = 366 individual tumors. f, g Immunoblot in 2D (f) H460, A549, BEAS-KP, and g Ezh2 WT UK803 cells treated with or without 500 μM SAM for 6 days. h Normalized protein abundance in H460, A549, and BEAS-KP cells treated with or without 500 μM SAM for 6 days plotted as mean values +/− SEM. * indicates P = 0.012, ** indicates P = 0.009 for EZH2 and P = 0.0045 for FOXP2 by two-tailed t test. n = 3 different cell lines. Genotypes are signified as: WT, +/+; heterozygous, ∆/+; null, ∆/∆. SAM, S-adenosyl methionine. Source Data are provided for this figure.
Given that EZH2 is rarely mutated in NSCLC, one possible reason for observing decoupling in tumors is that in highly proliferative tumors EZH2 is upregulated, but there may be a limited availability of the methyl donor, S-adenosyl methionine (SAM)10. To examine the hypothesis that SAM levels influence PRC2 stability and activity, 500 μM SAM was added to the culture media of human lines H460, A549, and BEAS-KP for 6 days. As expected, these human lines treated with SAM had higher H3K27me3, and interestingly had significantly lower EZH2 and FOXP2 compared to the controls (Fig. 7f, h). Treating an Ezh2 WT murine 2D culture with SAM led to a similar decoupling of EZH2 and the H3K27me3 mark (Fig. 7g and Supplementary Fig. 7e). In summary, EZH2 expression decouples from PRC2 activity, and a pattern of high EZH2, low H3K27me3 and high FOXP2 is enriched in poorly differentiated advanced lung cancers.
Discussion
Here, we uncovered un-appreciated roles for EZH2 heterozygosity and full deficiency in regulating tumor status, whereby Ezh2 heterozygous tumors had lower tumor burden, lower nuclear grade and were less lethal than Ezh2 null tumors. Although the Ezh2 heterozygous tumoroids had a 25% reduction in H3K27me3 marks, they had the strongest suppression of canonical PRC2 target genes, suggesting that increased canonical PRC2 activity, rather than partial loss, drives Ezh2 heterozygous phenotypes. These cells also had changes in cell cycle, increased type-I interferon response pathways, and in vivo had an enrichment for T cell infiltration, and these mechanisms may explain the lower tumor burden in Ezh2 heterozygous mice. In contrast, Ezh2 null tumors were more metastatic and higher grade. Despite the fact that EZH2 inhibitors may drive more aggressive phenotypes in cancer, there still may be clinical utility for EZH2 inhibitors due to effects on the tumor microenvironment, and the ability to drive additional therapeutic vulnerabilities of tumor cells. Our data suggest that in KRAS-driven lung cancers, a PRC2-deficient state, either occurring naturally or through EZH2 inhibition, may be vulnerable to BET inhibition. We observed that Ezh2 null tumoroids, subcutaneous grafts and autochthonous lung tumors were sensitive to the BET inhibitor JQ1. Mechanistically, we determined that FOXP2, which was previously identified as a key regulator of embryonic lung progenitor cells28, was downregulated in several cases where BET inhibition was effective, and in the H460 cell line, overexpression of FOXP2 conferred increased sensitivity to JQ1. We confirmed that FOXP2 was a direct target of PRC2-mediated gene silencing in malignant lung epithelial cells, and that enforced expression of FOXP2 increased stemness and migration of lung epithelial cells.
It is notable that our drug response findings and the identification of FOXP2 were discovered in vitro only due to the use of 3D tumoroid cultures. We found that murine 2D cultures had extremely high levels of S-phase and down-regulation of genes involved in lung lineage determination when compared to 3D tumoroids, which showed retention of lung specification transcriptional programs. This finding is consistent with the belief that 3D cultures are able to retain more characteristics of the native tissue from which they were derived40, and demonstrates the vast differences between the same cells cultured in 2D versus in 3D. The outcomes of loss of PRC2 activity in 3D and 2D cells were also divergent. In the Ezh2 null genotype, 2D cell lines had marked derepression of Hox genes, but in 3D tumor spheroids, Ezh2 null cells were able to continue to repress Hox genes and instead showed derepression of Foxp2. Together, these findings underscore the power of using 3D systems for genetic, epigenetic and drug validation experimentation in cancer research.
Finally, while the co-expression of FOXP2 and EZH2 in patient samples may seem contradictory to the results from our mouse model, we believe these results support each other. Several papers, including our own, have shown an inverse correlation between EZH2 and its catalytic mark, H3K27me3 in vivo6,7,8,9,41,42,43. This abundant EZH2 may be part of PRC2-independent complexes44,45,46, or unable to complete histone methylation for a variety of reasons (one reason postulated is lack of SAM). Therefore, high EZH2 can indicate lower PRC2 function, and our mouse model and in vitro data show that one result of lowering PRC2 function in KRAS-driven adenocarcinomas is to de-repress FOXP2. In this way, the mouse model, with EZH2 deletion to drive PRC2 dysfunction, and the human samples, with high EZH2 indicating PRC2 dysfunction, do match. Furthermore, retrospective analysis of patient sample RNA-sequencing data showed that high FOXP2 expression predicted significantly shorter relapse-free survival. Our data suggest that when SAM is abundant, which may be the case in more well-differentiated tumors and normal tissues, PRC2 activity is high and EZH2 is downregulated. However, in aggressive tumors that often have low levels of H3K27me37,47, EZH2 may be upregulated, but unable to perform its catalytic function due to lack of available SAM. Future research focuses on understanding the physiological or pathophysiological metabolic mechanisms by which the SAM cellular pool are regulated, especially in vivo, and whether this process can be reversed for therapeutic intervention.
