Chromatin proteomics reveals novel combinatorial histone modification signatures that mark distinct subpopulations of macrophage enhancers.

The integrated activity of cis-regulatory elements fine-tunes transcriptional programs of mammalian cells by recruiting cell type-specific as well as ubiquitous transcription factors (TFs). Despite their key role in modulating transcription, enhancers are still poorly characterized at the molecular level, and their limited DNA sequence conservation in evolution and variable distance from target genes make their unbiased identification challenging. The coexistence of high mono-methylation and low tri-methylation levels of lysine 4 of histone H3 is considered a signature of enhancers, but a comprehensive view of histone modifications associated to enhancers is still lacking. By combining chromatin immunoprecipitation (ChIP) with mass spectrometry, we investigated cis-regulatory regions in macrophages to comprehensively identify histone marks specifically associated with enhancers, and to profile their dynamics after transcriptional activation elicited by an inflammatory stimulation. The intersection of the proteomics data with ChIP-seq and RNA-seq analyses revealed the existence of novel subpopulations of enhancers, marked by specific histone modification signatures: specifically, H3K4me1/K36me2 marks transcribed enhancers, while H3K4me1/K36me3 and H3K4me1/K79me2 combinations mark distinct classes of intronic enhancers. Thus, our MS analysis of functionally distinct genomic regions revealed the combinatorial code of histone modifications, highlighting the potential of proteomics in addressing fundamental questions in epigenetics.

INTRODUCTION

In eukaryotes, gene transcription is tightly controlled in time and space by cis-regulatory regions. These include both proximal regions, such as promoters located directly upstream of the transcription start sites (TSS), and distal regulatory elements, such as enhancers, silencers, insulators and locus control regions (LCR), which are platforms for recruiting transcription activators and suppressors that act synergistically to enforce specific transcriptional programs during development, differentiation and response to environmental stimuli.

In particular, enhancers play a crucial role in setting cell- and time-specific gene expression programs. They are capable of promoting transcription of their target genes at a considerable distance (up to several hundred kb), in an orientation-independent manner (1–6). Enhancers are characterized by DNA hypersensitivity to nucleases, nucleosome-depleted regions and clustering of DNA binding-motifs for TFs (7–12). Usually, enhancer activation requires the collaborative interaction of multiple TFs, including ubiquitous, pioneer or lineage-determining (LD) and signal-dependent (SD) TFs (13–17). Lineage-determining TFs can bind specific DNA motifs in compacted chromatin, creating a nucleosome-free region that can be bound by other TFs and co-activators (such as histone modifiers, ATP-dependent chromatin remodelers or bridging factors for long-range cross-talks with the basal transcriptional machinery at promoters) that ultimately control transcription (18,19).

Genome-wide studies indicate that enhancers exhibit a characteristic histone post-translational modification (PTM) ‘signature’, consisting of high levels of mono-methylation, with low levels of tri-methylation, at lysine 4 on histone H3 (H3K4me1Hi/K4me3Low) (20–22). Nevertheless, only a fraction of genomic regions with this signature are inherently active enhancers, which are also marked by H3K27 acetylation and actively transcribing RNA polymerase II (23,24).

Recent evidence suggests that chromatin modifications provide an important layer of enhancer regulation, by either serving as docking sites for signal-dependent TFs or promoting the access to other TFs; this would support the assembly of transcriptional complexes and long-range chromatin interactions, which altogether enforce specific transcriptional responses (14,25). For instance, H2B mono-ubiquitination stabilizes the histone variant H2A.Z at inducible enhancers regulated by the estrogen receptor alpha, which blocks binding of the chromatin remodeler INO80 and, consequently, the transcriptional activation of these enhancers (26).

Genome-wide characterization of chromatin features associated with cis-regulatory elements has been facilitated by ChIP-seq analysis. However, a comprehensive and unbiased picture of histone modification patterns associated with these regions is still lacking (27,28). This information would be invaluable to first identify genomic enhancers and then classify them into distinct subpopulations that regulate specific transcriptional programs, as well as to decipher the molecular mechanisms underlying their function.

Taking advantage of a native chromatin proteomics (N-ChroP) approach (29), we combined mass spectrometry (MS) and native ChIP (N-ChIP) against H3K4me1 and H3K4me3 to determine which histone PTMs are differentially associated with enhancers and promoters in macrophages. Macrophages provide an ideal system, as they have been well-studied at the genomic and epigenomic levels and are highly responsive to inflammatory stimuli, with a rapid and dynamic activation of hundreds of inflammatory genes involved in the innate and adaptive immune response. MS has been previously used to identify and quantify histone PTMs in an unbiased and comprehensive manner (30–35), and to characterize long-distance synergies between different PTMs when combined with histone digestion protocols that produce long polypeptides (36–40). However, MS analysis for epigenetic profiling has been so far largely limited to bulk chromatin, with only a few MS-based studies attempting to dissect histone PTMs in specific chromatin regions (41,42).

Using N-ChroP on resting or lipopolysaccharide (LPS)-stimulated macrophages over a time course, we profiled 42 distinct modified peptides from histones H3 and H4 and identified known and novel marks specifically enriched at enhancers. The unexpected enrichment on histone H3 of lysine 36 di- and tri-methylation (H3K36me2 and H3K36me3) and lysine 79 di-methylation (H3K79me2) in H3K4me1-rich regions suggested the existence of subcategories of enhancers marked by these PTMs. Indeed, by analyzing ChIP-seq data using an unsupervised learning strategy, we identified significant genome-wide co-associations between these marks. Intersecting ChIP-seq data with total and nascent transcript analyses, we also found that H3K4me1/K36me3 and H3K4me1/K79me2 signatures mark intronic enhancers involved in distinct inflammatory gene expression programs, while the dual H3K4me1/K36me2 mark mainly accumulates at transcribed enhancer regions.

By applying N-ChroP to dissect enhancers, we have provided further support to the histone code hypothesis, showing that distinct combinations of marks define functionally different regions of the genome. The novel epigenetic signatures that have emerged in this study will be valuable both for de novo identification of enhancer subpopulations and for dissecting the molecular mechanisms behind the activation of distinct enhancers in response to stimuli.

MATERIALS AND METHODS

Cell culture

RAW 264.7 macrophage-like cells were cultured in high-glucose Dulbecco's modified Eagle's medium (DMEM, Invitrogen) containing 10% low endotoxin FCS, 1% glutamine and 1% Pen/Strep. Cells were stimulated with lipopolysaccharide (LPS) from Escherichia coli serotype EH100 (Alexis) at 100 ng/ml for 1 or 4 h. Transcription was inhibited with 1 μM flavopiridol (Sigma, St. Louis MO, USA) 1 h prior to LPS stimulus, as previously described (43).

Native chromatin immunoprecipitation (N-ChIP)

