Binary outcomes of enhancer activity underlie stable random monoallelic expression

  1. Djem U Kissiov
  2. Alexander Ethell
  3. Sean Chen
  4. Natalie K Wolf
  5. Chenyu Zhang
  6. Susanna M Dang
  7. Yeara Jo
  8. Katrine N Madsen
  9. Ishan Paranjpe
  10. Angus Y Lee
  11. Bryan Chim
  12. Stefan A Muljo
  13. David H Raulet  Is a corresponding author
  1. Division of Immunology and Pathogenesis, Department of Molecular and Cell Biology, University of California, Berkeley, United States
  2. Cancer Research Laboratory, University of California, Berkeley, United States
  3. Laboratory of Immune System Biology, National Institute of Allergy and Infectious Diseases, National Institutes of Health, United States

Abstract

Mitotically stable random monoallelic gene expression (RME) is documented for a small percentage of autosomal genes. We developed an in vivo genetic model to study the role of enhancers in RME using high-resolution single-cell analysis of natural killer (NK) cell receptor gene expression and enhancer deletions in the mouse germline. Enhancers of the RME NK receptor genes were accessible and enriched in H3K27ac on silent and active alleles alike in cells sorted according to allelic expression status, suggesting enhancer activation and gene expression status can be decoupled. In genes with multiple enhancers, enhancer deletion reduced gene expression frequency, in one instance converting the universally expressed gene encoding NKG2D into an RME gene, recapitulating all aspects of natural RME including mitotic stability of both the active and silent states. The results support the binary model of enhancer action, and suggest that RME is a consequence of general properties of gene regulation by enhancers rather than an RME-specific epigenetic program. Therefore, many and perhaps all genes may be subject to some degree of RME. Surprisingly, this was borne out by analysis of several genes that define different major hematopoietic lineages, that were previously thought to be universally expressed within those lineages: the genes encoding NKG2D, CD45, CD8α, and Thy-1. We propose that intrinsically probabilistic gene allele regulation is a general property of enhancer-controlled gene expression, with previously documented RME representing an extreme on a broad continuum.

Editor's evaluation

Kissiov et al., show that enhancers can play an instructive role in controlling stable random monoallelic expression (RME). In order to do so, they initially focus on a limited set of natural killer (NK) receptor genes that are subject to RME, which they investigate using several in vivo genetic models. Furthermore, they show that RME can be considerably more prevalent than previously thought, applying to other gene function than receptors and independently of promoter CpG content. Finally, they provide evidence that enhancer strength and/or number might influence the extent of RME, by affecting the probability of promoter activation. Overall, this is a highly relevant manuscript with major implications in gene regulation and enhancer biology and, thus, of broad scientific interest.

https://doi.org/10.7554/eLife.74204.sa0

Introduction

In most cases, both alleles of autosomal genes are co-expressed. In recent years random monoallelic expression (RME) has emerged as an important exception that may apply to ~0.5–10% of expressed genes in a given tissue and has been characterized as the autosomal analog of X-inactivation (Gendrel et al., 2016). In RME, different cells of a given type express only one allele, both or neither, and this expression pattern is mitotically stable. Notably, RME genes do not share an overarching unifying feature or function (Deng et al., 2014; Eckersley-Maslin and Spector, 2014; Gendrel et al., 2016; Gimelbrant et al., 2007; Reinius et al., 2016; Xu et al., 2017), and the biological role for RME is in most cases not known. RME is distinct from X-inactivation, genomic imprinting, and allelic exclusion of antigen receptor and odorant receptor genes, in that biallelic expression occurs at an appreciable frequency, and expression is largely stochastic rather than imposed by strict feedback regulatory mechanisms (Gendrel et al., 2016).

The molecular determinants of RME are poorly understood, in part because of difficulties analyzing single primary cells (Reinius and Sandberg, 2015). Chromatin analysis of RME alleles in primary populations has not been possible due to the difficulty of isolating pure cell populations ex vivo with defined RME expression patterns (Gendrel et al., 2016; Xu et al., 2017). As a result, previous analyses have been limited to clonal cell lines derived from F1 hybrids, where allelic expression is known and clonally stable (Eckersley-Maslin et al., 2014; Gendrel et al., 2014; Levin-Klein et al., 2017; Xu et al., 2017). These analyses revealed that enhancers of RME genes are constitutively accessible irrespective of gene or allelic expression status, whereas promoters are accessible only at active alleles (Levin-Klein et al., 2017; Xu et al., 2017). Therefore, promoter accessibility, rather than enhancer opening and activation, might be the ‘gatekeeper’ of RME, whereas enhancers, being constitutively open, were proposed to permit rather than impose expression of RME alleles (Xu et al., 2017).

Enhancers may play more than a permissive role in RME, however, in light of evidence that enhancers primarily influence the probability of mitotically stable expression, rather than the amount of expression per cell (Walters et al., 1995; Weintraub, 1988). In fact, deletions of enhancers resulted in mitotically stable gene variegation in both cell lines and normal tissues—notably Igh in B cell hybridomas and Cd8a in primary thymocytes, among others (Ellmeier et al., 2002; Garefalaki et al., 2002; Ronai et al., 1999; Sleckman et al., 1997; Walters et al., 1995; Weintraub, 1988; Xu et al., 1996). Collectively, these data support the binary or ‘on/off’ model of enhancer action (Blackwood and Kadonaga, 1998; Fiering et al., 2000), where an increase in enhancer activity at a genetic locus results in an increase in the probability of gene expression, rather than an increase in expression per cell. Conversely, weak or reduced enhancer activity results in a lower likelihood of expression, but the cells that express the gene express a similar amount of gene product.

We reasoned that the binary action of enhancers—when limiting—might be a driving principle of RME, and sought to test this in examples of RME with a clear biological purpose: the Klra genes (which encode the Ly49 family receptors) and the Klrc1 gene (which encodes the NKG2A receptor) (Chess, 2012; Eckersley-Maslin and Spector, 2014; Gendrel et al., 2016). These genes, clustered in a ~ 1 Mb stretch of the NK gene complex (NKC) on chromosome 6 in mice, encode cell surface receptors expressed by NK cells that bind MHC I molecules. They are expressed in a variegated (Raulet et al., 1997; Yokoyama et al., 1990), monoallelic (Held et al., 1995), stochastic and largely mitotically stable fashion (Raulet et al., 2001), resulting in subpopulations of NK cells that express random combinations of the receptors and consequently exhibit distinct reactivities for cells expressing different MHC I molecules. Regulation of each gene allele is independent, and expression of one Klra gene has minimal effects on expression of others (Tanamachi et al., 2001). While a clear biological purpose for RME at many genes is lacking, RME at Klra genes underlies the basis of the ‘missing self’ mode of NK cell target recognition (Kärre et al., 1986). Furthermore, the system represents a powerful in vivo genetic model of RME, where allelic expression states can be easily assessed at the population level in primary cells using allele-specific antibodies that we and others previously generated (Tanamachi et al., 2001; Vance et al., 2002), circumventing previous technical limitations to studying RME in single primary cells.

Importantly, competition between Klra genes for interaction with a shared enhancer or locus control region is not required for variegation of Klra genes, as a Klra1 (encodes Ly49A) genomic transgene ectopically integrated in different genomic sites was usually expressed with a frequency similar to that of the native Klra1 gene (~17% of NK cells) (Tanamachi et al., 2004). We previously identified a key DNase I hypersensitive element, Klra1Hss1, ~ 5 kb upstream of the Klra1 gene that is conserved in other Klra genes and required for expression of the Klra1 transgene (Tanamachi et al., 2004).

Our central hypothesis is that enhancers, rather than simply being permissive for RME, both limit and directly control the probability of expression of Klra genes—and RME alleles generally—in a stochastic and binary fashion. Binary enhancer action, when limiting, may represent a causal mechanism of RME, explaining the pervasiveness of RME across genes and cell types. We have carried out genetic dissection and population analyses to demonstrate that enhancers control the probability of allelic expression and have provided a more general model of the role of enhancers in RME as well as in other developmentally regulated genes.

Results

Elements upstream of the Klra, Klrc, and Klrk family genes are transcriptional enhancers with activity in mature NK cells

Klra family genes are expressed in a mitotically stable RME fashion by NK cells. Each harbors an accessible chromatin site (Hss1) ~5 kb upstream of the transcription start site (TSS) (Figure 1A and B; Figure 1—figure supplement 1A). We noticed that related NK receptor genes, including the variegated Klrc1 gene (encodes NKG2A) and the Klrk1 gene (encodes NKG2D and is expressed by ~all NK cells), harbor similar elements which we named Klrc15′E and Klrk15′E, respectively (Figure 1A, B; Figure 1—figure supplement 1A). All Hss1 and 5′E elements are bound by a similar suite of factors including Runx3, T-bet, and the enhancer-associated acetyltransferase p300 (Figure 1—figure supplement 1A).

Figure 1 with 2 supplements see all
The Klra1Hss 1, Klrc15′E and Klrk15′E elements display chromatin features of enhancers.

(A) ATAC-seq and H3K4me1:me3 log2 ratio ChIP-seq data of relevant NKC genes in primary NK cells; red denotes positive me1:me3 ratios (enhancer-like) while blue indicates negative values (promoter-like). Approximate gene locations are indicated (bottom). Standard gene names (Klr nomenclature) are indicated followed by names derived from the gene products (Ly49 or Nkg2) for reference. Gray ovals represent additional undiscussed Klra genes. Vertical yellow bars and arrows denote the positions of the Hss1 and 5′E enhancers at the indicated genes. Data are sourced from ref (Lara-Astiaso et al., 2014). Normalized data ranges are indicated on the left. (B) Normalized ChIP-seq and ATAC-seq results (sourced from Lara-Astiaso et al., 2014), showing enhancer and promoter histone modifications at Klra1, Klrc1, and Klrk1. Approximate locations of sgRNAs used in this study to delete enhancers are shown. All datasets are presented with the same vertical scale across sub-panels, which are indicated in normalized signal per million reads (SPMR) in the left sub-panel. (C) Heatmaps depict 51,650 ATAC-seq peaks in primary NK cells (excluding peaks ranking in the bottom 5% for either H3K4me1 or H3K4me3) ranked according to H3K4me1:me3 ratio of average ChIP-seq signal calculated over a 2kb window centered on the ATAC-seq peak midpoint. The indicated data are displayed over these peaks in each heatmap. The locations of selected NK receptor gene Hss1, 5′E and promoter elements within the me1:me3 ranking are shown. H3K4 methylation data are sourced from ref (Lara-Astiaso et al., 2014) while p300 is sourced from ref (Sciumè et al., 2020).

Figure 2 with 1 supplement see all
Klrc15′E and Klra7Hss1 are constitutively accessible, while promoters are accessible only at expressed alleles.

(A) FACS plot depicting splenic NK cells from a (B6 x BALB/c)F1 hybrid mouse stained with allele-specific antibodies, allowing separation of NK cells expressing both, either, or neither NKG2A allele. (B) (left) Normalized ATAC-seq data generated from the 4 cell populations depicted in (A) aligned to the mm10 reference genome. (right) Allele-informative reads were binned according to chromosome of origin, and displayed as signal mapping to the B6 or BALB/c chromosome. The Klrc15′E enhancer and promoter (Pro.) are boxed (dotted line). Vertical data range in SPMR is indicated for each track. (C and D) Data are as in (A and B), but using an allele-specific staining protocol with respect to the Ly49G2 receptor. Klra7Hss1, Klra7Hss5 and the dominant TSS (Pro3, Gays et al., 2011) are boxed. (E) CUT&RUN data depicting each of 4 indicated histone modifications at the Klra7 gene in IL-2 expanded NK cells sorted to express neither ‘N’ or both ‘B’ alleles of Klra7. Negative control CUT&RUN data were generated using a mouse IgG2aκ (cIgG) antibody, and a 50:50 mixture of IL-2 expanded NK cells that expressed the B6 or BALB/c alleles. These data are displayed in the bottom track in each sub-panel. The ATAC-seq patterns are shown for reference above each analysis; Klra7Hss1 is denoted as ‘1’, Klra7Hss5 is denoted as ‘5’. Arrows depict the locations of the dominant Pro3 TSS. All ATAC-seq and CUT&RUN data within a sub-panel are presented with the same vertical scale.

The KlraHss1 elements were hypothesized to serve as upstream bidirectional promoters active only in immature, developing NK cells (which do not yet express Ly49s); it was proposed that the direction of transcription predetermines the subsequent on or off state of the gene in mature NK cells, which is driven by a distinct downstream promoter (Saleh et al., 2004). Recent analysis of Klra gene expression in cell lines suggested instead that the KlraHss1 elements are transcriptional enhancers (Gays et al., 2015), and that the Hss1 transcripts represent enhancer RNAs (eRNAs). Those conclusions were in turn contested in a subsequent paper (McCullen et al., 2016). To address this issue in vivo, we analyzed chromatin features associated with the Hss1 elements in primary NK cells with published ChIP-seq data generated in mature primary splenic NK cells, using the H3K4me1:me3 ratio as an indicator of regulatory element identity (Calo and Wysocka, 2013). The Hss1 and 5′E elements are all enriched in H3K4me1 relative to H3K4me3 (Figure 1A), indicating enhancer identity. The putative NK receptor gene enhancers all ranked in the top 32% of ATAC-seq accessible peaks with respect to the H3K4me1:me3 ratio. In contrast, known promoters of the respective genes ranked in the bottom 21% (Figure 1C). In a deeper analysis, we independently defined enhancers and promoters in mature NK cells. NK cell promoters were defined as previously annotated mouse promoters from the EDPNew database (Dreos et al., 2017) enriched in H3K27ac in NK cells, and enhancers were defined as ATAC-seq peaks bound by the p300 histone acetyltransferase that do not overlap with the promoter list. Enhancers defined in this manner were highly skewed to high H3K4me1:me3 ratios, and promoters to low ratios (Figure 1—figure supplement 1B). All Hss1 and 5′E elements were classified as enhancers based on the p300-bound enhancer dataset (Figure 1—figure supplement 1B). These findings support the conclusion that the Hss1 and 5′E elements elements are enhancers in primary mature NK cells.

To test whether the Klra7Hss1 (Klra7 encodes Ly49G2) and Klrc15′E enhancers are required in mature NK cells, we adapted a CRISPR/Cas9 nucleofection protocol developed to edit primary human T cells (Roth et al., 2018; Figure 1—figure supplement 2). We used NK cells from (B6 x BALB/c)F1 hybrid mice and sorted NKG2AB6+ or Ly49G2B6+ cells using B6-allele reactive monoclonal antibodies against each receptor (Tanamachi et al., 2001; Vance et al., 2002) in order to follow the fate of a single allele in each case (Figure 1—figure supplement 2A). Editing efficiencies of NK cells were lower than that of T cells, resulting in only 30% or fewer cells with disruption of the control Ptprc locus encoding CD45 (Figure 1—figure supplement 2B). Targeting Klrc15′E increased the percentage of NKG2AB6-negative cells from ~10% to ~20%–40%. (Figure 1—figure supplement 2C-E), in line with our theoretical maximum editing efficiency. Similarly, targeting Klra7Hss1 resulted in marked loss of Ly49G2B6 expression, with minimal (<5%) loss of expression in non-targeting or non-nucleofected (no zap) control conditions (Figure 1—figure supplement 2F-H). These data show that Klrc15′E and Klra7Hss1 play key roles in maintaining expression of active alleles in mature NK cells, arguing against the proposal that Hss1 elements are only required in immature NK cells (Saleh et al., 2004; Saleh et al., 2002). The results provide functional evidence that Hss1 functions as an enhancer in mature NK cells.

The Klra7Hss1 and Klrc15′E enhancers are constitutively accessible

Analysis of bulk NK cells did not reveal a correlation between the gene expression frequency of an NK receptor gene and the accessibility, TF occupancy, or H3K27ac modifications of Hss1 and 5′E enhancers (Figure 1—figure supplement 1A). This lack of concordance raised the possibility that these enhancers were similarly active and occupied by TFs upstream of both silent and active alleles, as has been observed for RME genes genome-wide in F1 clones (Levin-Klein et al., 2017; Xu et al., 2017). It has not previously been possible, however, to address this issue for an RME gene in freshly isolated ex vivo cell populations.

To purify populations of cells expressing different alleles of Klrc1, we stained (B6 x BALB/c)F1 NK cells with allele-specific antibodies (Vance et al., 2002), allowing us to sort and perform ATAC-seq on NK cell populations expressing all four configurations of alleles: expressing both alleles of Klrc1, only B6, only BALB/c, or neither (Figure 2A and B). While the cells expressing both alleles and those expressing only the B6 allele are closely juxtaposed in the cytometry plots, post sort analysis demonstrated that they could be efficiently separated (Figure 2—figure supplement 1A), consistent with previously published data where allelic expression was confirmed by mRNA analysis (Rogers et al., 2006). SNPsplit (Krueger and Andrews, 2016) analysis of reads demonstrated that the enhancer element Klrc15′E was accessible on both active and inactive alleles in all four populations, whereas the Klrc1 promoter was accessible only at active alleles (Figure 2B). We used a similar allele-specific staining protocol (Tanamachi et al., 2001) to sort and analyze cells expressing either, both or neither Ly49G2 allele (Figure 2C, D). The Klra7Hss1 enhancer was accessible on both active and inactive alleles in all four populations, whereas the dominant promoter Pro3 (Gays et al., 2011) was accessible only on the active allele (Figure 2D). Notably, the Klra7 gene harbors a second minor enhancer element, Klra7Hss5 (Figure 1C; Figure 1—figure supplement 1B), which was similarly accessible at all alleles (Figure 2D).

