Deconvolution of multiplexed transcriptional responses to wood smoke particles defines rapid aryl hydrocarbon receptor signaling dynamics

The heterogeneity of respirable particulates and compounds complicates our understanding of transcriptional responses to air pollution. Here, we address this by applying precision nuclear run-on sequencing and the assay for transposase-accessible chromatin sequencing to measure nascent transcription and chromatin accessibility in airway epithelial cells after wood smoke particle (WSP) exposure. We used transcription factor enrichment analysis to identify temporally distinct roles for ternary response factor–serum response factor complexes, the aryl hydrocarbon receptor (AHR), and NFκB in regulating transcriptional changes induced by WSP. Transcription of canonical targets of the AHR, such as CYP1A1 and AHRR, was robustly increased after just 30 min of WSP exposure, and we discovered novel AHR-regulated pathways and targets including the DNA methyltransferase, DNMT3L. Transcription of these genes and associated enhancers rapidly returned to near baseline by 120 min after exposure. The kinetics of AHR- and NFκB-regulated responses to WSP were distinguishable based on the timing of both transcriptional responses and chromatin remodeling, with induction of several cytokines implicated in maintaining NFκB-mediated responses through 120 min of exposure. In aggregate, our data establish a direct and primary role for AHR in mediating airway epithelial responses to WSP and identify crosstalk between AHR and NFκB signaling in controlling proinflammatory gene expression. This work also defines an integrated genomics-based strategy for deconvoluting multiplexed transcriptional responses to heterogeneous environmental exposures.

The heterogeneity of respirable particulates and compounds complicates our understanding of transcriptional responses to air pollution. Here, we address this by applying precision nuclear run-on sequencing and the assay for transposaseaccessible chromatin sequencing to measure nascent transcription and chromatin accessibility in airway epithelial cells after wood smoke particle (WSP) exposure. We used transcription factor enrichment analysis to identify temporally distinct roles for ternary response factor-serum response factor complexes, the aryl hydrocarbon receptor (AHR), and NFκB in regulating transcriptional changes induced by WSP. Transcription of canonical targets of the AHR, such as CYP1A1 and AHRR, was robustly increased after just 30 min of WSP exposure, and we discovered novel AHR-regulated pathways and targets including the DNA methyltransferase, DNMT3L. Transcription of these genes and associated enhancers rapidly returned to near baseline by 120 min after exposure. The kinetics of AHR-and NFκB-regulated responses to WSP were distinguishable based on the timing of both transcriptional responses and chromatin remodeling, with induction of several cytokines implicated in maintaining NFκB-mediated responses through 120 min of exposure. In aggregate, our data establish a direct and primary role for AHR in mediating airway epithelial responses to WSP and identify crosstalk between AHR and NFκB signaling in controlling proinflammatory gene expression. This work also defines an integrated genomics-based strategy for deconvoluting multiplexed transcriptional responses to heterogeneous environmental exposures.
Particulate matter (PM) air pollution is estimated to cause over 4 million premature deaths annually on a worldwide basis (1). In addition to premature death, air pollution has been associated with a wide variety of health effects including acute and chronic respiratory symptoms, inflammation, exacerbations, and morbidity (2)(3)(4). However, no specific therapies have been developed to address either acute or chronic health impacts of air pollution, and our understanding of individual versus population risk from exposure is fragmentary. Compounding this, respirable PM pollution is extremely heterogeneous (5). Thus, there is an important need to understand the effects of PM at the molecular level and to distinguish the cellular pathways that are regulated in response to different pollution exposures that are individually and in aggregate comprised of complex mixtures of chemical and particulates.
Respirable pollution is deposited along the airway and is known to exert a variety of effects on airway epithelial cells in a process that results in significant changes in gene expression (6,7). This has been investigated directly in human exposure studies (8); however, these processes are more frequently modeled using cell culture systems and standardized sources of pollution. For example, several pathways have been identified as regulating gene expression responses to laboratory-generated diesel exhaust particles (9)(10)(11), including activation of the inflammasome, NFκB, and aryl hydrocarbon receptor (AHR) signaling (12)(13)(14). Inflammasome signaling, which results in NFκB activation, is associated with induction of cytokines and inflammatory responses (15). AHR, which is a basic helix-loop-helix-Pas domain transcription factor, is activated by a variety of potential AHR ligands that are byproducts of combustion (16). Upon ligand binding, AHR translocates to the nucleus and controls gene expression through binding to xenobiotic response elements (XREs) found in regulatory regions of canonical target genes (17)(18)(19)(20), including detoxification enzymes such as CYP1A1, and the AHR repressor, AHRR. Functional XREs are found in regulatory regions for NFκB target genes, suggestive of a proinflammatory role in some contexts (21,22). However, AHR also induces the expression of interleukin-24 (IL-24), a cytokine that represses NFκB (22). Thus, in addition to heterogeneity within air pollution, cellular responses to air pollution can encompass multiplexed transcriptional cascades whose crosstalk potentially results in both proinflammatory and antiinflammatory effects.
One area of growing importance that further highlights challenges in contending with particle and chemical heterogeneity and cellular responses is wildfire pollution. As a consequence of climate changes, wildfires are increasing in size, prevalence, and impact on heavily populated areas (23,24). Wood smoke particles (WSPs) derived from controlled combustion of wood have been used as a model to study cellular effects of wildfire smoke exposure (25)(26)(27). These WSPs are composed of a heterogeneous mix of compounds including fine and coarse PM, heavy metals, polyaromatic hydrocarbons (PAHs), and volatile organic compounds, which depend on the source of combustible material, similar to heterogeneity within bona fide wildfire pollution (26). Experimental WSP inhalation in healthy individuals recruits inflammatory cells, stimulates oxidative damage, and impairs immune responses to viral infections. WSP exposure in cell culture models causes cytotoxicity, genotoxicity, cellular stress, formation of reactive oxygen species, and increased expression of proinflammatory cytokines (27)(28)(29). However, the molecular and transcriptional control of these effects has not been fully elucidated.
In this study, we developed a general approach to define the primary transcriptional mediators that result from exposing cultured human airway epithelial cells to a complex pollutant mix. We exposed Beas-2B airway epithelial cells to particles and compounds generated from controlled combustion of white oak. We assayed genome-wide nascent transcriptional responses and chromatin accessibility changes after 30 and 120 min of exposure to the white oak-derived WSP. We applied bioinformatics tools to cluster expression and epigenetic responses into distinct temporal patterns and to identify transcription factors that mediate these effects. We employed chromatin immunoprecipitation (ChIP)-quantitative PCR (qPCR), siRNA-mediated gene knockdown, and primary airway epithelial cells to validate transcriptional effectors discovered through this genomics-based discovery method.

