Research ArticleSystems Biology

Integrative analysis suggests cell type–specific decoding of NF-κB dynamics

See allHide authors and affiliations

Science Signaling  25 Feb 2020:
Vol. 13, Issue 620, eaax7195
DOI: 10.1126/scisignal.aax7195

Decoding NF-κB dynamics

The nuclear factor κB (NF-κB) subunits RelA and c-Rel are transcriptional regulators that mediate the responses of cells to stimulation of various receptors by different pathogenic stimuli. Martin et al. investigated how the dynamic patterns of RelA and c-Rel signaling might be translated into gene regulation programs appropriate for a given stimulus. The authors integrated single-cell imaging analysis of the trafficking of fluorescently tagged NF-κB subunits in macrophage and fibroblast cell lines exposed to various pathogenic stimuli with data from gene expression studies. This analysis revealed shared and cell type–specific gene clusters, which suggests that different cell types decode NF-κB subunit dynamics differentially to produce appropriate responses to specific stimuli.

Abstract

The complex signaling dynamics of transcription factors can encode both qualitative and quantitative information about the extracellular environment, which increases the information transfer capacity and potentially supports accurate cellular decision-making. An important question is how these signaling dynamics patterns are translated into functionally appropriate gene regulation programs. To address this question for transcription factors of the nuclear factor κB (NF-κB) family, we profiled the single-cell dynamics of two major NF-κB subunits, RelA and c-Rel, induced by a panel of pathogen-derived stimuli in immune and nonimmune cellular contexts. Diverse NF-κB–activating ligands produced different patterns of RelA and c-Rel signaling dynamic features, such as variations in duration or time-integrated activity. Analysis of nascent transcripts delineated putative direct targets of NF-κB as compared to genes controlled by other transcriptional and posttranscriptional mechanisms and showed that the transcription of more than half of the induced genes was tightly linked to specific dynamic features of NF-κB signaling. Fibroblast and macrophage cell lines shared a cluster of such “NF-κB dynamics–decoding” genes, as well as cell type–specific decoding genes. Dissecting the subunit specificity of dynamics-decoding genes suggested that target genes were most often linked to both RelA and c-Rel or to RelA alone. Thus, our analysis reveals the cell type–specific interpretation of pathogenic information through the signaling dynamics of NF-κB.

INTRODUCTION

Innate immune cells can sense diverse pathogen-associated molecular patterns (PAMPs) as ligands binding to a panel of Toll-like receptors (TLRs) found either on the cell surface or in endosomes (1, 2). Whereas the PAMP-TLR interactions have evolved to be specialized for each receptor, the downstream signaling adaptors and transcription factors (TFs) are extensively shared across the PAMPs (Fig. 1A). Because innate immune cells respond to each PAMP with a functionally tailored gene expression program, the shared usage of TFs, such as nuclear factor κB (NF-κB), suggests that biochemical specificity alone cannot adequately explain the basis of PAMP-specific responses. The activity dynamics of TFs can have profound effects in gene regulation and cell fate decisions (3). Therefore, it is possible that the signaling dynamics of TLR-activated TFs may be used to distinguish among the different PAMP stimuli and achieve functional specificity.

Fig. 1 Experimental design for imaging multiple NF-κB subunits in two cell types.

(A) Macrophages (RAW264.7 cells) and fibroblasts (NIH3T3 cells) were treated with an array of TLR ligands representing bacterial, viral, and fungal pathogens. The TLR ligands that produced detectable NF-κB activation in this study, as well as the targeted TLRs, are shown. LPS activates TLR4 (48), Pam3CSK4 activates TLR1/TLR2 heterodimers (49), R848 activates TLR7 (TLR7 and TLR8 in humans) (50, 51), zymosan activates TLR2/TLR6 heterodimers (52) and dectin-1 (53) (not depicted), poly(I:C) activates TLR3 (54), CpG-DNA activates TLR9 (55, 56), flagellin activates TLR5 (57) (not shown), and profilin activates TLR11/TLR12 heterodimers (not shown) (58, 59). (B) Macrophages or fibroblasts stably transfected with mEGFP-RelA or mKOκ–c-Rel were cultured in the same dish, treated with each of the TLR ligands, and then subjected to live-cell imaging. ssRNA, single-stranded RNA; dsRNA, double-stranded RNA; ssDNA, single-stranded DNA.

NF-κB dimers are activated by the TLR-MyD88 (myeloid differentiation marker 88) or TLR-TRIF (TIR domain–containing adaptor protein–inducing interferon-β) signaling axes. Stimulus-dependent degradation of inhibitor of NF-κB proteins (IκBs) induces the release of NF-κB dimers from IκBs and the accumulation of NF-κB in the nucleus for gene regulation. NF-κB heterodimers or homodimers containing RelA (also known as p65) or c-Rel subunits are critical mediators of the rapid transcriptional response in TLR-stimulated cells. These subunits contain a transactivation domain that mediates the recruitment of transcriptional coactivators, such as p300 and CBP (4). Given that NF-κB signaling dynamics (the temporal profile of nuclear abundance) can encode information about the dose and identity of extracellular stimuli (5–8), an important question is whether the signaling dynamics of NF-κB may also be used to decode the upstream PAMP information in terms of gene regulation. Pioneering single-cell imaging studies of RelA coupled with gene expression analysis have established a convincing link between RelA signaling dynamics and the expression of individual target genes (5, 7, 9–17). However, these findings were limited to RelA dynamics induced by a few cytokines or a single TLR ligand, lipopolysaccharide (LPS). Moreover, the exclusive focus on RelA is difficult to justify, because c-Rel is a highly abundant subunit in many immune cells and may exert a more potent functional role than RelA in certain situations (18).

Understanding the overall contribution of NF-κB in the functional decoding of signaling information requires data about the dynamics of all of the relevant NF-κB dimers, in different cell contexts, in response to diverse stimuli. The examination of diverse stimuli is important, especially because there are no precise and rapid pharmacological methods of perturbing specific features (for example, duration and amplitude) of NF-κB signaling dynamics. Widely used perturbation methods, such as small-molecule inhibitors, small interfering RNA (siRNA)–mediated knockdown, and CRISPR-based knockout systems, entail nonspecific effects and long-term compensatory mechanisms that hinder our ability to manipulate specific features of NF-κB dynamics. Therefore, more insight will be gained about the functional importance of NF-κB signaling dynamics from a greater number of naturally occurring distinct dynamic patterns.

Here, we used a comprehensive profiling approach to investigate the link between NF-κB dynamics, as reflected in RelA and c-Rel, and gene expression in two distinct cell types. Macrophages are immune cells specialized in detecting and responding to diverse pathogens and pathogen-derived substances. We contrasted the responses of a macrophage-like cell line with those of a nonimmune cell type, NIH3T3 fibroblasts, in our integrative analysis and identified shared and tissue-specific decoding of signaling dynamics induced by an array of TLR ligands representing bacterial, viral, and fungal pathogens (Fig. 1). The results suggest that different cell types decode TF signaling dynamics differentially to interpret the stimuli and produce cell type–specific gene expression responses.

RESULTS

Both RelA and c-Rel are activated by TLR ligands in RAW264.7 and NIH3T3 cell lines

To monitor the real-time dynamics of NF-κB dimers containing RelA or c-Rel by live-cell imaging, we stably transfected pools of macrophage-like cells (RAW264.7 cells) (herein referred to as macrophages) with plasmids encoding mEGFP-RelA or mKOκ–c-Rel (fig. S1). RelB, p50, and p52 subunits were not considered for imaging analysis because RelB is not robustly activated by canonical stimuli such as TLR ligands in most cell types and the other two subunits form transcriptionally repressive homodimers (4). We also generated stable cell lines of transfected NIH3T3 cells (hereafter referred to as fibroblasts) using the same expression plasmids (fig. S1). We observed a wide range of fusion protein abundance in individual cells and excluded single cells overexpressing the fluorescent fusion protein in our microscopy analysis to avoid effects from unphysiologically high amounts of RelA or c-Rel. The cells stably expressing mEGFP-RelA or mKOκ–c-Rel were treated with TLR ligands that represent widely different types of pathogens, spanning several major TLRs and the downstream intracellular signaling pathways (Fig. 1). Flagellin and profilin were not included in our comprehensive analysis because they did not visibly activate RelA or c-Rel in either cell line. Live-cell, time-lapse imaging was performed to quantify the single-cell time course of nuclear abundance (referred to as signaling dynamics) of RelA and c-Rel induced by each TLR ligand. Custom-written MATLAB code was used for cell tracking and quantification of nuclear and total cellular intensities of each fusion protein.

