Early transcriptional effects of inflammatory cytokines reveal highly redundant cytokine networks

The immunologists and computational biologists of the Immunological Genome (ImmGen) Project aim to generate an exhaustive definition of gene expression and regulatory networks of the mouse immune system. Here, we charted the primary transcriptional changes that result, in the major cell types of the immune system, from in vivo exposure to major inflammatory cytokines.

Inflammatory cytokines are pleiotropic signaling molecules secreted by immunocytes in response to external stimuli such as infections or toxins, or internal stimuli during sterile inflammation (Dinarello, 2009; Rock et al., 2010; van Loo and Bertrand, 2023; Murakami et al., 2019). These cytokines mediate defense mechanisms by coordinating communication between various immunocytes and stromal components. Their production and release trigger a cascade of events that promote inflammatory responses through both paracrine and autocrine effects. This leads to initiation of cell differentiation, proliferation, and recruitment, as well as broader organismal changes, such as fever or behavioral effects (Conti et al., 2004; Larson and Dunn, 2001). They also have more discrete effects on phenotypic differentiation of immunocytes, like Th17 cells (Korn et al., 2009). Dysregulation of IL6 signaling is implicated in many chronic inflammatory conditions (e.g., rheumatoid arthritis, inflammatory bowel disease), autoimmune disorders, and cancers, and they are the targets of some of the most effective immunotherapies developed over the past 30 years.

The main inflammatory cytokines include members of the IL1 (Dinarello et al., 1974), IL6 (Hirano et al., 1986), and TNF families (Carswell et al., 1975). Although their actions are mediated by different families of receptors, these cytokines exhibit functional overlaps due to shared signaling mediators (e.g., NF-κB and STAT transducers [Weber et al., 2010; Smale, 2010; Hayden and Ghosh, 2014; Webster and Vucic, 2020; Guo et al., 2024]), as illustrated in Fig. 1 A. Notably, inflammatory cytokines regulate each other’s production and release, resulting in interconnected signaling pathways (Oeckinghaus et al., 2011).

IL1β is initially synthesized as an intracellular inactive precursor and released in response to proinflammatory stimuli (TLR ligands, inflammasome, and caspase activation) (Fields et al., 2019). It activates cells by binding the IL1R1 receptor together with the IL1RAcP co-receptor (Dinarello, 2009). The IL1β:IL1R1:IL1RAcP complex recruits MyD88, initiating a cascade of proinflammatory signals mediated by oligomerized adapter protein complexes, including IRAK4 and TRAF6, and ultimately activates the NF-κB, JNK, and p38/MAPK pathways (Weber et al., 2010). TNFα is also produced as a proprotein, membrane-bound, and only released after cleavage by ADAM17, which is induced by diverse signals like bacterial ligands, IL1β, or cell activation. When TNFα binds its receptor as a homotrimer (Kucka and Wajant, 2021), it recruits the TRADD/TRAF/RIP adaptor complex, which then activates downstream NF-κB and MAPK pathways (Hayden and Ghosh, 2014), as well as programmed cell death (van Loo and Bertrand, 2023).

IL6 represents yet another structural family (Tanaka et al., 2014). A single protein chain, it binds IL6R together with the gp130 signaling partner shared with other cytokines (Tanaka et al., 2014; Rose-John et al., 2023). IL6 signaling is thought to primarily operate via JAK1/2, and TYK2 adaptors that activate STAT transducers, especially STAT3, for direct transcriptional activation. IL6 also elicits parallel activation of NF-κB, MAPK, and PI3K/Akt pathways, generating some overlap with IL1- and TNF-induced signaling (Oeckinghaus et al., 2011; Zegeye et al., 2018; Fan et al., 2013).

While inflammatory cytokines protect the host by activating the immune system, their dysregulation leads to several chronic inflammatory disorders and autoimmune diseases. Therefore, understanding the events elicited by these cytokines across the immune system is of distinct importance. Previous studies have investigated the transcriptional consequences of exposure to inflammatory cytokines in mixed populations or individual cell types (Jura et al., 2008; Defois et al., 2023), but with frequent emphasis on cells in culture or recurrent cell types (Jiang et al., 2021). However, in vivo investigations using comparative approaches to analyze transcriptomic responses across diverse cell types exposed to IL1β, IL6, and TNFα are incomplete.

Here, we purified and profiled immunocytes of the standard “ImmGen 14-cell set” from mice treated systemically with IL1β, IL6, and TNFα. Our findings reveal a large degree of redundancy between these cytokines, with surprising overlap with responses elicited by cytokines of other families, and shared responses between cell types that diverged between more distant myeloid or lymphoid populations.

The study design followed that used in our previous study of transcriptional signatures of γc cytokines. We determined the response across all immunocyte lineages, as represented in the ImmGen 14-cell set, which ranges from mast (MCs) to regulatory T cells (Tregs), after i.v. injection of inflammatory cytokines (eschewing confounders from tissue culture) and leaving cells in their normal environment of inter-cellular signals. In line with our previous determinations in this systematic profiling of cytokine signatures in immunocytes, we continued to use population profiling from carefully sorted cell populations with ImmGen’s “ULI” (ultra-low input) RNA sequencing (RNAseq) protocol, as opposed to single-cell RNAseq. Although it lacks information at the single-cell level, population profiling yields deep signatures, less confounded by the sparsity of single-cell data and by the imbalanced representation of different cell types (and their potentially ambiguous integration on 2D projections).