These data demonstrated that enhancers within the Klra and Klrc gene families behave similarly to those of other RME genes analyzed in F1 hybrid clones (Levin-Klein et al., 2017; Xu et al., 2017), exhibiting an accessible configuration whether or not the gene was expressed. Importantly, this analysis further validated the NK receptor genes as a model for RME. While initially surprising, the decoupling of enhancer and promoter accessibility seen at NK receptor genes and other RME loci is consistent with a binary model of enhancer action, where enhancer activation occurs in all cells of a given type and acts stochastically to raise the binary “on or off” probability of gene expression, rather than regulate the per-cell amount of expression (Blackwood and Kadonaga, 1998; Walters et al., 1995).

We extended these observations by analyzing the pattern of active enhancer associated marks at silent and active alleles of Klra7 (which encodes Ly49G2). In order to obtain sufficient cell numbers of rare populations, we expanded NK cells from (B6 x BALB)F1 mice with IL-2 containing medium and then sorted cells that expressed neither (N) or both (B) Ly49G2 alleles and performed CUT&RUN for the enhancer-associated H3K4me1/2/3 and H3K27ac modifications. The Klra7 promoter and gene body displayed striking enrichment of H3K4me2/3 and H3K27ac in cells that expressed both Klra7 alleles, and as predicted lacked these modifications in cells where neither allele was expressed (Figure 2E). Notably, the Klra7Hss1 and Klra7Hss5 enhancers displayed equal enrichment of H3K27ac in cells expressing both alleles or neither (Figure 2E). As H3K27ac delineates active as opposed to poised enhancers (Calo and Wysocka, 2013), these data suggest constitutive enhancer activation on both silent and active alleles. These data are consistent with previous results demonstrating that enhancer-derived transcripts are produced at Hss1 elements in cells that do not express the target Klra gene (Gays et al., 2015). Thus, whether measured by accessibility, H3K27ac enrichment or eRNA production, Hss1 activity may be decoupled from target Klra expression status.

Klra1Hss1 and Klrc15′E are required for gene expression in vivo, and act in cis

We tested the requirement for Klra1Hss1 in the endogenous locus by deleting it in the B6 germline via CRISPR/Cas9 editing (Figure 3A-C, Figure 3—figure supplement 1A, B). Klra1Hss1Δ/Hss1Δ mice completely lacked Ly49A expression (Figure 3B, C; Figure 3—figure supplement 1B), but importantly expression of other Ly49 receptors was unaffected (Figure 3—figure supplement 1B), supporting the notion that the variegated NK receptor genes are regulated proximally and independently of each other. Klra1+/Hss1Δ mice displayed an intermediate percentage of Ly49A+ cells (Figure 3B and C; Figure 3—figure supplement 1B).

Figure 3 with 2 supplements see all
The Klra1Hss1 and Klrc15′E enhancers are required for gene expression.

(A) Locations of sgRNAs used to delete Klra1Hss1 in the B6 germline; NK cell ATAC-seq are displayed for reference with the vertical data scale in SPMR indicated. (B–C) Ly49A staining of the indicated Klra1Hss1 deletion littermates. MFI of staining ± SEM are depicted in gray. In (C), data are combined from two independent experiments (means ± SEMs, n = 5–12). (D–F) Data as in (A–C) for Klrc15′E. Data in (F) are combined from two experiments with the Klrc15′E(B3Δ) allele (Figure 3—figure supplement 1C) and were recapitulated in analysis of the B1Δ allele (means ± SEMs, n = 6–18). ****p < 0.0001; ***p < 0.0001 using one-way ANOVAs with Tukey’s multiple comparisons.

As with Klra1Hss1, deletion of both allelic copies of Klrc15′E in the germline eliminated NKG2A expression, and heterozygous mice displayed a reduced frequency of NKG2A+ NK cells (Figure 3D-F; Figure 3—figure supplement 1C, D).

Whether the activity of constitutively accessible enhancers of RME genes is coordinated in trans via a dedicated epigenetic mechanism that enacts RME is not known. We addressed the cis vs trans activity of Klrc15′E and Klra1Hss1 in F1 hybrids using allele-discriminating antibodies. F1 hybrid mice between B6-Klrc15′EΔ mice and BALB/c mice were generated (Klrc1B6-5′EΔ/BALB/c+ heterozygotes) along with WT F1 littermates. The percentages of cells expressing all four combinations of NKG2A alleles were determined using allele-specific NKG2A antibodies we previously generated. Assuming that the Klrc15′E acts only in cis we calculated the expected changes in the frequencies of these cells in the heterozygotes. The experimental data closely mirrored the predictions (Figure 3—figure supplement 2A-C). Therefore, the constitutively accessible Klrc15′E acted in cis and independently of the activity of the other copy. Similarly, in Klra1B6-Hss1Δ/BALB/c+ heterozygotes, the BALB/c allele was unaffected when expression of the B6 allele was abrogated (Figure 3—figure supplement 2D-F).

A cis-acting secondary enhancer in the Klra7 locus contributes to the high expression frequency of Klra7

Both the Klra1 and Klrc1 gene loci harbor only a single prominent proximal enhancer-like site (Figure 1A and B), and are completely dependent on those enhancers for expression (Figure 3), complicating analysis of the role of Klra1Hss1 and Klrc15′E in the RME of their target genes. We reasoned that analysis of an RME NK receptor gene with multiple proximal enhancers could reveal the role of overall enhancer strength in regulating expression frequency. We hypothesized, in accordance with the binary model, that despite the presence of multiple (weak) enhancers, enhancer activity at such loci is limiting resulting in RME. We predicted that limiting enhancer activity further by deleting a secondary enhancer in a natural RME gene would reduce, but not abrogate, gene expression probability.

The Klra7 locus is expressed by ~50% of NK cells and contains both an Hss1 element and another constitutively accessible enhancer, Klra7Hss5 (Figure 2D). Interestingly, the corresponding region of the highly related Klra1 gene, which is expressed by only ~17% of NK cells, is much less accessible and presumably less active (Figure 4A).

Figure 4 with 1 supplement see all
A minor cis-acting enhancer amplifies Ly49G2 expression frequency.

(A) Normalized ATAC-seq tracks of Klra1 and Klra7 in bulk NK cells; the vertical data range in SPMR is displayed on each track. Hss1 and Hss5 enhancers and the Pro3 promoter are highlighted. sgRNAs used to generate Klra7Hss5Δ alleles are shown (arrows). (B) Ly49G2 staining of NK cells in the indicated Klra7Hss5 deletion littermates (B2Δ allele, Figure 4—figure supplement 1A). (C–D) Ly49G2 percentages (C), and mean fluorescence intensities of the positive populations (D) (n = 4–11). Similar results were obtained with the B1Δ allele (Figure 4—figure supplement 1B). (*p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001 using One-way ANOVA with Tukey’s multiple comparisons). (E) Flow cytometry plots of gated Klra7B6-Hss5Δ/BALB/c+ NK cells using Ly49G2B6-specific and Ly49G2B6+BALB/c-specific antibodies (right) and a wildtype littermate (left). (F) Expected and observed percentages of populations depicted in ‘E’ in F1 mice with the Klra7Hss5Δ (hatched bar is expected, white bar is observed) or wildtype (black) Klra7 allele. Expected frequencies were calculated assuming stochastic cis regulation of alleles (see Materials and methods; note effect of genetic background). Data are representative of two experiments. All error bars represent SEM.

Germline deletion of Klra7Hss5, in homozygous configuration, resulted in a depressed percentage of Ly49G2+ cells (35%) compared to WT mice (50%), with only a minor change in expression per cell (measured by mean fluorescence intensity of staining, MFI) (Figure 4B-D; Figure 4—figure supplement 1A and B). Heterozygous mice displayed an intermediate percentage of Ly49G2+ cells (Figure 4B, C; Figure 4—figure supplement 1B). To test whether Klra7Hss5 acts entirely in cis, we crossed Klra7Hss5Δ to BALB/c mice. The NK cell populations in the F1 mice expressing the Ly49G2B6 alleles were reduced, and the populations expressing neither allele or only Ly49G2BALB/c were increased, in the proportion expected under probabilistic action of Klra7Hss5 in cis (Figure 4E, F). Thus, the constitutively active enhancer Klra7Hss5 is directly involved in regulating Ly49G2 expression frequency, and explains, at least in part, the high expression frequency of Ly49G2 in relation to other receptors including Ly49A.

Deletion of Klrk15′E is sufficient to recapitulate stable RME in Klrk1

Our hypothesis that RME of NK receptor genes is imparted by limiting binary enhancer activity predicts that a receptor expressed by all NK cells may be converted into a variegated receptor by weakening enhancer activity, for example by deleting one of multiple associated enhancer elements. We tested this for the Klrk1 gene encoding the NKG2D immunostimulatory receptor, which is expressed by all NK cells (Wensveen et al., 2018), is distantly related to the Klrc1 and Klra genes, and is flanked on both sides by enhancer-like chromatin, suggesting possible regulation by multiple enhancers (Figure 1B). Deletion of the enhancer-like ATAC-accessible site ~5 kb upstream of the Klrk1 gene (Klrk15′E) (Figure 1B; Figure 5—figure supplement 1A, B), resulted in variegated NKG2D expression in Klrk15′EΔ/5′EΔ animals (Figure 5A and B). Only ~65% of NK cells expressed NKG2D in Klrk15′EΔ/5′EΔ animals. The expression level per cell was only modestly affected and to an extent consistent with largely monoallelic expression vs expression in some cells of both alleles (Figure 5C and Figure 5—figure supplement 1C). These data suggest that the primary role of Klrk15′E is to increase the probability rather than the degree of Klrk1 expression, in line with the binary model of enhancer action.

Figure 5 with 2 supplements see all
Klrk15′E deletion results in mitotically stable RME, fully recapitulating natural variegation.

(A–C) NKG2D staining of splenocytes from Klrk15′E deletion littermates (B1Δ allele), and an Klrk1-/- mouse (n = 3–5). Results are representative of four experiments with two deletion alleles (Figure 4—figure supplement 1, C and D). (D) Splenocytes from Klrk15′EΔ/5′EΔ mice were cultured with IL-2 for 2–3 days before sorting NKG2D+ and NKG2D- NK cells, which were expanded in fresh IL-2 medium for 8–10 days before analysis (white fill). Expanded, unsorted NK cells are shown in gray. (E) Staining of splenic NK cells from mice of six genotypes. “+”, “-” and “Δ” refer to wildtype, gene knockout, and Klrk15′E deletion alleles, respectively. (F) Quantified results in (E) compiled from two experiments. (G) Expected and observed percentages of NKG2D+ NKcells in Klrk1-/5′EΔ mice. Expected expression is calculated based on observed NKG2D+ percentages in Klrk15′EΔ/5′EΔ mice, assuming stochastic expression (see Materials and methods). Data are comprised of selected groups displayed in (G). (H) Stochastic co-expression of NKG2D and NKG2A, Ly49I or Ly49G by NKp46+ NKcells in Klrk15′EΔ/5′EΔ mice. WT (+/+) mice are shown for comparison. (I) Expected (‘E’) and observed (‘O’) percentages of cells coexpressing the indicated receptors in Klrk15′EΔ/5′EΔ mice. Expected percentages were calculated by mutiplying percentages of cells in each mouse expressing each receptor individually (n = 4). Data are representative of two experiments. (J) NKG2D staining of presorted gated NK cells from Klrk15′EΔ/5′EΔ mice (bottom), compared to wildtype and Klrk1-/- NK cells. (K) Normalized ATAC-seq tracks generated from NKG2D+ and NKG2D- cells sorted from the Klrk15′EΔ/5′EΔ mouse shown in (J) and are presented on the same vertical scale. ATAC-seq results for WT splenic NK cells were sourced from Lara-Astiaso et al., 2014 and auto-scaled to match the data generated from the Klrk15′EΔ/5′EΔ mouse. Vertical data scale in SPMR is displayed for each track. Error bars represent SEM. **p < 0.01; ***p < 0.001; ****p < 0.0001, computed using One-way ANOVAs with Tukey’s multiple comparisons.

Significantly, the expression state of Klrk15′EΔ alleles was mitotically stable. NK cells from the enhancer knockouts were stimulated with IL-2 for 2–3 days before sorting NKG2D+ and NKG2D- populations, and further expanded in IL-2 for an additional 8–10 days, where they underwent an ~10–100 fold expansion. The NKG2D+ and NKG2D- phenotypes were highly stable despite extensive proliferation (Figure 5D).

In heterozygotes with the Klrk15′EΔ allele on one chromosome and a Klrk1 exon replacement knockout allele (Guerra et al., 2008) on the other (-/5′EΔ) the percentage of NKG2D+ cells was lower than in 5′EΔ/5′EΔ mice (Figure 5E and F). This nearly matched the expected percentage under the assumption that the 5′EΔ alleles are independently regulated in the heterozygotes, that is, the positive cells include cells expressing both alleles with a frequency that is the product of the individual frequencies (Raulet et al., 1997; Figure 5G). NKG2D expression per NKG2D+ cell in Klrk15′EΔ /5′EΔ animals appeared slightly higher than in Klrk1+/-, consistent with a proportion of cells expressing both Klrk1 alleles, a feature characteristic of natural RME (Eckersley-Maslin and Spector, 2014; Gendrel et al., 2016; Figure 5E and Figure 5—figure supplement 1C). Together, these data strongly argue that expression of Klrk15′EΔ alleles follows a stochastic RME pattern.

The Klrk15′EΔ allele mimics the expression and accessibility features of naturally variegated NK receptor genes

The stable RME of Klrk15′EΔ alleles recapitulated the stochastic expression pattern of naturally variegated NK receptor genes. Expression of NKG2D in Klrk15′EΔ/5′EΔ mice was approximately randomly distributed with respect to the expression of the RME genes encoding NKG2A, Ly49G2 or Ly49I (Figure 5H). Indeed, the co-expression frequencies approximated the products of the separate frequencies of the receptors studied (the ‘product rule’ Raulet et al., 1997; Figure 5I).

To examine the chromatin accessibility of the Klrk1 locus in Klrk15′EΔ/5′EΔ mice, we performed ATAC-seq with NKG2D+ and NKG2D- cells sorted from Klrk15′EΔ/5′EΔ mice (Figure 5J). Robust promoter accessibility was detected in NKG2D+ cells but not in NKG2D- cells (Figure 5K). The 5′E element was deleted in both populations and therefore not accessible, but a 3′ enhancer-like element was equally accessible in both populations. This accessibility pattern mirrors that of the naturally variegated NK receptor genes (Figure 2) and RME broadly (Xu et al., 2017). We speculate that this 3′ element works in concert with 5′E in WT cells to increase the expression frequency of Klrk1 in NK cells, and represents the residual enhancer that drives variegated Klrk1 expression after 5′E deletion. These speculations remain to be tested.

The results of our experiments with the Klrk1 gene established that stable RME and the stochastic and variegated NK receptor expression pattern could be recapitulated in full by weakening enhancer activity at a gene normally expressed in ~all NK cells. Furthermore, the similarity of Klrk15′EΔ/5′EΔ and natural NK receptor gene variegation suggests that previous examples of enhancer deletion-associated variegation such as that seen in the Cd8a locus (Ellmeier et al., 2002; Garefalaki et al., 2002) are rooted in similar mechanisms as naturally-occuring RME.

Silent Klrk1 alleles in Klrk15′EΔ/5′EΔ NK cells are CpG methylated at the promoter

Bisulfite conversion, PCR amplification, cloning and Sanger sequencing were carried out to examine methylation of the promoter region of Klrk1 in sorted NKG2D+ and NKG2D- NK cells from Klrk15′EΔ/5′EΔ mice, or similar NK cell populations that had been expanded in IL-2-containing medium for 8 days. Silent alleles in NKG2D-negative NK cells from Klrk15′EΔ/5′EΔ mice were hypermethylated at five CpG sites in an ~150 bp stretch surrounding the promoter, while these CpG sites were largely unmethylated in NKG2D+ NK cells, whether the latter cells came from Klrk15′EΔ/5′EΔ mice or WT mice (Figure 5—figure supplement 2A, B).

Correlations between expression and the absence of methylation of the corresponding promoters of the Klra and Klrc1 genes have been previously documented (Rouhi et al., 2007; Rouhi et al., 2006; Rouhi et al., 2009). The methylation of these promoters, or of Klrk1 promoters as documented here, is not necessarily a sufficient cause of silenced gene expression, however. For example, DNA methylation of promoters is not associated with silent alleles in the case of numerous other RME genes, nor is it necessary for maintaining the silenced state of the large majority of tested RME genes (da Rocha and Gendrel, 2019; Eckersley-Maslin et al., 2014; Gupta et al., 2022; Marion-Poll et al., 2021). In addition, inhibitors of DNA methyltransferases did not result in activation of silenced Klra genes or the Klrc1 gene (Rouhi et al., 2007; Rouhi et al., 2006; Rouhi et al., 2009). Treatments with 5-azacytidine proved to be toxic to primary NK cells, which precluded our ability to test the role of methylation in the maintenance of Klrk1 silencing.

A notable finding was that the same five CpG sites in the Klrk1 promoter were usually methylated in hematopoietic stem cells, which are progenitors of mature NK cells. Specifically, 13/16 of these sites were methylated based on analysis of whole genome bisulfite sequencing (WGBS) reads (Li et al., 2021). The promoter of the Klrc1 gene was also largely methylated in hematopoietic stem cells (Rogers et al., 2006), and we extended this analysis to the Klra1 promoter and Klra7 promoters. Using the dataset generated in Li et al., 2021, we found that nearly all the CpGs analyzed in 500 bp windows surrounding each promoter were methylated in hematopoietic stem cells (18 of 18 CpGs in reads mapping to Klra1, and 15/19 CpGs in reads mapping to Klra7). These findings suggest that activation of RME alleles is associated with demethylation, as opposed to a model where unmethylated promoters of RME genes are selectively methylated as a means to impose RME.

Silent NK receptor gene alleles lack repressive histone modifications associated with polycomb and heterochromatic repression