LPS, the prototypical immune stimulant, induced RelA responses in macrophages, consisting of a high-amplitude initial peak followed by two or more peaks (Fig. 2), as previously reported (7, 19). LPS-induced c-Rel signaling dynamics closely mirrored those of RelA (Fig. 2). Poly(I:C) and CpG-DNA induced lower-amplitude and variable RelA and c-Rel dynamics (Fig. 2). In both cell types, c-Rel signaling dynamics in response to the various TLR ligands were generally similar to those of RelA, with somewhat reduced magnitudes (Fig. 2). Only half of the TLR ligands (LPS, Pam3CSK4, and zymosan) induced appreciable RelA and c-Rel responses in fibroblasts (Fig. 2). The lack of response to poly(I:C), R848, and CpG-DNA in fibroblasts was presumably due to little or no expression of the cognate TLRs (fig. S2). The ligands that activated both macrophages and fibroblasts elicited cell type–specific signaling dynamics. The LPS and Pam3CSK4 responses in macrophages were characterized by marked oscillations, whereas the fibroblast responses largely consisted of a single peak of nuclear activity, which quickly returned to baseline with minor fluctuations thereafter (Fig. 2). Fibroblast responses to zymosan were delayed and heterogeneous, in comparison to those of macrophages. Note that we chose ligand concentrations that were sufficiently high to induce robust responses but were considered nonsaturating for both cell types. A higher concentration of zymosan could not be used because its autofluorescence interferes with signal quantification upon cell binding and internalization.

Fig. 2 TLR ligand–induced NF-κB dynamics in single cells.

Macrophages and fibroblasts stably expressing fluorescent RelA or c-Rel fusion proteins were treated with the indicated TLR ligands and imaged for approximately 12.5 hours. Ratios of nuclear-to-total (nuc:total) cellular fluorescence of RelA and c-Rel were calculated and plotted. Ligand concentrations were carefully chosen from the reference literature, where available. The concentrations of TLR ligands used were as follows: LPS (10 ng/ml), Pam3CSK4 (40 ng/ml), R848 (1 μM; 350 ng/ml), poly(I:C) (20 μg/ml), CpG-DNA (25 μg/ml), zymosan (5 μg/ml), and flagellin and profilin (250 ng/ml). n ≥ 20 cells for each ligand and NF-κB subunit in each cell type. Results were obtained from at least two independent experiments per ligand, NF-κB subunit, and cell type.

To capture quantitative features of NF-κB signaling dynamics, we focused on the following four measures, which essentially define a single-cell RelA or c-Rel response: (i) the delay until the nuclear signal reached a maximum amplitude (time to peak), (ii) the amplitude of the initial peak (peak amplitude), (iii) the duration of the initial peak (duration), and (iv) the time-integrated nuclear occupancy [area under the (nuclear:total ratio) curve (AUC)] (Fig. 3A). We also prepared a complete summary of the quantified results including statistical differences and similarities among the dynamics (table S1). The quantitative analysis showed that some dynamic features were ligand specific or cell type specific (Fig. 3B). First, the time to RelA or c-Rel peaks was ligand specific, with the shortest time observed for Pam3CSK4. Among the three ligands that activated both macrophages and fibroblasts, the rank order of the time to peak was preserved between the two cell types, with a general delay in fibroblasts. However, the rank orders of the other dynamic features across the ligands were not the same between the cell types. For example, the relative strengths of the responses to LPS and Pam3CSK4 were reversed in the two cell types (LPS and Pam3CSK4 resulted in higher values in macrophages and fibroblasts, respectively). Regarding c-Rel, the amplitude, duration, and AUC were slightly weaker than those of RelA, particularly in fibroblasts. Furthermore, each dynamic feature had distinct profiles across the ligands, particularly in macrophages, which suggests that no single feature can hold the information content of NF-κB responses (Fig. 3B). Each feature of NF-κB dynamics may have a differential effect in the transcriptional control of a specific subset of target genes, providing cells with multiple channels to decode the signaling information into functional cell type–specific gene expression responses.

Fig. 3 Quantification of TLR ligand–induced NF-κB dynamics.

(A) Features of NF-κB signaling dynamics that were analyzed. (B) Quantification of the single-cell RelA and c-Rel dynamics induced by the six TLR ligands that induced responses in either cell type (from the data in Fig. 2). The median value of each dynamic feature is depicted with quartile ranges. One-way ANOVA was performed to measure statistical significance among groups, whereas the Mann-Whitney test was performed for pairwise comparisons. Statistical test results can be found in table S1. AUC, area under the curve; A.U., arbitrary units.

TLR ligands induce cell type–specific patterns of gene expression

To explore whether distinct NF-κB dynamics translated into functional genome regulation, we next investigated the transcriptome-wide changes induced by each of the TLR ligands, as assessed by total RNA sequencing (RNA-seq). Parental lines of macrophages and fibroblasts were treated with the six TLR ligands for 6 or 12 hours (Fig. 4A). Because NF-κB is primarily a transcriptional activator, we focused our analysis on 1633 genes whose expression was induced [log2 fold change (log2FC) ≥ 1] by any of the TLR ligands in either cell type. We chose 6 and 12 hours of ligand stimulation as time points that correspond to the outcome of the gene regulatory responses. Although these are not early time points, a substantial proportion (46 to 65%) of the differentially expressed genes were previously reported to be immediate early genes (20, 21). The gene induction by TLR ligands largely reflected the strengths of the RelA and c-Rel responses from the respective ligands. In macrophages, LPS, Pam3CSK4, R848, and zymosan generated the greatest increases in gene expression, whereas poly(I:C) and CpG-DNA induced only moderate changes (Fig. 4A). In fibroblasts, LPS, Pam3CSK4, and zymosan also induced robust increases in gene expression, with overall FCs that were not as large as those observed in macrophages.

Fig. 4 Dynamic features of NF-κB signaling are linked to gene expression.

(A) Heatmap of the 1633 differentially expressed genes [log2 fold change (log2FC) ≥ 1 and false discovery rate (FDR) < 0.05 in at least one treatment] obtained from whole-transcriptome analysis of TLR ligand–induced changes in macrophage and fibroblast gene expression. The log2FC values of the treated cells relative to unstimulated cells are shown. Each row represents one gene. (B) Pearson correlation values were determined for the log2FC value (mean of the 6- and 12-hour treatments) of each gene with respect to the medians of the single-cell NF-κB dynamics in each cell type and for each NF-κB subunit. Clustering was performed on the basis of the correlation values using partitions around medoids. Values beyond the legend scales were set to extremes. The ordering of the genes is the same as that in (A). (C) Correlations for representative genes derived from each of the four clusters in (A) and (B). The genes’ expression level changes relative to the median of the single-cell dynamics features for each of the NF-κB subunits in each cell type are shown. The dynamics feature that most strongly correlated with gene expression for each example gene is depicted. (D) Gene Ontology Biological Process (GOBP) enrichment analysis of dynamics-correlated genes in each of the clusters presented in (A) and (B). The three most statistically significant GOBP terms for each cluster were merged, and the combined GOBP terms are shown. The number of genes with GOBP annotations is indicated in parentheses. Significant terms are marked as dots. Dot size represents the proportion of cluster genes annotated for each term. P values were adjusted using the default Benjamini-Hochberg method.

As expected, the TLR ligands induced numerous NF-κB–regulated immune response genes. Genes whose expression was increased in both cell types included Nfκbia, Tnfaip3, and Ccl2. Other genes, including Il1b and Cd40, were only increased in expression in macrophages, whereas genes such as Cxcl1 and Cxcl5 were only increased in expression in fibroblasts. The expression of Nfkbie, which encodes a feedback inhibitor of NF-κB (IκBε) (22, 23), was induced by TLR ligands to a greater degree in fibroblasts than in macrophages (table S2), despite macrophages exhibiting stronger NF-κB signaling (Figs. 2 and 3). The higher induction of Nfkbie expression and possibly that of other feedback genes in fibroblasts may contribute to the shorter-lived NF-κB dynamics that were observed in fibroblasts, in comparison to the sustained dynamics in macrophages. Our transcriptome-wide profiles show the extent of the cell type specificity in gene induction patterns (Fig. 4A and table S2).

Features of NF-κB dynamics relate to some, but not all, gene clusters

We sought to discover genes potentially involved in the decoding of NF-κB signaling dynamics by integrating the live-cell imaging data with the total RNA-seq data. We calculated the correlations between dynamic features of RelA and c-Rel signaling and gene induction across the TLR ligand treatments (Fig. 4, B and C; see Materials and Methods). If the correlation was strong, then the gene was operationally termed as “decoding” a dynamic feature of the NF-κB subunit. The analysis revealed a large number of genes whose expression was strongly correlated with the RelA and c-Rel dynamics across the TLR ligands (Fig. 4B). Because short response times are generally associated with stronger responses, time to peak was often anticorrelated with gene induction. We prepared a complete list of correlation coefficients and gene expression values (table S2).

