The NFκB subunit RELA is a master transcriptional regulator of the committed epithelial-mesenchymal transition in airway epithelial cells

The epithelial-mesenchymal transition (EMT) is a multistep dedifferentiation program important in tissue repair. Here, we examined the role of the transcriptional regulator NF-κB in EMT of primary human small airway epithelial cells (hSAECs). Surprisingly, transforming growth factor β (TGFβ) activated NF-κB/RELA proto-oncogene, NF-κB subunit (RELA) translocation within 1 day of stimulation, yet induction of its downstream gene regulatory network occurred only after 3 days. A time course of TGFβ-induced EMT transition was analyzed by RNA-Seq in the absence or presence of inducible shRNA-mediated silencing of RELA. In WT cells, TGFβ stimulation significantly affected the expression of 2,441 genes. Gene set enrichment analysis identified WNT, cadherin, and NF-κB signaling as the most prominent TGFβ-inducible pathways. By comparison, RELA controlled expression of 3,138 overlapping genes mapping to WNT, cadherin, and chemokine signaling pathways. Conducting upstream regulator analysis, we found that RELA controls six clusters of upstream transcription factors, many of which overlapped with a transcription factor topology map of EMT developed earlier. RELA triggered expression of three key EMT pathways: 1) the WNT/β-catenin morphogen pathway, 2) the JUN transcription factor, and 3) the Snail family transcriptional repressor 1 (SNAI1). RELA binding to target genes was confirmed by ChIP. Experiments independently validating WNT dependence on RELA were performed by silencing RELA via genome editing and indicated that TGFβ-induced WNT5B expression and downstream activation of the WNT target AXIN2 are RELA-dependent. We conclude that RELA is a master transcriptional regulator of EMT upstream of WNT morphogen, JUN, SNAI1–ZEB1, and interleukin-6 autocrine loops.

The respiratory tract represents one of the largest surfaces of the human body exposed to the external environment (1,2). Under normal conditions, this mucosal surface is a semi-impermeable barrier formed by highly differentiated epithelial cells that allow fluid conservation, permit gas exchange, and provide physical protection from pathogen invasion. Disruption of the epithelial barrier releases epithelial growth factors, including transforming growth factor (TGF), 3 sequestered within the extracellular matrix that coordinate expression of gene regulatory networks (GRNs) to promote repair. Phenotypically, TGF orchestrates a dedifferentiation program known as the type II epithelial-mesenchymal transition (EMT). Type II EMT involves extensive cytosolic restructuring, resulting in the loss of apical-basal polarity, dissolution of adherens junctions, enhanced motility, and expression of fibrotic genes (for an indepth review, see Ref. 3). As a consequence of this process, epithelial cells acquire stem cell-like properties, permitting the transitioned mesenchymal cell to repopulate regions of denuded epithelium, promoting tissue repair and extracellular matrix remodeling (3).
Transition to the mesenchymal state is not binary but rather the product of sequential cell-state changes beginning from the differentiated epithelial state transitioning into uncommitted "partial EMT (pEMT)" state(s) (4). Transition into pEMT is initiated by TGF␤ binding to TGF␤ receptor type II, signaling intracellularly via the SMAD-dependent, "canonical" and This work was supported by NIAID, National Institutes of Health Grant AI062885; National Center for Advancing Translational Sciences, National Institutes of Health Grant ULTR00002373; and National Science Foundation Grant DMS-1361411. The authors declare that they have no conflicts of interest with the contents of this article. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health. This article contains Figs. S1 and S2 and Tables S1-S3. The data discussed in this publication have been deposited in NCBI's Gene Expression Omnibus and are accessible through GEO Series accession number GSE84135. 1 Both authors contributed equally to this work. 2  SMAD-independent, "noncanonical" pathways (5,6). These two arms function coordinately to activate a cascade of transcription factors that result in ordered patterns of gene expression and inhibition, resulting in a metastable pEMT state (7). Gene expression patterns of SMAD are controlled in a cell type-specific manner by interaction with master transcription factors that function in coordinately regulated complexes to form acetylated histone-rich superenhancers on mesenchymal genes, promoting chromatin remodeling (8,9). Persistent growth factor and microenvironmental stimulation is required to produce a stable mesenchymal transition. Time-resolved RNA profiling coupled with ChIP-sequencing studies have elucidated the role of the JUNB complex as an example of a "master" transcription factor of the mesenchymal state (9). Another master regulator is zinc finger E-box-binding homeobox 1 (ZEB1), a bifunctional transcriptional regulator activated by SNAI1 that down-regulates E-cadherin and the EMT-suppressing microRNA-200 family in a complex with the C-terminal binding protein repressor (10 -12). Although the posttranslational regulation of the SNAI1-ZEB1 autoamplification loop by microRNAs is well-described (10), the factor(s) controlling the initial transcriptional spike in SNAI1 expression required to trigger the amplification loop is not known in type II EMT.
In addition to master transcription factors, autocrine and paracrine growth factor signaling plays an important role in the maintenance of the mesenchymal state (9). For example, up-regulation of the WNT morphogen ligands that activate intracellular ␤-catenin (CTNNB)-transcription factor (TCF3) complex in neighboring cells to maintain SNAI1 and TGF␤ expression is an important component of EMT (13,14). Local up-regulation of IL6 promotes growth factor independence and expression of stem cell-like features in EMT via autocrine/ paracrine signaling (15,16). The identity of autocrine pathways, their timing of expression, and their relationship to the master transcription factors in type II EMT are not known.
Our previous work identified that persistent NF-B action is central to the fully committed mesenchymal regulatory network in type II EMT (17,18). However, the role of RELA in the initiation of type II EMT has not been determined. To address this gap, we conducted time course experiments during the initiation of TGF␤-induced mesenchymal transition in primary human small airway epithelial cells (hSAECs). Although nuclear translocation of RELA was detected within 1 day of TGF␤ stimulation, the activating Ser-276 phosphorylation only significantly occurred after 3 days of stimulation, an event associated with activation of the NF-Bdependent GRN. hSAECs expressing a doxycycline (Dox)-regulated shRNA were profiled by RNA-Seq in the absence or presence of RELA knockdown (KD). We observed that RELA controls the autocrine growth factors WNT and IL6 as well as the JUN and SNAI1 master transcription factors. These data implicate RELA as a master transcription regulator of type II EMT.