Approximately 20–30 × 106 of RAW cells were homogenized in lysis buffer (10% sucrose, 0.5 mM EGTA [pH 8.0], 15 mM NaCl, 60 mM KCl, 15 mM HEPES, 0.5% Triton, 0.5 mM PMSF, 1 mM DTT, 5 mM NaF, 5 mM Na3VO4, 5 mM sodium butyrate, 5 μg/ml aprotinin, 5, μg/ml pepstatin A and 5 μg/ml leupeptin), and nuclei were separated from cytoplasm by centrifugation at 3750 rpm (4°C) for 30 min on sucrose cushions. Nuclear pellets were washed twice with phosphate buffered saline (PBS), re-suspended in digestion buffer (0.32 M sucrose, 50 mM Tris–HCl [pH 7.6], 4 mM MgCl2, 1 mM CaCl2 and 0.1 mM PMSF) and digested with micrococcal nuclease (MNase, Roche) at a final concentration of 0.01 U/μl, for 60 min at 37°C. The reaction was stopped by adding 1 mM EDTA and chilling on ice. After centrifugation of re-suspended nuclei at 10 000 rpm (4°C) for 10 min, the soluble fraction of chromatin (S1) comprising smaller fragments (mono- and di-nucleosomes), was collected as the supernatant. Pellets were re-suspended in dialysis buffer (10 mM Tris–HCl [pH 7.6], 1 mM EDTA, 0.5 mM PMSF, 5 mM NaF, 5 mM Na3VO4, 5 mM sodium butyrate and 1 × protease inhibitor cocktail) and dialyzed overnight at 4°C in a dialysis tube (cut off 3.5 kDa). The second soluble fraction of chromatin (S2), comprising large fragments (tri- to hepta-nucleosomes), was likewise obtained after centrifugation at 10 000 rpm (4°C) for 10 min. DNA was extracted using QIAquick columns (QIAGEN) and then run on 1% agarose gel to evaluate the chromatin fractions. The S1 fraction was combined with a small aliquot of the S2 fraction and diluted in ChIP dilution buffer (50 mM Tris–HCl [pH 7.6], 50 mM NaCl and 5 mM EDTA). Chromatin was incubated overnight at 4°C with H3K4me1 (ab8895, 25 μg), H3K4me3 (active motif 39159, 30 μl) or H3K9me3 (ab8898, 25 μg) antibodies. In parallel, 200–250 μl of G protein–coupled magnetic beads (Dynabeads, Invitrogen 100.04D) was blocked in 0.5% BSA in PBS overnight at 4°C. Blocked beads were washed, added to chromatin and incubated for 3 h at 4°C on a rotating wheel. Beads were washed four times (with 50 mM Tris–HCl [pH 7.6] and 10 mM EDTA) using increasing salt concentration for the last three washes (of 75, 125 and 175 mM NaCl). LDS sample buffer (Invitrogen NP0007) supplemented with 50 mM DTT was added to the beads for 5 min at 70°C to elute the immunoprecipitated proteins from the beads. Proteins were resolved on 4–12% Bis-Tris acrylamide SDS-PAGE pre-cast gels (Invitrogen NP0335BOX) on an Invitrogen system and visualized using Colloidal Blue Staining kit (Invitrogen LC6025) to ensure the immunopurification of intact nucleosomes containing all histones at the same stochiometry. To analyze the basal histone PTMs, three biological replicates (n = 3) for each ChIP, using H3K9me3 as a ‘negative’ control for cis-regulatory regions enrichment, were performed while for the dynamic experiment we performed three biological replicates (n = 3) for each time point (1 and 4 h after LPS treatment).

In-gel digestion of core histones

Bands corresponding to histones H3 and H4 were excised and in-gel digested, as previously described (29). Briefly, gel bands were cut in pieces and de-stained with H2O 100% alternated with dehydration steps in 50% acetonitrile (in H2O) and 100% acetonitrile. Gel pieces were then chemically alkylated in-gel by incubation with D6-acetic anhydride (Sigma 175641) 1:9 in 1M NH4HCO3 (60–70 μl final volume) and saturated CH3COONa solution as a catalyzer. After 3 h at 37°C with vigorous shaking in a thermomixer, chemically-modified gel slices were washed with NH4HCO3 and alternating with ACN at increasing % (from 50 to 100). In-gel digestion was performed with 100 ng/μl trypsin (Promega V5113) in 50 mM NH4HCO3 at 37°C overnight, in order to obtain an ‘in-gel’ ArgC-like digestion, which cleaves at the amide bond C-terminal to arginine residues, producing peptides with an optimal length for MS analysis. Finally, digested peptides were collected and extracted using 5% formic acid alternated with ACN 100%.

Digested peptides were desalted and concentrated using a combination of reverse-phase C18/carbon ‘sandwich’ system and ion-exchange (SCX) chromatography, on hand-made nano-columns (StageTips) (44); digested peptides loaded on C18/carbon and SCX StageTips were then eluted with high-organic solvent (80% ACN/0.5% acetic acid) and 5% NH4OH/30% methanol, respectively. Eluted peptides were lyophilized, re-suspended in 1% TFA in ddH2O, pooled and subjected to LC–MS/MS.

Nanoflow liquid chromatography and tandem-mass spectrometry (nLC-MS/MS)

Peptide mixtures were analyzed by online nanoflow liquid chromatography tandem-mass spectrometry (LC–MS/MS) using an EASY-nLC™ 1000 (Thermo Fisher Scientific, Odense, Denmark) connected to a Q-Exactive instrument (Thermo Fisher Scientific) through a nanoelectrospray ion source. The nano LC system was operated in one column set-up with a 25-cm analytical column (75 μm inner diameter, 350 μm outer diameter) packed with C18 resin (ReproSil, Pur C18AQ 1.9 μm, Dr Maisch, Germany) configuration. Solvent A was 0.1% FA in ddH2O, and solvent B was 80% ACN with 0.1% FA. Samples were injected in an aqueous 1% TFA solution at a flow rate of 500 nl/min. Peptides were separated with a gradient of 0–40% solvent B over 100 min followed by a gradient of 40–60% in 5 min and 60–95% over 5 min at a flow rate of 250 nl/min in the EASY-nLC 1000 system. The Q-Exactive instrument was operated in the data-dependent mode (DDA) to automatically switch between full-scan MS and MS/MS acquisition. Survey full scan MS spectra (from m/z 300–1350) were analyzed in the Orbitrap detector with resolution R = 60 000 at m/z 200. The ten most intense peptide ions with charge states ≥2 were sequentially isolated to a target value of 3e6 and fragmented by higher energy collision dissociation (HCD) with a normalized collision energy setting of 28%. The maximum allowed ion accumulation times were 20 ms for full scans and 100 ms for MS/MS, and the target value for MS/MS was set to 1e5 (R = 15 000 at m/z 200). The dynamic exclusion time was set to 20 s for Q-Exactive instrument. Standard mass spectrometric conditions for all experiments were: spray voltage, 2.4 kV; no sheath and auxiliary gas flow.

Analysis of mass spectrometric data and histone PTM profiling

Acquired raw data were analyzed by the integrated MaxQuant software v.1.5.2.8, using the Andromeda search engine (45,46). The MOUSE 1401 database (51 195 entries) was used for peptide identification. Enzyme specificity was set to Arg-C, because the chemical derivatization of lysines with deuterated acetic anhydride prior to trypsin digestion results in an Arg-C–like digestion. In MaxQuant, the estimated false discovery rate (FDR) of all peptide identifications was set to a maximum of 1%. The main search was performed with a mass tolerance of 6 ppm. A maximum of 2 missed cleavages was permitted, and the minimum peptide length was set at six amino acids.