Using an unsupervised clustering analysis, we identified four distinct gene clusters based on how the gene expression was correlated to NF-κB dynamic features (Fig. 4B). Genes in cluster A were distinguished by strong correlations with RelA and c-Rel dynamics in both macrophages and fibroblasts (Fig. 4, A and B; “decoding NF-κB dynamics”), with Ccl2 as an example of a duration-decoding gene. This cluster includes many genes previously established to be direct targets of NF-κB, such as Nfκbia, Nfkbie, Tnfaip3, and Myc (fig. S3). Cluster B genes had their expression changes strongly correlated to RelA and c-Rel dynamics only in macrophages (Fig. 4, A and B; “macrophage-specific decoding”), for example, Tnf, Il1b, Cd40, Icam1, Nlrp3, and Marcks. In fibroblasts, cluster B genes either had negligible correlations or were not expressed (41% with mean absolute correlation of <0.6 and 42% with undetectable expression; Fig. 4, A and B). Gene Ontology (GO) analysis showed that only these dynamics-decoding clusters (A and B) were enriched for “NF-κB signaling” genes (Fig. 4D).

The potential relevance of clusters C and D was less clear. Expression of genes in cluster C was mostly anticorrelated with NF-κB dynamics in both cell types (Fig. 4, A to C), suggesting that these genes may be regulated by other TFs that may function in a mutually exclusive manner with respect to NF-κB. For example, despite the inability of CpG-DNA to induce detectable RelA or c-Rel responses in fibroblasts, CpG-DNA produced substantial gene induction in cluster C (Fig. 4A). Cluster D genes exhibited a mixture of heterogeneous correlation patterns. These findings raised the possibility that the linkage of NF-κB dynamics to gene regulation may not have been accurate in this integrative approach, which prompted us to seek alternative methods.

Nascent transcript-based reclustering reveals cell type–specific decoding of NF-κB dynamics

Hypothesizing that the initial integrative analysis was potentially missing information, we wondered whether we could obtain deeper information from our total RNA-seq data. Because nascent transcripts are a more direct readout of active TFs, we surmised that they would produce a clearer relationship between NF-κB dynamics and gene regulation patterns. Furthermore, our bulk RNA-seq libraries had sufficient sequencing depths to enable the detection and quantification of nascent transcripts (using intron-mapped reads) that generally exist at lower abundance (Fig. 5A). Therefore, instead of using the default RNA-seq quantification (of mostly mature transcripts), we were able to quantify the expression of nascent transcripts using reads that map to introns of the mouse reference genome.

Fig. 5 NF-κB dynamics are better correlated to nascent transcription.

(A) The 1633 differentially expressed genes from Fig. 4A (left) were reanalyzed and reclustered on the basis of nascent-level gene expression changes (center) and correlations with NF-κB dynamics (right). The log2FC values of cells treated with the indicated TLR ligands relative to unstimulated cells are shown. Pearson correlation values were determined for the log2FC value (mean of the 6- and 12-hour treatments) of each gene at the nascent level with respect to the medians of single-cell NF-κB dynamics in each cell type and for each NF-κB subunit. Clustering was performed on the basis of the correlation values using partitions around medoids. Values beyond the legend scales were set to extremes. (B) GOBP enrichment analysis of dynamics-correlated genes in each cluster presented in (A). The three most statistically significant GOBP terms for each cluster were merged, and the combined GOBP terms are shown. The number of genes with GOBP annotations is indicated in parentheses. Significant terms are marked as dots. Dot size represents the proportion of cluster genes annotated for each term. P values were adjusted using the default Benjamini-Hochberg method.

Reanalysis of the RNA-seq data and reclustering of the genes based on correlations between NF-κB dynamics and nascent transcripts delineated the genes presumably regulated by transcriptional (clusters E, F, G, and H) versus alternative mechanisms (cluster I) (Fig. 5A). The display of expression changes in mature versus nascent transcripts alongside the gene clusters highlighted that the genes in cluster I did not have a well-defined correlation with NF-κB dynamics because only their mature transcripts, but not their nascent transcripts, were increased in abundance in ligand-treated cells. Some genes within this cluster had undetectable nascent transcripts, and the correlation to NF-κB dynamics could not even be defined (generating blank entries in the heatmap; Fig. 5A). The gene Zfp36, which encodes the posttranscriptional regulator tristetraprolin (TTP), is in cluster I (24, 25). TTP promotes the degradation of its target mRNAs and is regulated by p38 mitogen-activated protein kinase (MAPK) activity (25). This cluster also helps to explain why cluster D was difficult to characterize: 51% of the genes in cluster D were assigned to cluster I (fig. S4). These results suggest that some genes in cluster I may be regulated by alternative mechanisms, including posttranscriptional control, and that analysis of mature transcript abundances may have been inappropriate for correlation with NF-κB dynamics (as was initially done to generate clusters A to D).

Clusters E, F, and H generally recaptured the patterns observed in the previous clusters A, B, and C (Figs. 4 and 5). Cluster E represents a cell type–independent “NF-κB dynamics-decoding” gene cluster, which includes Nfkbia, Nfkbib, Nfkbid, Nfkbie, Nfkbiz, Relb, Nfkb1, Nfkb2, Bcl3, Tnfaip3, Ikbke, Ccl2, Ccl7, and Myc. A macrophage-specific decoding gene cluster was also redefined as cluster F, which included inflammation-related genes, such as Tnf, Nlrp3, Cd40, Tlr1, Tlr6, Ccl5, Il1a, Il1b, Lcn2, and H2-Q6. The macrophage-specific cluster F had genes with either undetectable expression (44%) or negligible correlation (41%) in fibroblasts. The reanalysis also resulted in ‘‘NF-κB signaling’’ being the top enriched GO term among the NF-κB dynamics–correlated genes in clusters E and F, confirming the specificity of these clusters (Fig. 5B). A “no decoding” cluster H, consisting of cell cycle–related genes showing little correlation to NF-κB dynamics, was reduced with fewer constituent genes in comparison to its previous counterpart, cluster C. These genes may be largely influenced by other TLR-activated factors [for example, interferon regulatory factors (IRFs) and activator protein-1 (AP-1)] that do not involve NF-κB.

Note that the new integrative analysis gave a fibroblast-specific gene cluster G, which was missing from the previous analysis. Cluster G contains Cxcl5, Vcam1, Isg20, Ifitm3, Ccl8, Ccl25, Irf1, Irf2, Numb (mistakenly thought to be anticorrelated in the previous analysis), H2-Q4, Brd4, and Ifit3. About half of the fibroblast-specific decoding genes are from cluster D, which was of unknown relevance, but their nascent transcript abundances strongly correlated with the dynamic features of both RelA and c-Rel (Fig. 5A and fig. S4). Unlike the macrophage-specific clusters described so far, the fibroblast specificity of cluster G was, in large part, due to the lack of correlation (rather than lack of expression) in macrophages: Only 25% of the genes were not expressed, and 69% had negligible correlation in macrophages.

The strength of correlation to the dynamic features of NF-κB may be gene specific, and individual genes may depend more on certain features over others for their regulation. For example, some genes may better decode RelA amplitude, whereas others mostly decode c-Rel AUC. Although such tantalizing examples could be found, we examined the genome-wide influence of each dynamic feature. In macrophages, the strongest absolute correlation between NF-κB dynamics and global changes in gene expression was from the duration of nuclear c-Rel signaling (fig. S5, A and B). Correlation to RelA duration was similarly strong, which is consistent with a study of RelA dynamics and target gene expression using single-cell RNA-seq (scRNA-seq) (16). In fibroblasts, duration and time to peak of RelA displayed the highest absolute correlations. Overall, when nascent transcript expression was used for analysis, most TLR-induced genes were linked to NF-κB dynamics: 91 and 64% were correlated to at least one feature of NF-κB dynamics in macrophages and fibroblasts, respectively (absolute correlation of ≥0.8), with RelA showing a broader and stronger effect than c-Rel (fig. S5C). These findings support the functional importance of NF-κB signaling dynamics in transcriptional responses to diverse TLR ligands, particularly in macrophages.

Transcription of dynamics-decoding genes is linked to RelA or to both RelA and c-Rel

Last, having observed a similar overall correlation of transcriptional responses to RelA versus c-Rel dynamics, we were curious to examine more closely whether there were any genes that seemed to track the dynamics of one subunit over the other. Because many NF-κB dynamics–decoding (positively correlated to the first three features and negatively correlated to time to peak) genes were cell type specific (Fig. 5A), we performed this analysis separately for fibroblasts and macrophages with genes from clusters E, F, and G. In each cell type, the largest set was from TLR-induced genes whose transcription was linked to dynamic features of both RelA and c-Rel (Fig. 6). Slightly smaller in numbers were genes showing transcriptional responses linked to RelA but only weakly to c-Rel. For example, the strongest fibroblast gene induction from the zymosan treatment was from the RelA dynamics–decoding cluster (Fig. 6B), which reflects the weak c-Rel responses relative to those of RelA in zymosan-treated fibroblasts (Figs. 2 and 3B). Few genes showed c-Rel–specific correlations. We note that these results should be interpreted with caution, given the unavoidable difference between the fusion protein abundances in the two reporter cell lines. Further validations will be necessary to ascertain the roles of RelA or c-Rel in specific regulation of these genes using rapid-acting tools that overcome the limitations of conventional methods such as knockdown, knockout, or pulsing of stimuli (11, 22, 26, 27). Nevertheless, GO analysis indicated that some functional categories related to LPS responses and NF-κB signaling were differentially enriched in these gene subsets (Fig. 6 and table S3).