WSP exposure causes rapid and temporally dynamic changes in gene transcription
As an initial step in developing a system to identify transcriptional mechanisms resulting from exposure to air pollution, we exposed Beas-2B cells, a transformed human bronchial epithelial cell line, to WSP in submerged culture at varying concentrations of 10, 100, and 1000 μg/ml (Fig. 1A). WSP is an established model of wildfire pollution that is known to include significant chemical and particulate heterogeneity. We analyzed changes in gene expression after 2, 4, and 24 h of exposure using qRT-PCR, focusing on previously identified targets of PM in airway cells and NFκB target genes (30)(31)(32)(33). We observed rapid induction of three candidate exposure-regulated genes (Fig. 1B), with peak expression observed after just 2 h of treatment. Thus, airway epithelial cells exhibit prompt, dynamic, and dose-dependent changes in gene expression in response to WSP.
To explore the earliest transcriptional responses to WSP and determine underlying mechanisms, we conducted nascent transcript analysis of Beas-2B cell responses to WSP using precision nuclear run-on sequencing (PRO-Seq) (34). PRO-Seq measures gene transcription with high temporal resolution based on RNA polymerase II activity and also allows for quantification of enhancer activity based on bidirectional signatures of enhancer RNA (eRNA) transcription (35,36). For these studies, we selected the highest tested dose of WSP (1000 μg/ml) in order to elicit maximal transcriptional responses for further bioinformatics analyses. To define direct transcriptional responses to WSP and their initial dynamics, we analyzed nascent transcription after 30-and 120-min exposure times relative to treatment with vehicle (PBS). We initially focused our analysis on changes in gene transcription in relationship to exposure time. As visualized in volcano plots ( Fig. 2A), at the 30-min time point, increased transcription was the dominant response to WSP exposure, with many of these rapidly induced targets exhibiting reduced expression at 120 min relative to 30 min. We subsequently defined three distinct temporal clusters of transcriptional induction in response to WSP, as indicated by Venn diagrams in Figure 2B and detailed in the Experimental procedures section, encompassing early peak, early plateau, and late clusters of differentially regulated gene transcripts. Complete gene lists for each cluster are presented in Tables S3-S5, respectively. PRO-Seq data were visualized in the Integrative Genomics Viewer (IGV) genome browser (37), and representative data tracks for each temporal cluster are shown in Figure 2C.
Next, to determine whether the different clusters represent distinct biologic processes regulated by WSP as a function of time, we performed functional annotation of the clusters using DAVID (38). Although there were annotation differences between all the three clusters ( Fig. 2D; complete output for each cluster is listed in Tables S6-S8), comparison of the late cluster to the early peak and early plateau annotations showed significant divergence in enriched pathways. In addition, annotation of the early peak cluster identified enrichment of AHR-related pathways including cytochrome P450 metabolism (p = 6.438e-5) and PAS domain-containing transcription factors (p = 3.399e-5). Functional annotation of early plateau cluster genes also identified enrichment of PAS domain-containing transcription factors (p = 7.424e-4) and NFκB signaling pathways (p = 9.775e-3). In aggregate, our analysis of nascent transcription indicates that WSP exposure causes a rapid remodeling of airway epithelial gene expression and associated cellular processes. These changes are temporally dynamic, and ontology analysis implicates AHR and NFκB as potential transcriptional mediators of early effects of WSP exposure on transcription.