Results
We previously characterized a reproducible model of type II EMT using a continuously replicating line of hSAECs immortalized with human telomerase (6,17). Unbiased genomic and proteomic profiling studies indicate that hSAEC expression patterns closely mimic those of primary small airway epithelial cells (19,20). In the absence of stimulation, hSAECs exhibit cuboidal epithelial morphology, express cytokeratin 17, and form pseudostratified columnar epithelium in air-liquid interfaces (21).
Upon TGF␤ stimulation, hSAECs exhibit a time-dependent mesenchymal dedifferentiation program associated with formation of filamentous (F) actin stress fibers, detected by phalloidin staining. Significantly enhanced phalloidin staining was seen initially after 12 h and dramatically increased by 48 h (Fig.  1A).

Kinetics of NF-B/RELA activation
Our earlier RNA-Seq analysis of hSAECs in a fully committed EMT state after 15 days of TGF␤ stimulation indicated the presence of activated RELA and a robust NF-B GRN signature (17). To understand the early kinetics of NF-B activation, we examined the effects of TGF␤ on RELA nuclear translocation by confocal fluorescence microscopy (CFM) over a detailed time course from 0 to 72 h. In unstimulated control hSAECs, RELA was predominately cytoplasmic (Fig. 1B). A small fraction of nuclei were RELA-positive by 12 h, and RELA was predominately nuclear after 48 h (Fig. 1, B and D).
Earlier work has shown that activation of the RELA transcriptional activity requires phosphorylation of Ser-276 by ribosomal S6 kinases (22,23). Ser-276 phosphorylation is a rate-limiting step required for p300-dependent acetylation and complex formation with bromodomain-containing protein 4 (BRD4) (18,24,25), a coactivator required for stable EMT (9,18). In contrast to the kinetics of nuclear RELA translocation, the formation of activating phospho-Ser-276 RELA is initially detected after 2 days and dramatically induced by 3 days of TGF␤ stimulation where a 15-and 42-fold increase is seen relative to unstimulated hSAECs, respectively ( Fig. 1, C and E).
The changes in RELA nuclear translocation were validated by Western immunoblotting of nuclear extracts where TGF␤ produced a time-dependent increase in nuclear RELA abundance after 2 and 3 days of TGF␤ stimulation (Fig. 1, F and G). Collectively, these data indicate that TGF␤ induces RELA translocation followed thereafter by RELA activating phosphorylation.

Validation of Dox-inducible RELA silencing
To systematically identify the direct genomic targets of RELA in the initiation of mesenchymal transition, we engineered hSAECs to stably express a Dox-inducible RELA shRNA. The rationale for selecting a Dox-inducible system was that this approach enables study of the effect of acute RELA depletion without the well-described compensation for NF-B deficiency in knockout cell lines (26). We first confirmed the effect of Dox treatment on RELA abundance where 5 days of Dox treatment reduced endogenous RELA mRNA to 15% that of untreated hSAECs ( Fig. 2A) and depleted total cellular protein to similar levels (Fig. 2B). To test whether this level of depletion functionally impaired RELA signaling, we examined the effect of TNF on the NF-Bdependent GRN. For this purpose, we examined the response of TNFAIP3/A20, NFKBIA/ NFB is a master regulator of EMT autocrine loops IB␣, and IL6, representative genes directly regulated by NF-B (27). We observed that an acute (1-h) TNF␣ stimulus induced TNFAIP3/A20 mRNA by 44-fold. By contrast, in RELA KD, TNFAIP3/A20 mRNA was only induced by 8-fold, representing a significant inhibition (p Ͻ 0.01, ANOVA, post hoc Tukey's t test). Similar findings were observed for the NF-B/RELAdependent NFKBIA/IB␣ and IL6 genes where the 22-fold induction of NFKBIA/IB␣ in WT cells was reduced to 5-fold after RELA silencing. The 80-fold induction of IL6 in WT cells was reduced to 5-fold in the RELA KD cells (p Ͻ 0.05, ANOVA, post hoc Tukey's test; Fig. 2A). These data indicate that the RELA KD functionally interfered with the NF-B-dependent GRN and would be suitable for testing the role of NF-B in the initiation of TGF␤-induced EMT.

RELA KD blocks the initiation of TGF␤-induced mesenchymal transition
We next examined the response of the RELA KD hSAECs to TGF␤-induced mesenchymal transition. WT or RELA KD cells were subjected to a short time course of TGF␤ induction (0, 1, and 3 days). Redistribution of F-actin was measured using phalloidin staining in CFM. In unstimulated controls, phalloidin staining is largely located in the cellular periphery (Fig. 2D, top). By contrast, TGF␤ stimulation produced organized F-actin formation throughout the cytoplasm and nucleus after 3 days. RELA KD cells maintained a low level of phalloidin staining localized to the cellular periphery (Fig. 2D, top).
To confirm the persistence of RELA silencing over the 3-day time course, RELA abundance was measured by CFM and quantitative real-time RT-PCR (Q-RT-PCR) (primer sequences are given in Table S1). The cytoplasmic-to-nuclear redistribution of RELA and enhanced total RELA abundance were evident in the WT cells but significantly reduced in the RELA KD cells (Fig. 2D, bottom). Similarly, RELA mRNA was reduced to Ͻ15% that of WT cells over the 3 days of TGF␤ stimulation (Fig. 2E). Compared with WT cells, RELA KD cells showed a significant reduction in TGF␤-induced fibronectin (FN1) expression from 900-fold in WT cells to 200-fold in RELA KD cells (Fig. 2E). Similarly, RELA silencing produced a significant inhibition of the 42-fold increase in matrix metalloproteinase 9 (MMP9) and the 5.5-fold increase in VIM expression to less than 2-fold (Fig. 2, E and F). We interpret these data to show that RELA plays an important role in initiating the early stages of type II EMT.