Fig. 6 Subunit-specific correlation analysis.

(A and B) For macrophages (A) and fibroblasts (B), genes from clusters E to G were classified on the basis of the average strength of correlation to RelA and c-Rel dynamic features (see Materials and Methods). Display schemes are the same as those used in Fig. 5. Right: GOBP enrichment analysis of dynamics-correlated genes in each gene subset in the heatmap. The three most statistically significant GOBP terms for each set were merged, and the combined GOBP terms are listed in the plot. The number of genes with GOBP annotations is indicated in parentheses. Significant terms are marked as dots. Dot size represents the proportion of cluster genes annotated for each term. P values were adjusted using the default Benjamini-Hochberg method.

DISCUSSION

Our study enabled a side-by-side comparison of an immune cell type and a nonimmune cell type responding to a panel of pathogen-related stimuli and revealed the extent of cell type–specific decoding of complex NF-κB signaling dynamics. Previous studies of NF-κB dynamics in single cells have largely focused on one cell type each (6, 7, 12, 13, 15, 16, 26, 28–30). Integrative analyses of the extensive transcriptomic profiles and NF-κB signaling dynamics in matching conditions suggest that distinct features of NF-κB dynamics are decoded by different target genes in a cell type–specific manner. Furthermore, a previously uncharacterized relationship between NF-κB dynamics and gene regulation emerged when nascent transcripts were considered instead of mature transcripts. The marked difference between information gained from mature versus nascent transcript abundance also serves as a cautionary note for the utility of single-cell RNA-seq methods for studying how TF signaling dynamics are decoded by target genes, because they generally measure mature transcripts.

We characterized c-Rel signaling dynamics by live-cell imaging and found that they closely resembled those of RelA (Figs. 2 and 3). The overall similarity of the responses suggests that RelA and c-Rel activities are likely regulated by the same or comparable signaling mechanisms. However, c-Rel signaling appeared to be somewhat less robust, relative to that of RelA, consistent with the expectation that RelA is the predominant NF-κB subunit in most settings. We note that the RelA and c-Rel fluorescent fusion proteins in this study were expressed from a constitutive promoter, reporting the abundance and localization of the existing pool of NF-κB containing these subunits. Therefore, the signaling dynamics of endogenous RelA and c-Rel, which are themselves induced by NF-κB, may diverge from those observed here, especially in the late phase of responses. Our findings will be useful in dissecting endogenous signaling dynamics, which are a superimposition of existing versus newly synthesized RelA or c-Rel proteins.

The cell type–specific decoding of NF-κB signaling dynamics may occur partly through gene-specific chromatin requirements for productive transcription. For example, some target genes may need the sustained presence of NF-κB at key regulatory sites, thereby decoding the duration of RelA or c-Rel signaling. Others may depend particularly on a high concentration of nuclear NF-κB, thereby decoding the amplitude of RelA or c-Rel. AUC-decoding genes may produce a transcriptional output that is generally proportional to the total NF-κB signaling activity over time. Genes decoding time to peak may diminish their transcription if the NF-κB arrives at their regulatory sites too late, such as when the enhancers are occupied or chromatin accessibility is limited by other early-acting TFs. Such tissue-specific genomic compartmentalization in the decoding of NF-κB signaling dynamics is an interesting but difficult topic to be investigated with improved tools in the future. Such a regulatory architecture is also consistent with the finding that some target genes are co-regulated by NF-κB dynamics (16). Last, it will be important in the future to obtain additional insight from varying the dose of activating ligands for each of the TLRs. In summary, our results highlight the need for more in-depth investigations of the potential effects of specific features of NF-κB dynamics on transcriptional control in both immune and nonimmune cells.

MATERIALS AND METHODS

Reagents

The plasmids mEGFP-N1 (54767), mKOκ (amplified from Chicken Mermaid S188) (53617) (31), RelA cFlag pcDNA3 (20012), and c-Rel cFlag pcDNA3 (20013) (32) were purchased from Addgene. The pSF-EF1α-Ub-Neo plasmid (OG606) was purchased from Oxford Genetics. Purchased TLR ligands included LPS (Enzo Life Sciences), Pam3CSK4 (Tocris, 4633), R848 (Invivogen, tlrl-R848), poly(I:C) (Invivogen, tlrl-picw), CpG-DNA (Invivogen, ODN 2395), zymosan (Sigma, 24250), profilin (Sigma, SRP8050), and flagellin (Sigma, SRP8029). The primary antibodies polyclonal rabbit anti-RelA (SC-372) and polyclonal rabbit anti–c-Rel (SC-71) were purchased from Santa Cruz Biotechnology. Monoclonal rabbit anti-Gapdh (14C10) was purchased from Cell Signaling Technology. Donkey anti-rabbit antibody (GE Healthcare, NA934V) was used as the secondary antibody.

Plasmid generation

Briefly, complementary DNA from plasmids purchased from Addgene was amplified by polymerase chain reaction (PCR) using Phusion polymerase mix [New England BioLabs (NEB), M0532S] and the appropriate primers (table S4). PCR products were gel-purified using a QIAquick Gel Extraction kit (Qiagen, 28115). The constructs and the pSF-EF1α-Ub-Neo plasmid were digested with restriction enzymes from NEB (Eco R1–HF, R3101S; Eco RV–HF, R3195S; Bsi W1, R0553S; Bsr GI–HF, R3575S), where applicable. The plasmid was dephosphorylated with Antarctic phosphatase (NEB, M02809S). Digested PCR products were purified with a QIAquick PCR purification kit (Qiagen, 28104) and ligated to each other and then into the pSF-EF1α-Ub-Neo plasmid using Promega LigaFast Rapid DNA Ligation kit (Fisher Scientific, PR-M8221). Portions of the ligation reactions were used to transform DH5α cells (Thermo Fisher Scientific, 18265017). Plasmids from single antibiotic-resistant clones were purified using Mini Prep kit(s) (Qiagen, 27104) and Maxi Prep kit(s) (Qiagen, 12362). Construction of the plasmids was verified by DNA sequencing.

Cell culture and transfections

NIH3T3 cells (CRL-1658) and RAW264.7 cells (TIB-71) were newly purchased from the American Type Culture Collection. All cell lines were cultured and maintained at 37°C in a 5% CO2/95% air environment in growth medium composed of phenol red–free Dulbecco’s modified Eagle’s medium (Gibco, 21063-029) supplemented with 10% heat-inactivated fetal bovine serum (Gemini Bio-Products, 100-500) and penicillin (100 U/ml), streptomycin (100 μg/ml), and 2 mM l-glutamine (1% PSG) (Gibco, 15140-122). Cells were passaged by brief exposure to 0.25% trypsin-EDTA (Gibco, 25200-056) (NIH3T3) or by gently scraping (RAW264.7). Cells in growth medium lacking 1% PSG were transfected with expression plasmids using Fugene HD (Promega, E2311) according to the manufacturer’s instructions. Pools of cells with stable integration of the EEF1A1 promoter-driven fusion constructs were obtained by selection with G418 sulfate solution (500 μg/ml) (Mirus, MIR5920). Stable expression and function of the fusion proteins were tested through a combination of Western blotting and microscopy.

Western blotting

Cells were lysed in ice-cold radioimmunoprecipitation assay cell lysis buffer (Millipore, 20-188) containing Complete Ultra Mini protease inhibitors (Roche, 05892970001). Lysates were vortexed, centrifuged at 18,000g for 15 min, and homogenized with insulin syringes. Sample protein concentrations were equalized using the Protein Assay Dye Reagent Concentrate (Bio-Rad, 5000006) and heated at 95°C for 5 min in lithium dodecyl sulfate (LDS) sample buffer (Thermo Fisher Scientific, NP0007) containing 10% 2-mercaptoethanol (Gibco, 31350-010). Proteins were separated by SDS–polyacrylamide gel electrophoresis (PAGE) using 4 to 12% NuPAGE bis-tris gels (Thermo Fisher Scientific, NP0322BOX) and Mops buffer (Thermo Fisher Scientific, NP0001). Proteins were transferred using transfer buffer (Thermo Fisher Scientific, NP00061) onto 0.45-μm polyvinylidene difluoride membranes (Millipore, IPVH00010). Membranes were blocked in 5% (w/v) non-fat milk for 30 min and then sequentially incubated with primary and horseradish peroxidase (HRP)–conjugated secondary antibodies. HRP activity was detected using SuperSignal West Pico Chemiluminescent Substrate (Thermo Fisher Scientific, 34580). Where indicated in the figure legends, membranes were stripped for 15 min with Restore Western Blot Stripping Buffer (21059).