eRNA transcription is regulated by WSP
To further study mechanisms underlying regulation of gene expression by WSP, we next analyzed changes in enhancer expression by utilizing Tfit to identify bidirectional signatures of RNA polymerase II activity associated with nascent eRNA and promoter transcription (39). We then used DEseq2 to determine differential expression of these bidirectional transcripts following vehicle or WSP exposure at both time points. Through this analysis, we defined clusters of bidirectional transcriptional regulation in response to WSP analogous to the clusters observed for gene transcription (Fig. 3A). Examples of these distinct bidirectional transcript response patterns to WSP as visualized in IGV are shown in Figure 3B. Complete lists of differentially transcribed bidirectionals from early peak, early plateau, and late temporal clusters, as well as genomic coordinates for all identified bidirectional transcripts, are provided in Tables S9-S12, respectively. Taken together, these data indicate that dynamic changes in airway epithelial cell Demultiplexing rapid effects of wood smoke on transcription gene expression in response to WSP are associated with rapid changes in enhancer activity.
Dynamic changes in enhancer activity generally result from differential transcription factor binding within the enhancer. Therefore, to identify transcription factors that control the primary response to WSP exposure, we interrogated enhancer regions whose activity, based on PRO-Seq signatures, changed dynamically in response to WSP exposure. To accomplish this, we identified the origin of bidirectional transcription and applied the transcription factor enrichment analysis (TFEA) tool to these regions to identify differential motif enrichment for specific transcription factors (40). TFEA quantifies the degree of colocalization of transcription factor motif instances within the center of specific genomic regions, which in this case are sites of bidirectional transcription. We therefore compared changes in enrichment scores between vehicle and 30, 30, and 120 min, and vehicle and 120 min, analogous to the three temporal clusters we previously defined based on induction of transcription. Motifs for the ternary response factor-serum response factor (TCF-SRF) complex, which is activated in response to extracellular signal-regulated kinase and mitogen-activated protein kinase signaling (41,42), were identified as significantly enriched in the 30-min dataset in comparison to vehicle (e.g., Fig. 3C, top), suggesting a rapid kinase-mediated transcriptional response to WSP. Comparison of 30 to 120 min showed that enrichment of TCF and SRF motifs returned to baseline by 120 min, whereas EGR family motifs were enriched at 120 min in comparison to vehicle (Fig. 3C, bottom). These data suggest that dynamic changes in transcription factor utilization underlie different kinetics in transcriptional responses to WSP. Significantly enriched and depleted motifs for each treatment comparison are summarized in Table S13.
Changes in chromatin structure identify AHR as driving early transcriptional responses to WSP Some transcription factors are more associated with chromatin remodeling than directly altering RNA polymerase A, Venn diagrams indicating number of differentially regulated bidirectional transcripts (based on p adj < 0.05) meeting indicated temporal clustering criteria, as described for Figure 2B. B, IGV-visualized PRO-Seq tracks of Tfit-called bidirectionals (boxed in black and magnified below each panel) and nearest genes exhibiting a similar pattern of regulation as representative examples of each temporal cluster characterized in (A). C, motif displacement distributions of significantly enriched transcription factor-binding motifs within early peak (top) versus late (bottom) clusters of differentially transcribed bidirectionals, as identified using TFEA. Each column represents the frequency of the indicated motif instance at the specified distance from the bidirectional center (labeled 0), where color specifies a normalized frequency relative to the entire 3-kb region. Darker colors indicate greater enrichment on a 0 to 1 scale. IGV, Integrative Genomics Viewer; PRO-Seq, precision nuclear run-on sequencing; TFEA, transcription factor enrichment analysis; WSP, wood smoke particle. activity, and TFEA can also be applied to regions based on changes in chromatin structure to identify these regulatory pathways (40). Therefore, to determine effects of WSP exposure on chromatin status, we performed the assay for transposase-accessible chromatin using sequencing (ATAC-Seq) in Beas-2B cells following 30 and 120 min of WSP exposure (43,44). Based on MACS2-defined peaks (45), we partitioned the ATAC-Seq data temporally into the timedependent clusters we had previously defined for nascent transcription (Fig. 4A), noting that the number of new peaks occurring after 120 min of WSP was significantly greater than the number of new peaks at 30 min. Examples of loci that fall into each cluster are shown in Figure 4B. These data indicate that WSP causes significant chromatin remodeling, and the extent of chromatin remodeling in response to WSP increases over time.
As a further approach to define transcriptional mediators of genomic responses to WSP, we interrogated 3 kb regions centered on MACS2-defined ATAC-Seq peaks using TFEA (46) to compare motif enrichment based on treatment and time. This analysis identified strong enrichment for the AHR nuclear translocator motif, the obligate AHR dimerization partner. This enrichment was evident at 30 min but decreased to baseline by 120 min (Fig. 4C, top). In contrast, differential motif enrichment for RELA, a member of the NFκB complex, was significantly increased at 120 min but not 30 min (Fig. 4C,  bottom). These data strongly implicate AHR signaling as mediating early transcriptional responses to WSP in association with rapid chromatin remodeling, whereas effects of NFκB appear to peak at a later time.