Systematic identification of NF-Bdependent genes early in the mesenchymal transition
To systematically identify the RELA-induced GRN in the initial stages of the mesenchymal transition, WT or RELA KD cells were stimulated in the absence or presence of TGF␤ for 1 and 3 days. Each experimental condition was reproduced in triplicate, representing 18 samples subjected to RNA-Seq analysis (Fig.  S1). The quality metrics and numbers of mapped reads for the data set are given in Table S2.
Principal component analysis was used as an unsupervised approach to visualize sample reproducibility and gene expression similarity of biological treatments. The first principal component separated cells on the time of TGF␤ treatment (Fig. 3B). By contrast, the second principal component separated cells on the basis of RELA activity. We observed that each biological replicate closely clustered with the other samples receiving the same treatment, indicating a high concordance of the biological replicates. Prior to TGF␤ stimulation, RELA KD cells were widely separated from WT cells. After 1 day of treatment, both WT and RELA-silenced cells closely clustered; these populations then separated after 3 days of TGF␤ treatment. We interpret these data to indicate that RELA controls global gene expression in the absence of stimulation, and its major effect is seen after 3 days of stimulation.

Analysis of the TGF␤-induced transcriptional program
We first considered the analysis of the TGF␤-induced GRN in the absence of shRNA expression followed by identification of the network component that is RELA-dependent. TGF␤ treatment changed the expression of 2,441 unique genes after filtering for a 4-fold change in gene expression and an adjusted p value of Ͻ0.01.
To ensure that we have captured the initial steps in mesenchymal transition, we examined the patterns of 28 well-established epithelial and mesenchymal "signature genes" (28) by hierarchical clustering. TGF␤ induced at least two temporally coordinated waves of gene expression, identifying three major clusters: a group high in the untreated cells (Fig. 3A, WTd0), a transiently up-regulated group at 1 day of stimulation (Fig. 3A,  WTd1), and a group later up-regulated after 3 days of stimulation (Fig. 3A, WTd3). The highly expressed genes in the WT untreated cells included characteristic epithelial genes, including keratin (KRT), NOTCH1, and TP53, whose expression fell after TGF␤ stimulation. A smaller group of genes was transiently up-regulated at 1 day with expression falling by 3 days. These included integrin (ITGA)-1, the transcription factor YY1, and the S100 calcium-binding protein (S100A14). Also, a group of genes was up-regulated at 1 day and sustained for 3 days, included the master mesenchymal transcription factors SNAI1, TWIST1, and JUN.
Finally, the genes induced after 3 days of treatment included the characteristic mesenchymal genes FN1, VIM, secreted protein, acidic, cysteine-rich (SPARC), MMP, and others (Fig. 3). Collectively, these data indicate that we have captured components of the epithelial-pEMT-mesenchymal transition, consis-

NFB is a master regulator of EMT autocrine loops
tent with the three-state transition observed in type III EMT in malignant cells (10). We also noted that the expression of the master transcription factors SNAI1, TWIST1, and JUN preceded that of ZEB1, consistent with the sequential expression of the SNAI1-ZEB1 regulatory network (10). We interpret these data to indicate that we have a robust data set that captures the early expression events in the initiation and pEMT-mesenchymal transition of primary cells.
To more systematically understand the temporal expression patterns, we compared the overlap of the differentially expressed genes on day 1 versus that of day 3. 417 genes were uniquely changed after 1 day of treatment, 915 genes were common, and 1,119 genes were regulated 3 days after treatment (Fig. 3B). To obtain information about the molecular processes underlying the TGF␤-induced gene regulatory network, we first performed gene set enrichment analysis (GSEA) on the highly differentially expressed genes in WT cells after 1 and 3 days. GSEA assumes that biological processes occur by coregulation of genes within common pathways and statistically identifies enriched pathways based on prior knowledge (29). The five most highly significant "hallmark" pathways included the EMT itself, "estrogen response," "TNF signaling via the NF-B pathway," and the "p53" and "coagulation" pathways ( Fig. 3C). Interestingly, 54 genes, over 25% of the hallmark EMT reference set, were identified in our RNA-Seq analysis; over half of these EMT genes contributed to the "estrogen response" and "TNF␣ signaling via NF-B" gene sets. A separate functional classification analysis on the WT day 1 gene set was performed, examining for enrichment of signaling pathways, and resulted in the identification of "integrin," "WNT signaling," "angiogenesis," and "chemokine/cytokine signaling" (Fig. 3D). We noted that the classifications by these two different analyses were highly concordant, with the genes contributing to the integrin and WNT signaling pathways also being contained within the GSEA/hallmark EMT gene set, and the chemokine/cytokine signaling within the hallmark GSEA "TNF␣ signaling" gene sets. Importantly, these data confirm our earlier observations that the NF-B-dependent GRN is a core component of the epithelial-to-pEMT state transition (17,30).
After 3 days of TGF␤ stimulation in WT cells, GSEA again identified "EMT," estrogen response, and TNF␣ signaling via NF-B as dominant molecular processes; additional pathways that appeared included the cell cycle-regulatory pathways "G 2 M checkpoint" and "E2F targets" (Fig. 3E). Signaling pathway enrichment was predicted for WNT, "cadherin," angiogenesis, and integrin pathways along with the chemokine/cytokine signaling pathway observed after 1 day of TGF␤ treatment (Fig. 3F).