Microscopy

Pools of transfected NIH3T3 and RAW264.7 cells stably expressing RelA or c-Rel fusion protein were seeded individually in Ibidi four-well culture inserts (80466) at medium confluence on 35-mm glass-bottomed dishes (MatTek, P35G-1.5-20-C) and cultured for up to 36 hours as described earlier. The Ibidi insert was then removed, and the medium was changed to growth medium containing 2% serum (1:5 dilution in serum-free growth medium). Cells were then placed in a Zeiss LSM880 AxioObserver confocal microscope incubation system and maintained at 37°C in a humidified environment of 5% CO2 for 3 to 4 hours, at which time recording was begun. Immediately before the second frame, cells were treated with 3× concentrations of each TLR ligand, also in growth medium containing 2% serum, to achieve the indicated 1× final concentrations. Time-lapse images were acquired at 5-min intervals for a time course of ~12.5 hours (150 frames) using laser powers of 0.12% (for NIH3T3 cells) and 0.8% (for RAW264.7 cells) and excitation wavelengths of 488 nm (for mEGFP) and 514 nm (for mKOκ). Additional imaging conditions are as follows: a 40×/1.4 numerical aperture (NA) Plan Apochromat oil objective, a fully open pinhole (600 μm), 16-bit, 1024 × 1024, LineSequential Scan Mode, ChS1 channel, 0.52-μs pixel dwell, 0.63-s frame time, 0.52-μs pixel time, 0.03-ms line time, 750 detector (master) gain, 1 digital gain, and 1.0 zoom. Acquired LSM files were exported as 16-bit TIFF files for further analysis.

Image analysis

Using custom MATLAB code that uses the Image Processing toolbox, responding cells in each time lapse field that did not divide, stayed in view and focus for at least 10 hours, had little overlaps with other cells, and exhibited no signs of cellular stress were subject to further analysis. Cellular boundaries were identified automatically by intensity thresholding. Boundaries were manually verified, adjusted, and separated if necessary. Because no nuclear markers were used, to avoid perturbing cell physiology and potentially NF-κB regulation, nuclear boundaries were carefully identified and segmented manually in every frame. This strict analysis precluded analyzing many more cells. Background-subtracted intensity values were used to calculate the ratios of mean nuclear fluorescence intensity to mean total cellular fluorescence intensity (nuc:total). Where necessary, fibroblasts that did not visibly respond to TLR ligands [R848, poly(I:C), and CpG-DNA] were analyzed.

Total RNA-seq

Parental NIH3T3 and RAW264.7 cells were seeded at medium confluence on 24-well glass-bottomed plates (MatTek, P24G-1.5-13-F) and cultured overnight as described earlier. The next day, the medium was changed to growth medium containing 2% serum (1:5 dilution in serum-free growth medium), and the cells were cultured for an additional 3 to 4 hours. As described earlier, cells were then treated with 3× concentrations of each TLR ligand to achieve final 1× concentrations. For each replicate, nontreated cells, as well as cells treated with TLR ligands for 6 and 12 hours, were lysed in RLT buffer and frozen at −80°C. RNA was later purified using QIAshredder kit(s) (Qiagen, 79656) and RNeasy kit(s) (Qiagen, 74136). RNA concentrations were determined using a NanoDrop spectrophotometer (Thermo Fisher Scientific). Determination of RNA integrity and other quality control (QC) was performed by the Center for Cancer Research (CCR) Sequencing Facility at the National Cancer Institute (Frederick, MD). Only high-quality RNA samples [RIN (RNA integrity number) ≥ 9] were used for RNA-seq library preparation with TruSeq Stranded Total RNA kit(s) (Illumina). Three independent biological replicates, from three independently thawed batches of cells, were performed to obtain two QC-passing biological replicates for every sample. In total, 56 RNA-seq libraries were generated [2 cell types × 6 ligands × 2 time points × 2 biological replicates, as well as 8 untreated control samples taken at 0 hour (two from each of the two replicates per cell type)]. The libraries were multiplexed and sequenced on two NextSeq flow cells at the CCR Sequencing Facility. From each sample, we collected from 16.1 to 55.5 million (median = 30.3 million) sequencing reads per library, with 79 to 90% nonduplicate reads and at least 91% reads with 90% bases having the Phred score quality Q30 and above. The percentage of uniquely mapped reads ranged from 72 to 86% (median = 84%).

RNA-seq analysis

The stranded paired-end 75–base pair (bp) reads were mapped to the mouse reference genome [GRCm38, mm10 (33)] with the STAR v2.5.3a mapper (34) using the Gencode M14 primary assembly annotation (35). For the exonic-based mature expression, the reads were counted per gene using RSEM (36). For the intronic-based nascent expression, we first extracted unique intronic regions, as suggested previously (37), with BEDtools v2.25.0 (38). The reads were counted with the featureCounts function of the Subread package v1.4.6-p3 (39) using the following options: -p –primary -B -C -s 2. The count data were then used for differential expression analysis with the edgeR v3.18.1 (40), R v3.4.2. Before the analysis, we filtered out the following genes: (i) lowly expressed genes [<1 FPKM (fragments per kilobase of transcript per million mapped reads) for exonic and <0.5 FPKM for intronic in >60% of samples in each cell type separately]; (ii) short genes (<200 bp); (iii) ribosomal, transporter RNA genes and other small RNA genes; and (iv) genes from unplaced chromosomes. Read counts were normalized for the library size and common, trended, and tagwise negative binomial dispersions by weighted likelihood empirical Bayes method using the calcNormFactors and estimateDisp functions of the edgeR package. The counts were then fitted into gene-wise negative binomial generalized linear models with the glmFit function. Differential expression was determined pairwise for each ligand and time point compared to untreated samples by likelihood ratio tests using the glmLRT function with contrasts accordingly. A detectable batch effect was caused by using multiple flow cells in the sequencing; therefore, this factor was included in our gene-wise negative binomial general linear model fit (~ flow cell + group) before testing for differential expression. For each cell type, the data were analyzed separately, but the P value correction was performed only once for all differential expression results. Among the differentially increased genes, we excluded 334 that were not annotated as protein-coding genes per Gencode M14 gene biotype. A gene was considered to be differentially expressed if |log2FC| ≥ 1 and false discovery rate (FDR) < 0.05.

Correlation analysis

We first filtered genes whose expression was statistically significantly induced by any of the ligands at any time point and by any expression measure (mature or nascent transcripts). We then used the mean of the edgeR-calculated log2FC values from the two time points (6 and 12 hours) for each of the ligands and calculated the Pearson’s correlations between gene expression (measured by log2FC upon stimulation compared to untreated cells) and each of the four dynamics features: the amplitude of the first peak (amplitude), the duration of the first peak (duration), the time-integrated area under the curve (AUC), and time to the first peak (time to peak) for each of the NF-κB subunits (RelA and c-Rel) and each cell type (macrophages and fibroblasts) separately. We used the median of measurable dynamics values across the number of single cells analyzed (table S1). We calculated the correlations using the corAndPvalue function of the WGCNA R package (41) for both mature and nascent expression (measured as log2FC relative to untreated cells). We then performed a cluster analysis with the partition around medoids (PAM) algorithm (42) based on the correlation values calculated between the expression changes and all four dynamics features for both cell types and both NF-κB subunits. For the purpose of clustering, missing values were set to zero and later reverted to missing values before plotting or any other analyses. The identified clusters were then used to generate a heatmap of correlated and anticorrelated genes with the ComplexHeatmap R package (43). The number of clusters was increased until no new patterns were detected. Genes within clusters were ordered by hierarchical clustering. The correlation heatmap was plotted next to the expression heatmap with the same gene order (based on PAM clustering of the correlation values). For the subunit-specific correlation, we considered genes to be correlated if the average of the absolute correlation coefficients was at least 0.7, whereas genes were considered to be uncorrelated if the average of the absolute correlation coefficients was below 0.6. Despite the panel of TLR ligands analyzed, our sample size was not large enough to support stringent statistical filtering. Hence, we did not correct the correlation coefficients for any multiple testing and instead used their absolute value as a measure of confidence.

Gene set enrichment analysis

The Ensembl gene identifiers for each of the clusters of genes were used for the gene set enrichment analysis using the ClusterProfiler R package (44) using all mouse genes (Gencode MV14) as a reference set. The tool uses a hypergeometrical test with the Benjamini-Hochberg correction for multiple testing (adjusted P values) (45). We used GO Biological Process terms (46, 47) (GO.db R package version 3.5.0) and performed compareCluster analysis with enrichGO function. The statistical significance threshold for the adjusted P value was set to 0.05. Only the three most statistically significantly enriched terms were included in all the GO analyses.

Statistical analysis

For the single-cell dynamics, statistically significant differences among groups were calculated using one-way analysis of variance (ANOVA) and Mann-Whitney tests (pairwise comparisons). P values of <0.05 were considered to be statistically significant. For RNA-seq and related analyses, statistical methods are described in the relevant sections of Materials and Methods.