Novel targets of AHR signaling defined through integrated genomics
In light of the known association of other forms of PM pollution with AHR signaling (12), and the identification of AHR through applying unbiased bioinformatics approaches to our genome-wide data, we scrutinized the early peak cluster for canonical AHR target genes. We found a significant increase in transcription of numerous AHR targets at 30 min of Figure 4. ATAC-Seq demonstrates consistent temporal dynamics between chromatin accessibility and transcriptional responses to WSP and predicts key regulatory factors and pathways. A, ATAC-Seq was performed following exposure of Beas-2B cells to vehicle or WSP for 30 or 120 min. Venn diagrams indicate number of differentially regulated ATAC-Seq peaks (based on p adj < 0.05) meeting the indicated temporal clustering criteria, as described previously. B, IGV-visualized ATAC-Seq tracks with MACS2-called peaks (boxed in black and magnified below each panel) and nearest genes exhibiting a similar pattern of regulation as representative examples of each temporal cluster characterized in (A). C, motif displacement distributions, as described for Figure 3C, from TFEA of differentially regulated MACS2-called ATAC-Seq peaks representative of early peak (top) and late (bottom) differential motif enrichment. ATAC-Seq, assay for transposase-accessible chromatin using sequencing; IGV, Integrative Genomics Viewer; WSP, wood smoke particle.
WSP exposure that sharply decreased after 120 min ( Table 1). Examples of integrated PRO-Seq and ATAC-Seq data for several established AHR targets, such as CYP1B1, and their associated enhancers visualized in IGV are shown in Figure 5A. MatInspector analysis of each of these dynamically regulated enhancer regions revealed matches for the canonical AHR/ AHR nuclear translocator response element (47), implicating AHR as likely mediating these rapid transcriptional effects. We then identified additional genes of interest demonstrating early peak kinetics (Table 1) and interrogated these genes for canonical AHR binding sites in regulatory regions. DNMTL3, which is reported to function within a protein complex to induce DNA methylation (48), is among the putative novel AHR targets we identified in this manner.
To definitively establish that AHR directs rapid and transient effects of WSP on airway epithelial gene expression, we performed ChIP-qPCR assays for AHR occupancy at several putative target loci. First, we confirmed by Western blot that AHR is expressed in this cell type (Fig. 5B). Then, we performed AHR ChIP-qPCR in Beas-2B cells following WSP exposure for 30 and 120 min. This demonstrated a gain in AHR occupancy in regulatory regions for canonical targets such as CYP1B1, as well as several novel targets, including DNMT3L, following exposure to WSP for 30 min (Fig. 5C). AHR occupancy was decreased relative to the 30-min time point after 120 min, consistent with the enhancer transcriptional signatures observed in the PRO-Seq datasets (Fig. 5A). Thus, AHR directly and rapidly activates gene expression in response to WSP, however, AHR occupancy and direct effects on gene expression peak before 2 h.
2,3,7,8-Tetrachlorodibenzo-p-dioxin exposure confirms novel AHR targets in Beas-2B and primary airway epithelial cells WSP contains a range of potential AHR ligands and whether the rapid AHR signaling response observed at 1 mg/ml of WSP is relevant to established AHR agonists is not clear. Therefore, we exposed Beas-2B cells to 2,3,7,8-tetrachlorodibenzo-pdioxin (TCDD), a prototypical ligand for AHR (20). We compared AHR induction between two doses of WSP and a standard 10 nM TCDD concentration (49). This demonstrated that expression of several, but not all, early peak target genes is dependent on AHR activation and occurs in response to standard doses of TCDD and WSP (Fig. 6A). To determine whether this also occurs in primary cells, we exposed primary small airway epithelial cells (Sm18) to TCDD and WSP. These primary cells manifested similar expression responses to these stimuli, including with respect to novel AHR target genes (Fig. 6B). AHR occupancy at regulatory regions neighboring several early peak AHR target genes was also increased with exposure to standard doses of TCDD. Our findings thus establish a physiologically diverse set of AHR targets in response to both combustion-derived AHR ligands and TCDD, with targets implicated in regulating DNA methylation (DNMT3L), inflammation (IL-24), cell-cell adhesion (LHX4), glucocorticoid responses (NR3C1), and transcriptional regulation (MAFF).

Crosstalk between AHR and NFκB signaling in response to WSP
Our data, in accordance with the literature on other sources of PM pollution, suggest that both AHR and NFκB respond to WSP, but whether there is direct crosstalk between these pathways is unclear. To study the role of AHR in inflammatory transcriptional responses to WSP, we analyzed gene expression responses to WSP in the setting of siRNA-mediated knockdown of AHR or the RELA subunit of NFκB. Using qRT-PCR, we quantified the expression of several canonical AHR targets as well as several established targets of NFκB (Fig. 7A). After WSP exposure for 2 h, AHR knockdown cells (verified by Western blotting, Fig. 7B) demonstrated decreased expression of AHRR and CYP1A1 in comparison to cells treated with si-Ctrl. There was no effect of AHR knockdown on expression of CXCL8, a prototypical NFκB target. Surprisingly, WSP-mediated increases in IL1A and IL1B expression were attenuated with AHR knockdown, indicating a direct or secondary role for AHR in inducing these cytokines in response to WSP. Knockdown of RELA resulted in an essentially reciprocal pattern, with a complete abrogation of CXCL8 induction in response to WSP, a minimal effect on the induction of AHRR and CYP1A1, and an intermediate effect on IL1A and IL1B (Fig. 7, C and D). Thus, AHR and NFκB are both required for maximal increases in IL1A and IL1B expression in response to WSP, indicating functional crosstalk between AHR and inflammatory responses to WSP. Moreover, as IL1A and IL1B induce NFκB activity, these data indicate convergence of AHR and inflammatory signaling in a positive feedback circuit that may promote further inflammation.

