CCR6− group 3 innate lymphoid cells (ILC3s) are mediators of intestinal immunity and barrier function that possess the capacity to acquire type 1 effector features and fully convert into ILC1s. The molecular mechanisms governing such plasticity are undefined. Here, we identified c-Maf as an essential regulator of ILC3 homeostasis and plasticity that limits physiological ILC1 conversion. Phenotypic analysis of effector status in Maf-deficient CCR6− ILC3s, coupled with evaluation of global changes in transcriptome, chromatin accessibility, and transcription factor motif enrichment, revealed that c-Maf enforces ILC3 identity. c-Maf promoted ILC3 accessibility and supported RORγt activity and expression of type 3 effector genes. Conversely, c-Maf antagonized type 1 programming, largely through restraint of T-bet expression and function. Mapping of the dynamic changes in chromatin landscape accompanying CCR6− ILC3 development and ILC1 conversion solidified c-Maf as a gatekeeper of type 1 regulatory transformation and a controller of ILC3 fate.
Introduction
Innate lymphoid cells (ILCs) are tissue-resident immune cells that lack antigen specificity and contribute to tissue immunity, homeostasis, and inflammation (Artis and Spits, 2015; Vivier et al., 2018). ILCs are preprogrammed for effector function during ontogeny, allowing for rapid cytokine secretion upon stimulation. In this regard, ILC progenitors diversify into specialized functional subclasses comprising cytotoxic natural killer (NK) cells and three groups of IL-7R–dependent helper ILCs characterized by the expression of lineage-defining transcription factors (TFs) and cytokines. These include ILC1s defined by T-bet (Tbx21) and IFNγ production; ILC2s defined by high levels of GATA-3 and IL-5 and IL-13 secretion; and ILC3s defined by RORγt and IL-22 production (Artis and Spits, 2015; Vivier et al., 2018).
ILC3s localize predominantly to the intestinal lamina propria; as the major source of IL-22 production (Sawa et al., 2011), they promote epithelial barrier function and protection against extracellular bacteria (Melo-Gonzalez and Hepworth, 2017). While defined by their dependence on RORγt, ILC3s are heterogeneous, encompassing two subclasses distinguished by the expression of CCR6 and thought to derive from distinct progenitors (Constantinides et al., 2014; Klose et al., 2014; Melo-Gonzalez and Hepworth, 2017). CCR6+ ILC3s include fetal lymphoid tissue inducer (LTi) and adult LTi-like cells that are essential for lymphoid tissue formation (Eberl et al., 2004; Melo-Gonzalez and Hepworth, 2017). In contrast, CCR6− ILC3s constitute a remarkably plastic subset characterized by the capacity to express T-bet, IFNγ, and the natural cytotoxicity receptor NKp46, thereby exhibiting a poised type 1 effector state. Within this group, CCR6− NKp46− double-negative (DN) ILC3s are precursors of the NKp46+ ILC3 subset in a developmental pathway that is driven by T-bet induction and dependent on Notch, IL-23, and microbiota signals (Klose et al., 2013; Rankin et al., 2013; Viant et al., 2016). Reminiscent of CD4+ T helper 17 (Th17)→Th1 plasticity, genetic fate-mapping studies demonstrated that NKp46+ ILC3s can further up-regulate T-bet, lose RORγt expression, and fully convert to IFNγ-producing ILC1s, known as ex-ILC3s (Vonarbourg et al., 2010; Artis and Spits, 2015). Notably, human ILC3s undergo similar effector transitions (Bernink et al., 2015; Cella et al., 2019). Although providing flexibility during infection, ILC3-to-ILC1 plasticity is also associated with immune-related pathology (Klose et al., 2013) and is implicated in human autoimmune conditions such as Crohn’s disease (Takayama et al., 2010; Fuchs et al., 2013; Bernink et al., 2015). Thus, plasticity is an inherent feature of physiological CCR6− ILC3 development, function, and dysfunction.
The molecular mechanisms that control ILC3 plasticity are largely unknown. ILC3 identity reflects the integration of both intrinsic TF and extrinsic environmental signals. Indeed, IL-7 and microbiota signals stabilize CCR6− ILC3 identity, whereas IL-12, IL-18, and IL-15 promote conversion to an ILC1 fate (Cella et al., 2010; Vonarbourg et al., 2010; Bernink et al., 2013, 2015). The effector flexibility of NKp46+ ILC3s is facilitated by a unique regulatory state characterized by the balanced activity of RORγt and T-bet lineage-defining TFs that supports a type 3 identity program poised for rapid type 1 functionality and conversion (Vonarbourg et al., 2010; Klose et al., 2013). Interestingly, the type 2 regulator, GATA-3, restrains RORγt expression, allowing for the balance of lineage-defining TFs compatible with development to the NKp46+ ILC3 stage (Zhong et al., 2016). Beyond this, the TF regulators that program the capacity of CCR6− ILC3s to dynamically control their functional fate in response to shifting contexts are largely undefined.
The AP-1 TF, c-Maf, is a pleiotropic regulator of effector programming in adaptive and innate-like lymphocytes. Indeed, c-Maf is essential for activation or repression of key cytokine loci (e.g., Il4, Il21, Il10, and Il22) in various functional subsets of CD4+ T cells (Ho et al., 1996; Bauquet et al., 2009; Rutz et al., 2011; Ciofani et al., 2012) and invariant NKT cells (Yu et al., 2017). c-Maf regulates antiinflammatory and tissue-residency programs in IL-17–producing CD4+ Th17 cells (Ciofani et al., 2012; Aschenbrenner et al., 2018) and is required for the adoption of specialized effector phenotypes by regulatory T cells (Wheaton et al., 2017; Xu et al., 2018) and for type 17 programming of γδ T cells (Zuberbuehler et al., 2019). Transcriptomic profiling of ILC subsets identified c-Maf as most highly expressed in ILC3 cells (Gury-BenAri et al., 2016), suggesting a function for c-Maf in this lineage. Moreover, a recent ILC network analysis revealed c-Maf as a regulator of ILC3/ILC1 balance in the intestine (Pokrovskii et al., 2019). Here, we define c-Maf as an essential regulator of CCR6− ILC3 homeostasis and plasticity that controls the balance of DN and NKp46+ ILC3s and limits physiological conversion to an ILC1 state. c-Maf directly restrained T-bet expression, which is essential for proper CCR6− ILC3 development and lineage stability. Accordingly, c-Maf constrained expression of the type 1 program, including high levels of IFNγ production, while conversely supporting RORγt activity and the expression of key type 3 effectors (e.g., Il22). Underlying global ILC3 transcriptional programming, c-Maf is required to both promote ILC3 chromatin accessibility and restrain the acquisition of an ILC1-like regulatory landscape characterized by T-bet, Runx3, and TCF-1 function. Taken together, our findings define c-Maf as a gatekeeper of type 1 effector acquisition that regulates plasticity and stabilizes NKp46+ ILC3 cell identity.
Results
c-Maf is highly expressed in T-bet+ ILC3
Within the ILC lineage, Maf expression is first induced during ILC specification at the early innate lymphoid progenitor stage immediately downstream of the common lymphoid progenitor (Harly et al., 2018). To explore the role of c-Maf in the ILC3 lineage, we evaluated the expression of c-Maf within ILC subsets located in the small intestine lamina propria (SILP), a tissue where all ILC lineages are present. Here, ILC3s express higher levels of c-Maf protein compared with ILC1s or ILC2s (Fig. 1 A), indicating that c-Maf is dynamically regulated during ILC specialization. c-Maf expression also varies within the ILC3 compartment; it is expressed most weakly in LTi-like CCR6+ ILC3s, at intermediate levels in NKp46− CCR6− DN ILC3s, and at the highest levels in NKp46+ ILC3s (Fig. 1 A). Interestingly, c-Maf expression in ILC3 cells was highly correlated with that of T-bet, with high levels of c-Maf restricted to T-bet+ DN and NKp46+ ILC3 subsets. This pattern suggested that c-Maf may be a relevant ILC3 regulator.
c-Maf is required for ILC3 homeostasis
To study the role of c-Maf in ILCs, we bred mice harboring a conditional Maf allele with mice expressing Cre recombinase from the Il7r locus (Il7rCre) for specific deletion in the lymphoid lineage (Schlenner et al., 2010; Wende et al., 2012). This deletes Maf in the common lymphoid progenitor cells that give rise to all ILC progenitors. Maffl/flIl7rCre mice were healthy, displaying no overt signs of inflammatory disease. Although c-Maf expression is induced in the earliest ILC precursors (Harly et al., 2018), the total number and proportion of lineage marker–negative (Lin−) IL-7Rα+ (also known as CD127) ILCs among lymphocytes was unaltered in the SILP of Maffl/fl Il7rCre compared with control Maf+/+Il7rCre mice (Fig. S1 A), indicating that c-Maf is not obligatory for ILC development. However, the relative proportion and number of individual ILC subsets were significantly affected in the absence of c-Maf (Figs. 1 B and S1 A). Specifically, ILC3 lineage cells were modestly decreased, while ILC2s were concomitantly increased in the Maffl/fl Il7rCre SILP (Figs. 1 B and S1 A), implying a selective requirement for c-Maf in the ILC3 lineage.
We found that ILC3 subset homeostasis was substantially altered in the SILP of Maffl/flIl7rCre mice. This was marked by a significant increase in the proportion of NKp46+ ILC3 and LTi-like CCR6+ ILC3 subsets, and a striking loss in the proportion and number of DN ILC3 cells (Figs. 1 C and S1 A). As DN cells include precursors for NKp46+ ILC3s (Klose et al., 2013; Rankin et al., 2013), this dysregulation suggested enhanced differentiation of DN cells into the NKp46+ ILC3 compartment in the absence of c-Maf. Accordingly, there was a 3.5-fold reduction in the percentage of precursor T-bet− DN ILC3s and a near-absence of transitional T-bet+ DN ILC3s in Maffl/flIl7rCre versus Maf+/+Il7rCre (Fig. 1 C), suggesting that in the absence of c-Maf, ILC3s transit to the NKp46+ stage immediately upon T-bet expression. Thus, c-Maf is critical in maintaining ILC3 subset balance.
To evaluate the role of c-Maf within the ILC3 lineage, we conditionally deleted Maf following the induction of ILC3 differentiation using the Rorc-Cre transgenic strain (Eberl and Littman, 2004). Despite lacking obvious signs of intestinal inflammatory disease, isolation of lymphocytes from the SILP of cohoused Maf+/+ and Maffl/flRorc-Cre mice was inconsistent, potentially due to dysbiosis; therefore colon lamina propria (CoLP) ILCs were primarily analyzed. There were similar numbers of ILC3s in the CoLP of Maffl/flRorc-Cre and control mice (Fig. 1 D). Notably, specific loss of c-Maf in ILC3s resulted in a selective significant increase in the number of ILC1s in Maffl/flRorc-Cre CoLP (Fig. 1 D), suggesting that c-Maf expression in ILC3s regulates the ILC1/ILC3 balance. ILC3s from Maffl/flRorc-Cre also displayed developmental dysregulation, characterized by a severely diminished transitional T-bet+ DN ILC3 population, and significantly elevated NKp46+ ILC3 cells compared with Maf+/+Rorc-Cre control mice (Figs. 1 E and S1 B). For comparison to the Il7rCre model, we observed a similar alteration in ILC3 subset proportions in the SILP of Maffl/flRorc-Cre mice that we were able to analyze (Fig. S1 C). These findings confirm the important role for c-Maf in regulating ILC3 homeostasis.
c-Maf constrains type 1 features and stabilizes type 3 identity in ILC3s
NKp46+ ILC3s possess a unique regulatory state characterized by expression of both lineage-defining TFs, T-bet and RORγt (Sciumé et al., 2012; Klose et al., 2013; Rankin et al., 2013); however, the factors that limit execution of the type 1 program in this type 3 cell are unknown. CCR6− ILC3s exhibit a low-to-high T-bet expression gradient that regulates the transition of DN ILC3s to the NKp46+ stage and subsequent conversion into ILC1s (Klose et al., 2013). Notably, this T-bet expression gradient is disrupted in CCR6− ILC3s from Maffl/flIl7rCre mice, with cells displaying uniformly high levels of T-bet compared with Maf+/+Il7rCre controls (Fig. 2 A). In the absence of c-Maf, T-bet expression is significantly elevated in NKp46+ ILC3s to the maximal level observed in ILC1 cells from control mice, indicating that c-Maf is required to limit T-bet expression in CCR6− ILC3s. Deletion of Maf also resulted in a low-level de-repression of T-bet in CCR6+ ILC3 cells that do not otherwise express T-bet (Fig. 2 A), revealing c-Maf–dependent repression of Tbx21 in type 3 ILCs. In contrast, we observed a modest, but reproducible, decrease in RORγt expression in CCR6− and NKp46+ ILC3 populations but not CCR6+ ILC3 cells from Maffl/flIl7rCre compared with Maf+/+Il7rCre SILP (Fig. 2 A). In light of the antagonistic relationship between T-bet and RORγt in ILC3s (Klose et al., 2013; Zhong et al., 2016), elevated T-bet levels may contribute to the selective reduction of RORγt expression in T-bet+ Maf-deficient ILC3 subsets. While small, similar differences in RORγt expression have been reported to affect ILC3 homeostasis (Zhong et al., 2016). Thus, c-Maf is required for the T-bet expression gradient that defines CCR6− ILC3 subsets and to regulate the balance of lineage-defining TF expression in NKp46+ ILC3s.
We next evaluated effector cytokine expression in Maf-deficient ILC3s. In line with elevated T-bet expression, we observed a striking and selective de-repression of IFNγ production in CCR6− ILC3s from Maffl/flIl7rCre mice compared with controls when stimulated with PMA/ionomycin (Fig. 2 B), and similarly with type 3 conditions (IL-23 and IL-1β; Fig. S2 A). A significant elevation of IFNγ expression was also observed when Maf was deleted in ILC3 cells using Rorc-Cre (Fig. S2 B). Notably, despite T-bet expression, control ILC3s produced minimal IFNγ at steady state, whereas Maf deficiency released the potential for high-level IFNγ among T-bet+ ILC3 cells (Figs. 2 B and S2 B). Intestinal ILC3s are otherwise a major source of homeostatic IL-22 (Satoh-Takayama et al., 2008; Luci et al., 2009; Sawa et al., 2010). While overall IL-22 production from total SILP ILC3s was unaltered in Maffl/flIl7rCre versus Maf+/+Il7rCre mice, there was a moderate but significant 24% decrease in the proportion of IL-22+ cells among CCR6− ILC3s, and interestingly, a reciprocal increase in IL-22 percentage for CCR6+ ILC3s in the absence of c-Maf (Fig. 2 B). Thus, c-Maf functions to constrain type 1 effector function and stabilize type 3 identity in T-bet–expressing ILC3 cells.
Next, we established the timing of c-Maf–dependent effects on ILC3 identity during the initial waves of ILC3 development. Before birth, from embryonic days (E) 15–18, there are few T-betlo and little to no NKp46+ ILC3 cells detected in the fetal intestine (Figs. 2 C and S2 C), as both NKp46 and T-bet expression are up-regulated in ILC3s in the first weeks of life, after microbial colonization (Sawa et al., 2010, 2011; Klose et al., 2013). Notably, in the E15 Maffl/flRorc-Cre fetal intestine, the proportion and level of T-bet expression were precociously elevated in Maf-deficient CCR6− ILC3s compared with control cells (Figs. 2 C and S2 D). This indicates that c-Maf also restrains T-bet expression in CCR6− ILC3s before the availability of T-bet–inductive signals and limits development of progenitor DN ILC3 cells into the T-bet+ ILC3 compartment. Of note, CCR6+ ILC3s are reciprocally reduced in frequency in the E15 Maffl/flRorc-Cre fetal intestine (Fig. 2 C), suggesting a stage-specific role for c-Maf in fetal LTi cells. Accordingly, we observed smaller and disorganized Peyer’s patches in adult mice in which c-Maf was deleted with either Rorc-Cre or Il7rCre (Fig. S2 E and data not shown).
The fetal analysis also reveals that T-bet repression by c-Maf occurs in a cell-intrinsic manner, as other potentially c-Maf–expressing RORγt+ CD3+ cells are absent from the E15 fetal intestine (Fig. S2 D). This is further supported by the observation that features of c-Maf deficiency are recapitulated in Maffl/flIl7rCre mice bred onto a Rag1−/− background lacking T and B cells, including a significant up-regulation of T-bet and IFNγ and down-regulation of RORγt and IL-22 production in CCR6− T-bet+ ILC3s (Fig. 2 D). Moreover, mixed bone marrow chimeras of CD45.2 Maf+/+Rorc-Cre or Maffl/flRorc-Cre with CD45.2+CD45.1+ Maf+/+Rorc-Cre cells similarly demonstrated a cell-intrinsic role for c-Maf in regulating CCR6− ILC3 homeostasis and limiting T-bet expression (Figs. 2 E and S2 F). Thus, c-Maf restrains type 1 features in steady-state ILC3s independently of adaptive immune cells and potential differences in microbiota.
c-Maf globally maintains type 3 and restrains type 1 features
To uncover c-Maf–dependent genes regulating ILC3 identity, we performed RNA-sequencing (RNA-seq) and differential expression analysis for NKp46+ ILC3 cells isolated from the SILP of Maffl/flIl7rCre and Maf+/+Il7rCre mice. For comparison, we determined c-Maf–dependent expression profiles for SILP ILC2 and CCR6+ ILC3 cells. Importantly, both mutant and control ILC3 subsets purified based on cell surface phenotype were uniformly RORγt+ (Fig. S3 A). While NKp46+ ILC3s shared some c-Maf–dependent regulation with CCR6+ ILC3 or ILC2 cells, a relatively large fraction of c-Maf–dependent genes (∼66%) were unique to each subset, indicating that c-Maf plays a distinct role in NKp46+ ILC3 cells (Fig. S3 B).
In the absence of c-Maf, a similar number of genes were significantly up- and down-regulated in NKp46+ ILC3 cells (Fig. 3 A). Interestingly, several down-regulated c-Maf–dependent genes were shared with CCR6+ ILC3 cells (Fig. S3 C), IL-17–producing γδ T cells (Tγδ17; Fig. S3 D), and RORγt+ regulatory T cells (Wheaton et al., 2017; e.g., Blk, Vdr, Ltb4r1, Icos, Rnase4, and Matn2), suggestive of canonical c-Maf targets. Consistent with the phenotype of c-Maf–deficient ILC3s, numerous type 3 effector genes were also down-regulated, including Rorc, Il22, Il23r, and Il1r1 (Fig. 3, A and B). Conversely, up-regulated genes displayed a striking abundance of loci encoding both type 1 effectors (e.g., Tbx21, Ifng, Bcl6, and Il21r), and NK cell–related receptors and factors (e.g., Gzmc, Fasl; Fig. 3, A and B). Among these, we confirmed that NK1.1 (Klrb1c) and NKG2A/C/E (Klrc1, Klrc2, and Klrc3) protein expression was substantially elevated on c-Maf–deficient NKp46+ ILC3s, with levels intermediate between control NKp46+ ILC3s and ILC1s (Fig. 3 C). Notably, Tbx21 was also significantly differentially expressed (DE) in Maf-deficient CCR6+ ILC3s; however, the resulting transcript is minimally expressed, with levels 23-fold lower than that of NKp46+ ILC3s (Fig. S3 C). In line with these trends, characterization of the dysregulated transcriptome in c-Maf–deficient NKp46+ ILC3s using gene set enrichment analysis (GSEA) revealed that the most down-regulated genes were significantly enriched in ILC3 signature genes, whereas the most up-regulated targets were significantly enriched for ILC1 signature genes, as well as for top-ranking genes positively regulated by T-bet (Fig. 3 B). These global shifts in effector program reveal a role for c-Maf in both stabilizing type 3 identity and constraining type 1 features in NKp46+ ILC3s.
c-Maf directly regulates the Tbx21 locus
We next evaluated regulation of Tbx21 by c-Maf. For this, we purified NKp46+ cells from the ILC3-like MNK-3 cell line, which recapitulates the effector features and the c-Maf, RORγt, and T-bet expression pattern of intestinal ILC3s (Fig. S3 E; Allan et al., 2015). Overexpression of c-Maf in NKp46+ MNK-3 cells via retroviral transduction resulted in decreased T-bet expression relative to empty virus control cells (Fig. 3 D), indicating that c-Maf is sufficient to restrict T-bet expression. We identified several Maf recognition elements (MAREs) within conserved noncoding sequences (CNSs) located within, and proximal to, the Tbx21 locus that corresponded to accessible regions in NKp46+ ILC3s as measured by the assay for transposase-accessible chromatin using sequencing (ATAC-seq; Fig. 3 E). As commercial c-Maf antibodies are not compatible with chromatin immunoprecipitation (ChIP) of rare primary cells, ChIP assay was performed using NKp46+ MNK-3 cells. This revealed significant binding of c-Maf relative to a negative control region (3′ UTR) at both CNS+10, located 10 kb from the transcription start site (TSS), and at exon 1 (Fig. 3 E), indicating that c-Maf directly attenuates Tbx21 expression. To identify potential locus silencer elements, Tbx21 CNS regions were cloned upstream of the active Rorc CNS+10 enhancer in a minimal promoter-driven luciferase reporter construct (pEnh), and assayed in NKp46+ MNK-3 cells for the ability to diminish enhancer activity. Interestingly, none of the candidate CNSs displayed silencing activity; rather, two regions, Tbx21 CNS+10 and +19, significantly enhanced reporter pEnh activity (Fig. 3 F). Focusing on the c-Maf–bound region, we confirmed that Tbx21 CNS+10 possesses robust enhancer activity, with a 15-fold increase in activity over that of the minimal promoter (Fig. 3 F). To assess whether c-Maf modulates the activity of CNS+10, we performed additional enhancer reporter assays in the context of c-Maf overexpression. Notably, Tbx21 CNS+10 activity was significantly attenuated by 30% in the presence of exogenous c-Maf expression compared with empty control plasmid, whereas the activity of a Rorc control enhancer was unaffected (Fig. 3 F). This implicates Tbx21 CNS+10 as a c-Maf–responsive enhancer through which c-Maf directs repressive activity at the Tbx21 locus and fine-tunes T-bet expression.
As c-Maf is a known regulator of Rorc (Tanaka et al., 2014; Zuberbuehler et al., 2019), we evaluated whether c-Maf directly activates Rorc expression in ILC3s. c-Maf ChIP of CNS regions that are accessible in NKp46+ ILC3s and harbor MAREs (Fig. 3 G) revealed a modest but significant binding of c-Maf at Rorc CNS-1.5 (Fig. 3 G). This suggests that c-Maf may contribute to Rorc(t) expression in NKp46+ ILC3s.
T-bet–dependent and –independent regulation of ILC3 identity by c-Maf
To determine the contribution of elevated T-bet to the c-Maf–deficient transcriptome, we performed differential expression analysis using RNA-seq data reported for NKp46+ ILC3 cells isolated from Tbx21+/− heterozygous versus Tbx21+/+ Rorc(t)GFP/+ reporter mice (Rankin et al., 2016). Comparing the c-Maf– and T-bet–dependent changes in transcription revealed a subset of genes showing opposing regulation by c-Maf and T-bet (Fig. 4 A), consistent with c-Maf antagonism of T-bet. As expected, this set was dominated by known T-bet targets including type 1 effector genes such as Ifng, Fasl, Gzma, and members of the killer-cell lectin-like receptor (KLR) family. Together, these findings suggest that the type 1 effector switch in Maf-deficient NKp46+ ILC3s is driven largely by T-bet up-regulation. In support of this, reversing the elevation of T-bet expression in Maffl/flTbx21fl/+Il7rCre genotype mice resulted in a near-restoration of SILP ILC3 subset homeostasis and a major reduction in IFNγ production compared with Maffl/flIl7rCre mice (Fig. 4 B). Nevertheless, IFNγ levels remained modestly higher than that of NKp46+ ILC3s from Maf+/+Il7rCre control mice, suggesting minor contributions of c-Maf to regulation of IFNγ expression.
Notwithstanding the influence of T-bet, the majority of c-Maf–dependent genes were not coregulated by T-bet (Fig. 4 A), revealing T-bet–independent functions for c-Maf in NKp46+ ILC3s. To delineate the distinct roles of c-Maf, we performed pathway analysis of genes that were significantly dysregulated in Maffl/flIl7rCre NKp46+ ILC3s and classified these genes into the following groups: (i) up-regulated, (ii) down-regulated, (iii) T-bet–dependent, or (iv) T-bet–independent (Fig. 4 C). Interestingly, hierarchical clustering of the enriched pathway results demonstrated that down-regulated and T-bet–independent gene sets were more closely related, as were up-regulated and T-bet–dependent subsets. The latter cluster was selectively enriched for type 1–related Th1, NK cell, and Granzyme A signaling pathways, in addition to Notch and mTor signaling pathways that are highly enriched in the ILC1 transcriptome (Gury-BenAri et al., 2016; Fig. 4 C). Other pathways related to vascular endothelial growth factor, STAT3, and T cell receptor signaling were represented in all c-Maf–dependent classes (Fig. 4 C). Unique to the down-regulated and T-bet–independent cluster of c-Maf–regulated loci, we identified a striking enrichment for metabolic pathways related to cholesterol biosynthesis, Ahr signaling, xenobiotic metabolism, and several nuclear receptor–related pathways (e.g., LXR, RXR, and TR; Fig. 4, C and D). The cholesterol biosynthesis pathway is particularly interesting, as cholesterol biosynthetic intermediates function as endogenous activating ligands for RORγ (Santori et al., 2015). In the absence of c-Maf, many genes in this pathway were significantly down-regulated (Fig. 4 D), including Srebf2 (SREBP-2), the activating master regulator TF for the cholesterol biosynthesis pathway (Yamauchi and Rogers, 2018). Moreover, genes encoding many biosynthetic enzymes that function upstream of RORγ ligand production were selectively down-regulated, potentially leading to reduced cholesterol biosynthetic intermediates and therefore decreased RORγt activity. Accordingly, Il22, a direct RORγt target (Qiu et al., 2012), was significantly down-regulated in NKp46+ ILC3s lacking c-Maf. Therefore, c-Maf may stabilize ILC3 identity by both restricting T-bet expression and promoting RORγt function.
c-Maf limits ILC3→ILC1 plasticity
In certain cytokine environments, NKp46+ ILC3 cells convert to RORγt− ex-ILC3s, which express high levels of T-bet and IFNγ and are phenotypically identical to ILC1 cells (Vonarbourg et al., 2010; Klose et al., 2013). The dysregulated state of c-Maf–deficient NKp46+ ILC3s is reminiscent of such conversion, suggesting that c-Maf functions to limit ILC3 plasticity. To evaluate this, we employed an ILC3 fate-mapping strategy. Specifically, Maffl/fl mice were bred with the Rorc-Cre Rosa26Zs-Green strain, which activates Zs-Green expression from a conditional allele in the ubiquitous Rosa26 locus, permanently marking cells with a history of Rorc expression, and allowing us to distinguish between bona fide Zs-Green− ILC1s and Zs-Green+ ex-ILC3s (Fig. 5 A). Interestingly, tracking ILC3 plasticity in WT Rorc-Cre Rosa26Zs-Green mice, we found that c-Maf expression is significantly reduced in ex-ILC3s relative to NKp46+ ILC3s, to the level observed for ILC1s, suggesting that down-regulation of c-Maf activity represents a physiological mechanism facilitating ILC3→ILC1-like effector transitions (Fig. 5 B). Consistent with this hypothesis, compared with WT controls, the CoLP from Maffl/flRorc-Cre Rosa26Zs-Green mice had an increased frequency of ex-ILC3s among total ILC1s (Zs-Green+ and Zs-Green−), indicative of enhanced conversion in the absence of c-Maf (Fig. 5 C). Importantly, the increased total ILC1s in Maffl/fl Rorc-Cre mice was attributable to elevated numbers of ex-ILC3s, as numbers of fate-mapped negative ILC1s (ZS−) were not significantly altered (Fig. 5 C). Interestingly, accelerated conversion was already observed in the E17 fetal intestine, where a substantial proportion of ex-ILC3s were detected in Maffl/flRorc-Cre Rosa26Zs-Green but not control embryos (Fig. 5 D). Taken together, c-Maf restrains physiological ILC3 plasticity and conversion to a type 1 ILC effector at steady state and during ontogeny.
Cytokine signals are key regulators of ILC3 identity. During transition into ex-ILC3s, NKp46+ ILC3s down-regulate IL-7Rα (CD127; Klose et al., 2014). Interestingly, CD127 expression was significantly reduced on Maffl/flRorc-Cre Rosa26Zs-Green Nkp46+ ILC3s compared with controls (Fig. 6 A). c-Maf ChIP assay revealed significant c-Maf occupancy at the Il7r promoter and CNS-5.4 (Fig. 6 B), implicating c-Maf in direct regulation of Il7r expression. Similarly, c-Maf–deficient ex-ILC3s also possessed lower levels of CD127 than WT ex-ILC3s, uniformly matching that of control ILC1s (ZS−; Fig. 6 A). These shifts in CD127 expression are consistent with the advanced acquisition of a type 1–like profile by NKp46+ ILC3 cells in the absence of c-Maf. However, since IL-7 signals functionally maintain ILC3 identity by promoting RORγt expression (Vonarbourg et al., 2010), c-Maf–dependent changes in IL-7Rα expression on NKp46+ ILC3s may also contribute to ILC1 conversion.
In light of the involvement of cytokines such as IL-7, IL-23, IL-12, IL-15, and IL-18 in promoting ILC3 identity or ILC1 conversion (Vonarbourg et al., 2010; Bernink et al., 2013, 2015; Fuchs et al., 2013; Klose et al., 2013), we next explored whether these signals can modulate c-Maf expression. For this, we cultured NKp46+ ILC3s from Rag1−/− mice with the OP9-DL1 stromal cell line for 24 h. Control cells were maintained in the presence of stem cell factor (SCF) and Flt3L, whereas treatment groups received additional test cytokines. Among these conditions, IL-7 and IL-23 enhanced c-Maf expression relative to the control, whereas IL-12 and IL-15 had little or no effect (Fig. 6, C and D). Moreover, addition of IL-12 to IL-7 was also not sufficient to down-regulate c-Maf to control levels (Fig. 6 D). Importantly, this was not due to lack of IL-12 activity, as IL-12 was able to synergize with IL-18 for IFNγ production (Fig. S3 F). Unexpectedly, IL-18 was a strong inducer of both c-Maf expression and IFNγ, which supports a pleiotropic role for this cytokine in promoting both type 1 and type 3 features, as reported by others (Bernink et al., 2015; Victor et al., 2017). Together, these trends imply that lineage-relevant cytokines promote c-Maf rather than target its reduction in NKp46+ ILC3s. In particular, c-Maf enhancement by IL-7, combined with the down-regulation of CD127 in c-Maf–deficient cells suggests a positive feedback loop model, whereby IL-7 enhances c-Maf expression that in turn supports CD127 expression to maintain NKp46+ ILC3 identity.
c-Maf supports an ILC3 chromatin accessibility landscape
To uncover the molecular mechanisms of c-Maf–mediated gene regulation of ILC3 identity, we used ATAC-seq to compare chromatin accessibility in NKp46+ ILC3s or CCR6+ ILC3s from Maffl/flIl7rCre versus control Maf+/+Il7rCre SILP. Among the total 106,251 accessible peaks in NKp46+ ILC3 and CCR6+ ILC3 subsets, ∼61% were shared (Fig. S4 A). However, among regions that were differentially accessible (DA) in the absence of c-Maf, 3,430 were unique in NKp46+ ILC3s compared with 1,177 unique in CCR6+ ILC3s, and only 1,014 shared, suggesting that c-Maf has distinct roles in each subset (Fig. S4 B). In this regard, 8.5% of ATAC regions were significantly DA in Maffl/flIl7rCre NKp46+ ILC3s (Fig. 7 A). Regions with diminished accessibility included those associated with canonical c-Maf targets or ILC3-related loci (e.g., Blk, Icos, Vdr, Podnl1, Cxcl2; Fig. 7 A), several of which were also DA in CCR6+ ILC3s (Fig. S4 C). Conversely, Tbx21 and T-bet–dependent loci such as Ifng, Klri2, Plac8, and Fasl encompassed regions displaying increased accessibility in Maf-deficient NKp46+ ILC3s to levels resembling that of ILC1 cells (Fig. 7 B and S4 D). In particular, we observed a significant increase in accessibility in the Tbx21 CNS+10 region that is occupied by c-Maf. This was correlated with a Tbx21 locus-wide significant increase in the active enhancer histone mark H3K27ac, particularly between c-Maf–bound regions at exon 1 and CNS+10, with no change in repressive H3K27me3 as measured by CUTRUN (cleavage under targets and release using nuclease) chromatin profiling (Skene and Henikoff, 2017; Fig. 7, B and C). This indicates that c-Maf suppresses Tbx21 by influencing histone acetylation in addition to chromatin accessibility. Moreover, c-Maf–dependent gains and losses in ATAC accessibility at promoter regions were well-correlated (r = 0.68) with differential expression of the corresponding loci in Maffl/flIl7rCre NKp46+ ILC3s (Fig. S4 E). Notably, among these DE loci, 29% had at least one DA region within 10 kb of the TSS (compared with 7.25% of nondependent loci, P = 4.80 × 10−41, Fisher’s exact test), signifying that c-Maf–dependent changes in chromatin accessibility contribute to shifts in gene expression associated with ILC3 identity and plasticity.
To identify factors that mediate c-Maf–dependent changes in chromatin accessibility, we performed de novo motif analysis of DA ATAC regions that were either shared or uniquely affected in NKp46+ or CCR6+ ILC3 subsets (Fig. 7, D and E; and Fig. S4 F). This identified MARE, NF-κB, and ROR response element consensus sites as the three most enriched motifs for NKp46+ ILC3-specific DA regions losing accessibility in the absence of c-Maf (Fig. 7 E). Here, DA regions containing MAREs in combination with NF-κB or RORγt motifs, or both, were more highly represented compared with non-DA regions (Fig. 7 F). Together, this implies a role for c-Maf in direct promotion of accessibility, potentially in collaboration with NF-κB and RORγt. Conversely, the MARE consensus was not identified in regions increasing in accessibility in NKp46+ ILC3s; instead, Runx, Tbox, and HMG motifs were highly enriched (Fig. 7 E). Notably, like Tbx21, c-Maf is bound to CNSs within the Runx3 and Tcf7 (TCF-1) loci (Fig. 7 H), which also display up-regulated RNA in the absence of c-Maf (Fig. 7 G), indicating that c-Maf constrains the expression of these TFs in addition to limiting their function at regulatory elements. Although the Runx motif is highly represented in ATAC regions, those that contained combinations of the Runx motif with T-bet, Tcf7, or both were selectively enriched in increasing DA regions, compared with regions that harbor only one of these motifs (Fig. 7 F). This finding is consistent with the described role for Runx as a T-bet coregulator in type 1 or plastic type 17 gene regulation (Djuretic et al., 2007; Wang et al., 2014) and implies a similar function in ILC3s. Importantly, while the MARE, Runx, and Tcf motifs were enriched in DA regions encompassed in both ILC3 subsets, enrichment of motifs for lineage-defining TFs (RORγt and T-bet) was distinct to the NKp46+ ILC3 DA regions (Figs. 7 E and S4 F). Therefore, c-Maf uses both direct and indirect mechanisms to regulate the changes in chromatin accessibility associated with identity in NKp46+ ILC3 cells.
c-Maf globally constrains transition to a type 1 chromatin accessibility landscape
We next aimed to place the c-Maf–regulated chromatin state in the context of the changes in the chromatin landscape that occur during the transition of ILC3s to the ILC1 state. For this, we generated additional ATAC-seq profiles from developmentally associated DN and NKp46+ ILC3s isolated from the SILP of C57Bl/6 mice (Fig. 8 A), and RORγthi NKp46+ ILC3s, transitioning RORγtlo NKp46+ ILC3s, and RORγt− ILC1s identified based on GFP reporter expression from RORγtGFP/+ mice (Eberl and Littman, 2004; Fig. 8 A). Consistent with the close relationship between populations, 56% of the 107,680 total accessible regions were shared by all five subsets, with many additional regions shared among two or more subsets, and few regions uniquely accessible in only one subset (Fig. S5 A). Accordingly, and in line with the ILC chromatin landscape being established early in development (Shih et al., 2016), only a small fraction of regions were significantly DA between developmental stages in the ILC3→ILC1 transition (Fig. 8 B). Among these, as expected, the ILC3-specific locus, Il22, showed progressively diminishing accessibility, while the locus for the ILC1 lineage-defining TF, Tbx21, displayed increasing accessibility with conversion (Fig. 8 C). Motif analysis of DA ATAC regions losing or gaining accessibility in consecutive transitions between NKp46+ ILC3 and ILC1 cells implicate the concerted loss of RORγt, NF-κB, and Maf activity, and enhanced T-bet and Ets function, respectively, in the conversion of ILC3 to ILC1 (Figs. 8 D and S5 C).
Lineage relatedness was assessed using principal component analysis (PCA) of the ATAC-seq accessibility profiles to visualize the main components of variability among the ILC3 to ILC1 lineage series. The CCR6+ ILC3, DN ILC3, WT NKp46+ ILC3, and ILC1 population samples clustered separately in four distinct quadrants (Fig. 8 E). Interestingly, the c-Maf–deficient NKp46+ ILC3 accessibility profile clustered closely with those of transitional NKp46+ RORγtlo ILC3s from RORγtGFP/+ reporter mice and showed greater similarity to ILC1s than to Maf+/+ NKp46+ ILC3s (Fig. 8 E), indicative of a shift toward an ILC1-like chromatin landscape. The resemblance and type 1 shift of both c-Maf–deficient NKp46+ ILC3s and RORγt heterozygous NKp46+ ILC3s from RORγtGFP/+ mice are notable, as RORγtGFP/+ ILC3s phenocopy Maffl/flIl7rCre ILC3s with respect to increased frequency of NKp46+ subset cells and enhanced expression of T-bet (Fig. S5 D; Zhong et al., 2016). Thus, c-Maf–deficient NKp46+ ILC3s have altered accessibility resembling that of a transitional state between ILC3 and ILC1 cells.
Hierarchical clustering of the 500 most variable ATAC regions for the DN ILC3→ILC1 transition series revealed two distinct clusters representing the balanced loss of ILC3 and acquisition of ILC1 identity (Figs. 8 F and S5 B). Specifically, cluster 1 encompassed regions decreasing in accessibility during the transition and was most significantly enriched for the RORγt motif, whereas cluster 2 included regions gaining accessibility from DN ILC3 to ILC1 and was dominated by the T-bet motif (Fig. 8 F). Accordingly, cluster 1 ATAC regions were associated with loci positively regulated by c-Maf and those more highly expressed in ILC3 compared with ILC1 (Fig. 8 G), whereas cluster 2 ATAC regions localize to T-bet–dependent and ILC1 signature genes (Fig. 8 G). In line with the PCA, the Maffl/flIl7rCre NKp46+ ILC3 ATAC signature clustered with that of the transitional NKp46+ ILC3 and ILC1 populations isolated from RORγtGFP/+ mice, unlike the Maf+/+Il7rCre NKp46+ ILC3 profile that clustered with that of DN ILC3 and NKp46+ ILC3s from C57Bl/6 mice (Fig. 8 E). Thus, the c-Maf–deficient regulatory landscape resides within the continuum of transitional regulatory states between NKp46+ ILC3 and ILC1. This is driven by the precocious and predominant acquisition of ILC1-like accessibility (cluster 2), and to a much lesser extent, the loss of ILC3 accessibility (cluster 1). Interestingly, RORγt heterozygosity (Rorc(t)GFP/+) alone was sufficient to more uniformly shift NKp46+ ILC3 accessibility to an ILC1-like profile, suggesting that loss of c-Maf activity precedes loss of RORγt function, and that global gains in ILC1-like chromatin states precede the silencing of ILC3 regions during plastic conversion. This later trend is further supported by the number of DA regions observed between transitional states, whereby the majority of gains in accessibility are acquired early on, coincident with T-bet up-regulation, and proportionally more regions lose accessibility in the final transition to ILC1 (Fig. 8 B). Taken together, these data show that c-Maf sustains type 3 identity by constraining a global regulatory transition to a type 1 chromatin accessibility landscape.
Discussion
Here, we identified c-Maf as an essential regulator of CCR6− ILC3 homeostasis and plasticity. c-Maf controls development at two distinct T-bet–dependent stages by limiting the initial DN→NKp46+ ILC3 transition, and subsequent NKp46+ ILC3→ILC1 effector conversion. Mechanistically, restraint of T-bet expression by c-Maf is required to regulate the T-bet expression gradient that is essential for proper CCR6− ILC3 developmental progression and to maintain ILC3 identity. Accordingly, c-Maf constrained expression of IFNγ, NK cell–related loci, and other type 1 loci (e.g., Tcf7 and Runx3), while conversely supporting RORγt function and the expression of Il22 and key type 3 effectors. Underlying lineage-defining transcriptional programming, c-Maf is required to both promote ILC3 chromatin accessibility and restrain the acquisition of an ILC1-like regulatory landscape. Taken together, our findings define c-Maf as a gatekeeper of type 1 effector acquisition that regulates physiological plasticity and stabilizes NKp46+ ILC3 cell identity.
c-Maf plays distinct roles in adaptive versus innate type 17 programming. In Th17 cells, c-Maf is a core regulator of effector genes, directly controlling cytokine, anti-inflammatory, and tissue residency gene expression programs, yet notably, not Rorc master regulator expression (Xu et al., 2009; Ciofani et al., 2012; Aschenbrenner et al., 2018; Gabryšová et al., 2018). In contrast, in innate programming of RORγt+ lymphoid effectors, c-Maf directly controls the expression of lineage-defining TFs to influence lineage identity and plasticity. Indeed, here, we show that c-Maf directly binds Tbx21 and limits expression in NKp46+ ILC3s to constrain acquisition of type 1 features and conversion to an ILC1 fate. In addition, c-Maf is essential for the activation of Rorc and type 17 program genes during effector commitment of innate-like Tγδ17 cells (Zuberbuehler et al., 2019). The context-dependent function of c-Maf is highlighted by distinctions in epigenetic programming, whereby c-Maf contributes to the establishment of an effector-associated regulatory landscape in innate ILC3 and Tγδ17 cells (Zuberbuehler et al., 2019), but not in adaptive Th17 cells (Gabryšová et al., 2018; Miraldi et al., 2019). Rather, in Th17 cells, regulation of lineage-defining TF expression and chromatin remodeling is mediated by the TFs IRF4, BATF, and STAT3 (Ciofani et al., 2012; Vahedi et al., 2012). Thus, c-Maf has a more prominent position in the regulatory hierarchy of factors controlling the specialization of innate compared with adaptive RORγt+ lymphocyte effectors.
NKp46+ ILC3s are poised for type 1 functionality. While environmental factors that promote type 1 conversion have been reported (e.g., IL-12, IL-18, and IL-15; Vonarbourg et al., 2010; Bernink et al., 2013, 2015), the molecular mechanisms that govern ILC3 plasticity are mostly undefined. Flexibility is facilitated by the simultaneous expression of type 1 and type 3 lineage-defining TFs: T-bet and RORγt (Sciumé et al., 2012; Klose et al., 2013). The balance between these mutually antagonistic regulators is critical for CCR6− ILC3 development and homeostasis (Klose et al., 2013; Zhong et al., 2016). Accordingly, we show that reducing levels of RORγt via genetic heterozygosity is sufficient to shift the accessibility of regions that are altered during conversion in the direction of an ILC1-like profile. While loss of RORγt accompanies the ILC3→ILC1 transition, its down-modulation by GATA-3 is also required to permit development of ILC3s into the NKp46+ ILC3 compartment (Zhong et al., 2016). What regulates T-bet expression to control the lineage-defining TF balance? Here, we show that c-Maf binds Tbx21 and restrains T-bet expression in CCR6− ILC3s. c-Maf is essential for the graded expression of T-bet that underlies the homeostasis of populations in the continuum of the DN ILC3→NKp46+ ILC3→ILC1 (ex-ILC3) progression. In this way, c-Maf is a rheostat that tunes T-bet levels, limits the acquisition of type 1 features, and facilitates the poised regulatory status of NKp46+ ILC3s.
c-Maf also limits the expression of other type 1–associated TFs, Runx3 and Tcf7 (TCF-1), that influence NKp46+ ILC3 and/or ILC1 development (Mielke et al., 2013; Ebihara et al., 2015). Interestingly, c-Maf also antagonizes TCF-1 in Tγδ17 cells (Zuberbuehler et al., 2019), revealing a conserved c-Maf–TCF-1 regulatory axis in innate-like RORγt+ effector lymphocytes. Enrichment of the motif trio, T-bet, Runx, and TCF-1, within regions that gain accessibility in the absence of c-Maf suggests that c-Maf also constrains Runx3 and TCF-1 function at suppressed regulatory elements that potentially drive type 1 conversion in WT NKp46+ ILC3s. In support of this notion, TCF-1 promotes the type 1 program in γδ T cells, while counteracting type 17 effector differentiation of Tγδ17, Th17, and Tc17 cells (Yu et al., 2011; Malhotra et al., 2013; Mielke et al., 2019), implying a similar role for TCF-1 in ILC3 plastic conversion. In addition, Runx3 is a cooperative binding partner of T-bet in Th1 gene regulation (Djuretic et al., 2007) and an essential coregulator of T-bet in Th17 cell conversion to an IFNγ+ pathogenic state (Wang et al., 2014). Thus, it is tempting to speculate that a similar T-bet–Runx3 collaboration mediates the ILC3→ILC1 transition when c-Maf activity is dampened. Taken together, these findings show that c-Maf prohibits the activation of a type 1 regulatory module in NKp46+ ILC3s.
c-Maf directly regulates ILC3 identity. Indeed, the ILC3 signature is globally down-regulated in the absence of c-Maf. Moreover, the enrichment of the MARE motif selectively at cis regions that lose accessibility with Maf deficiency, and that are proximal to c-Maf–dependent loci, indicates that c-Maf activates genes in the NKp46+ ILC3 program. The colocalization of motifs for ROR response element and NF-κB p65 with MAREs in these c-Maf–dependent ATAC regions suggests that positive regulation by c-Maf is integrated with that of lineage-defining TF RORγt and environmental signals (potentially IL-1β) that result in NF-κB activation. Thus, c-Maf globally regulates accessibility associated with ILC3 identity and function. In addition, c-Maf supports lineage-promoting metabolic pathways, such as the AhR and RXR pathways that are critical to ILC3 development (Qiu et al., 2012; van de Pavert et al., 2014). c-Maf also activates the expression of Srebf2, the master regulator of the cholesterol biosynthesis pathway, and in particular, selectively supports expression of enzymes involved in the generation of biosynthetic intermediates that are endogenous RORγ-activating ligands (Santori et al., 2015). Notably, c-Maf similarly activates cholesterol biosynthesis genes in Schwann cells, where it is required for proper myelination and cellular function (Kim et al., 2018). Together, the data support a model in which c-Maf balances suppression of T-bet with metabolic programming that supports RORγt activity and maintenance of NKp46+ ILC3 identity.
Environmental signals drive CCR6− ILC3 development and plasticity. T-bet up-regulation in CCR6− ILC3 is controlled by cues from the commensal microbiota and IL-23 (Klose et al., 2013). Thereafter, Notch signaling is essential for transition to the NKp46+ ILC3 stage (Rankin et al., 2013). Notably, these pathways intersect with known regulation of c-Maf, as both Notch and cytokine-mediated STAT3 signals are activators of c-Maf in CD4+ T cells (Xu et al., 2009; Auderset et al., 2013). Moreover, the microbiota supports IL-7 production by intestinal epithelial cells (Shalapour et al., 2010; Vonarbourg et al., 2010), a cytokine that we demonstrate also enhances c-Maf expression by NKp46+ ILC3 cells. This implies a model in which these environmental signals sequentially boost c-Maf levels during the DN→NKp46+ ILC3 transition, facilitating graded T-bet expression. Thereafter, physiological down-regulation of IL-7R expression by T-bet (Klose et al., 2013) in NKp46+ ILC3s may allow negative feedback down-modulation of c-Maf and release of the type 1 program, resulting in full conversion to an ILC1 state. Taking these findings together, we have identified a core regulator of ILC3 identity and plasticity that globally restrains the ILC1-like program and controls the switch from type 3 to type 1 immunity in T-bet+ ILC3 cells.
Materials and methods
Mice
Mice were housed under specific pathogen–free conditions and used in accordance with the Duke University Institutional Animal Care guidelines. Maf floxed allele mice (Maffl/fl) were obtained from C. Birchmeier (Max-Delbrück Center for Molecular Medicine, Berlin, Germany; Wende et al., 2012) and were backcrossed to C57BL/6 for at least five generations. Il7rCre mice were obtained from H.R. Rodewald (German Cancer Research Center, Heidelberg, Germany; Schlenner et al., 2010). The following commercially available strains were bred in our facility: C57Bl/6 mice (Taconic), Rag2−/−Il2rg−/− mice (stock 4111; Taconic), Rorc-Cre mice (stock 022791; Jackson), Rorc(t)GFP mice (stock 007572; Jackson), Rosa26ZS-green (stock 007906; Jackson), Tbx21 floxed allele mice (stock 022741; Jackson), and Rag1-deficient mice (stock 002216; Jackson). Adult mice were used at the ages of 8–16 wk. Littermate controls were used.
Generation of bone marrow chimeras
Bone marrow cells were harvested from either CD45.2+ Maf+/+Rorc-Cre or Maffl/flRorc-Cre mice and were mixed in a 1:1 ratio with CD45.2+CD45.1+ Maf+/+Rorc-Cre bone marrow cells. This mixture was retroorbitally injected into sublethally irradiated (450 rad) Rag2−/−Il2rg−/− mice. SILP was harvested 6 wk after reconstitution for analysis.
Cell isolation from intestines
Lamina propria lymphocytes from the small and large intestine were isolated as previously described (Carr et al., 2017), except the 30-min digestion occurred in the presence of 2% FBS/10 mM Hepes for the small intestine or 10% FBS/10 mM Hepes for the large intestine. Fetal intestines were harvested, minced, and digested for 60 min in a water bath at 37°C in HBSS supplemented with collagenase D (1 mg/ml; Roche), DNaseI (0.1 mg/ml; Sigma-Aldrich), dispase (0.1 U/ml; Worthington), 10% FBS, and 10 mM Hepes. The remaining isolation steps were as described for the adult intestine protocol (Carr et al., 2017).
Cell stimulation and flow cytometry
To evaluate cytokine expression, ex vivo lamina propria lymphocytes were stimulated for 4 h at 37°C in complete IMDM (IMDM supplemented with 10% FBS, 10 U ml−1 penicillin, 10 µg ml−1 streptomycin, 2 mM glutamine, 50 µg ml−1 gentamycin, and 55 µM β-mercaptoethanol) with phorbol 12-myristate 13-acetate (100 ng ml−1; Sigma-Aldrich) and ionomycin (375 ng ml−1; Sigma-Aldrich) or IL-23 (20 ng ml−1) and IL-1β (10 ng ml−1) in the presence of GolgiStop (BD) for the last 3 h.
For both stimulated and unstimulated samples, Fc receptor blocking, surface stain, fixation/permeabilization, and TF staining was performed as described (Carr et al., 2017). When sorting, cells were stained in PBS for 30 min at 4°C, washed, and resuspended in Ca2+/Mg2+-free PBS buffer containing 2 mM EDTA and 0.5% BSA. Sorting was performed using a MoFlo Astrios or MoFlo XDP cell sorter (Beckman Coulter). In all experiments, a fixable viability dye (eBioscience) was used to exclude dead cells. Antibodies were purchased from BioLegend (CD127 and CCR6), BD Biosciences (NK1.1, NKG2a/c/e, CCR6, IFNγ, and CD45.2), and eBioscience (CD45, RORγt, CD3, CD11b, Ter1119, B220, CD19, T-bet, CD90.2, Streptavidin, IL-22, NKp46, KLRG1, NK1.1, CD45.1, and c-Maf). All data were acquired on a FACSCanto II (BD Biosciences) or Fortessa X20 (BD Biosciences), and analyzed using FlowJo software (TreeStar). All flow cytometry data were pregated for singlets and live cells.
Cell culture
OP9-DL1 cells (provided by J.C. Zúñgia-Pflücker, Sunnybrook Research Institute, Toronto, Canada) were maintained in complete IMDM (Sigma-Aldrich). Rag1−/− NKp46+ ILC3s were sorted (CD45loCD90hi, KLRG1−, Nkp46+CCR6−) and cultured on OP9-DL1 stromal cells in complete IMDM in the presence of SCF (10 ng ml−1; eBioscience) and Flt3L (10 ng ml−1; eBioscience) for 24 h with additional cytokines as indicated. IL-7 (10 ng ml−1; eBioscience), IL-23 (50 ng ml−1; eBioscience), IL-12 (10 ng ml−1 for Fig. 6 D and 50 ng ml−1 for Fig. 6 C; eBioscience), IL-15 (10 ng ml−1; BioLegend), and IL-18 (50 ng ml−1; BioLegend) were used as well. MNK cells were cultured as described (Allan et al., 2015). MNK-3 cells were sorted to select for NKp46+ cells and are referred to as NKp46+ MNK-3. NKp46+ MNK-3 were cultured in IL-7 (10 ng ml−1; eBioscience) and IL-15 (10 ng ml−1; BioLegend).
Retroviral gene transfer
The retroviral construct was generated by cloning complementary DNA for c-Maf into murine stem cell virus–Thy1.1 5′ of the internal ribosomal entry site, allowing bicistronic expression with cell surface Thy1.1. Retroviral supernatants were generated by transfection of retroviral constructs into the Plat-E producer cell line (Morita et al., 2000) using Lipofectamine 2000 reagent (Thermo Fisher Scientific) and collection of culture medium after 48 h. For gene transfer, 500,000 NKp46+ MNK-3 cells were resuspended in 0.45-µm-filtered viral supernatant containing 6.7 µg ml−1 of hexadimethrine bromide (Sigma-Aldrich), transferred to 48-well plates, and spun at 2,400 rpm for 2 h at 30°C. Cells were cultured for 48 h in fresh medium containing IL-7 and IL-15 (10 ng ml−1) and harvested for flow cytometry analysis.
Dual luciferase reporter assay
Select CNS regions associated with Tbx21 were assessed by dual luciferase reporter assays for enhancer and silencer activity by cloning the sequences upstream of a minimal promoter driving a luciferase gene (pGL4.23 [luc2/minP]) or an experimentally determined enhancer from the Rorc locus (CNS+10), respectively. For this assay, 1 million NKp46+ MNK-3 cells were nucleofected with Lonza Amaxa Cell Line Kit V using program X-005. 1 µg of pCMV-RL was used in combination with 5 µg of the luciferase construct being tested. For c-Maf overexpression luciferase experiments, 5 µg of murine stem cell virus-Thy1.1 c-Maf or empty vector control were conucleofected with pCMV-RL and the luciferase construct being tested. Cells were transferred to prewarmed medium containing IL-7 and IL-15 and were harvested 18 h later. Luciferase assays were done using Promega’s Dual Luciferase Reporter Assay System (E1910). Firefly luciferase measurements were normalized to renilla luciferase readings for each sample. Data were presented as fold change over empty pGL4.23minP or Rorc CNS+10-pGL4.23minP, as appropriate. Regions assayed include the following (mm10): Tbx21 CNS+2.6 chr11:97,112,536–97,113,131; Tbx21 CNS+5 chr11:97,109,626–97,110,193; Tbx21 CNS+10 chr11:97,103,152–97,104,646; Tbx21 CNS+11 chr11:97,103,152–97,103,791; and Tbx21 CNS+19 chr11:97,096,037–97,096,518.
ChIP
20 million NKp46+ MNK-3 cells were used for c-Maf ChIP-qPCR. ChIP was performed as described (Carr et al., 2017) using a commercial c-Maf–specific antibody (Bethyl Laboratories; A300-613A). Primers specific for test regions containing MAREs as well as for a control region within the locus lacking MAREs were used to amplify the ChIP-enriched and input DNA by qPCR. The data are represented as percentage of input DNA. Primers for c-Maf ChIP-qPCR at various regions in the 5′ to 3′ orientation are as follows: Tbx21 CNS+8.4 F: GCTGGGTAATTAGCCTTCGCTCT, R: CACCCTTGTCAGTTTGCGTTGG; Tbx21 Exon 1 F: CCCGGCTCGGGATAGAAGAAAC, R: ATGCTGACCGGCACCGA; Tbx21 CNS+2.6 F: TGACAAGGCATTGAGAGGACTGG, R: GGGCCTCTAACCATGGGACTTT; Tbx21 CNS+10 F: TACGCCTCCGCCATGACCA, R: AGACTAGCGTAATGGCAACTGTGAGT; Tbx21 CNS+11 F: CAGGGCACCAGCAGACACAGTAAAT, R: GATGTGGCTCAGTGGTAGAGTGCT; Tbx21 CNS+19 F: GTGCAGCCTGGGTTGGTCTC, R: CAAGCCAGTTGAGGTGGCAAGT; Tbx21 3′ UTR F: CGACCCCTCCACCTCGCAGAA, R: CAGTCACGAACCTGGTGCTGCTT; Rorc(t) CNS-1.5 F: CATCCCTCCTCCCTGAGCTTGG, R: GGCTAGTGCACAGGAGGCTGAT; Rorc(t) CNS+2.5 F: CTGCCTCTGTTGCTGACTCAAAC, R: TGTCTTCCTTAGAGAGCACTCAGC; Rorc(t) CNS+6 F: CTCCCAGGCAACAGTTTGAGTAATCC, R: CGGAGAAACATCTTCCACAGCAAAG; Rorc(t) CNS+10 F: GGGCATGAGGTGAGGAAGACCCAG, R: CCCTGCATGTTTCCTTTCGCTTTGG; Rorc(t) Neg Control F: CCACAGAGACACCACCGGACATCT, R: CTGACAGTCCCTTCCACCTGCCT; Il7r pro F: TCTTGGTATTCCCTGAAGCTCCTGA, R: TGGACCCTGAGATCATCTGAAGCA; Il7r CNS-5.4 F: GCGAGGGTGGTTAAAGGAACCCAT, R: TTCAAAATGCTGGGGCTCCCCCT; Il7r Neg (Int 3) F: TGTCACATCATTCAGCCCAAACCCA, R: AAACTGTCTGGGACCAAAGCTGGA; Tcf7 CNS+16 F: TAGCCTCCCACTTCCTGCAGGTTC, R: GGGAGTTTCTCTGCTAGTGGTGACAGG; Tcf7 CNS-17 F: GCATGTGGGATTGCTCACGAGATGG, R: CAAAGGTTGGGCCCAGATGGACAAG; Tcf7 Neg F: AAAGAGTGGCTTGCAGTTGGTCACAG, R: CCTGCTCTCTGTGGGTTCCCTCTC; Runx3 Pro F: GGTGAGCCTCGTTCATTCATTCATT, R: GCAGAGTATCATTAAATGGTGGGAAGG; and Runx3 Neg F: GCCCTATGTCTGTCTTGTTCCTGGG, R: TGAACCAGACCCCAAAGCAGCAATC.
Pathway analysis
Ingenuity Pathway Analysis (Qiagen) was performed on DE genes (Benjamini–Hochberg false discovery rate [FDR] 0.05) based on DESeq2 differential expression analysis of RNA-seq data from Maf+/+Il7rCre versus Maffl/flIl7rCre NKp46+ ILC3. Pathways with P value 0.05 were considered significant.
Statistical analysis
Besides RNA-seq, ATAC-seq, and CUTRUN analysis, all data were analyzed using GraphPad Prism 7 or 8. Two-tailed unpaired Student’s t test was used for all flow cytometry phenotyping, luciferase reporter assay, and ChIP-qPCR data.
RNA preparation and RNA-seq
5,000–15,000 sort-purified Nkp46+ ILC3 (CD3−CD19−CD127+CD90hiKLRG1−CCR6−NKp46+) and CCR6+ ILC3 (CD3−CD19−CD127+CD90hiKLRG1−CCR6+NKp46−) and ∼17,000 ILC2 (CD3−CD19−CD127+KLRG1+) from Maf+/+Il7rCre and Maffl/flIl7rCre SILP were frozen in Trizol. RNA was purified from the aqueous phase using the RNeasy Plus Micro Kit (Qiagen). Samples were submitted to the Duke Sequencing and Genomic Technologies Shared Resource facility for library construction using the Clontech SMARTer v3/v4 Ultralow Input RNA-seq Kit (Takara Biosciences). Libraries for three biological replicates per subset were sequenced on the Illumina HiSeq 4000 using 50-bp single-end mode at a depth of 30–44 million reads per sample.
RNA-seq differential expression analysis
All RNA-seq samples were first validated for consistent quality using FastQC v0.11.2 (Babraham Institute). Raw reads were trimmed to remove adapters and bases with average quality score (Q; Phred33) of 20 using a 4-bp sliding window (SlidingWindow:4:20) with Trimmomatic v0.32 (Bolger et al., 2014). Trimmed reads were aligned to the primary assembly of the GRCm38 mouse genome using STAR v2.4.1a (Dobin et al., 2013) removing alignments containing noncanonical splice junctions (--outFilterIntronMotifs RemoveNoncanonical).
Aligned reads were assigned to genes in the GENCODE vM13 comprehensive gene annotation (Harrow et al., 2012) using featureCounts (v1.4.6-p4) with default settings (Liao et al., 2013). Differential expression analysis was performed using DESeq2 (v1.22.0; Love et al., 2014) running on R (v3.5.1). Briefly, raw counts were imported and filtered to remove genes with low or no expression (that is, keeping genes having two or more counts per million in two or more samples). Filtered counts were then normalized using the DESeq function, which internally uses estimated size factors accounting for library size, estimated gene, and global dispersion. To find significant DE genes, the nbinomWaldTest was used to test the coefficients in the fitted Negative Binomial GLM using the previously calculated size factors and dispersion estimates. Genes having an Benjamini–Hochberg FDR 0.05 were considered significant, unless otherwise indicated. Log2 fold-change values were shrunk toward zero using the adaptive shrinkage estimator from the ashr R package (Stephens, 2017). For estimating transcript abundance, transcripts per million (TPM) were computed using the rsem-calculate-expression function in the RSEM v1.2.21 package (Li and Dewey, 2011). Data were visualized using pandas v0.23.3 and seaborn v0.9.0 packages in Python 2.7.11.
ATAC sample preparation and ATAC-seq
Nkp46+ ILC3 and CCR6+ ILC3 from Maf+/+Il7rCre and Maffl/flIl7rCre SILP were used for ATAC-seq (see RNA-seq methods for gating strategy). For the ATAC-seq transitions, we sorted DN ILC3s (CD3−CD19−NK1.1−/CD45loCD90.2hi/KLRG1−/NKp46−CCR6−) and NKp46+ ILC3s (CD3−CD19−NK1.1−/CD45loCD90.2hi/KLRG1−/NKp46+CCR6−) from C57bl/6 mice. Additionally, we sorted NKp46+ ILC3s (CD3−CD19−/CD45+CD127+ NKp46+RORγthi), NKp46+ RORγtlo (CD3−CD19−/CD45+CD127+ NKp46+RORγtlo), and ILC1s (CD3−CD19−/CD45+CD127+NKp46+RORγt−) from Rorc(t)GFP/+ mice. Cells from three to seven mice were pooled before sorting. Omni-ATAC was performed on two biological replicates per subset using Tn5 transposase from the Nextera DNA Library Prep Kit (FC-121-1030). The transposition reaction was scaled according to the cell number. After transposition, Qiagen MinElute Kit (28204) was used to purify the DNA. Based on qPCR amplification curves, samples were amplified for five to seven additional cycles after the five-cycle preamplification (10–12 total cycles). For size selection and clean-up of material, three rounds of AMPure beads were performed. First, 1× AMPure beads were used to remove primer dimers, followed by 0.5× to remove large fragments 1,000 kb, and finally 1.8× to further clean up. ATAC libraries were submitted to the Duke Sequencing and Genomic Technologies Shared Resource facility for sequencing on an Illumina NextSeq500 using a 42-bp paired-end mode at a depth of 60–80 million reads per sample.
ATAC-seq preprocessing and alignment
As with RNA-seq data, quality was assessed with FastQC, and adapters were trimmed with Trimmomatic. Trimmed reads were aligned to the GRCm38 genome using Bowtie (v1.0.0; Langmead et al., 2009), reporting only alignments having no more than two mismatches, discarding multimapping reads (-v 2 --best-strata -m 1). Reads mapping to the ENCODE mm10 blacklisted regions (https://www.encodeproject.org/files/ENCFF547MET; manually curated regions with anomalous signal across multiple genomic assays and cell types) were removed using bedtools2 intersect (v2.25.0; Quinlan and Hall, 2010). Properly paired reads were then filtered to exclude presumed PCR duplicates using Picard MarkDuplicates (v1.130). Reads were used to generate RPM count bigWig files for visualization using deeptools bamCoverage (v3.0.1; Ramírez et al., 2014).
Differential ATAC analysis
Reads processed as above were used to call peaks for each sample individually using MACS2 (v2.1.1.20160309; Zhang et al., 2008). Peaks were called using two different settings previously reported in the literature: in paired-end mode, callpeak --nomodel -f BAMPE -q 0.1 --s; using only the first mate in single-end mode, --extsize 200 -f BED --nomodel -q 0.1 --shift -100. The single-end mode achieved a larger pool of peaks which were combined into a master peakset using bedtools merge for downstream analyses. The number of reads in each sample mapping to each peak was calculated using FeatureCounts. Counts were then processed using DESeq2, removing any peaks having fewer than two counts per million in all but one sample. Counts were normalized using the full library size (number of mapped reads as opposed to number of counts in peaks) and the dispersion estimates as described for RNA-seq samples. Regions having a Benjamini–Hochberg FDR 0.05 were considered to have significantly different accessibility between groups.
Hierarchical clustering
Hierarchical clustering was performed over regularized log-transformed (peak and gene) counts using the single method from the linkage function in the scipy package unless otherwise noted (Jones et al., 2001). The single method, also known as Nearest Point Algorithm, uses the minimum Euclidean distance between all points across pairs of clusters to iteratively build the tree of grouped elements. When possible, cell subpopulations were also hierarchically clustered using the Ward variance minimization algorithm implemented in the ward method. Clusters were visualized as heatmaps using the clustermap function in the seaborn package.
PCA
The top 500 regularized log transformed values sorted by variance were used to perform PCA analysis using the R package prcomp. Briefly, the calculation is performed by singular value decomposition of the data matrix. The results were visualized using ggplot2 in R, summarizing the amount of variance explained by each component.
Motif enrichment analysis
To determine the TF binding motifs that are enriched in DA ATAC-seq peaks, we used HOMER de novo and known motif enrichment analysis (HOMER findMotifsGenome.pl v4.9). Full-length differential ATAC-seq peaks with an FDR 0.05 were provided as foreground, with nondifferential peaks in any comparison selected as background.
GSEA
Gene sets used in the GSEA were generated by taking the top 100 up-regulated or down-regulated genes from DE analyses of RNA-seq data comparing ILC3 to ILC1 (ILC3 and ILC1 signatures; Gury-BenAri et al., 2016). For the T-bet–dependent signature, the top 100 down-regulated genes in NKp46+ ILC3 (Rankin et al., 2016) Tbx21+/− Rorc(t)GFP/+ compared with Tbx21+/+ Rorc(t)GFP/+ RNA-seq data were used for the gene set. Each gene in the NKp46+ ILC3 Maf+/+Il7rCre and Maffl/flIl7rCre DE dataset was assigned a rank by taking the −log10 of the DE P value and multiplying by the sign of the log2 fold-change for that gene. The option “number of permutations = 1,000” was used. For the enrichment statistic in Fig. S4 C, “classic” was chosen with the Hallmarks Cholesterol Biosynthesis pathway and “weighted” for Fig. 3 B ILC3, ILC1, and T-bet–dependent signatures.
CUTRUN sample preparation and sequencing
CUTRUN was performed on ∼13,000 Nkp46+ ILC3s from Maf+/+Il7rCre and Maffl/flIl7rCre SILP as previously described (Skene et al., 2018; see RNA-seq methods for gating strategy). In short, after sorting, the cells were permeabilized with 0.005% digitonin freshly made, bound to concanavalin A beads, incubated with H3K27ac, H3K27me3, or IgG antibodies overnight at 4°C with rotation. The remaining steps were performed as described in the original protocol. The following antibodies were used: H3K27ac (Abcam; ab4729), H3K27me3 (Cell Signaling; 9733S), and rabbit anti-mouse IgG (Abcam; ab46540). CUTRUN libraries were submitted to the Duke Sequencing and Genomic Technologies Shared Resource facility for sequencing on an Illumina NextSeq500 using a 42-bp paired-end mode.
CUTRUN data processing
Paired-end reads were mapped using Bowtie2 (v2.3.4.3; Langmead et al., 2009) against both Saccharomyces cerevisiae (sacCer3) and mouse (mm10) genomes in very-sensitive-local mode, allowing for fragments of ≤700 bp and discarding all but paired concordant alignments (parameters for bowtie2: -I 10 -X 700 --local --very-sensitive-local --no-discordant --no-mixed --no-unal --phred33). Duplicates were marked using Picard MarkDuplicates (v1.130) with standard parameters. Signal files were generated with deeptools bamCoverage (v3.0.1; Ramírez et al., 2014) ignoring duplicates, extending reads to match the fragment size defined by paired mates, and scaling the values by 1e6 divided by the total number of reads mapping against sacCer3 (parameters: --extendReads --ignoreDuplicates --scaleFactor 1e6/sacCer3_TOTAL_NREADS). Finally, to quantify histone signal over specific genomic regions, FeatureCounts was used, assigning fully aligned fragment counts to regions (parameters -p -B).
Data availability
All sequencing data that support the findings of this study have been deposited in the Gene Expression Omnibus of the National Center for Biotechnology Information under accession no. GSE137322.
Online supplemental material
Fig. S1 shows the frequencies of various ILC and ILC3 subsets from Maf+/+Il7rCre and Maffl/flIl7rCre mice as well as Maf+/+Rorc-Cre and Maffl/flRorc-Cre mice. Fig. S2 shows additional phenotyping data for Maf-deficient mice, including cytokine production, fetal intestine ILC subsets, Peyer’s patches cellularity, and summary data for mixed bone marrow chimeras. Fig. S3 shows the sort strategy for RNA-seq and ATAC-seq for Maf+/+Il7rCre and Maffl/flIl7rCre mice and additional RNA-seq–related plots. Fig. S4 shows ATAC-seq–related plots and tracks for loci of interest, correlation of DA and DE, and motif analysis. Fig. S5 shows the shared and unique ATAC peaks for the transitional populations, population clustering, and motif analysis.
Acknowledgments
We thank C. Birchmeier (Max Delbrück Center for Molecular Medicine, Berlin, Germany) for providing Maf conditional mice; J.C. Zúñiga-Pflücker (University of Toronto, Toronto, Canada) for providing OP9-DL1 cells; R. DePooter and J. Espinosa for critical reading of the manuscript; and Eunchong Park and Geoffrey Houtz for technical assistance. We acknowledge the expert assistance of N. Martin and L. Martinek with flow cytometry sorting.
This work was funded by a Whitehead Charitable Foundation Scholar Award (to M. Ciofani). M.E. Parker, J.D. Wheaton, and M. Ciofani were supported by National Institutes of Health grant R01 GM115474.
T.E. Reddy is co-founder and consultant for Element Genomics. The remaining authors declare no competing interests.
Author contributions: M.E. Parker, M.K. Zuberbuehler, and M. Ciofani designed, performed, and analyzed experiments. A. Barrera and J.D. Wheaton performed computational analyses. T.E. Reddy provided computational expertise. D.S.J. Allan and J.R. Carlyle provided reagents. M.E. Parker and M. Ciofani wrote the manuscript. M. Ciofani conceived the study.