SUPPLEMENTARY MATERIALS

stke.sciencemag.org/cgi/content/full/13/620/eaax7195/DC1

Fig. S1. Expression of NF-κB fusion proteins in stably transfected cells.

Fig. S2. TLR expression in macrophages and fibroblasts.

Fig. S3. TLR ligand–induced expression of NF-κB target genes and correlations with NF-κB dynamics.

Fig. S4. Reclassification of TLR-induced genes based on the correlation of changes in mature versus nascent transcripts to NF-κB dynamics.

Fig. S5. Correlation of TLR ligand–induced transcriptional changes with NF-κB dynamics.

Table S1. Single-cell dynamics characteristics.

Table S2. Complete list of differentially expressed genes and correlations.

Table S3. List of RelA versus c-Rel specificity in correlations.

Table S4. Primers used to generate expression plasmids containing mEGFP-RelA and mKOκ–c-Rel coding sequences.

REFERENCES AND NOTES

  1. ↵
    1. S. Akira,
    2. S. Uematsu,
    3. O. Takeuchi
    , Pathogen recognition and innate immunity. Cell 124, 783–801 (2006).
    OpenUrlCrossRefPubMedWeb of Science
  2. ↵
    1. T. Kawasaki,
    2. T. Kawai
    , Toll-like receptor signaling pathways. Front. Immunol. 5, 461 (2014).
    OpenUrlCrossRefPubMed
  3. ↵
    1. J. E. Purvis,
    2. K. W. Karhohs,
    3. C. Mock,
    4. E. Batchelor,
    5. A. Loewer,
    6. G. Lahav
    , p53 dynamics control cell fate. Science 336, 1440–1444 (2012).
    OpenUrlAbstract/FREE Full Text
  4. ↵
    1. M. S. Hayden,
    2. S. Ghosh
    , NF-κB, the first quarter-century: Remarkable progress and outstanding questions. Genes Dev. 26, 203–234 (2012).
    OpenUrlAbstract/FREE Full Text
  5. ↵
    1. S. Tay,
    2. J. J. Hughey,
    3. T. K. Lee,
    4. T. Lipniacki,
    5. S. R. Quake,
    6. M. W. Covert
    , Single-cell NF-κB dynamics reveal digital activation and analogue information processing. Nature 466, 267–271 (2010).
    OpenUrlCrossRefPubMedWeb of Science
  6. ↵
    1. D. A. Turner,
    2. P. Paszek,
    3. D. J. Woodcock,
    4. D. E. Nelson,
    5. C. A. Horton,
    6. Y. Wang,
    7. D. G. Spiller,
    8. D. A. Rand,
    9. M. R. White,
    10. C. V. Harper
    , Physiological levels of TNFalpha stimulation induce stochastic dynamics of NF-κB responses in single living cells. J. Cell Sci. 123, 2834–2843 (2010).
    OpenUrlAbstract/FREE Full Text
  7. ↵
    1. M.-H. Sung,
    2. N. Li,
    3. Q. Lao,
    4. R. A. Gottschalk,
    5. G. L. Hager,
    6. I. D. C. Fraser
    , Switching of the relative dominance between feedback mechanisms in lipopolysaccharide-induced NF-κB signaling. Sci. Signal. 7, ra6 (2014).
    OpenUrlAbstract/FREE Full Text
  8. ↵
    1. T. K. Lee,
    2. E. M. Denny,
    3. J. C. Sanghvi,
    4. J. E. Gaston,
    5. N. D. Maynard,
    6. J. J. Hughey,
    7. M. W. Covert
    , A noisy paracrine signal determines the cellular NF-κB response to lipopolysaccharide. Sci. Signal. 2, ra65 (2009).
    OpenUrlAbstract/FREE Full Text
  9. ↵
    1. D. E. Nelson,
    2. A. E. C. Ihekwaba,
    3. M. Elliott,
    4. J. R. Johnson,
    5. C. A. Gibney,
    6. B. E. Foreman,
    7. G. Nelson,
    8. V. See,
    9. C. A. Horton,
    10. D. G. Spiller,
    11. S. W. Edwards,
    12. H. P. McDowell,
    13. J. F. Unitt,
    14. E. Sullivan,
    15. R. Grimley,
    16. N. Benson,
    17. D. Broomhead,
    18. D. B. Kell,
    19. M. R. H. White
    , Oscillations in NF-κB signaling control the dynamics of gene expression. Science 306, 704–708 (2004).
    OpenUrlAbstract/FREE Full Text
    1. S. Zambrano,
    2. I. De Toma,
    3. A. Piffer,
    4. M. E. Bianchi,
    5. A. Agresti
    , NF-κB oscillations translate into functionally related patterns of gene expression. eLife 5, e09100 (2016).
    OpenUrl
  11. ↵
    1. R. A. Kellogg,
    2. S. Tay
    , Noise facilitates transcriptional control under dynamic inputs. Cell 160, 381–392 (2015).
    OpenUrlCrossRefPubMedWeb of Science
  12. ↵
    1. M.-H. Sung,
    2. L. Salvatore,
    3. R. De Lorenzi,
    4. A. Indrawan,
    5. M. Pasparakis,
    6. G. L. Hager,
    7. M. E. Bianchi,
    8. A. Agresti
    , Sustained oscillations of NF-κB produce distinct genome scanning and gene expression profiles. PLOS ONE 4, e7163 (2009).
    OpenUrlCrossRefPubMed
  13. ↵
    1. A. Adamson,
    2. C. Boddington,
    3. P. Downton,
    4. W. Rowe,
    5. J. Bagnall,
    6. C. Lam,
    7. A. Maya-Mendoza,
    8. L. Schmidt,
    9. C. V. Harper,
    10. D. G. Spiller,
    11. D. A. Rand,
    12. D. A. Jackson,
    13. M. R. H. White,
    14. P. Paszek
    , Signal transduction controls heterogeneous NF-κB dynamics and target gene expression through cytokine-specific refractory states. Nat. Commun. 7, 12057 (2016).
    OpenUrl
    1. R. E. C. Lee,
    2. S. R. Walker,
    3. K. Savery,
    4. D. A. Frank,
    5. S. Gaudet
    , Fold change of nuclear NF-κB determines TNF-induced transcription in single cells. Mol. Cell 53, 867–879 (2014).
    OpenUrlCrossRefPubMedWeb of Science
  15. ↵
    1. V. C. Wong,
    2. V. L. Bass,
    3. M. E. Bullock,
    4. A. K. Chavali,
    5. R. E. C. Lee,
    6. W. Mothes,
    7. S. Gaudet,
    8. K. Miller-Jensen
    , NF-κB-chromatin interactions drive diverse phenotypes by modulating transcriptional noise. Cell Rep. 22, 585–599 (2018).
    OpenUrl
  16. ↵
    1. K. Lane,
    2. D. Van Valen,
    3. M. M. DeFelice,
    4. D. N. Macklin,
    5. T. Kudo,
    6. A. Jaimovich,
    7. A. Carr,
    8. T. Meyer,
    9. D. Pe'er,
    10. S. C. Boutet,
    11. M. W. Covert
    , Measuring signaling and RNA-seq in the same cell links gene expression to dynamic patterns of NF-κB activation. Cell Syst. 4, 458–469.e5 (2017).
    OpenUrl
  17. ↵
    1. U. Rand,
    2. M. Rinas,
    3. J. Schwerk,
    4. G. Nöhren,
    5. M. Linnes,
    6. A. Kroger,
    7. M. Flossdorf,
    8. K. Kály-Kullai,
    9. H. Hauser,
    10. T. Hofer,
    11. M. Köster
    , Multi-layered stochasticity and paracrine signal propagation shape the type-I interferon response. Mol. Syst. Biol. 8, 584 (2012).
    OpenUrlAbstract/FREE Full Text
  18. ↵
    1. M. N. Shokhirev,
    2. J. Almaden,
    3. J. Davis-Turak,
    4. H. A. Birnbaum,
    5. T. M. Russell,
    6. J. A. Vargas,
    7. A. Hoffmann
    , A multi-scale approach reveals that NF-κB cRel enforces a B-cell decision to divide. Mol. Syst. Biol. 11, 783 (2015).
    OpenUrlAbstract/FREE Full Text
  19. ↵
    1. J. Bagnall,
    2. C. Boddington,
    3. H. England,
    4. R. Brignall,
    5. P. Downton,
    6. Z. Alsoufi,
    7. J. Boyd,
    8. W. Rowe,
    9. A. Bennett,
    10. C. Walker,
    11. A. Adamson,
    12. N. M. X. Patel,
    13. R. O'Cualain,
    14. L. Schmidt,
    15. D. G. Spiller,
    16. D. A. Jackson,
    17. W. Muller,
    18. M. Muldoon,
    19. M. R. H. White,
    20. P. Paszek
    , Quantitative analysis of competitive cytokine signaling predicts tissue thresholds for the propagation of macrophage activation. Sci. Signal. 11, eaaf3998 (2018).
    OpenUrlAbstract/FREE Full Text
  20. ↵
    1. K. S. Oh,
    2. H. Patel,
    3. R. A. Gottschalk,
    4. W. S. Lee,
    5. S. Baek,
    6. I. D. C. Fraser,
    7. G. L. Hager,
    8. M. H. Sung
    , Anti-Inflammatory chromatinscape suggests alternative mechanisms of glucocorticoid receptor action. Immunity 47, 298–309.e5 (2017).
    OpenUrlCrossRef
  21. ↵
    1. D. M. Bhatt,
    2. A. Pandya-Jones,
    3. A. J. Tong,
    4. I. Barozzi,
    5. M. M. Lissner,
    6. G. Natoli,
    7. D. L. Black,
    8. S. T. Smale
    , Transcript dynamics of proinflammatory genes revealed by sequence analysis of subcellular RNA fractions. Cell 150, 279–290 (2012).
    OpenUrlCrossRefPubMedWeb of Science
  22. ↵
    1. A. Hoffmann,
    2. A. Levchenko,
    3. M. L. Scott,
    4. D. Baltimore
    , The IκB-NF-κB signaling module: Temporal control and selective gene activation. Science 298, 1241–1245 (2002).
    OpenUrlAbstract/FREE Full Text
  23. ↵
    1. S. Basak,
    2. A. Hoffmann
    , Crosstalk via the NF-κB signaling system. Cytokine Growth Factor Rev. 19, 187–197 (2008).
    OpenUrlCrossRefPubMedWeb of Science
  24. ↵
    1. S. Hao,
    2. D. Baltimore
    , The stability of mRNA influences the temporal order of the induction of genes encoding inflammatory molecules. Nat. Immunol. 10, 281–288 (2009).
    OpenUrlCrossRefPubMedWeb of Science
  25. ↵
    1. F. Kratochvill,
    2. C. Machacek,
    3. C. Vogl,
    4. F. Ebner,
    5. V. Sedlyarov,
    6. A. R. Gruber,
    7. H. Hartweger,
    8. R. Vielnascher,
    9. M. Karaghiosoff,
    10. T. Rulicke,
    11. M. Müller,
    12. I. Hofacker,
    13. R. Lang,
    14. P. Kovarik
    , Tristetraprolin-driven regulatory circuit controls quality and timing of mRNA decay in inflammation. Mol. Syst. Biol. 7, 560 (2011).
    OpenUrlAbstract/FREE Full Text
  26. ↵
    1. R. E. C. Lee,
    2. M. A. Qasaimeh,
    3. X. Xia,
    4. D. Juncker,
    5. S. Gaudet
    , NF-κB signalling and cell fate decisions in response to a short pulse of tumour necrosis factor. Sci. Rep. 6, 39519 (2016).
    OpenUrl
  27. ↵
    1. L. Ashall,
    2. C. A. Horton,
    3. D. E. Nelson,
    4. P. Paszek,
    5. C. V. Harper,
    6. K. Sillitoe,
    7. S. Ryan,
    8. D. G. Spiller,
    9. J. F. Unitt,
    10. D. S. Broomhead,
    11. D. B. Kell,
    12. D. A. Rand,
    13. V. Sée,
    14. M. R. H. White
    , Pulsatile stimulation determines timing and specificity of NF-κB-dependent transcription. Science 324, 242–246 (2009).
    OpenUrlAbstract/FREE Full Text
  28. ↵
    1. M. V. Gutschow,
    2. J. C. Mason,
    3. K. M. Lane,
    4. I. Maayan,
    5. J. J. Hughey,
    6. B. T. Bajar,
    7. D. N. Amatya,
    8. S. D. Valle,
    9. M. W. Covert
    , Combinatorial processing of bacterial and host-derived innate immune stimuli at the single-cell level. Mol. Biol. Cell 30, 282–292 (2019).
    OpenUrl
    1. R. A. Kellogg,
    2. C. Tian,
    3. M. Etzrodt,
    4. S. Tay
    , Cellular decision making by non-integrative processing of TLR inputs. Cell Rep. 19, 125–135 (2017).
    OpenUrl
  30. ↵
    1. R. A. Kellogg,
    2. C. Tian,
    3. T. Lipniacki,
    4. S. R. Quake,
    5. S. Tay
    , Digital signaling decouples activation probability and population heterogeneity. eLife 4, e08931 (2015).
    OpenUrlCrossRef
  31. ↵
    1. Z. Han,
    2. L. Jin,
    3. J. Platisa,
    4. L. B. Cohen,
    5. B. J. Baker,
    6. V. A. Pieribone
    , Fluorescent protein voltage probes derived from ArcLight that respond to membrane voltage changes with fast kinetics. PLOS ONE 8, e81295 (2013).
    OpenUrlCrossRefPubMed
  32. ↵
    1. S. Sanjabi,
    2. K. J. Williams,
    3. S. Saccani,
    4. L. Zhou,
    5. A. Hoffmann,
    6. G. Ghosh,
    7. S. Gerondakis,
    8. G. Natoli,
    9. S. T. Smale
    , A c-Rel subdomain responsible for enhanced DNA-binding affinity and selective gene activation. Genes Dev. 19, 2138–2151 (2005).
    OpenUrlAbstract/FREE Full Text
  33. ↵
    1. D. M. Church,
    2. V. A. Schneider,
    3. T. Graves,
    4. K. Auger,
    5. F. Cunningham,
    6. N. Bouk,
    7. H. C. Chen,
    8. R. Agarwala,
    9. W. M. McLaren,
    10. G. R. Ritchie,
    11. D. Albracht,
    12. M. Kremitzki,
    13. S. Rock,
    14. H. Kotkiewicz,
    15. C. Kremitzki,
    16. A. Wollam,
    17. L. Trani,
    18. L. Fulton,
    19. R. Fulton,
    20. L. Matthews,
    21. S. Whitehead,
    22. W. Chow,
    23. J. Torrance,
    24. M. Dunn,
    25. G. Harden,
    26. G. Threadgold,
    27. J. Wood,
    28. J. Collins,
    29. P. Heath,
    30. G. Griffiths,
    31. S. Pelan,
    32. D. Grafham,
    33. E. E. Eichler,
    34. G. Weinstock,
    35. E. R. Mardis,
    36. R. K. Wilson,
    37. K. Howe,
    38. P. Flicek,
    39. T. Hubbard
    , Modernizing reference genome assemblies. PLOS Biol. 9, e1001091 (2011).
    OpenUrlCrossRefPubMed
  34. ↵
    1. A. Dobin,
    2. C. A. Davis,
    3. F. Schlesinger,
    4. J. Drenkow,
    5. C. Zaleski,
    6. S. Jha,
    7. P. Batut,
    8. M. Chaisson,
    9. T. R. Gingeras
    , STAR: Ultrafast universal RNA-seq aligner. Bioinformatics 29, 15–21 (2013).
    OpenUrlCrossRefPubMedWeb of Science
  35. ↵
    1. J. M. Mudge,
    2. J. Harrow
    , Creating reference gene annotation for the mouse C57BL6/J genome assembly. Mamm. Genome 26, 366–378 (2015).
    OpenUrlCrossRefPubMed
  36. ↵
    1. B. Li,
    2. C. N. Dewey
    , RSEM: Accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics 12, 323 (2011).
    OpenUrlCrossRefPubMed
  37. ↵
    1. J. G. S. Madsen,
    2. S. F. Schmidt,
    3. B. D. Larsen,
    4. A. Loft,
    5. R. Nielsen,
    6. S. Mandrup
    , iRNA-seq: Computational method for genome-wide assessment of acute transcriptional regulation from total RNA-seq data. Nucleic Acids Res. 43, e40 (2015).
    OpenUrlCrossRefPubMed
  38. ↵
    1. A. R. Quinlan,
    2. I. M. Hall
    , BEDTools: A flexible suite of utilities for comparing genomic features. Bioinformatics 26, 841–842 (2010).
    OpenUrlCrossRefPubMedWeb of Science
  39. ↵
    1. Y. Liao,
    2. G. K. Smyth,
    3. W. Shi
    , featureCounts: An efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 30, 923–930 (2014).
    OpenUrlCrossRefPubMedWeb of Science
  40. ↵
    1. M. D. Robinson,
    2. D. J. McCarthy,
    3. G. K. Smyth
    , edgeR: A Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26, 139–140 (2010).
    OpenUrlCrossRefPubMedWeb of Science
  41. ↵
    1. P. Langfelder,
    2. S. Horvath
    , WGCNA: An R package for weighted correlation network analysis. BMC Bioinformatics 9, 559 (2008).
    OpenUrlCrossRefPubMed
  42. ↵
    1. A. P. Reynolds,
    2. G. Richards,
    3. B. de la Iglesia,
    4. V. J. Rayward-Smith
    , Clustering rules: A comparison of partitioning and hierarchical clustering algorithms. J. Math. Model. Algorithms 5, 475–504 (2006).
    OpenUrlCrossRef
  43. ↵
    1. Z. Gu,
    2. R. Eils,
    3. M. Schlesner
    , Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics 32, 2847–2849 (2016).
    OpenUrlCrossRefPubMed
  44. ↵
    1. G. Yu,
    2. L.-G. Wang,
    3. Y. Han,
    4. Q.-Y. He
    , clusterProfiler: An R package for comparing biological themes among gene clusters. OMICS 16, 284–287 (2012).
    OpenUrlCrossRefPubMedWeb of Science
  45. ↵
    1. Y. Benjamini,
    2. Y. Hochberg
    , Controlling the false discovery rate: A practical and powerful approach to multiple testing. J. R. Stat. Soc. B. Methodol. 57, 289–300 (1995).
    OpenUrl
  46. ↵
    1. M. Ashburner,
    2. C. A. Ball,
    3. J. A. Blake,
    4. D. Botstein,
    5. H. Butler,
    6. J. M. Cherry,
    7. A. P. Davis,
    8. K. Dolinski,
    9. S. S. Dwight,
    10. J. T. Eppig,
    11. M. A. Harris,
    12. D. P. Hill,
    13. L. Issel-Tarver,
    14. A. Kasarskis,
    15. S. Lewis,
    16. J. C. Matese,
    17. J. E. Richardson,
    18. M. Ringwald,
    19. G. M. Rubin,
    20. G. Sherlock
    , Gene ontology: Tool for the unification of biology. The Gene Ontology Consortium. Nat. Genet. 25, 25–29 (2000).
    OpenUrlCrossRefPubMedWeb of Science
  47. ↵
    1. Gene Ontology Consortium
    , Expansion of the Gene Ontology knowledgebase and resources. Nucleic Acids Res. 45, D331–D338 (2017).
    OpenUrlCrossRefPubMed
  48. ↵
    1. A. Poltorak,
    2. X. He,
    3. I. Smirnova,
    4. M. Y. Liu,
    5. C. Van Huffel,
    6. X. Du,
    7. D. Birdwell,
    8. E. Alejos,
    9. M. Silva,
    10. C. Galanos,
    11. M. Freudenberg,
    12. P. Ricciardi-Castagnoli,
    13. B. Layton,
    14. B. Beutler
    , Defective LPS signaling in C3H/HeJ and C57BL/10ScCr mice: Mutations in Tlr4 gene. Science 282, 2085–2088 (1998).
    OpenUrlAbstract/FREE Full Text
  49. ↵
    1. M. S. Jin,
    2. S. E. Kim,
    3. J. Y. Heo,
    4. M. E. Lee,
    5. H. M. Kim,
    6. S. G. Paik,
    7. H. Lee,
    8. J. O. Lee
    , Crystal structure of the TLR1-TLR2 heterodimer induced by binding of a tri-acylated lipopeptide. Cell 130, 1071–1082 (2007).
    OpenUrlCrossRefPubMedWeb of Science
  50. ↵
    1. H. Hemmi,
    2. T. Kaisho,
    3. O. Takeuchi,
    4. S. Sato,
    5. H. Sanjo,
    6. K. Hoshino,
    7. T. Horiuchi,
    8. H. Tomizawa,
    9. K. Takeda,
    10. S. Akira
    , Small anti-viral compounds activate immune cells via the TLR7 MyD88-dependent signaling pathway. Nat. Immunol. 3, 196–200 (2002).
    OpenUrlCrossRefPubMedWeb of Science
  51. ↵
    1. M. Jurk,
    2. F. Heil,
    3. J. Vollmer,
    4. C. Schetter,
    5. A. M. Krieg,
    6. H. Wagner,
    7. G. Lipford,
    8. S. Bauer
    , Human TLR7 or TLR8 independently confer responsiveness to the antiviral compound R-848. Nat. Immunol. 3, 499 (2002).
    OpenUrlCrossRefPubMedWeb of Science
  52. ↵
    1. A. Ozinsky,
    2. D. M. Underhill,
    3. J. D. Fontenot,
    4. A. M. Hajjar,
    5. K. D. Smith,
    6. C. B. Wilson,
    7. L. Schroeder,
    8. A. Aderem
    , The repertoire for pattern recognition of pathogens by the innate immune system is defined by cooperation between toll-like receptors. Proc. Natl. Acad. Sci. U.S.A. 97, 13766–13771 (2000).
    OpenUrlAbstract/FREE Full Text
  53. ↵
    1. B. N. Gantner,
    2. R. M. Simmons,
    3. S. J. Canavera,
    4. S. Akira,
    5. D. M. Underhill
    , Collaborative induction of inflammatory responses by dectin-1 and Toll-like receptor 2. J. Exp. Med. 197, 1107–1117 (2003).
    OpenUrlAbstract/FREE Full Text
  54. ↵
    1. L. Alexopoulou,
    2. A. C. Holt,
    3. R. Medzhitov,
    4. R. A. Flavell
    , Recognition of double-stranded RNA and activation of NF-κB by Toll-like receptor 3. Nature 413, 732–738 (2001).
    OpenUrlCrossRefPubMedWeb of Science
  55. ↵
    1. S. Bauer,
    2. C. J. Kirschning,
    3. H. Häcker,
    4. V. Redecke,
    5. S. Hausmann,
    6. S. Akira,
    7. H. Wagner,
    8. G. B. Lipford
    , Human TLR9 confers responsiveness to bacterial DNA via species-specific CpG motif recognition. Proc. Natl. Acad. Sci. U.S.A. 98, 9237–9242 (2001).
    OpenUrlAbstract/FREE Full Text
  56. ↵
    1. M. Jurk,
    2. B. Schulte,
    3. A. Kritzler,
    4. B. Noll,
    5. E. Uhlmann,
    6. T. Wader,
    7. C. Schetter,
    8. A. M. Krieg,
    9. J. Vollmer
    , C-Class CpG ODN: Sequence requirements and characterization of immunostimulatory activities on mRNA level. Immunobiology 209, 141–154 (2004).
    OpenUrlCrossRefPubMed
  57. ↵
    1. F. Hayashi,
    2. K. D. Smith,
    3. A. Ozinsky,
    4. T. R. Hawn,
    5. E. C. Yi,
    6. D. R. Goodlett,
    7. J. K. Eng,
    8. S. Akira,
    9. D. M. Underhill,
    10. A. Aderem
    , The innate immune response to bacterial flagellin is mediated by Toll-like receptor 5. Nature 410, 1099–1103 (2001).
    OpenUrlCrossRefPubMedWeb of Science
  58. ↵
    1. F. Yarovinsky,
    2. D. Zhang,
    3. J. F. Andersen,
    4. G. L. Bannenberg,
    5. C. N. Serhan,
    6. M. S. Hayden,
    7. S. Hieny,
    8. F. S. Sutterwala,
    9. R. A. Flavell,
    10. S. Ghosh,
    11. A. Sher
    , TLR11 activation of dendritic cells by a protozoan profilin-like protein. Science 308, 1626–1629 (2005).
    OpenUrlAbstract/FREE Full Text
  59. ↵
    1. A. A. Koblansky,
    2. D. Jankovic,
    3. H. Oh,
    4. S. Hieny,
    5. W. Sungnak,
    6. R. Mathur,
    7. M. S. Hayden,
    8. S. Akira,
    9. A. Sher,
    10. S. Ghosh
    , Recognition of profilin by Toll-like receptor 12 is critical for host resistance to Toxoplasma gondii. Immunity 38, 119–130 (2013).
    OpenUrlCrossRefPubMedWeb of Science