We investigated histone modifications associated with gene repression to search for clues regarding the maintenance of the active and silent epigenetic states. We assayed the polycomb-associated marks H3K27me3 and H2AK119Ub1 (H2AUb1) and the heterochromatin-associated H3K9me3, which have previously been found at inactive alleles of some other monoallelically expressed genes, notably the odorant receptors and protocadherins (Eckersley-Maslin and Spector, 2014; Eckersley-Maslin et al., 2014; Gendrel et al., 2014; Gendrel et al., 2016; Xu et al., 2017). CUT&RUN analysis of repressive modifications in IL-2 expanded primary NK cells that were sorted as Ly49G2-negative (expressing neither allele, designated ‘N’), revealed that all three modifications were prevalent in the Hoxa gene cluster, as expected, but the entire NKC lacked appreciable signal for any of the modifications (Figure 6A). None of the 3 marks were enriched above background on silent Klra7 alleles (Figure 6B). Importantly, many other genes associated with non-NK cell lineages (e.g. Cd19 and Mstn, expressed in B cells and myocytes, respectively) were also not enriched for these repressive modifications (Figure 6B). In contrast, other genes such as Pdcd1 (encodes PD-1) and Spi1 (encodes the macrophage and B-cell lineage-regulating transcription factor PU.1) displayed all 3 marks. Therefore, with respect to repressive chromatin marks, silent Klra7 alleles resembled several genes normally expressed in other hematopoietic cell lineages but not in NK cells, rather than known repressed genes.

Figure 6 with 1 supplement see all
Silent NK receptor gene alleles resemble inactive genes expressed in non-NK lineages, rather than repressed genes.

(A) Repressive histone modification CUT&RUN data generated with primary IL-2 expanded NK cells sorted to express neither allele of Ly49G2 (‘N’ cells). IGV screenshots depicting the indicated histone modification or analyses with control mouse IgG2aκ (cIgG), which binds protein A. The Hoxa gene cluster (left) serves as a positive control. The entire NKC gene cluster is displayed on the right. The vertical scales in SPMR, indicated on the left of the panels, were matched for each type of mark for all samples analyzed and were chosen to provide strong signals for the positive control Hoxa cluster. The cIgG data were scaled the same as the H3K27me3 data, which had the weakest signal of the marks analyzed. (B) Data are displayed as in (A), at Klra7 (left), and gene loci belonging to the following classes: NK cell lineage-appropriate, NK cell lineage non-specific, and loci repressed in NK cells.

As silent Klra7 alleles appeared similar to lineage non-specific genes in our analysis of repressive chromatin marks, we extended the analysis of chromatin states using ChromHMM, which integrates multiple datasets to classify the genome into subdomains based on their chromatin signatures (Ernst and Kellis, 2012). Using data from cells expressing neither or both Ly49G2 alleles, we constructed a minimal 3 state model corresponding to transcriptionally active chromatin (high levels of H3K27ac and H3K4me3, both active marks), repressed chromatin (H2AUb1 and H3K9me3, both repressive marks), and inactive chromatin (lacking these active or repressive marks) (Figure 6—figure supplement 1A). As expected, the promoters of lineage-appropriate genes expressed in NK cells (e.g. Ncr1, Klrb1c, Ifng) fell into the ‘active’ chromatin state 1 (Figure 6—figure supplement 1B). Notably, genes commonly regarded as markers of non-NK cell hematopoietic lineages (e.g. Cd3e, Cd19, Ly6g, Siglech) fell into the ‘inactive’ chromatin state 2 (Figure 6—figure supplement 1C). Finally, promoters of other genes, often encoding transcription factors that promote non-NK cell fates such as Bcl11b, Batf3, and Pax5, fell into the ‘repressed’ state 3 (Figure 6—figure supplement 1D). These data suggest that many genes encoding immune effector molecules associated with non-NK lineages are not actively repressed but are inactive and stably silent, whereas genes promoting non-NK cell fates are actively repressed.

In cells expressing both copies of Klra7, the enhancer, promoter and gene body all fell within the active state 1 (Figure 6—figure supplement 1E), whereas in cells expressing neither copy, the enhancer remained in the active state 1 but the promoter and gene body became inactive (state 2) rather than repressed (state 3). Indeed, it was striking that the NKC as a whole lacked repressive state 3 chromatin (Figure 6A; Figure 6—figure supplement 1E). The lack of the repressive chromatin state at silent NK receptor genes suggests that repressive chromatin may not be required for stable RME generally, potentially explaining why repressive chromatin signatures are not a consistent feature of silent RME alleles in other instances (Eckersley-Maslin and Spector, 2014; Eckersley-Maslin et al., 2014; Gendrel et al., 2014; Gendrel et al., 2016). In lieu of active repression, other mechanisms must be invoked for the maintenance of RME patterns through cell division.

The lineage-defining receptor genes Ptprc/Cd45, Klrk1, Cd8a and Thy1 are monoallelically expressed, suggesting RME may be ubiquitous

Our findings that RME is rooted in broad and probabilistic properties of gene activation raised the possibility that these principles might apply to many and perhaps all genes. We hypothesized that many genes exhibit a minor extent of RME but escaped previous detection due to the methods used, which were limited by clone numbers and therefore lacked the resolution to detect very rare monoallelic expression (Eckersley-Maslin et al., 2014; Gendrel et al., 2014; Gimelbrant et al., 2007; Reinius et al., 2016). We employed flow cytometry to analyze millions of primary cells ex vivo for rare monoallelic expression of several membrane proteins, starting with NKG2D. We noticed that ~2.5% of NK cells in Klrk1+/- mice lacked expression of NKG2D, whereas the percentage in Klrk1+/+ mice was close to 0% (Figure 5F; Figure 7A, B). These data suggested that rare monoallelic expression of WT Klrk1 was obscured by expression of at least one allele in nearly all NK cells. Indeed, assuming allelic independence, the failure of each allele to be expressed in a random 2.5% of all NK cells (from here on referred to as an ‘allelic failure rate’) translated to only 0.063% of cells lacking both alleles. Importantly, in both Klrk15′EΔ/5′EΔ and Klrk1+/- mice, NKG2D-negative cells were as likely as NKG2D+ cells to express NKG2A, Ly49G2, or Ly49I (Figure 7—figure supplement 1A-C). These data suggest that in both 5’E knockout and WT mice, NKG2D is randomly expressed with respect to the other RME receptor genes. Thus, even the WT Klrk1 gene is expressed in an RME fashion.

Figure 7 with 2 supplements see all
The lineage-defining Klrk1, Ptprc, Cd8a and Thy1 genes are RME genes.

(A, B) Flow cytometry (A) and quantification of % NKG2D-negative cells (B) of selected Klrk1 genotypes. p = 0.0021, student’s t-test. (C) Monoallelic CD45 expression. Flow cytometry of gated Thy-1+ cells pooled from 2 Ptprca/b mice (left). The mean percentages ± SEM of each monoallelic population, combined from three experiments, are depicted within the plot. Right panel: a mixture of cells from Ptprca/a and Ptprcb/b mice. (D) CD45 allele single positive and double positive T cell populations were sorted from Ptprca/b mice using gates in panel C, expanded for 1 week in vitro, resorted to purity and expanded an additional ~5–8 fold. Histograms show CD45.1 and CD45.2 staining for the sorted populations after expansion. (E–F) Monoallelic expression of CD8α in (B6 x CBA)F1 mice presented as in (C) and (D). (F) shows CD8β+ cells from F1 mice sorted and expanded twice as in (D) (~5–10 fold expansion in the second stumulation). (G–H) Data are displayed as in (C–F) but with respect to Thy-1 allelic expression on CD3+CD4+ T cells in (B6 x AKR)F1 hybrid mice (6–8 fold expansion in the second stimulation). All experiments are representative of 2–3 performed. Error bars represent SEM.

To extend this approach, we sought to analyze allelic expression of ‘lineage-defining’ receptor genes for which (A) allele-specific antibodies exist and, (B) expression is normally considered to be universal in defined lymphohematopoietic lineages. We first analyzed expression of allelic variants of the Ptprc gene encoding the membrane phosphatase CD45, expression of which defines all lymphohematopoietic cells. The Ptprca allele, encoding CD45.1, and the Ptprcb allele, encoding CD45.2, are easily discriminated by flow cytometry with monoclonal antibodies in congenic mice. Heterozygous B6-Ptprca/b mice are expected to express both alleles on all B or T cells, but we were able to detect clearly defined, albeit very rare (~0.01%) subpopulations of B or T cells expressing only one allele or the other (Figure 7C; Figure 7—figure supplement 2A, B). The monoallelic cells exhibited similar staining intensity as homozygous B6-Ptprca/a and B6-Ptprbcb/b cells analyzed in parallel. Given the very low frequency of monoallelic Ptprc expression, cells lacking CD45 altogether would be predicted to be extremely rare and indeed were not detected. Sorted CD45.1 and CD45.2 single positive T cells from B6-Ptprca/b mice retained monoallelic expression over five- to eightfold expansion after stimulation with anti-CD3/CD28 beads in vitro, demonstrating that RME of Ptprc is mitotically stable (Figure 7D). Sanger sequencing of reverse-transcribed and amplified RNA isolated from the expanded cells displayed in Figure 7D revealed that the monoallelic populations expressed only the allele detected by cell surface staining, demonstrating that the rare observed RME of Ptprc reflected transcriptional differences (Figure 7—figure supplement 2C). These data argued strongly against the possibility that the monoallelic cells arose due to somatic mutations in one or the other Ptprc allele, since most such mutations would not be predicted to disrupt transcription.

Expression of Cd8a defines the cytotoxic lineage of T cells. Similar analysis of Cd8a for RME employed allele-specific CD8α antibodies. For Cd8a, we gated on cells expressing CD8β, the partner chain of CD8α. In (B6 x CBA)F1 mice, approximately 0.1% of CD8β+ cells lacked one or the other of the CD8α alleles, CD8.1 or CD8.2 (Figure 7E; Figure 7—figure supplement 2D, E). Sorted single positive cells that retained expression of CD8β after stimulation and expansion retained expression of the initially selected allele of CD8α, demonstrating mitotically stable RME (Figure 7F). Finally, analysis of Thy1, a marker of T cells in mice, also employed allele-specific antibodies to discriminate the allelic Thy-1.1 and Thy-1.2 proteins, expressed by AKR and B6 strains, respectively. In (B6 x AKR)F1 hybrids ~5% of CD4+ T cells expressed only one allele or the other (Figure 7G; Figure 7—figure supplement 2F, G). Again, the monoallelic populations displayed an impressive degree of mitotic stability (Figure 7H). Hence, Thy1 represents an RME gene with a remarkably high allelic failure rate. Notably, the Thy1 promoter contains a CpG island, indicating that the RME we identifiy in lineage-defining genes is not restricted to promoters with low CpG content. In conclusion, RME was detectable for all four genes we examined, all of which were previously considered to be expressed in all cells of the lineages analyzed.

These findings supported the notion that RME is characteristic of many genes and is a natural consequence of enhancer-promoter interactions, rather than a specialized form of gene expression. We quantified allelic failure rates of the various genes examined in this study (Figure 7—figure supplement 2H). We propose that in the absence of selection for biallelic expression, most genes exist along a continuum of allelic failure rates, and that RME and ‘non-RME’ genes differ quantitatively with respect to allelic failure rates rather than qualitatively with respect to dedicated, RME-specific regulatory programs.

Discussion

The Klra1Hss1 element was previously reported to be a ‘switch’ element active only in immature, Ly49-negative NK cells (Saleh et al., 2004). Our extensive analysis herein demonstrated that Hss1 displays properties of enhancers in mature cells, consistent with the conclusions of others based on reporter analysis (Gays et al., 2015). Furthermore, the loss of Ly49G2 expression after deletion of Klra7Hss1 that we have documented in mature Ly49G2+ cells is inconsistent with a solely developmental role of Hss1, and further supports its enhancer identity in mature NK cells. Finally, our results showed that variegation arises and is modulated by enhancer deletion (including Klra7Hss5 and Klrk15′E) rather than introduction of variegating switch elements.

The main significance of our results is to link RME with the binary model of enhancer action, placing previous observations of enhancer deletion-associated variegation in the context of the pervasive and naturally-occurring RME phenomenon. Remarkably, deletion of an enhancer upstream of the Klrk1 gene imparted an RME expression pattern that fully recapitulated the stochasticity, mitotic stability and promoter accessibility features of naturally variegated NK receptor genes. In this instance, enhancer-like elements downstream of Klrk1 may suffice to impart the lower frequency of expression. The striking commonalities of enhancer deletion variegation and natural variegation of NK receptor genes and other RME genes argues that RME is an extreme manifestation of the inherent probablistic nature of stable gene activation rather than a specialized mechanism to impose a variegated expression pattern.

The data also reveal the quantitative impact of enhancer strength on allelic expression frequencies. Deletion of Klra7Hss5, a relatively minor and constitutively accessible enhancer, reduced the frequency of expression of Klra7, a natural RME gene, directly tying the enhancer deletion-associated variegation phenomenon to RME. This result powerfully argues that enhancers are not simply permissive for expression at RME genes, but are also instructive regarding expression probability. We hypothesize that the broad range of frequencies with which different Klra genes are naturally expressed (~5–60%) in large part reflects differences in enhancer strength. Probabilistic enhancer action has also been documented in Drosophila, where regulation of the gap genes by multiple enhancers has been suggested to ensure a high probability of gene expression, and removal of one of multiple enhancers was shown to increase gene ‘failure rate’ (Perry et al., 2011). We propose that the binary decision to express a gene is regulated by quantitatively varying enhancer activity, which is comprised of both the number of enhancers acting upon a gene and the strength of individual enhancers within that set. Our results are consistent with recent findings that enhancers are probabilistic regulators of transcription burst frequency rather than burst size (Bartman et al., 2016; Larsson et al., 2019). How enhancer control of the probability of stable gene expression interfaces with the control of transcription burst frequency is an exciting area for future investigation.

Our data showing RME of the lineage-defining Klrk1, Ptprc, Cd8a, and Thy1 genes support the generality of probabilistic gene expression, and suggest that RME is even more prevalent than the previous estimates (Eckersley-Maslin et al., 2014; Gendrel et al., 2014; Nag et al., 2013; Reinius et al., 2016). Although RME has been associated previously with poorly expressed genes (Gendrel et al., 2014; Reinius et al., 2016), our results extend the phenomenon to relatively highly expressed genes. We propose that genes lie along a spectrum of allelic failure rates that are largely controlled by enhancer strength, with documented RME genes on the highest end of that spectrum. RME applies to genes encoding cell surface receptors (Klra, Klrc1, Klrk1 etc.) and to a gene encoding a transcription factor, Bc11b (Ng et al., 2018). It applies as well to CpG poor promoters (the aforementioned receptor genes) and to genes that harbor promoter proximal CpG islands (Thy1 and the gene encoding Bcl11b Ng et al., 2018). High-resolution genome-wide approaches will eventually provide a comprehensive picture of the full extent of RME.

We predict that the mechanism of mitotic stability of active and silent RME alleles is likely related to maintenance of gene expression states broadly, and may involve bookmarking of promoters (Xu et al., 2017), rather than repressive chromatin at silent alleles. Indeed, we found that expressed and silent Klra7 alleles are distinguished only by the presence of active marks at the promoters and gene bodies of active alleles. Enhancers are constitutively accessible and activated on expressed and silent alleles alike, consistent with the decoupling of enhancer and promoter accessibility seen at RME loci generally. The absence of a repressive chromatin state on silent NK receptor alleles suggests that the RME ‘off state’ can reflect a stably inactive, as opposed to repressed, chromatin state. Numerous lineage non-specific genes that are also silent in NK cells (e.g. Cd19, Cd3e) also lacked traditional repressive chromatin modifications, yet are generally not subject to subsequent activation after an initial failure to be activated. In the case of the NK receptor genes and RME genes broadly it appears that the inactive state is maintained in mature cells in spite of continued enhancer activation, suggesting silent promoters are no longer competent for activation—perhaps due to lack of critical promoter-activating pioneer factor activity or reduced nucleosome remodeling or ‘promoter opening’ capacity that is only present at sufficient levels for developmentally regulated gene activation during differentiation. The relevant biochemical activity is likely a property of general rather than gene-specific factors, given the apparent pervasiveness of RME even among lineage-defining genes.

The results, in concert with previous studies, provide strong evidence that the principles underlying enhancer deletion-associated variegation and RME are broadly applicable to the regulation of many if not all genes. A possible mechanistic explanation of the probabilistic nature of gene expression in RME is that the initially inactive promoters present an energetic threshold that must be overcome to stably activate gene expression, and that enhancer strength determines the probability of overcoming this threshold. Once overcome, the active state is largely stable even in the absense of the initial signal or enzymatic activity, which may be restricted to a cellular differentiation window. The energy barrier may be due to the initially chromatinized state of the promoter and the energy required for nucleosome remodeling of promoters. It may be this property that is the measure of ‘enhancer strength’, that is, the capacity of the enhancer to overcome the energy threshold presented by the promoter within a limited window during differentiation. The number of enhancer elements in a gene, the availability of relevant transcription factors that bind the enhancers and the affinity with which the enhancer binds the factors are all likely relevant in determining enhancer strength. It remains possible that in some cases of RME, additional modifications of the promoter may contribute to the energy barrier as well.

Differential DNA methylation of alleles is not broadly associated with all or even most RME genes (da Rocha and Gendrel, 2019; Eckersley-Maslin et al., 2014), nor does it appear to be responsible for maintenance of silent alleles in the majority of tested RME genes (Eckersley-Maslin and Spector, 2014; Gupta et al., 2022; Marion-Poll et al., 2021). Silent Klra and Klrc1 alleles in mouse NK cells were reported to be hypermethylated relative to active alleles (Rouhi et al., 2007; Rouhi et al., 2006; Rouhi et al., 2009), and our data show that silent Klrk1 alleles in 5′EΔ/5′EΔ NK cells are also hypermethylated. It is doubtful that DNA methylation is sufficient to maintain the silent state, however, since inhibitors of DNA methyltransferases failed to derepress inactive Klra or Klrc1 genes in cell lines (Rouhi et al., 2007; Rouhi et al., 2006; Rouhi et al., 2009), or other silent RME genes that were hypermethylated (Eckersley-Maslin and Spector, 2014; Gupta et al., 2022; Marion-Poll et al., 2021). In light of the collective data on DNA methylation and silent alleles in RME we suspect it is not generally causative in maintaining RME. Furthermore, the Klra, Klrc1, and Klrk1 genes are all hypermethylated on both alleles in progenitor hematopoietic stem cells, suggesting that demethylation is associated with gene activation as opposed to a model where alleles are randomly silenced by de novo DNA methylation. The possibility remains that for these genes, preexisting DNA methylation contributes to the energy barrier that must be overcome by enhancer action, and that demethylation accompanies de novo gene activation. DNA methylation is not a requisite feature of RME promoters, however, as many RME genes are hypomethylated on both alleles.