The NF-B-dependent transcription program controls master EMT regulators
3,138 unique differentially expressed genes were identified as NF-Bdependent by a 4-fold difference in comparing RELA KD cells versus WT for each sampling point in the experiment (Fig. 4A). The largest number of genes, 1,870, was identified in the RELA KD versus WT cells prior to TGF␤ treatment, consistent with the distinct mapping of the control cells in the principal component analysis (cf. Figs. 4A and S1). After TGF␤ stimulation, 671 genes were affected by RELA KD after 1 day of TGF␤ treatment, and 1,055 genes were affected after 3 days of TGF␤ treatment. Many are unique to each time point (Fig. 4A). Importantly, the number of NF-B-dependent genes represented over half of those regulated by TGF␤ at the same time point, indicating to us that NF-B has a global effect on the EMT transcriptional program. The 3,138 unique NF-B-dependent genes were subjected to analysis for molecular function by GO (Fig. 4B); the major functions identified included nucleic acid-binding transcriptional activity (10%) and binding (26%). GSEA and signaling pathway analyses were conducted for each KD versus WT comparison, day 0 ( Fig. 4, C and D), day 1 (Fig. 4, E and F), and day 3 (Fig. 4, G and H). Here the rankordered hallmark gene sets underwent a time-dependent shift from TNF␣ signaling via NF-B being the top ranked gene set in unstimulated cells (Fig. 4C) to EMT becoming the top ranked gene set in cells stimulated with TGF␤ for 3 days (Fig. 4G). Similarly, the rank order of the signaling pathways changed where the most predominant signaling pathways in control (day 0),"gonadotropin-releasing hormone," "epidermal growth factor receptor," and "chemokine and cytokine signaling" (Fig. 4D), were replaced by the WNT and cadherin signaling pathways after 3 days of TGF␤ treatment (Fig. 4H).

Inference of transcription factor networks under NF-B/RELA control
To obtain more insight into the role of RELA in control of transcription factors in type II EMT, we conducted upstream regulator analysis in Ingenuity Pathway Analysis. Upstream regulator analysis uses the observed expression change in our data sets and prior knowledge of expected effects between transcriptional regulators and their target genes to identify the cascades of upstream regulators responsible (31). We predicted 73 upstream transcription regulators that show a significant activation or inhibition effect (activation z-score Ͼ or Ͻ2, p value of overlap Ͻ0.01) that clustered into six expression patterns (Fig. 5). Many master transcription factors of the EMT were identified in Clusters 2, 5, and 6. Of these, we noted that Cluster

NFB is a master regulator of EMT autocrine loops
2 contained the master transcription factor JUN/Fos (9). Similarly, Cluster 5 included TWIST2/SNAI2 and TCF3. Cluster 6 included JUNB, SNAI1, SMADs 2/3/4, CTNNB1, and NF-B itself. From inspection of the effect of RELA KD, the actions of the transcription factors in Clusters 2 and 6 are the most significantly reduced. We also note that RELA KD results in the acti-

NFB is a master regulator of EMT autocrine loops
vation of mesenchymal-to-epithelial regulators SMAD7, MYC, and TBX2 (in Clusters 3 and 4).
Collectively, our data indicate that NF-B signaling functions as a master regulator of type II EMT upstream of the mesenchymal state-promoting transcription factors and negatively regulates epithelial state-promoting factors. We next sought to validate the NF-B dependence of the WNT pathway.