Acknowledgments: We thank M. Davidson for providing mEGFP-N1, V. Pieribone for Chicken Mermaid S188, and S. Smale for RelA cFlag pcDNA3 and c-Rel cFlag pcDNA3. We thank the National Cancer Institute CCR Sequencing Facility for RNA-seq QC, library preparation, and sequencing. This work used the computational resources of NIH High Performance Computing Systems, including the Biowulf cluster. Funding: This study was entirely funded by the NIA Intramural Research Program at the NIH. Author contributions: E.W.M. and M.-H.S. designed the study. E.W.M. performed the experiments. A.P. developed the bioinformatics workflow and performed the data analysis. H.P., H.D., and M.-H.S. generated MATLAB scripts for the quantification of imaging data. E.W.M., A.P., and M.-H.S. interpreted the results. E.W.M., A.P., and M.-H.S. wrote the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper or the Supplementary Materials. The materials and data in this study are available upon reasonable request or through appropriate material transfer agreements. Raw RNA-seq data from the 56 samples were deposited at the European Nucleotide Archive (ENA) and are available at www.ebi.ac.uk/ena/data/view/PRJEB24055.

Article Tools

My saved folders

Stay Connected to Science Signaling

Related Content

Similar Articles in:

Citing Articles in:

Navigate This Article