Discussion
Transcriptional regulations of airway epithelial responses to air pollution are challenging to define because of heterogeneity Demultiplexing rapid effects of wood smoke on transcription of the constituent respirable particulates and compounds. Through analyzing nascent transcription and performing ATAC-Seq, here we show that airway epithelial cells undergo temporal dynamic changes in gene expression and chromatin structure in response to WSP, which models both wildfire smoke and exposure-related heterogeneity. Using unbiased bioinformatics approaches, we identified distinct roles for three transcription factor families, TCF-SRF, AHR, and NFκB, in mediating responses to WSP. We subsequently defined a set of canonical AHR targets that exhibit increased transcription after 30 min of WSP exposure in association with increased activity of nearby enhancers. Using ChIP, we established AHR occupancy at canonical XREs within enhancers and promoters for several target genes, and we further confirmed a role for AHR in mediating gene expression responses through knockdown experiments and exposure to TCDD, a wellestablished AHR ligand. In aggregate, our data establish a framework for identifying transcription factors that respond to air pollution and define a direct role for AHR in inducing rapid and transient transcriptional effects in response to WSP. PAHs found within various forms of air pollution are known to function as AHR ligands (16); however, direct induction of AHR by wood smoke particulates had not been previously reported. By virtue of the high resolution afforded by integrating PRO-Seq and ATAC-Seq data, we identified canonical AHR target genes that are rapidly induced in Beas-2B cells following WSP exposure, thereby expanding our understanding of AHR signaling with respect to air pollution in airway epithelial cells. Rapidly induced targets included CYP1A1 and CYP1A2, which metabolize AHR-activating PAHs that are found in WSP; the AHR repressor, AHRR; and TIPARP, an ADP-ribosyltransferase that represses AHR activity though ribosylation (50,51). Thus, in airway epithelial cells, at least three distinct negative feedback mechanisms regulated by AHR are implicated in the rapid decrease in transcription of primary AHR targets following WSP exposure. We also Figure 5. Integrated analysis of chromatin remodeling and nascent eRNA profiling datasets identifies sites of direct AHR occupancy and functional regulation of early peak transcriptional targets of WSP. A, ATAC-Seq and PRO-Seq tracks aligned in IGV at early peak WSP response genes that are known targets of AHR signaling. Associated early peak eRNAs and ATAC-Seq peaks are boxed and magnified below each panel, with locations of matches to the consensus AHR binding motif, determined using MatInspector, indicated by light blue bars. B, Western blot verification of AHR protein expression in Beas-2B cells under basal culture conditions; GAPDH was a loading control. C, quantitative PCR primers were designed for the three regions examined in (A) plus one additional canonical (CCR7) and two novel (DNMT3L and LHX4) early peak targets of AHR, then ChIP-quantitative PCR was performed in Beas-2B cells following exposure to vehicle or WSP for 30 or 120 min. Bars represent AHR occupancy on a log 2 scale (±SD), expressed as the mean C T value at each target region relative to the geometric mean of C T values at three negative control regions (n = 4/group, *p < 0.05 for indicated comparison). AHR, aryl hydrocarbon receptor; ATAC-Seq, assay for transposase-accessible chromatin using sequencing; ChIP, chromatin immunoprecipitation; eRNA, enhancer RNA; IGV, Integrative Genomics Viewer; PRO-Seq, precision nuclear run-on sequencing; WSP, wood smoke particle. defined a novel set of genes regulated by AHR based on several criteria, including: (1) rapid and transient transcriptional induction patterns similar to that observed for canonical AHR targets with WSP exposure, (2) the discovery of WSPresponsive enhancers for these genes harboring centrally located canonical XREs, and (3) confirmation of target gene activation following TCDD exposure.
In contrast to the powerful negative feedback system that limits the duration of direct effects of AHR on gene transcription, a number of genes we identified are candidates to mediate prolonged secondary responses to WSP-mediated activation of AHR signaling. For example, transcription of the MEF2 family of transcription factors was induced after 30 min of WSP exposure. Based on TFEA indicating enrichment for motifs for each of these factors among active enhancers at 120 min versus vehicle (false discovery rate <0.05), our data implicate the MEF2 family in mediating secondary gene expression responses to WSP. In addition to this family, the bZip transcription factor, MAFF, LHX4, a lim homeobox factor, and the cytokine, IL-24, which is implicated in lung and airway remodeling in response to different inflammatory stimuli, were also among the set of AHR targets we identified (52)(53)(54)(55). Thus, the short-lived primary transcriptional response to WSP mediated by AHR includes targets that likely contribute to more sustained changes in gene expression, possibly in collaboration with other rapid transcriptional effectors of the response to WSP, such as the TCF-SRF family. Although we observed induction of AHR targets in primary human small airway epithelial cells and similar kinetics at early time points, future work is needed to more fully establish the relevance and temporal properties of AHR signaling responses to WSP in primary and differentiated airway epithelial cells cultured at air-liquid interface. How AHR signaling in response to a heterogenous exposure that results in the multiplexed transcriptional cascade activity characterized here compares with induction of AHR signaling in response to a single ligand, such as TCDD, also remains to be determined.
Persistent epigenetic changes such as DNA methylation, which can exert a memory effect on cell phenotype through relative to vehicle-treated controls (n = 4/group, *p < 0.05 versus vehicle). C, ChIP-quantitative PCR was performed in Beas-2B cells for novel early peak targets of AHR following stimulation with TCDD (10 nM). Bars represent AHR occupancy on a log 2 scale (±SD), expressed as the mean C T value at each target region relative to the geometric mean of C T values at three negative control regions (n = 4/group, *p < 0.05 for indicated comparison). AHR, aryl hydrocarbon receptor; ChIP, chromatin immunoprecipitation; TCDD, 2,3,7,8-tetrachlorodibenzo-p-dioxin; WSP, wood smoke particle.
modifying gene expression, are believed to contribute to long-term health effects of exposure to pollutants (56,57). Through PRO-Seq, we identified rapid and transient induction in response to WSP of the DNA methylase, DNMTL3 (58). Although PAH-containing pollutants have been previously reported to cause increased DNA methylation (59, 60), we are unaware of any reports that have identified specific methylation pathways that might mediate this epigenetic response to WSP or other pollutants containing AHR ligands. Whether this pathway is responsible for teratogenic and other long-term effects of AHR signaling, including altered DNA methylation, remains to be confirmed in future work.
In addition to AHR, we also identified NFκB as a mediator of transcriptional and chromatin responses to WSP in airway epithelial cells. Canonical targets of NFκB that were activated in response to WSP include TNFAIP3, CXCL8, and IL1B, among others. Knockdown of RELA confirmed that induction of TNFAIP3 and CXCL8 was largely dependent on the NFκB complex. Although other studies have implicated inflammasome activation in response to particulates as directly inducing NFκB activity (61), the mechanisms underlying NFκB activation in response to WSP remain to be elucidated. The timing of AHR and NFκB transcriptional motif utilization based on TFEA, however, suggests that NFκB-mediated responses to WSP occur in a temporally distinctive process from AHR signaling. Specifically, we observed significantly increased central NFκB motif enrichment within ATAC-Seq peaks at the 120-min time point relative to 30 min, whereas AHR motif utilization peaked at 30 min. Given the rapid induction of IL1B and CXCL8 in response to WSP, both of which are known activators of NFκB signaling, NFκB activity in response to WSP exposure is likely partially attributable to positive feedback mediated through these and other cytokines that activate NFκB.
The identification of novel AHR targets and AHR crosstalk with NFκB are two important aspects of the transcriptional response to WSP that we captured in this study. Our approach also comprehensively mapped the diverse gene targets and transcriptional regulators that respond to WSP over two time points, including defining rapid and transient transcriptional responses mediated by several families of transcription factors in response to this exposure. We employed a WSP exposure in this study both to model wildfire particulates and to represent the multiplexed transcriptional programming that arises in consequence to respirable particle and compound heterogeneity, which is a defining characteristic of many forms of air pollution. Our study thus describes a novel and generalizable approach to partition heterogeneous pollutant exposures into discrete transcriptional response pathways (Fig. 8). The strength of this method is supported by validation experiments we performed using conventional assays, such as qRT-PCR and ChIP-PCR, and a single chemical ligand, TCDD, to activate AHR directly. Although further experiments are necessary to investigate the pathways identified by our approach and their physiological consequences, in aggregate our data both define novel pathways that respond to WSP and inform a general strategy to deconvolute multiplexed transcriptional responses that arise from air pollution exposure.