The peptide search focused on lysine methylation and acetylation, with the variable modifications of lysine D3-acetylation (+45.0294 Da), lysine mono-methylation (masses of D3-acetylation [+45.0294] plus mass of mono-methylation [+14.016 Da]), lysine di-methylation (+28.031 Da), lysine tri-methylation (42.046 Da) and lysine acetylation (+42.010 Da) were included. To reduce the rate of false positives, which increases with increasing the number of variable modifications included in the database search (47), raw data were analyzed with multiple and parallel MaxQuant jobs, setting different combinations of variable modifications: (i) D3-acetylation, lysine mono-methylation, lysine di-methylation and lysine acetylation; (ii) D3-acetylayion with lysine mono-methylation and lysine tri-methylation and (iii) D3-acetylayion with tri-methylation and lysine acetylation. MaxQuant search results were exported, and peptides scoring lower than 60 (corresponding to a Mascot score of 15 (46)) and with a localization probability score <0.75 (following previous phospho-PTM studies (48,49)) were removed, as they were identified with low confidence (peptides parameters and output from MQ analysis are given in Supplementary Tables S1 and S2). Filtered data were subjected to manual inspection and validation using QualBrowser version 2.0.7 (ThermoFisher Scientific). Extracted ion chromatograms (XIC) were constructed for each precursor based on the m/z value, using a mass tolerance of 10 ppm and a mass precision of up to four decimals. For each histone-modified peptide isoform, its relative abundance percentage (RA%) was estimated by dividing the area under the curve (AUC) of each modified peptide by the sum of all the AUC values of all observed isoforms of that peptide. Hence, the RA % value for each specific peptide represents the relative proportion of that peptide among all quantified peptide forms sharing the same amino acid sequence (50,51). The relative enrichment of each PTM, both at the basal state and after LPS stimulus, was calculated as a ratio between its RA % in the ChIP and in the input. Visualization of histone PTM enrichment and unsupervised hierarchical clustering were performed with Perseus, setting correlation distances and complete linkage as parameters (52). Differences among histone PTMs associated to a specific chromatin region were assessed by t-test analysis, followed by Bonferroni's post-hoc test. Visualisation of MS/MS spectra for each modification was performed using the Viewer tool (Supplementary Figure S9) in MaxQuant.

ChIP-seq and RNA-seq analyses

For ChIP-seq, ∼10 × 106 of untreated or stimulated (1 or 4 h after LPS) RAW cells were used for each experiment. Cells were crosslinked with 0.75% formaldehyde for 10 min at room temperature, and chromatin was sonicated as described previously (29). Each chromatin input was immunoprecipitated with 10 μg of the following antibodies: H3K4me1 (Abcam 8895), H3K4me3 (Active Motif 39159), H3K36me2 (Abcam 9049), H3K79me2 (Abcam ab3594), H3K36me3 (Abcam 9050), H3K27ac (Abcam 4729), H3K27me3 (Cell Signalling 9733) and Pu.1 (Santa Cruz sc-352). After immunoprecipitation, beads were washed three times in buffer A (20 mM Tris–HCl [pH 7.6], 2 mM EDTA, 0.1% SDS, 1% Triton-100 and 150 mM NaCl), once in buffer B (20 mM Tris–HCl [pH 7.6], 2 mM EDTA, 0.1% SDS, 1% Triton-100 and 300 mM NaCl) and once in TE containing 50 mM NaCl. DNA was eluted in TE containing 2% SDS and de-crosslinked overnight at 65°C. DNA was purified by QIAquick columns (QIAGEN) and quantified with Picogreen. Sequencing libraries were generated as previously described (53,54) and sequenced on an Illumina HiSeq2000. Total RNA was extracted from 2 × 106 cells using RNeasy Kit (QIAGEN) with DNase I treatment. Libraries were then prepared using TruSeq RNA sample preparation Kit (Illumina) after depleting ribosomal RNA and sequenced on an Illumina HiSeq2000.

Computational methods

ChIP-seq analysis

Reads derived from internal ChIP-seq experiments and H3K9me3 ChIP-seq reads taken from (54) were aligned to the mouse genome (assembly mm9) using BowTie (55), allowing up to three mismatches and keeping only reads mapping uniquely (for histone PTMs) or at most at 10 different positions of the genome (for Pu.1). MACS2 was used to identify regions of reads enrichment compared to a control DNA input, either by using default parameters (Pu.1) or by allowing the identification of broad regions of enrichment in the analysis of histone PTMs (56).

To analyze composite combinatorial interactions of the different histone PTMs—in both basal and inflamed states—we employed the automated computational system for learning chromatin states ChromHMM (57). ChromHMM segments the genome into regions of 200 bp, such that it was possible to evaluate the presence/absence for each histone PTM according to the reads enrichment against a control. The resulting data was used to train a multivariate hidden Markov model for unsupervised learning of chromatin states, where each state corresponded to a combination of presence/absence of each modification; training used default parameters, and models up to 35 states were evaluated. In each analysis, the model with the maximum number of non-redundant states was chosen. When analyzing chromatin states across different conditions (untreated or induced with LPS for 1 or 4 h), models were learned using all histone PTMs under analysis, with a distinct cell-type label assigned to all modifications belonging to the same treatment, in order to obtain a chromatin state model applicable to all conditions. For the subsequent analysis on localization of specific chromatin states, regions identified as K4me1-positive (K4me1/K36me2, K4me1/K36me3 or K4me1/K79me2) that were separated on the genome by less than 2 kb were merged, and overlapping regions of H3K4me1-only containing states were removed, in order to obtain less ambiguous calls for the overlap with genomic features.

Genomic annotations of peaks or chromatin states were performed either using HOMER or with the bedIntersect function of BedTools (14,58).

Read distributions around Pu.1 summits were calculated using in-house–developed scripts. Statistical analysis and plots were performed using R. Gene ontology (GO) term enrichment was performed with the latest version of the online tool DAVID (59). Every analysis used standard parameters, with a specific background superset of genes for every group of analyzed genes; for instance, for genes harboring K79me2- or H3K36me3-marked enhancers, a background of the respective population of all H3K79me2 or H3K36me3 positive genes was used, while for the H3K79me2-marked active enhancers gene subpopulation, a background of all genes with an intergenic enhancer marked by H3K79me2 was used.

To identify intronic enhancers transitioning from class II to class I, class II was first defined as only those regions that showed co-localization of both K4me1 and K27me3 in the untreated condition, based on the 27-state ChromHMM model. Overlap between these regions, and emissions marked by both K4me1 and K27ac identified after LPS treatment (1 and 4 h) were determined, and regions were mapped onto RefSeq genes. GO enrichment of the resulting lists was computed with DAVID, using the whole mm9 RefSeq gene annotation as background.

Enrichment analysis of known TF binding motifs was carried out using the corresponding function of the HOMER suite, using the entire length of the input regions to calculate the enrichments.

RNA-seq analysis

Reads resulting from total RNA-seq and nascent RNA-seq (taken from (60), GSE66955) were aligned and quantified over the mouse genome (mm9) RefSeq gene annotation using RSEM (61), with default parameters and the addition of the strand-specific information for the nascent RNA dataset. BedIntersect was used to map chromatin states on RefSeq genes. Nascent RNA was aligned mm9 on the using TopHat (62). Intergenic emissions were considered as transcribed only when the mapped read yielded a transcription value passing a cut-off of RPKM >0.5. Genomic localizations of intergenic transcripts were evaluated as in (60); briefly, overlap was first evaluated for lincRNAs (as annotated by Ensembl (63)), then for super-enhancers (annotation taken from (64)), and finally with regions within 2 kb from a TSS or a TTS. Remaining regions were annotated as ‘active’ or ‘poised’ based on their K27ac status. All statistical analyses were performed using the corresponding R functions and packages.

RESULTS

