Altered Transcriptional Control Networks with Trans-Differentiation of Isogenic Mutant-KRas NSCLC Models

Background: The capacity of cancer cells to undergo epithelial mesenchymal trans-differentiation has been implicated as a factor driving metastasis, through the acquisition of enhanced migratory/invasive cell programs and the engagement of anti-apoptotic mechanisms promoting drug and radiation resistance. Our aim was to define molecular signaling changes associated with mesenchymal trans-differentiation in two KRas mutant NSCLC models. We focused on central transcription and epigenetic regulators predicted to be important for mesenchymal cell survival. Experimental design: We have modeled trans-differentiation and cancer stemness in inducible isogenic mutant-KRas H358 and A549 non-small cell lung cell backgrounds. As expected, our models show mesenchymal-like tumor cells acquire novel mechanisms of cellular signaling not apparent in their epithelial counterparts. We employed large-scale quantitative phosphoproteomic, proteomic, protein–protein interaction, RNA-Seq, and network function prediction approaches to dissect the molecular events associated with the establishment and maintenance of the mesenchymal state. Results: Gene-set enrichment and pathway prediction indicated BMI1, KDM5B, RUNX2, MYC/MAX, NFκB, LEF1, and HIF1 target networks were significantly enriched in the trans-differentiation of H358 and A549 NSCLC models. Physical overlaps between multiple networks implicate NR4A1 as an overlapping control between TCF and NFκB pathways. Enrichment correlations also indicated marked decrease in cell cycling, which occurred early in the EMT process. RNA abundance time course studies also indicated early expression of epigenetic and chromatin regulators within 8–24 h, including CITED4, RUNX3, CMBX1, and SIRT4. Conclusion: Multiple transcription and epigenetic pathways where altered between epithelial and mesenchymal tumor cell states, notably the polycomb repressive complex-1, HP1γ, and BAF/Swi-Snf. Network analysis suggests redundancy in the activation and inhibition of pathway regulators, notably factors controlling epithelial cell state. Through large-scale transcriptional and epigenetic cell reprograming, mesenchymal trans-differentiation can promote diversification of signaling networks potentially important in resistance to cancer therapies.