Cell culture
Beas-2B immortalized human bronchial epithelial cells were obtained from American Type Culture Collection and cultured in Dulbecco's modified Eagle's medium (Corning) with L-glutamine and 4.5 g/l glucose supplemented with 10% fetal bovine serum (VWR) and 1% penicillin/streptomycin Figure 7. Airway epithelial inflammatory responses to WSP are mediated in part through complex crosstalk between AHR and NFKB. A, Beas-2B cells transiently transfected with siRNA targeting AHR (si-AHR) or a scrambled control construct (si-Ctrl) were treated with vehicle or WSP for 2 h and then assayed for indicated gene expression using quantitative RT-PCR. Bars represent mean normalized C T values on a log 2 scale (±SD) relative to si-Ctrl + vehicle-treated controls (n = 4/group, *p < 0.05 for indicated comparison). B, Western blot verification of AHR knockdown by si-AHR transfection. C, si-RELAor si-Ctrl-transfected Beas-2B cells were treated and assayed as described for (A). D, Western blot verification of RELA knockdown by si-RELA. AHR, aryl hydrocarbon receptor; WSP, wood smoke particle.
(Corning). Deidentified primary human small airway epithelial cells (smAECs; <2 mm diameter) were obtained from the National Jewish Health Biobank and plated onto an irradiated National Institutes of Health (NIH)/3T3 (American Type Culture Collection) fibroblast feeder layer in F-medium containing 1 μM Y-27632 (APEX Bio). Upon visible colony formation (7-10 days), smAECs were removed with 0.25% trypsin (Corning), plated on tissue culture dishes doublecoated with type I bovine collagen solution (Advanced Biomatrix), and grown to confluence in BronchiaLife Epithelial Medium supplemented with the complete LifeFactors Kit from Lifeline Cell Technology. All cells were maintained in 5% CO 2 at 37 C.