Characterization of histone PTMs associated with cis-regulatory regions in macrophages

ChIP followed by MS-based proteomics (ChroP) enables in-depth characterization of the proteomic composition of spatially and functionally defined chromatin regions (29). Here, this approach was further implemented to carry out a comprehensive and unbiased mapping of histone PTMs specifically co-enriched at proximal and distal cis-regulatory regions of mouse macrophage-derived cells RAW 264.7. These regions were enriched prior to MS using as baits H3K4me1 and H3K4me3 as markers of enhancers and promoters, respectively, as well as H3K9me3 as a marker of inactive chromatin regions (as a negative control). Antibody specificity for specific methylation-states and residues within histones was confirmed by using MODifiedTM histone peptide arrays and dot blot analysis (Supplementary Figure S1A, B and Table S3).

Nuclei from asynchronously growing cells were enzymatically digested with micrococcal nuclease (MNase) and used as input for chromatin immunoprecipitations, with each individual N-ChroP experiment performed in three biological replicates. After immunoprecipitation, proteins were separated by SDS-PAGE and bands corresponding to histones H3 and H4 were excised and digested (65). Peptides were analyzed through ultra-high-performance liquid chromatography coupled to high-resolution MS/MS analysis (UHPLC-MS/MS), to enhance the separation and detection of hyper-modified histone peptides (66). Enrichment of each histone modification, and combinations thereof, in the distinct chromatin fractions was evaluated as the ratio of its percentage of relative abundance (RA%) between the ChIP and the input (Figure 1A).

Figure 1.

Distinct chromatin regions isolated by the N-ChroP strategy. (A) Scheme of the N-ChroP experimental approach. (B) Heatmap visualisation of histone PTM relative enrichment (log2) in H3K4me1, H3K4me3 and H3K9me3-immunopurified nucleosomes; modifications are clustered in peptides belonging to histone H3 (up) and H4 (down). Data from three biological independent replicates (I, II and III) are shown. (C) Volcano plots of the histone PTMs that were identified and quantified. P-values and t-test differences were plotted against each other in the volcano plot, comparing ChIPs of H3K4me1 (left) and H3K4me3 (right) with that of H3K9me3. Histone PTMs with P ≤ 0.05 were selected and highlighted as significantly enriched, using two-side t-test and Bonferroni correction. (D) Western blot validation of H3K36me2 (left) and H3K79me2 (right) enrichment in H3K4me1 ChIPs, with intensity signals quantified by ImageJ. The enrichment for each modification in the different ChIPs was calculated as follows: the signal of each PTM in the ChIPs and in the input was divided by the signal of the corresponding unmodified H3, after which the obtained ratio in each ChIP was normalized to the input value.

A significant proportion of the H3K4me1-nucleosomes were captured by ChIP (Supplementary Figure S1C). A similar efficiency was achieved with H3K4me3 and H3K9me3 ChIPs, as previously shown (29). H3K4me1 was chosen for enhancer enrichment on the basis of the well-known H3K4me1Hi/ H3K4me3Low signature commonly associated with enhancers. However, H3K4me1 enrichment is often found associated with other genomic regions (e.g. regions flanking TSS of transcribed genes). Thus, to assess the selectivity of our bait, we compared the ChIP-seq profiles of H3K4me1 and Pu.1, a pioneer TF highly expressed in macrophages that binds to macrophage-specific enhancers and is essential for establishing and maintaining macrophage identity and viability (67). Pu.1 is also essential for stably maintaining H3K4me1 at these regulatory regions and is a well-known marker of enhancers in macrophages, with a much lower occurrence at TSS (14,54,68,69). H3K4me1-positive regions significantly overlapped with Pu.1 binding sites (∼70%; Supplementary Figure S1D). Moreover, only 5% of the fraction of TSS overlapping with H3K4me1-positive regions was excluded when Pu.1 was used as a further discriminant (Supplementary Figure S1E), indicating that the specificity of H3K4me1 for enhancers was comparable to that of Pu.1, with only minor cross-contamination from TSS.

We next profiled the enrichment/depletion of 42 differentially modified histone peptides from H3K4me1, H3K4me3 and H3K9me3 N-ChroP experiments. Unsupervised hierarchical clustering segregated three distinct groups that matched the three baits, with biological replicates of each individual experiment clustering together, underscoring the high experimental and biological reproducibility of the assay (Figure 1B). Notably, an enrichment of H3K4me1 at putative enhancers was complemented by an overall H3K4me3 slight depletion in the same regions, confirming the existence of the H3K4me1Hi/K4me3Low signature (Figure 1B, Supplementary Figure S2A–C) (20). The H3K4me3 ChIP, in contrast, was characterized by a strong enrichment of H3K4me3, as well as of H3K4me1 (albeit to a much lower extent), indicating promoter enrichment (Figure 1B, Supplementary Figure S2A–C). In the H3K9me3 ChIP, the bait was specifically enriched, both alone and in combination with acetylated H3K14, as expected (Figure 1B, Supplementary Figure S2A and B). The observed MS-based enrichment of distinct chromatin regions pertaining to promoters and enhancers was further supported by results from western blots (Supplementary Figure S2D) and quantitative PCR analyses (Supplementary Figure S2E).

The significant depletion of H3K9me2 and H3K9me3 in both H3K4me1 and H3K4me3 ChroP experiments (Figure 1C) is in line with their known accumulation at inactive chromatin regions. Strikingly, both enhancers and promoters were enriched in hyper-acetylated H3 and H4 peptides; more specifically the tetra-acetylated N-terminal tail of H4 was significantly enriched at both regions, while acetylated H3K18, H3K23 and H3K14 (alone or together with H3K9 acetylation) were more enriched at promoters (t-test P ≤ 0.05) (Figure 1C). These data are consistent with the well-described functional link between histone acetylation and transcriptionally active chromatin, suggesting different specificity of the corresponding HATs at the two regulatory regions.

Interestingly, H3K79me1 and H3K79me2 were significantly enriched both at enhancers and promoters (Figure 1C), in line with available evidence indicating that these marks accumulate at actively transcribed genes (70). Moreover, N-ChroP revealed a specific co-association of H3K4me1 with H3K36me2/K36me3 marks (Figure 1C), which are usually associated with transcription elongation due to their accumulation at the 3′ of active genes (71,72). A more in-depth evaluation of the MS data revealed that all H3 peptides containing di- and tri-methylated K36, including K27me1/K36me2, K27me2/K36me2 and K27me1/K36me3, were significantly enriched at H3K4me1 enhancers (Figure 1C). As independent experiments in which ChIP was followed by western blot confirmed these MS results, we conclude that H3K36me2 is specifically enriched at enhancers and that H3K79me2 is found at both enhancers and TSS (Figure 1D).

Unexpectedly, we found that acetylation of H3K27 (K27ac), a known marker of active enhancers (23), was not enriched in chromatin domains containing H3K4me1 (Figure 1B, Supplementary Figure S3A). As N-ChroP selects the entire pool of H3K4me1 enhancers, this lack of enrichment could be due to not having a sufficiently abundant fraction of K27ac-labelled ones for reliable MS detection. However, the observation that this modification is detectable per se by N-ChroP (with a slight enrichment in H3K4me3-positive regions; Supplementary Figure S3B) challenges the detection limit hypothesis. Recent results showing that K27ac associates significantly to TSS may explain this quite unexpected result (73).