NF-B/RELA regulates the WNT pathway in EMT
The WNT signaling pathway mediates cell-cell interactions controlling pattern formation in embryogenesis. The pathway consists of the WNT ligands encoded by 19 distinct human genes (14) whose products bind to a heterodimeric plasma membrane receptor consisting of Frizzled (FZ) and LRP5/6 proteins. WNT ligands activate two opposing pathways, depending on the ligand and context of stimulation. WNT3 signaling stabilizes CTNNB, functioning as a cofactor of T cellregulatory factor (TCF), nuclear factor of activated T cells (NFAT) 5, and JUN. In contrast, WNT5B activates a noncanonical pathway through receptor tyrosine kinase-like orphan receptor 2 (ROR2) important in oncogenesis (13).
Experimentally observed changes in mRNA expression of the WNT pathway components are shown visually by hierarchical clustering of the RNA-Seq data. In WT cells, we observed that TGF␤ stimulation down-regulated FZ component 9 receptor and WNT3a and -4 ligands and induced a strong up-regulation of WNT3A, PRICKLE2, and the NFATC isoforms within 1 day of stimulation (Fig. 6A). A later induction of WNT5B and -7A ligands was observed after 3 day of TGF␤ stimulation. Interestingly, these major components of the WNT pathway required RELA; in the RELA KD cells, the TGF␤-induced up-regulation of WNT3, -5B, NFAT, and PRICKLE2 was blocked (Fig. 6A).
To confirm this observation, we conducted independent Q-RT-PCR assays of WNT3, -5B, and PRICKLE2 expression. Both WNT3 and PRICKLE2 were highly induced by TGF␤ in WT cells, peaking after 1 day of stimulation at 4.5-and 9-fold over that of untreated cells, respectively (Fig. 6B). By contrast, in the RELA KD cells, the TGF␤ induction of WNT3 was significantly reduced at all of time points (p Ͻ 0.01), whereas PRICKLE2 expression was reduced after 1 and 3 days of TGF␤ stimulation. Similarly, WNT5B was highly induced in WT cells by 3 days of TGF␤ treatment, increasing to 80-fold, and significantly inhibited after 1 and 3 days of stimulation in the RELA KD cells. These data confirm NF-B/RELA dependence of these members of the WNT signaling pathway.
To provide a mechanism for their RELA dependence, we used two-step cross-linking ChIP (XChIP) assays to determine whether RELA directly interacted with these genes within their native chromatin environment. Q-gPCR assays were developed for high-affinity NF-Bbinding sites in the proximal regions of WNT3, WNT5B, and PRICKLE2 promoters (Table S3), cor-responding to RELA peaks identified in ENCODE ChIP-Seq data (Fig. S2). We observed that RELA was strongly bound to these promoters in unstimulated cells by comparison of the signal from RELA immunoprecipitation with that of control IgG where a 10-fold increase in promoter signal was observed by Q-gPCR (p Ͻ 0.01; Fig. 6C). A 128-fold difference in signal intensity was observed when we examined the enrichment of an NF-B-binding site relative to a distant region of the gene lacking an NF-B-binding site (not shown). The effect of TGF␤ stimulation did not significantly change the amount of RELA binding to the promoter over that of control cells, although an increase in RELA binding was observed by concomitant stimulation with TNF, a ligand that produces a much more robust RELA translocation than does TGF␤ and potentiates TGF␤induced EMT (17,18,32). For each treatment condition, less RELA binding was observed in the RELA KD versus WT cells, further indicating assay specificity (Fig. 6C). The inability to detect the 2-fold change in RELA abundance by 3 days of TGF␤ stimulation reflects a limitation of the XChIP assay to small differences in RELA abundance/binding. These data indicate that the mechanisms of TGF␤-induced RELA activation early in the course of EMT are distinct from the robust nuclear translocation characteristic of the TNF-induced canonical pathway.
To further support the conclusion that RELA was required for WNT signaling in mesenchymal transition, we performed CFM experiments in control or TGF␤-stimulated WT or RELA KD cells. In the absence of stimulation, no WNT5B antigen was detected; however, TGF␤ induced a dramatic up-regulation of WNT5B immunoreactivity within the nucleus (Fig. 6D, top  row). This nuclear localization is consistent with WNT5B distribution within the kidney epithelium undergoing EMT (33). By contrast, in RELA KD cells, WNT5B immunoreactivity was undetectable in the absence or presence of TGF␤ stimulation (Fig. 6D). Similarly, CTNNB1 immunoreactivity was strongly induced in the cytoplasm of TGF␤-stimulated WT cells, whereas its expression was also completely inhibited in RELAsilenced cells (Fig. 6D, middle row). Finally, PRICKLE2 immunoreactivity was TGF␤-inducible in WT cells in a cytoplasmic vesicle-bound distribution, consistent with its localization to Golgi bodies (Fig. 6D, bottom row). No PRICKLE2 immunostaining was detected in the RELA KD cells (Fig. 6D).

Independent validation of RELA dependence
To independently confirm that the WNT pathway is RELAdependent, we performed validation experiments using CRISPR/ Cas9 -mediated RELA silencing by genome editing. In these cells, RELA protein and mRNA are undetectable (34) and fail to induce MMP9 after TGF␤ stimulation (see Figs. 2E and 7A). Because MMP9 is a hallmark gene of EMT, these data confirm the RELA dependence for the initiation of EMT. Next, WT or RELA Ϫ/Ϫ hSAECs were TGF␤-stimulated, and the expression Note that the major GO term is transcription factor activity. C, E, and G, GSEA for RELA-dependent genes in control (day 0 (d0)) and 1 and 3 days after stimulation. D, F, and H, pathway enrichment analysis. Shown are the top 10 overrepresented pathways for the differentially expressed genes in RELA KD cells in control (day 0) and at 1 and 3 days as indicated. CCKR, cholecystokinin receptor; GHRH, growth hormone-releasing hormone; EGFR, epidermal growth factor receptor; PDGF, platelet-derived growth factor; ACHR, acetylcholine receptor; FDR, false discovery rate.

NFB is a master regulator of EMT autocrine loops
of WNT ligand and downstream genes was measured by Q-RT-PCR. In a pattern similar to that observed in the RELA KD cells, the TGF␤-induced expression of WNT5B, AXIN2, and Dickkopf (DKK) was absent in RELA Ϫ/Ϫ hSAECs (Fig. 7B). Impor-tantly, previous work has shown that AXIN2 expression is induced by WNT signaling, indicating that TGF␤ functionally activates the WNT pathway (35). Finally, we confirmed that the up-regulation of the ZEB1-and -2 genes is RELA-dependent as

NFB is a master regulator of EMT autocrine loops
basal and TGF␤-induced expression is reduced in RELA Ϫ/Ϫ hSAECs (Fig. 7C).

NF-B/RELA is upstream of the master EMT transcriptional regulators JUN and SNAI1
The pEMT-to-full mesenchymal transition depends on the activation of JUN and the SNAI1-ZEB1 master transcription factors (9,11). In this regard, SNAI1 is one of the early expressed core EMT regulators important in triggering ZEB1 (10). SNAI1 mRNA expression was rapidly activated 7-fold within 1 day of TGF␤ stimulation, peaking at 9-fold in WT cells, and was only transiently induced 5-fold in the RELA KD cells, falling to Ͻ2-fold after 3 days (Fig. 8A). A strong RELA signal was found associated with the endogenous SNAI1 promoter in control cells; TGF␤ treatment slightly, but statistically signifi-cantly, increased this binding (Fig. 8B). In contrast, RELA binding was reduced for each treatment condition in the RELA KD cells (Fig. 8B).
Because we noted that JUN was significantly affected by RELA silencing, we also sought to confirm its direct regulation by NF-B. We confirmed a dramatic 30-fold up-regulation of JUN mRNA expression within 1 day of TGF␤ stimulation (Fig.  8C). The TGF␤ induction of JUN was significantly blunted in the RELA KD cells (Fig. 8C). We found a pattern of basal and TGF␤-increased RELA binding on the JUN promoter similar to that observed for SNAI1 in XChIP (Fig. 8D).

