Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 7-Day Trial for You or Your Team.

Learn More →

Synergistic action of master transcription factors controls epithelial-to-mesenchymal transition

Synergistic action of master transcription factors controls epithelial-to-mesenchymal transition 2514–2527 Nucleic Acids Research, 2016, Vol. 44, No. 6 Published online 28 February 2016 doi: 10.1093/nar/gkw126 Synergistic action of master transcription factors controls epithelial-to-mesenchymal transition 1,2,† 1,† 1,† 2,3 1,4 Hongyuan Chang , Yuwei Liu , Mengzhu Xue , Haiyue Liu , Shaowei Du , 1,2,5,6 1,2,6,* Liwen Zhang and Peng Wang Laboratory of Systems Biology, Shanghai Advanced Research Institute, Chinese Academy of Sciences, 100 Haike Road, Zhangjiang High-Tech Park, Shanghai 201210, PR China, University of Chinese Academy of Sciences, No. 19A Yuquan Road, Beijing 100049, PR China, Shanghai Institute of Biochemistry and Cell Biology, Shanghai Institutes for Biological Sciences, Chinese Academy of Sciences, 320 Yueyang Road, Shanghai 200031, PR China, 4 5 School of Life Sciences, Shanghai University, 333 Nanchen Road, Shanghai 200444, PR China, Shanghai Institute of Materia Medica, Chinese Academy of Sciences, 555 Zu Chong Zhi Road, Zhangjiang High-Tech Park, Shanghai 201203, PR China and School of Life Science and Technology, ShanghaiTech University, 100 Haike Road, Zhangjiang High-Tech Park, Shanghai 201210, PR China Received November 06, 2015; Revised February 07, 2016; Accepted February 20, 2016 metastasis, in which carcinoma cells lose their epithelial ABSTRACT features and acquire invasive properties (1–4). Although Epithelial-to-mesenchymal transition (EMT) is a com- metastasis is the main cause of cancer-related death, cur- plex multistep process in which phenotype switches rent EMT markers have not been consistently associated are mediated by a network of transcription factors with a poor clinical outcome (5). This disparity is partly be- (TFs). Systematic characterization of all dynamic cause the molecular mechanism of EMT is strongly depen- TFs controlling EMT state transitions, especially dent on the cellular context in which it happens (2,3)and the associated gene regulatory networks (GRNs) are often for the intermediate partial-EMT state, represents a only partially characterized. For example, TGF- signaling highly relevant yet largely unexplored task. Here, we can induce EMT in numerous cell lines and during the dif- performed a computational analysis that integrated ferentiation of multiple tissues and organs, but very differ- time-course EMT transcriptomic data with public ent transcriptional responses are observed in different cell cistromic data and identified three synergistic mas- types (6,7) and the underlying mechanisms are largely un- ter TFs (ETS2, HNF4A and JUNB) that regulate the known. Recently, it has been shown that cell-type-specific transition through the partial-EMT state. Overexpres- master transcription factors (TFs) are responsible for di- sion of these regulators predicted a poor clinical out- recting SMADs to gene targets thus can determine cell- come, and their elimination readily abolished TGF- type-specific responses of TGF-  signaling (8). This obser- -induced EMT. Importantly, these factors utilized vation is consistent with the extensively characterized tran- a clique motif, physically interact and their cumula- scriptional regulation of the E-cadherin gene during EMT tive binding generally characterized EMT-associated where more than ten TFs have been reported to repress E- cadherin expression in numerous cellular models (2,3), sug- genes. Furthermore, analyses of H3K27ac ChIP-seq gesting that additional undiscovered factors could regulate data revealed that ETS2, HNF4A and JUNB are asso- EMT. ciated with super-enhancers and the administration Another limitation in our current understanding of EMT of BRD4 inhibitor readily abolished TGF--induced is that the established transcriptional regulators are typ- EMT. These findings have implications for system- ically studied individually. Hence, the synergistic action atic discovery of master EMT regulators and super- among cooperative TFs, which has emerged as a general enhancers as novel targets for controlling metasta- characteristic of enhancers (9), has not been explored in sis. the transcriptional programs of EMT. Moreover, super- enhancers, which are characterized by unusually high cumu- INTRODUCTION lative binding of multiple TFs, show more sensitivity than regular enhancers toward inhibitors of BRD4, a transcrip- Epithelial-to-mesenchymal transition (EMT) is a funda- tional coactivator that is generally involved in enhancer mental developmental process that has been associated with To whom correspondence should be addressed. Tel: +86 21 20350913; Fax: +86 21 20350912; Email: [email protected] These authors contributed equally to this work. C The Author(s) 2016. Published by Oxford University Press on behalf of Nucleic Acids Research. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact [email protected] Nucleic Acids Research, 2016, Vol. 44, No. 6 2515 function (10). Because master TFs that dictate transcrip- Biosystems). The experiments were performed in triplicate, tional programs in numerous cellular models have been as- and the values were normalized to that of GAPDH. For sociated with super-enhancers (11,12), it is of particular in- absolute qPCR, the standard curves were generated using terest to examine whether master TFs regulating EMT are coding regions of SNAI1, SNAI2 and JUNB as described also associated with super-enhancers and consequently, the by Denzler et al. (22). The primer sequences are as follows: utility of BRD4 inhibitors to abolish EMT. ETS2 qPCR fwd, CCCCTGTGGCTAACAGTTACA; In addition to employing complex molecular mecha- ETS2 qPCR rev, AGGTAGCTTTTAAGGCTTGACTC; nisms, EMT also exhibits phenotype dynamics that is char- HNF4A qPCR fwd, CACGGGCAAACACTACGGT; acterized by progressive loss of epithelial hallmarks coupled HNF4A qPCR rev, TTGACCTTCGAGTGCTGATCC; with gradual gain of mesenchymal properties (3,4). Emerg- JUNB qPCR fwd, ACAAACTCCTGAAACCGAGCC; ing evidence suggests the existence of a partial-EMT (P) JUNB qPCR rev, CGAGCCCTGACCAGAAAAGTA; state that serves as an intermediate state between the ep- FOXP1 qPCR fwd, TGGCATCTCATAAACCATCAGC; ithelial state (E) to mesenchymal state (M) transition (13– FOXP1 qPCR rev, GGTCCACTCATCTTCGTCTCAG; 16). The three-state model of EMT suggests that a clear un- CDH1 qPCR fwd, TCAGGCGTCTGTAGAGGCTT; derstanding of the partial-EMT state could provide key in- CDH1 qPCR rev, ATGCACATCCTTCGATAAGACTG; sights into the metastatic process. Thus far, EMT has been CDH2 qPCR fwd, ACAGTGGCCACCTACAAAGG; largely studied by analyzing the beginning epithelial state CDH2 qPCR rev, CCGAGATGGGGTTGATAATG; or the ending mesenchymal state, mainly due to the tremen- SNAI1 qPCR fwd, GGCCCACCTCCAGACCCACT; dous difficulty in isolating partial-EMT cells ( 17). However, SNAI1 qPCR rev, GCGGGGACATCCTGAGCAGC; a number of cellular models, such as TGF--induced EMT, SNAI2 qPCR fwd, TGTGACAAGGAATATGTGAGCC; have been established (18–20). Thus, a time-course analysis SNAI2 qPCR rev, TGAGCCCTCAGATTTGACCTG; of cells undergoing EMT could reveal the essential charac- GAPDH qPCR fwd, GAGTCAACGGATTTGGTCGT; teristics of partial-EMT cells. GAPDH qPCR rev, GATCTCGCTCCTGGAAGATG. In this study, we addressed the hypothesis that the tem- poral activation of multiple cooperative master TFs me- Transient siRNA knockdown diated by super-enhancers might characterize the partial- EMT state. To this aim we conducted a systematic anal- Cells were seeded in six-well plates at a density of 1 × 10 ysis of the time-course gene expression profile of TGF- - cells/well and were cultured until they reached 40% con- induced EMT and of public cistromic data. We identified fluence. Cells were then transfected with a negative control three novel synergistic master TFs that are associated with siRNA (Silencer , Life Technologies) or with siRNAs super-enhancers. The cumulative binding of the master TFs targeting JUNB, ETS2 and HNF4 at a final concentration characterizes EMT-associated genes and the synergistic ac- of 100 nM using Lipofectamine 2000 reagent (Invitrogen) tion of the master regulators is essential for the partial- according to the manufacturer’s instructions. Twenty-four EMT state. Our analysis demonstrated the utility of system- hours after transfection, TGF-1 was added to the cell atic approaches to identify master regulators of EMT and culture to induce EMT. The mRNAs and proteins were suggested that super-enhancer-associated master regulators extracted for analysis of the effects of the siRNAs at 24 may be valuable for the diagnosis and therapeutic interven- h into EMT. The following are the siRNA sequences: tion of metastasis. ETS2, sense 5 -GUCUGGUGAACGUGAAUCUTT- 3 , antisense 5 -AGAUUCACGUUCACCAGACTT-3 ; HNF4A, sense 5 -ACACCACCCUGGAAUUUGATT-3 , MATERIALS AND METHODS antisense 5 -UCAAAUUCCAGGGUGGUGUTT-3 ; Cell culture JUNB, sense 5 -ACAAGGUGAAGACGCUCAATT-3 , antisense 5 -UUGAGCGUCUUCACCUUGUTT-3 . A549 cells were maintained in Dulbecco’s Modified Eagle’s Medium:Nutrient Mixture F-12 (DMEM/F-12) (Corning) supplemented with 10% fetal bovine serum (FBS). To in- Transwell migration and invasion assay duce EMT, A549 cells were seeded 24 h before induction The in vitro cell migration assay was performed using Tran- to reach ∼70% confluence. The induction of EMT was swell chambers (8-m pore size; Costar). A549 cells trans- performed following published protocols (21). Specifically, fected with siRNAs were plated 24 h post-transfection in TGF-1 (Pepro Tech) was added to a final concentration serum-free medium (5 × 10 cells per Transwell). TGF-1 of 4 ng/ml and medium with TGF-1 was replaced every was added to upper chamber at a final concentration of 5 other day. ng/ml. Medium containing 10% FBS in the lower cham- ber served as a chemoattractant. In the case of JQ1, 500 RNA isolation and quantitative real-time PCR nM compound was added to both upper and lower cham- Total RNA was isolated from cells using TRIzol reagent bers 6 h after the addition of TGF-1. After 24 h, the non- (TaKaRa). Reverse transcription was performed using migrating cells were removed from the upper face of the fil- TM the PrimeScript RT Reagent Kit (TaKaRa) and gDNA ters using cotton swabs and the migratory cells located on Eraser (Perfect Real Time) according to the manufacturers’ the lower side of the chamber were stained with crystal vio- instructions. Quantitative real-time PCR (qPCR) was let, air dried, photographed and counted. Images of six ran- performed using the SYBR Premix Ex Taq (Tli RNase H dom fields at 10 × magnification were captured from each Plus) system on an ABI Step One Plus machine (Applied membrane, and the number of migratory cells was counted. 2516 Nucleic Acids Research, 2016, Vol. 44, No. 6 Similar inserts coated with Matrigel were used to determine using the appropriate antibodies followed by detection with TM cell’s invasive potential in the invasion assay. Clean-Blot IP Detection Kit (HRP) (Thermo Fisher Sci- entific). All of the co-IPs were repeated two or more times. Fluorescence microscopy and staining for ETS2, HNF4A, JUNB and SMAD3 ChIP-qPCR 2×10 A549 cells were seeded in a 30-mm confocal dish. The ChIP assay for BRD4, ETS2, HNF4A and JUNB 24 h into TGF--induced EMT, cells were fixed in 4% were performed using HighCell ChIP kit (Diagenode) paraformaldehyde, permeabilized in 0.2% Triton X-100 according to the manufacturer’s instructions. A549 cells and probed with specific primary antibodies for ETS2 were harvested at 48 h into TGF--induced EMT and 10 (Santa Cruz, sc-365666), HNF4A (Santa Cruz, sc-6557), cells were used per ChIP assay. Chromatin was sonicated to JUNB (CST, 3753) or SMAD3 (CST, C67H9). The pri- an average fragment size of 800 bp using a Qsonica Q700 mary antibodies were then detected using matching fluo- sonicator. A fraction (1%) of the sonicated chromatin was rescent secondary antibodies. To detect nuclei, cells were used as ‘input’ DNA and the qPCR results were analyzed co-stained with 4 -6-diamidino-2-phenylindole (DAPI; In- using the Percent Input Method. The antibodies used for vitrogen). Cells were observed on a Leica TCS SP8 Confo- ChIP assay are BRD4 (Bethyl, A301-985A); ETS2 (Santa cal Laser Scanning Microscope. Images were analyzed with cruz, sc-365666); HNF4A (Abcam, ab41898) and JUNB Zeiss ZEN software. (CST, 3753). The primer sequences for qPCR are as follows: EGFR qPCR fwd, GGTTCAAAGAGCAGACGTGG; EGFR qPCR rev, CTGAGAAGCCACGGAAGAGA; Immunoblotting analyses ETS2 qPCR fwd, TCTAGGGAGAGGGAACAGCT; Protein lysates were prepared in the presence of PIC ETS2 qPCR rev, CATTCTCAACACAGCAGCCA; (protease–inhibitor complex) and PMSF. Twenty- HNF4A qPCR fwd, GGTTTTCCCACTACTTCCTGC; microgram aliquots were separated on 8% SDS- HNF4A qPCR rev, GAACCCAGAGCCAGGTGTAT; polyacrylamide electrophoresis gels, and the proteins JUNB qPCR fwd, TTTACAAGGACACGCGCTTC; were transferred onto a polyvinylidene difluoride mem- JUNB qPCR rev, CTTTCCTGGCGTCGTTTCC; VIM brane (Merck Millipore). The membrane was incubated qPCR fwd, GTTTCCTCGTTCCCCTTTGG; VIM qPCR for 1 h in blocking buffer (Tris-buffered saline containing rev, GAATTGCTCGTGGGTTGTGT. 0.1% Tween (TBS-T), and 5% non-fat dry milk) followed by incubation overnight at 4 C with the primary anti- RNA sequencing and data analysis bodies for ETS2 (Santa Cruz, sc-365666), HNF4A (Santa Cruz, sc-6557), JUNB (CST, 3753), FOXP1 (Abcam, The sequencing libraries were constructed according to ab16645), CDH1 (CST, 3195S) and CDH2 (CST, 14215S). the protocol for the Illumina TruSeq Sample preparation After washing with TBS-T, the blot was incubated with kit. Sequencing was performed on the Illumina HiSeq horseradish peroxidase (HRP)-conjugated secondary 2500 sequencer. Library construction and sequencing were antibody and the signals were visualized using an enhanced performed at the Genergy Biotech Co., Ltd. (Shang- chemiluminescence system according to the manufacturer’s hai). Paired-end RNA-seq reads with 101 bp at each end instructions (Kodak). were aligned to the ENSEMBL version 75 transcriptome database using the Bowtie program (23). The levels of gene expression were quantified using RSEM software ( 24). Dif- Immunoprecipitations (IPs) ferential expression analysis was conducted using maSigPro Twenty-four hours into TGF--induced EMT, A549 cells (25), and the significant genes were selected using an FDR were washed, and whole cell lysates were prepared in ly- cutoff value of 0.05. sis buffer containing 20 mM HEPES–KOH (pH 7.9), 150 mM NaCl, 0.5% Nonidet P-40, 1 mM EDTA, and 10% Compiling the list of EMT-related genes glycerol and freshly supplemented with 10 mM NaF, 1 mM Na VO , 1 mM PMSF, and a protease inhibitor mixture We collected the genes with a well-documented role in 3 4 (Sigma). Antibodies directed against ETS2 (Santa Cruz, sc- EMT from several recent reviews (1–3). Those genes were 351), HNF4A (Santa Cruz, sc-6557) or JUNB (CST, 3753) then combined with the four regulators identified in this were added to 1 mg of a whole cell lysate and this mix- study (ETS2, HNF4A, JUNB and FOXP1) and genes ture was incubated for 3 h with rotation at 4 C. Prior to in the QIAGEN EMT PCR Array to generate a non- adding the antibodies, an aliquot was removed for use as the redundant list of 120 genes. The gene list is then partitioned input sample. Pre-equilibrated Protein A (for rabbit anti- into E-markers (overexpressed in epithelial state) and M- body), Protein G (for mouse antibody) or Protein A+G (for markers (overexpressed in partial-EMT or mesenchymal goat antibody) Dynabeads in cell-lysis buffer were added to state). The list of E-markers is as follows: AREG, CDH1, the whole cell lysate/antibody mixtures, and these prepara- COL1A2, CXCL5, ELF3, ELF5, ERBB3, FGFBP1, tions were incubated overnight with rotation at 4 C. The FOXA1, FOXA2, FOXO3A, FZD7, GCLC, GDF15, next day, the Dynabeads were immunomagnetically cap- GSC, IL1RN, KRT19, MITF, MST1R, OCLN, PTP4A1, tured and were washed three times using PBS. The pulled- RAB27B, SIP1, SMAD3, SNAI3, STEAP1, TM4SF18, down proteins were then extracted in sample buffer, sepa- TSPAN13. The list of M-markers is as follows: AH- rated using SDS-PAGE, and subjected to immunoblotting NAK, AKT1, BMP7, CALD1, CAMK2N1, CAV2, CD59, Nucleic Acids Research, 2016, Vol. 44, No. 6 2517 CDH2, COL3A1, COL5A2, CTGF, CTNNB1, DSC2, (±3 kb of the ETS2 peak center) and a 25-bp bin size. The DSP, E47, EGFR, ESR1, ETS2, F11R, FN1, FOXC2, tag-density matrix was created by counting the tags using a FOXD3, FOXF1, FOXP1, FOXQ1, GATA4, GATA6, 3-kb window (±1.5 kb of the ETS2 peak center) and a 25- GNG11, GSK3B, HMGA2, HNF4A, IGFBP4, ILK, bp bin size. Super-enhancers were identified as described by ITGA5, ITGAV, ITGB1, JAG1, JUNB, KLF8, KRT14, Whyte et al. All of the plots were generated using R version KRT7, MAP1B, MMP2, MMP3, MMP9, MSN, NODAL, 3.1.1 software. NOTCH1, NUDT13, PDGFRB, PDK4, PLEK2, PP- PDE2, PRDX1, PRX1, PTK2, RAC1, RGS2, SERPINE1, Statistical and survival analysis SIX1, SMAD2, SNAI1, SNAI2, SOX10, SOX4, SOX9, Survival analysis was performed using the following pub- SPARC, SPP1, STAT3, TCF3, TCF4, TFPI2, TGFB1, lic lung cancer data sets: TCGA Lung Adenocarcinoma, TGFB2, TGFB3, TIMP1, TMEFF1, TMEM132A, TPM1, GSE50081 and GSE31210. In the case of the TCGA RNA- TWIST1, TWIST2, VCAN, VIM, VPS13A, WNT11, seq data, samples overexpressing target gene were identified WNT5A, WNT5B, ZEB1, ZEB2, ZNF703. using the procedure outlined by cBioPortal (29,30). Briefly, the gene expression levels in diploid samples (samples har- Motif enrichment analysis boring no amplification or deletion of the target gene) were We used the Clover algorithm with the default parame- used to calculate the mean value and standard deviation for ter settings for the motif enrichment analysis (26). Human the analyzed gene. In the case of the microarray gene ex- chromosome 20 sequences were used as the background. pression values, the R package mclust was used to identify The enhancer sequences were extracted from the hg19 bimodal expression patterns and the lower-expressing clus- genome using the coordinates derived from the ENCODE ter was used to calculate the mean value and standard de- A549 H3K27ac ChIP-seq data that was also enriched with viation. Samples overexpressing the target gene were then H3K4me1 but not H3K4me3 signals. Enhancers were clas- identified using a z -score cutoff of 2. All of the statistical sified into three categories: proximal enhancers (located tests and the survival analysis were performed using R ver- within 10 kb upstream of transcription start sites), distal en- sion 3.1.1 software. hancers (located greater than 10 kb upstream of transcrip- tion start sites) and enhancers within gene body and were RESULTS analyzed separately with Clover. Time-course gene expression analysis revealed a three-state EMT transcriptome signature DREM dynamic regulatory network analysis To determine the transcriptomic dynamics over the course We used DREM version 2.0 software (27) to model, ana- of EMT, we used an in vitro model in which the human lyze and visualize dynamic GRN during EMT, using the small cell lung cancer cell line A549 underwent EMT when following parameters: Filtering Options = Minimum Ab- treated with TGF- (18). Whereas the control cells main- solute Expression Change 2 (difference from 0); Search Op- tained the characteristics of an epithelial phenotype, A549 tions = Maximum number of paths out of a split: 2; Model cells treated with TGF- progressively developed a spindle- Selection Options = Penalized Likelihood, Node Penalty: like morphology and lost their tight connections (Supple- 40; Expression Scaling Options = Incorporate expression in mentary Figure S1A). Immunoblotting analysis revealed regulator data for TF; Minimum TF expression after scal- the progressive loss of an epithelial marker (E-cadherin) and ing: 5; Default values were set for other parameters. The key gain of a mesenchymal marker (N-cadherin) (Supplemen- regulators were selected according to the split scores (score tary Figure S1B), confirming that TGF-  induced EMT in threshold < 0.001). A549 cells. We then conducted RNA sequencing (RNA-Seq) ex- ChIP-seq data from published sources periments using TGF- treated A549 cells with biologi- cal replicates at eight time points (Supplementary Figure We used published ChIP-seq data from the Gene Expres- S2). Global analysis of the EMT transcriptomic data re- sion Omnibus database for ETS2 (GSE49402), HNF4A vealed 1,633 differentially expressed genes, which fell into (GSE25021 and GSE23436), JUNB (GSE51142), FOXP1 three major clusters (Figure 1A, Supplementary Table S1). (GSE30992 and GSE51142), SMAD3 (GSE41580 and The expression of known epithelial-associated genes, such GSE51510), H3K27ac (GSE36204, GSE29611, GSE17312, as SMAD3, ELF3 and E-cadherin (also known as CDH1), GSE26320, GSE52658, GSE16256, GSE27823, GSE40129 was, as expected, rapidly downregulated upon TGF- and GSE36354) and RNAPII ChIA-PET (GSE33664). treatment. In contrast, the expression of mesenchymal- associated genes, such as VIM, FN1 and N-cadherin (also ChIP-seq data analysis and visualization known as CDH2), was induced early and was progressively We downloaded the processed ChIP-seq data from the upregulated throughout the course of EMT. Whereas those Cistrome data browser (28) and visualized the ChIP-seq two clusters largely recapitulated the classic gene expres- peaks using the UCSC genome browser and the WashU sion changes that occur during EMT, a smaller gene clus- EpiGenome browser. The ChIP-seq peak gene assignment ter (196 genes) showed transient upregulated expression 12– and peak overlap were analyzed using HOMER soft- 36 h into EMT. At 12–36 h into EMT, both epithelial and ware (http://homer.salk. edu/homer/). The heatmap matri- mesenchymal markers were expressed at intermediate levels ces were created by counting the tags using a 6-kb window (Figure 1A), which is consistent with the current definition Color key 2518 Nucleic Acids Research, 2016, Vol. 44, No. 6 A B -log10(Benjamini) 010 20 30 40 Cell cycle SMAD3 Cell cycle process Cell division Organelle organization Chromosome segregation Epithelial Cellular response to stimulus high genes Cellular metabolic process Nitrogen compound metabolic process Response to stress DNA packaging Cell adhesion Multicellular organismal development Anatomical structure development Response to external stimulus -2 Mesenchymal Cell motion high genes Anatomical structure morphogenesis Cell motility Localization of cell Cellular developmental process Regulation of developmental process CDH1 ELF3 SOX4 FN1 -log10(Benjamini) CDH2 0.0 0.5 1.0 1.5 2.0 Regulation of localization VIM Regulation of locomotion Response to chemical stimulus Negative regulation of locomotion Regulation of biological quality Partial-EMT TCF4 Regulation of molecular function high genes Vesicle−mediated transport Positive regulation of locomotion Membrane organization Positive regulation of developmental process TGFB2 Figure 1. Time-course transcriptome analysis reveals dynamics of EMT gene expression program. (A) Heatmap of the time-course gene expression data for A549 cells undergoing TGF--induced EMT. The differentially expressed genes are separated into three clusters. Representative EMT genes are indicated. (B) Gene Ontology terms enriched in the three clusters of genes that were differentially expressed during TGF--induced EMT. of partial-EMT. Overall, the transcriptome of A549 cells ulus’, indicating that our clustering captured the key gene stimulated with TGF- to induce EMT exhibited three tem- expression dynamics of TGF--induced EMT. poral waves of gene expression, early (epithelial-high), inter- mediate (partial-EMT-high) and late (mesenchymal-high). ETS2, HNF4A and JUNB are putative master TFs of the This observation is consistent with both the time-course partial-EMT state EMT transcriptomic results reported by other groups (31) and the current prevailing regulatory model of EMT, the The gene expression dynamics confirmed that A549 cells Cascading Bistable Switches (CBS) model, which predicts progressively transited through three distinctive states dur- three stable states (21). ing TGF--induced EMT. These coordinated transcrip- Next, we analyzed the Gene Ontology (GO) enrichments tomic changes may be the result of differential TF activi- in each cluster (Figure 1B). The epithelial-high genes were ties. In accordance with the results of previous studies (18), highly enriched with GO terms related to cell cycle. Re- we found that the expression of canonical EMT TFs, such duced proliferation is a well-known phenotype associated as TWIST1/2 and ZEB1/2, was not significantly altered with TGF- signaling (6), and the observed enrichment in in A549 cells undergoing TGF--induced EMT (Supple- cell cycle functions for the early-wave genes was consistent mentary Figure S3A). Although we observed upregulated with reports demonstrating a reduced proliferation rate in SNAI1/2 expression, as reported previously (18), the max- cells undergoing EMT (32). In contrast, the top GO terms imal TPM (transcripts per million) values for SNAI1 (5.58 enriched in the partial-EMT-high gene cluster were regu- TPM) and SNAI2 (4.11 TPM), which were confirmed by lation of localization and locomotion, and cell adhesion absolute qPCR (Supplementary Figure S3B), were within was the top GO term associated with the mesenchymal- the bottom 40% of the values for expressed genes accord- high genes. These results are consistent with the increased ing to the definition of Hebenstreit et al. (33), which sug- cell motility observed upon TGF- treatment, another hall- gested that these genes are unlikely to be master TFs for mark of EMT. Despite these differences, all three clusters the EMT of A549 cells. Thus additional unknown factors were enriched with GO terms related to ‘response to stim- may be involved. We therefore calculated the binding mo- tif enrichment of selected TFs, with the following require- 0 h 6 h 12 h 24 h 36 h 48 h 72 h 96 h Partial-EMT Mesenchymal high genes Epithelial high genes high genes Nucleic Acids Research, 2016, Vol. 44, No. 6 2519 ments for each factor: (i) has a known position weight ma- EMT regulators, we assembled a dynamic GRN for EMT trix; (ii) with a maximal TPM of >10; (iii) shows differen- using the DREM (Dynamic Regulatory Events Miner) soft- tial expression during TGF--induced EMT of A549 cells; ware (27). For this analysis, we integrated the time-course (iv) is found within enhancers associated with epithelial- gene expression data with the computationally predicted high, partial-EMT-high and mesenchymal-high genes. We genome-wide binding data for 52 TFs that were differen- applied the TPM cutoff to focus on the TFs with relatively tially expressed during TGF--induced EMT. DREM iden- high expression levels, an intrinsic characteristic of master tified a dynamic GRN with prominent splitting points at TFs that provides advantages in competing for co-factors 6 and 48 h into EMT, which agreed well with the timing to dominate the transcriptional landscape (12). To identify of the putative state transitions from E to P and P to M, potential key regulators responsible for the state transitions, respectively (Figure 3A). Reassuringly, DREM identified we calculated the fold changes in the expression of each fac- ETS2, HNF4A and JUNB as the top-ranked key regula- tor during state transition, using their 0-h level of mRNA tors at the E to P transition. Similar to the E to P split- expression as the reference. We then nominated those TFs ting points, ETS2, HNF4A and JUNB were nominated by with significant binding motif enrichment ( P < 0.05) and DREM as the top-ranked regulators at the P to M transi- fold expression change (absolute fold change ≥4) as the po- tion. In addition to ETS2, HNF4a and JUNB, DREM also tential key regulators (Figure 2A, Supplementary Figure identified known EMT regulators such as SMAD3, ELF3 S4). This analysis confirmed the role of known key regula- and HMGA2 as key TFs operating at the splitting points. tor of EMT, such as SMAD3 and ELF3, in maintaining the Importantly, the expression levels of ETS2, HNF4A and epithelial phenotype. Notably, four TFs, ETS2, HNF4A, JUNB were upregulated at the splitting points compared JUNB and FOXP1, emerged as the top putative regulators to the levels at the 0-h time point, suggesting that they may of the partial-EMT state in proximal enhancers, distal en- play more active roles during state transitions than regu- hancers and enhancers within gene body (Figure 2A, Sup- lators with downregulated expression, such as SMAD3 or plementary Figure S4). Interestingly, whereas the expres- ELF3. Overall, the predictions obtained using the two in- sion level of FOXP1 increased further upon the transition dependent computational approaches (DREM and binding into the mesenchymal state, the expression levels of ETS2, motif enrichment) agreed well, further establishing ETS2, HNF4A and JUNB fell below the 4-fold cutoff at 72 h into HNF4A and JUNB as candidate master regulators of the EMT (Figure 2B). We validated the protein levels of the partial-EMT state. four TFs during EMT (Figure 2C), confirming that ETS2, To gain further insight into the dynamic GRN associated HNF4A and JUNB were transiently overexpressed 12–36 h with EMT, we extracted the core module from the DREM- into EMT and that FOXP1 displayed monotonic upregu- predicted model, retrieving only the interactions connecting lated expression. Consistent with the CBS model of EMT, the putative key regulators and the key markers of EMT the E-low, P-high and M-low expression pattern suggested (CDH1 and CDH2) (Figure 3B). This analysis revealed that ETS2, HNF4A and JUNB are likely regulators of the a highly modular topology, with 22 of 30 possible intra- −16 partial-EMT state and that FOXP1 is a candidate driver of module interactions implemented (modularity P < 2.2 , the mesenchymal state. Fisher’s exact test). This module was characterized by feed- To further validate the predicted binding motif enrich- forward loops and negative feedback loops, key network ment data, we queried the public chromatin immunopre- motifs to generate bimodality, and was consistent with the cipitation followed by sequencing (ChIP-seq) data for the canonical EMT regulatory model (21). Specifically, four of known (SMAD3) and putative (ETS2, HNF4A, JUNB and four putative EMT regulators had an autoregulatory loop, FOXP1) EMT regulators identified in the motif enrich- a common feature of master TFs (12). Finally, the three pu- ment analysis. Because those ChIP-seq analyses were per- tative partial-EMT state regulators used a clique motif, the formed in cells other than A549, we used ENCODE A549 most abundant three-protein interaction motif and one of H3K27ac ChIP-seq data as a reference and retained only the most conserved GRN motifs (34,35), suggesting coop- the TF ChIP-seq peaks that fell within the established A549 erativeness among those regulators. cell active enhancers. We then analyzed the ChIP-seq peaks associated with the 1633 genes that were differentially ex- Validation of the partial-EMT transcriptional regulatory pressed during EMT. Reassuringly, each of the vfi e TFs oc- module cupied the enhancers of a substantial number of the 1633 genes, ranging from 512 genes for JUNB to 784 genes for According to our dynamic EMT GRN model, the three HNF4A (Figure 2D). Furthermore, their target EMT genes putative partial-EMT state regulators, ETS2, HNF4A that were inferred from ChIP-seq data showed a significant and JUNB, autoregulate and positively regulate each −16 overlap (pair-wise overlap P < 10 , Fisher’s exact test), other, repress the expression of the epithelial state mark- which was consistent with the binding motif prediction re- ers and upregulate the expression of mesenchymal state sults (Figure 2D). regulators/markers. To validate the core partial-EMT reg- ulatory module, we first mined public ChIP-seq data for evidence of direct regulation. Binding events of the puta- ETS2, HNF4a and JUNB control the transition through the tive master TFs on themselves were readily supported by partial-EMT state in a dynamic GRN model the ChIP-seq data (Figure 3C–E). To validate the ChIP- To examine whether the same set of master regulators could seq results, we performed ChIP-qPCR on selected sites and be derived using an independent algorithm and to obtain a confirmed that ETS2, HNF4A and JUNB directly colocal- dynamic, systems-level overview of the roles of the putative ized to their own enhancers in A549 cells (Supplementary 2520 Nucleic Acids Research, 2016, Vol. 44, No. 6 Figure 2. Integrative analysis identified ETS2, HNF4A, JUNB and FOXP1 as the putative master TFs of EMT. ( A) Enrichment P -values of TF binding motifs in the proximal enhancers of partial-EMT high genes (x-axis) were plotted against the fold changes in gene expression during EMT (y-axis; 24 h versus 0 h, left panel; 72 h versus 0 h, right panel). Each circle represents a TF. The dotted lines represent the cutoff for the enrichment (P -value < 0.05) and the expression change (absolute fold change ≥ 4). (B) Levels of ETS2, HNF4A, JUNB and FOXP1 mRNA during EMT as estimated using RNA-seq. (C) Immunoblots showing the levels of ETS2, HNF4A, JUNB and FOXP1 proteins during EMT. (D) Venn diagram showing the overlap of the targets of the vfi e EMT TFs derived from published ChIP-seq data. Figures S5A and S6). Moreover, gene targets supported by function during the partial-EMT state, siRNAs were in- the ChIP-seq data and gene targets supported by the motif troduced into A549 cells 24 h prior to TGF- treatment, enrichment data showed extensive overlap (Supplementary and their effects on mRNA and protein expression were Figure S5B), suggesting that the majority of the predicted determined 24 h into TGF--induced EMT. Compared regulations were the result of direct TF-DNA interactions. to the effects of the control siRNAs, siRNAs for ETS2, To validate the predicted functional gene regulation by HNF4A or JUNB significantly reduced the levels of ETS2, the three putative partial-EMT master regulators, we per- HNF4, JUNB, FOXP1 and CDH2 mRNAs and proteins formed siRNA-mediated silencing of each regulator in and markedly increased the levels of CDH1 mRNA and A549 cells undergoing TGF--induced EMT. To focus on protein (Figure 4A–D), which was consistent with the pre- Nucleic Acids Research, 2016, Vol. 44, No. 6 2521 Figure 3. ETS2, HNF4A and JUNB are key regulators in a dynamic EMT Gene Regulatory Network model. (A) DREM output showing a dynamic gene regulatory map of EMT. The y-axis shows the levels of gene expression normalized to those at 0 h. Each line represents a path (cluster) of genes with a similar expression pattern, and nodes represent the hidden states of a hidden Markov model. The green nodes represent splitting points. For each splitting point, TFs with a split score < 0.001 were listed in ranked order of scores and were colored according to the expression level changes (blue for upregulation; red for downregulation). (B) Diagram illustrating the putative core regulatory network derived from DREM analysis. (C) Genome browser representation of ETS2, HNF4A, JUNB and FOXP1 binding events near ETS2 gene. (D) Same as (C) for HNF4A gene. (E) Same as (C) for JUNB gene. dictions of the putative EMT GRN model. Importantly, A unique feature of our predicted partial-EMT regula- knocking down the expression of any one of the ETS2, tory circuit was that three master regulators were identified. HNF4A and JUNB trio had significant effects, suggesting These master regulators displayed significant co-expression that their cooperative activity is critical for EMT. More- and formed a fully connected module connected by positive over, treating A549 cells undergoing EMT with siRNAs for feedback loops and autoregulations (a clique). Those fea- ETS2, HNF4A or JUNB substantially reduced their abil- tures are indicative of cooperative activity. Because a hall- ity to migrate and invade (Figure 4E, F, P < 0.001) with- mark of enhancer activation is the cooperative assembly out a detectable influence on cell’s proliferation or viabil- of a large trans-factor complex consisting of multiple TFs ity (data not shown), supporting the critical roles of these and their co-factors (9), we proceeded to test the hypothesis TFs in modulating key phenotypes of EMT. We next tested that ETS2, HNF4A and JUNB are part of the same trans- whether overexpression of the three master TFs was as- activation complex. Consistent with our ChIP-qPCR data, sociated with a poor clinical outcome by examining three which showed colocalization of ETS2, HNF4A and JUNB independent lung adenocarcinoma data sets. We classified on numerous enhancers (Supplementary Figures S5A and the samples into two categories, those that overexpressed at S6), examining ChIP-seq data revealed that the binding sites least one of the three master TFs and those in which none of occupied by ETS2 were also highly occupied by HNF4A the master TFs was overexpressed. Kaplan-Meier survival and JUNB and that the genomic locations of their maxi- analysis revealed that patients with elevated expression of at mal peak intensities were within 50 bp of each other (Figure least one of the three genes had a significantly worse chance 5A, B). We next evaluated the physical interaction of ETS2, of relapse-free survival (Figure 4G), demonstrating that the HNF4A and JUNB using immunoprecipitation. Compared overexpression of ETS2, HNF4a or JUNB is a robust prog- to the control antibodies, the antibodies specific for each nostic indicator. of the three TFs readily pulled down the other two fac- 2522 Nucleic Acids Research, 2016, Vol. 44, No. 6 Figure 4. ETS2, HNF4A and JUNB play critical roles in partial-EMT. (A) Quantitative reverse transcription polymerase chain reaction results showing the levels of gene expression in A549 cells during TGF--induced EMT after silencing ETS2 expression using a specific siRNA. n = 3; error bars indicate mean ± SD. * P < 0.05; ** P < 0.01, determined using the two-tailed Student’s t-test. (B) Same as (C) for silencing HNF4A. (C) Same as (A) for silencing JUNB. (D) Immunoblotting analysis of the protein abundance of indicated genes in A549 cells undergoing EMT treated with siRNAs targeting ETS2, HNF4A or JUNB. (E) A549 cells undergoing TGF--induced EMT were treated with siRNAs targeting ETS2, HNF4A or JUNB and were subjected to a migration assay. The migratory cells were quantified (bar charts). Scale bars: 100 m. n = 6; error bars indicate mean ± SD. * P < 0.05; ** P < 0.01, determined using the two-tailed Student’s t-test. (F) Same as (E) for the invasion assay. (G) Kaplan–Meier survival analysis based on the ETS2, HNF4A and JUNB expression levels in three independent Lung Adenocarcinoma data sets for disease-free survival. tors (Figure 5C). Moreover, immunofluorescence staining To determine whether super-enhancers also characterize showed that ETS2, HNF4A and JUNB were colocalized in the identified EMT master TFs, we analyzed ENCODE the nucleus of A549 cells 24 h into TGF--induced EMT A549 H3K27ac ChIP-seq data and identified 1050 super- (Figure 5D). Taken together, these data suggested that the enhancers (Figure 6A). Whereas the enhancers associated three factors indeed are part of the same regulatory complex with ETS2 barely missed the cutoff value to be classified as and cooperatively regulate EMT-associated genes. a super-enhancer, HNF4A and JUNB were found to be as- sociated with putative super-enhancers (Figure 6A and B). Super-enhancer characterizes master EMT transcriptional To test whether other key EMT genes were also asso- regulators ciated with super-enhancers, we compiled a list of 120 well-established EMT genes, classified them into epithelial- Recent studies suggested that master TFs dominate the associated markers (E-markers) and mesenchymal- transcriptional programs in various cell types and super- associated markers (M-markers) and evaluated their enhancers drive the expression of master TFs (10–12,36). Nucleic Acids Research, 2016, Vol. 44, No. 6 2523 (Figure 6E, Supplementary Figure S8). The lack of sig- nificant association between super-enhancers in mesenchy- mal cells with M-markers together with the lack of signifi- cant association between super-enhancers in epithelial cells with E-markers raised the interesting possibility that EMT marker genes, regardless of their association with epithelial or mesenchymal phenotype, are always marked with active enhancers. To test this hypothesis, we examined the associ- ation of all 120 EMT marker genes with active enhancers in the 9 tested cell lines and found that the vast majority of both E-marker and M-marker genes are consistently associ- ated with active enhancers (Figure 6F). This observation is consistent with EMT being a ubiquitous biological process and suggested that the core set of EMT-associated genes are generally active in diverse cell types. BRD4 inhibitor abolishes TGF--induced EMT BRD4 is a transcriptional coactivator from the bromod- omain and extraterminal (BET) subfamily of human bro- modomain proteins. Recent studies suggested that BRD4 is generally associated with active enhancers (10). Impor- tantly, JQ1, a recently developed BRD4 inhibitor, prefer- entially abolishes BRD4 functionality at super-enhancers (10), indicating that JQ1 may abolish TGF--induced EMT by inhibiting the master regulators ETS2, HNF4A and JUNB through the disruption of BRD4 functionality. We first analyzed the colocalization of BRD4 to sites occu- pied by EST2, HNF4A and JUNB by performing ChIP- qPCR assays utilizing the same set of primers used in ETS2, HNF4A and JUNB ChIP-qPCRs. Reassuringly, BRD4 binding is readily detected at tested enhancers near ETS2, Figure 5. ETS2, HNF4A and JUNB are synergistic master regulators HNF4A, JUNB, EGFR and VIM, confirming that BRD4 of EMT. (A) Heatmap showing the colocalization of ETS2, HNF4A and colocalized to the same sites occupied by ETS2, HNF4A JUNB at EMT-associated enhancers. The presence of ETS2, HNF4A and JUNB ChIP-seq peaks are displayed within a 6-kb window centered on the and JUNB (Supplementary Figure S9). Interestingly, while ETS2-bound site. (B) Summary plot for the ETS2, HNF4A and JUNB knocking down HNF4A or ETS2 with specific siRNAs sig- ChIP-seq peak enrichment across the ETS2-binding site associated with nificantly reduced BRD4 binding at some of the tested sites, EMT genes. (C) Endogenous association of ETS2, HNF4A and JUNB. At knocking down JUNB consistently reduced BRD4 binding 24 h into TGF--induced EMT, A549 cell lysates were prepared and were at all the tested sites (Supplementary Figure S9). The result immunoprecipitated using the indicated antibody. The presence of associ- ated proteins was then analyzed using immunoblotting. (D) Immunoflu- is consistent with a recent proteome-scale protein-protein orescence staining for ETS2, HNF4A or JUNB in A549 cells stimulated interaction analysis showed that JUNB makes physical con- with TGF- for 24 h. tact with BRD4 (37), suggesting that the master regulator JUNB recruits BRD4 to those sites. To functionally validate the roles of BRD4 inhibitors in intersection with the A549 cell super-enhancer-associated EMT, we treated A549 cells undergoing TGF--induced genes. Interestingly, 28 genes from the list of 120 EMT- EMT with JQ1, which has been extensively utilized to dis- associated genes were associated with the super-enhancers rupt super-enhancers (10,36). We tested the effect of JQ1 of A549 cells (Figure 6A, B, Supplementary Figures S7A– on the expression of key EMT genes by adding JQ1 to −14 C). This enrichment was highly significant ( P < 10 , the cell culture medium at 6 h into TGF--induced EMT Fisher’s exact test), suggesting that super-enhancers may and examined the cells 18 h later. This setup allowed us indeed drive the high expression of numerous key EMT to specifically target ETS2, HNF4A and JUNB associated genes. Importantly, binding motifs for ETS2, HNF4A and super-enhancers at the partial-EMT state without disrupt JUNB were highly enriched in those super-enhancers (P SMAD3 functionality, which is critical for mediating TGF- −6 < 10 ), suggesting that those master regulators drive the  signaling. Indeed, the addition of JQ1 did not interrupt formation of super-enhancers associated with EMT genes SMAD3 nuclear localization following TGF- treatment (Figure 6C and D). (Supplementary Figure S10). The short JQ1 treatment also We next extended the super-enhancer analyses to 8 addi- excluded molecular and functional impacts from JQ1 in- tional cell lines. While all examined cell lines demonstrated duced apoptosis, which occurs after 48 h of JQ1 expo- an enriched association of super-enhancers with EMT sure (38). As expected, the mRNA and protein levels of markers, neither mesenchymal nor epithelial cell lines dis- ETS2, HNF4A and JUNB were significantly reduced com- played biased enrichment toward E-markers or M-markers pared with those of the DMSO controls (Figure 7A, B), 2524 Nucleic Acids Research, 2016, Vol. 44, No. 6 Figure 6. Super-enhancers are associated with master EMT regulators. (A) The distribution of A549 H3K27ac tag intensities revealed 1050 super- enhancers. The red circles indicate EMT genes associated with super-enhancers. (B) Genome browser tracks showing the super-enhancers associated with ETS2, HNF4A, JUNB and SMAD3. The black bars denote super-enhancers. (C) Genome browser tracks showing super-enhancers associated with the EGFR gene and the binding events of ETS2, HNF4A and JUNB around EGFR gene. (D) Table depicting the enrichment of ETS2, HNF4a and JUNB binding motifs in EMT-associated super-enhancers relative to the genomic background. (E) Bar plots showing the enriched association of EMT genes with super-enhancers. Bars represent the observed number of super-enhancer-associated EMT genes and circles represent the expected number of super- enhancer-associated genes assuming no enrichment. * P < 0.001, determined using Fisher’s Exact Test. (F) Bar plots showing the universal association of EMT genes with active enhancers. Bars represent the observed percentages of active-enhancer-associated EMT genes. Nucleic Acids Research, 2016, Vol. 44, No. 6 2525 Figure 7. BRD4 inhibitor abolishes TGF--induced EMT. (A) Quantitative reverse transcription polymerase chain reaction results showing the expression levels of indicated genes in A549 cells undergoing EMT that were treated with various concentrations of JQ1. n = 3; error bars indicate mean ± SD. * P < 0.05; ** P < 0.01, as determined using the two-tailed Student’s t-test. DMSO, dimethylsulphoxide. (B) Immunoblotting analysis of the expression levels of indicated genes in A549 cells undergoing EMT that were treated with various concentrations of JQ1. (C) A549 cells undergoing TGF--induced EMT were treated with 500 nM JQ1 (left panels) or the DMSO control (right panel) and were subjected to a migration assay (top panels) or an invasion assay (bottom panels). (D) The migratory or invasive cells were quantified (bar charts). Scale bars: 100 m. n = 6; error bars indicate mean ± SD. * P < 0.05; ** P < 0.01, determined using the two-tailed Student’s t-test. consistent with the association of these genes with super- of published ChIP-seq data. We successfully recovered and enhancers. Consequently, the expression of FOXP1 and the validated novel partial-EMT regulators using the small cell mesenchymal marker gene (CDH2) regulated by those mas- lung cancer A549 cell line. In principle, our approach can ter TFs was also downregulated. Furthermore, cells under- be extended to any system in which time-course gene expres- going EMT that were treated with JQ1 displayed signifi- sion data is accessible and could be valuable for deciphering cantly reduced abilities to migrate and invade (Figure 7C, EMT mechanisms in other cancers. D, P < 0.01) and thus were losing key EMT characteristics. In sharp contrast to the canonical EMT model, in which a single TF such as SNAI1 of ZEB1 controls the EMT pro- gram, we identified multiple master TFs that synergistically DISCUSSION control the EMT transcriptional program in A549 cells. The master TFs physically interact and demonstrate cumulative EMT has been implicated as an integral part of metasta- binding at enhancers of EMT-associated genes, suggesting sis, and a clear understanding of the underlying molecular that the synergistic action of multiple master TFs could rep- mechanism of the latter is crucial for successful clinical ther- resent an intrinsic property of the EMT transcriptional pro- apies. Unfortunately, EMT is a phenomenon that is highly gram. Similar to the master TFs identified in other cellu- dependent on the cellular context in which it occurs. Despite lar models, the master TFs identified in our study are as- extensive efforts, the regulatory mechanism of EMT in most sociated with super-enhancers. In addition to A549 cells, cell/tissue types is not fully resolved. Our understanding is JQ1 inhibited migration and invasion of multiple breast and particularly poor in the case of partial-EMT, the interme- prostate cancer cell lines (data not shown). Moreover, the diate state that was recently proposed to hold key insights disruption of BRD4 function by JQ1 effectively inhibits tu- into metastasis (16). In this study, we employed a system- mor formation and metastasis in numerous mouse models atic approach to identify all potential key regulators and (39,40). Although it is not possible to deconvolute the role corresponding GRNs of EMT using an integrative analysis of BRD4 inhibitors in EMT using vivo assays owing to the of time-course gene expression data and a large repertoire 2526 Nucleic Acids Research, 2016, Vol. 44, No. 6 massive apoptosis effects from JQ1 administration (36,38), 2. Lamouille,S., Xu,J. and Derynck,R. (2014) Molecular mechanisms of epithelial-mesenchymal transition. Nat. Rev. Mol. Cell Biol., 15, these observations supported that super-enhancers, the crit- 178–196. ical role of which in other aspects of tumorigenesis has been 3. Thiery,J.P., Acloque,H., Huang,R.Y. and Nieto,M.A. (2009) demonstrated in several recent studies, are key players in Epithelial-mesenchymal transitions in development and disease. Cell, EMT and may together represent a general target for coin- 139, 871–890. 4. Yang,J. and Weinberg,R.A. (2008) Epithelial-mesenchymal hibiting multiple hallmarks of cancers. transition: at the crossroads of development and tumor metastasis. Consistent with the mechanism of enhancer formation, Dev. Cell, 14, 818–829. we identified 3 synergistic master TFs regulating TGF- - 5. Garber,K. (2008) Epithelial-to-mesenchymal transition is important induced EMT in A549 cells. Interestingly, each of the three to metastasis, but questions remain. J. Natl. Cancer Inst., 100, 232–233. factors, ETS2, HNF4A and JUNB, has been previously as- 6. Ikushima,H. and Miyazono,K. (2010) TGFbeta signalling: a complex sociated with EMT, albeit in different cellular backgrounds web in cancer progression. Nat. Rev. Cancer, 10, 415–424. (41–43). Aberrant reactivation of lineage specific TFs is a 7. Massague,J. and Gomis,R.R. (2006) The logic of TGFbeta signaling. well-known phenomena during tumorigenesis (44,45). Cur- FEBS Lett., 580, 2811–2820. rently, more than 30 TFs have been implicated in EMT pro- 8. Mullen,A.C., Orlando,D.A., Newman,J.J., Loven,J., Kumar,R.M., Bilodeau,S., Reddy,J., Guenther,M.G., DeKoter,R.P. and grams underlying a variety of tissue remodeling events (1,2). Young,R.A. (2011) Master transcription factors determine It is conceivable that many of them could be reactivated dur- cell-type-specific responses to TGF-beta signaling. Cell, 147, 565–576. ing metastasis. Thus, systematic discovery and characteriza- 9. Spitz,F. and Furlong,E.E. (2012) Transcription factors: from tion of their synergistic interactions may offer substantial enhancer binding to developmental control. Nat. Rev. Genet., 13, insights into the molecular mechanisms underlying EMT 613–626. 10. Loven,J., Hoke,H.A., Lin,C.Y., Lau,A., Orlando,D.A., Vakoc,C.R., and may identify novel prognostic markers as well as new Bradner,J.E., Lee,T.I. and Young,R.A. (2013) Selective inhibition of therapeutic targets for controlling metastasis. tumor oncogenes by disruption of super-enhancers. Cell, 153, 320–334. 11. Whyte,W.A., Orlando,D.A., Hnisz,D., Abraham,B.J., Lin,C.Y., AVAILABILITY OF SUPPORTING DATA Kagey,M.H., Rahl,P.B., Lee,T.I. and Young,R.A. (2013) Master transcription factors and mediator establish super-enhancers at key The data reported in this study have been deposited in GEO cell identity genes. Cell, 153, 307–319. under accession number GSE69667. 12. Hnisz,D., Abraham,B.J., Lee,T.I., Lau,A., Saint-Andre,V., Sigova,A.A., Hoke,H.A. and Young,R.A. (2013) Super-enhancers in the control of cell identity and disease. Cell, 155, 934–947. SUPPLEMENTARY DATA 13. Revenu,C. and Gilmour,D. (2009) EMT 2.0: shaping epithelia through collective migration. Curr. Opin. Genet. Dev., 19, 338–342. Supplementary Data are available at NAR Online. 14. Chao,Y., Wu,Q., Acquafondata,M., Dhir,R. and Wells,A. (2012) Partial mesenchymal to epithelial reverting transition in breast and prostate cancer metastases. Cancer Microenviron., 5, 19–28. ACKNOWLEDGEMENTS 15. Theveneau,E. and Mayor,R. (2012) Neural crest delamination and The authors would like to thank members of the Labora- migration: from epithelium-to-mesenchyme transition to collective cell migration. Dev. Biol., 366, 34–54. tory of Systems Biology for participating in helpful discus- 16. Jordan,N.V., Johnson,G.L. and Abell,A.N. (2011) Tracking the sions and National Center for Protein Science Shanghai for intermediate stages of epithelial-mesenchymal transition in epithelial assistance with confocal microscopy. stem cells and cancer. Cell cycle, 10, 2865–2873. Authors’ contributions: H.C., Y.L., H.L., S.D. and L.Z. 17. Yu,M., Bardia,A., Wittner,B.S., Stott,S.L., Smas,M.E., Ting,D.T., Isakoff,S.J., Ciciliano,J.C., Wells,M.N., Shah,A.M. et al. (2013) performed the experiments. M.X. conducted the com- Circulating breast tumor cells exhibit dynamic changes in epithelial putational analysis. H.C. and H.L. contributed to the and mesenchymal composition. Science, 339, 580–584. computational analysis. P.W. designed and supervised the 18. Liu,J., Hu,G., Chen,D., Gong,A.Y., Soori,G.S., Dobleman,T.J. and study, contributed to the computational analysis and wrote Chen,X.M. (2013) Suppression of SCARA5 by Snail1 is essential for the manuscript. All authors have read and approved the EMT-associated cell migration of A549 cells. Oncogenesis, 2, e73. 19. Zavadil,J. and Bottinger,E.P. (2005) TGF-beta and manuscript. epithelial-to-mesenchymal transitions. Oncogene, 24, 5764–5774. 20. Xu,J., Lamouille,S. and Derynck,R. (2009) TGF-beta-induced epithelial to mesenchymal transition. Cell Res., 19, 156–172. FUNDING 21. Zhang,J.Y., Tian,X.J., Zhang,H., Teng,Y., Li,R.Y., Bai,F., National Natural Science Foundation of China (NSFC) Elankumaran,S. and Xing,J.H. (2014) TGF-beta-induced epithelial-to-mesenchymal transition proceeds through stepwise [31271413]; Science and Technology Commission of Shang- activation of multiple feedback loops. Sci. Signal, 7, ra91. hai Municipality [12DZ1910800]; Chinese Academy of Sci- 22. Denzler,R., Agarwal,V., Stefano,J., Bartel,D.P. and Stoffel,M. (2014) ences ‘Hundred Talent Program’ (to P.W.). Funding for Assessing the ceRNA hypothesis with quantitative measurements of open access charge: National Natural Science Founda- miRNA and target abundance. Mol. Cell, 54, 766–776. tion of China (NSFC) [31271413]; Science and Technology 23. Langmead,B., Trapnell,C., Pop,M. and Salzberg,S.L. (2009) Ultrafast and memory-efficient alignment of short DNA sequences to the Commission of Shanghai Municipality [12DZ1910800]; human genome. Genome Biol., 10, R25. Chinese Academy of Sciences ‘Hundred Talent Program’. 24. Li,B. and Dewey,C.N. (2011) RSEM: accurate transcript Conflict of interest statement. None declared. quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics, 12, 323. 25. Conesa,A., Nueda,M.J., Ferrer,A. and Talon,M. (2006) maSigPro: a REFERENCES method to identify significantly differential expression profiles in time-course microarray experiments. Bioinformatics, 22, 1096–1102. 1. De Craene,B. and Berx,G. (2013) Regulatory networks defining EMT during cancer initiation and progression. Nat. Rev. Cancer, 13, 97–110. Nucleic Acids Research, 2016, Vol. 44, No. 6 2527 26. Frith,M.C., Fu,Y., Yu,L., Chen,J.F., Hansen,U. and Weng,Z. (2004) integrated cellular networks of transcription-regulation and Detection of functional DNA motifs via statistical protein-protein interaction. Proc. Natl. Acad. Sci. U.S.A., 101, over-representation. Nucleic Acids Res., 32, 1372–1381. 5934–5939. 27. Schulz,M.H., Devanny,W.E., Gitter,A., Zhong,S., Ernst,J. and 36. Asangani,I.A., Dommeti,V.L., Wang,X., Malik,R., Cieslik,M., Bar-Joseph,Z. (2012) DREM 2.0: Improved reconstruction of Yang,R., Escara-Wilke,J., Wilder-Romans,K., Dhanireddy,S., dynamic regulatory networks from time-series expression data. BMC Engelke,C. et al. (2014) Therapeutic targeting of BET bromodomain Syst. Biol., 6, 104. proteins in castration-resistant prostate cancer. Nature, 510, 278–282. 28. Liu,T., Ortiz,J.A., Taing,L., Meyer,C.A., Lee,B., Zhang,Y., Shin,H., 37. Hein,M.Y., Hubner,N.C., Poser,I., Cox,J., Nagaraj,N., Toyoda,Y., Wong,S.S., Ma,J., Lei,Y. et al. (2011) Cistrome: an integrative Gak,I.A., Weisswange,I., Mansfeld,J., Buchholz,F. et al. (2015) A platform for transcriptional regulation studies. Genome Biol., 12, R83. human interactome in three quantitative dimensions organized by 29. Cerami,E., Gao,J., Dogrusoz,U., Gross,B.E., Sumer,S.O., stoichiometries and abundances. Cell, 163, 712–723. Aksoy,B.A., Jacobsen,A., Byrne,C.J., Heuer,M.L., Larsson,E. et al. 38. Lockwood,W.W., Zejnullahu,K., Bradner,J.E. and Varmus,H. (2012) (2012) The cBio cancer genomics portal: an open platform for Sensitivity of human lung adenocarcinoma cell lines to targeted exploring multidimensional cancer genomics data. Cancer Discov., 2, inhibition of BET epigenetic signaling proteins. Proc. Natl. Acad. Sci. 401–404. U.S.A., 109, 19408–19413. 30. Gao,J., Aksoy,B.A., Dogrusoz,U., Dresdner,G., Gross,B., 39. Alsarraj,J. and Hunter,K.W. (2012) Bromodomain-Containing Sumer,S.O., Sun,Y., Jacobsen,A., Sinha,R., Larsson,E. et al. (2013) Protein 4: A Dynamic Regulator of Breast Cancer Metastasis Integrative analysis of complex cancer genomics and clinical profiles through Modulation of the Extracellular Matrix. Int. J. Breast using the cBioPortal. Sci. Signal., 6, pl1. Cancer, 2012, 670632. 31. Javaid,S., Zhang,J., Anderssen,E., Black,J.C., Wittner,B.S., 40. Hu,Y., Zhou,J., Ye,F., Xiong,H., Peng,L., Zheng,Z., Xu,F., Cui,M., Tajima,K., Ting,D.T., Smolen,G.A., Zubrowski,M., Desai,R. et al. Wei,C., Wang,X. et al. (2015) BRD4 inhibitor inhibits colorectal (2013) Dynamic chromatin modification sustains cancer growth and metastasis. Int. J. Mol. Sci., 16, 1928–1948. epithelial-mesenchymal transition following inducible expression of 41. Yang,M., Li,S.N., Anjum,K.M., Gui,L.X., Zhu,S.S., Liu,J., Snail-1. Cell Rep., 5, 1679–1689. Chen,J.K., Liu,Q.F., Ye,G.D., Wang,W.J. et al. (2013) A 32. Evdokimova,V., Tognon,C., Ng,T., Ruzanov,P., Melnyk,N., Fink,D., double-negative feedback loop between Wnt-beta-catenin signaling Sorokin,A., Ovchinnikov,L.P., Davicioni,E., Triche,T.J. et al. (2009) and HNF4alpha regulates epithelial-mesenchymal transition in Translational activation of snail1 and other developmentally hepatocellular carcinoma. J. Cell Sci., 126, 5692–5703. regulated transcription factors by YB-1 promotes an 42. Gervasi,M., Bianchi-Smiraglia,A., Cummings,M., Zheng,Q., epithelial-mesenchymal transition. Cancer Cell, 15, 402–415. Wang,D., Liu,S. and Bakin,A.V. (2012) JunB contributes to Id2 33. Hebenstreit,D., Fang,M., Gu,M., Charoensawan,V., van repression and the epithelial-mesenchymal transition in response to Oudenaarden,A. and Teichmann,S.A. (2011) RNA sequencing transforming growth factor-beta. J. Cell Biol., 196, 589–603. reveals two major classes of gene expression levels in metazoan cells. 43. Polydorou,C. and Georgiades,P. (2013) Ets2-dependent trophoblast Mol. Syst. Biol., 7, 497. signalling is required for gastrulation progression after primitive 34. Stergachis,A.B., Neph,S., Sandstrom,R., Haugen,E., Reynolds,A.P., streak initiation. Nat. Commun., 4, 1658. 44. Darnell,J.E. Jr (2002) Transcription factors as targets for cancer Zhang,M., Byron,R., Canfield,T., Stelhing-Sun,S., Lee,K. et al. therapy. Nat. Rev. Cancer, 2, 740–749. (2014) Conservation of trans-acting circuitry during mammalian 45. Vaquerizas,J.M., Kummerfeld,S.K., Teichmann,S.A. and regulatory evolution. Nature, 515, 365–370. 35. Yeger-Lotem,E., Sattath,S., Kashtan,N., Itzkovitz,S., Milo,R., Luscombe,N.M. (2009) A census of human transcription factors: Pinter,R.Y., Alon,U. and Margalit,H. (2004) Network motifs in function, expression and evolution. Nat. Rev. Genet., 10, 252–263. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Nucleic Acids Research Oxford University Press

Synergistic action of master transcription factors controls epithelial-to-mesenchymal transition

Download PDF
14 pages

Loading next page...
 
/lp/oxford-university-press/synergistic-action-of-master-transcription-factors-controls-epithelial-HofIlGAFUI

References (46)

Publisher
Oxford University Press
Copyright
The Author(s) 2016. Published by Oxford University Press on behalf of Nucleic Acids Research.
ISSN
0305-1048
eISSN
1362-4962
DOI
10.1093/nar/gkw126
pmid
26926107
Publisher site
See Article on Publisher Site

Abstract

2514–2527 Nucleic Acids Research, 2016, Vol. 44, No. 6 Published online 28 February 2016 doi: 10.1093/nar/gkw126 Synergistic action of master transcription factors controls epithelial-to-mesenchymal transition 1,2,† 1,† 1,† 2,3 1,4 Hongyuan Chang , Yuwei Liu , Mengzhu Xue , Haiyue Liu , Shaowei Du , 1,2,5,6 1,2,6,* Liwen Zhang and Peng Wang Laboratory of Systems Biology, Shanghai Advanced Research Institute, Chinese Academy of Sciences, 100 Haike Road, Zhangjiang High-Tech Park, Shanghai 201210, PR China, University of Chinese Academy of Sciences, No. 19A Yuquan Road, Beijing 100049, PR China, Shanghai Institute of Biochemistry and Cell Biology, Shanghai Institutes for Biological Sciences, Chinese Academy of Sciences, 320 Yueyang Road, Shanghai 200031, PR China, 4 5 School of Life Sciences, Shanghai University, 333 Nanchen Road, Shanghai 200444, PR China, Shanghai Institute of Materia Medica, Chinese Academy of Sciences, 555 Zu Chong Zhi Road, Zhangjiang High-Tech Park, Shanghai 201203, PR China and School of Life Science and Technology, ShanghaiTech University, 100 Haike Road, Zhangjiang High-Tech Park, Shanghai 201210, PR China Received November 06, 2015; Revised February 07, 2016; Accepted February 20, 2016 metastasis, in which carcinoma cells lose their epithelial ABSTRACT features and acquire invasive properties (1–4). Although Epithelial-to-mesenchymal transition (EMT) is a com- metastasis is the main cause of cancer-related death, cur- plex multistep process in which phenotype switches rent EMT markers have not been consistently associated are mediated by a network of transcription factors with a poor clinical outcome (5). This disparity is partly be- (TFs). Systematic characterization of all dynamic cause the molecular mechanism of EMT is strongly depen- TFs controlling EMT state transitions, especially dent on the cellular context in which it happens (2,3)and the associated gene regulatory networks (GRNs) are often for the intermediate partial-EMT state, represents a only partially characterized. For example, TGF- signaling highly relevant yet largely unexplored task. Here, we can induce EMT in numerous cell lines and during the dif- performed a computational analysis that integrated ferentiation of multiple tissues and organs, but very differ- time-course EMT transcriptomic data with public ent transcriptional responses are observed in different cell cistromic data and identified three synergistic mas- types (6,7) and the underlying mechanisms are largely un- ter TFs (ETS2, HNF4A and JUNB) that regulate the known. Recently, it has been shown that cell-type-specific transition through the partial-EMT state. Overexpres- master transcription factors (TFs) are responsible for di- sion of these regulators predicted a poor clinical out- recting SMADs to gene targets thus can determine cell- come, and their elimination readily abolished TGF- type-specific responses of TGF-  signaling (8). This obser- -induced EMT. Importantly, these factors utilized vation is consistent with the extensively characterized tran- a clique motif, physically interact and their cumula- scriptional regulation of the E-cadherin gene during EMT tive binding generally characterized EMT-associated where more than ten TFs have been reported to repress E- cadherin expression in numerous cellular models (2,3), sug- genes. Furthermore, analyses of H3K27ac ChIP-seq gesting that additional undiscovered factors could regulate data revealed that ETS2, HNF4A and JUNB are asso- EMT. ciated with super-enhancers and the administration Another limitation in our current understanding of EMT of BRD4 inhibitor readily abolished TGF--induced is that the established transcriptional regulators are typ- EMT. These findings have implications for system- ically studied individually. Hence, the synergistic action atic discovery of master EMT regulators and super- among cooperative TFs, which has emerged as a general enhancers as novel targets for controlling metasta- characteristic of enhancers (9), has not been explored in sis. the transcriptional programs of EMT. Moreover, super- enhancers, which are characterized by unusually high cumu- INTRODUCTION lative binding of multiple TFs, show more sensitivity than regular enhancers toward inhibitors of BRD4, a transcrip- Epithelial-to-mesenchymal transition (EMT) is a funda- tional coactivator that is generally involved in enhancer mental developmental process that has been associated with To whom correspondence should be addressed. Tel: +86 21 20350913; Fax: +86 21 20350912; Email: [email protected] These authors contributed equally to this work. C The Author(s) 2016. Published by Oxford University Press on behalf of Nucleic Acids Research. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact [email protected] Nucleic Acids Research, 2016, Vol. 44, No. 6 2515 function (10). Because master TFs that dictate transcrip- Biosystems). The experiments were performed in triplicate, tional programs in numerous cellular models have been as- and the values were normalized to that of GAPDH. For sociated with super-enhancers (11,12), it is of particular in- absolute qPCR, the standard curves were generated using terest to examine whether master TFs regulating EMT are coding regions of SNAI1, SNAI2 and JUNB as described also associated with super-enhancers and consequently, the by Denzler et al. (22). The primer sequences are as follows: utility of BRD4 inhibitors to abolish EMT. ETS2 qPCR fwd, CCCCTGTGGCTAACAGTTACA; In addition to employing complex molecular mecha- ETS2 qPCR rev, AGGTAGCTTTTAAGGCTTGACTC; nisms, EMT also exhibits phenotype dynamics that is char- HNF4A qPCR fwd, CACGGGCAAACACTACGGT; acterized by progressive loss of epithelial hallmarks coupled HNF4A qPCR rev, TTGACCTTCGAGTGCTGATCC; with gradual gain of mesenchymal properties (3,4). Emerg- JUNB qPCR fwd, ACAAACTCCTGAAACCGAGCC; ing evidence suggests the existence of a partial-EMT (P) JUNB qPCR rev, CGAGCCCTGACCAGAAAAGTA; state that serves as an intermediate state between the ep- FOXP1 qPCR fwd, TGGCATCTCATAAACCATCAGC; ithelial state (E) to mesenchymal state (M) transition (13– FOXP1 qPCR rev, GGTCCACTCATCTTCGTCTCAG; 16). The three-state model of EMT suggests that a clear un- CDH1 qPCR fwd, TCAGGCGTCTGTAGAGGCTT; derstanding of the partial-EMT state could provide key in- CDH1 qPCR rev, ATGCACATCCTTCGATAAGACTG; sights into the metastatic process. Thus far, EMT has been CDH2 qPCR fwd, ACAGTGGCCACCTACAAAGG; largely studied by analyzing the beginning epithelial state CDH2 qPCR rev, CCGAGATGGGGTTGATAATG; or the ending mesenchymal state, mainly due to the tremen- SNAI1 qPCR fwd, GGCCCACCTCCAGACCCACT; dous difficulty in isolating partial-EMT cells ( 17). However, SNAI1 qPCR rev, GCGGGGACATCCTGAGCAGC; a number of cellular models, such as TGF--induced EMT, SNAI2 qPCR fwd, TGTGACAAGGAATATGTGAGCC; have been established (18–20). Thus, a time-course analysis SNAI2 qPCR rev, TGAGCCCTCAGATTTGACCTG; of cells undergoing EMT could reveal the essential charac- GAPDH qPCR fwd, GAGTCAACGGATTTGGTCGT; teristics of partial-EMT cells. GAPDH qPCR rev, GATCTCGCTCCTGGAAGATG. In this study, we addressed the hypothesis that the tem- poral activation of multiple cooperative master TFs me- Transient siRNA knockdown diated by super-enhancers might characterize the partial- EMT state. To this aim we conducted a systematic anal- Cells were seeded in six-well plates at a density of 1 × 10 ysis of the time-course gene expression profile of TGF- - cells/well and were cultured until they reached 40% con- induced EMT and of public cistromic data. We identified fluence. Cells were then transfected with a negative control three novel synergistic master TFs that are associated with siRNA (Silencer , Life Technologies) or with siRNAs super-enhancers. The cumulative binding of the master TFs targeting JUNB, ETS2 and HNF4 at a final concentration characterizes EMT-associated genes and the synergistic ac- of 100 nM using Lipofectamine 2000 reagent (Invitrogen) tion of the master regulators is essential for the partial- according to the manufacturer’s instructions. Twenty-four EMT state. Our analysis demonstrated the utility of system- hours after transfection, TGF-1 was added to the cell atic approaches to identify master regulators of EMT and culture to induce EMT. The mRNAs and proteins were suggested that super-enhancer-associated master regulators extracted for analysis of the effects of the siRNAs at 24 may be valuable for the diagnosis and therapeutic interven- h into EMT. The following are the siRNA sequences: tion of metastasis. ETS2, sense 5 -GUCUGGUGAACGUGAAUCUTT- 3 , antisense 5 -AGAUUCACGUUCACCAGACTT-3 ; HNF4A, sense 5 -ACACCACCCUGGAAUUUGATT-3 , MATERIALS AND METHODS antisense 5 -UCAAAUUCCAGGGUGGUGUTT-3 ; Cell culture JUNB, sense 5 -ACAAGGUGAAGACGCUCAATT-3 , antisense 5 -UUGAGCGUCUUCACCUUGUTT-3 . A549 cells were maintained in Dulbecco’s Modified Eagle’s Medium:Nutrient Mixture F-12 (DMEM/F-12) (Corning) supplemented with 10% fetal bovine serum (FBS). To in- Transwell migration and invasion assay duce EMT, A549 cells were seeded 24 h before induction The in vitro cell migration assay was performed using Tran- to reach ∼70% confluence. The induction of EMT was swell chambers (8-m pore size; Costar). A549 cells trans- performed following published protocols (21). Specifically, fected with siRNAs were plated 24 h post-transfection in TGF-1 (Pepro Tech) was added to a final concentration serum-free medium (5 × 10 cells per Transwell). TGF-1 of 4 ng/ml and medium with TGF-1 was replaced every was added to upper chamber at a final concentration of 5 other day. ng/ml. Medium containing 10% FBS in the lower cham- ber served as a chemoattractant. In the case of JQ1, 500 RNA isolation and quantitative real-time PCR nM compound was added to both upper and lower cham- Total RNA was isolated from cells using TRIzol reagent bers 6 h after the addition of TGF-1. After 24 h, the non- (TaKaRa). Reverse transcription was performed using migrating cells were removed from the upper face of the fil- TM the PrimeScript RT Reagent Kit (TaKaRa) and gDNA ters using cotton swabs and the migratory cells located on Eraser (Perfect Real Time) according to the manufacturers’ the lower side of the chamber were stained with crystal vio- instructions. Quantitative real-time PCR (qPCR) was let, air dried, photographed and counted. Images of six ran- performed using the SYBR Premix Ex Taq (Tli RNase H dom fields at 10 × magnification were captured from each Plus) system on an ABI Step One Plus machine (Applied membrane, and the number of migratory cells was counted. 2516 Nucleic Acids Research, 2016, Vol. 44, No. 6 Similar inserts coated with Matrigel were used to determine using the appropriate antibodies followed by detection with TM cell’s invasive potential in the invasion assay. Clean-Blot IP Detection Kit (HRP) (Thermo Fisher Sci- entific). All of the co-IPs were repeated two or more times. Fluorescence microscopy and staining for ETS2, HNF4A, JUNB and SMAD3 ChIP-qPCR 2×10 A549 cells were seeded in a 30-mm confocal dish. The ChIP assay for BRD4, ETS2, HNF4A and JUNB 24 h into TGF--induced EMT, cells were fixed in 4% were performed using HighCell ChIP kit (Diagenode) paraformaldehyde, permeabilized in 0.2% Triton X-100 according to the manufacturer’s instructions. A549 cells and probed with specific primary antibodies for ETS2 were harvested at 48 h into TGF--induced EMT and 10 (Santa Cruz, sc-365666), HNF4A (Santa Cruz, sc-6557), cells were used per ChIP assay. Chromatin was sonicated to JUNB (CST, 3753) or SMAD3 (CST, C67H9). The pri- an average fragment size of 800 bp using a Qsonica Q700 mary antibodies were then detected using matching fluo- sonicator. A fraction (1%) of the sonicated chromatin was rescent secondary antibodies. To detect nuclei, cells were used as ‘input’ DNA and the qPCR results were analyzed co-stained with 4 -6-diamidino-2-phenylindole (DAPI; In- using the Percent Input Method. The antibodies used for vitrogen). Cells were observed on a Leica TCS SP8 Confo- ChIP assay are BRD4 (Bethyl, A301-985A); ETS2 (Santa cal Laser Scanning Microscope. Images were analyzed with cruz, sc-365666); HNF4A (Abcam, ab41898) and JUNB Zeiss ZEN software. (CST, 3753). The primer sequences for qPCR are as follows: EGFR qPCR fwd, GGTTCAAAGAGCAGACGTGG; EGFR qPCR rev, CTGAGAAGCCACGGAAGAGA; Immunoblotting analyses ETS2 qPCR fwd, TCTAGGGAGAGGGAACAGCT; Protein lysates were prepared in the presence of PIC ETS2 qPCR rev, CATTCTCAACACAGCAGCCA; (protease–inhibitor complex) and PMSF. Twenty- HNF4A qPCR fwd, GGTTTTCCCACTACTTCCTGC; microgram aliquots were separated on 8% SDS- HNF4A qPCR rev, GAACCCAGAGCCAGGTGTAT; polyacrylamide electrophoresis gels, and the proteins JUNB qPCR fwd, TTTACAAGGACACGCGCTTC; were transferred onto a polyvinylidene difluoride mem- JUNB qPCR rev, CTTTCCTGGCGTCGTTTCC; VIM brane (Merck Millipore). The membrane was incubated qPCR fwd, GTTTCCTCGTTCCCCTTTGG; VIM qPCR for 1 h in blocking buffer (Tris-buffered saline containing rev, GAATTGCTCGTGGGTTGTGT. 0.1% Tween (TBS-T), and 5% non-fat dry milk) followed by incubation overnight at 4 C with the primary anti- RNA sequencing and data analysis bodies for ETS2 (Santa Cruz, sc-365666), HNF4A (Santa Cruz, sc-6557), JUNB (CST, 3753), FOXP1 (Abcam, The sequencing libraries were constructed according to ab16645), CDH1 (CST, 3195S) and CDH2 (CST, 14215S). the protocol for the Illumina TruSeq Sample preparation After washing with TBS-T, the blot was incubated with kit. Sequencing was performed on the Illumina HiSeq horseradish peroxidase (HRP)-conjugated secondary 2500 sequencer. Library construction and sequencing were antibody and the signals were visualized using an enhanced performed at the Genergy Biotech Co., Ltd. (Shang- chemiluminescence system according to the manufacturer’s hai). Paired-end RNA-seq reads with 101 bp at each end instructions (Kodak). were aligned to the ENSEMBL version 75 transcriptome database using the Bowtie program (23). The levels of gene expression were quantified using RSEM software ( 24). Dif- Immunoprecipitations (IPs) ferential expression analysis was conducted using maSigPro Twenty-four hours into TGF--induced EMT, A549 cells (25), and the significant genes were selected using an FDR were washed, and whole cell lysates were prepared in ly- cutoff value of 0.05. sis buffer containing 20 mM HEPES–KOH (pH 7.9), 150 mM NaCl, 0.5% Nonidet P-40, 1 mM EDTA, and 10% Compiling the list of EMT-related genes glycerol and freshly supplemented with 10 mM NaF, 1 mM Na VO , 1 mM PMSF, and a protease inhibitor mixture We collected the genes with a well-documented role in 3 4 (Sigma). Antibodies directed against ETS2 (Santa Cruz, sc- EMT from several recent reviews (1–3). Those genes were 351), HNF4A (Santa Cruz, sc-6557) or JUNB (CST, 3753) then combined with the four regulators identified in this were added to 1 mg of a whole cell lysate and this mix- study (ETS2, HNF4A, JUNB and FOXP1) and genes ture was incubated for 3 h with rotation at 4 C. Prior to in the QIAGEN EMT PCR Array to generate a non- adding the antibodies, an aliquot was removed for use as the redundant list of 120 genes. The gene list is then partitioned input sample. Pre-equilibrated Protein A (for rabbit anti- into E-markers (overexpressed in epithelial state) and M- body), Protein G (for mouse antibody) or Protein A+G (for markers (overexpressed in partial-EMT or mesenchymal goat antibody) Dynabeads in cell-lysis buffer were added to state). The list of E-markers is as follows: AREG, CDH1, the whole cell lysate/antibody mixtures, and these prepara- COL1A2, CXCL5, ELF3, ELF5, ERBB3, FGFBP1, tions were incubated overnight with rotation at 4 C. The FOXA1, FOXA2, FOXO3A, FZD7, GCLC, GDF15, next day, the Dynabeads were immunomagnetically cap- GSC, IL1RN, KRT19, MITF, MST1R, OCLN, PTP4A1, tured and were washed three times using PBS. The pulled- RAB27B, SIP1, SMAD3, SNAI3, STEAP1, TM4SF18, down proteins were then extracted in sample buffer, sepa- TSPAN13. The list of M-markers is as follows: AH- rated using SDS-PAGE, and subjected to immunoblotting NAK, AKT1, BMP7, CALD1, CAMK2N1, CAV2, CD59, Nucleic Acids Research, 2016, Vol. 44, No. 6 2517 CDH2, COL3A1, COL5A2, CTGF, CTNNB1, DSC2, (±3 kb of the ETS2 peak center) and a 25-bp bin size. The DSP, E47, EGFR, ESR1, ETS2, F11R, FN1, FOXC2, tag-density matrix was created by counting the tags using a FOXD3, FOXF1, FOXP1, FOXQ1, GATA4, GATA6, 3-kb window (±1.5 kb of the ETS2 peak center) and a 25- GNG11, GSK3B, HMGA2, HNF4A, IGFBP4, ILK, bp bin size. Super-enhancers were identified as described by ITGA5, ITGAV, ITGB1, JAG1, JUNB, KLF8, KRT14, Whyte et al. All of the plots were generated using R version KRT7, MAP1B, MMP2, MMP3, MMP9, MSN, NODAL, 3.1.1 software. NOTCH1, NUDT13, PDGFRB, PDK4, PLEK2, PP- PDE2, PRDX1, PRX1, PTK2, RAC1, RGS2, SERPINE1, Statistical and survival analysis SIX1, SMAD2, SNAI1, SNAI2, SOX10, SOX4, SOX9, Survival analysis was performed using the following pub- SPARC, SPP1, STAT3, TCF3, TCF4, TFPI2, TGFB1, lic lung cancer data sets: TCGA Lung Adenocarcinoma, TGFB2, TGFB3, TIMP1, TMEFF1, TMEM132A, TPM1, GSE50081 and GSE31210. In the case of the TCGA RNA- TWIST1, TWIST2, VCAN, VIM, VPS13A, WNT11, seq data, samples overexpressing target gene were identified WNT5A, WNT5B, ZEB1, ZEB2, ZNF703. using the procedure outlined by cBioPortal (29,30). Briefly, the gene expression levels in diploid samples (samples har- Motif enrichment analysis boring no amplification or deletion of the target gene) were We used the Clover algorithm with the default parame- used to calculate the mean value and standard deviation for ter settings for the motif enrichment analysis (26). Human the analyzed gene. In the case of the microarray gene ex- chromosome 20 sequences were used as the background. pression values, the R package mclust was used to identify The enhancer sequences were extracted from the hg19 bimodal expression patterns and the lower-expressing clus- genome using the coordinates derived from the ENCODE ter was used to calculate the mean value and standard de- A549 H3K27ac ChIP-seq data that was also enriched with viation. Samples overexpressing the target gene were then H3K4me1 but not H3K4me3 signals. Enhancers were clas- identified using a z -score cutoff of 2. All of the statistical sified into three categories: proximal enhancers (located tests and the survival analysis were performed using R ver- within 10 kb upstream of transcription start sites), distal en- sion 3.1.1 software. hancers (located greater than 10 kb upstream of transcrip- tion start sites) and enhancers within gene body and were RESULTS analyzed separately with Clover. Time-course gene expression analysis revealed a three-state EMT transcriptome signature DREM dynamic regulatory network analysis To determine the transcriptomic dynamics over the course We used DREM version 2.0 software (27) to model, ana- of EMT, we used an in vitro model in which the human lyze and visualize dynamic GRN during EMT, using the small cell lung cancer cell line A549 underwent EMT when following parameters: Filtering Options = Minimum Ab- treated with TGF- (18). Whereas the control cells main- solute Expression Change 2 (difference from 0); Search Op- tained the characteristics of an epithelial phenotype, A549 tions = Maximum number of paths out of a split: 2; Model cells treated with TGF- progressively developed a spindle- Selection Options = Penalized Likelihood, Node Penalty: like morphology and lost their tight connections (Supple- 40; Expression Scaling Options = Incorporate expression in mentary Figure S1A). Immunoblotting analysis revealed regulator data for TF; Minimum TF expression after scal- the progressive loss of an epithelial marker (E-cadherin) and ing: 5; Default values were set for other parameters. The key gain of a mesenchymal marker (N-cadherin) (Supplemen- regulators were selected according to the split scores (score tary Figure S1B), confirming that TGF-  induced EMT in threshold < 0.001). A549 cells. We then conducted RNA sequencing (RNA-Seq) ex- ChIP-seq data from published sources periments using TGF- treated A549 cells with biologi- cal replicates at eight time points (Supplementary Figure We used published ChIP-seq data from the Gene Expres- S2). Global analysis of the EMT transcriptomic data re- sion Omnibus database for ETS2 (GSE49402), HNF4A vealed 1,633 differentially expressed genes, which fell into (GSE25021 and GSE23436), JUNB (GSE51142), FOXP1 three major clusters (Figure 1A, Supplementary Table S1). (GSE30992 and GSE51142), SMAD3 (GSE41580 and The expression of known epithelial-associated genes, such GSE51510), H3K27ac (GSE36204, GSE29611, GSE17312, as SMAD3, ELF3 and E-cadherin (also known as CDH1), GSE26320, GSE52658, GSE16256, GSE27823, GSE40129 was, as expected, rapidly downregulated upon TGF- and GSE36354) and RNAPII ChIA-PET (GSE33664). treatment. In contrast, the expression of mesenchymal- associated genes, such as VIM, FN1 and N-cadherin (also ChIP-seq data analysis and visualization known as CDH2), was induced early and was progressively We downloaded the processed ChIP-seq data from the upregulated throughout the course of EMT. Whereas those Cistrome data browser (28) and visualized the ChIP-seq two clusters largely recapitulated the classic gene expres- peaks using the UCSC genome browser and the WashU sion changes that occur during EMT, a smaller gene clus- EpiGenome browser. The ChIP-seq peak gene assignment ter (196 genes) showed transient upregulated expression 12– and peak overlap were analyzed using HOMER soft- 36 h into EMT. At 12–36 h into EMT, both epithelial and ware (http://homer.salk. edu/homer/). The heatmap matri- mesenchymal markers were expressed at intermediate levels ces were created by counting the tags using a 6-kb window (Figure 1A), which is consistent with the current definition Color key 2518 Nucleic Acids Research, 2016, Vol. 44, No. 6 A B -log10(Benjamini) 010 20 30 40 Cell cycle SMAD3 Cell cycle process Cell division Organelle organization Chromosome segregation Epithelial Cellular response to stimulus high genes Cellular metabolic process Nitrogen compound metabolic process Response to stress DNA packaging Cell adhesion Multicellular organismal development Anatomical structure development Response to external stimulus -2 Mesenchymal Cell motion high genes Anatomical structure morphogenesis Cell motility Localization of cell Cellular developmental process Regulation of developmental process CDH1 ELF3 SOX4 FN1 -log10(Benjamini) CDH2 0.0 0.5 1.0 1.5 2.0 Regulation of localization VIM Regulation of locomotion Response to chemical stimulus Negative regulation of locomotion Regulation of biological quality Partial-EMT TCF4 Regulation of molecular function high genes Vesicle−mediated transport Positive regulation of locomotion Membrane organization Positive regulation of developmental process TGFB2 Figure 1. Time-course transcriptome analysis reveals dynamics of EMT gene expression program. (A) Heatmap of the time-course gene expression data for A549 cells undergoing TGF--induced EMT. The differentially expressed genes are separated into three clusters. Representative EMT genes are indicated. (B) Gene Ontology terms enriched in the three clusters of genes that were differentially expressed during TGF--induced EMT. of partial-EMT. Overall, the transcriptome of A549 cells ulus’, indicating that our clustering captured the key gene stimulated with TGF- to induce EMT exhibited three tem- expression dynamics of TGF--induced EMT. poral waves of gene expression, early (epithelial-high), inter- mediate (partial-EMT-high) and late (mesenchymal-high). ETS2, HNF4A and JUNB are putative master TFs of the This observation is consistent with both the time-course partial-EMT state EMT transcriptomic results reported by other groups (31) and the current prevailing regulatory model of EMT, the The gene expression dynamics confirmed that A549 cells Cascading Bistable Switches (CBS) model, which predicts progressively transited through three distinctive states dur- three stable states (21). ing TGF--induced EMT. These coordinated transcrip- Next, we analyzed the Gene Ontology (GO) enrichments tomic changes may be the result of differential TF activi- in each cluster (Figure 1B). The epithelial-high genes were ties. In accordance with the results of previous studies (18), highly enriched with GO terms related to cell cycle. Re- we found that the expression of canonical EMT TFs, such duced proliferation is a well-known phenotype associated as TWIST1/2 and ZEB1/2, was not significantly altered with TGF- signaling (6), and the observed enrichment in in A549 cells undergoing TGF--induced EMT (Supple- cell cycle functions for the early-wave genes was consistent mentary Figure S3A). Although we observed upregulated with reports demonstrating a reduced proliferation rate in SNAI1/2 expression, as reported previously (18), the max- cells undergoing EMT (32). In contrast, the top GO terms imal TPM (transcripts per million) values for SNAI1 (5.58 enriched in the partial-EMT-high gene cluster were regu- TPM) and SNAI2 (4.11 TPM), which were confirmed by lation of localization and locomotion, and cell adhesion absolute qPCR (Supplementary Figure S3B), were within was the top GO term associated with the mesenchymal- the bottom 40% of the values for expressed genes accord- high genes. These results are consistent with the increased ing to the definition of Hebenstreit et al. (33), which sug- cell motility observed upon TGF- treatment, another hall- gested that these genes are unlikely to be master TFs for mark of EMT. Despite these differences, all three clusters the EMT of A549 cells. Thus additional unknown factors were enriched with GO terms related to ‘response to stim- may be involved. We therefore calculated the binding mo- tif enrichment of selected TFs, with the following require- 0 h 6 h 12 h 24 h 36 h 48 h 72 h 96 h Partial-EMT Mesenchymal high genes Epithelial high genes high genes Nucleic Acids Research, 2016, Vol. 44, No. 6 2519 ments for each factor: (i) has a known position weight ma- EMT regulators, we assembled a dynamic GRN for EMT trix; (ii) with a maximal TPM of >10; (iii) shows differen- using the DREM (Dynamic Regulatory Events Miner) soft- tial expression during TGF--induced EMT of A549 cells; ware (27). For this analysis, we integrated the time-course (iv) is found within enhancers associated with epithelial- gene expression data with the computationally predicted high, partial-EMT-high and mesenchymal-high genes. We genome-wide binding data for 52 TFs that were differen- applied the TPM cutoff to focus on the TFs with relatively tially expressed during TGF--induced EMT. DREM iden- high expression levels, an intrinsic characteristic of master tified a dynamic GRN with prominent splitting points at TFs that provides advantages in competing for co-factors 6 and 48 h into EMT, which agreed well with the timing to dominate the transcriptional landscape (12). To identify of the putative state transitions from E to P and P to M, potential key regulators responsible for the state transitions, respectively (Figure 3A). Reassuringly, DREM identified we calculated the fold changes in the expression of each fac- ETS2, HNF4A and JUNB as the top-ranked key regula- tor during state transition, using their 0-h level of mRNA tors at the E to P transition. Similar to the E to P split- expression as the reference. We then nominated those TFs ting points, ETS2, HNF4A and JUNB were nominated by with significant binding motif enrichment ( P < 0.05) and DREM as the top-ranked regulators at the P to M transi- fold expression change (absolute fold change ≥4) as the po- tion. In addition to ETS2, HNF4a and JUNB, DREM also tential key regulators (Figure 2A, Supplementary Figure identified known EMT regulators such as SMAD3, ELF3 S4). This analysis confirmed the role of known key regula- and HMGA2 as key TFs operating at the splitting points. tor of EMT, such as SMAD3 and ELF3, in maintaining the Importantly, the expression levels of ETS2, HNF4A and epithelial phenotype. Notably, four TFs, ETS2, HNF4A, JUNB were upregulated at the splitting points compared JUNB and FOXP1, emerged as the top putative regulators to the levels at the 0-h time point, suggesting that they may of the partial-EMT state in proximal enhancers, distal en- play more active roles during state transitions than regu- hancers and enhancers within gene body (Figure 2A, Sup- lators with downregulated expression, such as SMAD3 or plementary Figure S4). Interestingly, whereas the expres- ELF3. Overall, the predictions obtained using the two in- sion level of FOXP1 increased further upon the transition dependent computational approaches (DREM and binding into the mesenchymal state, the expression levels of ETS2, motif enrichment) agreed well, further establishing ETS2, HNF4A and JUNB fell below the 4-fold cutoff at 72 h into HNF4A and JUNB as candidate master regulators of the EMT (Figure 2B). We validated the protein levels of the partial-EMT state. four TFs during EMT (Figure 2C), confirming that ETS2, To gain further insight into the dynamic GRN associated HNF4A and JUNB were transiently overexpressed 12–36 h with EMT, we extracted the core module from the DREM- into EMT and that FOXP1 displayed monotonic upregu- predicted model, retrieving only the interactions connecting lated expression. Consistent with the CBS model of EMT, the putative key regulators and the key markers of EMT the E-low, P-high and M-low expression pattern suggested (CDH1 and CDH2) (Figure 3B). This analysis revealed that ETS2, HNF4A and JUNB are likely regulators of the a highly modular topology, with 22 of 30 possible intra- −16 partial-EMT state and that FOXP1 is a candidate driver of module interactions implemented (modularity P < 2.2 , the mesenchymal state. Fisher’s exact test). This module was characterized by feed- To further validate the predicted binding motif enrich- forward loops and negative feedback loops, key network ment data, we queried the public chromatin immunopre- motifs to generate bimodality, and was consistent with the cipitation followed by sequencing (ChIP-seq) data for the canonical EMT regulatory model (21). Specifically, four of known (SMAD3) and putative (ETS2, HNF4A, JUNB and four putative EMT regulators had an autoregulatory loop, FOXP1) EMT regulators identified in the motif enrich- a common feature of master TFs (12). Finally, the three pu- ment analysis. Because those ChIP-seq analyses were per- tative partial-EMT state regulators used a clique motif, the formed in cells other than A549, we used ENCODE A549 most abundant three-protein interaction motif and one of H3K27ac ChIP-seq data as a reference and retained only the most conserved GRN motifs (34,35), suggesting coop- the TF ChIP-seq peaks that fell within the established A549 erativeness among those regulators. cell active enhancers. We then analyzed the ChIP-seq peaks associated with the 1633 genes that were differentially ex- Validation of the partial-EMT transcriptional regulatory pressed during EMT. Reassuringly, each of the vfi e TFs oc- module cupied the enhancers of a substantial number of the 1633 genes, ranging from 512 genes for JUNB to 784 genes for According to our dynamic EMT GRN model, the three HNF4A (Figure 2D). Furthermore, their target EMT genes putative partial-EMT state regulators, ETS2, HNF4A that were inferred from ChIP-seq data showed a significant and JUNB, autoregulate and positively regulate each −16 overlap (pair-wise overlap P < 10 , Fisher’s exact test), other, repress the expression of the epithelial state mark- which was consistent with the binding motif prediction re- ers and upregulate the expression of mesenchymal state sults (Figure 2D). regulators/markers. To validate the core partial-EMT reg- ulatory module, we first mined public ChIP-seq data for evidence of direct regulation. Binding events of the puta- ETS2, HNF4a and JUNB control the transition through the tive master TFs on themselves were readily supported by partial-EMT state in a dynamic GRN model the ChIP-seq data (Figure 3C–E). To validate the ChIP- To examine whether the same set of master regulators could seq results, we performed ChIP-qPCR on selected sites and be derived using an independent algorithm and to obtain a confirmed that ETS2, HNF4A and JUNB directly colocal- dynamic, systems-level overview of the roles of the putative ized to their own enhancers in A549 cells (Supplementary 2520 Nucleic Acids Research, 2016, Vol. 44, No. 6 Figure 2. Integrative analysis identified ETS2, HNF4A, JUNB and FOXP1 as the putative master TFs of EMT. ( A) Enrichment P -values of TF binding motifs in the proximal enhancers of partial-EMT high genes (x-axis) were plotted against the fold changes in gene expression during EMT (y-axis; 24 h versus 0 h, left panel; 72 h versus 0 h, right panel). Each circle represents a TF. The dotted lines represent the cutoff for the enrichment (P -value < 0.05) and the expression change (absolute fold change ≥ 4). (B) Levels of ETS2, HNF4A, JUNB and FOXP1 mRNA during EMT as estimated using RNA-seq. (C) Immunoblots showing the levels of ETS2, HNF4A, JUNB and FOXP1 proteins during EMT. (D) Venn diagram showing the overlap of the targets of the vfi e EMT TFs derived from published ChIP-seq data. Figures S5A and S6). Moreover, gene targets supported by function during the partial-EMT state, siRNAs were in- the ChIP-seq data and gene targets supported by the motif troduced into A549 cells 24 h prior to TGF- treatment, enrichment data showed extensive overlap (Supplementary and their effects on mRNA and protein expression were Figure S5B), suggesting that the majority of the predicted determined 24 h into TGF--induced EMT. Compared regulations were the result of direct TF-DNA interactions. to the effects of the control siRNAs, siRNAs for ETS2, To validate the predicted functional gene regulation by HNF4A or JUNB significantly reduced the levels of ETS2, the three putative partial-EMT master regulators, we per- HNF4, JUNB, FOXP1 and CDH2 mRNAs and proteins formed siRNA-mediated silencing of each regulator in and markedly increased the levels of CDH1 mRNA and A549 cells undergoing TGF--induced EMT. To focus on protein (Figure 4A–D), which was consistent with the pre- Nucleic Acids Research, 2016, Vol. 44, No. 6 2521 Figure 3. ETS2, HNF4A and JUNB are key regulators in a dynamic EMT Gene Regulatory Network model. (A) DREM output showing a dynamic gene regulatory map of EMT. The y-axis shows the levels of gene expression normalized to those at 0 h. Each line represents a path (cluster) of genes with a similar expression pattern, and nodes represent the hidden states of a hidden Markov model. The green nodes represent splitting points. For each splitting point, TFs with a split score < 0.001 were listed in ranked order of scores and were colored according to the expression level changes (blue for upregulation; red for downregulation). (B) Diagram illustrating the putative core regulatory network derived from DREM analysis. (C) Genome browser representation of ETS2, HNF4A, JUNB and FOXP1 binding events near ETS2 gene. (D) Same as (C) for HNF4A gene. (E) Same as (C) for JUNB gene. dictions of the putative EMT GRN model. Importantly, A unique feature of our predicted partial-EMT regula- knocking down the expression of any one of the ETS2, tory circuit was that three master regulators were identified. HNF4A and JUNB trio had significant effects, suggesting These master regulators displayed significant co-expression that their cooperative activity is critical for EMT. More- and formed a fully connected module connected by positive over, treating A549 cells undergoing EMT with siRNAs for feedback loops and autoregulations (a clique). Those fea- ETS2, HNF4A or JUNB substantially reduced their abil- tures are indicative of cooperative activity. Because a hall- ity to migrate and invade (Figure 4E, F, P < 0.001) with- mark of enhancer activation is the cooperative assembly out a detectable influence on cell’s proliferation or viabil- of a large trans-factor complex consisting of multiple TFs ity (data not shown), supporting the critical roles of these and their co-factors (9), we proceeded to test the hypothesis TFs in modulating key phenotypes of EMT. We next tested that ETS2, HNF4A and JUNB are part of the same trans- whether overexpression of the three master TFs was as- activation complex. Consistent with our ChIP-qPCR data, sociated with a poor clinical outcome by examining three which showed colocalization of ETS2, HNF4A and JUNB independent lung adenocarcinoma data sets. We classified on numerous enhancers (Supplementary Figures S5A and the samples into two categories, those that overexpressed at S6), examining ChIP-seq data revealed that the binding sites least one of the three master TFs and those in which none of occupied by ETS2 were also highly occupied by HNF4A the master TFs was overexpressed. Kaplan-Meier survival and JUNB and that the genomic locations of their maxi- analysis revealed that patients with elevated expression of at mal peak intensities were within 50 bp of each other (Figure least one of the three genes had a significantly worse chance 5A, B). We next evaluated the physical interaction of ETS2, of relapse-free survival (Figure 4G), demonstrating that the HNF4A and JUNB using immunoprecipitation. Compared overexpression of ETS2, HNF4a or JUNB is a robust prog- to the control antibodies, the antibodies specific for each nostic indicator. of the three TFs readily pulled down the other two fac- 2522 Nucleic Acids Research, 2016, Vol. 44, No. 6 Figure 4. ETS2, HNF4A and JUNB play critical roles in partial-EMT. (A) Quantitative reverse transcription polymerase chain reaction results showing the levels of gene expression in A549 cells during TGF--induced EMT after silencing ETS2 expression using a specific siRNA. n = 3; error bars indicate mean ± SD. * P < 0.05; ** P < 0.01, determined using the two-tailed Student’s t-test. (B) Same as (C) for silencing HNF4A. (C) Same as (A) for silencing JUNB. (D) Immunoblotting analysis of the protein abundance of indicated genes in A549 cells undergoing EMT treated with siRNAs targeting ETS2, HNF4A or JUNB. (E) A549 cells undergoing TGF--induced EMT were treated with siRNAs targeting ETS2, HNF4A or JUNB and were subjected to a migration assay. The migratory cells were quantified (bar charts). Scale bars: 100 m. n = 6; error bars indicate mean ± SD. * P < 0.05; ** P < 0.01, determined using the two-tailed Student’s t-test. (F) Same as (E) for the invasion assay. (G) Kaplan–Meier survival analysis based on the ETS2, HNF4A and JUNB expression levels in three independent Lung Adenocarcinoma data sets for disease-free survival. tors (Figure 5C). Moreover, immunofluorescence staining To determine whether super-enhancers also characterize showed that ETS2, HNF4A and JUNB were colocalized in the identified EMT master TFs, we analyzed ENCODE the nucleus of A549 cells 24 h into TGF--induced EMT A549 H3K27ac ChIP-seq data and identified 1050 super- (Figure 5D). Taken together, these data suggested that the enhancers (Figure 6A). Whereas the enhancers associated three factors indeed are part of the same regulatory complex with ETS2 barely missed the cutoff value to be classified as and cooperatively regulate EMT-associated genes. a super-enhancer, HNF4A and JUNB were found to be as- sociated with putative super-enhancers (Figure 6A and B). Super-enhancer characterizes master EMT transcriptional To test whether other key EMT genes were also asso- regulators ciated with super-enhancers, we compiled a list of 120 well-established EMT genes, classified them into epithelial- Recent studies suggested that master TFs dominate the associated markers (E-markers) and mesenchymal- transcriptional programs in various cell types and super- associated markers (M-markers) and evaluated their enhancers drive the expression of master TFs (10–12,36). Nucleic Acids Research, 2016, Vol. 44, No. 6 2523 (Figure 6E, Supplementary Figure S8). The lack of sig- nificant association between super-enhancers in mesenchy- mal cells with M-markers together with the lack of signifi- cant association between super-enhancers in epithelial cells with E-markers raised the interesting possibility that EMT marker genes, regardless of their association with epithelial or mesenchymal phenotype, are always marked with active enhancers. To test this hypothesis, we examined the associ- ation of all 120 EMT marker genes with active enhancers in the 9 tested cell lines and found that the vast majority of both E-marker and M-marker genes are consistently associ- ated with active enhancers (Figure 6F). This observation is consistent with EMT being a ubiquitous biological process and suggested that the core set of EMT-associated genes are generally active in diverse cell types. BRD4 inhibitor abolishes TGF--induced EMT BRD4 is a transcriptional coactivator from the bromod- omain and extraterminal (BET) subfamily of human bro- modomain proteins. Recent studies suggested that BRD4 is generally associated with active enhancers (10). Impor- tantly, JQ1, a recently developed BRD4 inhibitor, prefer- entially abolishes BRD4 functionality at super-enhancers (10), indicating that JQ1 may abolish TGF--induced EMT by inhibiting the master regulators ETS2, HNF4A and JUNB through the disruption of BRD4 functionality. We first analyzed the colocalization of BRD4 to sites occu- pied by EST2, HNF4A and JUNB by performing ChIP- qPCR assays utilizing the same set of primers used in ETS2, HNF4A and JUNB ChIP-qPCRs. Reassuringly, BRD4 binding is readily detected at tested enhancers near ETS2, Figure 5. ETS2, HNF4A and JUNB are synergistic master regulators HNF4A, JUNB, EGFR and VIM, confirming that BRD4 of EMT. (A) Heatmap showing the colocalization of ETS2, HNF4A and colocalized to the same sites occupied by ETS2, HNF4A JUNB at EMT-associated enhancers. The presence of ETS2, HNF4A and JUNB ChIP-seq peaks are displayed within a 6-kb window centered on the and JUNB (Supplementary Figure S9). Interestingly, while ETS2-bound site. (B) Summary plot for the ETS2, HNF4A and JUNB knocking down HNF4A or ETS2 with specific siRNAs sig- ChIP-seq peak enrichment across the ETS2-binding site associated with nificantly reduced BRD4 binding at some of the tested sites, EMT genes. (C) Endogenous association of ETS2, HNF4A and JUNB. At knocking down JUNB consistently reduced BRD4 binding 24 h into TGF--induced EMT, A549 cell lysates were prepared and were at all the tested sites (Supplementary Figure S9). The result immunoprecipitated using the indicated antibody. The presence of associ- ated proteins was then analyzed using immunoblotting. (D) Immunoflu- is consistent with a recent proteome-scale protein-protein orescence staining for ETS2, HNF4A or JUNB in A549 cells stimulated interaction analysis showed that JUNB makes physical con- with TGF- for 24 h. tact with BRD4 (37), suggesting that the master regulator JUNB recruits BRD4 to those sites. To functionally validate the roles of BRD4 inhibitors in intersection with the A549 cell super-enhancer-associated EMT, we treated A549 cells undergoing TGF--induced genes. Interestingly, 28 genes from the list of 120 EMT- EMT with JQ1, which has been extensively utilized to dis- associated genes were associated with the super-enhancers rupt super-enhancers (10,36). We tested the effect of JQ1 of A549 cells (Figure 6A, B, Supplementary Figures S7A– on the expression of key EMT genes by adding JQ1 to −14 C). This enrichment was highly significant ( P < 10 , the cell culture medium at 6 h into TGF--induced EMT Fisher’s exact test), suggesting that super-enhancers may and examined the cells 18 h later. This setup allowed us indeed drive the high expression of numerous key EMT to specifically target ETS2, HNF4A and JUNB associated genes. Importantly, binding motifs for ETS2, HNF4A and super-enhancers at the partial-EMT state without disrupt JUNB were highly enriched in those super-enhancers (P SMAD3 functionality, which is critical for mediating TGF- −6 < 10 ), suggesting that those master regulators drive the  signaling. Indeed, the addition of JQ1 did not interrupt formation of super-enhancers associated with EMT genes SMAD3 nuclear localization following TGF- treatment (Figure 6C and D). (Supplementary Figure S10). The short JQ1 treatment also We next extended the super-enhancer analyses to 8 addi- excluded molecular and functional impacts from JQ1 in- tional cell lines. While all examined cell lines demonstrated duced apoptosis, which occurs after 48 h of JQ1 expo- an enriched association of super-enhancers with EMT sure (38). As expected, the mRNA and protein levels of markers, neither mesenchymal nor epithelial cell lines dis- ETS2, HNF4A and JUNB were significantly reduced com- played biased enrichment toward E-markers or M-markers pared with those of the DMSO controls (Figure 7A, B), 2524 Nucleic Acids Research, 2016, Vol. 44, No. 6 Figure 6. Super-enhancers are associated with master EMT regulators. (A) The distribution of A549 H3K27ac tag intensities revealed 1050 super- enhancers. The red circles indicate EMT genes associated with super-enhancers. (B) Genome browser tracks showing the super-enhancers associated with ETS2, HNF4A, JUNB and SMAD3. The black bars denote super-enhancers. (C) Genome browser tracks showing super-enhancers associated with the EGFR gene and the binding events of ETS2, HNF4A and JUNB around EGFR gene. (D) Table depicting the enrichment of ETS2, HNF4a and JUNB binding motifs in EMT-associated super-enhancers relative to the genomic background. (E) Bar plots showing the enriched association of EMT genes with super-enhancers. Bars represent the observed number of super-enhancer-associated EMT genes and circles represent the expected number of super- enhancer-associated genes assuming no enrichment. * P < 0.001, determined using Fisher’s Exact Test. (F) Bar plots showing the universal association of EMT genes with active enhancers. Bars represent the observed percentages of active-enhancer-associated EMT genes. Nucleic Acids Research, 2016, Vol. 44, No. 6 2525 Figure 7. BRD4 inhibitor abolishes TGF--induced EMT. (A) Quantitative reverse transcription polymerase chain reaction results showing the expression levels of indicated genes in A549 cells undergoing EMT that were treated with various concentrations of JQ1. n = 3; error bars indicate mean ± SD. * P < 0.05; ** P < 0.01, as determined using the two-tailed Student’s t-test. DMSO, dimethylsulphoxide. (B) Immunoblotting analysis of the expression levels of indicated genes in A549 cells undergoing EMT that were treated with various concentrations of JQ1. (C) A549 cells undergoing TGF--induced EMT were treated with 500 nM JQ1 (left panels) or the DMSO control (right panel) and were subjected to a migration assay (top panels) or an invasion assay (bottom panels). (D) The migratory or invasive cells were quantified (bar charts). Scale bars: 100 m. n = 6; error bars indicate mean ± SD. * P < 0.05; ** P < 0.01, determined using the two-tailed Student’s t-test. consistent with the association of these genes with super- of published ChIP-seq data. We successfully recovered and enhancers. Consequently, the expression of FOXP1 and the validated novel partial-EMT regulators using the small cell mesenchymal marker gene (CDH2) regulated by those mas- lung cancer A549 cell line. In principle, our approach can ter TFs was also downregulated. Furthermore, cells under- be extended to any system in which time-course gene expres- going EMT that were treated with JQ1 displayed signifi- sion data is accessible and could be valuable for deciphering cantly reduced abilities to migrate and invade (Figure 7C, EMT mechanisms in other cancers. D, P < 0.01) and thus were losing key EMT characteristics. In sharp contrast to the canonical EMT model, in which a single TF such as SNAI1 of ZEB1 controls the EMT pro- gram, we identified multiple master TFs that synergistically DISCUSSION control the EMT transcriptional program in A549 cells. The master TFs physically interact and demonstrate cumulative EMT has been implicated as an integral part of metasta- binding at enhancers of EMT-associated genes, suggesting sis, and a clear understanding of the underlying molecular that the synergistic action of multiple master TFs could rep- mechanism of the latter is crucial for successful clinical ther- resent an intrinsic property of the EMT transcriptional pro- apies. Unfortunately, EMT is a phenomenon that is highly gram. Similar to the master TFs identified in other cellu- dependent on the cellular context in which it occurs. Despite lar models, the master TFs identified in our study are as- extensive efforts, the regulatory mechanism of EMT in most sociated with super-enhancers. In addition to A549 cells, cell/tissue types is not fully resolved. Our understanding is JQ1 inhibited migration and invasion of multiple breast and particularly poor in the case of partial-EMT, the interme- prostate cancer cell lines (data not shown). Moreover, the diate state that was recently proposed to hold key insights disruption of BRD4 function by JQ1 effectively inhibits tu- into metastasis (16). In this study, we employed a system- mor formation and metastasis in numerous mouse models atic approach to identify all potential key regulators and (39,40). Although it is not possible to deconvolute the role corresponding GRNs of EMT using an integrative analysis of BRD4 inhibitors in EMT using vivo assays owing to the of time-course gene expression data and a large repertoire 2526 Nucleic Acids Research, 2016, Vol. 44, No. 6 massive apoptosis effects from JQ1 administration (36,38), 2. Lamouille,S., Xu,J. and Derynck,R. (2014) Molecular mechanisms of epithelial-mesenchymal transition. Nat. Rev. Mol. Cell Biol., 15, these observations supported that super-enhancers, the crit- 178–196. ical role of which in other aspects of tumorigenesis has been 3. Thiery,J.P., Acloque,H., Huang,R.Y. and Nieto,M.A. (2009) demonstrated in several recent studies, are key players in Epithelial-mesenchymal transitions in development and disease. Cell, EMT and may together represent a general target for coin- 139, 871–890. 4. Yang,J. and Weinberg,R.A. (2008) Epithelial-mesenchymal hibiting multiple hallmarks of cancers. transition: at the crossroads of development and tumor metastasis. Consistent with the mechanism of enhancer formation, Dev. Cell, 14, 818–829. we identified 3 synergistic master TFs regulating TGF- - 5. Garber,K. (2008) Epithelial-to-mesenchymal transition is important induced EMT in A549 cells. Interestingly, each of the three to metastasis, but questions remain. J. Natl. Cancer Inst., 100, 232–233. factors, ETS2, HNF4A and JUNB, has been previously as- 6. Ikushima,H. and Miyazono,K. (2010) TGFbeta signalling: a complex sociated with EMT, albeit in different cellular backgrounds web in cancer progression. Nat. Rev. Cancer, 10, 415–424. (41–43). Aberrant reactivation of lineage specific TFs is a 7. Massague,J. and Gomis,R.R. (2006) The logic of TGFbeta signaling. well-known phenomena during tumorigenesis (44,45). Cur- FEBS Lett., 580, 2811–2820. rently, more than 30 TFs have been implicated in EMT pro- 8. Mullen,A.C., Orlando,D.A., Newman,J.J., Loven,J., Kumar,R.M., Bilodeau,S., Reddy,J., Guenther,M.G., DeKoter,R.P. and grams underlying a variety of tissue remodeling events (1,2). Young,R.A. (2011) Master transcription factors determine It is conceivable that many of them could be reactivated dur- cell-type-specific responses to TGF-beta signaling. Cell, 147, 565–576. ing metastasis. Thus, systematic discovery and characteriza- 9. Spitz,F. and Furlong,E.E. (2012) Transcription factors: from tion of their synergistic interactions may offer substantial enhancer binding to developmental control. Nat. Rev. Genet., 13, insights into the molecular mechanisms underlying EMT 613–626. 10. Loven,J., Hoke,H.A., Lin,C.Y., Lau,A., Orlando,D.A., Vakoc,C.R., and may identify novel prognostic markers as well as new Bradner,J.E., Lee,T.I. and Young,R.A. (2013) Selective inhibition of therapeutic targets for controlling metastasis. tumor oncogenes by disruption of super-enhancers. Cell, 153, 320–334. 11. Whyte,W.A., Orlando,D.A., Hnisz,D., Abraham,B.J., Lin,C.Y., AVAILABILITY OF SUPPORTING DATA Kagey,M.H., Rahl,P.B., Lee,T.I. and Young,R.A. (2013) Master transcription factors and mediator establish super-enhancers at key The data reported in this study have been deposited in GEO cell identity genes. Cell, 153, 307–319. under accession number GSE69667. 12. Hnisz,D., Abraham,B.J., Lee,T.I., Lau,A., Saint-Andre,V., Sigova,A.A., Hoke,H.A. and Young,R.A. (2013) Super-enhancers in the control of cell identity and disease. Cell, 155, 934–947. SUPPLEMENTARY DATA 13. Revenu,C. and Gilmour,D. (2009) EMT 2.0: shaping epithelia through collective migration. Curr. Opin. Genet. Dev., 19, 338–342. Supplementary Data are available at NAR Online. 14. Chao,Y., Wu,Q., Acquafondata,M., Dhir,R. and Wells,A. (2012) Partial mesenchymal to epithelial reverting transition in breast and prostate cancer metastases. Cancer Microenviron., 5, 19–28. ACKNOWLEDGEMENTS 15. Theveneau,E. and Mayor,R. (2012) Neural crest delamination and The authors would like to thank members of the Labora- migration: from epithelium-to-mesenchyme transition to collective cell migration. Dev. Biol., 366, 34–54. tory of Systems Biology for participating in helpful discus- 16. Jordan,N.V., Johnson,G.L. and Abell,A.N. (2011) Tracking the sions and National Center for Protein Science Shanghai for intermediate stages of epithelial-mesenchymal transition in epithelial assistance with confocal microscopy. stem cells and cancer. Cell cycle, 10, 2865–2873. Authors’ contributions: H.C., Y.L., H.L., S.D. and L.Z. 17. Yu,M., Bardia,A., Wittner,B.S., Stott,S.L., Smas,M.E., Ting,D.T., Isakoff,S.J., Ciciliano,J.C., Wells,M.N., Shah,A.M. et al. (2013) performed the experiments. M.X. conducted the com- Circulating breast tumor cells exhibit dynamic changes in epithelial putational analysis. H.C. and H.L. contributed to the and mesenchymal composition. Science, 339, 580–584. computational analysis. P.W. designed and supervised the 18. Liu,J., Hu,G., Chen,D., Gong,A.Y., Soori,G.S., Dobleman,T.J. and study, contributed to the computational analysis and wrote Chen,X.M. (2013) Suppression of SCARA5 by Snail1 is essential for the manuscript. All authors have read and approved the EMT-associated cell migration of A549 cells. Oncogenesis, 2, e73. 19. Zavadil,J. and Bottinger,E.P. (2005) TGF-beta and manuscript. epithelial-to-mesenchymal transitions. Oncogene, 24, 5764–5774. 20. Xu,J., Lamouille,S. and Derynck,R. (2009) TGF-beta-induced epithelial to mesenchymal transition. Cell Res., 19, 156–172. FUNDING 21. Zhang,J.Y., Tian,X.J., Zhang,H., Teng,Y., Li,R.Y., Bai,F., National Natural Science Foundation of China (NSFC) Elankumaran,S. and Xing,J.H. (2014) TGF-beta-induced epithelial-to-mesenchymal transition proceeds through stepwise [31271413]; Science and Technology Commission of Shang- activation of multiple feedback loops. Sci. Signal, 7, ra91. hai Municipality [12DZ1910800]; Chinese Academy of Sci- 22. Denzler,R., Agarwal,V., Stefano,J., Bartel,D.P. and Stoffel,M. (2014) ences ‘Hundred Talent Program’ (to P.W.). Funding for Assessing the ceRNA hypothesis with quantitative measurements of open access charge: National Natural Science Founda- miRNA and target abundance. Mol. Cell, 54, 766–776. tion of China (NSFC) [31271413]; Science and Technology 23. Langmead,B., Trapnell,C., Pop,M. and Salzberg,S.L. (2009) Ultrafast and memory-efficient alignment of short DNA sequences to the Commission of Shanghai Municipality [12DZ1910800]; human genome. Genome Biol., 10, R25. Chinese Academy of Sciences ‘Hundred Talent Program’. 24. Li,B. and Dewey,C.N. (2011) RSEM: accurate transcript Conflict of interest statement. None declared. quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics, 12, 323. 25. Conesa,A., Nueda,M.J., Ferrer,A. and Talon,M. (2006) maSigPro: a REFERENCES method to identify significantly differential expression profiles in time-course microarray experiments. Bioinformatics, 22, 1096–1102. 1. De Craene,B. and Berx,G. (2013) Regulatory networks defining EMT during cancer initiation and progression. Nat. Rev. Cancer, 13, 97–110. Nucleic Acids Research, 2016, Vol. 44, No. 6 2527 26. Frith,M.C., Fu,Y., Yu,L., Chen,J.F., Hansen,U. and Weng,Z. (2004) integrated cellular networks of transcription-regulation and Detection of functional DNA motifs via statistical protein-protein interaction. Proc. Natl. Acad. Sci. U.S.A., 101, over-representation. Nucleic Acids Res., 32, 1372–1381. 5934–5939. 27. Schulz,M.H., Devanny,W.E., Gitter,A., Zhong,S., Ernst,J. and 36. Asangani,I.A., Dommeti,V.L., Wang,X., Malik,R., Cieslik,M., Bar-Joseph,Z. (2012) DREM 2.0: Improved reconstruction of Yang,R., Escara-Wilke,J., Wilder-Romans,K., Dhanireddy,S., dynamic regulatory networks from time-series expression data. BMC Engelke,C. et al. (2014) Therapeutic targeting of BET bromodomain Syst. Biol., 6, 104. proteins in castration-resistant prostate cancer. Nature, 510, 278–282. 28. Liu,T., Ortiz,J.A., Taing,L., Meyer,C.A., Lee,B., Zhang,Y., Shin,H., 37. Hein,M.Y., Hubner,N.C., Poser,I., Cox,J., Nagaraj,N., Toyoda,Y., Wong,S.S., Ma,J., Lei,Y. et al. (2011) Cistrome: an integrative Gak,I.A., Weisswange,I., Mansfeld,J., Buchholz,F. et al. (2015) A platform for transcriptional regulation studies. Genome Biol., 12, R83. human interactome in three quantitative dimensions organized by 29. Cerami,E., Gao,J., Dogrusoz,U., Gross,B.E., Sumer,S.O., stoichiometries and abundances. Cell, 163, 712–723. Aksoy,B.A., Jacobsen,A., Byrne,C.J., Heuer,M.L., Larsson,E. et al. 38. Lockwood,W.W., Zejnullahu,K., Bradner,J.E. and Varmus,H. (2012) (2012) The cBio cancer genomics portal: an open platform for Sensitivity of human lung adenocarcinoma cell lines to targeted exploring multidimensional cancer genomics data. Cancer Discov., 2, inhibition of BET epigenetic signaling proteins. Proc. Natl. Acad. Sci. 401–404. U.S.A., 109, 19408–19413. 30. Gao,J., Aksoy,B.A., Dogrusoz,U., Dresdner,G., Gross,B., 39. Alsarraj,J. and Hunter,K.W. (2012) Bromodomain-Containing Sumer,S.O., Sun,Y., Jacobsen,A., Sinha,R., Larsson,E. et al. (2013) Protein 4: A Dynamic Regulator of Breast Cancer Metastasis Integrative analysis of complex cancer genomics and clinical profiles through Modulation of the Extracellular Matrix. Int. J. Breast using the cBioPortal. Sci. Signal., 6, pl1. Cancer, 2012, 670632. 31. Javaid,S., Zhang,J., Anderssen,E., Black,J.C., Wittner,B.S., 40. Hu,Y., Zhou,J., Ye,F., Xiong,H., Peng,L., Zheng,Z., Xu,F., Cui,M., Tajima,K., Ting,D.T., Smolen,G.A., Zubrowski,M., Desai,R. et al. Wei,C., Wang,X. et al. (2015) BRD4 inhibitor inhibits colorectal (2013) Dynamic chromatin modification sustains cancer growth and metastasis. Int. J. Mol. Sci., 16, 1928–1948. epithelial-mesenchymal transition following inducible expression of 41. Yang,M., Li,S.N., Anjum,K.M., Gui,L.X., Zhu,S.S., Liu,J., Snail-1. Cell Rep., 5, 1679–1689. Chen,J.K., Liu,Q.F., Ye,G.D., Wang,W.J. et al. (2013) A 32. Evdokimova,V., Tognon,C., Ng,T., Ruzanov,P., Melnyk,N., Fink,D., double-negative feedback loop between Wnt-beta-catenin signaling Sorokin,A., Ovchinnikov,L.P., Davicioni,E., Triche,T.J. et al. (2009) and HNF4alpha regulates epithelial-mesenchymal transition in Translational activation of snail1 and other developmentally hepatocellular carcinoma. J. Cell Sci., 126, 5692–5703. regulated transcription factors by YB-1 promotes an 42. Gervasi,M., Bianchi-Smiraglia,A., Cummings,M., Zheng,Q., epithelial-mesenchymal transition. Cancer Cell, 15, 402–415. Wang,D., Liu,S. and Bakin,A.V. (2012) JunB contributes to Id2 33. Hebenstreit,D., Fang,M., Gu,M., Charoensawan,V., van repression and the epithelial-mesenchymal transition in response to Oudenaarden,A. and Teichmann,S.A. (2011) RNA sequencing transforming growth factor-beta. J. Cell Biol., 196, 589–603. reveals two major classes of gene expression levels in metazoan cells. 43. Polydorou,C. and Georgiades,P. (2013) Ets2-dependent trophoblast Mol. Syst. Biol., 7, 497. signalling is required for gastrulation progression after primitive 34. Stergachis,A.B., Neph,S., Sandstrom,R., Haugen,E., Reynolds,A.P., streak initiation. Nat. Commun., 4, 1658. 44. Darnell,J.E. Jr (2002) Transcription factors as targets for cancer Zhang,M., Byron,R., Canfield,T., Stelhing-Sun,S., Lee,K. et al. therapy. Nat. Rev. Cancer, 2, 740–749. (2014) Conservation of trans-acting circuitry during mammalian 45. Vaquerizas,J.M., Kummerfeld,S.K., Teichmann,S.A. and regulatory evolution. Nature, 515, 365–370. 35. Yeger-Lotem,E., Sattath,S., Kashtan,N., Itzkovitz,S., Milo,R., Luscombe,N.M. (2009) A census of human transcription factors: Pinter,R.Y., Alon,U. and Margalit,H. (2004) Network motifs in function, expression and evolution. Nat. Rev. Genet., 10, 252–263.

Journal

Nucleic Acids Research – Oxford University Press

Published: Apr 7, 2016

There are no references for this article.