H3K27me1 and H3K18me1 were significantly depleted at enhancers, and H3K18me1 was also depleted at TSS (Figure 1B and C). The H3K18me1 enrichment in silent domains—mirrored by its depletion at active domains—had been already observed by N-ChroP (29) and is in line with data supporting its role in silent chromatin, such as its antagonism to H3K18ac and having a maximal half-life lower than that of all histone lysine mono-methylations with active functions (74,75).

Genome-wide mapping of H3K36me2, H3K36me3 and H3K79me2 corroborates N-ChroP data

The most unexpected observation from systematic histone PTM analysis by N-ChroP was the co-association of H3K36me2/me3 and K79me2 with H3K4me1, which we therefore further investigated with genome-wide localization analyses. We validated the association of H3K36me2, H3K36me3 and H3K79me2 with H3K4me1-nucleosomes by ChIP-seq using antibodies whose specificity was carefully assessed (Supplementary Figures S1A, B and S3B, Supplementary Table S3): genomic distribution of these marks showed that H3K36me2 (in a total of 50 748 peaks) frequently associates with intronic and extragenic regions (64% and 34%, respectively), with a pattern very similar to H3K4me1. In contrast, H3K36me3 (in a total of 28 167 peaks) and H3K79me2 (in a total of 13 412 peaks) accumulate mainly at introns (83% and 84%, respectively). As a control, ChIP-seq profiles of the three baits showed the well-characterized genomic localizations (Figure 2A).

Figure 2.

ChIP-seq analysis validates the MS-based histone PTM profiling. (A) Global genomic distribution of peaks for H3K36me2, H3K36me3, H3K79me2 and histone PTMs used as baits in N-ChroP, at the basal state. TSS are defined as regions ±1 kb around known transcription starting sites in the RefSeq database. (B) Heatmap representing the chromatin state model generated by ChromHMM for histone PTMs at the basal state. (C) Relative quantification of the overlap between H43K4me1 and other histone modifications based on chromatin state identification. Overlaps between H3K4me1 and H3K36me2 (26% of H3K4me1-containing chromatin states), H3K36me3 (10%) or H3K79me2 (24%) are highlighted.

We next used the ChromHMM segmentation approach, which identifies significant genome-wide co-associations of histone marks (chromatin states) through unsupervised clustering of different histone PTM profiles (57). Analysis of all ChIP-seq profiles in resting macrophages generated known chromatin states—such as the H3K4me1/K4me3 combination that marks regions flanking active TSS—and confirmed that H3K4me1/me3 and H3K9me3 do not overlap. Moreover, we identified regions in which H3K4me1 was co-enriched with H3K36me2/K36me3 and H3K79me2 (Figure 2B), confirming the patterns that emerged from N-ChroP. In addition, frequency analysis of chromatin states shows that H3K4me1-rich regions had 26%, 10% and 24% overlap with H3K36me2-, H3K36me3- and H3K79me2-positive regions, respectively, further highlighting the sensitivity of N-ChroP in detecting PTM co-associations with relatively low abundance (Figure 2C).

To assess how these marks distribute within H3K4me1-enriched chromatin states, we next assessed their cumulative distributions around Pu.1 bound regions (which are in H3K4me1-enriched chromatin states) (Figure Supplementary S4A). While H3K4me1 is enriched in the nucleosomes flanking the nucleosome-depleted regions bound by Pu.1 (14), H3K36me2/K36me3 and H3K79me2 have a bimodal pattern that centers around the nucleosome-free region, with a maximum enrichment located further downstream and upstream of Pu.1 bound regions (Figure 3A, B, Supplementary S4B). Hence, we concluded that these modifications can co-occur with H3K4me1 to define a broader region surrounding Pu.1 peaks.

Figure 3.

Genome-wide distribution of H3K36me2, H3K36me3 and H3K79me2 around Pu.1. (A) Reads distribution around Pu.1 binding sites that overlap with enhancers marked by H3K36me2, H3K36me3 or H3K79me2. Chromatin states #2, #11 and #5 indicate regions in which H3K4me1 co-localized with H3K36me2, H3K36me3 or H3K79me2, respectively. (B) Representative ChIP-seq signal snapshots showing the distribution of H3K36me2, H3K36me3 and H3K79me2 relative to the H3K4me1 and Pu.1 profiles.

Dynamic profiling of histone modification changes at cis-regulatory regions upon inflammatory stimulation

In macrophages, inflammatory stimuli elicit extended transcriptional activation, during which synchronized changes in gene expression are caused by the combined activity of cell-specific pioneer and stimulus-dependent chromatin determinants at cis-regulatory regions of inflammatory genes (16,76,77). The expression profiles of a panel of inflammatory genes at different times after LPS stimulation showed that most of them were already significantly transcriptionally activated at 1 h, with peaks at 4 h after LPS (Supplementary Figure S5A, B). To cover the acute phase of the inflammatory transcriptional response, we used the time points of 0, 1 and 4 h after LPS treatment for our dynamic chromatin proteomic experiments in RAW 264.7 cells (Figure 4A). Consistent with previous ChIP-seq evidence, the level of H3K4me1 in the immunoprecipitated chromatin remained constant across the three time points (Figure 4B and C), confirming that it is an ideal bait for profiling co-associated modifications over time.

Figure 4.

Dynamic profiling of histone PTMs at enhancers during the inflammatory response. (A) Experimental design for the N-ChroP time-course using H3K4me1 as bait. UT, untreated cells; 1h or 4h, cells treated with lipopolysaccharide (LPS) for 1 or 4 h, respectively. (B) Heatmap visualization of the dynamic relative enrichment of histone PTMs during LPS treatment (log2). Data are shown as a mean of three independent biological replicates (n = 3). (C) Line graphs for a subset of histone PTMs, representing the fold enrichment over time. UT, 1 h and 4 h (x-axis) correspond to cells untreated or treated with LPS for 1 h or 4 h, respectively. (D) Left histograms represent the log2 relative enrichment of peptides containing K36me2 and K36me3 with K27me in the H3 27–40 aa peptide. The right histogram indicates the relative enrichment of K79 modifications in the H3 73–83 aa peptide. Values are expressed as mean ± s.e. (standard error) from three independent biological experiments (n = 3), before and after LPS treatment. (E) Global genomic distribution of peaks for H3K36me2, H3K36me3 and H3K79me2 after 1 h or 4 h (upper and lower panels, respectively) of LPS treatment. TSS are defined as regions ±1 kb around known TSS in the RefSeq database. (F) Representative ChIP-Seq signal snapshots showing the distribution of H3K4me1, H3K36me2, H3K36me3 and H3K79me2, before and after LPS treatment.

The analysis of 42 differentially modified histone peptides at the three time points showed a mild increase of global histone H3 and H4 acetylation, in line with the model of hyper-acetylation leading to chromatin relaxation, which in turn increases the accessibility of stimulus-specific TFs (Figure 4B and C) (78,79). A more composite picture, however, emerged after an in-depth dissection of distinct acetylation sites: some acetylated histone peptides (e.g. tri- and tetra-acetylated isoforms of the H4 tail) increased during inflammation, whereas others (e.g. H3K18ac/K23ac) remained constant (Figure 4B and C). This result could indicate that the constitutively acetylated state is characterized by marks different than the LPS-induced acetylated state, as a consequence of the activity of distinct HATs. H3K4me3, a prototypic mark of promoters of actively transcribed genes (20,80), also increased after LPS treatment (Figure 4B and C, Supplementary Figure S6A). Interestingly, H3K4me2 followed the same trend, suggesting that, in this context, this mark may have a role more similar to H3K4me3 than to H3K4me1.

