Combined Metabolic and Chemical (CoMetChem) Labeling Using Stable Isotopes—a Strategy to Reveal Site-Specific Histone Acetylation and Deacetylation Rates by LC–MS

Histone acetylation is an important, reversible post-translational protein modification and a hallmark of epigenetic regulation. However, little is known about the dynamics of this process, due to the lack of analytical methods that can capture site-specific acetylation and deacetylation reactions. We present a new approach that combines metabolic and chemical labeling (CoMetChem) using uniformly 13C-labeled glucose and stable isotope-labeled acetic anhydride. Thereby, chemically equivalent, fully acetylated histone species are generated, enabling accurate relative quantification of site-specific lysine acetylation dynamics in tryptic peptides using high-resolution mass spectrometry. We show that CoMetChem enables site-specific quantification of the incorporation or loss of lysine acetylation over time, allowing the determination of reaction rates for acetylation and deacetylation. Thus, the CoMetChem methodology provides a comprehensive description of site-specific acetylation dynamics.

(B) Chemical structure of MS-275 (also known as Entinostat). (C) Workflow for the metabolic labeling of the CoMetChem approach to study histone acetylation dynamics upon HDAC inhibitor (HDACi) treatment. RAW264.7 cells were first cultured in a [U-12C]-Glucose containing medium for 22 hours, followed by a preincubation of 16 hours with either a histone deacetylase inhibitor (HDACi, MS-275 (c= 1 µM) or SAHA (c= 0.41 µM)) or a carrier (DMSO, c= 0.01%). After pre-incubation, the culture medium was replaced by a [U-13C]-Gluc containing medium including either HDACi or carrier. The nuclei were isolated from the samples at different time points, followed by histone extraction, chemical derivatization of unmodified lysine residues at the protein level using 13 C6,D6-acetic anhydride, tryptic digestion and quantitative LC-MS analysis.   Linear fit of the measurements from samples taken within 8 hours after addition of [U-13C]-Glucose used to calculate the acetylation rates of the single-acetylated H3(18-26) species with an acetyllysine at K18ac (A) and K23ac (B), the deacetylation rates at K23ac (C) and K18ac (D); the acetylation rates of the double-acetylated H3(18-26) species with an acetyllysine at K23ac (E) and K18ac (F), and the deacetylation rate of the doubleacetylated H3(18-26) species (G).  The table lists the reaction rates and the corresponding standard error and R 2 values calculated using the initial flux fitting method described in the Materials and Methods section. The naming of acetylation sites and isotopologues has been done according to Figure 1-3. The reaction names are as indicated in Figure 3 and Supplementary Figures S6 and S7. The fourth column refers to the experiment the data are derived from (C for control, MS-275 and SAHA for the inhibitor treatments). The table lists the reaction rates and the corresponding standard error and R 2 values calculated using the initial flux fitting method described in the Materials and Methods section. The naming of acetylation sites and reactions are as indicated in Figure S6 A. The fourth column refers to the experiment the data are derived from (C for control, MS-275 and SAHA for the inhibitor treatments).

Histone extraction
For histone extraction, the cell nuclei were isolated first. The cells were collected and washed three times with PBS. For lysis, 500 µL ice-cold nuclear isolation buffer (15mM Tris-HCl, 15 mM NaCl, 60 mM KCl, 5 mM MgCl2, 1 mM CaCl2, 250 mM sucrose, 1mM dithiothreitol (DTT), 1mM phenylmethylsulfonyl fluoride (PMSF), 10mM nicotinamide, and 10mM sodium butyrate, 0.01% sodium deoxycholate(SDC)) was used. The homogenized cells were incubated for 10 min on ice and washed three times with the nuclear isolation buffer without SDC. Histones were extracted from the nuclei using 0.2 M hydrogen chloride (HCl) at a cell density of 8 × 103 cells/μL and samples were incubated on a rotator at 4°C for 30 min. Afterwards, the samples were centrifuged (16 000× g, 10 min, 4°C). The supernatants were transferred to a new reaction tube and trichloroacetic acid was slowly added to a final concentration of 25% (w/v). The reaction tubes were inverted five times and incubated for 35 min on ice. The histones were pelleted by centrifugation (16 000 × g, 10 min, 4°C) and washed three times with ice-cold acetone. The samples were air-dried in a fume hood for 30 min at RT. The histone pellets were dissolved in 50 μL H2O (HPLC grade). Protein concentration was determined using the microplate BCA protein assay kit following the manufacturer's instructions. A BSA standard (2 mg/mL) was used to calibrate the assay across the concentration range of 2-40 µg/mL. The absorbance was measured at 580 nm using a plate reader.