INTRODUCTION
Cellular plasticity in epithelial cancers is associated with a progression to a metastatic state (1,2) and resistance to anti-cancer therapies (3,4). Over 90% of cancer patient deaths are attributable to complications arising from metastatic dissemination of cancer cells to distant organ sites. Tumor plasticity associated with epithelial mesenchymal transition (EMT) (5,6) contributes to metastasis, drug resistance and is correlated with poor prognosis (7,8). The epigenetic reprograming associated with EMT promotes the disassembly of epithelial cell-junctions, a loss of epithelial polarity (4,9), and the formation of molecular assemblies allowing cell migration and invasion (10). In multiple model systems, EMT-derived mesenchymal cells can show the properties of cancer stem cells including tumor initiation at low cell inoculation in vivo, sphere formation in vitro, and CD44 high , CD24 low , and ALDH1 active pluripotent stem cell markers (11,12). Despite advances in the treatment of non-small cell lung cancers (NSCLC), for example, the development of EGFR tyrosine kinase inhibitors (EGFR TKIs) for patients with activating EGFR mutations, survival rates for patients with mutant-KRas lung cancer are poor. Mutant-KRas is observed in~20% of NSCLC cases and is generally independent of EGFR mutations (13). In a recent clinical study, the median survival time of patients treated with standard of care was 7.7 months in NSCLC patients with mutant-KRas, in marked contrast to 38 months in patients with mutant EGFR (14). Several studies suggest tumor cells harboring mutant-KRas may be primed to undergo epithelial mesenchymal trans-differentiation, for example (15,16). These mesenchymal-like carcinoma cells have been shown to be resistant to many conventional lung cancer therapies, including taxanes, pemetrexed, gemcitabine, and EGFR TKIs (8,12,17). Chemotherapy has been shown to promote selection for EMT-derived cells in a number of solid tumor types (7,17,18) and conversely markers of EMT appears to contribute to chemotherapy resistance (19) and predict response to EGFR and PI3K inhibitors (20). Similar data show EMT-derived cells can serve as reservoirs for cancer recurrence in relevant genetically engineered models (21,22). Therapies are needed that specifically target drug resistant mesenchymal-like tumor cells. However, survival signaling networks in EMT-derived cells appear heterogeneous and critical dependencies common to mesenchymal tumor cells remain ill defined. Current inhibitors of mesenchymal stemlike tumor cells have been largely restricted to specific tumor and driver types, for example, JAK or PDGFR inhibitors. Since mesenchymal trans-differentiation involves an epigenetic reprograming, the pharmacological use of epigenetic modulators would appear attractive. However, initial therapeutic successes with single agent DNA methyltransferase inhibitors, HDAC inhibitors, and EZH2 inhibitors have been more pronounced in the hematologic malignancies, as opposed to epithelial-derived carcinomas, which can undergo EMT.
In order to understand and identify transcriptional and epigenetic networks in a mutant-KRas mesenchymal cell context, we have modeled metastable reversible EMT in two NSCLC cell backgrounds. "Metastable" refers to a reversible EMT, achieved for example by withdrawal of inducer (23), where in intermediate stages cells can express both epithelial and mesenchymal markers (24). Here, we compared CD44 low A549 cells with isogenic CD44 high A549/transforming growth factor beta (TGFβ) cells and CDH1/E-cadherin high H358 cells with isogenic CDH1 low H358/dox-TGFβ cells. As expected both models reversibly undergo EMT and are less responsive to paclitaxel, gemcitabine, and erlotinib in the mesenchymal state (data not shown). These models confirm the clinical observations and show that tumor cells that have undergone EMT show increased resistance to standard of care cancer treatments. Importantly, upon undergoing EMT tumor cells acquire novel mechanisms of cellular signaling and resistance to apoptosis not apparent in their epithelial counterparts.
We sought to define at signaling network changes distinguishing epithelial and mesenchymal tumor states, with an aim of identifying key nodes important in mesenchymal cell survival. Pharmacological blockade of mesenchymal survival pathways may overcome a limitation of current therapies, which preferentially impact tumor cells with an epithelial phenotype (17). We used a combination of RNA-Seq, phosphoproteomic, and bioinformatics approaches to identify transcriptional and epigenetic networks modulated as a consequence of epithelial mesenchymal transition (EMT). We employed TGFβ to induce EMT in the A549 and H358 NSCLC models, where cells were preselected for a starting epithelial state. TGFβ is a physiologically relevant inducer of EMT (25) associated with a macrophage/monocyte rich inflammatory microenvironment (26). Detailed RNA, protein, and phosphopeptide abundance datasets were obtained. Using pathway prediction and gene-set enrichment approaches, and protein-protein interaction data we observed altered regulation of transcriptional, epigenetic, and chromatin modulators and explored intersections between key pathways.

RNA-Seq
Duplicate cell samples (~3 × 10 6 cells) were lyzed in 600 µl Qiagen RLT buffer, 1% β-mercaptoethanol, followed by column isolation (qiagen.com). RNA isolation, library construction, and RNA-Seq essentially followed Illumina protocols (illumina.com). cDNA was prepared from total RNA with globin RNA reduction, followed by library generation (Illumina True Seq). cDNA libraries were immobilized (Illumina, Standard Cluster Generation Kit) and sequenced on an Illumina HySeq instrument. Read depths for H358 and A549 were 30 million and 25 million, respectively. RNA-Seq was performed on Illumina HySeq instrument 1 . FASTQ file reads passing Illumina purity filter (PF reads) were aligned and quantitation performed by RSEM (29), generating files where normalized counts for each detected gene and isoform (UCSC Known Gene). Aligned RNAs (27976 genes; 48241 isoforms) passing QC thresholds were used to calculate mesenchymal-epithelial transcript abundance ratios followed by log2 linear scaling. Additional analysis was performed using Cufflinks (30) on BAM files aligned to hg19. For RNA-Seq data were selected where log2 fold change (FC) was ≥1.0 and where rsem reads were >5 in any sample. In difference analysis where either numerators or denominators were equal to zero, values were generated (specified as >5 or <5) in gene lists if differences >10 rsem reads were observed in both H358 and A549 models with concordant log2 signs. Concordance required >2 FC in three out of four samples (Mes-1 and Mes-2 over Epi-ave) for both A549 and H358 cell models.