The RME of the NK receptor genes resembles the monoallelic expression pattern of cytokine genes including Il2, Il4, Il5, Il10, and Il13 (Bix and Locksley, 1998; Gendrel et al., 2016; Kelly and Locksley, 2000; Rivière et al., 1998). The cytokine genes are inducible in response to TCR stimulation and therefore expression is inherently unstable, but impressive stability over several mitotic divisions was observed for Il4 (Bix and Locksley, 1998), and the probability of Il4 allelic activation and biallelic expression correlated with the strength of the inducing signal (Rivière et al., 1998). Intriguingly, the Il4 and Il13 genes are closely linked and are co-regulated by an enhancer, CNS1, that was found to be constitutively acetylated at histone H3, and thus permissive for expression (Guo et al., 2005). Expression of the Il4 and Il13 gene alleles was independent, however, in a manner strikingly similar to the NK receptor genes. It is probable that the general principles uncovered by us and others apply to genes that are induced via stimulation as well as genes whose expression is acquired during differentiation. Importantly, enhancer deletions have been observed to reduce the stochastic expression probabilities of the monoallelically expressed odorant receptor genes (Khan et al., 2011) and trace amine-associated receptors of the olfactory epithelium (Fei et al., 2021). It therefore seems likely that the probabilistic nature of gene activation and binary model of enhancer-promotor communication accounts for the initial stochastic step across many or even all types, even if they might differ with respect to feedback and allelic exclusion mechanisms.

Apparently, RME often occurs at such a low rate that it is both beneath ready detection and presumably irrelevant for the function of a cell lineage. In the case of NK receptors, we propose that evolution has exploited RME to generate a complex combinatorial repertoire of NK cell specificities. More speculatively, by regulating expression of fate-determining mediators, RME may underlie stochastic cell fate decisions in some instances of cellular development (Ng et al., 2018). From an evolutionary perspective, appreciable RME of a gene could arise by mutation of strong enhancers of a precursor gene, by providing a new gene with a weak enhancer, or by diminishing the concentration of relevant enhancer-binding transcription factors in a given lineage of cells (as has been shown for several relevant TFs for the Klra genes) (Ioannidis et al., 2003; Ohno et al., 2008).

Finally, our results suggest that allelic failure rates may in some cases dwarf the rates of null alleles generated by somatic mutation. As a novel mechanism of genetic haploinsufficiency at the cellular level, RME might have broad implications in genetic disease etiology and penetrance of disease phenotypes in heterozygous individuals.

Materials and methods

Animals and animal procedures

Request a detailed protocol

All mice were maintained at the University of California, Berkeley. Klrk1-/- mice are available at the Jackson Laboratory (JAX Stock No. 022733). C57BL/6 J (B6), BALBc/J and B6;129-Ncr1tm1Oman/J (Ncrgfp) were purchased from the Jackson Laboratory and bred at UC Berkeley. BALB/cJ, CBA/J and AKR/J and B6.SJL-Ptprca Pepcb/BoyJ mice were purchased from the Jackson Laboratory. F1 hybrid mice were purchased from the Jackson Laboratory or were generated at UC Berkeley from inbred parents.

For the generation of CRISPR edited mice, Cas9 RNP was delivered to single-cell embryos either through microinjection or CRISPR-EZ electroporation, both of which are described in reference (Modzelewski et al., 2018). Klra1Hss1Δ mice were generated by microinjection, while Klrc15′EΔ, Klrk15′EΔ and Klra7Hss1Δ mice were generated by CRISPR-EZ electroporation. Whether through microinjection or electroporation, we used paired sgRNAs flanking the enhancer to generate enhancer deletion mice. sgRNAs were selected using the GPP web portal from the Broad Institute. Guides with highest predicted editing efficiencies were prioritized, while also minimizing for predicted off-target cutting in protein-coding genes. sgRNAs were generated using the HiScribe T7 Quick High Yield RNA Synthesis Kit (New England Biolabs). Founder mice (F0) harboring deletion alleles were backcrossed to C57BL/6 J (B6) mice to generate heterozygous F1 mice, and were then intercrossed to generate WT, heterozygous and homozygous littermates for experiments. All sgRNAs used for the generation of enhancer deletion mice are listed in Supplementary file 1. Primers used to PCR identify edited founders and genotype subsequent filial generations are listed in Supplementary file 1. All animals were used between 8 and 32 weeks of age, and all experiments were approved by the UC Berkeley Animal Care and Use Committee (ACUC).

Flow cytometry

Request a detailed protocol

Single-cell splenocyte suspensions were generated by passing spleens through a 40 μm filter. Red blood cells were lysed with ACK buffer. Fresh splenocytes, or where indicated cells cultured with 1000 U/ml recombinant human IL-2 (National Cancer Institute) were stained for flow cytometry in PBS containing 2.5% FCS (FACS Buffer). Before staining with antibodies, FcγRII/III receptors were blocked for 15 min at 4°C using 2.4G2 hybridoma supernatant. Cells were washed with FACS buffer and then stained with antibodies directly conjugated to fluorochromes or biotin at 4°C for 15–30 min. In order to differentiate between alleles of a receptor in (B6 x BALB/c)F1 hybrid NK cells, the B6-specfic clone was used first in order to block epitopes in competition with the clone recognizing both alleles. For example, to discriminate Ly49G2 alleles, cells were stained for at least 15 min with 3/25 which recognizes Ly49G2B6, and then 4D11 (which recognizes both alleles) was added. For discriminating alleles of NKG2A, cells were stained first with the NKG2AB6-specific 16a11, followed by 20d5, which binds to both alleles. Anti-Ly49AB6 (A1) was added before the non-discriminating JR9 mAb, but in this case, cells expressing only the B6 allele did not resolve from the population of cells expressing both alleles. When necessary, cells were washed and then stained with secondary antibody or fluorochrome-conjugated streptavidin. Near-IR viability dye (Invitrogen L34975) or DAPI (Biolegend 422801) were used to discriminate live cells. Flow cytometry was carried out using an LSR Fortessa or X20 from BD Biosciences, and data were analyzed using FlowJo software. In all cases, NK cells were defined as CD3-NKp46+ splenocytes. For sorting on a BD FACSAria II sorter, the samples were prepared nearly identically as they were for flow cytometric analysis with the exception that the medium used was sterile RPMI 1640 (ThermoFisher) with 5% FCS.

Antibodies used in flow cytometry

Request a detailed protocol

From Biolegend: anti-CD3ε (145–2C11) PE-Cy5, anti-CD19 (6D5) PE-Cy5, anti-F4/80 (BM8) PE-Cy5, anti-Ter119 (TER-119) PE-Cy5, anti-NKp46 (29A1.4) BV421, anti-NKG2AB6 (16a11) PE, anti-Ly49AB6 (A1) PE, anti-NKG2D (CX5) PE/Dazzle 594, anti-CD8β (YTS156.7.7) PE-Cy7, anti-CD45.1 (A20) APC, anti-CD45.2 (104) FITC, anti-CD90.2 (53–2.1) PE or FITC, goat-anti-mouse IgG (Poly4053) PE. From eBioscience/ThermoFisher: anti-NKG2A (20d5) PerCP, anti-Ly49I (YLI-90) FITC, anti-Ly49G2 (4D11) PerCP-eFlour 710 or PE-Cy7, anti-CD90.1 (HIS51) FITC, anti-rat IgG F(ab’)2 (polyclonal, lot 17-4822-82) APC. From BD Biosciences: anti-CD4 (GK1.5) BUV737. From BioXCell: anti-CD8.1 (116–13.1) unconjugated, anti-CD8.2 (2.43) unconjugated. Prepared in-house: anti-Ly49A (JR9) (Roland and Cazenave, 1992) biotin, anti-Ly49G2B6 (3/25) (Tanamachi et al., 2001), anti-NKG2D (MI-6) biotin. Biotin conjugated antibodies were used in conjunction with streptavidin (conjugated with APC or PE) from Biolegend.

Ex vivo NK cell cultures

Request a detailed protocol

NK cells were prepared from spleens by passage through a 40 μm filter. Red blood cells were lysed with ACK. Splenocytes were cultured in RPMI 1640 (ThermoFisher) with 1000 U/mL IL-2 (National Cancer Institute) and 5% FCS. In all cases, media was supplemented with 0.2 mg/mL glutamine (Sigma), 100 U/mL penicillin (ThermoFisher), 100 μg/mL streptomycin (Thermo Fisher Scientific), 10 μg/mL gentamycin sulfate (Fisher Scientific), 50 μM β-mercaptoethanol (EMD Biosciences), and 20 mM HEPES (ThermoFisher).

Analysis of the stability of monoallelic expression of NKG2D

Request a detailed protocol

NKG2D+/- NK cells were sorted from WT or Klrk15′EΔ/5′EΔ mice on day 2 or 3 of ex vivo NK cell culture in IL-2 medium as described above. Cells were cultured in vitro in IL-2 containing media for a further 8–10 days, during which cells expanded ~10–100 fold based on hemocytometer counts. Cells were analyzed for NKG2D expression by flow cytometry. In all cases medium contained 5% FCS (Omega Scientific), 0.2 mg/mL glutamine (Sigma), 100 U/mL penicillin (ThermoFisher), 100 μg/mL streptomycin (Thermo Fisher Scientific), 10 μg/mL gentamycin sulfate (Fisher Scientific), 50 μM β-mercaptoethanol (EMD Biosciences), and 20 mM HEPES (ThermoFisher). Cells were incubated at 37°C in 5% CO2.

Ex vivo assay for the stability of monoallelic expression in T cells

Request a detailed protocol

Cells from the spleens and a collection of lymph nodes (brachial, axial, inguinal, mesenteric) from F1 hybrid mice and parental inbred line controls were combined and passed through a 40 μm filter, and red blood cells were lysed with ACK buffer. Cells were prepared for sorting as described above, staining with the relevant allele-specific antibodies. For CD45 monoallelic expression, Thy1+ cells were further gated according to CD45 allelic expression. For CD8α monoallelic expression, CD3+CD8β+ cells were analyzed for CD8α allelic expression. For Thy1 monoallelic expression, CD4+MHC II- cells were analyzed for Thy-1 allelic expression. Cells expressing either the paternal or maternal allele (or both) of the receptor studied were sorted and expanded for 1 week in RPMI 1640 (ThermoFisher) containing 200 U/mL recombinant IL-2, Dynabeads mouse T-activator CD3/CD28 (ThermoFisher) beads at a 1:1 cells to beads ratio, 10% FCS, and supplemented as described for ex vivo NK cell cultures. After 1 week of expansion, cells were harvested, counted by hemocytometer and prepared for a second sort. After sorting for expression of the relevant receptor allele again in order to ensure purity, cells were once again expanded in a restimulation, this time with a cells to beads ratio of 10:1. After the second expansion, cells were again counted, stained and prepped for final analysis of monoallelic receptor expression by flow cytometry.

In analysis of Ptprc monoallelic expression, RNA was isolated from expanded T cells expression either or both Ptprc alleles as displayed in Figure 7D using the iScript cDNA synthesis kit (BioRad), from 10,000 to 40,000 cells. Half of the reaction volume (10 μL out of 20 μL) was used to PCR amplify a region of the Ptprc transcript using intron-spanning PCR primers (Supplementary file 1).

Enhancer deletion in primary NK cells via Cas9-RNP nucleofection

Request a detailed protocol

