Abstract
The histone variant H2A.Z is an essential and conserved regulator of eukaryotic gene transcription. However, the exact role of this histone in the transcriptional process remains perplexing. In vertebrates, H2A.Z has two hypervariants, H2A.Z.1 and H2A.Z.2, that have almost identical sequences except for three amino acid residues. Due to such similarity, functional specificity of these hypervariants in neurobiological processes, if any, remain largely unknown. In this study with dissociated rat cortical neurons, we asked if H2A.Z hypervariants have distinct functions in regulating basal and activity-induced gene transcription. Hypervariant-specific RNAi and microarray analyses revealed that H2A.Z.1 and H2A.Z.2 regulate basal expression of largely nonoverlapping gene sets, including genes that code for several synaptic proteins. In response to neuronal activity, rapid transcription of our model gene Arc is impaired by depletion of H2A.Z.2, but not H2A.Z.1. This impairment is partially rescued by codepletion of the H2A.Z chaperone, ANP32E. In contrast, under a different context (after 48 h of tetrodotoxin, TTX), rapid transcription of Arc is impaired by depletion of either hypervariant. Such context-dependent roles of H2A.Z hypervariants, as revealed by our multiplexed gene expression assays, are also evident with several other immediate early genes, where regulatory roles of these hypervariants vary from gene to gene under different conditions. Together, our data suggest that H2A.Z hypervariants have context-specific roles that complement each other to mediate activity-induced neuronal gene transcription.
Significance Statement
Epigenetic regulation of activity-induced gene transcription is pivotal in mediating neuronal responses that underlie several transcription-dependent downstream brain processes, yet it remains poorly understood. Thus, understanding roles of such epigenetic processes and their core components, such as variant histone H2A.Z, is necessary to comprehend brain development and function. In vertebrates, H2A.Z has two hypervariants, H2A.Z.1 and H2A.Z.2, which are encoded by different genes and differ by only three amino acids. Despite such similar sequences, we provide evidence suggesting that they regulate nonoverlapping gene cohorts in neurons and play context-specific, nonredundant roles in activity-induced transcription of immediate early genes. Together, our findings represent a substantial departure from the current H2A.Z biology by shifting the focus to independent contributions of H2A.Z hypervariants.
Introduction
Activity-induced gene transcription, an integral component of neural plasticity, is conducted at the chromatin level by enabling the transcription machinery to access and “read-out” DNA (Greer and Greenberg, 2008). Several activity-induced neuronal immediate early genes (IEGs) are accessed independent of signal by a transcriptionally engaged RNA polymerase II (Pol II) that pauses proximal to the promoter until an activity signal is received (Saha et al., 2011; Adelman and Lis, 2012). Activity-induced “turn-on” signals release the paused Pol II, which then efficiently transcribes the gene by running through spools of nucleosomes (productive elongation). Nucleosomes, the fundamental chromatin unit, regulate Pol II elongation by acting as either barriers or facilitators (Owen-Hughes and Gkikopoulos, 2012; Teves et al., 2014). Such activity-dependent nucleosomal role selection primarily depends on post-translational modifications of the histone tails and nucleosomal remodeling, including turnover of different histone variants (Michod et al., 2012).
In mammals, three of the four canonical histones, H2A, H2B, and H3, have diversified into several variant histones with distinct amino acid sequences that can influence transcriptionally relevant nucleosomal properties and dynamics (Weber and Henikoff, 2014). Recently, a few studies have unveiled important roles of histone variants H2A.Z and H3.3 in neuronal gene transcription underlying brain development and function. For example, H2A.Z, which is conserved from yeast to human and has been implicated in both transcriptional activation and repression (Soboleva et al., 2014), plays a role in memory consolidation and cerebellar development in the brain (Zovkic et al., 2014; Yang et al., 2016). Similarly, H3.3, which differs from canonical H3 histones by only four to five amino acids (Filipescu et al., 2013), regulates embryonic and adult neuronal gene expression patterns, thereby playing an essential role in plasticity and cognition (Maze et al., 2015). This example of H3.3 illustrates how very small dissimilarities in sequences of two closely related histone proteins may produce large differences in neuronal gene transcription (Maze et al., 2014).
Similar small differences also exist within the H2A.Z histone class. In vertebrates, H2A.Z has two very closely related paralogs that have been referred to as “hypervariants”: H2A.Z.1 (formerly H2A.Z) and H2A.Z.2 (formerly H2A.V; Talbert et al., 2012). Two nonallelic genes, H2afz and H2afv, which are driven by independent promoters on different chromosomes, encode these H2A.Z hypervariants (Dryhurst et al., 2009; Matsuda et al., 2010). H2A.Z.1 and H2A.Z.2 differ by only three amino acids, which are located far apart on the polypeptide (Fig. 1A) and are structurally similar to their respective counterpart in that there are no usable difference in their antigenicity (e.g., hydrophobic vs hydrophilic, ring structure vs linear, or large vs small in residue size). Thus, there are currently no antibodies available to distinguish between H2A.Z.1 and H2A.Z.2. Due to this limitation, our understanding of H2A.Z is most likely a portmanteau of both hypervariant functions. Such a composite picture of their function would not be a concern if these hypervariants were functionally redundant and homogeneous (Soboleva et al., 2014). However, a few structural and functional features suggest the contrary. Recently resolved crystal structures of H2A.Z-containing nucleosomes revealed subtle differences in the H2A.Z.1 and H2A.Z.2 L1 loop (Horikoshi et al., 2013), but it remains unknown whether such differences manifest any functional consequences. In H2A.Z.1−/− mice, which are embryonic lethal (Faast et al., 2001), H2A.Z.2 is apparently unable to compensate for the loss of its twin isoform, thereby suggesting a unique role for H2A.Z.1 during early embryonic development. Similar functional inequality of H2A.Z.1 and H2A.Z.2 has also been recently reported in melanoma cells (Vardabasso et al., 2015). Taken together, there remains an untested possibility that these two H2A.Z hypervariants may not be functionally identical, but instead have discrete nonoverlapping roles (functional specificity) in gene transcription. We tested this hypothesis in neurons and report here that H2A.Z.1 and H2A.Z.2 have independent and context-specific roles in basal and neuronal activity-induced gene transcription.
Validation of H2A.Z hypervariant-specific RNAi. A, H2A.Z.1 and H2A.Z.2 amino acid sequences are depicted where three different amino acid residues are indicated in red. The target sequences of shH2afz and shH2afv are noted in blue. B, left, Representative images of smFISH of H2afv and H2afz in E18 rat cortex. Right, Quantification of H2afz and H2afv molecules; particle count per 212 μm2. N = 3, *p < 0.05 (one-tailed unpaired t test). Scale bar, 10 μm. C, Cortical neurons in culture were infected with lentiviruses delivering either shH2afz or shH2afv and total RNA was collected after indicated number of days. Compared to control neurons, changes in H2afz and H2afv mRNA levels were determined by qPCR, normalized to Gapdh mRNA levels (internal control) and are depicted here as fold change. NS, not significant. *p < 0.05. N = 3–6. D, Timeline of infection and treatment for all assays performed in dissociated cortical neurons. Cells were infected with lentiviruses between 5 and 7 d after plating (gray line), and assays were performed between 4 and 7 d after infection (between DIV10 and DIV14; blue line). E, Total acid-extracted histone from neurons infected with lentiviruses delivering either shH2afz or shH2afv or both for 6–7 d, resolved by electrophoresis, and blotted for indicated histone. N = 3. Note the cumulative effect of depleting both hypervariants when probed with total H2A.Z antibody. F, IP of nucleosomal H2A.Z from nuclear extract of neurons infected with lentiviruses delivering either shH2afz or shH2afv or both for 6 d resolved by electrophoresis. N = 3. G, H, Neurons were coinfected with indicated constructs for 4–5 d and nuclear extracts were resolved by electrophoresis. IS, insensitive; shRNA target regions were swapped in these constructs. N = 3.
Materials and Methods
Plasmids and subcloning
Several web-based software programs were used to design eight (four each) hypervariant-specific shRNA for RNAi. The shRNA were inserted into pLKO.1-puro (designed by RNAi consortium or TRC; obtained from Addgene) following the protocol on the Addgene website. Self-inactivating HIV lentivirus particles were produced by transfecting 293T cells with the shRNA vector, envelope (pMD2.G; Addgene), and packaging plasmids (psPAX2; Addgene). For construction of exogenous tagged constructs, DNA sequences encoding HA-H2A.Z.1, HA.IS1, HA-ISB, MYC-H2A.Z.2, MYC-IS2, and MYC-ISB were generated each with 5′ EcoRI and 3′ NotI restriction sites, using Life Technologies GeneArt Gene Synthesis “strings.” These EcoRI-NotI DNA fragments were cloned into the System Biosciences lentiviral plasmid CD811A-1, using the corresponding restriction sites on the vector. Subsequently, DNA constructs were characterized by restriction mapping using EcoRI and NotI enzymes (New England Biolabs), and inserted elements were validated by sequencing. These constructs were packaged into lentiviruses as mentioned above.
Dissociated neuronal culture, RNAi, and cell treatment
Cultures of cortical neurons were prepared from embryonic day 18 (E18) Sprague Dawley rats of either sex. All animal procedures were performed in accordance with the National Institute of Environmental Health Sciences (NIEHS) and the University of California Merced animal care committee’s regulations [NIEHS Institutional Animal Care and Use Committee (IACUC) approval: ASP#01-21; and University of California Merced IACUC approval: ASP#13-0007 and ASP#16-0004]. Dissociated cortical neurons were plated in Neurobasal medium (Invitrogen) supplemented with 25 μM glutamate (Sigma-Aldrich), 0.5 mM L-glutamine (Sigma-Aldrich), either B27 (Invitrogen) or NS21, and maintained in a similar medium without the glutamate. NS21 was prepared in the laboratory as previously described (Chen et al., 2008). Neurons were routinely used between 10 and 14 d in vitro (DIV), a stage when synapses have formed and are effective. For infection with recombinant lentiviruses, the viral supernatant was diluted in neuronal media and cells were infected between DIV5 and DIV7 at a multiplicity of infection ranging from 2 to 5 for 4–7 d. To induce gene transcription under resting conditions using synaptic circuits, we cotreated neurons with 50 μM bicuculline (Bic; Sigma-Aldrich) and 75 μM 4-aminopyridine (4AP; Acros Organics; Papadia et al., 2005). To induce gene transcription extrasynaptically, we blocked neuronal activity with 1–2 μM tetrodotoxin (TTX; Calbiochem) and activated the MAP kinase (MAPK) pathway (via PKC) with 1 μM phorbol ester myristate (PMA; Sigma-Aldrich; Schultz et al., 1997). To induce gene transcription at the intranuclear level, neurons were treated with 100 μM 5, 6-dichloro-1-β-D-ribofuranosyl-benzimidazole (DRB; Sigma-Aldrich) alone, or in combination with 250–500 nM triptolide (Trip; Tocris) for 1 h followed by washout with preconditioned media. To induce genes after homeostatic changes, a previously described protocol using prolonged treatment of TTX for 48 h followed by its washout was used (Saha et al., 2011).
Electrophysiological recordings
Whole-cell patch-clamp recordings were performed on dissociated cortical neurons coexpressing GFP and shRNA by an experimenter blinded to the identity of the treatment. Miniature EPSCs (mEPSCs) were assessed for frequency and amplitude as detailed previously (Saha et al., 2011). For visualization of infected neurons using GFP coexpression, the shRNA component was subcloned from the pLKO.1 to the pLKO.3G backbone as detailed previously (Saha et al., 2011). To ensure adequate knockdown, PCR was performed in sister neuronal cultures with the same batch of viruses to confirm their efficacy.
RNA extraction and gene transcription quantification
Total RNA was isolated from dissociated neurons using the RNeasy Mini kit (QIAGEN) with in-column DNase (QIAGEN) digestion. cDNA was synthesized using MuLV reverse transcriptase (Promega), random primers (Promega), oligo dTs (Promega), and RNase inhibitors (Thermo Scientific). qPCR was performed using iTaq Universal Sybr Green Supermix (Bio-Rad) and the Bio-Rad CFX Connect real-time PCR Detection System. Pre-mRNA was estimated as described previously (Saha, et al., 2011). Rat PCR primers used in this study are listed in Table 1.
List of primers used in the study
Expression microarray and data analysis
Gene expression analysis was performed using RNA obtained from four biological replicates. The Agilent Whole Rat Genome 4x44 multiplex format oligo arrays (014879, Agilent Technologies) were used following the Agilent 1-color microarray-based gene expression analysis protocol. Using 500 ng of total RNA, Cy3-labeled cRNA was produced according to the manufacturer’s protocol. For each sample 1.65 μg of Cy3-labeled cRNAs were fragmented and hybridized for 17 h in a rotating hybridization oven. Slides were washed and then scanned with an Agilent Scanner. Data were obtained using the Agilent Feature Extraction software (v9.5), using the 1-color defaults for all parameters. The Agilent Feature Extraction Software performed error modeling, adjusting for additive and multiplicative noise. Then, “gProcessedSignal” values were extracted for each probe and sample and used as the baseline expression measurements and quantile normalization was performed to normalize the samples. Probe-based expression data were collapsed down to a single measurement per gene by taking the median of all probes associated with each gene. Next, the average expression was calculated across all samples for each gene and the bottom 10% of genes with the lowest average expression was removed. Differentially expressed genes (DEGs) were identified by first performing a two-sample Student’s t test for each gene followed by Benjamini Hochberg correction. Genes were defined as differentially expressed if they had p < 0.001 and a fold change of 1.25 or greater. For the heat map depicting expression of these DEGs, the expression values were normalized to have a mean of zero and a standard deviation of one, and hierarchical clustering of these genes was performed. GO terms that were enriched for these DEGs were identified using the online tool DAVID (http://david.ncifcrf.gov/). GEO accession number for these datasets is GSE96886.
Sample preparation for electrophoresis
Neurons were lysed in ice-cold 1× RIPA buffer (25 mM Tris, pH 7.5, 150 mM NaCl, 1% Na-deoxycholate, 0.1% SDS, and 1% NP-40) and supplemented with 1:100 protease inhibitor cocktail (Sigma-Aldrich). Lysed neurons were sheared by sonication, cell debris pelleted at 15000 rpm for 5 min at 4°C, and the clarified supernatant transferred to a prechilled 1.7 ml microcentrifuge tube. Total cell extracts were denatured at 95°C for 5 min, using either home-made 5× Laemmli buffer, BOLT 4× Sample buffer (Life Technologies) and BOLT 10× reducing agent (Life Technologies), 2×-, or 4× Laemmli sample buffer (both from Bio-Rad). To obtain synaptic fractions, cell pellets were treated with 1 ml ice-cold Syn-PER reagent (Thermo Scientific catalogue number 87793) and supplemented with 1:100 protease inhibitor cocktail (Sigma-Aldrich). Neurons were incubated on ice for 10 min, gently triturated 10 times, centrifuged at 1200 × g for 10 min at 4°C, and supernatant was transferred to ice-cold microcentrifuge tubes, then recentrifuged at 15,000 × g for 20 min at 4°C. The pellet was then resuspended in 20 µl 2× Laemmli sample buffer (Bio-Rad). Acid-extracted histone fractions were prepared using a protocol described previously (Shechter et al., 2007). Histone pellets were denatured using 2× Laemmli sample buffer (Bio-Rad).
Immunoprecipitations (IPs)
All IPs were performed overnight, at 4°C in 1× IP buffer (0.5% Triton X-100, 0.002 M EDTA, 0.02 M Tris pH 7.75, 0.15 M NaCl, and 10% glycerol), or in 1× IP buffer from the SimpleChIP Enzymatic Chromatin IP kit (Cell Signaling Technology catalogue number 9003). Proteins were isolated from overnight IPs by incubating with Pierce Protein A/G Magnetic Beads (catalogue number 88803) for 1 h at 4°C. Subsequently, the beads were washed three times with IP buffer, denatured, separated by electrophoresis, and immune blotted as noted below.
Western blotting and imaging
Denatured protein samples were resolved on 4–20% (Bio-Rad catalogue number 4568095), 4–15% (Bio-Rad catalogue number 456–1083), 8–16% (Bio-Rad catalogue number 456–1103) Mini PROTEAN gels, or on BOLT 4–12% Bis-Tris Plus Gels (Life Technologies catalogue number NW04122BOX), in Tris/Glycine/SDS (Bio-Rad catalogue number 1610772), or 1× MOPS buffer (Life Technologies catalogue number B0001), respectively. Resolved proteins were transferred onto LF PVDF membrane, using the Bio-Rad TBT RTA kit and protocol (catalogue number 1704272). PVDF membranes were incubated at 4°C overnight with appropriate primary antibodies in 1× TBS-T with 0.5% BSA. The next day, membranes were washed three times in 1× TBS-T, probed with appropriate Alexa Fluor secondary antibodies (Life Technologies) for 30–45 min at room temperature, washed three times with 1× TBS-T, and imaged using Bio-Rad ChemiDoc Imaging System.
Chromatin digestion
Micrococcal nuclease (MNase) digests of rat dissociated neurons were performed using the SimpleChIP Enzymatic Chromatin IP kit (Cell Signaling Technology catalogue number 9003). Cells were fixed using Thermo Scientific Pierce methanol-free formaldehyde (catalogue number PI28906). The manufacturer’s protocol was used with the exception that Thermo Scientific Pierce Micrococcal Nuclease (catalogue number PI88216) was used in place of CST MNase. Postdigest, samples were briefly sonicated at low settings to rupture the nucleus and then were either subjected to chromatin IP (ChIP), DNA extraction, or protein co-IP following the CST kit protocol.
ChIP
Neurons were fixed with formaldehyde for 5–10 min, lysed in the sonication buffer (20 mM Tris, pH 7.8, 2 mM EDTA, 0.5 mM EGTA, 0.5% SDS, and inhibitors for proteases and phosphatases), and then sonicated for sixteen cycles (bursts and intervals of 30 s each) with the Bioruptor (Diagenode) to obtain 200- to 1000-bp genomic DNA fragments. A fraction of the sonicated sample was then immunoprecipitated overnight at 4°C in the 1× IP buffer with 2- to 5-μg antibody. Antigen-antibody complexes were immunoprecipitated with Pierce Protein A/G Magnetic Beads, washed once with low salt buffer, three times with high salt buffer, once with LiCl buffer, and once with Tris-EDTA buffer. Samples were reverse cross-linked at 65°C overnight, and chromatin DNA was eluted using the QIAquick Nucleotide Removal kit (QIAGEN). Eluted chromatin was quantified by qPCR.
Immunocytochemistry and microscopy
Infected neurons were washed twice with 1× ice-cold PBS (Fisher Scientific). The cells were then incubated with 4% paraformaldehyde (Sigma-Aldrich) in 1× PBS for 15 min at room temperature and then washed twice with 1× PBS, permeabilized at room temperature for 20 min with 0.5% Triton X-100 (Fisher Scientific), washed twice with 1× PBS, and blocked for 30 min with 10% goat serum (Gibco) in 1× PBS. Cells were incubated at 4°C overnight in 3% goat serum in 1× PBS with primary antibodies at 1:500 dilution. The next day, the primary antibody solution was removed, and cells were washed three times with 0.05% Tween (Fisher) in 1× PBS (0.05% PBS-T), and incubated with appropriate secondary antibody in 3% goat serum in 1× PBS for 45 min, washed three times with 0.05% PBS-T, and cured overnight using ProLong Anti-Fade Gold with DAPI (Life Technologies). Images were captured with a Keyence BZ9000-E microscope at 20× magnification.
Antibodies
Sheep anti-H2A.Z (Millipore, 09-862; RRID:AB_1587118), mouse anti-MYC (ThermoFisher Scientific, R95025; RRID:AB_2556560), rabbit anti-H3 (Cell Signaling Technology, 2650S; RRID:AB_2115124), mouse anti-β-ACTIN (ThermoFisher Scientific, AM4302; RRID:AB_2536382), rabbit anti-ARC (Synaptic Systems, 156003; RRID:AB_887694), rabbit anti-p-p44/42 MAPK (T202/Y204; Cell Signaling Technology, 4370S; RRID:AB_2315112), mouse anti-p44/42 MAPK (Erk1/2; L34F12; Cell Signaling Technology, 4696S; RRID:AB_390780), rabbit anti-HA (Cell Signaling Technology, 3724S; RRID:AB_1549585), mouse anti-MYC (Cell Signaling Technology, 2276S; RRID:AB_331783), rabbit anti-H2A (Cell Signaling Technology, 12349S; RRID:AB_2687875), mouse anti-H2B (Cell Signaling Technology, 2934S; RRID:AB_2295301), mouse anti-H4 (Cell Signaling Technology, 2935S; RRID:AB_1147658), rabbit anti-H2A.Z (Cell Signaling Technology, 2718S; RRID:AB_10694716), mouse anti-PSD95 (Antibodies Incorporated, 75-028; RRID:AB_2292909), rabbit anti-SH3 and multiple ankyrin repeat domains 3 (SHANK3) (Santa Cruz Biotechnology, SC30193; RRID:AB_2301759), goat anti-SHANK3 N16 (Santa Cruz Biotechnology, SC23547; RRID:AB_2187727), mouse anti-GluR2 (anti-GluA2; Millipore, MAB397; RRID:AB_2113875), mouse anti-HOMER1 (Synaptic Systems, 160011; RRID:AB_2120992), mouse anti-Synapsin1 (Synaptic Systems, 106001; RRID:AB_887805), rabbit anti-RPB1 (Cell Signaling Technology, 14958S; RRID:AB_2687876), rabbit anti-acetyl-H2A.Z Lys 5,7,11 (Millipore, ABE1363; RRID:AB_2687877), and rabbit anti-pSer5-RPB1 (Active Motif, 61085; RRID:AB_2687451).
Single molecule fluorescent in situ hybridization (smFISH)
Rat brains were harvested at E18 and individually embedded in 22 × 22 mm cryomolds using optimal cutting temperature compound (OCT) and flash frozen in 2-methylbutane, which was cooled to −20°C on a dry ice/ethanol slurry. Brains were sectioned on a Leica cryostat at 20 μm and processed for single-molecule FISH according to the RNAscope Fluorescent Multiplex kit instructions (Advanced Cell Diagnostics). The following probes were used with the RNAscope fluorescent multiplex reagent kit: Rn-H2afv (20 ZZ probes targeting 414–1815 of NM_001106019.1; part# 495451), Rn-Hz1fz-C2 (13 ZZ probes targeting 12–792 of NM_022674.1; part# 495441-C2), 3-plex Positive Control Probe, rat (part# 320891), and 3-plex Negative Control Probe (part# 320871). Images were acquired on a Zeiss 880 confocal microscope using a 40× oil-immersion lens. The image acquisition parameters were set using 3-plex negative controls (cDNA probes against bacterial RNAs not present in rat tissue) in each of the three channels (fluorescein isothiocyanate, CY3 and CY5) so that any signal above the level of background was acquired. At least three 212 × 212 μm z-stack images (1-μm steps) of somatosensory cortex were acquired per rat. Images were processed and analyzed using FIJI software (NIH v.2). Maximum intensity projection images were thresholded to an average value per channel and particles were counted using the analyze particle function in FIJI. Background particle counts were subtracted from each value before taking the average particle count per channel per rat (N = 3).
Nanostring assay and data analysis
NanoString probes were designed for indicated pre-mRNAs by NanoString technologies and assays were performed following the manufacturer’s protocol. Data were processed in R (version 3.3.2), using the limma package (version 3.30.8) from Bioconductor (version 3.4). Data were log2 transformed, then batch-adjusted within each experiment and time point to account for sample processing errors. The pre-mRNA measurements were normalized using pre-mRNA housekeeper genes. Sample quality and processing were reviewed and confirmed using MA-plots and heat maps at each step. Pairwise comparisons were performed using the limma-voom methodology as previously described. Statistical hits were defined by requiring one sample group mean log2 signal at least 2, absolute fold change above 1.5, and Benjamini Hochberg adjusted p < 0.01. Heat maps were generated using group mean log2 normalized abundances, centered to the indicated control group mean. Violin plots were generated using the same values, centered to the treated scrambled for the respective group of treatments. Gene class changes per experiment and time point were tested using group mean values in a nonparametric, paired Wilcoxon test, with gene as the pairing factor, to test for coordinated gene changes between treated shH2A.Z.1 or treated shH2A.Z.2 and the corresponding treated scrambled of the same experiment time point.
Experimental replicates and statistics
Biological and technical replicates for each figure are indicated as N and n, respectively, in corresponding figure legends. Error bars represent standard error of mean throughout this article. Statistical comparison of datasets was performed with the two-tailed Student’s t test (with Bonferroni corrections for multiple comparisons), two-way ANOVA, or Wilcoxon paired nonparametric test, unless otherwise stated in the figure legend.
Results
H2A.Z hypervariants and strategy to target them independently
First, we asked if both H2A.Z isoforms, H2A.Z.1 and H2A.Z.2 (Fig. 1A), are expressed in the brain. Using the smFISH approach and specific probes that bind to H2afz or H2afv, we detected both H2A.Z isoforms in the E18 rat brain (Fig. 1B). While there were significantly more H2afz particles compared to H2afv, we observed that the majority of DAPI-positive cells expressed both RNAs. Further studies were performed in dissociated cells obtained from E18 rat brains.
Next, we used RNAi to target these hypervariants independently. Both genes have been duplicated during evolution and have multiple mammalian “gene copies” that have variable noncoding regions despite having the same coding region. For example, in rat, when compared to the annotated H2afz on chromosome 2, an unnamed “pseudogene” on chromosome 16 (rat rn6 location: 71438461-71439260) codes for an mRNA that retains the H2afz coding sequence in totality amid dissimilar noncoding regions. Therefore, to ensure knockdown of these probable pseudogene mRNAs as well, the coding region was targeted for RNAi. Initially, we designed four shRNA against each of the H2afz and H2afv mRNA coding regions and evaluated their knockdown efficacy by qPCR. Only one of these four shRNA pairs efficiently knocked down each hypervariant independently; others knocked down both hypervariants to some degree. Interestingly, these efficient shRNAs against H2afz (shH2afz) and H2afv (shH2afv) targeted the same region on the respective mRNA (Fig. 1A). These 24-bp target regions on H2afz and H2afv mRNAs differ by only 5 bp. Because of such scant difference, we tested whether shH2afz and shH2afv were indeed specific for their respective hypervariant targets.
Lentivirus-mediated delivery of shH2afz to dissociated rat cortical neurons significantly knocked down H2afz, but not H2afv, mRNA at all tested times after viral infection (Fig. 1C). Similarly, H2afv, but not H2afz, mRNA was knocked down by shH2afv (Fig. 1C). We also noted that loss of one hypervariant was not compensated by significantly enhanced expression of the other hypervariant. Because efficient knockdown was noted between days 4 and 7 after infection, we decided to use this timeline for the rest of our assays (Fig. 1D). At the protein level, as detected by a commercial anti-H2A.Z antibody that does not discriminate between hypervariants, infection with shH2afz or shH2afv viruses partially depleted endogenous H2A.Z levels without affecting the expression of other histones (Fig. 1E). Although we were unable to completely abolish H2A.Z signal on infection with both shH2afz and shH2afv viruses, H2A.Z protein depletion was more pronounced on coinfection. We suspect that the protein remaining after coinfection is nucleosomal H2A.Z protected in the heterochromatin (Fig. 1E,F). The cumulative effect of such coinfection suggests that neurons express both hypervariants, and that they are sensitive to their respective shRNA only.
To test the above stated possibility further, we coexpressed shH2afz or shH2afv with either tagged HA-H2A.Z.1 or tagged Myc-H2A.Z.2. Protein expression of HA-H2A.Z.1 was attenuated by shH2afz, but not shH2afv. Similarly, MYC-H2A.Z.2 protein expression was sensitive to shH2afv, but not shH2afz (Fig. 1G,H, left panels). To test shH2afz and shH2afv specificity even further, we exchanged their 24-bp target sequence in our exogenous constructs with the expectation that shRNA sensitivity of these constructs will also be reversed. As expected, HA-H2A.Z.1-IS, the H2afz mRNA with the shH2afv 24-bp target sequence, is sensitive to shH2afv, but not shH2afz. Similarly, Myc-H2A.Z.2-IS is sensitive to shH2afz, but not shH2afv (Fig. 1G,H, right panels). Taken together, we show that shH2afz and shH2afv are specific for H2afz and H2afv mRNAs respectively, and that we can efficiently knockdown each targeted hypervariant independently.
H2A.Z hypervariant depletion alters expression of nonoverlapping sets of genes
To examine the impact of individual H2A.Z hypervariant knockdown on basal level of neuronal gene transcription, we performed a genome-wide microarray analysis of RNA obtained from dissociated cortical neurons depleted of H2A.Z.1 or H2A.Z.2. Gene expression levels in these samples were compared to samples from neurons treated with a scrambled shRNA (referred to as control neurons hereafter). The degree of change in the mRNA levels for every gene, whose transcription was altered due to H2A.Z.1 or H2A.Z.2 deficiency (>1.25-fold), is plotted in Figure 2 as a heat map. A total of 338 genes were differentially expressed due to depletion of H2A.Z.1 (154 genes were upregulated and 184 were downregulated), while H2A.Z.2 depletion altered expression of 271 genes (94 genes were upregulated and 177 were downregulated). Only 24 genes were sensitive to depletion of either H2A.Z isoform (Fig. 2B,C; Tables 2–5). Such relative dearth of common genes affected by depletion of either hypervariant suggests that H2A.Z.1 and H2A.Z.2 regulate, directly or indirectly, independent gene cohorts. Functional gene ontology analysis of H2A.Z.1 and H2A.Z.2 target genes revealed that targets of the former are enriched for ionotropic postsynaptic receptors (Fig. 2D). Three of the top 10 molecular GO pathways belonged to this category. This finding contrasts with H2A.Z.2 target genes, where none of the hits were related to synaptic function. Interestingly, knockdown of H2A.Z.1, but not H2A.Z.2, led to altered expression of several candidate genes linked to psychiatric disorders such as autism and schizophrenia (Neurexin1, Shank3, CaMK2G, Vamp7, Gsk3b, Tuba8, Grm8, etc.).
H2A.Z hypervariant depletion alters baseline expression of mostly nonidentical gene cohorts. A, Heat map showing changes in gene expression as detected by microarray analysis after depletion of H2A.Z.1 and H2A.Z.2 compared to control neurons. N = 4. Genes showing >1.25-fold change with p < 0.001 are shown here. B, C, Venn diagrams depicting total number of genes that are either upregulated (B) or downregulated (C) after depletion of H2A.Z.1 or H2A.Z.2. Note the nonoverlapping nature of such gene sets. D, Top five GO terms derived from molecular function gene ontology analysis (DAVID) for H2A.Z.1 and H2A.Z.2 target genes.
List of DEGs after H2A.Z.1 and H2A.Z.2 KD-part A (upregulated after H2A.Z.1 KD)
H2A.Z hypervariant depletion changes synaptic protein expression
Because glutamatergic synaptic function was indicated as a hypervariant-specific GO term in the above analysis, and to test whether H2A.Z hypervariant depletions have any impact on the physiologic health of neurons, we performed electrophysiological analyses of membrane and excitatory synaptic properties. Resting membrane potentials and input resistances were not significantly different between control and H2A.Z.1- or H2A.Z.2-depleted neurons. Presynaptic function and/or synaptic density, as inferred by mEPSC frequency, was also not significantly different between conditions (mEPSC frequency per 3 s for control, H2A.Z.1-, and H2A.Z.2-depleted neurons were 0.31 ± 0.31, 0.14 ± 0.06, and 0.13 ± 0.07, respectively). However, mEPSC amplitude, typically indicative of postsynaptic function, was significantly diminished in neurons depleted of either H2A.Z hypervariant (Fig. 3A,B). To better understand mechanisms of such mEPSC amplitude reduction in H2A.Z hypervariant-depleted neurons, we assessed mRNA and protein levels of several postsynaptic molecules, including many from our microarray screen. A few interesting cases representing diverse effects of hypervariant depletion are presented in Figure 3C,D.
H2A.Z hypervariant-specific RNAi has synaptic effects. A, B, Amplitude of mEPSCs decrease after H2A.Z.1 or H2A.Z.2 depletion. Shown are example traces in A with quantification of mEPSC amplitude is in B (N = 3; n = 15 neurons for Sc, 10 for H2A.Z.1 KD, and 8 for H2A.Z.2 KD; *p < 0.05 and **p < 0.001). Neurons were infected for 5–7 d. C, mRNA levels of indicated genes, normalized to Gapdh, represented as fold change of their abundance in control neurons (broken line at 1 on y-axis). N = 3–5. *p < 0.05. D, Representative Western blottings of synaptic fraction from control or H2A.Z hypervariant-depleted neurons as indicated. Loading control (β-ACTIN) for individual blots are provided and molecular weights of the nearest protein ladder band is indicated. Arrow indicates quantified SHANK3 band (see text). Neurons were infected for 7 d. N = 3. Note: levels of synapsin, a presynaptic molecule, remain unaltered after hypervariant depletions.
A contrasting effect of hypervariant depletion was noted in two instances: postsynaptic density protein 95 (PSD-95) and SHANK3. PSD-95 is the anchoring protein encoded by Dlg4 and found exclusively in the postsynaptic density, while SHANK3 is another member of the postsynaptic density that is encoded by the high-confidence autism candidate gene Shank3. These two gene products were strongly reduced in neurons depleted of H2A.Z.1, but not H2A.Z.2 (Fig. 3D). Compared to control, PSD-95 band intensities were 48.58 ± 12.4%, 102.34 ± 26.72%, and 65.95 ± 17.42% in H2A.Z.1, H2A.Z.2, and both KD, respectively (N = 3). Such an effect of H2A.Z.1 KD on PSD-95 must be indirect (not regulated at the gene transcription level), as Dlg4 mRNA levels remained unchanged in H2A.Z.1-depleted neurons (Fig. 3C). For SHANK3, compared to control, band intensities were 20.25 ± 7.6%, 47.74 ± 20.35%, and 15.34 ± 12.6% in H2A.Z.1, H2A.Z.2, and both KD, respectively (N = 3). However, the effect on SHANK3 is interesting, as Shank3 mRNA levels are upregulated specifically in response to H2A.Z.1 depletion, perhaps as a response to cellular feedback pathways triggered by loss of the protein.
We also tested levels of the ligand-gated ion channel GluA2 (AMPA receptor; coded by Gria2, also called GluR2) and postsynaptic density protein HOMER1 (coded by Homer1) after H2A.Z hypervariant depletion. As shown in Figure 3D, we observe no notable difference in signal of either GluA2 or HOMER1B/C after H2A.Z.1 or H2A.Z.2 knockdown. However, the activity-induced short Homer1a mRNA isoform was downregulated by H2A.Z.1 loss and upregulated by H2A.Z.2 depletion (Fig. 3C). Further studies are required to confirm if such effects are direct influences of H2A.Z hypervariants on Homer1a transcription.
Together, we find that depletion of both H2A.Z.1 and H2A.Z.2 leads to decreased amplitude of mEPSCs, which may be explained by altered expression of postsynaptic proteins. The effects of H2A.Z.1 and H2A.Z.2 knockdown on the tested synaptic proteins, however, are apparently mediated by transcription-independent means in most cases given the lack of change at the pre-mRNA level.
H2A.Z.1 and H2A.Z.2 depletion differentially impact Arc/Arg3.1 transcription
To study direct roles of H2A.Z hypervariants on transcriptional regulation of a model gene, we chose to study the IEG Arc/Arg3.1 (referred to as Arc from here on), which responds rapidly to neuronal activity (Guzowski et al., 1999; Saha et al., 2011). Arc was chosen for several reasons: (1) it is a neuron-specific IEG, which has been extensively studied previously by us and others, (2) its promoter and transcriptional start site (TSS) regions are enriched with H2A.Z-containing nucleosomes (our unpublished observations), and (3) it is one of the genes recently shown to be regulated by H2A.Z in cognitive functions (Zovkic et al., 2014). Neuronal activity was stimulated in dissociated neurons by GABAA receptor antagonist (Bic)-induced disinhibition and K+ channel blocker (4AP)-induced increased burst frequency (Bic + 4AP; Papadia et al., 2005). Using the Bic + 4AP protocol, we induced gene transcription and evaluated Arc mRNA abundance at 15 min after treatment. Basal and induced Arc mRNA levels were unaffected by H2A.Z.1 depletion but were significantly diminished in H2A.Z.2-depleted cells (Fig. 4A).
H2A.Z hypervariant-specific RNAi has differential effects on activity-induced Arc transcription. A, Graphical representation of Arc mRNA level at 15 min after treatment with Bic + 4AP in indicated neuronal groups as detected by qPCR and normalized to Gapdh. N = 3–4 for each group; *p < 0.05. NS, not significant. B, C, Representation of Arc pre-mRNA level at different time points after treatment with Bic + 4AP (B) or TTX + PMA (C) in indicated groups of neurons (5–6 d postinfection) as detected by qPCR and normalized to Gapdh. N = 3–4 for each time point; *p < 0.05 and **p < 0.01. NS, not significant. Two sets of data in statistical comparison are represented by the position and the color of asterisk(s) or a NS. For example, in B, the green * on top of the 15 min H2A.Z.2 KD time point (red triangle) represents a statistical difference between the means of H2A.Z.2 and H2A.Z.1 KD group (green square) at the same time point. D, E, Representative Western blottings of whole-cell lysate from neurons infected as indicated for 5–6 d. Neurons were treated with Bic + 4AP (D) or TTX + PMA (E) for either 60 min to detect ARC, or for 15 min to detect phosphorylated ERK1/2. F, Representation of Arc pre-mRNA level at 15 min after treatment with Bic + 4AP in indicated neuronal groups as detected by qPCR and normalized to Gapdh. N = 3; **p < 0.01.
Because mRNA abundance is not necessarily a reflection of transcription per se but rather the net result of both synthesis and decay rates, we next evaluated the direct Arc transcription product: its pre-mRNA (Wu and Brewer, 2012). We detected no significant differences in the basal level of Arc pre-mRNA with H2A.Z hypervariant depletions (Fig. 4B, 0-min time point). After Bic + 4AP treatment, Arc pre-mRNA levels increased significantly within 15 min and continued to increase with time in control neurons. However, in agreement with the mRNA data (Fig. 4A), the pre-mRNA response in H2A.Z.2-depleted neurons was attenuated at the early time point (15 min) and was also significantly reduced at later points (Fig. 4B). In contrast, H2A.Z.1 depletion had no effect on Arc pre-mRNA levels at 15 min. However, at later time points, H2A.Z.1-depleted neurons had significantly less Arc pre-mRNA than in control neurons and yet had significantly more than in neurons depleted of H2A.Z.2 (Fig. 4B). These results suggest that although the loss of H2A.Z.1 impaired Arc transcription significantly at later time points, loss of H2A.Z.2 had a much stronger effect at all time points. At the protein level, hypervariant-specific effects on activity-induced ARC abundance were generally reflective of Arc pre-mRNA expression patterns (Fig. 4D). When both hypervariants were depleted, pre-mRNA results were not significantly greater than H2A.Z.2 depletion alone, and were in fact intermediate between the two hypervariants. While trying to knockdown both H2A.Z.1 and H2A.Z.2, we had to titer our viral infection and infect with less of each virus kind to avoid causing viral overload-induced cell death. Therefore, less knockdown was achieved in the dual KD cells compared to their single KD counterparts. Such reduced knockdown may explain the intermediate effect of knocking down both H2A.Z.1 and H2A.Z.2 on Arc transcription.
Because some synaptic alterations were noted due to H2A.Z hypervariant depletion and because Bic + 4AP treatment is a synapse-dependent activity induction protocol, we next tried to induce Arc nonsynaptically by jump-starting the MAPK pathway, otherwise employed by activity-induced gene transcription signals (Rosen et al., 1994; Wiegert and Bading, 2011). The MAPK pathway was extrasynaptically activated via PKC using PMA (Schultz et al., 1997) in the presence of TTX to prevent any extracellular input via neuronal activity, and this treatment induced rapid transcription of Arc (Fig. 4C). When H2A.Z hypervariant-depleted neurons were stimulated using this protocol, the outcome was similar to what was seen with the Bic + 4AP protocol, both at the pre-mRNA and protein levels (Fig. 4C,E). These observations suggest that consequences of H2A.Z hypervariant depletion on activity-induced Arc transcription are most likely not limited to effects of changing synaptic components. Supporting this assertion is the observation that H2A.Z hypervariant depletion had no effect on activity-induced activation of the MAPK pathway as indicated by phospho-ERK1/2 levels (Fig. 4D,E).
To further test the specificity of the H2A.Z.2 depletion-induced phenotype, we attempted to rescue rapid Arc induction by combining shRNA-mediated knockdown with coexpression of shRNA-insensitive H2A.Z.2 (Myc-H2A.Z.2-IS). Several such experimental designs were unsuccessful, likely due to lack of regulatory post-translational modifications on exogenous constructs in the absence of endogenous proteins and/or perturbation of a critical balance of endogenous levels of these hypervariants and their chaperones. The latter possibility is based on a recent study where it was demonstrated that phenotype loss due to depletion of ANP32E, an H2A.Z chaperone that removes nucleosomal H2A.Z, could be rescued by codepletion of ANP32E and H2A.Z (Alatwi and Downs, 2015). Accordingly, we codepleted ANP32E (92.0 ± 0.1% knockdown on day 5 after infection) with H2A.Z.2 and observed a partial rescue of rapid Arc transcription (Fig. 4F). Taken together, our data show nonredundant roles of H2A.Z hypervariants in the regulation of Arc transcription and hint at the importance of a fine balance between histones and their chaperones in mediating such responses.
H2A.Z hypervariants are present near the Arc TSS and regulate its transcription
H2A.Z is enriched near Arc TSS (Zovkic et al., 2014). To identify specific hypervariant(s) in this enrichment, we exogenously expressed distinctly tagged hypervariants (HA-H2A.Z.1 and MYC-H2A.Z.2) in neurons and probed for their expression using antibodies against the HA- or MYC-epitope tags (Fig. 5A). Usually, protein overexpression is concerning, because it may cause, among other things, artifacts that include ectopic subcellular localizations and inaccurate protein complex formations. Despite such limitations, we reasoned that increased dosage of tagged-H2A.Z hypervariants in the free-floating nucleoplasmic pool may not have any negative consequences as H2A.Z incorporation into its functionally relevant complex (nucleosomes) is tightly regulated by chaperone-dependent mechanisms (Billon and Cote, 2013), and so this overexpression will not necessarily have unwanted nucleosomal incorporation of H2A.Z despite its artificial abundance.
Distinct properties of HA-H2A.Z.1- and MYC-H2A.Z.2-containing nucleosomes. A, Representative image of immunocytochemical analysis of HA-H2A.Z.1 and MYC-H2A.Z.2 expression in unstimulated neurons. Approximately 70% neurons expressed both hypervariants. Scale bar, 20 μm. B, Co-IP assay from nuclear lysate preparation of neurons exogenously expressing a control vector or HA-H2A.Z.1 and MYC-H2A.Z.2. N = 3. C, Representative image of genomic DNA obtained from untreated or MNase-treated nuclei of neurons treated as described above. The band size of ∼150 bp indicates that MNase digestion predominantly produced mononucleosomes. N = 3. D, Co-IP assay from single nucleosome preparation of neurons exogenously expressing HA-H2A.Z.1 and MYC-H2A.Z.2. Antigen-antibody complex was resolved by electrophoresis and blotted for indicated histones. The 15-kDa protein ladder band is indicated in all panels. Arrow indicates unique band found only in anti-MYC precipitate. N = 3.
To verify nucleosomal incorporation of tagged-H2A.Z hypervariants, we performed co-IP assays. These assays, with anti-IgG, anti-H2A.Z, anti-HA and anti-MYC from nuclear lysates, or anti-HA and anti-MYC from mono-nucleosome preparations (Fig. 5B–D), detected stable association of both hypervariants with other nucleosomal histones H2B, H2A, H3 and H4 suggesting that tagged-H2A.Z molecules are incorporated into nucleosomes. Interestingly, (1) HA and MYC IPs probed for their respective hypervariants revealed very little association between HA-H2A.Z.1 and MYC-H2A.Z.2, suggesting that very few H2A.Z nucleosomes contain both hypervariants; and (2) we detected a unique band in H2B, presumably representing an unidentified post-translationally modified state of H2B, in the precipitate with MYC-H2A.Z.2, but not HA-H2A.Z.1 (Fig. 5B,D). These data suggest that certain H2A.Z.2-containing nucleosomes may have unique post-translational signatures and possibly unique functions.
To confirm that functional incorporation of tagged-H2A.Z molecules into the Arc promoter region did not interfere with transcriptional activity, we expressed HA-H2A.Z.1 and MYC-H2A.Z.2 in mature neurons and induced Arc transcription with PMA + TTX. Overexpression of one or both hypervariants did not alter Arc transcription at pre-mRNA or protein level (Fig. 6A).
Both H2A.Z hypervariants are present near the Arc TSS and mediate its transcription. A, Representative Western blotting of whole-cell lysate from neurons infected as indicated for 5–6 d. Neurons were treated with TTX + PMA for 60 min to assess ARC expression. N = 3. B, top, Graphical map (not to scale) depicting relative position of ChIP primers near Arc TSS. Bottom, Quantification of H2A.Z, MYC-H2A.Z.2, and HA-H2A.Z.1 binding to the Arc TSS region in untreated neurons as determined by ChIP with antibodies against H2A.Z, MYC, and HA epitopes. N = 4. C, Quantification of H2A.Z, Ac-H2A.Z, MYC-H2A.Z.2, and HA-H2A.Z.1 binding to Arc TSS region in neurons treated with Bic + 4AP for indicated time periods. N = 4. Data are normalized as indicated and expressed as fold change compared to untreated samples (denoted by the broken line at 1 on y-axis).
Next, we performed ChIP assays to detect H2A.Z and its hypervariants near the Arc TSS. We tested several anti-H2A.Z antibodies and found the ChIP grade antibody from Millipore to be most efficient. Using this antibody, we detected significant enrichment of H2A.Z near the Arc TSS (primer pair 1, 2, and 3) compared to Arc gene body (primer 4; Fig. 6B). Primer pairs 1 and 3 roughly represent the −1 and +2 (often referred to as +1) nucleosomes across the Arc TSS (Kim et al., 2010), while primer pair 2 denotes a region across the TSS, which instead of being nucleosome-free, houses a nucleosome (+1) in Arc (see H3K4.me data in the study by Kim et al., 2010). Like total H2A.Z, both MYC-H2A.Z.2 and HA-H2A.Z.1 were also enriched near the TSS compared to the gene body (Fig. 6B). Although strong signals were obtained near the TSS with both anti-HA and anti-MYC ChIP compared to empty vector controls, these signals cannot be directly compared as their apparent difference is likely to be due to dissimilarity in antibody efficiency. For example, this trend is reversed when similar ChIP assays were performed after overexpressing hypervariants with swapped tags (data not shown). Taken together, we conclude that both H2A.Z hypervariants are present near the Arc TSS.
To distinguish the potentially different roles of H2A.Z hypervariants in activity-induced Arc transcription at the gene level, we treated neurons with Bic + 4AP for various time periods and performed ChIP with anti-H2A.Z, anti-MYC, anti-HA, and anti-acetylated H2A.Z (Ac-H2A.Z) antibodies (Fig. 6C). Acetylation of H2A.Z is an important epigenetic modification associated with active or primed promoters (Soboleva et al., 2014). On Bic + 4AP treatment, H3-normalized total H2A.Z levels were unchanged at any of the three nucleosomes across the Arc TSS at all time points tested (see Discussion). However, ChIP signal for Ac-H2A.Z drops significantly at all three nucleosomes at 15 min after treatment and remains reduced in the +1 and +2 nucleosomes after 30 min. Although total level of H2A.Z remained unchanged, we detected a significant drop of H2A.Z-normalized MYC-H2A.Z.2 and HA-H2A.Z.1 at the +2 nucleosome 15 and 30 min after treatment, suggesting ongoing activity-induced nucleosomal turnover of H2A.Z (nucleosomal remodeling). Interestingly, at the +1 nucleosome, only MYC-H2A.Z.2 levels dropped significantly at 15 min after treatment, indicating a possible H2A.Z.2-specific role for this nucleosome at this early time point. This inference is corroborated by our data in that the most severe effect on Arc induction was found after the depletion of H2A.Z.2 (Fig. 4B–E).
H2A.Z.2 is required for RNA Pol II pausing near the Arc TSS
To study the underlying mechanism of H2A.Z.2’s role in rapid Arc transcription (Fig. 4), we tested the effect of its knockdown on histone acetylation and promoter proximal RNA Pol II pausing based on recent reports (Saha et al., 2011; Vardabasso et al., 2015). In melanoma cells, knockdown of H2A.Z.2, but not H2A.Z.1, attenuates acetylation of histone H3 and H4 globally (Vardabasso et al., 2015). We observed no such effect of H2A.Z.2 depletion on H3 or H4 acetylation in neurons (data not shown). This effect is presumably a cell-type-specific phenomenon relevant to dividing cells. Instead, we noted that H2A.Z.2 depletion led to a decrease in Ac-H2A.Z levels globally (Fig. 7A) and near the Arc TSS (Fig. 7B,C). Note that the data presented in Figure 7C are normalized by total H2A.Z, the level of which remained unaltered after hypervariant knockdown. The effect was not mimicked by H2A.Z.1 depletion. This suggests that H2A.Z.2 depletion results in selective loss of acetylation, not total H2A.Z, near the TSS.
H2A.Z.2 facilitates priming of the Arc promoter. A, Representative Western blotting of whole-cell lysate showing Ac-H2A.Z levels in neurons depleted of H2A.Z hypervariants as indicated. N = 3. B, Graphical map (not to scale) depicting relative position of ChIP primers near Arc TSS. C, Quantification of Ac-H2A.Z (normalized by total H2A.Z) binding to Arc TSS region in neurons treated as indicated, determined by ChIP. N = 4. Two sets of data in statistical comparison are represented by the position and the color of asterisk(s). D, ChIP data demonstrating binding of RPB1 to the Arc TSS region in control neurons or H2A.Z hypervariant-depleted neurons as indicated. N = 4. Note the enrichment of RPB1 at the TSS (primer 2) in control neurons (black) indicating promoter-proximal RNA Pol II pausing. E, ChIP data demonstrating binding of pSer5-RPB1 to the Arc TSS in control and hypervariant-depleted neurons. N = 3. *p < 0.05 and **p < 0.01.
Presence of TSS proximal H2A.Z highly correlates with another chief component of rapid IEG transcription, RNA Pol II pausing (Day et al., 2016). We have previously shown that RNA Pol II is loaded and paused near TSS of Arc and other rapidly induced IEGs, and underlies the ability of these genes to respond to activity within a few minutes (Saha et al., 2011). Here, we tested whether there is a relation between Pol II pausing and H2A.Z hypervariants at the Arc TSS. Using ChIP with an antibody against RPB1 (N-terminus), the largest Pol II subunit, we found that RNA Pol II pausing at the Arc TSS is diminished significantly in neurons depleted of H2A.Z.2, but not H2A.Z.1. While it is unlikely that our anti-RPB1 antibody precipitates nonspecific proteins or is less efficient after hypervariant depletion, we sought to corroborate our findings with another antibody to further ensure that H2A.Z.2 depletion was affecting paused Pol II. We performed ChIP with an antibody specific to RPB1 phosphorylated at serine 5 (pSer5) residue of C-terminal domain heptad repeats. pSer5-RPB1 represents the pool of RPB1 that is transcriptionally engaged, but has not entered productive elongation (paused). Our data indicate that actively engaged Pol II at the +1 nucleosome of Arc is significantly reduced by the loss of H2A.Z.2, but not H2A.Z.1 (Fig. 7E). This hypervariant-specific loss of Pol II pausing provides an explanation for our observations at 15-min post-Bic + 4AP treatment (Fig. 4B,C), where Arc pre-mRNA levels were minimal in H2A.Z.2-, but not in H2A.Z.1-, depleted neurons.
H2A.Z.1 plays context-dependent roles in Arc transcription
To determine a potential role for H2A.Z.1 in Arc transcription, we hypothesized that it may facilitate Pol II recruitment or elongation during active transcription. To test these possibilities, we induced Arc transcription using the DRB washout protocol with or without Trip. DRB is a nucleoside analog that inhibits mammalian transcription by interfering with the DRB sensitivity-inducing factor complex that enables RNA Pol II pausing in the presence of negative elongation factor (NELF), and then facilitates productive elongation once NELF falls off in response to external signals (Kwak and Lis, 2013). Because DRB can be washed out, DRB washout assays have been used to achieve rapid and synchronized transcription elongation without additional extranuclear signals (Lee et al., 2015). When control neurons were subjected to this assay, we observed a rapid and robust Arc induction (Fig. 8A). DRB washout resulted in significant upregulation of Arc pre-mRNA levels within three minutes, which then continued to increase with time to reach >100-fold at 15 min. Depletion of H2A.Z.2 resulted in, as we saw with Bic + 4AP and TTX + PMA assays, an attenuated response at the early time point that then partially recovered at later time points while remaining significantly less than the normal response (Fig. 8A). In contrast, depletion of H2A.Z.1 had no effect on DRB washout-triggered robust Arc transcription (Fig. 8A), indicating that H2A.Z.1 is likely not involved in the productive elongation phase.
H2A.Z.1 plays context-dependent roles in Arc transcription. A, Graphical representation of the level of Arc pre-mRNA at different time points after 1-h treatment with DRB followed by its washout for indicated time periods in the presence or absence of Trip in indicated groups of neurons (infected for 5–6 d) as detected by qPCR and normalized to Rn18s (18s rRNA; DRB inhibition-insensitive, Pol III-dependent transcription). Gray, green and red asterisk(s) or NS represent statistical comparison between control (Scrambled) responses and responses in Scrambled w/Trip, H2A.Z.1 KD w/Trip and H2A.Z.2 KD, respectively. Connecting brackets delineate other comparisons, shown in black asterisk(s) or NS. B, Similar dataset as in A, except that control neurons or hypervariant-depleted neurons were treated with TTX for 48 h followed by its washout for indicated time periods. Two sets of data in statistical comparison are represented by the position and the color of asterisk(s) or NS. For both A, B, N = 4–5 for each time point; *p < 0.05 and **p < 0.01. NS, not significant.
Next, to test for any role of H2A.Z in RNA Pol II recruitment, a phenomenon that sustains robust responses by providing the gene promoter with a chain of Pol II complexes, we washed out DRB in the presence of Trip, a diterpene triepoxide isolated from natural sources. Trip inhibits Pol II recruitment to promoters by inhibiting TFIIH, the pre-initiation complex component that mediates Pol II binding to DNA (Titov et al., 2011). When DRB washout was performed in the presence of a low concentration of Trip, Arc transcription was significantly reduced at later time points (Fig. 8A) suggesting that the slightest interference with Pol II recruitment compromises the robustness of the response. Had H2A.Z.1 played a role in Pol II recruitment, its depletion would have had an effect comparable to washout with Trip in control cells. However, that does not seem to be the case as the Arc response is significantly different between these two groups (Fig. 8A, compare broken black line with green line). Moreover, control and H2A.Z.1-depleted neurons responded similarly to induction in the presence of Trip (Fig. 8A, compare broken black and green lines), suggesting that H2A.Z.1 is unlikely to be required for Pol II recruitment.
To test the roles of H2A.Z hypervariants in another neurobiological context, we next performed the TTX washout assay (Rao et al., 2006; Saha et al., 2011). This assay involves treating neurons with TTX for 48 h to induce homeostatic synaptic changes, followed by TTX washout to trigger rapid and robust network activity resulting in Arc transcription (Fig. 8B). Consistent with other induction assays so far, H2A.Z.2 depletion had an attenuating effect on early transcriptional response, followed by partial recovery at later time points (Fig. 8B). Surprisingly, depletion of H2A.Z.1 significantly reduced the Arc response throughout, even the early time points. This is in contrast with the insensitivity of the Arc response at early time points in all induction assays, including 15-min postinduction with Bic + 4AP under resting (control) conditions. Together, our data suggest that H2A.Z.1 exerts context-dependent effects on activity-induced Arc transcription.
Gene- and context-specific roles of H2A.Z hypervariants in rapid IEG transcription
In a previous study, we identified two groups of IEGs classified based on the presence or absence of paused Pol II near the promoter (Saha et al., 2011). We refer to these two groups as rapid IEGs (genes with paused Pol II) and delayed IEGs (genes without paused Pol II; Saha and Dudek, 2013). Because H2A.Z.2 likely regulates Arc transcription via RNA Pol II pausing, we next asked if the roles played by these hypervariants are universal across other rapid IEGs. Several rapid IEGs, identified in our previous screen, were then studied under two neurobiological contexts: (1) activity induced in previously untreated neurons (Bic + 4AP assays); and (2) activity induced after homeostatic changes (TTX washout assays). Responses of a few IEGs in control and H2A.Z hypervariant-depleted neurons are shown in Figure 9, and salient inferences drawn from these data are noted in the following paragraph.
H2A.Z.1 and H2A.Z.2 play context-dependent roles in activity-induced transcription of rapid IEGs. Graphical representation of pre-mRNA levels for noted IEGs at different time points after two treatment regimens. Left, Bic + 4AP treatment. Right, 48-h treatment with TTX followed by its washout for specified time periods in indicated groups of neurons as detected by qPCR and normalized to Gapdh. N = 4–7 for each time point; *p < 0.05 and **p < 0.01. NS, not significant. Two sets of data in statistical comparison are represented by the position and the color of asterisk(s).
All IEGs in Figure 9 respond to activity under both normal and homeostatic conditions. With some genes, such responses are more robust at early time points after 48 h of TTX treatment (Npas4, cFos, Cyr61, etc.). In general though, most rapid IEGs maintained similar temporal dynamics of expression in both conditions; Btg2 was an exception in that its response peaked at different times under different conditions. In most of these instances, depletion of either H2A.Z.1 or H2A.Z.2 had either no effect or reduced gene transcription. Fbxo33 was an exception to this observation, with H2A.Z.1 depletion resulting in a significantly stronger transcriptional response at the early time point in the Bic + 4AP condition only. Several genes (Cyr61, Btg2, and Dusp1) resembled Arc in that H2A.Z.2 depletion caused the most severe reductions. In contrast to Arc and these genes, H2A.Z.1 or H2A.Z.2 knockdown exerted statistically similar effect on several rapid IEGs (Npas4, cFos, Nr4a3, etc.) under both contexts. Interestingly, we noticed that H2A.Z hypervariant depletions induced different outcomes under different contexts in a few cases. For example, H2A.Z.2 depletion did not alter Btg2 and Dusp1 responses significantly at the early time point under normal conditions, but impaired them at the same time point after 48 h of TTX. Interestingly, these two genes are insensitive to H2A.Z.1 depletion for both contexts at all tested time points. Collectively, our data in Figure 9 suggest that H2A.Z.1 and H2A.Z.2 do not have generalized roles in rapid IEG transcription; instead, they play gene-specific roles, which may vary further under different neurobiological contexts.
To further investigate the role of H2A.Z hypervariants in IEG transcription, we conducted NanoString analyses, a multiplexed gene expression analysis system (Geiss et al., 2008), on control and hypervariant-depleted cells under the different neurobiological contexts described above using all IEGs identified in our previous screen (Saha et al., 2011). Figure 10A,B shows the results for rapid IEGs and non-IEGs (control genes) displayed as heat maps. Notably, results of the NanoString analyses were consistent with our initial qPCR analyses of rapid IEGs in Figure 9. Although there appeared to be gene-specific regulation due to loss of H2A.Z.1 and/or H2A.Z.2, overall comparison of all rapid IEGs tested in control neurons to those tested in hypervariant-depleted neurons suggests that loss of either hypervariant generally resulted in diminished transcription of rapid IEGs at each time point under both neurobiological contexts (Fig. 10C,D).
Effect of H2A.Z hypervariant depletion on transcription of rapid IEGs and nonactivity-induced genes. A, Nanostring multiplexed gene expression data (heat maps) representing pre-mRNA levels in control and hypervariant-depleted neurons treated with Bic + 4AP for indicated times. B, Heat maps representing pre-mRNA levels in control and hypervariant-depleted neurons treated TTX for 48 h followed by washout and collected at indicated times. C, D, Violin plots depicting pre-mRNA abundance of all rapid IEGs from hypervariant-depleted neurons compared to that of control neurons treated with Bic + 4AP (C) or TTX for 48 h followed by washout (D) for indicated times. N = 3 for each treatment all time points. *p < 0.05 and **p < 0.01 (Wilcoxon paired nonparametric test).
Gene- and context-specific roles of H2A.Z hypervariants in transcription of delayed IEGs
Because we saw gene-to-gene variation in H2A.Z hypervariant regulation of rapid IEG transcription, we asked if these hypervariants might also be involved in regulating transcription of delayed IEGs. Based on the lack of a role for paused Pol II in the delayed IEG transcription, we expected to find no effect of hypervariant depletion. However, we found a context-dependent regulation of delayed IEGs, and some trends are noted along the following lines: at 60-min postinduction, depletion of H2A.Z.1 appears to result in downregulation of Cox2, Klf4, and Nur77 with Bic + 4AP induction, yet results in upregulation of the same genes on induction after TTX-induced changes. Differential regulation between neurobiological contexts was also evident after loss of H2A.Z.2 in Cartpt and Fam46a (Fig. 11A). Interestingly, group statistical comparisons of delayed IEG responses in hypervariant-depleted neurons to control neurons show that loss of H2A.Z.1 only has a significant effect at the early time point on induction after 48 h of TTX. However, on induction in normal neurons (Bic + 4AP assays), H2A.Z.1 depletion only significantly affects the group at the late time point (Fig. 11B). Overall, these data suggest that the roles of H2A.Z hypervariants in regulating transcription of delayed IEGs may vary depending on specific genes and cellular states.
Effect of H2A.Z hypervariant depletion on transcription of delayed IEGs. A, Nanostring multiplexed gene expression data (heat maps) representing pre-mRNA levels in control and hypervariant-depleted neurons treated with Bic + 4AP for indicated times. B, Heat map representing pre-mRNA levels in control and hypervariant-depleted neurons treated TTX for 48 h followed by washout and collected at indicated times. C, D, Violin plots depicting pre-mRNA abundance of all delayed IEGs from hypervariant-depleted neurons compared to that of control neurons treated with Bic + 4AP (C) or TTX for 48 h followed by washout (D) for indicated times. N = 3 for each treatment at all time points. *p < 0.05 and **p < 0.01 (Wilcoxon paired nonparametric test).
Discussion
H2A.Z is enriched in nucleosomes near promoters and enhancers of many vertebrate genes and therefore is postulated to play a critical role in gene transcription (Farris et al., 2005; Raisner and Madhani, 2006; Zlatanova and Thakar, 2008). However, roles of the hypervariants H2A.Z.1 and H2A.Z.2 in these processes have often been presumed to be redundant, and have not been categorically questioned (Soboleva et al., 2014). In the current study, we have demonstrated several instances of H2A.Z.1 and H2A.Z.2 nonredundancies: (1) H2A.Z.1 and H2A.Z.2 depletion in cortical neurons affects basal expression of predominantly nonoverlapping sets of genes; (2) depletion of these hypervariants has differential effects on important synaptic molecules like PSD-95 and SHANK3; (3) H2B, perhaps when modified in a certain way, interacts with MYC-H2A.Z.2, but not HA-H2A.Z.1; (4) H2A.Z.1 and H2A.Z.2 depletion differently influence activity-induced Arc transcription kinetics; (5) in normal neurons, H2A.Z.2, but not H2A.Z.1, is required near the Arc TSS for RNA Pol II pausing; (6) activity-induced transcription of certain IEGs is sensitive to depletion of only one of the two hypervariants (e.g., Btg2 and Dusp1 response is less after H2A.Z.2, but not H2A.Z.1, knockdown); and (7) certain IEGs appear to be differentially sensitive to one or both hypervariants in a context-dependent manner (Fbxo33 and Maff after H2A.Z.1 depletion, Fam46a after loss of H2A.Z.2). Taken together, our data strongly suggest that H2A.Z hypervariants have important, nonredundant roles and/or functional specificity that should be considered for a full understanding of the part played by H2A.Z in gene transcription.
Three lines of evidence inspired our hypothesis about functional specificity of H2A.Z hypervariants. First, several instances have been observed in nature where genes duplicated during evolution encode very similar protein products that have achieved functional diversity. An important example is found in yeast where several paralogous ribosomal proteins are functionally dissimilar to their nearly identical counterparts (Komili et al., 2007; Steffen et al., 2012). Because mechanisms that work are often preserved in nature, we reasoned that the three amino acid differences in these hypervariants might be enough to allow functional diversity. Second, all three variant amino acid residues reside in functionally pertinent regions thereby making them suitable to mediate alternative functionalities. H2A.Z.114S/H2A.Z.214A is located in the H2A.Z N-terminal tail, which is heavily acetylated and is linked to H2A.Z transcriptional activity (Bruce et al., 2005; Valdés-Mora et al., 2012). H2A.Z.138T/H2A.Z.238S is located at the edge of the L1 loop, which plays an essential part in intranucleosomal interactions of H2A.Z (Horikoshi et al., 2013). Such intranucleosomal interactions are also facilitated by the unstructured region of H2A.Z C-terminus (Wratting et al., 2012), which houses the H2A.Z.112128V/H2A.Z.2128A variation. Third, knockout of H2A.Z.1 results in embryonic lethality in mice (Faast et al., 2001) indicating that H2A.Z.2 cannot compensate for the loss of its paralog. Together, these observations suggest that H2A.Z hypervariants have nonredundant functions that have not yet been adequately explored.
Our study into functional specificity of H2A.Z hypervariants in activity-induced Arc transcription revealed a few interesting aspects of H2A.Z biology. Overall, Arc transcription requires both H2A.Z.1 and H2A.Z.2. Depletion of either hypervariant had unique effects on Arc transcriptional kinetics that were not compensated by the other paralog. In normal previously untreated neurons (i.e., without homeostatic changes induced by silencing), H2A.Z.2, but not H2A.Z.1, was found to be necessary for the rapid phase of activity-induced Arc transcription (within 15 min of treatment) by facilitating RNA Pol II promoter proximal pausing (Fig. 7). Pol II pausing near promoters is facilitated, among other factors, by the position of the +1 nucleosome (Kwak and Lis, 2013; Jimeno-González et al., 2015). Although exact mechanisms remain unknown, it is presumed that the +1 nucleosome acts as a physical barrier to the transcriptionally engaged Pol II. In mammals, the presence of H2A.Z in the +1 nucleosome is positively correlated with Pol II pausing (Ku et al., 2012; Day et al., 2016). In our study, considering that (1) H2A.Z.2 is required for Pol II pausing near the Arc TSS and (2) Bic + 4AP treatment removes MYC-H2A.Z.2 from the +1 nucleosome (Fig. 6), we suggest that H2A.Z.2-containing nucleosomes play a dual role in mediating activity-induced transcription. In the resting state before activity, the +1 H2A.Z.2-containing nucleosome would maintain a rigid barrier to Pol II, thereby facilitating its pausing. After activity, dynamic turnover of this nucleosome allows Pol II to commence productive elongation. If true, such a dual-role model could fuel further queries about possible activity-induced signals that convert a rigid nucleosomal barrier into a dynamic gate for Pol II productive elongation. Although this model holds true for activity induced under both normal and prolonged TTX treatment-enabled conditions for Arc, further investigation into other IEGs revealed the potential for both gene- and context-specific roles of H2A.Z hypervariants in gene transcription. However, future studies will be required to determine if these hypervariants are directly or indirectly involved in the transcription of other IEGs investigated in this study.
One of our observations was puzzling and therefore difficult to interpret: although there was turnover of MYC-H2A.Z.2 in response to acute increases in neuronal activity, the total levels of H2A.Z (as reported by the Millipore anti-H2A.Z antibody) near the Arc TSS did not change over the entire duration of testing. This could be due to weak efficiency of the commercial H2A.Z antibody (Fig. 6B, compare signals), which was the best anti-H2A.Z ChIP antibody in our hands. We believe that, as has been previously reported (Zovkic et al., 2014) and is also shown by our anti-MYC ChIP data, activity does indeed induce H2A.Z turnover at the Arc TSS.
In comparison to the Pol II pausing-related defined role of H2A.Z.2 in Arc transcription, the exact role of H2A.Z.1 in the process is less apparent. In normal, unsilenced neurons, lack of H2A.Z.1 impairs only the late segment of activity-induced Arc transcription (30 min and later). As suggested by our DRB washout assays, this is most likely via indirect mechanisms (discussed further below). However, after subjecting neurons to homeostatic plasticity-related changes with TTX exposure for 48 h, H2A.Z.1 depletion also impairs early responses of Arc transcription. This repurposing of H2A.Z.1 after prolonged lack of activity is likely mediated by a context-dependent newfound ability to facilitate Pol II pausing. Alternatively, these data could also be explained by heterogeneity of nucleosomal composition in the cell population (discussed further below).
Considering that we could not detect strong interactions of tagged H2A.Z.1 and H2A.Z.2 in our co-IPs (Fig. 5), indicating that they likely do not often coexist in the same nucleosome, and yet show H2A.Z.1 and H2A.Z.2 enrichment at the same nucleosomal positions (Fig. 6B), it is possible that these Arc TSS nucleosomal loci are populated by different H2A.Z hypervariants in different alleles and/or cells in the population. In normal neurons, only the fraction of these alleles and/or cells that have H2A.Z.2, and therefore paused Pol II, responds to activity. After prolonged TTX treatment, additional alleles and/or cells with H2A.Z.1 attain Pol II pausing via unknown mechanisms and facilitate an enhanced response after TTX washout. This explanation aligns with the noticeably more robust expression of Arc after prolonged TTX treatment and washout in comparison to disinhibition-induced activity under previously untreated conditions (compare Figs. 4B, 8B). Also, the effect of H2A.Z.1 depletion on late, not early, Arc transcriptional response (Fig. 4B,D) falls within the scope of this explanation. On induction of activity in resting conditions, H2A.Z.1-containing alleles and/or cells undertake Arc transcription in a Pol II pausing-independent “delayed” fashion, as delayed IEGs do (Saha et al., 2011), thereby accounting for impairments at later time points in H2A.Z.1-depleted neurons. These are exciting possibilities that will require further probing.
Additional experiments will also be required to understand the relevance of the H2A.Z chaperone ANP32E in regulating H2A.Z functions (Mao et al., 2014; Obri et al., 2014). In our hands, several traditional “rescue” attempts involving expression of shRNA-insensitive H2A.Z.2 failed to reverse the loss of Arc transcription in H2A.Z.2-depleted neurons. However, the phenotype loss could be partially rescued by codepletion of ANP32E. Because ANP32E removes H2A.Z from nucleosomes (Mao et al., 2014; Obri et al., 2014; Gursoy-Yuzugullu et al., 2015), it is likely that the presence of ANP32E in H2A.Z.2-depleted neurons facilitates continued removal of nucleosomal H2A.Z.2 and subsequent fallouts thereof, including loss of Pol II pausing and inability to respond rapidly. However, codepletion of the chaperone along with the histone preserves nucleosomal H2A.Z.2, at least in part, enabling Arc to rapidly respond to some extent. A similar precedence of H2A.Z and ANP32E codepletion partially rescuing phenotype loss due to depletion of one of these factors exists in non-neuronal cells (Alatwi and Downs, 2015). Our findings, combined with this previous report, advocate for the need to approach the complex biology of H2A.Z hypervariants and their chaperones together.
List of DEGs after H2A.Z.1 and H2A.Z.2 KD-part B (downregulated after H2A.Z.1 KD)
List of DEGs after H2A.Z.1 and H2A.Z.2 KD-part C (upregulated after H2A.Z.2 KD)
List of DEGs after H2A.Z.1 and H2A.Z.2 KD-part D (downregulated after H2A.Z.2 KD)
Collectively, our current studies indicate that H2A.Z hypervariants mediate activity-induced IEG transcription in a context-dependent (H2A.Z.1) or context-independent (H2A.Z.2) fashion, where the extent of their involvement in the transcriptional process varies from gene to gene. Our data also suggest spatial and temporal cooperation of these hypervariants in manifesting the maximum gene response under a given circumstance. Natural extensions of the current work will be to test the above hypothesis about allele/cell-specificity of these hypervariants, and to investigate molecular mechanisms connecting H2A.Z hypervariants and paused Pol II. Additionally, sequentially mutating H2A.Z.1 to H2A.Z.2 could provide insight into which amino acid differences are necessary for functional specificity. Particular attention may be paid to the H2A.Z.114S/H2A.Z.214A variation, where H2A.Z.1, but not H2A.Z.2, could be potentially phosphorylated. It will also be of interest to learn about hypervariant-specific H2A.Z chaperones, if any, and understand how these hypervariants differentially regulate expression of several crucial synaptic molecules. Finally, it is our hope that the data presented here will increase appreciation for functional specificity of H2A.Z hypervariants and promote tool-building efforts to enable a comprehensive understanding of diverse roles played by these paralogs at different genes and under different neurobiological contexts.
Acknowledgments
Acknowledgements: We thank Ms. A. Funk and Dr. B. Bennett for project assistance, the Saha lab members for critical reading of this manuscript, and Dr. Mike Cleary and Dr. Sergei Nechaev. We also thank the personnel in the Viral Vector and Molecular Genomics cores at National Institute of Environmental Health Sciences, and the Stem Cell Instrumentation Foundry at University of California Merced for use of their facility.
Footnotes
The authors declare no competing financial interests.
This work was supported by the United States Department of Health and Human Services (HHS) National Institutes of Health (NIH) National Institute of Mental Health Grant R00MH096941 (to R.N.S.), the HHS NIH National Institute of Environmental Health Sciences Grant ES100221 (to S.M.D.), and by start-up funds from University of California Merced (R.N.S.). For the K99 phase of the R00 Award, this research was supported by the Intramural Research Program of the National Institute of Environmental Health Sciences, National Institutes of Health.
This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International license, which permits unrestricted use, distribution and reproduction in any medium provided that the original work is properly attributed.
References
In this issue