PEPTIDE IDENTIFICATION AND QUANTIFICATION BY LC-TANDEM MS
LC-MS/MS was performed essentially as previously described (33,34). Two (A549) or three (H358) LC-MS/MS experiments were performed for each fraction (1 fraction for pY, 10 fractions for SCX/TiO2, and 10 fractions for SCX/total peptide). Proteins were identified from survey and product ion spectra data, using the Paragon algorithm of ProteinPilot [v4.0; (35)] and GPM [(36) v2.2.1]. Two missed tryptic cleavages were allowed and posttranslational modifications considered included cysteine derivitization, STY phosphorylation, deamidation, carbamylation, oxidation, and SILAC labels. Database searches used the human UniProt FASTA database (10-2012; 70,391 sequences including common contaminants). When multiple protein isoforms were identified, Protein Pilot allowed only peptides specific to each detected isoform to be used, which factored in ion counts for weighting in the protein ratio calculation (37). Parsimony of protein results was assured by rigorous protein inference with the ProGroup algorithm. Protein identification complied with the guidelines of Bradshaw et al. (38) where two or more unique isoform-specific peptides were required for inclusion. False discovery rates of phosphotyrosine peptide capture experiments were <1%. False positive rates for the high complexity TiO2 and total peptide assignments ranged between 1.4 and 0.8%. For statistical analysis we required four or more peptides with individual peptide assignments at >95% confidence. Phosphopeptide peak areas were normally distributed by log10 + 1 conversion followed by paired t -test. Phosphopeptides sites identified four or more times are listed in Table S2 in Supplementary Material and coordinate changes in phosphopeptide abundance between H358 and A549 models were observed for 73 unique sites are listed in Table S3 in Supplementary Material. Fifty-eight protein abundance changes were coordinately associated with H358 and A549 trans-differentiation, with two or more unique peptides also at 95% peptide ID confidence are listed in Table S4 in Supplementary Material. Co-correlation of overlapping peptides from skipped cleavage and multiple charge states were used to further reduce false discovery rates and correctly bin previously defined benchmark peptides from E-cadherin, vimentin, and fibronectin to their respective epithelial or mesenchymal cell states.

FUNCTIONAL ANALYSIS AND PERFORMANCE PARAMETERS
A schema for the statistical and categorical analysis of protein, phosphopeptide, and RNA transcript data is provided in Figure 1. Changes in the abundance of proteins, phosphopeptides, and RNA transcripts were compared between fractions, experiments, and models where membership within specific signaling networks only was established using two or more independent lines of evidence. Cross-correlation between two independent isogenic KRas models was used to reduce noise and

CHARACTERIZATION AND VALIDATION OF H358 AND A549 MUTANT-KRas MODELS
Two KRas mutant adenocarcinoma NSCLC cell lines H358 and A549 were used as model systems to molecularly define transcriptional and epigenetic reprograming following mesenchymal trans-differentiation. H358 and A549 cells can spontaneously generate populations of CDH1 high /CD44 low and CDH1 low /CD44 high , with epithelial and mesenchymal-like phenotypes, respectively (41). Spontaneous inter-conversion has been previously reported (42,43). H358 contain relatively rare (estimated at~2-4%) CDH1 low /CD44 high cells while A549 are typically a more mixed population of each phenotype. As such, all H358 experiments were initiated from epithelial CDH1 high /VIM low clones with predominant epithelial cell-junctions, from which subsequent H358/dox-TGFβ clones were derived. Doxycycline (0.5 µg/ml) induction of transgene expression, a constitutively active form of TGFβ1 (37), was verified by immunoblot and was shown to correctly modulate EMT marker expression (CDH1 low , CD44 high , and VIM high ) as shown in Figure 1A. Fluorescence microscopy ( Figure 1B) showed loss of E-cadherin membrane localization and gain of CD44 expression in H358/dox-TGFβ cells relative to the minus dox control. Multiple H358/dox-TGFβ clones exhibited correct isogenic mesenchymal trans-differentiation and in contrast vector control cells remained epithelial in the presence of doxycycline (24). Similarly, the percentage of aldefluor positive cells, a marker of aldehyde dehydrogenase activity and putative stemness, was increased after mesenchymal trans-differentiation with TGFβ for 14 days. The percentage of aldefluor positive cells was 6.2% in the −dox control H358 cells and 16.8% in the +dox H358/TGFβ producing cells ( Figure 1C). A549 cells were anti-CD44 antibody counter selected using magnetic beads into a CD44 low starting epithelial population, which was subsequently induced with TGFβ (10 ng/ml) to yield a uniformly mesenchymal cell population after 14 days. This allowed a more direct comparison of the CD44 low starting population with the TGFβ-induced mesenchymal population, which became CD44 high . Immunoblot for fibronectin, E-cadherin, and vimentin confirmed an EMT transition after 7 and 14 days exposure to TGFβ in A549 cells (10 ng/ml; Figure 1D). EMT occurs in H358 and A549 epithelial CDH1 + /CD44 low subclones exposed to TGFβ over a prolonged period.