Ex vivo editing of primary mouse NK cells was carried out according to a modified version of the protocol used to modify primary human T cells described in reference (Roth et al., 2018). Cas9 was purchased from the UC Berkeley Macro Lab core (40 µM Cas9 in 20 mM HEPES-KOH, pH 7.5, 150 mM KCl, 10% glycerol, 1 mM DTT), and sgRNAs were transcribed in vitro according to the Corn lab online protocol (https://www.protocols.io/view/in-vitro-transcription-of-guide-rnas-and-5-triphos-bqjbmuin). NK cells were prepared by sorting day 5 IL-2 cultured NK cells from (B6 x BALB/c)F1 hybrids. CD3-NKp46+ cells were sorted to be positive for either NKG2AB6 using the 16a11 clone or Ly49G2B6 using the 3/25 clone, and cells were further cultured overnight in RPMI 1640 media containing 5% FCS and 1000 U/mL IL-2 (National Cancer Institute). On day 6, 1 million sorted NK cells were prepared for nucleofection using the Lonza 4D-Nucleofector per condition. Cas9 and sgRNAs were complexed at a molar ratio of 1:2 (2.5 μL of 40 μM Cas9 was added to 2.5 μL of sgRNA suspended at 80 μM (6.5 μg) in nuclease-free H2O). If two flanking guides were used, 1.25 μL of each were used, maintaining the Cas9 to sgRNA molar ratio. Cas9-RNP was complexed for 15 min at 37°C and transferred to a single well of a 96-well strip nucleofection cuvette from Lonza for use with the Nucleofector 4D. 1 million sorted day 6 IL-2 cultured NK cells were resuspended in 18 μL of supplemented Lonza P3 buffer from the P3 Primary Cell kit, and added to the Cas9-RNP complex. Cells were nucleofected using the CM137 nucleofection protocol and 80 μL pre-warmed RPMI 1640 with 5% FCS was immediately added. After a 15-min recovery period at 37°C, cells were returned to culture in 1 mL of RPMI 1640 with 5% FCS and 1000 U/mL IL-2. After 5–7 days in culture maintaining a maximum density of 1 million cells/mL, receptor expression was assayed by flow cytometry. In order to validate enhancer flanking guides (Supplementary file 1) an identical protocol was followed with either day 5 IL-2 cultured splenocytes, or day 5 IL-2 cultured NK cells isolated using the MojoSort NK isolation kit from Biolegend, but instead of analysis by flow cytometry, gDNA was prepared and used as a template for PCR to detect the expected deletion.

F1 hybrid genetics and calculations of expected changes in receptor-expressing NK cell populations

Request a detailed protocol

F1 hybrid genetics were carried out by breeding WT or CRISPR/Cas9-edited males on the B6 background to females from the following backgrounds: BALBc/J, CBA/J, AKR/J. Edited alleles were crossed only to BALBc/J, while CBA/J and AKR/J were used in the F1 hybrid analysis of monoallelic expression of CD8α and Thy-1, respectively.

We estimated the expected frequencies of NK cells in (Klrc1B6-5′EΔ/BALB/c+) F1 mice by assuming independence of allelic expression. That assumption leads to the following predictions:

The percentage of cells expressing neither allele in the mutant will equal the sum of the percentages of the two NK cell populations that lack NKG2ABALB/c in WT (B6 x BALB/c)F1 hybrids, that is the cells that express neither allele, and cells expressing only the B6 allele.

The percentage of cells expressing NKG2ABALB/c only in the mutant will equal the sum of the percentages of the NK cell populations in WT (B6 x BALB/c)F1 hybrids that express NKG2ABALB/c, that is the cells that express only the BALB/c allele, and the cells expressing both alleles.

The percentages of cells expressing NKG2AB6 only or both NKG2AB6 and NKG2ABALB/c will be 0, since NKG2AB6 is not expressed.

The expected frequency of cells expressing NKG2D in Klrk1-/5′EΔ mice was calculated assuming stochastic expression of alleles, and was based on the frequency of cells expressing NKG2D, or not, in Klrk15′EΔ/5′EΔ mice. The frequency of cells lacking expression of a given allele is the square root of the frequency of cells expressing neither allele. Subtraction of this proportion from 1 yields the predicted frequency of cells expressing NKG2D in Klrk1-/5′EΔ mice. For example, an observed NKG2D expression frequency of ~67% in a Klrk15′EΔ/5′EΔ mouse would result in an expected frequency datapoint of ~43%.

The expected changes in populations with respect to Ly49G2 alleles in Klra7B6-Hss5Δ/BALB/c+ mice were calculated with the same assumption of independent regulation of alleles.

We started by calculating the overall percentage of cells expressing Ly49G2B6 in the F1 mice with the mutation, which averaged 47.7% of that in WT F1 mice.

The predicted percentage of cells expressing only Ly49AB6 in the mutant F1 was then 47.7% of the percentage of cells expressing only Ly49AB6 in WT mice.

And the predicted percentage of cells expressing both alleles in the mutant F1 was 47.7% of the percentage of cells expressing both alleles in WT mice.

The predicted percentage of cells expressing neither allele in the mutant F1 was calculated as the percentage of cells expressing neither allele in WT mice plus 52.3% (100%–47.7%) of the precentage of NK cells that express only Ly49G2B6 in WT mice.

Finally, the predicted percentage of NK cells expressing only Ly49G2BALB/c in the mutant was calculated as the percentage expressing only Ly49G2BALB/c in WT mice plus 52.3% of the NK cells expressing both alleles in WT mice.

Note that the genetic background of the mice significantly influences Klra7 expression even in WT mice, presumably reflecting trans-acting events (e.g. each WT Klra7B6 allele is expressed on ~31% of NK cells in B6 mice, but only ~19% in F1 hybrid mice). Therefore expected data are calculated using Ly49G2B6 expression frequencies in Klra7B6-5′E+/BALB/c+ mice.

ATAC-Seq

Request a detailed protocol

ATAC-seq was performed as previously described in reference (Buenrostro et al., 2013). Briefly, 50,000 sorted NK cells were washed in cold PBS and resuspended in lysis buffer (10 mM Tris-HCl, pH 7.4; 10 mM NaCl; 3 mM MgCl2; 0.1% (v/v) Igepal CA-630). The crude nuclear prep was then centrifuged and resuspended in 1 x TD buffer containing the Tn5 transposase (Illumina FC-121–1030). The transposition reaction was incubated at 37°C for 30 min and immediately purified using the Qiagen MinElute kit. Libraries were PCR amplified using the Nextera complementary primers listed in reference (Buenrostro et al., 2013) and were sequenced using an Illumina NextSeq 500 or a HiSeq 4000.

CUT&RUN

Request a detailed protocol

CUT&RUN was performed essentially as previously described (Skene et al., 2018). Briefly, 50,000–500,000 NK cells sorted as indicated were washed and immobilized on Con A beads (Bangs Laboratories) and permeabilized with wash buffer containing 0.05% w/v Digitonin (Sigma-Aldrich). Cells were incubated rotating for 2 hr at 4°C with antibody at a concentration of 10–20 μg/mL. Permeabilized cells were washed and incubated rotating at room temperature for 10 min with pA-MNase (kindly provided by the Henikoff lab) at a concentration of 700 ng/mL. After washing, cells were incubated at 0°C and MNase digestion was initiated by addition of CaCl2 to 1.3 mM. After 30 min, the reaction was stopped by the addition of EDTA and EGTA. Chromatin fragments were released by incubation at 37°C for 10 min, purified by overnight proteinase K digestion at a concentration of 120 μg/mL with 0.1% wt/vol SDS at 55°C. DNA was finally purified by phenol/chloroform extraction followed by PEG-8000 precipitation (final concentration of 15% wt/vol) using Sera-mag SpeedBeads (Fisher) (https://ethanomics.files.wordpress.com/2012/08/serapure_v2-2.pdf).

Libraries were prepared using the New England Biolabs Ultra II DNA library prep kit for Illumina as described online (https://www.protocols.io/view/library-prep-for-cut-amp-run-with-nebnext-ultra-ii-bagaibse?version_warning=no) with the following specifications and modifications. The entire preparation of purified CUT&RUN fragments from a reaction were used to create libraries. For histone modifications, end repair and dA-tailing were carried out at 65°C. NEB hairpin adapters (From NEBNext Multiplex Oligos for Illumina) were diluted 25-fold in TBS buffer and ligated at 20°C for 15 min, and hairpins were cleaved by the addition of USER enzyme. Size selection was performed with AmpureXP beads (Agencourt), adding 0.4 X volumes to remove large fragments. The supernatant was recovered, and a further 0.6 X volumes of AmpureXP beads were added along with 0.6 X volumes of PEG-8000 (20% wt/vol PEG-8000, 2.5 M NaCl) for quantitative recovery of smaller fragments. Adapter-ligated libraries were amplified for 15 cycles using NEBNext Ultra II Q5 Master Mix using the universal primer and an indexing primer provided with the NEBNext oligos. Amplified libraries were further purified with the addition of 1.0 X volumes of AmpureXP beads to remove adapter dimer and eluted in 25 μL H2O. Libraries were quantified by Qubit (ThermoFisher) and Bioanalyzer (Agilent) before sequencing on an Illumina HiSeq 4,000 or MiniSeq as paired-ends to a depth of 10–32 million.

The following antibodies were used for CUT&RUN: Abcam: anti-H3K4me1(ab8895), anti-H3K4me2 (ab7766), anti-H3K4me3 (ab8580), anti-H3K27ac (ab4729), anti-H3K9me3 (ab8898). Cell Signaling: anti-H3K27me3 (C36B11), anti-H2AUb1 (D27C4). Control IgG (cIgG) from Biolegend: Mouse IgG2aκ (MOPC-173).

Bisufite conversion and analysis of promoter CpG methylation

Request a detailed protocol

Genomic DNA was extracted from 50,000 to 100,000 FACS-sorted cells by addition of 10 μl of water, 10 μl PBS, 5 μl proteinase K, and 15 μl lysis buffer. The lysate was incubated for 30 min at 58°C. Bisulfite conversion was performed using the EpiTect Fast LyseAll Bisulfite Kit (Qiagen) according to the manufacturer’s protocol. Cleanup was also performed according to the manufacture’s protocol using MinElute columns. Bisulfite-treated DNA was amplified by PCR using EpiTaq HS (Takara) using the primers targeting a region spanning 481bp and 5 CpG sites in the the Klrk1 promoter indicated in Figure 5—figure supplement 2A: (Forward: 5’- ATAGAGATAGTAGAAAAAAAATTTGTTAGAAT – 3’; Reverse: 5’-AAAACTTTCCACAATCTCAAAAACTAAATT - 3’). 35 amplification cycles were performed: 10s at 98°C, 30s at 55°C, and 60s at 72°C. PCR products were cloned into the TOPO TA plasmid from Invitrogen, and transformed into XL-1 competent cells (UC Berkeley Macro Lab). Clones were selected by blue/white screening, used for colony PCR, and analyzed by Sanger sequencing of the colony PCR product.

We further analyzed the CpG methylation status of the Klrk1 promoter in hematopoietic stem cells (HSCs) by mining published whole genome bisulfite sequencing (WGBS) data from Li et al., 2021. We downloaded .cov files from the associated GEO series (GSE167237), and scored the frequency of methylated CpGs at the same five sites that we analyzed in the bisulfite analysis described above. This series contains 4 WGBS datasets generated using E14 fetal liver HSCs and 2 generated using adult bone marrow HSCs; datapoints were compiled across all six datasets. We performed a similar analysis at the Klra1 and Klra7 promoters in a region spanning 500 bp around each promoter.

Datasets and processing and visualization

Request a detailed protocol

Raw mined datasets were downloaded from NCBI Gene Expression Omnibus (GEO) or the European Bioinformatics Institute (EBI). NK cell ATAC-seq and histone modification (H3K4me1, H3K4me2, H3K4me3, H3K27ac) were from reference (Lara-Astiaso et al., 2014) under GEO accession numbers GSE59992 and GSE60103. Runx3 ChIP-seq data and non-immune serum control in NK cells were sourced from reference (Levanon et al., 2014) (GSE52625) and T-bet ChIP-seq data and input control were sourced from reference (Shih et al., 2016) (GSE77695). p300 ChIP-seq raw data was sourced from reference (Sciumè et al., 2020) (GSE145299). p300 ChIP-seq peaks were called in reference (Sciumè et al., 2020) and downloaded in .csv format.

Raw data from all datasets (mined or generated in this study) were processed using a pipeline assembled in-house. Datasets were tested with FastQC. Paired-end reads were then aligned to the mm10 reference genome using Bowtie2 with the --sensitive parameter. Paired-end CUT&RUN libraries were tested and aligned with the same pipeline. All reads aligned to the mitochondrial chromosome were removed with samtools. Aligned reads were then sorted, indexed, and filtered for a mapping quality of ≥10 with samtools. PCR duplicates were removed with Picard (Broad Institute). Reads covering blacklisted regions (ENCODE mm10 database) were removed with bedtools. Data were then normalized to signal per million reads (SPMR) when calling narrow peaks with MACS2. Resultant bedgraph files were converted to bigwigs with the bedGraphToBigWig program from the UCSC Genome Browser toolkit for visualization on Integrative Genomics Viewer (IGV) (Thorvaldsdóttir et al., 2013). Data in Figure 1—figure supplement 1A were plotted using the Bioconductor package SeqPlots (Stempor and Ahringer, 2016).

Ranking of accessible sites in NK cells according to H3K4me1:me3 ratio

Request a detailed protocol

Reads from duplicate ChIP-seq datasets (for both H3K4me1 and H3K4me3) from reference (Lara-Astiaso et al., 2014) were merged to ensure robust signal, and the resultant files were processed and normalized as above. NK cell ATAC-seq peaks were called in the Ly49G2B6+BALB+ NK cell ATAC-seq dataset using macs2 narrowpeaks. Before ranking, ATAC-seq peaks were filtered such that only peaks that fell within the top 95% of both H3K4me1 and H3K4me3 signal computed over a 2 kb window from the peak midpoint using pandas and numpy in Python 3.7.4, resulting in 51,650 usable peaks. H3K4me1:me3 raw ratio and log2 ratio bigwigs were generated with the bamCompare utility from deepTools (v2.5.4). The log2 ratio track was visualized on IGV, and the raw ratio was used to rank ATAC-seq peaks. Heatmaps were generated with the computeMatrix and plotHeatmap utilities from deepTools (v2.5.4). Heatmaps were sorted by the mean H3K4me1:me3 ratio signal over a 2 kb window centered at the midpoint of the 51,650 ATAC-seq peaks. Hss1 and 5′E enhancer regions and corresponding promoters at NKC genes were individually predefined and the position of each was then marked on the heatmap.

Definition of NK cell promoters and enhancers and ranking of regulatory elements according to H3K4me1:me3 ratio

Request a detailed protocol

Annotated mouse promoters (defined as the TSS at a single nucleotide) in the mm10 genome assembly were downloaded as a BED file from the EDPNew database (Dreos et al., 2017). To identify likely active promoters in NK cells, broad regions of H3K27ac were called based on ChIP-seq data sourced from reference (Lara-Astiaso et al., 2014) using the “macs2 callpeak --broad” command. Mouse EDPNew promoters falling within broad H3K27ac domains were identified using the “bedtools intersect -wa” command, resulting in a set of 9901 active promoters in mouse NK cells (Source data 1-mouse NK cell promoters).

Enhancers in naïve mouse NK cells were defined as the intersection of ATAC-seq and p300 peaks not found at the promoters as defined above. p300 ChIP-seq peaks in resting NK cells were previously defined and downloaded from reference (Sciumè et al., 2020). ATAC-seq peaks that were enriched in p300 binding were identified using the “bedtools intersect -wa” command. To define enhancers that do not overlap annotated promoters, EDPNew promoters were subtracted from p300-enriched ATAC-seq peaks using the “bedtools subtract” command resulting in 10,246 NK cell enhancers (Source data 1-mouse NK cell enhancers).

SNPsplit chromosome of origin reads analysis

Request a detailed protocol

Delineation of allele-informative reads was performed similarly as in reference (Xu et al., 2017). SNPs between the C57BL/6 (B6) and BALB/cJ (BALB) mouse strains were sourced from the Wellcome Sanger Institute Mouse Genomes Project dbSNP (v142). In order to perform unbiased alignment of reads originating from both the B6 and BALB genomes, SNPs marked by the database were replaced by ‘N’ in the mm10 reference genome that we use for alignment using SNPsplit (Babraham Institute) (Krueger and Andrews, 2016). ATAC-seq datasets generated in (B6 x BALB) F1 hybrid NK cells were then aligned to the N-masked genome using bowtie2 and further processed and normalized as above. Reads that overlapped the annotated SNPs were marked as allelically informative reads after alignment and quality control using SNPsplit. Allele-informative reads were then processed and normalized as described above. ~4% of ATAC-seq reads across the dataset were allele-informative.

ChromHMM construction of three-state model

Request a detailed protocol

CUT&RUN data for four histone modifications (H3K4me3, H3K9me3, H3K27ac, H2AK119ub1) generated in cells expressing neither allele (DN) or both alleles (DP) of Ly49G2 were separately used to construct chromatin states using ChromHMM (v1.22) (Ernst and Kellis, 2012). The genome was segmented into three distinct states: state 1 (active chromatin; enriched in H3K4me3 and H3K27ac), state 2 (inactive chromatin; lacking enrichment of all four marks), and state 3 (repressed chromatin enriched in H3K9me3 and H2AK119ub1). The resultant .bed file outputs were visualized with IGV. Chromatin states for both Ly49G2 DN and DP cells are provided in Source data 4.

Statistical analysis

Request a detailed protocol

In vivo germline-edited mouse data were compared with one-way ANOVAs with Tukey’s multiple comparisons (when three or more genotypes were compared) or student’s t-tests (when only two groups were compared). Ex vivo edited NK cell experiments were analyzed by ratio paired t-tests comparing experimental and control samples within a single experiment. In all cases, *p < 0.05; **p < 0.01; ***p < 0.001; **** < 0.0001.

Appendix 1

Appendix 1—key resources table
Reagent type (species) or resourceDesignationSource or referenceIdentifiersAdditional information
Genetic reagent (Mus musculus)C57BL/6J backgroundJackson LaboratoryStrain #:000664 RRID:IMSR_JAX:000664
Genetic reagent (Mus musculus)BALB/cJ backgroundJackson LaboratoryStrain #:000651 RRID:IMSR_JAX:000651
Genetic reagent (Mus musculus)AKR/JJackson LaboratoryStrain #:000648 RRID:IMSR_JAX:000648
Genetic reagent (Mus musculus)CBA/JJackson LaboratoryStrain #:000656 RRID:IMSR_JAX:000656
Genetic reagent (Mus musculus)CB6F1/JJackson LaboratoryStrain #:100007 RRID:IMSR_JAX:100007Purchased from the Jackson Laboratory or generated in-house
Genetic reagent (Mus musculus)B6CBAF1/JJackson LaboratoryStrain #:100011 RRID:IMSR_JAX:100011
Genetic reagent (Mus musculus)B6AKRF1/JJackson LaboratoryGenerated in-house by crossing C57BL/6J with AKR/J
Genetic reagent (Mus musculus)B6.Cg-Klrk1tm1Dhr/JJackson LaboratoryStrain #:022733 RRID:IMSR_JAX:022733
Genetic reagent (Mus musculus)B6-Klra1Hss1ΔThis paperGenerated by Cas9 RNP microinjection
Genetic reagent (Mus musculus)B6-Klrc15’EΔThis paperGenerated by Cas9 RNP electroporation (CRISPR-EZ)
Genetic reagent (Mus musculus)B6-Klrk15’EΔThis paperGenerated by Cas9 RNP electroporation (CRISPR-EZ)
Genetic reagent (Mus musculus)B6-Klra7Hss5ΔThis paperGenerated by Cas9 RNP electroporation (CRISPR-EZ)
Antibody(Armenian hamster monoclonal) anti-CD3e (clone 145–2C11) in PE-Cy5BiolegendRRID: AB_312667 (BioLegend Cat. No. 100302)FACS (1:400)
Antibody(rat monoconal) anti-CD4 (clone GK1.5) in BUV737BD BiosciencesBD Biosciences Cat# 612844, RRID:AB_2870166FACS (1:200)
Antibody(rat monoclonal) anti-CD19 (clone 6D5) in PE-Cy5BiolegendRRID: AB_313644 (BioLegend Cat. No. 115509)FACS (1:400)
Antibody(rat monoclonal) anti-F4/80 (clone BM8) in PE-Cy5BiolegendRRID: AB_893482 (BioLegend Cat. No. 123112)FACS (1:400)
Antibody(rat monoclonal) anti-Ter119 (clone TER-119) in PE-Cy5BiolegendRRID: AB_313711 (BioLegend Cat. No. 116210)FACS (1:400)
Antibody(rat monoclonal) anti-NKp46 (clone 29A1.4) in BV421BiolegendRRID: AB_10915472 (BioLegend Cat. No. 137611)FACS (1:100)
Antibody(mouse monoclonal) anti-NKG2AB6 (clone 16a11) in PEBiolegendRRID: AB_10959654 (BioLegend Cat. No. 142803)FACS (1:50)
Antibody(mouse monoclonal) anti-Ly49AB6 (clone A1) in PEBiolegendRRID: AB_2134787 (BioLegend Cat. No. 138703)FACS (1:50)
Antibody(rat monoclonal) anti-NKG2D (clone CX5) in PE-Dazzle 504BiolegendRRID: AB_2728147 (BioLegend Cat. No. 130213)FACS (1:100)
Antibody(rat monoclonal) anti-CD8β (clone YTS156.7.7) in PE-Cy7BiolegendRRID: AB_2562777 (BioLegend Cat. No. 126616)FACS (1:200)
Antibody(mouse monoclonal) anti-CD45.1 (clone A20) in APCBiolegendRRID: AB_313503 (BioLegend Cat. No. 110714)FACS (1:200)
Antibody(mouse monoclonal) anti-CD45.2 (clone 104) in FITCBiolegendRRID: AB_313443 (BioLegend Cat. No. 109806) Cat. No. 109806FACS (1:200)
Antibody(rat monoclonal) anti-CD90.2 (clone 53–2.1) in FITCBiolegendRRID: AB_10641145 (BioLegend Cat. No. 140308)FACS (1:200)
Antibody(rat monoclonal) anti-CD90.2 (clone 53–2.1) in PEBiolegendRRID: AB_10641145 (BioLegend Cat. No. 140308)FACS (1:200)
Antibody(goat polyclonal) anti-mouse IgG (Poly4053) in PEBiolegendRRID: AB_315010 (BioLegend Cat. No. 405307)FACS (1:200)
Antibody(mouse monoclonal) anti-NKG2A (clone 20d5) in PerCP-eFlour 710eBioscience/ThermoFisherRRID: AB_10853352 (Catalog # 46-5896-82)FACS (1:100)
Antibody(mouse monoclonal) anti-Ly49I (clone YLI-90) in FITCeBioscience/ThermoFIsherRRID: AB_2534426 Catalog # A15413FACS (1:100)
Antibody(mouse monoclonal) anti-Ly49G2 (clone 4D11) in PerCP-eFlour 710eBioscience/ThermoFIsherRRID: AB_1834437 Catalog # 46-5781-82FACS (1:100)
Antibody(mouse monoclonal) anti-CD90.1 (clone HIS51) in FITCeBioscience/ThermoFIsherRRID: AB_465151 Catalog # 11-0900-81FACS (1:100)
Antibody(donkey polyclonal) anti-rat IgG F(ab’)2 in APCeBioscience/ThermoFIsherRRID:AB_469453 polyclonal, lot 17-4822-82 (discontinued)FACS (1:200)
Antibody(mouse monoclonal) anti-CD8.1 (clone 116–13.1)BioXCellRRID: AB_10949065 Catalog # BE0118FACS (1:250)
Antibody(rat monoclonal) anti-CD8.2 (clone 2.43)BioXCellRRID: AB_1125541 Catalog # BE0061FACS (1:50)
Antibody(mouse monoclonal) anti-Ly49A (clone JR9) biotin conjugatedPurified in-house. Ref: Roland and Cazenave Int. Immunol. 1992 PMID: 1535510FACS (1:100)
Antibody(mouse monoclonal) anti-Ly49G2B6 (clone 3/25) unconjugatedUsed as ascites Ref: Tanamachi et al. J. Exp. Med. 2001 PMID: 11157051FACS (1:100)
Antibody(rat monoclonal) anti-NKG2D (clone MI-6) conjugated to biotin in-houseeBioscience/ThermoFIsherRRID: AB_494129 Catalog # 16-5880-86FACS (1:100)
Antibody(rabbit polyclonal) anti-H3K4me1 (ab8895)AbcamAbcam Cat# ab8895, RRID:AB_306847CUT&RUN (1:50)
Antibody(rabbit polyclonal) anti-H3K4me2 (ab7766)AbcamAbcam Cat# ab7766, RRID:AB_2560996CUT&RUN (1:50)
Antibody(rabbit polyclonal) anti-H3K4me3 (ab8580)AbcamAbcam Cat# ab8580, RRID:AB_306649CUT&RUN (1:50)
Antibody(rabbit polyclonal) anti-H3K27ac (ab4729)AbcamAbcam Cat# ab4729, RRID:AB_2118291CUT&RUN (1:50)
Antibody(rabbit polyclonal) anti-H3K9me3 (ab8898)AbcamAbcam Cat# ab8898, RRID:AB_306848CUT&RUN (1:50)
Antibody(rabbit monoclonal) anti-H3K27me3 (clone C36B11)Cell SignalingCell Signaling Technology Cat# 4395, RRID:AB_11220433CUT&RUN (1:50)
Antibody(rabbit monoclonal) anti-H2AUb1 (clone D27C4)Cell SignalingCell Signaling Technology Cat# 8240, RRID:AB_10891618CUT&RUN (1:50)
Antibody(mouse monoclonal) IgG2a k (clone MOPC-173)BiolegendCat # 400202CUT&RUN (1:50)
Commercial assay or kitHiScribe T7 Quick High Yield RNA Synthesis KitNew England BiolabsCat # E2050S
Commercial assay or kitiScript cDNA Synthesis KitBio-radCat # 1708890
Commercial assay or kitP3 Primary Cell 4D-Nucleofector X Kit SLonzaCat #: V4XP-3032
Commercial assay or kitMojosort Mouse NK Cell Isolation KitBiolegendCat # 480049
Commercial assay or kitNextera DNA Library Prep KitNexteraCat # FC-121–1030
Commercial assay or kitNEBNext Ultra II DNA Library Prep Kit for IlluminaNew England BiolabsCat # E7645S
Commercial assay or kitNEBNext Multiplex Oligos for Illumina (Index Primers Set 1)New England BiolabsCat # E7335S
Commercial assay or kitEpitect Fast DNA Bisulfite KitQiagenCat # 59,824
Commercial assay or kitInvitrogen TOPO TA Cloning Kit for Subcloning, without competent cellsThermoFisherCat # Invitrogen 450641
Commercial assay or kitLIVE/DEAD Fixable Near-IR Dead Cell Stain KitThermoFisherCat # L34975Used at (1:1000) in PBS for flow cytometry
Peptide, recombinant proteinRecombinant human IL-2 (teceleukin)National Cancer Institute (BRB Preclinical Biologics Repository)Used at 1,000U/mL for NK cell culture
Peptide, recombinant proteinCas9-NLS (40uM)UC QB3 MacroLab15.6mg used per nucleofection reaction
Peptide, recombinant proteinpA-MNase (Batch #6, 143mg/ml)Kindly provided by the Henikoff labUsed at 0.7mg/mL for CUT&RUN
Peptide, recombinant proteinEpiTaq HS (for bisulfite-treated DNA)TakaraCat # R110B
OtherDAPIBiolegendCat # 422801Used at (1:2000) for flow cytometry
OtherDynabeads mouse T-activator CD3/CD28ThermoFisherCat # 11,456D
OtherBioMagPlus Concanavalin ABangs LaboratoriesCat # BP531
OtherAmpureXP beads, 5mLBeckman CoulterCat # A63880
Software, algorithmFlowJo Version 10FlowJohttps://www.flowjo.com/ RRID:SCR_008520
Software, algorithmPython 3.7.4Pythonhttps://www.python.org/downloads/release/python-374/
Software, algorithmBowtie 2.1.1DOI: 10.1038/nmeth.1923
Software, algorithmPicardBroad Institutehttps://broadinstitute.github.io/picard/
Software, algorithmSAMtools 1.8DOI:10.1093/bioinformatics/btp352
Software, algorithmBEDToolsDOI: 10.1093/bioinformatics/btq033
Software, algorithmMACS2DOI: 10.1186/gb-2008-9-9-r137
Software, algorithmbedGraphToBigWigDOI: 10.1093/bioinformatics/btq351
Software, algorithmIGVDOI: 10.1093/bib/bbs017
Software, algorithmSeqPlotsDOI: 0.12688/wellcomeopenres.10004.1
Software, algorithmdeepToolsDOI: 10.1093/nar/gku365

Data availability

Sequencing data have been deposited in GEO under the accession code GSE181197.

The following data sets were generated
    1. Muljo SA
    2. Raulet DH
    (2021) NCBI Gene Expression Omnibus
    ID GSE181197. Binary outcomes of enhancer activity underlie stable random monoallelic expression.
The following previously published data sets were used
    1. Levanon D
    2. Negreanu V
    3. Lotem J
    4. Bone KR
    5. Brenner O
    6. Leshkowitz D
    7. Groner Y
    (2014) NCBI Gene Expression Omnibus
    ID GSE52625. Runx3 Regulates Interleukin-15-Dependent Natural Killer Cell Activation.
    1. Shih HY
    2. Sciume G
    3. Mikami Y
    4. Guo L
    5. Sun HW
    6. Brooks SR
    7. Urban JF
    8. Davis FP
    9. Kanno Y
    10. O'Shea JJ
    (2016) NCBI Gene Expression Omnibus
    ID GSE77695. Developmental Acquisition of Regulomes Underlies Innate Lymphoid Cell Functionality.
    1. Shih HY
    2. Sciume G
    3. Mikami Y
    4. Nagashima H
    5. Sun HW
    6. Brooks SR
    7. Jankovic D
    8. Yao C
    9. Villarino A
    10. Davis FP
    11. Kanno Y
    12. O'Shea JJ
    (2020) NCBI Gene Expression Omnibus
    ID GSE145299. Rapid remodeling of poised chromatin landscapes and transcription factor repurposing facilitate gene induction in natural killer cells.
    1. Zhou J
    2. Liu D
    (2021) NCBI Gene Expression Omnibus
    ID GSE167237. DNA methylation regulates haematopoietic development.

References

    1. Kubota A
    2. Kubota S
    3. Lohwasser S
    4. Mager DL
    5. Takei F
    (1999)
    Diversity of NK cell receptor repertoire in adult and neonatal mice
    Journal of Immunology (Baltimore, Md 163:212–216.
    1. Yokoyama WM
    2. Kehn PJ
    3. Cohen DI
    4. Shevach EM
    (1990)
    Chromosomal location of the Ly-49 (A1, YE1/48) multigene family: Genetic association with the NK 1.1 antigen
    Journal of Immunology (Baltimore, Md 145:2353–2358.

Decision letter

  1. Deborah Bourc'his
    Reviewing Editor; Institut Curie, France
  2. Tadatsugu Taniguchi
    Senior Editor; Institute of Industrial Science, The University of Tokyo, Japan
  3. Anne-Valerie Gendrel
    Reviewer; Institute Curie, France
  4. Alvaro Rada-Iglesias
    Reviewer; Instituto de Biomedicina y Biotecnologia de Cantabria, Spain

Our editorial process produces two outputs: i) public reviews designed to be posted alongside the preprint for the benefit of readers; ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Decision letter after peer review:

Thank you for submitting your article "Binary outcomes of enhancer activity underlie stable random monoallelic expression" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Tadatsugu Taniguchi as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Anne-Valerie Gendrel (Reviewer #1); Alvaro Rada-Iglesias (Reviewer #2).

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Overall, the reviewers found the manuscript to be relevant, the data clearly presented and the main conclusions well supported by the data. They appreciated the usage of an impressive battery of mouse models, and recognized novel insights brought by this study, into how RME can be regulated and how enhancers can work in a binary manner controlling the probability rather than the level of gene expression. Two slight critics may be (1) the limited mechanistic insights the manuscript provides onto how enhances contribute to RME, and (2) whether the enhancer-RME process described here applies to other types of genes and other cell types, i.e. for genes that are rather CpG-rich than CpG-poor, that are not expressed in the hematopoietic lineage and that do not encode surface markers.

Please refer to the two points further expanded below to prepare your revisions, and please consult the detailed points raised by the reviewers, for further rephrasing of the manuscript and discussion in a point-by-point rebuttal.

1 – Provide additional analyses to investigate the nature of the promoters studied here, before making the claim that RME may be more widespread than originally thought. Using published data, it should be feasible to study these peculiar promoters for their CpG content, and motifs for transcription factor binding sites. In comparison, are genes with CpG-rich promoters also subject to RME and do enhancers influence RME at CpG-rich genes? Finally, although the authors suggest in the discussion that DNA methylation is unlikely to play a major role in RME, it is definitely worth investigating, as an attempt to explain the lack of enhancer-gene communication and also to link with a potential influence on RME extent. Promoters of active and inactive genes should be compared for their DNA methylation states, either using public datasets (if they exist) or using targeted PCR-based DNA methylation analysis.

2 – Moderate the tone regarding the claims that HSSs identified by ATAC-seq are definitely enhancers and not promoters (solely based on epigenomic profiling), provided that this topics is quite controversial in the field.

Reviewer #1 (Recommendations for the authors):

Here are my specific comments and questions for the authors to improve the manuscript.

1. Throughout the manuscript, mature NK, immature NK, primary NK, IL-2 expanded NK cells are mentioned in different instances. It is unclear whether this matters or not? Is this important for the expression of the different genes analysed? This should be clarified, especially for non-specialists in these lineages.

What is the stimulation of NK cells with IL-2 important for; is this for proliferation, for activity? Why are cells sometimes expanded with IL-2 and sometimes not?

2. How were the % of expression frequencies on Figure 1—figure supplement 1A calculated?

3. There is no ATAC-seq peak visible for the promoter of the Ly49a gene (Figure 1A, Figure 3A, Figure 4A), unlike for the other genes of the cluster. Could the authors please comment on this? Doesn´t this argue against the claim that HSS1 is an enhancer and not a promoter?

4. On the FACS plot on Figure 2A, it is difficult to distinguish the population of cells expressing the B6 allele only from both alleles. It would help to draw a circle around each cell populations that were selected for the sorting and further analysed for ATAC-seq. Moreover, was there a control of populations after sorting?

5. On Figure 2E, for the control mouse IgG2ak (cIgG) antibody in IL-2 expanded NK cells, there is no ATAC-seq peak at the enhancer (nor CUT&RUN signal). Shouldn't an ATAC-seq peak at the enhancer expected in this population too? Isn't it similar to cells expressing neither alleles?

6. On Figure 3—figure supplement 2C and F. Isn't there a slight increase in the % of cells expressing the BALB only allele in heterozygous F1 mice for both Ly49a and Nkg2a?

7. On lines 344-345, the authors say: "the corresponding region of the highly related Ly49a gene, which is expressed by only ~17% of NK cells, is much less accessible and presumably less active (Figure 4A)". How is the level of accessibility seen/observed? This is not clear to me.

8. Figure 4E – similar to comment #4, isn't the % of cells expressing the BALB allele only increased? …which would be indicative of a slight dosage compensation mechanism?

9. The possible role of the HSS5 enhancer in increasing the expression probability/frequency of the Ly49g gene in NK cells is appealing. Given the similarity between the Ly49a and Ly49g loci regarding the HSS1 enhancer, would it be possible to test this by introducing the HSS5 element upstream of the Ly49a gene to see whether this would increase the probability/frequency of Ly49a expression? (maybe using a transgene construct, not necessarily a knockin).

10. Line 442; there is no introduction for the knockout model of Nkg2d. What is the KO allele exactly? What is deleted: promoter, exons, whole locus?

11. Line 464; The authors mention the 3'enhancer for the Nkg2d gene. Could they speculate on its function?

12. Line 573-575; the following sentence is not clear to me: "Importantly, NKG2D- cells in both Nkg2d5'ED/5'ED and Nkg2d+/- mice were as likely to express NKG2A, Ly49G2 or Ly49I as NKG2D+ cells, suggesting stochastic expression of the WT Nkg2d allele (Figure 7—figure supplement 1, A-C)." Could the authors please clarify?

13. The following two important points remain unanswered in this study and could be discussed further in the discussion part of the manuscript:

a – What define the strength of an enhancer for the Ly49 and Nkg2 genes. What is the difference between a weak and a strong enhancer? Related to the binding of specific transcription/chromatin factors? Level of enrichment of specific histone marks? 3D conformation?

b – Could this model apply to chromatin or nuclear factors, expressed in other cell lineages than the hematopoietic lineage?

Reviewer #2 (Recommendations for the authors):

Overall, I found the manuscript highly relevant and the main conclusions are well supported by the presented data. However, there are some aspects that could be modified and/or expanded in order to improve the overall relevance of the work:

My main concern is that, as presented, the work strongly suggest that enhancers play an important and instructive work in RME. However, there are limited insights into the molecular mechanisms behind such an instructive role. Given that the promoter but not the enhancer chromatin state correlates with RME, I think the authors should explore more how enhancer-gene communication might be altered in expressing and non-expressing alleles. Since the investigated enhancers are quite proximal to the target gene promoters, investigating topological aspects of enhancer-gene communication can be technically quite challenging. Moreover, it is highly intriguing that in alleles subject to RME the promoters stably lose their responsiveness to proximal active enhancers. In this regard, I have a few suggestions:

– Although H3K27ac is strongly associated with active enhancers, it can be uncoupled from expression of the target genes (Pachano et al., 2021), while the production of eRNAs seems tope strongly coupled with mRNA expression (Kim et al., 2010; Hirabayashi et al., 2019). Therefore, I think it could be quite relevant to investigate eRNA levels at the investigated enhancers in the context of expressing and non-expressing alleles.

– As the authors clearly show, the promoters of the NK receptor genes do not acquire histone modifications associated with gene repression. This is not surprising given that, as the authors mention in the discussion, these promoters are CpG-poor. Inactive CpG-poor promoters do not typically show enrichment for these histone marks but display high DNA methylation levels instead, in contrast with CpG-rich promoters, which are typically hypomethylated regardless of whether they are active or not. Therefore, although the authors mention in the Discussion that DNA methylation is unlikely to play a major role in RME, I still think that it would be quite insightful to compare the DNA methylation levels of active and inactive promoters for example by traditional locus-specific bisulfite sequencing. Moreover, the authors could investigate whether their findings still apply to genes with CpG-rich promoters (are CpG-rich genes subject to RME? Do enhancers influence RME at CpG-rich genes?) and/ or, alternatively, introduce a CpG island into one of the NK receptor gene promoters.

– If the authors find clear differences in the DNA methylation status of active and inactive promoters, one possibility is that epigenetic heterogeneity at gene promoters might occur before enhancer activation and thus influence, together with enhancer strength, the extent of RME. The authors could investigate DNA methylation at promoter regions at earlier cellular states, which could represent a critical temporal window at which DNA demethylating agents could have an impact on RME.

Reviewer #3 (Recommendations for the authors):

1 – In the first section of the results, the authors test the idea that Ly49Hss and NKg25'E regulatory elements are enhancers rather than promoters in NK cells as previously suggested. A full demonstration of this is not provided since the data presented is essentially epigenetic features and ATAC-seq. The provided evidence does not for example demonstrate that Ly49Hss can not behave as promoters in immature cells as proposed before since no RNA data is shown in the article for NKR analyses.

2 – In the same section, based on epigenetic profiles and the relative amount of H3K4me1/3 and other predictive methods (ATAC, K27ac) the authors formally conclude on 'enhancer identity' (lines 186-187) and 'these findings provide definitive support for the conclusion that Ly49Hss1 and Nkg25'E elements represent enhancers' (lines 196-197). It is important to keep in mind that epigenomic features never demonstrate enhancer identity but rather suggest of putative enhancer features. Thus, the phrasing used here should be toned down.

3 – For all ATAC-seq and ChIP-seq data shown, including examples the units should appear on the graphs as normalized equivalently to the read counts (for example in Figure 1B, Figure 1s2C and all the following). A table indicating all published and original data sets used with GSE numbers should be provided, including replicates used and sequencing depth. The presented novel data does not seem to have been submitted to GEO database. This should be done.

4 – If the model proposed by the authors is correct, it implies a level of stochasticity in promoter-enhancer dependent activation. Based on sequence elements (for example presence or absence of CpG islands), do the promoters of the NKR locus tend to be less constitutive? More generally how do the authors envision the critical parameters limiting promoter's action and their strong dependency to enhancer stochastic activation? These points could be discussed in the Discussion section.

https://doi.org/10.7554/eLife.74204.sa1

Author response

1 – Provide additional analyses to investigate the nature of the promoters studied here, before making the claim that RME may be more widespread than originally thought. Using published data, it should be feasible to study these peculiar promoters for their CpG content, and motifs for transcription factor binding sites. In comparison, are genes with CpG-rich promoters also subject to RME and do enhancers influence RME at CpG-rich genes? Finally, although the authors suggest in the discussion that DNA methylation is unlikely to play a major role in RME, it is definitely worth investigating, as an attempt to explain the lack of enhancer-gene communication and also to link with a potential influence on RME extent. Promoters of active and inactive genes should be compared for their DNA methylation states, either using public datasets (if they exist) or using targeted PCR-based DNA methylation analysis.

We address these concerns in detail in the response to the reviewers’ comments below. In summary, RME can be found at genes with both high and low CpG content (Eckersley-Maslin et al., Dev. Cell 2014 PMID: 24576421). As one example from our work, in Figure 7, we show random monoallelic expression of the Thy1 gene, which is expressed by T cells and is widely regarded as a T cell “marker”. As documented in the UCSC genome browser screenshot depicted in reviewer data, Figure 2 (embedded in the response to reviewer 2) the Thy1 gene includes a CpG island at its promoter. Furthermore, the Bcl11b gene, which encodes a transcription factor that commits thymocytes to the T cell fate, contains a promoter CpG-island and was recently shown to be an RME gene and whose expression probability is further reduced by enhancer-deletion (Ng et al., eLife 2018 PMID: 30457103). Therefore, we do not believe there is reason to believe that our results are unique to a class of genes with respect to CpG density; rather, the RME phenomenon (as well as the enhancer deletion-associated variegation phenomenon) appears to apply across cell types and gene classes.

Concerning DNA methylation, we have extensively revised the manuscript to address this issue and will discuss this below in response to reviewer 2.

2 – Moderate the tone regarding the claims that HSSs identified by ATAC-seq are definitely enhancers and not promoters (solely based on epigenomic profiling), provided that this topics is quite controversial in the field.

We softened the language, but note that our data provide evidence beyond epigenomic profiling that these elements are enhancers in mature NK cells: specifically, we showed that deletion of Hss1 in mature NK cells extinguishes expression, which is known to initiate at promoters well downstream of the enhancers in mature NK cells.

Reviewer #1 (Recommendations for the authors):

Here are my specific comments and questions for the authors to improve the manuscript.

1. Throughout the manuscript, mature NK, immature NK, primary NK, IL-2 expanded NK cells are mentioned in different instances. It is unclear whether this matters or not? Is this important for the expression of the different genes analysed? This should be clarified, especially for non-specialists in these lineages.

Where relevant, we have attempted to explain the characteristics of each type of NK cell. Briefly, immature NK cells are not yet terminally differentiated and importantly do not yet express Ly49 genes, expression of which is associated with mature, or fully differentiated NK cells. This is now indicated on page 10.

What is the stimulation of NK cells with IL-2 important for; is this for proliferation, for activity? Why are cells sometimes expanded with IL-2 and sometimes not?

Generally, we use IL-2 to expand NK cells when it is necessary obtain sufficient cells for analysis, which we now indicate in the text on page 15. To elaborate, NK cells are fairly rare among spleen cells (3%) and only 3 % of those co-express both alleles of Ly49G2 in (B6 x BALB/c)F1 hybrid mice, for a frequency of 0.09% of splenocytes. Generally, CUT&RUN experiments utilized IL-2 expanded cells, while ATAC-seq experiments utilized fresh cells. These data appear to cross-validate each other well (Figure 2). Furthermore, Ly49 expression states (both on and off) are stable in NK cells expanded in IL-2. For these reasons, we used IL-2 expanded NK cells when we needed larger numbers of primary cells for chromatin assays.

2. How were the % of expression frequencies on Figure 1—figure supplement 1A calculated?

The percentages are based on previously published historical data generated by both our group and others. We have added citations in the Figure 1—figure supplement 1 legend in the revised manuscript. The numbers depict approximations only.

3. There is no ATAC-seq peak visible for the promoter of the Ly49a gene (Figure 1A, Figure 3A, Figure 4A), unlike for the other genes of the cluster. Could the authors please comment on this? Doesn´t this argue against the claim that HSS1 is an enhancer and not a promoter?

The transcription start sites (TSSs) of the Ly49 genes were previously documented by 5’ oligo-capping RACE (Gays et al., PLos One 2011 PMID: 21483805), and for the Ly49 genes these range from ~5-6kb downstream of HSS1, at sites called “Pro2” and “Pro3,” hereafter referred to as the TSSs. The citation is now included in the Figure 2 legend. mRNAs in mature, Ly49+ NK cells do NOT initiate at the HSS1 elements, and this is true across the gene family. Thus, HSS1 is not the promoter for Ly49a or any of the Ly49 genes in mature NK cells, nor has any published report made the claim that it is. While the signals from the documented TSSs for Ly49a are weak, they are present (see Author response image 1). The weakness of these ATAC sites is likely due to (1) Multiple TSSs such that the ATAC-seq accessibility is diffused to several sites, and (2) the fact that less than 20% of cells express the gene. The signal would therefore be ~5 times higher (or more) if the cells were sorted to be Ly49A+.

Author response image 1
ATAC-seq tracks in NK cells at the Ly49a locus show weak signal at the Ly49aHss5 and Ly49aPro elements.

Here we display four ATAC-seq tracks generated from (B6 x BALB/c)F1 NK cells sorted according to Ly49G2 allelic expression status as shown in Figure 2D, but at the Ly49a locus. The locations of Hss1, Hss5, and the major annotated TSSs are indicated.

4. On the FACS plot on Figure 2A, it is difficult to distinguish the population of cells expressing the B6 allele only from both alleles. It would help to draw a circle around each cell populations that were selected for the sorting and further analysed for ATAC-seq. Moreover, was there a control of populations after sorting?

The NKG2AB6-specific mAb (16a11) was first generated in our lab (Vance et al. PNAS 2002 PMID 11782535), and it was subsequently shown by Dixie Mager’s lab that the use of this clone in combination with the NKG2AB6+BALB/c reactive clone 20d5 enables the sorting the populations expressing each separately, which was confirmed at the mRNA level (Rogers et al., J. Immunol. 2006 PMID: 16785537) and is now cited on p. 11. We added a new figure (Figure 2 figure supplement 1) that shows the gating and post sort analysis of the NKG2A-sorted cells and the Ly49G2-sorted cells used to generate ATAC-seq data and the CUT&RUN data presented in Figure 2.

5. On Figure 2E, for the control mouse IgG2ak (cIgG) antibody in IL-2 expanded NK cells, there is no ATAC-seq peak at the enhancer (nor CUT&RUN signal). Shouldn't an ATAC-seq peak at the enhancer expected in this population too? Isn't it similar to cells expressing neither alleles?

The control IgG “cIgG” refers to the antibody used as a control for the CUT&RUN itself, not as a control for sorting; we do not expect to observe a signal with an IgG that should not interact with the chromatin. We regret the omission of a clear explanation. The cells used for the cIgG control are a mixture of NK cells that express either the B6 or BALB/c allele, in approximately equal amounts. In order to clarify this, we changed the coloring of the cIgG label to black to distinguish it from the labeling of the allele expressed. We further explain in the revised legend for Figure 2 the nature of this control experiment and specify which cells were used to generate the cIgG track. The protein A in the pA-Mnase fusion binds mIgG2ak with high affinity, which is the reason for the selection of this control. In general, CnR produces remarkably low noise, and in all cIgG CnR datasets we have produced we have not seen appreciable signal enrichment at ATAC-seq accessible sites.

6. On Figure 3—figure supplement 2C and F. Isn't there a slight increase in the % of cells expressing the BALB only allele in heterozygous F1 mice for both Ly49a and Nkg2a?

In each of the referenced panels, we expect to see an increase in frequency of BALB/c only cells without invoking a compensation mechanism. This is not due to an increased likelihood of expression of the BALB/c allele, but rather because the cells that would have expressed both alleles are no longer able to express the B6 allele as they lack the critical enhancer on the B6 chromosome. In other words, while the proportion of cells expressing only the BALB/c allele increases, the total proportion of cells expressing the BALB/c allele remains unchanged.

If the reviewer is instead referring to the second sub-panel in Figure 3—figure supplement C (where there is a slight increase in the observed BALB/c-only cells over expected) this is a separate point. We expect an increase in the BALB/c only cells for reasons entirely independent of a compensatory effect, as described above. However, this slight increase in observed over expected may imply some slight compensatory mechanism. This increase is slight, and the overall result is that the total number of cells expressing the BALB/c allele is not drastically increased as we might expect if there were a trans-acting transcriptional “counting” mechanism to ensure a particular percentage of cells express NKG2A.

7. On lines 344-345, the authors say: "the corresponding region of the highly related Ly49a gene, which is expressed by only ~17% of NK cells, is much less accessible and presumably less active (Figure 4A)". How is the level of accessibility seen/observed? This is not clear to me.

Figure 4A shows that the region of Ly49a that aligns with Hss5 of Ly49g has only a very weak ATAC signal. This was true for all ATAC-seq datasets generated in NK cells that we have examined. The ATAC-seq data in Reviewer Data Figure 1 also shows a weak Hss5 signal that is barely distinguishable above background, and certainly much less detectable than Ly49g-Hss5 (Figure 2D; Figure 4A).

8. Figure 4E – similar to comment #4, isn't the % of cells expressing the BALB allele only increased? …which would be indicative of a slight dosage compensation mechanism?

See comment 6. Cells that would otherwise express both alleles in a WT mouse that fail to express the B6 allele in the mutant F1 are now identified as BALB only cells, resulting in an increased proportion without any compensation mechanism.

If, again, the reviewer is referring to the slight increase in observed BALB/c-only cells as compared to expected, the effect is even more slight here and is much more consistent with stochastic and independent regulation of alleles rather than a compensatory mechanism that acts in trans.

9. The possible role of the HSS5 enhancer in increasing the expression probability/frequency of the Ly49g gene in NK cells is appealing. Given the similarity between the Ly49a and Ly49g loci regarding the HSS1 enhancer, would it be possible to test this by introducing the HSS5 element upstream of the Ly49a gene to see whether this would increase the probability/frequency of Ly49a expression? (maybe using a transgene construct, not necessarily a knockin).

This is potentially a revealing experiment that we would love to do and plan to do. In order to do it properly, it should be a germline knockin. Transgenes suffer from various technical problems including variable copy number insertions and position effects. We previously successfully used a genomic transgene to establish the proximal regulation of Ly49a (Tanamachi et al., J. Immunol. 2004 PMID: 14707081). The data, however, were variable for different founders, to an extent that we believe would compromise this approach for the question at hand. Knockin approaches for these genes have proven to be extremely challenging in our hands likely due to the homology between the targeting sequence and multiple related Ly49 genes which appear to diminish the frequency of replacement of the desired gene. For these reasons we decided this experiment, while potentially revealing, will be subject of future efforts and, if successful, included in a future publication.

10. Line 442; there is no introduction for the knockout model of Nkg2d. What is the KO allele exactly? What is deleted: promoter, exons, whole locus?

The Nkg2d knockout model was generated and characterized by our group (Guerra et al., Immunity 2008 PMID: 18394936). The gene is comprised by 7 exons, of which exons 2-6 are deleted in the knockout model. We added a statement to this effect and a reference on p. 27 of the revised manuscript.

11. Line 464; The authors mention the 3'enhancer for the Nkg2d gene. Could they speculate on its function?

We added a sentence on p 28 to note the possibility that the 3’ element constitutes the residual enhancer that drives expression after 5’E deletion, and which acts in combination with 5’E to drive a very high probability of Nkg2d expression in NK cells. However, we have no direct evidence at present that the 3’ ATAC-seq peak is actually an enhancer of Nkg2d.

12. Line 573-575; the following sentence is not clear to me: "Importantly, NKG2D- cells in both Nkg2d5'ED/5'ED and Nkg2d+/- mice were as likely to express NKG2A, Ly49G2 or Ly49I as NKG2D+ cells, suggesting stochastic expression of the WT Nkg2d allele (Figure 7—figure supplement 1, A-C)." Could the authors please clarify?

We revised this sentence (p. 33) to lend more clarity. It now states:

“Importantly, in both Nkg2d5′ED/5′ED and Nkg2d+/- mice, NKG2D-negative cells were as likely as NKG2D+ cells to express NKG2A, Ly49G2 or Ly49I (Figure 7—figure supplement 1, A-C). These data suggest that in both knockout and WT mice, NKG2D is randomly expressed with respect to the other RME receptor genes.”

13. The following two important points remain unanswered in this study and could be discussed further in the discussion part of the manuscript:

a – What define the strength of an enhancer for the Ly49 and Nkg2 genes. What is the difference between a weak and a strong enhancer? Related to the binding of specific transcription/chromatin factors? Level of enrichment of specific histone marks? 3D conformation?

We provide a speculative answer in our response to the public reviews and in the revised manuscript on p. 43. On p 43 of the revised manuscript, after proposing that the randomness of RME may reflect the probability of overcoming the energy barrier required for nucleosome remodeling of initially chromatinized promoters, we added the following passage:

“It may be this property that is the measure of “enhancer strength”, i.e., the capacity of the enhancer to overcome the energy threshold presented by the promoter within a limited window during differentiation. The number of enhancer elements in a gene, the availability of relevant transcription factors that bind the enhancers and the affinity with which the enhancer binds the factors are all likely relevant in determining enhancer strength.”

b – Could this model apply to chromatin or nuclear factors, expressed in other cell lineages than the hematopoietic lineage?

We believe the answer is yes. As already mentioned, the gene encoding the transcription factor Bcl11b is expressed in an RME fashion during thymocyte development, and the probability of expression of this gene is further altered by enhancer deletion variegation (Ng et al., eLife 2018 PMID: 30457103). We add a reference to this in this discussion on page 41, indicating that this supports the generality of our model to genes not encoding cell surface receptors. With respect to non-hematopoietic lineages, deletion of enhancers regulating the odorant receptor genes (Khan et al., Cell 2011 PMID: 22078886) as well as trace amine-associated and trace-amine associated receptors (Fei et al., Nat. Commun. 2021 PMID: 34145235) in the olfactory epithelium reduced the expression probability of the target genes, demonstrating that RME rooted in binary enhancer action applies to non-hematopoietic lineage. On page 44 in the revised manuscript, we add statements and citations to these papers in order to clarify that the mechanisms we describe are likely not restricted to hematopoietic cells.

Reviewer #2 (Recommendations for the authors):

Overall, I found the manuscript highly relevant and the main conclusions are well supported by the presented data. However, there are some aspects that could be modified and/or expanded in order to improve the overall relevance of the work:

My main concern is that, as presented, the work strongly suggest that enhancers play an important and instructive work in RME. However, there are limited insights into the molecular mechanisms behind such an instructive role. Given that the promoter but not the enhancer chromatin state correlates with RME, I think the authors should explore more how enhancer-gene communication might be altered in expressing and non-expressing alleles. Since the investigated enhancers are quite proximal to the target gene promoters, investigating topological aspects of enhancer-gene communication can be technically quite challenging. Moreover, it is highly intriguing that in alleles subject to RME the promoters stably lose their responsiveness to proximal active enhancers. In this regard, I have a few suggestions:

– Although H3K27ac is strongly associated with active enhancers, it can be uncoupled from expression of the target genes (Pachano et al., 2021), while the production of eRNAs seems tope strongly coupled with mRNA expression (Kim et al., 2010; Hirabayashi et al., 2019). Therefore, I think it could be quite relevant to investigate eRNA levels at the investigated enhancers in the context of expressing and non-expressing alleles.

We appreciate the close relationship between eRNA and target mRNA production. While we have not directly assayed eRNAs ourselves, they have been analyzed extensively by RT-PCR in both primary NK cells and NK cell related cell lines (Gays et al., J. Immunol. 2015 PMID: 25926675) (cited on p. 9 of revised manuscript). This study showed that bidirectional enhancer-derived transcripts are present in cells that do not express the corresponding receptor gene. Therefore, even when enhancer-derived transcripts are used to measure enhancer activity, Hss1 activity can be decoupled from target gene expression status, corroborating our results. Taken together, these results suggest that the relationship between the enhancer activity acting on a locus is fundamentally probabilistic rather than deterministic (i.e., not 1:1). This decoupling likely takes place infrequently, when overall enhancer activity is below the threshold to activate the target gene in the large majority of cells. In many and perhaps most cases, enhancer activity as measured by eRNA production will correlate very strongly with mRNA production, especially when measured in bulk populations rather than single cells. We added a sentence on p. 14 addressing the relationship between eRNA production and target Ly49 gene expression status.

– As the authors clearly show, the promoters of the NK receptor genes do not acquire histone modifications associated with gene repression. This is not surprising given that, as the authors mention in the discussion, these promoters are CpG-poor. Inactive CpG-poor promoters do not typically show enrichment for these histone marks but display high DNA methylation levels instead, in contrast with CpG-rich promoters, which are typically hypomethylated regardless of whether they are active or not. Therefore, although the authors mention in the Discussion that DNA methylation is unlikely to play a major role in RME, I still think that it would be quite insightful to compare the DNA methylation levels of active and inactive promoters for example by traditional locus-specific bisulfite sequencing. Moreover, the authors could investigate whether their findings still apply to genes with CpG-rich promoters (are CpG-rich genes subject to RME? Do enhancers influence RME at CpG-rich genes?) and/ or, alternatively, introduce a CpG island into one of the NK receptor gene promoters.

We regret that we did not dedicate more space in the first version of the paper to the discussion of DNA methylation. We have added data and much discussion on the topic, including a section in the results beginning on p 28 and supplementary figure (Figure 5—figure supplement 2) showing differential methylation of the Nkg2d promoter on active and silent alleles in Nkg2d5’E deletion mice. We further revised the discussion to address the role of DNA methylation in RME beginning on p44. Briefly, many RME genes do not display differential promoter methylation (da Rocha, Gendrel Essays Biochem. PMID: 31782494). Recent studies have further indicated that only a small minority of RME genes are subject to derepression of silent alleles through pharmacological inhibition of DNA methylation (Gupta et al., G3 2021 https://doi.org/10.1093/g3journal/jkab428 ; Marion-Poll et al., Nat. Commun. 2021 PMID: 34504093). Furthermore, inhibition of DNA methylation did not derepress silent Igh alleles in a B cell hybridoma model of enhancer-deletion variegation (Ronai et al., J. Immunol. 2002 PMID: 12471125). The role of DNA methylation in RME was addressed in a recent review (da Rocha, Gendrel Essays Biochem. PMID: 31782494), and it appears to be accepted in the RME field that DNA methylation is unlikely to underly all or most cases of RME (Eckersley-Maslin, Spector Trends Genet. 2014 PMID: 24780084). As we mention in the discussion (page 43), there remains the possibility that DNA methylation may play a role in determining the probability of promoter activation for some RME genes, including the ones we studied in the paper.

Regarding the issue of CpG density: Ly49 promoters are indeed CpG-poor, but it is the case that CpG-rich genes are also subject to RME (Eckersley-Maslin et al., Dev. Cell 2014 PMID: 24576421). Indeed, the Thy1 gene, which we show is an RME gene despite being considered a marker of all T cells, harbors a CpG island annotated on the UCSC genome browser (Author response image 2). The RME Bcl11b (Ng et al., eLife 2018 PMID: 30457103) also harbors a CpG island (Author response image 3). We add a statement in the discussion beginning (p. 42) to clarify that RME (and enhancer-controlled RME) apply to genes with CpG-dense and CpG poor promoter regions.

Author response image 2
The Thy1 gene harbors a promoter proximal CpG island.

We provide a screenshot of the UCSC genome browser displaying the Thy1 locus with the CpG track. The red box highlights the region of interest. We show in Figure 7 of the manuscript that Thy1, which is widely regarded as a marker of all T cells, is actually an RME gene. CpG islands are defined as those with GC content of >50%, and the ratio of observed/expected CpG dinucleotides is >0.6, as described in the following link: https://genome.ucsc.edu/cgi-bin/hgc?db=mm10&c=chr9&l=44043001&r=44048197&o=44043488&t=44043701&g=cpgIslandExt&i=CpG%3A+16.

Author response image 3
The Thy1 gene harbors a promoter proximal CpG island.

Similar to Author response image 2, we provide a UCSC genome browser screenshot of the Bcl11b locus with the CpG track displayed. It was recently shown that the gene, which encodes a key transcription factor commiting thymocytes to the T cell lineage, is expressed in an RME fashion and is further subect to enhancer deletion-associated variegation (Ng et al., eLife 2018 PMID: 30457103). The promoter region of interest is highlight in the red box, and the CpG island is indicated with the red arrow.

– If the authors find clear differences in the DNA methylation status of active and inactive promoters, one possibility is that epigenetic heterogeneity at gene promoters might occur before enhancer activation and thus influence, together with enhancer strength, the extent of RME. The authors could investigate DNA methylation at promoter regions at earlier cellular states, which could represent a critical temporal window at which DNA demethylating agents could have an impact on RME.

With respect to the methylation of CpG sites in the Ly49 and Nkg2a promoters themselves, several studies have been published by Dixie Mager’s group demonstrating a clear correlation between DNA methylation and silent alleles (Rouhi et al., J. Immunol. 2006 PMID: 16493057; Rouhi et al., Mol. Immunol. 2007 PMID: 16750269; Rogers et al., J. Immunol. 2006 PMID: 16785537). Inhibition of DNA methylation on its own appears insufficient to derepress silent alleles, however. In the case of Nkg2a, the promoter region is heavily methylated both in immature cells that do not yet express NKG2A, and in the case of silent alleles in mature cells, suggesting that methylation is not specifically deposited on silent alleles, but rather that active alleles are de-methylated. This is mentioned in the revised discussion (page 43). Furthermore, the Ly49a promoter region (Pro2) was shown to be methylated in non-hematopoietic cells. All of these data support the notion that DNA methylation of the NK receptor gene promoters is associated with the default off state, and that demethylation takes place specifically on activated alleles, likely concomitant with activation. These considerations make it less likely that there exists an immature developmental state in which the two alleles might be differentially methylated, preparatory of subsequent RME. Another important point is that the enhancer deletion variegation (and associated allelic “failure”) was observed in Drosophila (Perry et al., PNAS 2011 PMID: 21825127), where DNA methylation is not thought to be present at detectable levels and thus does not play a major role in gene regulation. This consideration confirms that DNA methylation is not generally required for the binary enhancer-driven RME phenomenon.

As discussed in the revised Discussion section (p 42), we favor a mechanism in which enhancer action must overcome the energy barrier associated with the promoters being initially in an inactive state, and success occurs with a random probability, albeit dependent on enhancer strength. In our thinking, the inactive state that must be overcome could simply be the closed chromatin state. However, we also propose the possibility that promoter DNA methylation could contribute to the energy barrier in some cases.

Reviewer #3 (Recommendations for the authors):

1 – In the first section of the results, the authors test the idea that Ly49Hss and NKg25'E regulatory elements are enhancers rather than promoters in NK cells as previously suggested. A full demonstration of this is not provided since the data presented is essentially epigenetic features and ATAC-seq. The provided evidence does not for example demonstrate that Ly49Hss can not behave as promoters in immature cells as proposed before since no RNA data is shown in the article for NKR analyses.

As requested, we toned down the statements claiming the chromatin modification data provided “definitive” support for enhancer as opposed to promoter identity (p. 10). Nevertheless, we feel that the evidence is strong. Our argument depends not only on the epigenetic marks and ATACseq results, but also on our finding that the elements are necessary for maintaining gene expression in mature NK cells, where it is known that the mRNAs are initiated well downstream of Hss1 (Figure 1—figure supplement 2). The revised text on p. 10-11 emphasizes this point (see our response to reviewer 1, comment 3 for further details). The reviewer points out that we did not examine RNA, but others have, and we do not contest their findings that transcripts initiating at Hss1 and transcribing through the gene have been identified in immature NK cells. The difficulty is that transcripts (eRNAs) are known to initiate from enhancers, so the RNA data does not resolve the issue. Another group has made a good argument that the observed transcripts initiating at Hss1 are indeed eRNAs (Gays et al., J. Immunol. 2015 PMID: 25926675). We note that the existence of eRNAs in other genes has not led to theories that such RNAs are part of a developmental switch that turns those genes on and off. An additional argument in favor of our conclusions is that deletion of the 5’ enhancer of Nkg2d resulted in an expression pattern indistinguishable from that of Ly49 or Nkg2a genes. Since it is unlikely that enhancer deletion would CREATE a bidirectional promoter that drives variegated expression, the simplest hypothesis is that variegation in Ly49 and enhancer deletion variegation alike is the consequence of weak enhancer activity as opposed to the action of a bidirectional upstream promoter.

2 – In the same section, based on epigenetic profiles and the relative amount of H3K4me1/3 and other predictive methods (ATAC, K27ac) the authors formally conclude on 'enhancer identity' (lines 186-187) and 'these findings provide definitive support for the conclusion that Ly49Hss1 and Nkg25'E elements represent enhancers' (lines 196-197). It is important to keep in mind that epigenomic features never demonstrate enhancer identity but rather suggest of putative enhancer features. Thus, the phrasing used here should be toned down.

See our response to the previous comment, where we describe the modifications we made and their basis.

3 – For all ATAC-seq and ChIP-seq data shown, including examples the units should appear on the graphs as normalized equivalently to the read counts (for example in Figure 1B, Figure 1s2C and all the following). A table indicating all published and original data sets used with GSE numbers should be provided, including replicates used and sequencing depth. The presented novel data does not seem to have been submitted to GEO database. This should be done.

We have added the units on the graphs as requested. We have submitted to GEO and provided the accession number: GSE181197. The secure token avcpeismhpehlgd is required to access the series prior to publication. Furthermore, we have provided the metadata for the GEO series, with an added column indicating sequencing depth for each file as the file “220302_GEO submission metadata.”

4 – If the model proposed by the authors is correct, it implies a level of stochasticity in promoter-enhancer dependent activation. Based on sequence elements (for example presence or absence of CpG islands), do the promoters of the NKR locus tend to be less constitutive? More generally how do the authors envision the critical parameters limiting promoter's action and their strong dependency to enhancer stochastic activation? These points could be discussed in the Discussion section.

As noted in our response to Reviewer 2’s comments, we have now expanded the Discussion section to include discussion of our conception of the critical parameters that determine stochasticity. Some of these (enhancer strength) come from our work and others, such as the abundance of relevant transcription factors, come from the work of others, which we cite.

In conclusion, we have made substantial changes to the manuscript to address the comments of all the reviewers. We believe the manuscript is improved as a result.

https://doi.org/10.7554/eLife.74204.sa2

Article and author information

Author details

  1. Djem U Kissiov

    Division of Immunology and Pathogenesis, Department of Molecular and Cell Biology, University of California, Berkeley, Berkeley, United States
    Contribution
    Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Resources, Software, Validation, Visualization, Writing - original draft, Writing - review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-6279-342X
  2. Alexander Ethell

    Division of Immunology and Pathogenesis, Department of Molecular and Cell Biology, University of California, Berkeley, Berkeley, United States
    Contribution
    Data curation, Formal analysis, Methodology, Validation, Visualization
    Competing interests
    No competing interests declared
  3. Sean Chen

    Division of Immunology and Pathogenesis, Department of Molecular and Cell Biology, University of California, Berkeley, Berkeley, United States
    Contribution
    Investigation, Methodology, Resources
    Competing interests
    No competing interests declared
  4. Natalie K Wolf

    Division of Immunology and Pathogenesis, Department of Molecular and Cell Biology, University of California, Berkeley, Berkeley, United States
    Contribution
    Data curation, Investigation, Methodology, Writing - review and editing
    Competing interests
    No competing interests declared
  5. Chenyu Zhang

    Division of Immunology and Pathogenesis, Department of Molecular and Cell Biology, University of California, Berkeley, Berkeley, United States
    Contribution
    Investigation, Methodology, Writing - review and editing
    Competing interests
    No competing interests declared
  6. Susanna M Dang

    Division of Immunology and Pathogenesis, Department of Molecular and Cell Biology, University of California, Berkeley, Berkeley, United States
    Contribution
    Investigation, Methodology
    Competing interests
    No competing interests declared
  7. Yeara Jo

    Division of Immunology and Pathogenesis, Department of Molecular and Cell Biology, University of California, Berkeley, Berkeley, United States
    Contribution
    Investigation, Methodology
    Competing interests
    No competing interests declared
  8. Katrine N Madsen

    Division of Immunology and Pathogenesis, Department of Molecular and Cell Biology, University of California, Berkeley, Berkeley, United States
    Contribution
    Investigation, Methodology, Writing - review and editing
    Competing interests
    No competing interests declared
  9. Ishan Paranjpe

    Division of Immunology and Pathogenesis, Department of Molecular and Cell Biology, University of California, Berkeley, Berkeley, United States
    Contribution
    Investigation, Methodology
    Competing interests
    No competing interests declared
  10. Angus Y Lee

    Cancer Research Laboratory, University of California, Berkeley, Berkeley, United States
    Contribution
    Methodology, Resources, Writing - review and editing
    Competing interests
    No competing interests declared
  11. Bryan Chim

    Laboratory of Immune System Biology, National Institute of Allergy and Infectious Diseases, National Institutes of Health, Bethesda, United States
    Contribution
    Formal analysis, Investigation, Methodology, Writing - review and editing
    Competing interests
    No competing interests declared
  12. Stefan A Muljo

    Laboratory of Immune System Biology, National Institute of Allergy and Infectious Diseases, National Institutes of Health, Bethesda, United States
    Contribution
    Methodology, Project administration, Resources, Writing - review and editing
    Competing interests
    No competing interests declared
  13. David H Raulet

    Division of Immunology and Pathogenesis, Department of Molecular and Cell Biology, University of California, Berkeley, Berkeley, United States
    Contribution
    Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing - original draft, Writing - review and editing
    For correspondence
    raulet@berkeley.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-1257-8649

Funding

National Institute of Allergy and Infectious Diseases (R01-AI113041)

  • David H Raulet

National Institute of Allergy and Infectious Diseases (Intramural Research Program of NIH)

  • Stefan A Muljo

National Science Foundation (DGE 1752814)

  • Natalie K Wolf

Cancer Research Coordinating Committee (Predoctoral fellowship)

  • Djem U Kissiov

National Institute of Allergy and Infectious Diseases (NIAID)

  • Stefan A Muljo

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank Lily Zhang and Erik Seidel for assistance, Hector Nola and Alma Valeros in the Cancer Research Laboratory at UC Berkeley for expert assistance with flow cytometry and cell sorting, and M Wong for assistance with some Illumina sequencing experiments. We thank Drs. Kathleen Pestal, Han-Yu Shih and Giuseppe Sciume for discussion and feedback. We thank Drs. Jasper Rine, Russell Vance, Ellen Robey, Michel DuPage, and David Martin for critical evaluation of the manuscript.

Ethics

This study was performed strictly and in full accordance with the guidelines set forth by the University of California, Berkeley's Animal Care and Use Committee (ACUC). The protocol was approved by ACUC under the protocol number # AUP-2015-10-8058-1. All experiments were performed on mice euthanized according to the protocol.

Senior Editor

  1. Tadatsugu Taniguchi, Institute of Industrial Science, The University of Tokyo, Japan

Reviewing Editor

  1. Deborah Bourc'his, Institut Curie, France

Reviewers

  1. Anne-Valerie Gendrel, Institute Curie, France
  2. Alvaro Rada-Iglesias, Instituto de Biomedicina y Biotecnologia de Cantabria, Spain

Publication history

  1. Preprint posted: August 28, 2021 (view preprint)
  2. Received: September 24, 2021
  3. Accepted: March 21, 2022
  4. Version of Record published: May 26, 2022 (version 1)

Copyright

This is an open-access article, free of all copyright, and may be freely reproduced, distributed, transmitted, modified, built upon, or otherwise used by anyone for any lawful purpose. The work is made available under the Creative Commons CC0 public domain dedication.

Metrics

  • 1,373
    Page views
  • 110
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Djem U Kissiov
  2. Alexander Ethell
  3. Sean Chen
  4. Natalie K Wolf
  5. Chenyu Zhang
  6. Susanna M Dang
  7. Yeara Jo
  8. Katrine N Madsen
  9. Ishan Paranjpe
  10. Angus Y Lee
  11. Bryan Chim
  12. Stefan A Muljo
  13. David H Raulet
(2022)
Binary outcomes of enhancer activity underlie stable random monoallelic expression
eLife 11:e74204.
https://doi.org/10.7554/eLife.74204

Further reading

    1. Chromosomes and Gene Expression
    2. Genetics and Genomics
    Joseph V Geisberg, Zarmik Moqtaderi ... Kevin Struhl
    Research Advance

    Alternative polyadenylation yields many mRNA isoforms whose 3' termini occur disproportionately in clusters within 3' UTRs. Previously, we showed that profiles of poly(A) site usage are regulated by the rate of transcriptional elongation by RNA polymerase (Pol) II (Geisberg et., 2020). Pol II derivatives with slow elongation rates confer an upstream-shifted poly(A) profile, whereas fast Pol II strains confer a downstream-shifted poly(A) profile. Within yeast isoform clusters, these shifts occur steadily from one isoform to the next across nucleotide distances. In contrast, the shift between clusters from the last isoform of one cluster to the first isoform of the next - is much less pronounced, even over large distances. GC content in a region 13-30 nt downstream from isoform clusters correlates with their sensitivity to Pol II elongation rate. In human cells, the upstream shift caused by a slow Pol II mutant also occurs continuously at the nucleotide level within clusters, but not between them. Pol II occupancy increases just downstream of the most speed-sensitive poly(A) sites, suggesting a linkage between reduced elongation rate and cluster formation. These observations suggest that 1) Pol II elongation speed affects the nucleotide-level dwell time allowing polyadenylation to occur, 2) poly(A) site clusters are linked to the local elongation rate and hence do not arise simply by intrinsically imprecise cleavage and polyadenylation of the RNA substrate, 3) DNA sequence elements can affect Pol II elongation and poly(A) profiles, and 4) the cleavage/polyadenylation and Pol II elongation complexes are spatially, and perhaps physically, coupled so that polyadenylation occurs rapidly upon emergence of the nascent RNA from the Pol II elongation complex.

    1. Chromosomes and Gene Expression
    2. Structural Biology and Molecular Biophysics
    Yu Chen, Claudia Cattoglio ... Xavier Darzacq
    Research Article Updated

    Transcription factors (TFs) are classically attributed a modular construction, containing well-structured sequence-specific DNA-binding domains (DBDs) paired with disordered activation domains (ADs) responsible for protein-protein interactions targeting co-factors or the core transcription initiation machinery. However, this simple division of labor model struggles to explain why TFs with identical DNA-binding sequence specificity determined in vitro exhibit distinct binding profiles in vivo. The family of hypoxia-inducible factors (HIFs) offer a stark example: aberrantly expressed in several cancer types, HIF-1α and HIF-2α subunit isoforms recognize the same DNA motif in vitro – the hypoxia response element (HRE) – but only share a subset of their target genes in vivo, while eliciting contrasting effects on cancer development and progression under certain circumstances. To probe the mechanisms mediating isoform-specific gene regulation, we used live-cell single particle tracking (SPT) to investigate HIF nuclear dynamics and how they change upon genetic perturbation or drug treatment. We found that HIF-α subunits and their dimerization partner HIF-1β exhibit distinct diffusion and binding characteristics that are exquisitely sensitive to concentration and subunit stoichiometry. Using domain-swap variants, mutations, and a HIF-2α specific inhibitor, we found that although the DBD and dimerization domains are important, another main determinant of chromatin binding and diffusion behavior is the AD-containing intrinsically disordered region (IDR). Using Cut&Run and RNA-seq as orthogonal genomic approaches, we also confirmed IDR-dependent binding and activation of a specific subset of HIF target genes. These findings reveal a previously unappreciated role of IDRs in regulating the TF search and binding process that contribute to functional target site selectivity on chromatin.