Chemical acetylation using 13C4,D6-acetic anhydride and tryptic digestion to investigate chemical acetylation efficiency
To study the chemical acetylation efficiency, 6 µg of the histone extract were diluted in HPLC-H2O to a final volume of 50 µL. For reduction, the samples were incubated with 10 mM DTT (dissolved in 100 mM sodium borate (SB, pH 8.5), 100 mM ammonium bicarbonate (ABC, pH 8.5) or triethylammonium bicarbonate (TEAB, pH 8.5) for 10 min at 57°C on a shaker (600 rpm). Afterwards, the samples were alkylated with 20 mM iodoacetamide (IAA, dissolved in 100 mM SB, 100 mM ABC or 100 mM TEAB, pH 8.5) and incubated for 30 min at room temperature in the dark. Subsequently, free IAA was quenched by adding 10 mM of DTT (dissolved in 100 mM SB, 100 mM ABC or 100 mM SB, pH 8.5). Subsequently, 50 µL of 100 mM SB (pH 8.5), 100 mM ABC (pH 8.5) or 100 mM TEAB (pH 8.5) were added. Then 1 µL of 4.1 mM 13C4,D6-acetic anhydride (13C4,D6-AA, dissolved in water-free DMSO) was added, followed by incubation on a thermomixer (1000 rpm) at 25°C for 10 minutes. The chemical acetylation was repeated twice. O-acetylation was reverted by adding hydroxylamine to final concentration of 7.6 mM, followed by an incubation on a thermomixer (450 rpm) at 25°C for 120 minutes. For tryptic digestion, the samples were incubated with 1 µL of a trypsin solution (c = 0.2 µg/µL sequencing grade modified trypsin, dissolved in trypsin resuspension buffer, Promega, Walldorf, Germany) for 16 hours at 37°C. Afterwards, the samples were acidified by adding formic acid (FA) to final concentration of 0.1% FA. The searches were performed using the following parameters: precursor mass tolerance was set to 6 ppm and fragment mass tolerance was set to 20 ppm. For peptide identification, five missed cleavages were allowed. Carbamidomethylation of cysteine residues was considered as a fixed modification. Oxidation of methionine residues and an acetylation of lysine residues (light: H(2) C(2) O, heavy: D(3) C(2) O) were allowed as variable modifications. Peptides and proteins were identified with an FDR of 1%. For analysis of chemical acetylation efficiencies, peptide intensities of the "modification_specific_peptide" output file were used. Statistical analysis was performed with GraphPad Prism (version 5) using two-tailed unpaired t-test.

LC-MS/MS analysis of CoMetChem labeled histone species
For LC-MS/MS analysis, 100 ng of the tryptic peptide digests were injected on a nano-ultra pressure liquid chromatography system (Dionex UltiMate 3000 RSLCnano pro flow, Thermo Scientific, Bremen, Germany) coupled via an electrospray ionization (ESI) source to a tribrid orbitrap mass spectrometer (Orbitrap Fusion Lumos, Thermo Scientific, San Jose, CA, USA

MS2 fragment ion quantification using PASTAQ
Raw-files were converted to mzML-files using msconvert (http://proteowizard.sourceforge.net/tools.shtml) and MS2 spectra were extracted using PASTAQ's [1] Python bindings. The resulting MS/MS spectra were filtered to only contain those with precursor m/z and retention time within the given min/max retention time window for each inclusion list entry. Each MS/MS spectrum was searched for the theoretical fragments of interest by extracting a slice of data centered around the theoretical m/z. Finally, the fragments were quantified by performing a weighted Gaussian curve fitting using non-linear least squares and the area under the curve was calculated using the trapezoidal rule and exported together with the raw total intensity, fitted height, fitted m/z, and fitted width of the peak. For site-specific quantification of isobaric H3(18-26) isotopologue species, the relative MS2 abundance of the isotopologue species is multiplied with the corresponding MS1 abundance. For site-specific quantification of the single-acetylated H4(4-17) species, the intensities for the y5, y7, and y12 ions of the single-acetylated precursors were quantified and the relative abundances were calculated according to the approach of Feller et al. [2] as follows:

Supplemental Tutorial -CoMetChem data analysis pipeline
In the following section, a detailed tutorial on bioinformatics workflow ( Figure S11) and data processing can be found. The different steps of the pipeline are exemplarily explained by the H3(18-26) peptide KQLATKAAR.

Generation of extracted ion chromatograms (EICs) of the H3(18-26) isotopologue species
Extracted ion chromatograms (EICs) of the isotopologue species have to be generated for the analysis.  Figure S13) the following parameters are used in the "chromatogram ranges" tab: Detector Type: MS Trace Type: Mass Range Ranges: theoretical m/z of the corresponding isotopologue species Mass Tolerance: 3 ppm (Note: The mass tolerance should be adjusted according to the measurement accuracy). 1.4 Export the workspace of peak detection including peak heights and peak areas of the EICs as *.csv file (see example Table S4).

Natural stable isotope correction using PICor
To quantify the different isotopologue species, a correction of measured isotope intensities for naturally occurring stable isotopes is performed using PICor. The python PICor package and source code is available under https://github.com/MolecularBioinformatics/PICor. The detailed approach of PICor is described on bioRxiv [3] .
2.1 Read in *.csv file containing the intensities of the EICs of the isotopes of the H3(18-26) isotopologue species MS1 data (see section 1). The columns containing intensities have to be renamed so they contain information about the number of isotopically labeled atoms (e.g. 2C13 3H02 for the species containing one chemical acetylation).

S-20
2.2 Generate a *.tsv containing the name, sum formula and charge of your analyte of interest (Note: the sum formula has to contain all atoms including the acetyl group of the acetyllysine residues) (Table S5). 2.3 Run the python PICor package from the command line using the following command: "picor input.csv -o outfile.csv -x Rep -x Exp --molecules-file peptides.tsv -r --resolution 60000 -mz-calibration 200"¬ Where Rep and Exp are column names that do not contain intensity values and should be excluded. The input file containing intensity values is generated in section 2.1, the peptides.tsv file was generated in section 2.2. Figure S14 shows the heights of the different isotopologue species before (upper graph) and after (lower graph) PICor correction.

Raw data file conversion using msConvert
The single-acetylated H3(18-26) isotopologue species II*/II* and V*/VI*, as well as the doubleacetylated species VII*/VIII* are positional isomers and isobaric at MS1 level but are distinguishable based on their fragment ion spectra. To detect the isotopologue-specific fragment ions by PASTAQ whereas SEQUENCE_precursorSpecificIon.csv includes all ions capable for identification and quantification of a specific species. 4.4 The python script generates a list with all theoretically possible fragment ions with the following outputs: theoretical m/z (z= 2) of the intact isotopologue species(precursor m/z); type of fragment ion (b-ions, y-ions); theoretical m/z of the fragment ions; amino acid sequence (one letter code) of the fragment ion including the type of acetyl group of the acetyllysine residues. 4.5 Export isotopologue species-specific fragment ion information as a .csv file.

MS2 fragment ion quantification using PASTAQ
To determine site-specific acetylation levels for MS1 isobaric isotopologue species, the relative abundances of site-specific acetyllysine containing fragment ions were quantified using Pipelines and Systems for Threshold Avoiding Quantification of LC-MS/MS data (PASTAQ) [1] . PASTAQ can be downloaded under https://pastaq.horvatovichlab.com/.
5.1 For MS2 fragment quantification analysis with PASTAQ, a sample table must be generated (*.csv).
The sample table (see Table S6) has to include the following information: sample name: Sample name identical to the LC-MS raw data input file (*.mzML).    Figure S15). 5.6 Start the analysis to extract the MS2 spectra within the precursor m/z and retention time tolerance window as given in the sample table (see 5.1) (Note: It is not necessary to specify a tolerance range for m/z, as the width of peaks in the m/z dimension only depends on the resolution and type of the MS analyzer). Each MS/MS spectrum is searched for the theoretical fragments of interest obtained from the fragment table (see 5.2) by extracting a slice of data centered around the theoretical m/z. Finally, the fragments are quantified by performing a weighted Gaussian curve fitting using nonlinear least squares. If the fitting is successful, the area under the curve is calculated using the trapezoidal rule and this number is saved, along with the raw total intensity, fitted height, fitted m/z, and fitted width of the peak (sigma m/z) (Table S8). To ensure the fitting is correct, a plot for each raw data slice and its corresponding Gaussian fitting model is generated ( Figure S16). 5.7 Export the results as *.csv.

Python: read in and calculate site-specific quantification
For the site-specific quantification of isobaric isotoplogue species, a relative quantification of the isobaric species has to be performed at the MS2-level. The relative abundance of the isobaric isotopologue species is multiplied with the corresponding MS1 abundance. For this, we use the scripts MS2_abundances_H3.py and Site-specific_intensities_H3.py, which can be downloaded at https://github.com/functional-proteomics-and-metabolomics/CoMetChem. 6.1 For relative quantification of the site-specific isotopologue species at MS2 level run the Python script MS2_abundances_H3.py using the MS2 data from section 5.6 as an input. Specify the PATH were the PASTAQ outputs are stored and the dictionaries REPLICATES/TIMEPOINTS. 6.2 The script Site-specific_intensities_H3 takes the Results.csv file from the script MS2_abundances_H3.py and the file containing the MS1 data (see 2. 3) as input and corrects the results for the MS1 abundance. The SPECIES dictionary is necessary and has to be updated for the analysis of different isotopologue species. 6.3 The *.csv output file is used for further analysis

Calculation of site-specific half-lives and turnovers
7.1 Read in the file with the site-specific intensities into a python interpreter of your choice. 7.2 Using the SciPy package, fit an exponential curve over the ratio of heavy to heavy + light isotopologues over time using the following equation: S-25 The exponential factor k (turnover rate) can then be used to calculate the half-life of a species as 1 2 ⁄ = − (0.5) Figure S17 shows the corresponding for the single-acetylated K18ac H3(18-26) species. 7.3 To calculate abundance-corrected turnovers, the turnover rate (k) has to be multiplied with the abundance of the species (A).

Calculation of site-specific reaction rates of acetylation and deacetylation
8.1 To calculate the site-specific acetylation and deacetylation rates, a linear fit was performed over the abundance of different species in the first 8 hours using the python-based ecosystem SciPy and the function linregress ( Figure S18). The slope dx/dt and the intercept (Ax, abundance at T=0) are needed for further calculations. 8.2 Acetylation and deacetylation rates can be determined for the reactions outlined in Figure 3 in the main manuscript using the following equations: In the case of deacetylation, where Ax is the abundance of the acetylated species x at T=0 and x = × In the case of acetylation, where Asubstrate is the abundance of the non-acetylated species at T=0. Note: The data processing pipeline is currently operating in a targeted way. The extracted ion chromatograms of the isotopes of the different isotopologue species must be generated in FreeStyle (section 1) or a similar program that allows to quantify MS1 isotopes and provide EICs of the isotopes of the isotopologues as *.csv. In addition, the sequence and modification of the isobaric isotopologues species must be specified in the Python script to calculate isotopologue-specific fragment ions. The data processing modules FreeStyle (section 1) and Fragmentation.py (section 4) can be easily adjusted to analyze any isotopologue species of interest. For a more global analysis approach, the Fragmentation.py module can be linked to the output of database search engines if they specify peptide sequences in the single-letter amino acid code including the specification of the modified acetyllysine as described in Section 4, and the FreeStyle module can be exchanged by a program that allows to quantify MS1 isotopes and provide EICs of the isotopes of the isotopologue species as *.csv.