INTEGRATION OF RNA, PROTEIN, AND PHOSPHOPROTEIN EMT STATE-SPECIFIC MEASUREMENTS
In order the globally assess differences between isogenic epithelial and mesenchymal cell states, we measured RNA, protein, and phosphorylation changes as outlined in Figure 1E. RNA-Seq (44) was performed where non-zero ratios with RSEM reads ≥5 for any condition are listed in Table S1 in Supplementary Material (20,443 genes). Correlation between biological replicate RNA-Seq samples was r 2 = 0.90 and r 2 = 0.84 for H358 and A549 models, respectively (Figure 2A). Protein and phosphopeptide changes between epithelial and mesenchymal cell states were measured using SILAC labeling (31) in H358 and A549 TGFβ-treated or control. Analysis of RNA and protein abundance ( Figure 2B) confirmed the gain of vimentin (VIM), fibronectin (FN1), and loss of keratin-8 (KRT8) and S100A6, expected RNA and protein changes occurring with mesenchymal trans-differentiation. We focused on changes in RNAs encoding transcription regulators, where 82 transcripts were coordinately regulated in comparing mesenchymal and epithelial cell states for both isogenic H358 and A549 models ( Figure 2C).

FUNCTIONAL ANNOTATION OF TRANSCRIPTION AND EPIGENETIC NETWORKS
Functional association with RNA abundance changes was assessed by GSEA. As expected multiple signatures associated with KRas transformation, stemness, and mesenchymal trans-differentiation were significantly enriched (

FIGURE 2 | (A)
Regression analysis of duplicate mesenchymal RNA-Seq samples expressed at a ratio to the mean epithelial control for both H358 and A549 models (greater than twofold changes are indicated in black). (B) Concordant RNA and protein changes for both isogenic H358 and A549 models.
(C) Transcription associated RNA transcripts with differential abundance between epithelial and mesenchymal cell states in both H358 and A549 models.
Frontiers in Oncology | Molecular and Cellular Oncology prediction analysis predicted activation of TGFβ signaling (45), and inhibition of DKK1 and SMAD7 signaling. These findings further served as validation benchmarks for the GSEA and pathway prediction approaches. Similarly benchmark RNA expression values for transcriptional regulators, which can induce EMT [Snail, Slug, Twist, and Zeb (9, 46-48)] were markedly increased in a model dependent manner (Table S6 in Supplementary Material). Transcriptional signatures surrounding LEF1, NFκB, and BMI1 were highly enriched in datasets comparing H358 and A549 epithelial and mesenchymal-like cell states using GSEA and/or pathway prediction analysis ( Table 1). We considered other predicted activated or inhibited nodes common to both KRas transdifferentiation models. Pathway prediction analysis was used to model signaling activation or inhibition based on transcription factor response. Table 2 summarizes pathway state based on H358 and A549 RNA-Seq datasets. Expected activation or Snail (H358), Twist (A549), and STAT pathways [H358 and A549 via increased IL11 and IL6; (23,49)] was expected but highlight model heterogeneity, even within closely related NSCLC systems. Interestingly, RNA expression correlated the mesenchymal state with KDM5B, SMARCA4, and EGR1 pathway activation, and Myc/Max, SIN3A, and SPDEF pathway inhibition.

ALTERED TCF/LEF AND NFκB NETWORKS WITH THE MESENCHYMAL CELL STATE
Canonical Wnt signaling leads to the nuclear translocation of βcatenin and the activation of TCF and LEF family transcription factors, which in turn promote pro-survival gene expression programs. The activation of the Wnt/β-catenin/TCF-LEF axis has been implicated in epithelial mesenchymal transition and metastatic behavior (50,51). GSEA correlations predict LEF1 pathway activation in the mesenchymal state ( Table 1). Figure 3A shows www.frontiersin.org GSEA plots of genes down-regulated by overexpression of LEF1 in epithelial DLD1 cells (52), positively correlating with the epithelial cell state. Correlations between the H358 and A549 RNA abundance changes and LEF1 gene signatures were statistically significant (p < 0.0001 and FDR q-value <0.02; Table 1A). The top 20 positively and negatively correlated genes were identified and heat mapped for H358 ( Figure 3B; control and TGFβ duplicate samples) and for A549 (Figure 3C; control and TGFβ duplicate samples). The data suggested TCF/LEF1 signaling were up-regulated in EMT-derived mesenchymal lung tumor cells. We asked whether EMT state in H358 or A549 cells might alter the abundance of RNA transcripts encoding components of the Wnt signaling pathway itself, where log2 ratios are shown for duplicate samples comparing 44 Wnt pathway components between mesenchymal and epithelial states ( Figure 4A). Both positively acting and negatively acting Wnt pathway encoded RNAs were observed. The Wnt pathway positive regulators TCF4/TCF7L1/LEF1 and WNT5A/5B were found to be up-regulated in either H358 or A549 mesenchymal cell states, while the negative regulators DKK1 and NDK2 were down-regulated in mesenchymal cell states. Nuclear β-catenin localization has been shown to correlate with neoplastic transformation and cancer patient outcome (53). In normal epithelial cells β-catenin forms a complex with α-catenin and E-cadherin on the inner side of the plasma membrane, important for cell-cell interaction. By immunofluorescence microscopy, we observed that β-catenin was localized at cell-cell-junctions on the inner side of the plasma membrane in the epithelial cell state, as expected ( Figure 4B; top panels). In the mesenchymal cell state, β-catenin staining was lost from the cell membrane, with expected (54) diffuse cytoplasmic and speckled nuclear staining ( Figure 4B; bottom panels). We asked whether signaling through the Wnt/β-catenin pathway might show a change in activity between epithelial and mesenchymal lung states. H358 cells were transfected with TCF promoter-reporter plasmid (sTOP), measuring TCF/LEF dependent transcription, or control (FOP) plasmid with mutated TCF/LEF binding sites. A control TK-renilla plasmid was used to control for transfection efficiency. Triplicate biological experiments were performed, each run in triplicate, which were normalized to the TK-renilla control and averaged ( Figure 4C). We observed the signal from the TCF/LEF reporter was increased Frontiers in Oncology | Molecular and Cellular Oncology   (Table  S5 in Supplementary Material) for TGFβ, stem cell, and KRas pathways.
Gene-set enrichment and pathway prediction analysis of RNAs altered in abundance with mesenchymal trans-differentiation, indicated potential activation of the NFκB pathway, particularly in A549 (Table 1B; Figure 5A). The Hinata NFκB signature was derived from normal keratinocytes over-expressing NFKB1 and RELA (52  down-regulated in mesenchymal-like A375 melanoma cells treated with an NFκB inhibitor (55). The top 20 positively and negatively correlated genes form the Hinata signature correlation were identified and heat mapped for A549 ( Figure 5B). RNAs encoding components of the pathway, notably IL1, TLR5, and TLR9, were markedly increased in the mesenchymal state ( Figure 5C). It is beyond the scope of this study to fully define upstream signaling pathways, which promote the transcriptional and epigenetic changes observed in H358 and A549 models. However, network predictions from RNA sequencing datasets, phosphoproteomic analysis, and kinase activity immunoblots all suggest pronounced activation of ERK1/2 (56, 57) and casein kinase-2 (CK2) and the inhibition of glycogen synthase-3 (GSK3) protein kinase activities in both A549 and H358 mesenchymal states. Phosphorylation of GSK3beta on the kinase activating site Y216 is reduced in H358 mesenchymal cells (data not shown). The inhibition of

Frontiers in Oncology | Molecular and Cellular Oncology
GSK3 has been shown to stabilize Snail, which in turn increases Zeb1 expression, while CK2 kinase activity has been correlated with NFκB activation (58). Wnt5A is a transcriptional target of NFκB, which is increased in the mesenchymal state in both H358 and A549 models. Wnt5A driven TCF/LEF targets are known to include CD44 and Snail, further reinforcing the establishment of the mesenchymal state. Interaction nodes between TCF/LEF and NFκB pathways were mapped using protein-protein interaction data (BioGrid; Figure 5D). The orphan nuclear receptor NR4A1/TR3/Nur77 was identified as an up-regulated common component of the LEF1 and NFκB 1 protein-protein interaction networks. NR4A1 showed coordinate regulation with redox sensitive genes, e.g., TXNDC5, suggesting NR4A1 may regulate oxidative stress as has been observed in pancreas cancer models (59

ALTERED TRANSCRIPTION AND EPIGENETIC NETWORKS
BMI1 is a component of the polycomb repressive complex (PRC1), formed by complexes with heterochromatin adaptor family proteins (e.g., CBX8). PRC1 functions as an inhibitor of transcription during embryogenesis and in neoplastic transformation, as a tumor suppressor. BMI1 is required for Twist mediated EMT (64). Multiple GSEA signatures associated with changes in BMI1 target gene expression were observed (Table 2C), with the majority showing stringent statistical significance (p < 0.05 and FDR qvalue <0.05). The BMI_DN and Wiederschain signatures reflect genes up-regulated the stem-like medulloblastoma line DAOY by knockdown of BMI1 with or without co-knockdown of the polycomb ring finger gene PCGF2/Mel-18 (65). In contrast, the PRC1_BMI signature (66) is derived by knockdown of BMI1 in fibroblasts, cells of mesodermal origin, showed a reversed correlation. The GSEA plots correlating BMI1 knockdown with the H358 and A549 mesenchymal cell state are shown ( Figure 6A). The top 20 positively and negatively correlated genes were identified and heat mapped for H358 ( Figure 6B) and for A549 (Figure 6C), suggesting altered regulation of the BMI1 pathway with EMT in the two mutant-KRas NSCLC models. Interestingly, the GSEA indicated better correlation with inhibition of BMI1 target genes with the mesenchymal cell phenotype. However, direct analysis of up-regulated and down-regulated BMI1 target genes showed enrichment in both H358 and A549 epithelial and mesenchymal cell states (Table S7 in Supplementary Material), which also was observed by GSEA (anti-correlation) at lower statistical stringency. These findings suggest additional factor(s) likely contribute to BMI1 target gene recognition and/or the directionality of regulation between epithelial and mesenchymal cell states. respectively. RUNX3, previously thought to function as a tumor suppressor, has recently been associated with cancer progression (68) and cooperative induction with NFκB of the inflammatory cytokine IL23 (69). In NSCLC overexpression of RUNX2 has been observed in comparisons of tumor and normal tissues and was implicated with poor outcome (70). The overexpression of exogenous RUNX2 also has been shown to promote EMT (71,72), increase migratory and invasive behavior (73), and increase expression of Twist and Slug, both of which are markedly increased with EMT in the A549 background. We examined interaction networks around RUNX2 and RUNX3, where proteins forming direct contacts and potentially establishing signaling complexes, are shown in Figure 7. In H358 RUNX3 RNA expression was correlated with increased expression for potential protein-protein interactors SP110, NOTCH1, DLG4, JUN, KAT2B, HIVEP3, and HDAC5 and decreased expression and potential interaction with EZH2, SUV39H1, RBM14, CREBBP, and SMAD3. This also was associated with an increased SP110 phosphorylation of S256. In A549, Runx2 expression correlated with increased KAT2B, HIVEP3, ETS1, SP110, and HDAC5 expression and decreased expression of STAT3 and FOS. Multiple transcription factors associated with the maintenance of the epithelial state were markedly decreased with EMT in both models. Grainyhead-like 1 (GRHL1) RNA expression is decreased in both mesenchymal models and GRHL2 is reduced in the H358 www.frontiersin.org model. GRHL2 and GRHL1 share~70% homology and GRHL2 has been associated with the maintenance of an epithelial state (74,75), as a target repressed by Zeb1. Similarly transcription factor AP-2 gamma (TFAP2C), a transcription factor associated with estrogen receptor signaling in breast cancer (76), showed reduced expression in the mesenchymal state. Emerging data indicate TFAP2C is important for the normal luminal epithelial differentiation (77) where knockdown promotes mesenchymal transition (78). Conversely, overexpression of the TFAP2C cDNA reduced CD44 expression, and high TFAP2C expression was correlated with response to neoadjuvant chemotherapy (79). EHF, an ETS family transcription factor, also has been associated with an epithelial differentiation program and EHF RNA expression was decreased (74). Interestingly, the epithelial-specific splicing factors ESRP1 and ESRP2 (80) are markedly attenuated with EMT in the H358 model, but their expression is absent in both epithelial and mesenchymal A549 states, again indicating considerable heterogeneity in EMT programs within closely related NSCLC models. Finally, Odd-skipped related 2 (OSR2) is a Smad3/4 downregulated palate and limb developmental gene (81), reduced in the both mesenchymal models. Epithelial mesenchymal transition state was correlated with statistically significant changes in site-specific phosphorylation ( Table 3). Several members of the BAF SWI/SNF complex chromatin remodeling proteins were altered in their pattern of site-specific phosphorylation dependent on EMT state. Notably SMARCC1 showed increased phosphorylation at positions S328 and S330, and ARID1A was increased in phosphorylation at position S696. CBX3/HP1γ is a chromatin remodeling factor implicated in euchromatin silencing in embryonic stem cells (82) and in transcription elongation by RNA polymerase II on heterochromatic genes (83). CBX3/HP1γ has shown overexpression Frontiers in Oncology | Molecular and Cellular Oncology   in NSCLC (84) as compared with normal adjacent tissue, and expression was correlated with poor survival rate [p = 0.02; (85)]. Phosphorylation measurements also suggested EMT state may be regulated in part by increased phosphorylation of CBX3/HP1 at position S95 (p < 0.001) and likely S93 (p = 0.10; not shown) in the central region between chromo and chromoshadow domains. CBX3/HP1 functions as a transcriptional silencer through histone H3 K9 interaction, where phosphorylation at S93, likely by PKA, results in a reduction of CBX3/HP1 mediated silencing and transcription elongation (86). Change in phosphorylation of CBX8, a component of the PRC1 complex, at position S191 within the atrophin-1 domain also was observed ( Table 3). Phosphorylation of HDAC2 was increased at positions S422 and S424 in the mesenchymal state. S422 and S424 are known CK2a sites, associated with inhibition of deacetylase activity (87). S394 is also a CK2a phosphosite increased in H358, decreased in A549. The reported role S394 in activation (88) and inhibition (87) of HDAC2 may be context dependent through as yet undefined factors. The PAF1 complex protein CTR9/SH2BP1, which plays a role in the maintenance of ESC pluripotency showed phosphosite regulation at T925 (A549) and S970 (H358) ( Table 3). Additional transcriptional and epigenetic regulators were found to be modified by phosphorylation, including NCOR2/SMRT, MEF2D, FOXK1, BRD3, and the PRC1 component PHC3. Decreased phosphorylation of the bromodomain protein kinase BAZ1B/WSTF at position S1468, increased phosphorylation of EAF1 at sites S158 and S165 were observed. The phosphorylation of TRIM28/KAP1 S473, likely by CHK2, inhibits co-repression of CDKN1A/p21 (89), and may contribute the observed attenuation of cell cycling. ZC3H8 functions as a repressor of GATA3 and shows increased phosphorylation at position S77 in mesenchymal H358 cells, where GATA3 has been associated with the TGFβ growth suppression response (90).

EARLY TRANSCRIPTIONAL EVENTS IN THE TRANS-DIFFERENTIATION OF MUTANT-KRas NSCLC MODELS
We asked whether key transcription factors, highlighted by GSEA correlations, might show coordinate early expression during mesenchymal transition. RNA samples were isolated from the H358 model following TGFβ induction at 0, 1, 2, 4, 8, 16, 24, 72 h, 7 days (~168 h), 21 days (~500 h), and the long term/steady-state condition (>4500 h), and duplicate samples sequenced. TGFβ1 served as an internal control, where expression was observed within 1 h after addition of doxycycline.
The exchange of cell cycle inhibition with migratory and invasive gene activation programs is a hallmark of metastable mesenchymal trans-differentiation. Reduced RNA expression of cell cycle activators and increased expression of inhibitors was observed and served as a benchmark ( Table 4). Early in the EMT epigenetic reprograming, inhibition of E2F family, MCM family, CCNE1 and CDCA7L gene products, and an increase in BTG2 expression was observed. This was accompanied by reduced phosphorylation of Histone H1B and H1E at positions S18 and T18, respectively ( Table 3). Both sites are cell cycle sensitive, where T18 is a known CDK1 and CDK2 site.
We then examined changes in RNA abundance for encoded transcriptional and epigenetic regulators. We identified 81 transcriptional and epigenetic regulators expressed at early time points in primarily two time bins (6-8 and 18-24 h; Table 5). ELF3, ATOH8, EAF2, MYC, IRX5, CITED2, and ID2 showed relatively rapid repression~8 h after induction, which was continued at later time points. ELF3 encodes an E74-like domain transcription factor, which is epithelial-specific. FOXS1, STAT5A, and HIC1 were increased at the early, continuing on to later time points. By 18-24 h HR, FOXR2, SOX2, RORA, and VSX1 showed a decrease in RNA abundance, while FOXS1, DMBX1, www.frontiersin.org  CDMBX1, CITED4, RORC, VGLL3, RUNX3, and BHLHE40 were markedly increased. The homeodomain transcription factor DMBX1 has been shown to be required for reprograming of induced pluripotent stem cells, and for midbrain development (91). Several transcription regulators showed biphasic expression, where an early decrease in RNA abundance was followed by an increase at the >6 months steady-state time point for the helix-loop-helix dominant negative E-box binding inhibitor ID2 and RAR-related orphan receptor RORA, while the reverse was true for PBX/knotted 1 homeobox 2 (PKNOX2), which was increased at early time points and markedly attenuated by 21 days ( Table 5). Genes from the 6-8 and 18-24 h sets (16,058 genes each) were used for GSEA correlations using oncogene and curated signature datasets. The top ranked gene signatures, divided into functional families and correlating with early and middle time points are shown in Table 6. The LEF1 and NFκB activation signatures and BMI1 and EZH2 inhibition signatures were correlated with the mesenchymal state at early time points. Interestingly, HOXA9 target gene inhibition also was correlated with the mesenchymal state, where crosstalk between HOXA9 and BMI1 can impact cell cycling and senescence programs (92). The positive correlation of the mesenchymal state with HIF1 hypoxia signaling, embryonic stem cell programs, and STAT signaling were observed at the early time points. A mesenchymal correlation with Myc target inhibition also was an early event ( Table 6). Multiple cell cycle signatures correlated the mesenchymal state with cell cycle inhibition at early time points, consistent with previous complete EMT endpoint data (Tables 1 and 2).

DISCUSSION
A detailed investigation of molecular differences contrasting isogenic epithelial and stem-like mesenchymal tumor cell states has been undertaken in isogenic lung adenocarcinoma cell models harboring mutant-KRas. Our observations reinforce the important role that cancer stemness and EMT can have in driving drug resistance in tumor cells (7,8,12,17,93) and highlight the wide diversity of mechanisms (94) that can be used by tumor cells to evade targeted-and chemo-therapies. Agents and combinations successfully targeting mutant-KRas containing tumor cells, in an epithelial or mesenchymal cell context, would have a marked impact in the treatment of NSCLC and particularly pancreas cancer where trans-differentiation is a frequent event (22). An understanding of these mechanisms and molecular contexts will have important implications in driving combinatorial drug therapy in cancer patients in the future.
Mutation of the KRas oncoprotein leads to the interaction and activation of BRaf/CRaf and Mek1/2 kinases and phosphorylation of Erk1/2 on the activation segment. Erk1/2 phosphorylation on the loop TXY motif activated Erk1/2 catalytic activity and leads to downstream cytoplasmic and nuclear pro-survival substrate target phosphorylation. In both A549 and H358 models, a markedly increased Erk active site phosphorylation was observed in the mesenchymal state (Table S3 in Supplementary Material). Similar to the H358 and A549 models described here we have observed TGFβ mediated induction of EMT in mt-EGFR NSCLC adenocarcinoma lines HCC4006 and HCC827 (data not shown). A common feature in both mt-KRas and mt-EGFR models is an elevated activation of the Erk1/2 kinases. These data support a hypothesis that chronic Erk activation may contribute to the initiation of EMT, and once EMT has occurred, Erk is further activated in part through integrin/paxillin/FAK and TGFbetaR1 signaling networks. Thus we suggest that KRas, through Erk, may act as both to prime initiation of EMT and to maintain the mesenchymal state once a transition has occurred. Preliminary RNA interference studies indicate enhanced mesenchymal cell sensitivity to Erk1/2 knockdown. Interestingly, synthetic lethality in multiple mt-KRas backgrounds has been observed with knockdown of Snail2/Slug (95).
Though our principle interest is in identifying key survival targets required for the viability of cells in the mesenchymal state, we considered that druggable transcriptional or epigenetic regulators acting early in the EMT process might serve as useful candidates. Previous studies have highlighted the heterogeneity of distinct EMT programs and the plethora of potential targets, for example, FGFR (57), Axl (96)  changes occurring as cells transition from epithelial to mesenchymal states, we specifically defined networks, which can modulate epigenetic states. The rationale is that pharmacological modulation of epigenetic regulators can alter the expression of many potential mesenchymal drug targets, potentially overcoming the redundancy and heterogeneity of mesenchymal therapeutic targets. Here, we have identified multiple epigenetic regulatory networks, for example, the PRC1 complex, HP1γ, and BAF/Swi-Snf complex, which can contribute to the formation and maintenance of epithelial and mesenchymal tumor cell states. Overlapping transcriptional programs also were observed, for example, the redundant loss of epithelial state transcription (GRHL2, TRAPC2, EHF, and OSR2) and splicing (ESRP1 and ESRP2) factors. This theme of modulation of overlapping or redundant network components was also observed for enhancers of the mesenchymal cell state (e.g., overlapping Zeb1/2, Snail1/2, Twist) and chromatin reprograming machinery. Synthetic lethality and therapeutic reprograming of the mt-KRas NSCLC models can be investigated by future RNAi or CRISPR approaches and/or cDNA overexpression studies, using the transcriptional and epigenetic nodes defined.

ACKNOWLEDGMENTS
We thank Melissa Monaghan, Wei Li, and Jiang Wang for technical assistance.