Discussion
Understanding the cellular pathways that trigger EMT is fundamental to the pathogenesis of asthma (6), chronic obstructive Figure 6. Effect of RELA KD on the WNT signaling pathway. A, hierarchical clustering of temporal expression of WNT pathway members. Mean transcripts per million counts for each experimental condition were z-score-normalized by row and subjected to hierarchical clustering using Euclidean distance. Gene names are given by universal gene symbol. B, Q-RT-PCR assay was independently conducted to confirm RNA-Seq expression profiles for WNT3, WNT5B, and PRICKLE2 mRNA. Shown are individual mRNA abundance in replicates normalized to POLB. *, p Ͻ 0.05; **, p Ͻ 0.01, ANOVA, post hoc Tukey's test. C, XChIP assay was conducted to quantify RELA binding to the proximal promoter. Stimulation (Stim) conditions are control, TGF␤, or TGF␤ ϩ TNF. Antibodies used are normal anti-rabbit IgG or anti-RELA. Q-gPCR was conducted using gene-specific primer pairs (Table S3) corrected for input (70). Data are expressed as -fold change relative to control WT cells immunoprecipitated with IgG. Shown are individual replicates. D, expression of WNT5B, CTNNB1, and PRICKLE2 was determined by CFM. Ab detection was by secondary Alexa Fluor 488 (green); nuclei are stained with DAPI (blue). CFM was taken at 63ϫ. d, day.

NFB is a master regulator of EMT autocrine loops
pulmonary diseases (36), and pulmonary fibrosis (37), diseases that affect over 300 million people worldwide (38). In the differentiated vertebrate airway, type II EMT plays a central role in the response to injury. Mucosal injury releases growth factors triggering the mesenchymal dedifferentiation program that produces extensive cytosolic restructuring resulting in the loss of apical-basal polarity, dissolution of adherens junctions, enhanced motility, and expression of mesenchymal fibrotic genes (3). These cells acquire stem cell-like properties, become resistant to reactive oxygen species stress, and become motile. This phenotypic change permits the transitioned mesenchymal cell to repopulate regions of denuded epithelium, promoting tissue repair and fibrosis (3). Currently, we understand that EMT is regulated by sequential activation of autocrine loops that are responsible for stepwise transition of epithelial to pEMT and subsequently to the full mesenchymal phenotype. However, the GRNs of the EMT program are highly cell typedependent, initially determined through the interaction of TGF␤-induced SMAD3 with tissue-specific master transcription factors and modified by oncogenic GTPase signaling proteins (8,17).
Our studies of type II EMT in primary human airway cells and those of type III EMT in transformed cells have demonstrated that NF-B action is important for both types of mesenchymal transition (18, 39 -41). NF-B has multiple regulatory interfaces with EMT at the transcriptional and posttranscriptional levels (42). In this study, we used RNA-Seq to systematically identify the early transcriptional events under NF-B control. We found that RELA controls expression of a significant fraction of those genes regulated by TGF␤. We have discovered and validated in two independent models of RELA silencing that RELA coordinately regulates ligands, receptors, and downstream targets of the WNT pathway as well as JUN and SNAI1-ZEB1 master transcription factors.