WSP exposure
WSPs were generated as described and generously gifted by Dr Andrew Ghio at the Environmental Protection Agency Laboratory in North Carolina (25). Briefly, wood smoke was generated by heating white oak wood on an electric heating element (Brinkmann Corporation) in a Quadra-fire 3100 woodstove (Colville). The mean diameter of the freshly generated WSP was 0.14 micron. In the freshly generated particle, metal content was negligible. Particles were obtained by its mechanical disruption from the piping above the woodstove followed by sonication in water to disaggregate the particles. The ratio of elemental to organic carbon in the WSP was 0.004 (Sunset Laboratories) suggesting elemental carbon comprised approximately 0.4% of the particle. Diluted particle was injected into a gas chromatograph (Gerstel Thermal Desorption System/Agilent 6890 equipped with an Agilent 5973 mass selective detector). Identifiable compounds included levoglucosan (25.3%; m/m). After deposition and aggregation on the piping, there were measurable concentrations of metal including calcium (7.66 ppm), iron (0.76 ppm), magnesium (0.35 ppm), and aluminum (0.31 ppm). Particles were dehydrated for transportation and resuspended in PBS. Particles were disaggregated by sonication and sterilized under UV light for 20 min prior to addition directly to cell culture media at concentrations of 0.01 to 1 mg/ml as indicated in the article.

RNA purification and qRT-PCR
Beas-2B cells or smAECs were grown to confluence in 6-well tissue culture dishes (collagen coated for smAECs, as described previously) and treated with vehicle (PBS) or WSP for 2, 4, or 24 h. Cells were harvested in TRIzol (Life Technologies) and RNA purified by PureLink RNA Mini Kit (Life Technologies) prior to qRT-PCR, performed with normalization to RPL19 as previously detailed (62). Sequences of primers used for qRT-PCR analysis are provided in Table S1.

PRO-Seq
Beas-2B cells were plated on 3 × 15 cm tissue culture dishes per treatment and grown to confluence and then treated with vehicle (PBS) or WSP for 30 or 120 min. Cells were harvested Figure 8. Genomics-based deconvolution of multiplexed transcriptional responses to complex air pollutants. WSP exposure is predicted to result in complex signal transduction events that lead to activation of multiple transcription factor pathways. To deconstruct these molecular processes, we generated paired genome-wide chromatin accessibility and nascent transcription profiles in response to two WSP exposure time points. Unbiased transcription factor enrichment analysis (TFEA) of these data identified distinct families of transcription factors that mediate changes in gene expression, enhancer activity, and chromatin structure in response to WSP. Transcription factor activity could be further clustered based on timing, and in some cases, responses were rapid and transient. This system can be applied to identify and compare epigenetic and transcription factor responses to different sources of air pollution. WSP, wood smoke particle. and nuclei prepared as reported previously (62). Aliquots of 10E6 nuclei were subjected to 3-min nuclear run-on reactions in the presence of Biotin-11-CTP (PerkinElmer), and PRO-Seq libraries were constructed in duplicate as described (63). Uniquely indexed libraries were pooled and sequenced on an Illumina NextSeq instrument using 75 bp single-end reads by the BioFrontiers Sequencing Facility at the University of Colorado Boulder.