We observed an interesting decrease for the silencing marks H3K9me3 and H3K27me2/me3, and for H3K36me1 (Figure 4B and C), corroborating previous evidence indicating a synergy between H3K36me1 and H3K27me2 in establishing gene silencing (81).

The overall levels of H3K36me2/me3 and H3K79me2 remained stable during the early inflammatory response (Figure 4B and C), displaying a behavior similar to H3K4me1; this was confirmed by western blot analysis (Supplementary Figure S6A). ChIP-seq profiles of H3K36me2/K36me3 and H3K79me2 at the three time points further confirmed that the overall abundance and genomic localization of these marks was stable (Figure 4D and E, Supplementary Figure S6B and C), suggesting that these marks synergize with H3K4me1 to maintain specific epigenetic signatures within distinct enhancer subpopulations during the acute inflammatory response.

The co-association of H3K4me1 with H3K79me2 and H3K36me3 identifies distinct intronic enhancer subpopulations

H3K79me2 and H3K36me3 are typically enriched at transcribed genes, with the former mainly localized in proximity of TSS and the latter accumulating at the 3′ end of gene bodies (71,82). Accordingly, histone PTM analysis displayed the enrichment only of H3K79me2 in H3K4me3-enriched chromatin, while both were over-represented in H3K4me1-containing nucleosomes (Figure 1B). Chromatin state analysis identified three chromatin states: H3K4me1/K79me2, H3K4me1/K4me3/K79me2, and H3K4me3/K79me2 (chromatin states #5, #6 and #7, respectively) (Figure 2B). While chromatin states #6 and #7 accumulated around promoters (but located differentially relative to TSS), the chromatin state #5 was found mainly within introns (Figure 5A and B). Interestingly, the chromatin state characterized by H3K4me1/K36me3 (#11; Figure 2B) was over-represented at introns (Figure 5A and B), suggesting that multiple subpopulations of intronic enhancers exist and are marked by distinct combinations of H3K4me1 with other PTMs.

Figure 5.

Characterization of different intronic enhancer subpopulations. (A) Distribution and enrichment of chromatin states in a 10-kb window surrounding TSS. Chromatin states #5 and #6 correspond to regions in which H3K4me1 co-localized with H3K79me2 alone or with H3K4me3, respectively. Chromatin state #7 represents regions in which H3K79me2 was co-enriched with H3K4me3. The chromatin state #11 represents regions in which H3K4me1 and H3K36me3 co-localized. (B) Global distribution of chromatin states #5, #6, #7 and #11. (C) Heatmap representing the chromatin state model at basal level and after 1 h or 4 h of LPS treatment, including ChIP-seq profiles for H3K27ac and H3K27me3. (D) Distribution of RNA quantification values, either from total (left and central panel) or nascent (right panel) RNA analysis, associated with different sets of chromatin states. Chromatin states #7, #8, #14 and #17 represent regions in which H3K4me1 co-localized with K79me2, K27ac/K79me2, K27ac/K36me3 or K27ac, respectively. Chromatin states #16 and #18 represents instead regions in which H3K4me1 is present alone or in combination with H3K27me3. (E) TF binding motifs enriched for enhancers marked by H3K4me1/K36me3 and H3K4me1/K79me2.

To gain insight into the role of H3K79me2 and H3K36me3 in regulating the respective subpopulations of intronic enhancers, we repeated the chromatin state analysis but including all ChIP-seq profiles from both unstimulated macrophages and those stimulated with LPS for 1 or 4 h, to obtain a single chromatin model applicable to all conditions. We also added the H3K27ac and the H3K27me3 marks to this chromatin model, to help discriminating poised and active/non-active enhancers (Supplementary Figure S7A and B). The newly generated clusters (Figure 5C) revealed the presence of two groups of H3K4me1/K79me2 enhancers: those with or without H3K27ac, indicated as active or non-active H3K79me2 intronic enhancers (states #8 and #7, respectively). Intriguingly, the chromatin state marked by H3K4me1/K36me3 was identified only with H3K27ac (state #14), indicating that enhancers marked by K36me3 are more frequent in an active state. In contrast, H3K27me3 clustered with H3K4me1, but not with H3K79me2 or H3K36me3. After LPS treatment, chromatin states marked by H3K4me1/K79me2 or H3K4me1/K36me3 did not change their distribution (Supplementary Figure S7B), confirming that the inflammatory response did not cause major changes in genomic distribution of these PTM combinations.

We then intersected the newly generated chromatin states with RNA-seq expression data from resting macrophages to assess the contribution of these enhancer subpopulations to the expression of the underlying genes. As expected, transcript levels were significantly higher for the genes marked by either H3K79me2 or H3K36me3 with respect to all genes (P < 0.001; two-tailed Z-test) (Figure 5D, left panel), confirming them as markers of active transcription. As an additional control, we assessed the level of transcription of genes containing any of the H3K4me1-positive chromatin states. As expected, while they displayed significantly higher transcription compared to all genes, their overall expression was lower than that of genes marked by either H3K79me2 or H3K36me3 (Figure 5D, left panel). Transcript levels that associated with genes marked by H3K27me3 were significantly lower than the expression level of all genes (Figure 5D, left panel). When we assessed the effect of intronic enhancers on transcription, no significant differences (P > 0.2, two-tailed Z-test) were seen between genes containing H3K4me1/K79me2 or H3K4me1/K36me3 regions and those with H3K79me2 or H3K36me3 regions. Similarly, no significant differences were detected when taking into account the activation state of these intronic enhancers based on their co-association with H3K27ac (Figure 5D, middle panel).

Since steady-state transcripts do not necessarily reflect on-going transcription rates, we profiled nascent RNAs from 4-thiouridine (4sU) labeled cells (60) and discovered a significantly higher transcriptional activity (P < 0.001; two-tailed Z-test) at genes with H3K79me2- or H3K36me3-positive active enhancers (Figure 5D, right panel). This suggests that active intronic enhancers marked by H3K4me1/K79me2 and H3K4me1/K36me3 might not simply have a physical connection with their host genes but rather also play a functional role in their transcriptional regulation. The transcriptional activity of poised enhancers, in contrasts, was lower when both total and nascent transcript levels were compared (Figure 5D, middle and right panels), in line with previous observations (83).

To assess a potential role of these newly identified classes of intronic enhancers in controlling inflammatory gene expression, we interrogated the underlying genomic regions for enrichment of binding motifs for TFs. Interestingly, we found specific enrichment of motifs recognized by the NF-κB and IRF proteins (Figure 5E), whose co-occurrence at a subset of regulatory regions has been previously described (84,85). In addition, we retrieved motifs for various TFs known to regulate inflammatory gene expression (e.g. AP-1, Jun, Fra1 and Fosl2) or to enhance transcriptional activity (e.g. ETS) (86), and the enrichment of these TFs was more prominent after LPS stimulation (Figure 5E, Supplementary Table S5). Hence, our data indicated that these two subpopulations of enhancers have a role in the transcriptional activation of the associated genes and contain binding sites for relevant inflammatory TFs.

Prompted by these findings, we carried out a gene ontology (GO) enrichment analysis of the genes hosting the two enhancer subpopulations marked by H3K4me1/K27ac/K79me2 and H3K4me1/K27ac/K36me3, in both basal and inflamed states. Several GO terms associated to immunity and inflammatory response were significantly enriched initially and were further overrepresented after LPS treatment (Figure 6A), strengthening the hypothesis that H3K79me2- and H3K36me3-positive active enhancers increase the transcriptional activity of genes involved in inflammation-related pathways. In particular, while more generic GO categories such as ‘activation of MAPK activity’ and ‘immune response regulating signal transduction’ are common to both groups, more specific biological functions are associated with specifically with a mark; thus, activation and ‘positive regulation of JUN kinase activity’ was uniquely associated to H3K4me1/K27ac/K79me2 enhancers, while ‘response to lipopolysaccharide’ was over-represented in the H3K4me1/K27ac/K36me3 enhancers. None of these GO terms were enriched in control genes marked by H3K4me1/K27ac only (referred to as class I; (83)), corroborating the assumption that intronic enhancers marked by the H3K4me1/K79me2 and H3K4me1/K36me3 signatures enforce specific transcriptional programs within the inflammatory response. In contrast, genes marked by H3K4me1/K27me3 (referred to as class II; (83)) were linked to homeostasis (Figure 6A).

Figure 6.

Functional analysis of genes associated with different populations of intronic enhancers. (A) Gene ontology (GO) analysis and enriched pathways in genes hosting K4me1-, K4me1/K79me2- or K4me1/K36me3-marked active enhancers, at basal state or after LPS stimulus. Significantly enriched GO terms are indicated both by colour code and the P-value shown within each cell. (B) Conversion of intronic enhancers after LPS induction: schematic representation of class II enhancers (K4me1/K27me3) that convert to class I enhancers (K4me1/K27ac) after LPS stimulus. GO enrichment of the novel class I enhancers marked by H3K4me1/K27ac and H3K4me1/K27ac/K36me3 is shown.

Enhancers that convert genes from a poised (class II) to an active (class I) state play a crucial role during the transcriptional changes induced in response to stimuli or that occur during functional and developmental transitions (83); we therefore focused on the H3K4me1/K27me3 regions that acquired H3K27ac after 4 h of LPS treatment. We observed that about 2 000 enhancers underwent this conversion, with ∼5% of them displaying H3K36me2 or H3K79me2 as an additional mark (Figure 6B). In particular, GO analysis of the genes with class I H3K36me3 enhancers showed an enrichment of terms related to inflammatory pathways, including that for I-NF-κB. Although GO analysis of the genes associated to class I H3K79me2 enhancers revealed no specific enrichment of terms or categories, various interesting genes were included in this group, such as: MAPK14, a Ser/Thr kinase induced by LPS with a well-characterized role in inflammation (reviewed in (87)); SH3PXD2B, encoding a protein required for podosome formation (88); and Trim46 and Tpm1, the proteins of which appear to distinguish functionally different macrophage subpopulations (89).

The H3K36me2/K4me1 combination as a novel marker of enhancer transcription

The presence of short (<2 kb) non-coding RNAs transcribed bi-directionally from enhancers (eRNAs) has been recently correlated to enhancer activity and expression of nearby genes (90,91). The combination of different molecular features, such as eRNAs and histone modifications (e.g. H3K27ac) has been previously used to identify active enhancers (23,92). A recent study used a logistic regression model with selected histone PTMs to reveal modification patterns predictive of eRNA transcription (93). Nevertheless, the functional correlation between specific PTM patterns and eRNA transcription has not been systematically addressed yet.

Our N-ChroP data showed that H3K36me2 is enriched at H3K4me1-chromatin regions; as this modification is generally associated with transcriptional elongation (94–96), we hypothesized that transcriptional events occurred at regions marked by this pair of histone marks. We tested this hypothesis using our chromatin state model to identify extragenic regions marked either by the combination H3K4me1/K36me2 (#2 and #3; 8757 regions), or H3K4me1 alone (#16 and #17; 17 941 regions). Both groups were further divided as ‘active’ or ‘inactive’ based on the co-association of H3K27ac (Figure 7A and B). When these regions were interrogated for the presence of actively ongoing transcription, a significant increase was observed for nascent RNA reads from intergenic enhancers marked by H3K4me1/K36me2, as compared to reads from regions marked by H3K4me1 only (P < 0.0001, Mann–Whitney test) (Figure 7A). After applying a RPKM cutoff of 0.5, we classified the two sets of transcribed extragenic enhancers on the basis of annotated genomic features, focusing on super-enhancers (64) and active enhancers (positive for H3K4me1 and H3K27ac; (23)) as these elements are actively transcribed. Indeed, both elements showed a statistically significant enrichment in H3K36me2-marked enhancer regions, as compared to the H3K4me regions (Figure 7C).

Figure 7.

The combinatorial mark H3K36me2/K4me1 associates with transcribed enhancers. (A) Comparison between the number of intergenic enhancers marked by the H3K4me1/K36me2 combination (states #2 and #3) or by H3K4me1 alone (states #16 and #17), as identified by the 27-state ChromHMM model. Within these classes, active or non-active enhancers were also divided based on the presence or absence of H3K27ac, and the distribution of transcript levels associated with the two groups is shown. (B) Examples of representative genomic regions displaying H3K4me1, H3K36me2, H3K27ac and nascent RNA enrichment. (C) Distribution of transcribed extragenic enhancers belonging to the K4me1 and K4me1/K36me2 groups, as compared to annotated genomic features; the enrichment of active and super-enhancers in the H3K4me1/K36me2 group is highlighted. (D) Distribution of H3K4me1 (left) and H3K36me2 (right) normalized read counts in H3K4me1/K36me2-marked enhancers in cells untreated or LPS-stimulated, with or without flavopiridol (FP) pretreatment. (E) mRNA levels of NSD1 and NSD2 measured by RT-qPCR analyses in RAW 264.7 cells expressing the scrambled (shScramble) and NSD1- and NSD2-specific shRNAs (shNSD1 and shNSD2). Gene expression is normalized to TBP mRNA level. Values are expressed as mean ± s.e. (standard error) of two biological replicates (n = 2). (F) NSD1 protein level in cells expressing scrambled (shScramble) and NSD1-specific shRNA (upper panel). Protein LFQ ratio distribution from label-free quantitative (LFQ) MS-analysis of whole extracts from cells expressing NSD2-specific shRNAs (shNSD2) versus shScramble cells. NSD2 LFQ ratio value is indicated with the green dot and equals log2 = –1.4, confirming its down-regulation. (G) Western blot analysis of H3K36me2 from protein extracts from cells expressing scrambled (shScramble) and NSD1- and NSD2-specific shRNAs (shNSD1 and shNSD2). Vinculin and Ponceau staining are used as loading control. (H) qPCR of extragenic regions marked by H3K4me1/K36me2 after H3K36me2 ChIP in cells expressing scrambled (shScramble) or NSD1- and NSD2-specific shRNAs (shNSD1 and shNSD2), at 4 h after LPS treatment. Data are expressed as percent of the input. (I) Quantitative RT-qPCR using 4sU-labeled RNAs on extragenic regions marked by H3K4me1/K36me2 in cells expressing scrambled (shScramble) and NSD1- and NSD2-specific shRNAs (shNSD1 and shNSD2), at 4 h after LPS treatment. Values are expressed as mean ± s.e. (standard error) of two biological replicates (n = 2).

Recent evidence proposes that the RNA pol II actively contributes to establishing and maintaining characteristic histone PTM signatures at enhancers, through their transcription (43). To further investigate the emerging cross-talk between enhancer transcription and H3K36me2, we profiled H3K4me1 and H3K36me2 by ChIP-seq after flavopiridol treatment, a drug that inhibits elongation by blocking Cdk9 (97), both at basal conditions and after LPS stimulation (Supplementary Figure S8A). Transcriptional inhibition correlated with the impairment of H3K4me1 deposition (in line with previous reports; (43)) in regions marked by H3K4me1, either alone or in combination with H3K36me2 (Figure 7D, left panel, Supplementary Figure S8B). Interestingly, we observed that H3K36me2 decreased at intergenic regions, both in basal and LPS-induced states, when transcription elongation was inhibited (Figure 7D, right panel). Moreover, the reduced deposition of H3K4me1 was more prominent in extragenic enhancers when it co-localized with H3K36me2 (Figure 7D, left panel), supporting the idea that K36me2 is an additional mark of transcribed enhancers and that its deposition at these regulatory regions requires transcription.

As a complementary strategy, we also assessed the impact of H3K36me2 deposition on enhancer transcription during inflammation. In mammalians, SETD2 is known to tri-methylate H3K36, while various enzymes display mono- and di-methylase activity at this residue, such as NSD1, NSD2 and NSD3, which constitute the NSD family of K-methyl-transferases and were shown to mono- and di-methylate K36 in vivo (98–100). Based on this enzymatic activity and in order to tackle H3K36me2 level, we silenced NSD1 and NSD2, individually or in combination, in RAW 264.7 cells, unstimulated or stimulated with LPS for 4 h in two biological replicates. Enzymes knock-down was achieved with lentivirus-mediated short hairpin RNA interference (shRNA) (Figure 7E and F). Bulk H3K36me2 reduction was assessed by western blot (Figure 7G) and further confirmed by H3K36me2 ChIP, followed by qPCR, at a pool of extragenic regions previously identified as marked by H3K4me1/K36me2 (Figure 7H). Intriguingly, we observed that after LPS treatment nascent transcripts were reduced in samples depleted for NSD1 and NSD2 at about 40% of the regions assessed (Figure 7I). This reduction seems to be more pronounced after NSD1 depletion and in the combined knock-down than in NDS2-depleted cells, suggesting a more prominent role of NSD1 in the deposition of K36me2 at these transcribed enhancers.

Overall, our results indicate that -at least in a subset of enhancers- H3K36me2 deposition is interlinked with their transcription, possibly serving as docking site for transcriptional co-regulators.

DISCUSSION

Here, we performed a widespread characterization of histone PTMs associated to cis-regulatory regions in macrophages, focusing on enhancers that are poorly characterized from a molecular point of view despite their established role in sustaining and coordinating cell- and stimulus-specific transcriptional programs. As bait for our N-ChroP proteomics approach, we chose three histone marks with distinct genomic locations and functions: H3K4me1 and H3K4me3 mark enhancers and promoters, respectively (20), and H3K9me3 marks heterochromatin (used as a negative control). The enrichment of specific regulatory regions by multiple histone modifications has been recently used for the proteomics identification of interactors (101), but a systematic MS-characterization of histone modifications at these genomic elements has never previously been addressed.

We profiled 42 differentially modified peptides for H3 and H4, with the unique possibility of inferring functional PTM associations, not only within the same histone molecule, but also among different histones belonging to the same nucleosome, thanks to the enrichment of intact nucleosomes. Analyzing the modification patterns more specifically associated with H3K4me1 and H3K4me3 chromatin, we observed that H3K79me2 was enriched at both cis-regulatory regions, while H3K36me2 and H3K36me3 accumulated only at H3K4me1-enhancers. Because of their typical association with transcribed gene bodies, the presence of H3K36me2/me3 and H3K79me2 at enhancers was unexpected and prompted us to investigate it further. ChIP-seq analysis confirmed the co-association of H3K4me1 with these marks, and hinted at the existence of different enhancer subpopulations marked by distinct PTM signatures. The subsequent intersection of the resulting chromatin states with transcriptomics data permitted us to distinguish different enhancer classes that have distinct features and transcriptional activity. In particular, we identified two groups of intronic enhancers marked by either H3K4me1/K79me2 or H3K4me1/K36me3, which associated with highly transcribed genes involved in different transcriptional programs of macrophage inflammation.

We propose that these intergenic enhancers modulate the expression of the nearby genes. This still needs to be functionally corroborated with further experiments; for instance, DNA editing could be used to selectively delete putative intronic enhancers and then assess the transcriptional response upon inflammation; alternatively, chromatin conformation capture approaches to determine the three-dimensional organization analysis of the genomes around these regulatory regions may reveal interactions between putative enhancer sub-groups and promoters, thereby helping to identify regulated genes.

By integrating nascent transcript analysis with ChIP-seq and proteomics data, we discovered that the H3K4me1/K36me2 combination tags actively transcribed enhancers and, in particular, super-enhancers, and that H3K36me2 deposition at these regions is linked to their active transcription. Super-enhancers are defined as clusters of transcriptional enhancers densely bound by TFs and other co-regulators, which play a critical role in defining cell identity by regulating cell type–specific genes located in close proximity (64,102). Although H3K27ac accumulates at these elements and was initially considered as a hallmark of super-enhancer, this mark alone is in fact insufficient to identify these elements comprehensively and unambiguously. Given that the presence of eRNAs marks enhancer activation in macrophages (103), our finding suggests that the newly identified H3K4me1/K27ac/K36me2 signature may be used for de novo identification of active super-enhancers. Interestingly, we also found that H3K36me2 deposition seems to be linked to the transcription of a subset of extragenic regions marked by the dual K4me1/K36me2 signature, suggesting a possible role of this mark in enforcing active enhancer transcription.

Quite unexpectedly, H3K27 acetylation showed a decreasing trend in ChroP over the time course after LPS. This result may be explained by the fact that macrophages react to inflammation with a composite response comprising both activation of ‘poised’ enhancers that regulate inflammatory genes and silencing of active enhancers of housekeeping or survival genes; hence, downregulation of K27ac may result from an overall reduction in the number of active enhancers within the pool of the precipitated H3K4me1-nucleosomes (54). This reveals a possible limitation of the current N-ChroP setup when applied to the dynamic profiling of low abundance PTMs: because N-ChroP enriches the complete pool for all possible H3K4me1-positive enhancers, the MS-based readout ‘averages’ among histone modification patterns of the bulk of enhancers, which may hide the dynamic behavior of rare modifications that mark only a small fraction of enhancers. Re-ChIP experiments, in combination with MS, could offer a solution to overcome this constraint, although MS sensitivity remains a major challenge.

In sum, our results with N-ChroP have confirmed its strong potential for discovering novel combinatorial histone PTMs that mark functionally distinct genomic regions. While the MS findings still need to be complemented with next-generation sequencing data (e.g. ChIP-seq and RNA-seq) to define the genome-wide distribution of these combinatorial marks and to assess their functional outcome, integrating proteomics with ChIP-seq data can improve enhancer identification and provide hints about unexpected cross-talk between different chromatin determinants (such as between eRNAs and histone modifiers, as shown here).

DATA AVAILABILITY

Sequencing data have been deposited in the NCBI’s Gene Expression Omnibus (GEO) under GEO SuperSeries accession number GSE91009. The mass spectrometry proteomics data have been deposited to the ProteomeXchange Consortium (104) via the PRIDE partner repository with the dataset identifier PXD007582.

Click here for additional data file.