TGF␤ induces a mode of NF-B activation distinct from the canonical pathway
NF-B is a pleiotropic, cytokine-inducible transcription factor regulated by post-translational modifications (43,44). In resting cells, NF-B/RELA is retained in the cytoplasm, but a small fraction undergoes continuous nuclear-cytoplasmic shuttling (45,46). The well-understood canonical NF-B pathway is rapidly (within minutes) triggered by cytokines to phosphorylate and ubiquitylate IB␣, resulting in RELA cytoplasmic release. RELA release is followed by a second activating step involving serine phosphorylation and acetylation, triggering complex formation with BRD4/CDK9 coactivator (47). Interestingly, our studies suggest that TGF␤ activation of NF-B is distinctly slower than the fast activation by cytokines where TGF␤-stimulated RELA translocation occurs much more slowly with activating Ser-276 phosphorylation even later.
In addition to the differences in kinetics, the 2-fold magnitude of RELA translocation in response to TGF␤ is distinctly less than that produced by canonical NF-B translocation (see  (Table S3) corrected for input (70). *, p Ͻ 0.05; **, p Ͻ 0.01, ANOVA, post hoc Tukey's test. Refs. 17 and 30). These data suggest that the major effect of TGF␤ is post-translational regulation of the RELA phosphorylation state. Our previous work has suggested that RELA translocation is indirectly mediated by a paracrine factor secreted into the conditioned medium (18). The identity of this factor(s) is currently unknown and will require further study.

NF-B mediates the WNT morphogen signaling
The mesenchymal transition is the product of sequential cascades involving the transition from epithelial to an uncommitted "partial" mesenchymal state followed by commitment to an irreversible mesenchymal state (11). These cell-state transitions are achieved by activating secondary (downstream) autocrine signaling loops, including IL6 and fibroblast growth factor (15,16). In this study, we extend our understanding to the RELA dependence of the WNT pathway in TGF␤-induced type II mesenchymal transition.
The WNT pathway is important in the formation of cell fate determination, formation of cellular polarity in organ development, and control of cellular migration (14). The inducible WNT ligands are secreted, palmitoylated glycoproteins that function as morphogens by forming concentration gradients. All WNT ligands signal by binding heterodimeric FZ receptors, triggering intracellular signaling pathways, including the canonical and noncanonical pathways with the latter composed of both the WNT/Ca ϩ and planar cell polarity pathways (48). The canonical pathway stabilizes CTNNB1, which serves as a transcriptional coactivator inducing downstream homeodomain proteins important in pattern formation. By contrast, the noncanonical pathways are CTNNB1-independent, activating the tyrosine kinase ROR2 (49) and Rho GTPase, stimulating actin polymerization and cellular migration characteristic of the mesenchymal state (13,14). The planar cell polarity pathway induces differential localization of PRICKLE2. Our RNA-Seq studies indicate that TGF␤ suppresses expression of WNT3 and -4 but activates WNT5B and -7B ligands, suggesting that the mesenchymal state shifts the modes of WNT signaling. Consistent with this hypothesis, WNT5B partially inhibits canonical signaling (13), perhaps explaining the cytoplasmic CTNNB1 accumulation observed by CFM (Fig. 6D). A complete understanding of the functional effect of this differential ligand expression will require further study. Although RELA has been previously reported to activate WNT5A (50), our data are the first to show that NF-B controls expression of WNT5B and -7A during mesenchymal transition, suggesting a much broader impact of RELA on WNT signaling.
Although our studies focus on RELA-dependent GRN in the initiation of EMT, our data permit us to infer the activity of the WNT pathway. The observations of selective up-regulation of AXIN2, a direct target of WNT signaling (35), and down-regulation of DKK, an inhibitor of the WNT signaling pathway (51), suggest that WNT signaling is activated early in the initiation of the mesenchymal transition by multiple mechanisms.

NF-B/RELA is a master regulator of EMT upstream of SNAI1 and JUN
SNAI1 and ZEB1 are master transcriptional regulators of the mesenchymal transition whose expression are dramatically up-regulated by TGF␤ (52). Currently, it is understood that SNAI1 and ZEB1 are silenced in epithelial cells by a double-negative feedback loop mediated by the microRNAs miR34 and miR200 (10,53). Under basal conditions, miR34 post-transcriptionally inhibits SNAI1 expression, whereas miR200 blocks ZEB1 expression. Equilibrium of the double-negative feedback loops is perturbed by TGF␤-induced transcriptional up-regulation of SNAI1, which binds to the promoter and inhibits miR200 expression, releasing ZEB1. Our data here indicate that RELA plays a direct role in activating this SNAI1-ZEB1 negative feedback loop. Our data demonstrate that SNAI1 is directly activated by RELA binding to its proximal promoter. Initially, the TGF␤ induction of SNAI1 after 1 day is not affected, but continued SNAI1 expression after 3 days requires RELA. This finding suggests a complex mode of SNAI1 regulation where SNAI1 is initially up-regulated directly by SMAD3 (7, 54) but later is functionally replaced by RELA for sustained expression. Our earlier studies have shown that RELA activation of SNAI1 is mediated by a RELA-dependent coactivator recruitment of the atypical histone acetyltransferase BRD4 (18). BRD4 activates SNAI1 by phosphorylation of the C-terminal repeat domain of RNA polymerase II, promoting its transcriptional elongation.
In addition to the SNAI1-ZEB1 double-feedback loop, a recent study has identified the essential role for a clique of JUNB-HNF4 -ETS transcription factors in the pEMT-mesenchymal transition (9). JUNB-HNF4 -ETS are induced by TGF␤ stimulation, physically associate, and coordinately bind to the promoters of many mesenchymal genes (9). These factors are also associated with dense clusters of H3K27Ac loci, a.k.a. "superenhancers" in the A549 carcinoma cell line. Although JUNB is not induced by TGF␤ in hSAECs, we provide evidence here that JUND, a homologous isoform with identical DNA binding specificity, is directly regulated by RELA. Collectively, these data show that RELA is upstream of these two core transcription factor loops necessary for mesenchymal transition, qualifying RELA to be considered as a "master transcriptional regulator" of EMT.
Master regulators are engaged in the coordinate regulation of cliques of downstream transcription factors by maintaining their expression by the formation of superenhancers (9,55). Our data show that RELA directly binds to the WNT5B, JUND, SNAI, and ZEB1 promoters and is involved in their constitutive expression. Previous work has shown that NF-B/RELA controls the activity of superenhancers important in innate inflammation (56). More work will be required to identify whether superenhancers mediate type II EMT and how the activated RELA complex maintains mesenchymal transition genes in open chromatin domains.
In summary, our study identifies the role of NF-B/RELA as a master regulator of mesenchymal transition in primary human airway epithelial cells. Our studies demonstrate that RELA is required for the activation of the downstream WNT, JUN, and SNAI1-ZEB1 core regulatory factors required for transition from pEMT to the fully committed mesenchymal program.

Experimental procedures hSAEC culture and induction of the mesenchymal transition
An immortalized hSAEC line established by infecting primary hSAECs with human telomerase and cyclin-dependent kinase-4 retrovirus constructs was described (20). The immortalized hSAECs were grown in small airway epithelial cell growth medium (Lonza, Walkersville, MD) in a humidified atmosphere of 5% CO 2 . For induction of EMT, TGF␤ (10 ng/ml; PeproTech, Rocky Hill, NJ) was added to growth medium for the indicated times.

RELA silencing
hSAECs expressing Dox-regulated RELA shRNA were established using puromycin selection. Briefly, TRIPZ Tet-On inducible lentiviral RELA shRNA plasmid (Dharmacon, Thermo Fisher Scientific, Lafayette, CO) was reverse transfected with a packaging construct into BOS23 cell lines, and supernatants were harvested. Naïve hSAECs were infected and selected for puromycin resistance (4 g/ml). Puromycin-resistant hSAECs were pooled and characterized. RELA depletion was produced by addition of Dox (2 g/ml; Sigma-Aldrich) to the culture medium for 5 days.
RELA deletion was performed by CRISPR/Cas9 genome editing. A 20-nucleotide seed sequence in exon 4 of the RELA gene was used to introduce a stop codon that induces nonsensemediated RNA decay. RELA Ϫ/Ϫ hSAECs lack detectable RELA mRNA and protein (34,57).

Q-RT-PCR
Total RNA was extracted using acid guanidium phenol (TRI Reagent, Sigma-Aldrich). For gene expression analyses, 1 g of RNA was reverse transcribed in a 20-l reaction mixture using SuperScript III (58,59). Relative changes in gene expression were quantified relative to control ␤-polymerase (POLB) transcripts using the ⌬⌬Ct method (60). The forward and reverse gene-specific Q-RT-PCR primers are listed in Table S1.

Subcellular fractionation and Western immunoblot analyses
Nuclear and cytoplasmic proteins were fractionated as described previously (61). For Western blotting, equal amounts of nuclear protein were fractionated by SDS-PAGE and transferred to polyvinylidene difluoride membranes. The membranes were incubated with affinity-purified rabbit polyclonal antibodies to RELA (Santa Cruz Biotechnology). Washed membranes were then incubated with IRDye 800 -labeled antirabbit IgG antibodies (Rockland Immunochemicals, Gilbertsville, PA), and immune complexes were quantified using the Odyssey IR Imaging System (LI-COR Biosciences, Lincoln, NE).

Immunostaining and CFM
hSAECs incubated in the absence or presence of Dox (2 g/ml, 5 days) were stimulated in the absence or presence of TGF␤ (10 ng/ml, 3 days) and replated on rat-tail collagentreated coverslips (Roche Applied Sciences). Cells were fixed with paraformaldehyde in PBS (4%, v/v; 30 min) and incubated in 0.1 M NH 4 Cl in PBS for 10 min. For visualization of cytoplasmic F-actin distribution, fixed cells were stained with Alexa Fluor 488 -or Alexa Fluor 568 -phalloidin (Life Technologies) and counterstained with 4Ј,6-diamidino-2-phenylindole (DAPI) for visualizing the nucleus. The cells were then imaged with a Zeiss LSM510 fluorescence confocal microscope at a magnification of 63ϫ.

RNA extraction and qualification
Total cellular RNA was extracted using TRI Reagent according to the manufacturer's recommendations. RNA was quantified spectrophotometrically using a NanoDrop ND-1000 (NanoDrop Technologies). Quality of the purified RNA was assessed by visualization of 18S and 28S RNA bands using an Agilent BioAnalyzer 2100 (Agilent Technologies, CA). The resulting electropherograms were used in the calculation of the 28S/18S ratio and the RNA integrity number.

Library construction and sequencing
A RiboZero Gold kit (Illumina, San Diego CA) was used to remove rRNA from 1 g of total RNA. First and second strand synthesis, adapter ligation, and amplification of the library were performed using the Illumina TruSeq v2 Stranded RNA Sample Preparation kit under conditions prescribed by the manufacturer (Illumina). Samples were tracked using "index tags" incorporated into the adapters. Library quality was evaluated using an Agilent DNA-1000 chip on an Agilent 2100 Bioanalyzer. Quantification of library DNA templates was performed using Q-PCR and a known-size reference standard. Cluster formation of the library DNA templates was performed using the TruSeq PE Cluster kit v3 (Illumina) and the Illumina cBot work station using conditions recommended by the manufacturer. Template input was adjusted to obtain a cluster density of 700,000 -900,000/mm 2 . Paired-end 50-base sequencing by synthesis was performed using TruSeq SBS kit v3 (Illumina) on three lanes of an Illumina HiSeq 1500 using the manufacturer's protocols.

RNA-Seq analysis
The raw next-generation sequencing analysis passed 11 separate quality analyses by FastQC (http://www.bioinformatics. babraham.ac.uk/projects/fastqc/, 4 updated November 2013, accessed June 6, 2014). Reads were aligned to the human GRCh38 reference genome using STAR, an ultrafast splicingaware alignment program, version 2.5.1b (63). Read counts per gene were obtained with the featureCounts function (v1.5.0-p1) of the subread software package (64) using the GENCODE version 24 basic annotation file. Counts were input into the DESeq2 differential gene expression analysis program (version 1.12.1) (65). Factors were RELA level (KD versus control) and TGF␤ treatment (0, 1, and 3 days). Differentially expressed genes were called with a log 2 -fold change greater than 2 or less than Ϫ2 and a Benjamini-Hochberg-adjusted p value (false discovery rate) of less than 0.01.

Two-step XChIP
XChIP was performed as described previously (70). hSAECs (4 ϫ 10 6 -6 ϫ 10 6 /100-mm dish) were washed twice with PBS. Protein-protein cross-linking was first performed with disuccinimidyl glutarate (2 mM, 45 min; Thermo Fisher Scientific) followed by protein-DNA cross-linking with formaldehyde (1%, 15 min; Thermo Fisher Scientific). Chromatin was fragmented by sonication, and equal amounts were immunoprecipitated overnight at 4°C with 4 g of the indicated Ab in ChIP dilution buffer. Immunoprecipitates were collected with 40 l of protein A magnetic beads (Dynal Inc.), washed, and eluted in 250 l of elution buffer for 15 min at room temperature. Samples were de-cross-linked in 0.2 M NaCl at 65°C for 2 h. The precipitated DNA was phenol-chloroform-extracted, ethanolprecipitated, and dried.
Quantitation was performed by Q-gPCR using gene-specific primers (Table S3). These PCR primers were designed to amplify NF-B-binding sites identified by TRANSFAC analysis and confirmed by ChIP-Seq tracks from ENCODE. Data are presented as -fold change relative to RELA IPs of untreated (control) cells. For this calculation, for each gene-specific primer set, the median Ct value of input was subtracted from the Ct value produced by RELA IP for each cellular treatment (⌬Ct). For each cellular treatment, the -fold change relative to control cells was calculated by subtracting the ⌬Ct from the ⌬Ct of the RELA IP from control cells.

Statistical analysis
Differences among the means were tested for significance by ANOVA with Fisher's least significance difference test. Individual comparisons were made using post hoc Tukey's test. Significance was reached when the p value was Ͻ0.05. Values presented are means Ϯ S.D.