PRO-Seq computational analysis
PRO-Seq data were processed using a standardized Nextflow pipeline (https://github.com/Dowell-Lab/Nascent-Flow). A complete pipeline report detailing all software programs and versions utilized and a detailed quality control report including trimming, mapping, coverage, and complexity metrics are included in Supplemental File S1. TDF coverage files output by the pipeline, normalized by reads per million mapped, were visualized using the IGV genome browser (version 2.8.0). FStitch (version 1.0) and Tfit (version 1.0) were used to identify regions with bidirectional transcriptional activity as described (62). Counts were calculated for each sorted BAM file using multiBamCov from the BEDTools suite (version 2.28.0) (64) and RefSeq:National Center for Biotechnology Information Reference Sequences for hg38 downloaded from the University of California Santa Cruz genome browser (May 18, 2018). Genes and long noncoding RNAs were then filtered such that only the isoform with the highest number of reads per annotated length was kept, and DESeq2 (version 1.20.0; Bioconductor release, version 3.7) was used to determine differentially transcribed genes between treatments. Gene clustering was then performed by intersecting differentially transcribed gene lists to form early peak (upregulated in 30 min WSP versus vehicle and downregulated in 120 min versus 30 min WSP), early plateau (upregulated in 30 min WSP versus vehicle and up in 120 min WSP versus vehicle), and late (upregulated in 120 min WSP versus vehicle and upregulated in 120 min versus 30 min WSP) temporal clusters. Functional annotation of differentially regulated genes was performed using DAVID, version 6.8 (38). For bidirectional comparisons, all predicted bidirectional Tfit calls were aggregated using mergeBed (argument −d 60) from BEDTools (version 2.28.0) to generate an annotation file. Counts were then calculated for each sample using multicov (BEDTools, version 2.28.0), and DESeq2 was used to identify differentially transcribed bidirectionals between treatments.

Assay for ATAC-Seq
Beas-2B cells were grown to confluence in 6-well tissue culture dishes and treated with vehicle (PBS) or WSP for 30 or 120 min. Cells were rinsed and scraped in ice-cold PBS, then 50,000 cells from each treatment were pelleted and processed in duplicate for Omni-ATAC-Seq using the protocol developed by Corces et al. (43). Uniquely indexed libraries were pooled and sequenced on an Illumina NextSeq using 37 bp paired-end reads by the BioFrontiers Sequencing Facility at the University of Colorado Boulder.

ATAC-Seq computational analysis
ATAC-Seq reads were trimmed for adapters, minimum length, and minimum quality using the bbduk tool from the BBMap Suite (version 38.73) with arguments "ref=adapters.fa ktrim=r qtrim=10 k = 23 mink=11 hdist=1 maq=10 min-len=20." Quality control was monitored both pretrim and post-trim for all samples using FastQC (version 0.11.8). Trimmed reads were mapped to the human genome (hg38; downloaded from the University of California Santa Cruz genome browser on September 16, 2019, with corresponding hisat2 index files) using hisat2 (version 2.1.0). Resulting SAM files were converted to sorted BAM files using samtools (version 1.9) and to bedGraph coverage format using genomeCoverageBed from the BEDTools suite (version 2.29.2). Read coverage was then normalized to reads per million mapped using a custom python script, and files were converted to TDF format using igvtools (version 2.5.3) for visualization in IGV. MACS2 (version 2.1.4) callpeak with "-SPMR" argument was applied to sorted BAM files for each replicate pair of samples with parameters set to call peaks with log 2 fold change >1 above baseline with q < 1e-5 as significant. Peaks were then assigned to early peak, early plateau, and late temporal clusters based on simple presence/absence criteria within indicated treatment groups, illustrated by Venn diagrams in Figure 4.

Differential transcription factor motif enrichment
TFEA for both PRO-Seq and ATAC-Seq datasets was performed using TFEA (version 1.0; https://github.com/Dowell-Lab/TFEA; (40)). Sets of consensus regions of interest (ROIs) are first defined for all Tfit-called bidirectionals (for PRO-Seq) or all MACS2-called peaks (for ATAC-Seq) using the muMerge function of TFEA. TFEA then calculates read coverage for each replicate within each ROI and applies DESeq2 to determine differentially represented ROIs by treatment condition and assign rank as a function of p value. For each ranked ROI, the FIMO tool (65) from the MEME suite scans the 3-kb sequence surrounding the ROI center (bidirectional origin or ATAC-Seq peak summit) for positions of transcription factor consensus motif matches (with p value cutoff of 10 −5 ), represented by a curated set of position weight matrices (66). Enrichment scores for each motif are then calculated and corrected for sequence content to reduce known biases associated with local GC enrichment, and p values are determined using Z scores. Motif displacement distributions (as in Fig. 3C) are displayed as heat maps where the location of all motif instances within ±1500 bp is indicated relative to the ROI center. The intensity of color is a normalized fraction of motif instances within that base pair range relative to the total number of motifs in the entire 3-kb region.

ChIP-qPCR
Beas-2B cells were grown to confluence in 10-cm tissue culture dishes and treated with vehicle or WSP for 30 or 120 min. Cells were crosslinked by adding 16% methanol-free Demultiplexing rapid effects of wood smoke on transcription formaldehyde to a final concentration of 1% and incubating for 5 min at room temperature, then ChIP-qPCR was then performed as described (67) using 35 cycles of sonication. A mouse monoclonal antibody raised against AHR was purchased from Santa Cruz (sc-133088X) and used at a concentration of 5 μg/sample. Assays were performed in biologic quadruplicate. ChIP-qPCR primer sequences are listed in Table S2.

siRNA-mediated gene knockdown
Beas-2B cells were plated in 6-well tissue culture dishes in antibiotic-free medium. Approximately 24 h later, cells were transfected with ON-TARGETplus siRNA SMARTpools (Horizon/Dharmacon) against human AHR (si-AHR; L-004990-00), RELA (si-RELA; L-003533-00), or a scrambled control (si-Ctrl; D-001810-10), at a final concentration of 25 nM using Lipofectamine RNAiMAX transfection reagent (Life Technologies) as instructed by the manufacturer. Media were replaced with fresh complete 18 h post-transfection, then 24 h later, cells were treated with vehicle or WSP for 2 h for qRT-PCR assays or harvested without additional stimulation for verification of knockdown via Western blotting.

Statistics
Statistical comparisons for qRT-PCR and ChIP-qPCR assays were made by two-tailed t test using the Bonferroni correction when appropriate. These analyses were conducted using statistical software embedded in Open Office (Apache OpenOffice, version 4.1.7; http://www.openoffice.org/ welcome/credits.html).

Data availability
Raw and processed sequencing data have been submitted to the National Center for Biotechnology Information Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/) under accession number GSE167372.
Supporting information-This article contains supporting information.