Inflammatory cytokines, components of innate response pathways, elicit very fast responses (Cronkite and Strutt, 2018). Thus, for the primary phase of the screen, cells were harvested for profiling 2 h after injection, when primary responses should dominate and indirect effects be comparatively limited. Mice were injected with IL1β, IL6, or TNFα at non-lethal doses chosen to match those commonly used in functional experiments. Peritoneal cells and splenocytes were harvested 2 h later, and cell populations were purified by flow double-sorting for ULI RNAseq following standard ImmGen sorts (see Materials and methods for gating strategies and population acronyms). We noted that IL1β and TNFα treatment led to the rapid removal of CD155 (CSF-1R) from the surface of monocytes, and we switched to CX3CR1 as a sorting marker for monocytes (Fig. S1 A). In all, 185 RNAseq datasets passing quality control thresholds were assembled in biological triplicates for all cell/cytokine combinations (except IL1β in plasmacytoid dendritic cells [pDCs]). A dedicated interface of ImmGen website (https://www.immgen.org/Databrowser19/Cytokines.html) presents these results in an interactive manner.

As illustrated in Fig. 1 B by the number of transcripts up- or downregulated at an arbitrary threshold (FoldChange >2 and t test on log-transformed P value <10⁻⁴), the scale of responses to inflammatory cytokines differed between lineages: these cytokines mostly acted upon myeloid cells, particularly neutrophils (GNs), where they had the most extensive effects, affecting 200–2,000 transcripts (for comparison, similarly calculated effects of common-γ-chain (γc) cytokines were in the 50–800 range [Baysoy et al., 2023]). Responses were particularly narrow in B cells (10-fold fewer than in macrophages) and also relatively modest in T lymphocytes. Throughout cell types, IL1β and TNFα had comparatively similar effects, both generally fivefold stronger than IL6.

The clustered heatmap of Fig. 2 A represents 2,331 transcripts altered at P value <10⁻⁴ by one cytokine in at least one cell type (every replicate is shown as a FoldChange versus the mean of the controls, emphasizing the robustness of the data). This synthetic overview of these changes brings forth several major conclusions, which will be further detailed below: (1) There was extremely high similarity between responses to IL1β and TNFα, throughout all cell types. (2) Clusters of affected transcripts were mostly shared between myeloid cells, with quantitative shades, and few truly cell-type-specific effects, with the exception of Clusters 10 and 7, induced by all three cytokines in GNs, and either unaffected (Cluster10) or even repressed (Cluster7) elsewhere. (3) Only one cluster (Cluster6) was predominantly induced in lymphoid cells, which otherwise tended to show changes in the same clusters as myeloid cells, albeit weaker. This similarity was particularly marked for natural killer (NK) cells, reminiscent of notions that NK cells (and innate lymphoid cells more generally) are somewhat an intermediate cell type between adaptive lymphocytes and innate myeloid cells (Bezman et al., 2012). (4) Although one tends to consider responses to cytokines as mainly inductive (via STAT or NF-κB activation), transcript downregulation was extensive even at these early times, numerically equivalent to upregulation, and mirroring the cell/cytokine distribution (a paradox that we had noted earlier for γc cytokines [Baysoy et al., 2023]).

As might be expected, genes induced by inflammatory cytokines include a wide array of cellular functions, as revealed by enrichment analysis. A large number of over-represented pathways were found (544 BioCarta terms enriched at P_adj <10⁻²), with a diversity of functions exemplified in Fig. 2 B. Enrichment in Hallmark pathways brought forth TNFα and IL6 signaling, comfortingly, but also revealed an overlap with signaling pathways of other cytokine families, e.g., interferons (IFNs), IL2, or TGFβ (Fig. 2 C), which we further investigate below. In contrast, genes repressed by inflammatory cytokines showed hardly any enrichment in particular pathways or gene ontologies (Fig. 2 B). One might speculate that, unlike the induction facet which turns up entire functional pathways, this repression is quantitatively important to “clear out space” on ribosomes but does not shut down entire cell functions.

A more discriminating analysis of the different gene clusters defined in Fig. 2 A revealed clear partitions according to their transcriptional control (Fig. 2 D) and functional characteristics. From a regulatory standpoint, the different clusters formed a patchwork of enrichment for very distinct combinations of transcription factor binding motifs. These included several of the expected transcription factors (TFs), NF-κB, MYC, IRF, STAT (Weber et al., 2010; Smale, 2010; Hayden and Ghosh, 2014; Webster and Vucic, 2020; Guo et al., 2024), and CHD1 (Honma et al., 2020; Zhao et al., 2020). MYC stood out most uniquely as driving Cluster4, and its transcript was strongly induced in myeloid cells, marginal zone B cells (MZB), and Treg (but not in GNs). This involvement echoes the large representation of MYC targets in response to γc cytokines (Sarosiek et al., 2010; Sinclair et al., 2013; Villarino et al., 2022; Baysoy et al., 2023). Signaling through NF-κB is expected to represent a major component of the IL1β and TNFα-driven responses (Weber et al., 2010; Smale, 2010; Hayden and Ghosh, 2014; Webster and Vucic, 2020). NF-κB–binding motifs were indeed present, perhaps more narrowly than one might have expected, with the over-representation of NF-κB–binding sites in Clusters 1 and 9. Control of the former, which is almost exclusively myeloid-specific, also involved the canonical myeloid TFs IRF8 and PU.1, whose motifs were less represented in regulatory regions of the more broadly induced Cluster9. Interestingly, IRF8 appears dominant for Cluster6, one of the main response clusters in T cells, but here associated with RUNX and STAT. CHD1 is a chromatin remodeler that directly upregulates Il6 (Zhao et al., 2020), appears involved in regulating several clusters, in particular Cluster11, and it is interesting to note that Il6 belongs to the myeloid-specific Cluster11. Overall, the response to inflammatory cytokines thus seems to involve a combination of TFs.

Gene set enrichment analysis of the upregulated clusters unveiled broad pathways (production of cytokines and chemokines, metabolic processes, and response to stress) but each cluster was also characterized by a functional bias (list and statistics in Table S2). For instance, Cluster1 (264 transcripts) is one of the strongly myeloid-dominant clusters, with divergent effects of IL6 compared with IL1β and TNFα (especially in monocytes). It is enriched in genes associated with innate pattern recognition and includes Tnf itself along with several members of the NF-κB transcription factor family members (Rela, Relb, Nfkb1, and Nfkb2), suggesting a fast positive feedback loop. Cluster4 (174 transcripts) is also myeloid-dominant, yet mostly downregulated by IL1β and TNFα in GNs. This gene set is highly enriched in post-transcriptional processes such as ribosome biogenesis or translational regulation. Cluster9 (151 transcripts), on the other hand, is more ubiquitously induced by IL1β and TNFα and includes many transcripts associated with cytokine signaling (Il4ra, Il15ra, Jak3, Stat3, and Socs3).

Inflammatory cytokines are known to partake in positive feed-forward loops: cytokines induce their own transcripts, which helps amplify fast responses by the innate immune system. It was thus of interest to ask how these feed-forward elements are distributed among immunological lineages. Indeed, TNFα and IL1β reciprocally induced their transcripts in most myeloid cells, except for pDCs, and GNs where Il1β transcripts were extremely abundant even prior to treatment (Fig. 3 A). Feed-forward induction was absent in all lymphoid cell types, even for Tnf whose baseline expression levels were similar in myeloid and lymphoid cells. Il6 transcripts were also induced by IL1β and TNFα in most myeloid cells, but IL6 protein was completely unable to induce the other inflammatory cytokines (Fig. 3 A), indicating a unidirectional spreading of the inflammatory response, at least at this early point. More generally, inflammatory cytokines induce a wide array of other chemokines and other cytokines in myeloid cells but fewer in lymphoid cells, showing the importance of basal states in allowing feed-forward responses to develop (note that the opposite was true for γc cytokines, with their preponderant responses in lymphoid cells).

We also asked whether inflammatory cytokines upregulate the expression of their receptors. This was generally not the case (Fig. 3 B): except for slight induction of Il1r2, genes encoding IL6, TNFα, and IL1 receptors were impervious to these inflammatory cytokines. Thus, by inducing themselves but not their receptors, inflammatory cytokines depart operationally from γc cytokines: IL2 and IL4 do not induce themselves but strongly induce their key receptors, IL2RA and IL4RA (Baysoy et al., 2023).

To determine whether these early transcriptional changes are mirrored at the chromatin level, we analyzed changes in the chromatin histone induced by IL1β using CUT&RUN (Skene and Henikoff, 2017) datasets that we generated in the context of another ImmGen program (https://www.immgen.org/Databrowser19/Chromatin.html). Because of the nature of this dataset, it was not possible to fully match treated and untreated samples from the same organ, so we merged chromatin CUT&RUN data from GNs of several organs from mice left untreated (spleen, bone marrow, and blood) or treated with IL1β 3 h prior to harvest (lung and peritoneal cavity). We reasoned that, while some tissue-related patterns would be expected in the data, tissue-specific effects would be partially averaged out by the merging, and focusing the analysis on the vicinity of IL1β-responsive genes would allow the main traits associated with the response to IL1β to come through. The marks analyzed included H3K4me3 (promoter), H3K4me1 (active enhancer), H3K27ac (active enhancer), H3K27me3 (Polycomb-repressed chromatin), and H3K36me3 (active gene body transcription), and samples from male and female mice were used as replicates.

Normalized signal intensities from GN-specific peak atlas were generated using Epic2 (Stovner and Sætrom, 2019) from merged GN CUT&RUN bam files for each histone target, with peak-gene assignment based on the nearest gene with ChIPseeker (Yu et al., 2015). The comparison revealed marked changes in the vicinity of genes induced by IL1β for H3K4me3 in control and IL1β-stimulated GNs (Fig. 4 A). These findings are in agreement with previous observations (Markovics et al., 2011). Other histone markers showed more limited (H3K27Ac) or no changes (H3K4me1 and H3K27me3). For repressed genes, a more modest decrease in H3K4me3 signals was observed, with essentially no changes in other histones. Genomic plots centered around highly induced genes, such as Cd14 and Socs3, correspondingly showed a very strong increase in H3K4me3 at and downstream from the promoter (Fig. 4 B, left two panels). There was little change in the enhancer-associated marks. Together with the more general scatterplots, these profiles suggest that the early inductive response to IL1β primarily reflects action at the level of the promoter, more than activation of distal enhancers. Changes in the distal enhancer after exposure to cytokines have been reported (Ostuni et al., 2013), albeit in different timeframes, and a detailed comparative analysis over time would be of interest. For repressed genes, while most genes (Atf3; Fig. 4 B) displayed no changes in histone marks, others such as Trim65 (Fig. 4 B, far right) exhibited clear reductions in H3K4me3. These patterns suggest that downregulation by inflammatory cytokines results from post-transcriptional changes (i.e., mRNA degradation), but also some degree of transcriptional turn-down, as we had previously suggested for γc cytokines (Baysoy et al., 2023).

The heatmap of Fig. 2 A indicated a striking degree of similarity between the signatures of TNFα and IL1β. With only mild quantitative variation in their impact on most clusters, IL6 had lesser but also overlapping effects. Heatmap representations can be misleading and mask subtle divergence, so we compared the responses more closely on FoldChange/FoldChange plots (Fig. 5). In peritoneal macrophages (MF.PC), the response to IL1β and TNFα was virtually identical (Pearson R = 0.79), with mostly only quantitative variation in the degree of induction. As detailed in Fig. S1 B, these joint changes included Socs2, Ifn, Il1rn, and Vdr. Only a few genes responded only to TNFα (i.e., F3, Tnfsf9, or Ccl5; Fig. S1 B). TNFα and IL1β also had highly concordant effects in follicular B cells (Bfo) and naive CD4⁺ T cells (T4n) (Fig. 5, top), but more divergence was observed in GNs. A strong overlap between IL1β and TNFα signatures was expected to some extent since they both engage NF-κB signaling pathways.

Responses to IL6 also included components shared with IL1β and TNFα but were generally more divergent (Pearson R versus TNFα ranging from 0.1 to 0.52 in different cells) Transcripts induced by TNFα and IL1β but not IL6 included the NF-κB family (Rel, Relb, Nfkb2, and Nfkbia, Fig. S1 C), consistent with the comparatively limited role of NF-κB in IL6 signaling, as well as Irf1, Cd274, Icosl, or Il12rb2. Gene sets induced by both IL6 and TNFα/IL1β could be found in many cells (e.g., Ccr5, Il18rap, Sgk1, or Igfbp4, Fig. S1 C). Overall, these data indicate that inflammatory cytokines elicit intrinsically similar responses. Indeed, a systematic search for transcripts that could uniquely be induced by only one of these cytokines across the set of cell types tested yielded only 30, 30, and 237 transcripts for IL1β, IL6, and TNFα, respectively.

A strong similarity between the responses of different immunocytes to inflammatory cytokines was obvious from the shared clusters of Fig. 2 A, induced as well as repressed in all myeloid cells, macrophages, DCs, and GNs, albeit with some quantitative variation. The only exception to this widespread sharing of response between myeloid cells was in Cluster7, almost exclusively responsive in GNs (Fig. 2 A). Lymphoid cells shared their own, more restricted, general patterns. This widespread distribution echoed our prior observations after the challenge with γc cytokines, where essentially no cell-type-specific response could be detected at early times (Baysoy et al., 2023). While this manuscript was in preparation, a broad analysis of cytokine responses in vivo by single-cell RNAseq (Cui et al., 2024) concluded that most cytokines induce highly cell-type-specific responses. In particular, it was claimed that IL1β induces distinct gene programs in almost every cell type. As illustrated in Fig. S2, reanalysis of these data also supports our conclusion of similar responses between cell types, the misleading conclusions having resulted from scoring cell specificity by non-parametric intersection of lists of responsive genes called at a fixed statistical threshold, without considering actual FoldChanges. It was important, however, to assess more directly the cell specificity of the responses observed in the present data. For a general overview, we first showed Jaccard indices of the overlap between responses in different cell types (Fig. 6 A), which confirmed the visual impression from Fig. 2 A. In the volcano plots of Fig. 6 B, the signatures of IL1β treatment in GATA6⁺ MF.PC (at FoldChange >2, P value <0.01) were completely biased in monocytes (>90% bias). A strong bias, albeit less extreme, was found for responses to IL1β in NK and effector-memory CD8⁺ T cells (T8em) (Fig. 6 B, right). For a broader perspective, we plotted the responses to IL1β in MF.PC (x-axis) versus several other cell types of the same mice (Fig. 6 C, top row). Responses in other myeloid cells were highly related, while the milder responses in lymphoid cells were less correlated (except for NK cells, which clearly showed some similarity to MF.PC). Conversely, the responses of T4n (Fig. 6 C, bottom row) bore little similarity to those of GNs, bore a little more similarity to those of DC8, even more to B or NK cells, and were almost superimposable to those of Tregs. To generally quantify the degree of cell specificity in these responses, we simply calculated the number of responsive genes shared between any two cell types and those truly specific to one cell type (i.e., not induced by >1.5-fold in any other cell; Table S3). For IL1β, MCs and GNs exhibited 18% and 32% of specific responses, but for all other cell types, this proportion was very low (3–13%).

These large similarities notwithstanding, it is important to point out that there are fine but true differences between responses of closely related cells, and these may be functionally relevant. For instance, the “closeup” on responses to IL1β in T4n and Treg cells (Fig. 6 D) reveals that Dusp4 or Fos are only induced in T4n cells but not in Tregs (Jun is induced in both), while Il10 or Il1rl1 (encodes the ST2 receptor for IL33) only responds in Tregs. The latter is suggestive given the role of ST2 in Treg subsets with proreparative potential.

Thus, inflammatory cytokines induce very few cell-type-specific responses, but there is a degree of fine specificity of high functional relevance that can be missed by agglomerative analyses.

Redundancy between cytokine-elicited responses has long been recognized (Paul, 1989; Ozaki and Leonard, 2002), most immediately because cytokines within a family share receptor or signaling components, as is the case for γc cytokines, or for IL1β and TNFα. Aside from these immediate causes of overlap, activating secondary signal transducers can lead to unexpected similarities, as can the induction of shared negative feedback molecules of the Socs family. For instance, we described the crosstalk between signatures of Type-1 interferon and γc cytokines, likely tied to the activation of STAT1 by the latter. This overlap showed cell specificity, present in most cell types, but missing in NK cells (Baysoy et al., 2023). We thus asked whether similarities could be found between the signatures of inflammatory cytokines and of γc cytokines. Indeed, a simple count of transcripts responsive to IL1β or TNFα that were also responsive, in the same cell types, to each of the γc cytokines revealed an overlap involving dozens to hundreds of transcripts (Fig. 7 A). As a proportion, this similarity was most marked in NK cells treated with IL2 or IL15, involving almost half of the response to inflammatory cytokines (Fig. 7 A, right). A FoldChange comparison revealed a global correlation between the two responses (Fig. 7 B). For IL2, the key trophic cytokine of Treg cells, we found that >50% of its signature genes were shifted to some degree by IL1β or IL6 (Fig. 7 C).

This overlap between γc and inflammatory cytokines included a significant contribution of MYC-dependent transcripts, as illustrated in Fig. 7 D for IL4 and TNFα in monocytes, suggesting that several co-regulated clusters might be shared. This was indeed the case, as Cluster4 (defined in Fig. 2 A) included many of the transcripts of the γc-induced “Cluster5” defined in Baysoy et al. (2023) (Fig. 7 E). Similarly, Cluster6, which is preferentially induced in lymphoid cells, overlapped with γc-induced “Cluster6.” The major downregulated clusters (3, 7, and 8) also coincided with γc-repressed Clusters 11 and 12. Thus, the intersection between responses to inflammatory and γc cytokines does not merely include scattered transcripts and co-regulated gene modules, consistent with the notion that they share signaling pathways. This point was reinforced by Gene Ontology analysis, which showed that the intersection of transcripts upregulated by γc and inflammatory cytokines (all clusters) includes distinct and well-defined functional clusters (Fig. 7 F and Table S1 D). Most notable is a tight group of ribosome biogenesis transcripts, along with many generic changes related to the activation of transcription and translation. Thus, a common result of inflammatory and γc cytokines act to rapidly prepare immunocytes for major biosynthetic changes.

The analysis of IFN-sensitive genes (ISG) after TNFα stimulation with inflammatory cytokines proved paradoxical: it was only in the generally less-responsive lymphoid cells that the core ISG signature was fully induced, and not or much less consistently in myeloid cells (Fig. 8 A). The distribution of ISG induction was essentially similar for IL1β and only weak for IL6. The cell-specific overlap rules out an indirect effect of IFN that would be released after TNFα and suggests that signaling from TNFα and IL1β connects to the main STAT1/STAT2 signaling pathways in lymphocytes, but not in myeloid cells. We previously used network inference to dissect core ISG into clusters that are controlled by different groups of TFs (Mostafavi et al., 2016): a major “Cluster3” controlled by canonical STAT and IRF factors, and several smaller clusters controlled by other regulators (reprinted in Fig. 8 B). In response to TNFα, ISG cluster3 was strongly induced in lymphocytes but not in myeloid cells (exemplified by MF.PC and T4n in Fig. 8 C), Cluster5 responded in both, and other ISG clusters seemed non-responsive anywhere. Thus, the overlap in responses to different cytokines shows baroque intricacy in relation to the diverse signaling modules that are cross-activated in one cell type but not the other.

To focus on the immediate and mostly direct consequences of cytokine action, the data above depicted very early responses to inflammatory cytokines. On the other hand, it was of interest to follow the temporal dynamics of the effects observed at early times and ask whether later-appearing responses were elicited by the administration. We thus conducted a time-course experiment, purifying cells for RNAseq as above, at time points ranging from 2 to 24 h after systemic injection as above. Because running a time course on all cell types and the three cytokines would be prohibitive, we limited this survey to TNFα and splenic monocytes and MF.PC, those cell types with the strongest early responses. The results showed that the response dominated 2–4 h after injection (somewhat earlier for MF.PC than for monocytes), followed by attenuation by 8 h and further diminishment by 24 h (Fig. 9 A; note that the actual numbers are lower than in Fig. 1 B because of lower statistical power from having only two replicates). When these responses were aligned to the genes and clusters of the early response (2 h, ordered per Fig. 2 A), many persisted to the 4 h time point, decreased at 8 h, and finally resolved at 24 h. Indeed, the “amplifying induction” of inflammatory cytokines themselves followed that trajectory (Fig. 9 C). However, broader examination by cumulating transcripts induced at any one time point highlighted clusters of transcripts that persisted longer or that only appeared after 8 h (Fig. 9 D). Secondarily unfolding responses are common in many response systems, but perhaps less expected for inflammatory cytokines, which are thought to represent rapidly deployed responses to innate immune system triggers. These late-responsive transcripts were enriched in RNA metabolism, including several associated with the cell cycle, with enrichment in motifs for E2F factors (enrichR, hypergeometric P_adj <10⁻¹⁰). In terms of cell specificity, the high degree of similarity between early responses in monocytes and MF.PC tended to diffuse somewhat at later times, with an interesting discordance for some cell-cycle associated transcripts, induced in monocytes but repressed in macrophages (Fig. 9 E).

In the context of the systematic determination of cytokine signatures by the ImmGen consortium, this work provides a valuable compendium of the impact of the key inflammatory cytokines across all lineages of the mouse immune system. Some of the results conformed to expectations, such as the predominantly myeloid-centric distribution of the response and the feed-forward loops induced reciprocally by TNFα and IL1β. However, the results also brought several surprising aspects, leading to some reconsideration of the interconnectivity of cytokine networks.

The first surprise was the extreme similarity of early transcriptional responses elicited by TNFα and IL1β. These two cytokines belong to very different structural families, as do their receptors, and are secreted very differently. They are therapeutic targets of blocking antibodies in mostly different disease indications and are generally considered complementary yet different players in inflammatory responses. Yet the target transcripts they elicited were extremely similar, in both myeloid and lymphoid cells (Fig. 2 A and Fig. 5), as were the most strongly responding cell types. Some overlap was expected, as they both signal strongly through the NF-κB and MAPK cascades (Weber et al., 2010; Smale, 2010; Hayden and Ghosh, 2014; Webster and Vucic, 2020; Guo et al., 2024), but superimposition to such a degree was unexpected. Their early transcriptional signatures being so similar implies that the different biologies of TNFα and IL1β result from differences in availability (preferential release in different organismal contexts, or at different times, of inflammation) rather than of effects.

We should mention one caveat here: since TNFα and IL1β induce each other, this cross-induction may have contributed in part to the overlapping signatures. However, the same cross-induction would apply to normal responses in vivo. In parallel experiments, we similarly profiled immunocyte transcriptomes 2 h after LPS injection: those responses were highly similar to those provoked by TNFα and IL1β, Thus, there appears to be a stereotypical response to many innate immune triggers. One might consider TNFα and IL1β as endogenous endotoxins of sorts.

Second, motif enrichment analysis in the regulatory regions of the different response clusters suggested that NF-κB may not be as dominant a signaling conduit as often represented (Weber et al., 2010; Smale, 2010; Hayden and Ghosh, 2014; Webster and Vucic, 2020; Guo et al., 2024). Certainly, NF-κB appeared as the dominant TF controlling the induction of Clusters 1 and 9, each time in association with different TFs. But NF-κB motifs were conspicuously absent from other response clusters, which were instead dominated by other TFs (MYC, IRF8, STAT, and CHD1). MYC, the deep cellular reprogramming that it controls, was likely activated by the MAPK cascade known to be excited by both TNFα and IL1β. MYC activity can be induced by phosphorylation (via ERK) and by stimulated transcription. Indeed, its expression was also strongly induced by IL1β and TNFα. An involvement of PU.1 in myeloid cells is plausible as it has been associated with IL1 signaling (Marecki et al., 2001; Pietras et al., 2016). IRF8 were presumably brought into play by post-translational modification, as it was not appreciably induced by any cytokine. Overall, the early responses to inflammatory cytokines appeared to involve an interlaced patchwork of transcriptional controllers.

Also surprisingly, IL6 did not induce a clear response cluster of its own, instead mostly a qualitatively and quantitatively reduced component of the response to TNFα and IL1β. There is precedent for such similar effects, like the ability of IL6 and IL1β to induce IL17-producing T cells. Because IL6 belongs to another structural family, signals through a gp130-coupled receptor, and has some unique functions in the immune system (e.g., final maturation of B cells into plasma cells), one might have expected distinct response clusters unique to IL6. These were not apparent in any cell type, however. It is possible that biological consequences specific to IL6 map to only a very small set of genes: the sharing of a majority of signature genes does not preclude fine specificity linked to a few outlier targets.

Third, there was very little cell-type specificity in the responses (Figs. 2 and 6). When taken globally, large differences were present between the main branches of the immune system (e.g., between myeloid and lymphoid cells, Fig. 6 C), but responses were very similar within classes of cells (e.g., between CD4⁺ and CD8⁺ T cell lineages, or between macrophages and monocytes). However, even when the overall signatures were quasi-identical, it was also possible to find isolated but functionally significant differences: Il10 was only induced by IL1β in Treg but not in CD4⁺ T cells (Fig. 6 D). Thus, cell specificity needs to be considered both at the level of overall signatures and of the fine-grained detail of unique response.

A final surprise was the extensive overlap between signatures of inflammatory cytokines and those stimulated by other cytokines, in particular γc cytokines. γc cytokines are mostly involved in trophic support and differentiation and belong to different chapters of cytokine textbooks. Yet, almost 50% overlap was observed in some cells when compared with rough thresholds (Fig. 7 A) and were even more aligned when responses were compared in a more sensitive fashion (Fig. 7, B and C). This was particularly paradoxical for IL4, which is generally considered an anti-inflammatory cytokine. Some of the overlap could be ascribed to shared response modules (Fig. 7 E), in particular the MYC-dependent cluster (Fig. 7 F). The latter suggests that both classes of cytokines prepare cells for reprogramming by activating generic cellular functions like transcription and resetting translation with a strong increase in ribosomes. The overlap with IFN signature genes mirrored our previous observations (Baysoy et al., 2023), which we surmise operates via STAT1/2 (IFN transcripts were not induced). But, we remain baffled by the lineage specificity of this overlap, which was mostly restricted to lymphoid cells that were otherwise low responders to inflammatory cytokines. From a functional standpoint, ISGs provide a first line of antiviral defense, so why would inflammatory cytokines induce them in only T or B lymphocytes?

Overall, these results paint a picture of a multiply redundant cytokine network. Rather than a sharper vision in which cytokines have discrete and cell-specific actions with well-defined immunological consequences, cytokines of different classes can activate the same response modules, and those responses operate similarly in different cell types. This core organization of genomic responses to inflammatory cytokines that is observed here will be further tuned by interactions and modifications at the protein level and by differences in the accessibility to cytokines imposed by cellular architecture.

6–8-wk-old C57BL/6J mice obtained from the Jackson Laboratory were housed in specific pathogen–free conditions. IL1β (Cat# 575102; BioLegend, 5 µg), IL6 (Cat# 575702; BioLegend, 5 µg), and TNFα (Cat# 575202; BioLegend, 5 µg) were injected i.v. in 120 µl PBS. Dosage choices were guided by existing literature (Ramadoss et al., 2009; Chai et al., 1996; Putoczki et al., 2013; Libert et al., 1990; Biesmans et al., 2015). For LPS treatment, 1 ng E. coli LPS O55:B5 (cat.#L2880; Sigma Aldrich) was injected i.v. in 120 μl PBS. All animal experimentation was performed under protocols approved by the Harvard Medical School Institutional Animal Care and Use Committee, protocol IS00001257.

Cells from two mice were pooled for each biological replicate, three biological replicates were sorted for cytokine treatments, and five for the PBS control, unless indicated otherwise. Mice were euthanized at the specified time point after injection, followed by cell preparation per ImmGen standard operating procedure (https://www.immgen.org), identified as MC (mast cell, sorted as CD45⁺[CD11b⁻CD11c⁻][CD19⁻CD4⁻CD8⁻]CD117⁺FcEr1a⁺), GN (neutrophil, sorted as TCRB⁻B220⁻Ly6G⁺), MF.PC/MFpc (peritoneal macrophage, sorted as ICAM2⁺F4/80), MF.RP/MFrp (red pulp macrophage, sorted as TCRB⁻MHCIIint⁻F4/80hi), Mo (monocyte, sorted as TCRB⁻CX3CR1⁺Ly6C⁺), DC8 (CD8⁺ DC, sorted as TCRB⁻CD11chiCD8a⁺), pDC (plasmacytoid DC, sorted as TCRB⁻CD11cintB220⁺SiglecH⁺), Bfo (follicular B cell, sorted as CD19⁺TCRB⁻CD93⁻CD23⁺CD43⁻CD5⁻), MZB (marginal zone B cell, sorted as CD19⁺TCRB⁻CD93⁻CD23⁻CD21hi), NK (natural killer cell, sorted as CD19⁻CD8⁻TER119⁻GR1⁻TCRB⁻NK1.1⁺), Tgd (γδ T cell, sorted as CD19⁻CD8⁻TER119⁻GR1⁻TCRB⁻TCRgd⁺), T8em (CD8⁺ effector-memory T cell, sorted as CD19⁻TCRB⁺CD4⁻CD8⁺CD62⁻), T4n (CD4⁺ naïve T cell, sorted as CD19⁻TCRB⁺CD4⁺CD8⁻CD62Lhi), Treg (regulatory T cell, sorted as CD19⁻CD8⁻TER119⁻GR1⁻TCRB⁺CD4⁺CD25hi). Peritoneal cells (MF.PC, MC) were collected by injecting 10 ml medium (Phenol red-free DMEM, 2% FCS, 0.1% azide, and 10 mM HEPES, pH 7.9) peritoneally and subsequent collection. Spleens were homogenized through a 100-µm filter, and red blood cells were lysed in AKC lysing buffer at a ratio of 1 ml per spleen. The resuspended cell pellets were stained as per the standard 14-cell set protocol (https://www.immgen.org), except for CD115, which showed a downregulation in response to IL1β and TNFα treatments and was replaced by CX3CR1 to define monocytes (Fig. S1 A). 1,000 cells were double-sorted on a BD Aria II sorter into 1.5 ml DNA LoBind tubes containing 5 μl TCL lysis buffer (QIAGEN) supplemented to 1% vol/vol β-mercaptoethanol. Sorted samples were kept on ice for 5 min, centrifuged, and quick-frozen on dry ice.

RNAseq was conducted in accordance with the standard ImmGen low-input protocol at the Broad Institute, as detailed previously (Picelli et al., 2013, 2014). In summary, total RNA was purified using a 2.2X RNA-SPRI bead cleanup, and polyadenylated mRNA was selected using an anchored 3′ reverse transcription primer (5′-AAGCAGTGGTATCAACGCAGAGTACT30VN-3′) (IDT). The purified RNA was converted into cDNA through reverse transcription in a limited polymerase chain reaction (PCR) amplification of first-strand cDNA. Then, the first-strand cDNA was tagmented by Tn5-transposon to add adapter sequences, following the Nextera XT DNA Library Preparation Kit (Illumina). Libraries underwent PCR amplification for 18 cycles with Illumina P7 and P5 index adapters, were pooled, and underwent a 0.9X SPRI clean up. Paired-end sequencing was performed on an Illumina NextSeq500 using 2 × 38 bp reads.

FASTQ reads were preprocessed and filtered as previously described (Baysoy et al., 2023). Briefly, reads were aligned to the mouse genome (GRCm38/mm10 primary assembly and collapsed to one transcript per gene based on GENCODE gene annotations vM25; https://www.gencodegenes.org/mouse/release_M25.html) with STAR 2.7.3a (Dobin et al., 2013). Ribosomal RNA gene annotations were removed, and gene-level quantification was calculated by featureCounts from the Subread package. Then, reads were normalized using the standard median ratio method in DESeq2 from Bioconductor (Love et al., 2014) to a final range between 1 and 450,000. Sample RNA transcript integrity (TIN) scores were calculated using RSeQC 2.6.4 (Wang et al., 2016). After normalization, samples with a TIN score below 45, low maximum correlation coefficients (<0.9) between biological replicates, <8,000 mapped genes, and <1 million uniquely mapped reads were excluded. Additionally, genes with <10 aligned reads and a coefficient of variation exceeding 0.3 were excluded from further analysis. The data generated in this study are available in the GEO database under the accession number GSE239946.

To avoid noise created by genes with very low expression, robustly expressed genes were chosen as outlined in Baysoy et al. (2023). Briefly, for all cell types, genes with an expression level >20 in at least one treatment were retained. Some genes yielded intrinsically high variance even in untreated datasets. These “noisy” genes were removed from consideration (on a per-cell-type basis) if their coefficient of variation in that cell was >1, leaving a total of 16,195 genes for further analyses (Table S1, A and B). FoldChange and t test P values (computed on log-transformed data) were determined for each one of these genes in each cytokine treatment relative to PBS controls. For gene–cell pairs that did not pass the expression level or noise criteria, FoldChanges were set to 1 and −log₁₀(P value) to zero.

For robustness, while leveraging the power of confirmatory responses in different cells or different cytokines, signature genes were included in a given cytokine/cell-type pair if (1) (FoldChange >1.8 and t test −log₁₀[P value] >3 [nominal]) or (2) (FoldChange >1.5 and t test −log₁₀[P value] >2) if the same gene showed (FoldChange >3 and −log₁₀[P value] >4) in another cytokine/cell combination (and the inverse for downregulated genes). Note that these signature tables include all genes responsive in at least one of the cell types, but these genes are not necessarily responsive in all cell types (Table S3).

For a comprehensive cytokine comparison, transcripts with a cytokine/PBS t test P value <10⁻⁴, in at least one cytokine/cell-type pair, were selected. Subsequently, 2,331 genes were selected and grouped into 14 K-means clusters using Morpheus (https://software.broadinstitute.org/morpheus) (Table S1 C). Ontology/pathway and TF motif enrichment analyses (hypergeometric tests, Benjamini-Hochberg correction) were performed in STRING (https://string-db.org) (Szklarczyk et al., 2019) and EnrichR (https://maayanlab.cloud/Enrichr/) (Xie et al., 2021) (computed scores and adjusted metrics of variance listed in Table S2).

Overlapping responses in different cells were quantitated by computing the Jaccard index (length[intersect]/length[union]) between signature gene sets at a threshold FoldChange of 1.5 (0.66 for repressed genes).

γc cytokine signature gene sets were retrieved from the ImmGen website (Baysoy et al., 2023). Interferon signature genes were retrieved from the pan-immunocyte IFN responses (Mostafavi et al., 2016). Hallmark MYC genes were extracted from the Molecular Signatures database (“HALLMARK_MYC_TARGETS_V1” and “HALLMARK_MYC_TARGETS_V2” gene sets obtained from https://www.gsea-msigdb.org/) (Liberzon et al., 2015), and NF-κB target genes were retrieved from MotifMap (https://motifmap-rna.ics.uci.edu/) (Xie et al., 2009). Lastly, cell cycle signature genes were retrieved from Dominguez et al. (2016).

For the rough counts table of Fig. 7 A, we simply counted, on a per cell-type basis, the number of transcripts with FoldChanges >2 and nominal P value <0.01 after either IL1β of TNFα treatment (same data as in Fig. 1 B and Fig. 2 A) and also >2 after 2 h in vivo treatment with one of the γc cytokines, with no conditions on P values (data from Baysoy et al. [2023]).

Data from single-cell RNAseq profiling of cytokine responses (IL2 and IL1b) (Cui et al., 2024) were obtained from the authors as pseudo-bulk tables with per-gene reads-per-10K expression values, mean FoldChange, and adjusted P value between cytokine-treated and control cells. These values were used without further modification in Fig. S2, matching with the present data via GeneSymbol annotation. Responsive genes attributed as cell specific in Fig. 2 of Cui et al. (2024) were extracted from Table S4 from that paper.

CUT&RUN data were prepared and generated in the context of a broad ImmGen program to analyze the histone code across all immunological lineages in (B6xCast)F1 mice, treating samples from male and female mice as replicates, following a previously reported protocol (Baysoy et al., 2023). GNs (10,000 per reaction) from untreated or IL1β-treated mice were sorted into 2X Nuclei Extraction (NE) buffer (40 mM HEPES, 20 mM KCl, 0.2% Triton X-100, 40% Glycerol, 2 mM DTT, 1 mM Spermidine, 2X Roche Complete Protease Inhibitor [#11873580001; Millipore Sigma], freshly supplemented to 2X KDAC inhibitor cocktail [2 µM trichostatin A 1 mM sodium butyrate, 1 mM nicotinamide in 70% DMSO]), diluted to a final concentration of 1X. Samples were slow-frozen in isopropanol at −80°C and processed in batch mode for CUT&RUN. Upon thawing, samples were diluted to 10⁵ cells/ml in 1X NE Buffer and a mixture of 10 μl of activated Concanavalin A (ConA) beads, 2 μl of 1:50 (vol/vol) SNAP-CUTANA K-MetStat Panel, and 0.5 μg of primary antibody (rabbit IgG [13-0042; EpiCypher; Lot 20335004-04], H3K4me1 [701763; Thermo Fisher Scientific; Lot 2135869], H3K4me3 [13-0041; EpiCypher; Lot 210760004-01], H3K27me3 [MA5-11198; Thermo Fisher Scientific; Lot VL3152691], H3K36me3 [ab9050; Abcam; Lot GR3386101-1], H3K27ac [8173S; CST; Lot 8], and H3.3 [91191; Active Motif; Lot 25820004]) was added per reaction (10⁴ cells per target reaction) and incubated overnight.

ConA beads were then washed twice in 250 μl Digitonin Buffer (20 mM pH 7.5 HEPES, 150 mM NaCl, 0.5 mM Spermidine, 1X Roche Complete mini, 0.01% digitonin), incubated in 5 μl CUTANA pAG-MNase (20X) in 50 μl of Digitonin Buffer, washed twice again in 250 μl Digitonin Buffer, and resuspended in 50 μl Digitonin Buffer. MNase activity was activated with the addition of 2 mM CaCl₂ per reaction for chromatin digestion for 2 h at 4°C and terminated by adding 33 μl of High Salt Stop Buffer (750 mM NaCl, 26.4 mM EDTA, 5.28 mM EGTA, 66 μg/ml RNase A, and 66 μg/ml Glycogen). 20 pg CUTANA E. coli Spike-in DNA (18-1401; EpiCypher) was added per sample for experimental normalization. CUT&RUN-enriched DNA was released from ConA beads after a 10-min incubation at 37°C and cleaned up using 2:1 (Bead:DNA) ratio of Serapure beads. Libraries were then prepared using a CUTANA CUT&RUN Library Prep Kit (#14-1001; EpiCypher) and sequenced on Illumina NovaSeq 6000 SP (paired-end 2 × 75 bp read). All steps were optimized and performed on Tecan Freedom EVO robotics platforms with gentle rocking for incubation steps and magnetic capture for medium exchange/washing steps.

Fastq data were processed with the following pipeline: reads were adaptor-trimmed using Trim Galore v0.6.6 and aligned to C57BL/6J and Cast reference genome and pseudogenome sequences in parallel (https://csbio.unc.edu/CCstatus/index.py) using Bowtie2 v.2.3.4.3 with parameters -X 1000 -I 10 --very-sensitive-local --no-mixed --no-discordant --phred33. After converting the genomic coordinates of the aligned reads to mm10 coordinates, bam files containing allele-specific and non-allele-specific reads were merged for this study. Unmapped, mitochondrial, and duplicated reads were removed with SAMTools view and Picard MarkDuplicates, and reads overlapping with the ENCODE blacklist were filtered using BEDTools intersect. GN-lineage specific peaks were called using Epic2 (SICER implementation) (Stovner and Sætrom, 2019) with target-specific width and gap parameters (K4me3: -bin 200 -g 2, K4me1 and K27ac: -bin 200 -g 3 and K27me3: -bin 500 -g 10) using merged bam files of appropriate samples. These peak sets were then used to create a final GN-specific peak atlas using IDR. Reads in the peak data matrix were generated using deepTools multiBamSummary and normalized by counts per million. Peak annotation was performed with ChIPseeker (Yu et al., 2015) and GenomicRanges (Lawrence et al., 2013) packages in R. First, the genomic coordinates of the peaks were processed using GenomicRanges and then peaks were annotated with ChIPseeker by mapping their genomic location to the nearest genes in a reference genome (TxDb.Mmusculus.UCSC.mm10.knownGene, the Mus musculus genome assembly, mm10, from the UCSC Genome Browser) by identifying the closest transcriptional start sites to each peak. Genome coverage tracks were generated using deepTools bamCoverage with parameters --binSize 10 --normalizeUsing CPM --ignoreDuplicates --extendReads --smoothLength 50 and visualized using IGV v2.17.4.

Fig. S1 shows the downregulation of CD115 in splenic monocytes in response to IL-1β (A) and FoldChange/FoldChange analysis of inflammatory cytokine-unique genes in MF.PC (B and C) and T4n. Fig. S2 shows a reanalysis of single-cell RNAseq data from Cui et al. (2024) concluding that immune cell types exhibit similar transcriptional responses when treated with IL1β and IL2. Table S1 shows a list of Log2 FoldChange values (A and C) and P values (B) of induced genes in immunocytes in response to IL1β, IL6, and TNFα. Genes overlapping with γc-induced genes (Baysoy et al., 2023) are shown in Table S1 D. Log₂ FoldChange values of TNFα-responsive time-course genes of splenic monocytes and MF.PC are shown in Table S1, E–G. Table S2 lists a complete output of StringDB and EnrichR analyses of each cluster shown in Fig. 1 B. Table S3 shows counts and proportions of genes induced by inflammatory cytokines that are shared between immunocytes.

Results are displayed on the ImmGen website (https://rstats.immgen.org/Skyline/skyline.html?datagroup=ImmGen%20Cytokines%20Inflammatory). Primary data are available on NCBI Gene Expression Omnibus (GSE239946). Custom code is available from the corresponding author upon reasonable request.

We thank Drs J. Kagan, S. Vos, and N. Hacohen for discussion and advice; K. Hattori, A. Ortiz-Lopez, I. Magill, and A. Wood for help with mice and cell sorting; and the Broad Genetics Platform for profiling.

This study was funded by a resource grant from the National Institutes of Health/National Institute of Allergy and Infectious Diseases to the ImmGen consortium (AI072073).

Author contributions: J.J. Lee: Conceptualization, Data curation, Formal analysis, Investigation, Software, Validation, Visualization, Writing - original draft, Writing - review & editing, L. Yang: Formal analysis, Software, Visualization, Writing - review & editing, J.J. Kotzin: Investigation, Methodology, D. Ahimovic: Data curation, Methodology, M.J. Bale: Data curation, Formal analysis, Methodology, Resources, P.A. Nigrovic: Investigation, Resources, Writing - review & editing, S.Z. Josefowicz: Data curation, Methodology, Supervision, D. Mathis: Supervision, C. Benoist: Conceptualization, Formal analysis, Funding acquisition, Project administration, Supervision, Writing - original draft, Writing - review & editing.