a Identification of the BRD1 interaction network and its impact on mental disorder risk

Tue Fryland, Jane H. Christensen, Jonatan Pallesen, Manuel Mattheisen, Johan Palmfeldt and Mads Bak

Genome Medicine. 8.53 (May 3, 2016): From InfoTrac Health and Medical Collection 2017.

Author(s): Tue Fryland^1,2,3 , Jane H. Christensen^1,2,3 , Jonatan Pallesen^1,2,3 , Manuel Mattheisen^1,2,3, Johan Palmfeldt⁴ , Mads Bak⁵ , Jakob Grove^1,2,3,6 , Ditte Demontis^1,2,3 , Jenny Blechingberg¹ , Hong Sain Ooi^1,2,3 , Mette Nyegaard^1,2,3 , Mads E. Hauberg^1,2,3 , Niels Tommerup⁵ , Niels Gregersen⁴ , Ole Mors^2,3,7 , Thomas J. Corydon^1,2,3 , Anders L. Nielsen^1,2,3 and Anders D. Barglum^1,2,3,7

Background

Genetic variations in the bromodomain containing 1 (BRD1) gene located on 22q13.33 have repeatedly been associated with both schizophrenia and bipolar disorder [1-3]. The BRD1 locus approached genome-wide significance ( P = 3.31 x 10^-7 ) in the most recent schizophrenia genome-wide association study (GWAS) by the Psychiatric Genomics Consortium using conventional statistical methods [4]. Moreover, the locus was found genome-wide significant using an Empirical Bayes statistical approach and predicted to be highly replicable [5]. Furthermore, in a schizophrenia GWAS meta-analysis (11,185 cases/10,768 controls) and a family-based replication study (3286 cases from 1811 families) the SNP rs138880 in the promoter region of BRD1 showed the overall most significant association [6].

BRD1 is essential for normal brain development and inactivation of both alleles of Brd1 in mice leads to impaired neural tube closure [7]. Co-immunoprecipitation (co-IP) of epitope tagged and endogenous BRD1 and MYST2 from human K562 and HEK293 cells suggest that ING4, MEAF6, and MYST2 constitute the primary histone acetyltransferase complex of BRD1 [7]. Additionally, a focused promoter ChIP-on-chip (chromatin immunoprecipitation combined with microarray analysis) of co-expressed epitope tagged BRD1 and MYST2 in human K562 cells identified a large overlap in target genes between the two proteins suggesting a pivotal role of the BRD1/MYST2 complex in transcriptional regulation [7]. Equally, Brd1 ^-/- and Myst2 ^-/- knockout mouse embryos show severely decreased levels of the overall histone H3K14 acetylation, suggesting that the Brd1/Myst2 complex is responsible for the majority of histone H3K14 acetylation during embryonic development in mice [7, 8]. We have previously reported a divergence in immunohistochemical staining between the BRD1-S and BRD1-L isoforms in the brain [9] as well as differential transcriptional regulation of theBrd1-S and Brd1-L splice variants in prefrontal cortex and hippocampus following chronic restrained stress [10] and electroconvulsive seizures [11] in adult rats, indicating that BRD1 isoforms can perform separate functions dependent on the specific cell type and tissue. To gain more knowledge about the biological functions of BRD1 and how these might be involved in the pathogenesis of schizophrenia and related mental disorders, we sought in the present study to identify and analyze the BRD1 interaction network, encompassing BRD1-S and BRD1-L protein-protein interactions (PPIs) and chromatin interactions as well as genes being regulated upon up- or downregulation of BRD1. Moreover, we interrogated large GWAS datasets and found that the BRD1 interaction network is enriched for schizophrenia risk.

Methods

Cell work

The generation of cell lines stably expressing BRD1-S-V5 and BRD1-L-V5 have previously been described [9]. HEK293T cells (controls and stable BRD1-S-V5 and BRD1-L-V5 cell lines) were grown in DMEM medium (Invitrogen, San Diego, CA, USA) supplemented with 5 % fetal calf serum (FCS), 175 mg/L glutamine, 36 mg/L penicillin, and 60 mg/L streptomycine at 37 °C in 5 % CO₂ .

Co-immunoprecipitation (Co-IP)

Preparation of cell extract was performed according to the two-step procedure described in [12]. Experiments were carried out in 10 cm or 15 cm petri dishes with 1 x 10⁷ cells or 2 x 10⁷ cells plated, respectively. 1 x 10⁸ cells were used for each immunoprecipitation (IP). Cells were counted using a Nucleocounter (ChemoMetec A/S, Alleroed, Denmark) and plated 24 h before harvested using 1 mL per 10 x 10⁶ cells hypotonic Triton X-100 lysis buffer (20 mM Tris-HCl [pH 7.4], 10 mM KaCl, 10 mM MgCl₂ , 2 mM EDTA, 10 % glycerol, 1 % Triton X-100, 2.5 mM [beta]-glycerophosphate, 1 mM NaF, 1 mM DTT + protease inhibitors (Roche, Mannheim, Germany]) for 10 min on ice. Cell lysate was distributed to 15 mL tubes with 2 mL in each for sonication. DNA was fragmented by sonication (Bioruptor, settings: on 0.5, off 0.5) for 15 min at 6 °C. A total of 5 M NaCl was added to a final concentration of 420 mM, mixed and incubated on ice for 15 min after which the DNA fragmentation was repeated. Sonicated cell lysate was then cleared by centrifugation at maximum speed for 15 min and the supernatant was recovered for IP.

IP of V5 epitope tagged proteins was performed as follows: Anti-V5 and anti-HA antibody conjugated agarose beads (Sigma Aldrich, Steinheim, Germany) were washed twice in PBS before use and blocked in 1 % BSA. Cell lysates were pre-cleared for 30 min at 4 °C under rotation using 30 [mu]L anti-HA antibody conjugated bead solution before collecting the supernatant (400 rpm for 5 min) for IP. IP was performed using 60 [mu]L of the bead solution per 10 mL lysate. Beads were added to the pre-cleared supernatant and incubated at 4 °C overnight before they were collected at 400 rpm for 15 min. Beads were washed twice in 10 mL Triton X-100 buffer 250 (20 mM Tris-HCl [pH 7.4], 250 mM NaCl, 10 mM MgCl₂ , 2 mM EDTA, 10 % glycerol, 1 % Triton X-100, 2.5 mM [beta]-glycerophosphate, 1 mM NaF, 1 mM DTT, protease inhibitors [Roche]) then moved to Protein LoBind tubes (Eppendorf, Hamburg, Germany) and then washed six times in Triton X-100 buffer 150 (20 mM Tris-HCl [pH 7.4], 150 mM NaCl, 10 mM MgCl₂ , 2 mM EDTA, 10 % glycerol, 1 % Triton X-100, 2.5 mM [beta]-glycerophosphate,1 mM NaF, 1 mM DTT, protease inhibitors [Roche]). Proteins were eluted in 100 [mu]L Elution buffer (HENG) (10 mM HEPES-KOH [pH 9.0], 1.5 mM MgCl2, 0.25 mM EDTA, 20 % glycerol, 250 mM KCl, 0.3 % NP40 + 0.5 mg/mL V5 peptide [Sigma], protease inhibitors [Roche]). Cell lysate was kept cold at all times, between 1 °C and 6 °C. The process was carried out in one continuous procedure without any freezing/thawing of the cell lysate. IP with anti-V5 antibody and IP with anti-HA antibody was performed four times with either extracts from HEK293T, stable BRD1-S-V5, or BRD1-L-V5 expressing cell lines, adding to a total of 24 IP samples prepared for mass spectrometry.

IP of endogenous PBRM1 was performed using the same procedure as described above. Protein A and G beads (GE Healthcare, Uppsala, Sweden) were washed twice in PBS and blocked in 1 % BSA (final concentration). Pre-clearing and IP was performed in a 50/50 [mu]L mix of blocked protein A and G bead solutions. A total of 30 [mu]L of protein A and G bead solution was used for preclearing of cell extracts and 60 [mu]L of protein A and G bead solution was used for IP. Anti-PBRM1 antibody (Cat. no. A301-590A, Bethyl, Montgomery, AL, USA) was diluted in Triton X-100 buffer 150 according to antibody specifications and incubated with blocked protein A and G bead solution at 4 °C overnight. IP control reactions were performed in parallel by the same method but without the addition of antibody. Furthermore, elution was performed with 1 % SDS instead of the V5 peptide solution as described above.

Sample preparation for mass spectrometry

Proteins eluted upon IP were precipitated in six volumes of acetone (-20 °C), incubated at -20 °C overnight, and centrifuged at 15,000 x g at 4 °C. The supernatant was cautiously removed and the protein pellet was carefully washed in ice-cold 90 % acetone, centrifuged at 15,000 g and air-dried. Disulfide bonds of the proteins were reduced in 90 [mu]L buffer containing 50 mM tris-(2-carboxyethyl)phosphine (TCEP) (Sigma Aldrich, Steinheim, Germany) and 50 mM ammonium bicarbonate, and incubated at 60 °C for 10 min. The reduced cystein residues were blocked by 50 mM of the alkylation agent iodoacetamide (Sigma-Aldrich, Steinheim, Germany). Proteins were digested into peptides by addition of 1 [mu]g trypsin (Trypsin Gold from Promega, Madison, WI, USA) in 100 [mu]L 50 mM ammonium bicarbonate.

Mass spectrometry analysis

The peptide mixtures were analyzed by nanoliquid chromatography (Easy nLC from Proxeon, Denmark) coupled to mass spectrometry (MS) (Thermo Fisher, LTQ-Orbitrap), with separation on a reverse phase column (75 [mu]m, 100 mm, and 3.5 [mu]m C18 particles) at a flow rate of 300 nL/min using a 100 min gradient (5-35 %) of acetonitrile in 0.4 % acetic acid. The MS detection constituted a full scan (m/z 400 - 2000) in Orbitrap (<3 ppm mass accuracy) followed by up to four data dependent MS/MS fragmentation scans using collision induced dissociation (CID).

Database searches and MS statistics

The resulting MS files were processed essentially as previously described [13]. Briefly, Mascot version 2.2.04 (Matrix Science) was used for peptide identification and MaxQuant version 1.0.13.18 for protein identification and label-free quantification [14]. The MS data were searched against IPI protein database version 3.52 containing 73,928 sequences and the same number of reversed sequences for false discovery rate calculations (FDR). FDR was set to 0.01 for both identification of peptides and proteins. MS/MS mass tolerance was 5 ppm for peptide masses and 0.5 Da for fragment masses. Setting of trypsin digestion was cleavage at C-terminus of lysine and arginine except before proline and up to two missed cleavages were accepted. Carbamidomethylation at cysteine residues was set as fixed modification and oxidation of methionine was set as variable modification. Only peptides with a minimum length of six amino acid residues were accepted and at least two peptides (and one unique peptide) were required for protein identification.

Proteomics data analysis

The combined spectral intensity of each individual protein was normalized to the mean intensity of all the proteins for each sample or control. Intensities below the detection limit were assigned the value of mean of the minimum intensities from across all datasets, and then normalized with the mean intensity of all values in the sample or control. Extracts from stable BRD1-S-V5 and BRD1-L-V5 expressing cells and from HEK293T cells were used in co-IP experiments. Each co-IP experiment was performed four times followed by LC-MS/MS analysis, using beads with anti-V5 conjugated antibody (IP:V5) or anti-HA conjugated antibody (IP:HA) and extracts from stable BRD1-S-V5, BRD1-L-V5 expressing cells, or HEK293T cells, yielding: (4 x IP:V5 4 x BRD1-S-V5, 4 x IP:HA BRD1-S-V5, 4 x IP:V5 BRD1-L-V5, 4 x IP:HA BRD1-L-V5, 4 x IP:V5 HEK293T, 4 x IP:HA HEK293T). Four criteria (1-4) were raised to systematically identify repeatedly co-immunoprecipitated (co-IPed) proteins by calculating the enrichment ratio (log₂ (sample) - log₂ (control)), where control experiment is IP without proper BRD1 in the sample or without antibody against the proper BRD1 construct. The four criteria were: (1) the enrichment ratio IP:V5 BRD1-(S or L)-V5 - IP:HA BRD1-(S or L)-V5 of a protein is > cutoff value; (2) the enrichment ratio IP:V5 BRD1-(S or L)-V5 - IP:V5 HEK293T of a protein is > cutoff value; (3) the enrichment ratio IP:V5 BRD1-(S or L)-V5 - IP:HA HEK293T of a protein is > cutoff value; (4) the enrichment ratio IP:V5 BRD1-(S or L)-V5 - the mean of all controls of a protein is > cutoff value. Both two- and threefold cutoffs were applied to sort the identified protein-protein interactions according to co-IP strength (Table 1).

Table 1: Summary of BRD1 PPIs identified by co-IP mass spectrometry [see PDF for image]

Western blotting

Protein samples were prepared for SDS-PAGE using 5x loading buffer and 20x reducing agent (Fermentas, St. Leon-Rot, Germany) and then placed in a heating block for 5 min at 100 °C. Denatured proteins were loaded onto TGX or Tris-HCL 4-15 % linear gradient polyacrylamide gels (Bio-Rad Laboratories) and separated by electrophoresis at 100-120 V. Proteins were transferred to Hybond-P Membranes (GE Healthcare, Uppsala, Sweden) using a wet blotting apparatus (Bio-Rad Laboratories). Membranes were activated in 99 % alcohol for 20 s, washed in redistilled water for 1 min, and transferred to transfer buffer (Tris 50 mM Tris, 385 mM glycine, 7 mM SDS, 10 % Ethanol). Membranes were blocked for 1 h in blocking buffer (10 % skimmed milk powder, 1 % Tween 20, and PBS). Primary antibodies were: monoclonal anti-V5 antibody (Invitrogen, R960-25) diluted 1:2000 in wash buffer (0.05 % Skimmed milk powder, 0.005 % Tween 20, and PBS) and anti-BRD1 antibody (65/66 Sigma Genosys), previously described [9]. Secondary antibodies were: HRP conjugated anti-mouse IgG antibody (Dako Cytomation, P0447, Glostrup, Denmark) and HRP conjugated anti-rabbit IgG antibody (Dako Cytomation, P0448, Glostrup, Denmark). Dilutions were according to the manufacturer's protocol. All incubations were performed in 5-10 mL wash buffer. Membranes were incubated with primary antibodies for 24 h at 4 °C, followed by 3x 5 min wash and for 1 h with secondary antibody at 4 °C and washed 5x 5 min with wash buffer. Immunoreactivity was detected using BM Chemoluminescence substrate (POD) (Roche Applied Science, Mannheim, Germany) according to the manufacturer's protocol and the Image Reader LAS-3000 v2.2 program (Fujifilm, Minato (Tokio), Japan). Densitometric analysis was performed with the Multi Gauge v3.1 software (Fujifilm, Minato (Tokio), Japan).

Immunofluorescence

HEK293T cells stably expressing BRD1-S and BRD1-L and untransfected HEK293T cells were cultured for 48 h in 10 cm² slideflasks. Cells were fixed in 4 % freshly prepared paraformaldehyde, 0.1 % gluteraldehyde, 0.1 % Trition X-100 solution on ice for 30 min, followed by wash with PBS for 5 min. Slides were incubated with anti-V5 antibody (Invitrogen) for 1 h at 4 °C using a 1:200 dilution in PBS and washed 3x 5 min with PBS containing 0.1 % Triton X-100. Incubation with anti-mouse immunoglobulins/FITC F(ab')2 (Dako, Glostrup, Denmark) was performed for 1 h at 4 °C using a 1:200 dilution in PBS supplemented with 5 % FCS. Cell nuclei were stained by Hoechst staining (Hoechst stain 1 [mu]g/mL) for 10 min, slides were then washed 2x 5 min in PBS containing 1 % Triton X-100 and dried. Two drops of fluorescent mounting media were added to the fixed cells and stored in the dark at 4 °C until media were dried out and then viewed under a fluorescent microscope (Additional file 1).

Chromatin immunoprecipitation (ChIP)

ChIP of V5 epitope proteins were carried out as described in [15] with a few adjustments. Eight million cells per reaction were seeded 24 h before use in 10 mL medium. Formaldehyde, to the final concentration of 1 %, was added directly to the cell medium, on a shaker. The cross-linking reaction was stopped after 10 min by addition of glycine to the final concentration of 0.125 M, on a shaker. Cells were carefully washed twice in 10 mL cold PBS, harvested in 1 mL ChIP dilution buffer (0.01 % SDS, 1.1 % Triton X-100, 1.2 mM EDTA, 16.7 mM Tris-HCL, 167 mM NaCl) and protease inhibitors (Complete mini, Roche, Mannheim, Germany) and sonicated on ice (Bioruptor, settings: on 0.5, off 1.0) for 15 min at 6 °C. Cell lysate was cleared by 2 min maximum centrifugation at 6 °C. Supernatant was collected and stored at -80 °C.

Immunoprecipitations were carried out using 20 [mu]L of anti- V5 or anti- HA conjugated agarose beads (Sigma, A7345, A2095). First, the beads were washed twice in PBS and protease inhibitors (Roche, Mannheim, Germany) with intermediate low g centrifugation. Beads were blocked in 10 x volume PBS, 200 [mu]g/mL sonicated herring DNA, and 1.5 % BSA at RT and washed in ChIP dilution buffer. Extracts were pre-cleared with 20 [mu]L of blocked anti-HA beads at rotation for 30 min before immunoprecipitation with blocked anti-V5 and anti-HA beads over night at 6 °C. Beads were collected by centrifugation and washed twice in low salt buffer (0.1 % SDS, 1 % Triton-X 100, 2 mM EDTA, 20 mM Tris-HCL, pH 8.1, 150 mM NaCl, and proteinase inhibitors [Roche]), twice in high buffer (0.1 % SDS 1 % Triton-X 100, 2 mM EDTA, 20 mM Tris-HCL pH 8.1, 500 mM NaCl, and proteinase inhibitors [Roche, Mannheim, Germany]), once in LiCl buffer (0.25 M LiCl, 1 % IGEPAL-CA630, 1 % deoxycholixacid sodium salt, 1 mM EDTA, 10 mM Tris-HCl pH 8.1, and proteinase inhibitors [Roche, Mannheim, Germany]), and twice in TE buffer. Elution was performed in 500 [mu]L freshly made elution buffer (1 % SDS, 0.1 M NaHCO ₃ , H₂ O). A total of 20 [mu]L of 5 M NaCl was added to eluates and incubated for 4 h at 65 °C for reversing cross-links. A total of 10 [mu]L of 5 M EDTA, 20 [mu]L of 1 M Tris-HCl pH 6.5, and 2 [mu]L of proteinase K (Finnzymes, F-202S, 20 mg/mL) was added to eluates and incubated for 1 h at 45 °C. Phenol/chloroform and isopropanol (v/v) was added, mixed, and centrifuged. A total of 2 [mu]L of glycogen (Sigma Aldrich, cat. no. G1767) and 0.7 x volume isopropanol was added to supernatant and incubated at -20 °C overnight. Precipitated DNA was collected by centrifugation at maximum speed for 30 min at 4 °C. Pellet was carefully washed in 70 % ethanol, air-dried, and dissolved in 50 [mu]L redistilled water. ChIP DNA concentrations were measured on a microplate reader (Thermo Fisher, Flouroskan Ascent FL) by fluorescence of double-stranded DNA (dsDNA) (Lifetechnologies/Invitrogen, Qant-iT⢠picogreen, cat. no. P7589). All ChIP DNA was diluted to the same concentration before quantitative real-time PCR.

ChIP sequencing

DNA libraries were prepared using the ChIP-seq DNA Sample Prep kit (Illumina; IP-102-1001) and sequenced on a Genome Analyzer IIx. The sequence reads were aligned to the human genome (hg19), through the Galaxy project portal [16], using the Illumina Analysis Pipeline. Sequenced reads were mapped with Bowtie [17] allowing one mismatch. Peak calling was performed using the Model-based Analysis of ChIP-Seq v1.4 (MACS) [18], through the Galaxy project portal. Settings: distance = 100; and otherwise standard settings. Ensemble (GRCh37.p13) known genes and transcripts were used as reference gene annotations to identify promoter target genes (PTGs).

RNA and complementary DNA from cell culture

Cultured cells were lysed using Qiashredder columns (Qiagen, Hilden, Germany) and RNA was purified using the RNAeasy Mini kit (Qiagen) according to the manufacturer's protocol. RNA concentrations were determined using a NanoDrop 1000 version 3.7.1 (Thermo Fisher Scientific, Waltham, MA, USA). RNA integrity was assessed by 1 % agarose gel electrophoresis. A complementary DNA (cDNA) library was generated from 1 [mu]g of RNA using the iScript⢠cDNA Synthesis Kit (Bio-Rad Laboratories, Hercules, CA, USA). All steps were performed according to the manufacturer's protocol using a mixture of random hexamers and Oligo-dT primers provided with the kit.

Microarray expression profiling

Microarray (Affymetrix U219) and RNA quality control was performed at Aros Biotechnology, Skejby. Further information is provided elsewhere [19, 20]. A threshold of 1.5 was set for the fold change of microarray probes (Additional file 2).

Quantitative real-time PCR

Quantitative real-time PCR was performed on the LightCycler®480 (Roche Applied Science, Mannheim, Germany) system. A total of 96 or 384 well plates, with a respective total reaction volume of 10 [mu]L or 5 Îl, were used, respectively. Fluorescence of dsDNA was determined by addition of SYBR®Green (Roche Applied Science, Mannheim, Germany) to the reactions in the concentration described by the manufacturer. The standard settings on the LightCycler software were used for excitation and detection of fluorescence. Primer3 [21] was used in the design of primers (MWG Operon, Ebersberg, Germany).

Gene set enrichment analysis

We used publicly available summary statistics from single-marker GWASs [4, 22-31] considering only variants outside the broad MHC-region (chr6:25 M-35 M) [28] and filtered for info score [greater than or equai to]0.8 if available.

Data on schizophrenia, bipolar disorder, anorexia nervosa, autism, MDD, ADHD, and cross-disorders were accessed and downloaded via the PGC website [32]. Data on rheumatoid arthritis were downloaded from [33]. Data on coronary artery disease/myocardial infarction have been contributed by CARDIoGRAMplusC4D investigators and have been downloaded from [34]. Data on type 2 diabetes have been contributed by DIAGRAM and downloaded from [35]. Data on Crohn's disease have been contributed by IIBDG and downloaded from [36]. Data on Alzheimer's were downloaded from [37] - International Genomics of Alzheimer's Project (IGAP) is a large two-stage study based upon GWAS on individuals of European ancestry. In stage 1, IGAP used genotyped and imputed data on 7,055,881 single nucleotide polymorphisms (SNPs) to meta-analyze four previously published GWAS datasets consisting of 17,008 Alzheimer's disease cases and 37,154 controls (The European Alzheimer's Disease Initiative - EADI; the Alzheimer Disease Genetics Consortium - ADGC; The Cohorts for Heart and Aging Research in Genomic Epidemiology Consortium - CHARGE; The Genetic and Environmental Risk in AD Consortium - GERAD). In stage 2, 11,632 SNPs were genotyped and tested for association in an independent set of 8572 Alzheimer's disease cases and 11,312 controls. Finally, a meta-analysis was performed combining results from stages 1 and 2. For the risk gene enrichment analysis, we used MAGMA [38] and default settings. Genes were annotated using Ensemble (GRCh37.p13, same reference for identifying BRD1 PTGs). Information about the genetic correlation pattern in the data (linkage disequilibrium) was obtained using the 1000 Genomes European panel [39]. To assess whether or not the BRD1 protein network was enriched with de novo mutations with relevance to autism risk [40, 41] or rare mutations with relevance to schizophrenia risk [42], we compared the proportion of mutations in the genes among cases to controls (obtained from the referenced studies) using a one-sided binominal test correcting for the overall ratio of mutations in cases compared to controls.

Identification of spatiotemporal networks in the human brain

Human brain transcriptome (RNA-seq) data were obtained from www.brainspan.org [43]. We filtered the dataset including only genes with coefficient of variance more than 0.1 thereby removing genes with little or no information with regards to the spatiotemporal dynamics. Since the BRD1 network contained both protein-coding and non-protein-coding genes, both types of genes were included in the analysis yielding 21,396 uniquely identified genes. Data were segregated into 32 spatiotemporal intervals consisting of eight temporal intervals (P1-P8) and four brain regions (R1-R4). Brain regions (R) were grouped according to transcriptional similarities based on the hierarchical clustering described here [44]. Temporal intervals (P) were grouped as follows: P1 included post-conceptual week (pcw) 8-13 (first trimester); P2 included 16-26 pcw (second trimester); P3 included 35-37 pcw (third trimester), P4 included 4-10 months; P5 included 1-4 years; P6 included 8-13 years; P7 included 15-19 years; P8 included 21-40 years. The fraction of co-expressed genes (correlation coefficient >0.5) was calculated for each spatiotemporal interval and each BRD1 sub-network and compared to the background of genes co-expressed with BRD1 in the entire dataset. To determine whether or not the fraction of co-expressed genes was significantly higher than the background, a one-sided binominal test was performed while adjusting the P value for the number of tests performed (Bonferroni correction).

Software and statistical analysis

Ingenuity Pathway analysis was used for the bioinformatics analysis. Simulations, statistical analysis, and graphical illustrations were conducted in R, Python, Excel (Microsoft, Office 2010, USA), and GraphPad Prism 5.

Availability of data

The ChIPseq and microarray data discussed in this publication have been deposited in NCBI's Gene Expression Omnibus [45] and are accessible through GEO Series accession numbers GSE62811 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE62811) and GSE79255 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE79255). The protein-protein interactions from this publication have been submitted to the IMEx consortium [46, 47] through IntAct [48] and assigned the identifier IM-24990.

Results

Identification of BRD1 protein-protein interactions

In order to identify PPIs of the two BRD1 isoforms (BRD1-S and BRD1-L), we generated HEK293T cells stably expressing either V5 epitope tagged BRD1-S or BRD1-L (referred to as BRD1-S and BRD1-L cells). Both epitope tagged BRD1 isoforms localized to the nucleus of the HEK293T cells (Additional file 1), suggesting that the V5 epitope does not disrupt nuclear localization of the proteins. By co-immunoprecipitations, followed by identification of associated proteins using nanoLC coupled to a LTQ Orbitrap MS/MS system (Fig. 1a-c), we identified 24 and 16 PPIs for BRD1-S and BRD1-L, respectively (Table 1). The previously reported BRD1 interaction partners MEAF6, ING4, and ING5 [7, 49] were identified to interact with both BRD1-S and BRD1-L, whereas MYST2 [7] was found to interact solely with BRD1-S (Fig. 1d and Table 1). Likewise, our finding that histone H3 interacts with BRD1-L (Table 1) is in line with a previous study that identified the unmodified N-tail of histone H3 to interact with the PHD1 domain of BRD1 [50]. We cross-referenced the identified protein networks with public databases (BioGRID, CCSB Human Interaction Network Database, DIP, HPRD, IntAct, MINT, BIND, MPPI, and GNP) and discovered 20 novel BRD1-S PPIs and 13 novel BRD1-L PPIs which have not previously been reported, some of which were identified solely for one isoform (Fig. 1d and Table 1). A recent large-scale protein-protein interaction study found that the majority of protein isoforms share less than 50 % of their interaction partners [51]. In line with these results, 62.5 % of the identified PPIs were found to be common among the two BRD1 isoforms. As an example, the Histone-lysine N-methyltransferase protein suppressor of variegation 4-20 homologue 1 (SUV420H1) was identified as a specific PPI for BRD1-S (Fig. 1d and Table 1). As part of the list of common BRD1 PPIs (i.e. interacting with both BRD1-S and BRD1-L), we identified the 14-3-3 tyrosine monooxygenase proteins YWHAE, YWHAH, YWHAZ, and the poly bromo 1, PBRM1 (Table 1). These PPIs are of particular interest since genetic association studies have previously implicated YWHAE , YWHAH , and YWHAZ with schizophrenia and bipolar disorder [52-55] and thePBRM1 locus has surpassed genome-wide significance in both schizophrenia and bipolar disorder GWASs [4, 56-58]. We confirmed the interaction between BRD1 (for both isoforms) and PBRM1, by co-IP with antibodies against endogenous PBRM1 followed by western blot detection of BRD1 (Additional file 3).

Fig. 1: Co-immunoprecipitation (co-IP) coupled to mass spectrometry MS/MS analysis of BRD1-S and BRD1-L. a co-IPs were performed on cell extracts harvested from BRD1-S, BRD1-L, and HEK293T cell lines, using either anti-V5 antibody conjugated (IP-V5) or anti-HA antibody conjugated (IP-HA) beads. Cell culturing, extract preparation, and co-IP were repeated four times before preparation for mass spectrometry analysis (for further details, see "Methods"). b , c The combined spectral intensity of each individual protein was normalized to the mean intensity of all the proteins for each sample or control. In total, 319 proteins were identified from at least one sample or control in the analysis. The plots represent the normalized (Log₂ ) spectral intensities of each identified protein (black ) plotted as sample (y-axis ) against control (x-axis ). BRD1 (blue ) was identified as the most enriched protein in the BRD1-S (b ) and BRD1-L (c ) pull-downs. The plots , from the upper right to the lower left , illustrate a stepwise sorting of proteins that showed > two- or > threefold enrichment (red and orange , respectively, and otherwise black ) in the sample compared to controls. The last plots in (b ) and (c ) show the proteins that were identified as enriched in all controls, where the x-axis represents the mean of all the control experiments. These proteins were then further divided into three levels of co-IP strength taking into account the reproducibility of the co-IP within the repeated experiments (see Table 1). d BRD1 PPI network. Gray lines represent the interactions identified in this study by co-IP LC-MS/MS and black lines represent experimentally verified PPIs annotated in publicly available databases (BioGRID, CCSB Human Interaction Network Database, DIP, HPRD, IntAct, MINT, BIND, MPPI, and GNP). A high level of red in nodes indicates a high number of edges while a high level of blue indicates a low number of edges [see PDF for image]

Identification of chromatin-binding sites of BRD1 by ChIP sequencing

To identify BRD1 chromatin interactions and target genes, we conducted chromatin immunoprecipitation of epitope tagged BRD1-S and BRD1-L from the cell lines described above (Additional file 4) followed by next generation DNA sequencing (ChIP-seq) generating on average 24 million mapped reads (Additional file 5). We identified 2205 and 1722 ChIP-seq peaks for BRD1-S and BRD1-L, respectively. The identified peaks intersected more with promoter regions (defined as the region from a transcription start site and 5 kb upstream), transcribed regions (all exons, introns, and UTRs), exons (all exons), 5'UTR, and exon 1 of protein-coding genes compared to random genomic regions of the same sizes and chromosome distributions (Fig. 2a). This was not observed for 3'UTR regions and exon 2, providing evidence that the increased intersection with transcribed regions, is caused by frequent intersections with 5'UTR and/or exon 1 regions (Fig. 2a). Furthermore, BRD1-S and BRD1-L peaks were predominantly located near transcription start sites (TSSs) of genes (Fig. 2b, left panel). In order to gain more insights as to where BRD1-S and BRD1-L primarily bind at TSSs, we counted the number of sequencing reads obtained from the ChIP-seq data across a window of +/- 5 kb from TSS for all protein-coding genes in the human genome. This approach showed that the majority of sequencing reads aligned upstream (0 to -2 kb) from TSS of genes while another, smaller part of the reads aligned downstream (0 to 2 kb) from TSS of genes (Fig. 2c) suggesting that BRD1-S and BRD1-L bind mainly upstream but also downstream of TSS of genes. To further investigate whether BRD1 binds in a co-occurring manner at TSSs, we generated a heatmap of BRD1-S ChIP-seq read counts in TSSs from chromosome 1 and 10 (Fig. 2d). The analysis did not indicate that BRD1 binding was occurring twice at the same TSS. To evaluate the functional role of BRD1-S and BRD1-L, we compared all ChIP-seq regions to ChIP-seq regions of histone marks [59] and other chromatin binding proteins [60] and generated genomic binding profiles for both isoforms (Fig. 2e and Additional file 6). Both BRD1-S and BRD1-L were found to bind regions that highly overlap with histone H3K9ac genomic regions. Interestingly, the histone H3K9ac level was previously found to be reduced in erythroblasts (CD71⁺ Ter119^- ) isolated fromBrd1 ^-/- mouse embryos [7]. Furthermore, large subsets of BRD1-S and BRD1-L ChIP-seq regions were found in close proximity of or directly overlapping with RNA polymerase II binding sites (Fig. 2b, right panel and E). In view of the aforementioned results, predominant binding at RNA polymerase II binding sites, histone H3K9ac sites, promoter regions, 5'UTR, exon 1, and TSSs suggest that BRD1 primarily binds to actively transcribed regions of the genome. We observed a substantial proportion of the BRD1-S and BRD1-L peaks to be located at bi-directional/head-to-head genes that may share promoter or enhancer elements. Binding of BRD1 at these sites potentially affects transcription of both genes. Consequently, we defined BRD1 PTGs as genes having TSS within a window of +/- 5 kb from a BRD1 ChIP-seq peak. Overall, we identified sets of 1540 and 823 PTGs for BRD1-S and BRD1-L, respectively, with 251 PTGs being common to both isoforms (Fig. 2f). Among these, we validated five by ChIP followed by quantitative real-time PCR (Additional file 7 and Additional file 8).

Fig. 2: ChIP sequencing analysis of BRD1-S and BRD1-L. a The percentages of BRD1-S and BRD1-L ChIPseq peaks that intersect promoter regions of protein-coding genes (defined as 5 kb upstream from transcription start site), transcribed regions (all exons, introns, and UTRs), exons (all exons), 5'UTR, exon 1 exon 2, and 3'UTR, were identified and shown by the perforated blue and green linesin the graph. Random peak regions, comprising the same chromosome distribution and region sizes as the BRD1 ChIPseq regions, were generated and shown as Random 1 (2205 regions, red ) and Random 2 (1722 regions, gray ). b The minimum relative distance from a BRD1-S and BRD1-L ChIP-seq region to transcription start sites (left panel ) or to RNA Pol II binding sites (right panel ) was identified for all 2205 (BRD1-S) and 1722 (BRD1-L) regions and illustrated as the blue histograms . Random peak regions comprising the same chromosome distribution and region sizes as the BRD1 ChIPseq regions were generated and shown as red line histograms (Random). c In order to generate a profile of BRD1-S (blue ), BRD1-L (green ) binding to the TSS of genes, we identified and plotted the number of sequenced reads per total reads across a window of +/- 5 kb from the TSS of all protein-coding genes in the human genome. The control (Ctrl, black ) represents ChIP-seq using an antibody against the HA-epitope. d A heatmap of read counts from BRD1-S ChIP-seq across a +/- 2.5 kb window from the TSS of genes located on chromosome 1 and 10. Only regions with high read counts are shown here. e ChIPseq regions of chromatin-binding proteins and histone marks, identified in HEK293 cell lines, were obtained from the ENCODE database and from a histone H3K9ac dataset [59]. The percentages of BRD1-S and BRD1-L ChIPseq peak regions that intersected with the binding of other chromatin-binding proteins or the location of histone marks were identified and illustrated as heatmaps . f Genes comprising transcription start sites in a window of +/- 5 kb from a BRD1-S or BRD1-L ChIP-seq peak region were identified as promoter target genes (PTGs). Venn diagram illustrates the number of PTGs identified for BRD1-S and BRD1-L [see PDF for image]

Pathway analysis of the BRD1 interaction network

To further characterize the identified BRD1 interaction network, we used the Ingenuity Pathway Analysis (IPA) software to identify enrichment of specific pathways. Three datasets (BRD1-S PTGs + PPIs, BRD1-L PTGs + PPIs, and the combined BRD1 PTGs + PPIs) were analyzed according to gene categories (diseases and disorders, molecular and cellular functions, and physiological system development functions) and pathways as defined by IPA. After correction for multiple testing, we identified two, one, and 21 significant IPA categories, as well as, 14, six, and 27 significant canonical pathways among the three tested BRD1 interaction networks, respectively (categories, pathways, and genes within these can be found in Additional file 9). The most significant category in all three BRD1 interaction networks was the gene-expression category. Of neuropsychiatric relevance, we identified the abnormal morphology of neural arch and malformation of brain as significant sub-categories when analyzing the full set of BRD1 PTGs + PPIs. The most significant pathways identified for the subsets of BRD1-S and BRD1-L were IGF-1 Signaling and Cell Cycle: G1/S Checkpoint Regulation, respectively. Together the results indicate that BRD1 is a regulator of transcriptional processes and that it is involved in pathways important for brain development.

Identifying differentially expressed genes following changes in BRD1 expression levels

To identify differentially expressed genes following expression level changes in BRD1, we extracted RNA from BRD1-S and BRD1-L cells and parental HEK293T cells (Additional file 10) as well as from HEK293T cells transfected with siRNA directed against BRD1 or scrambled siRNA (Additional file 11). Microarray expression analyses covering 17,051 genes were performed on cDNA libraries constructed from the RNA samples (Additional file 2). Combining the results from both BRD1 up- and downregulation studies yielded a list of 4643 differentially expressed genes (DEGs), defined as having a probe fold change of more than 1.5, and corresponding to 27 % of the genes covered by the microarray. We subsequently examined whether the identified DEGs were enriched within BRD1 PTGs. For both BRD1-S and BRD1-L, DEGs accounted for 33 % of all PTGs, which represents a significant enrichment (P <0.001; Additional file 12). The observed enrichment supports that BRD1 binding at TSSs takes part in regulating gene expression.

Exploring the gene-regulatory potential of BRD1-S and BRD1-L isoforms

In order to explore the gene-regulatory potential of BRD1-S and BRD1-L isoforms, we performed integrative analyses of ChIP-seq and gene expression data. The number of up- and downregulated DEGs according to expression level from each study (BRD1-S or BRD1-L upregulation or BRD1 siRNA-mediated downregulation) is illustrated in Fig. 3a-d. Also, we investigated the overlap of DEGs identified after upregulating either BRD1-S or BRD1-L (Fig. 3e). We then created a window including +/- 200 kb from each BRD1 binding site, fragmented the window in 10 kb segments, and determined the number of DEGs having TSSs in these segments. The analysis was performed for DEGs identified after upregulating BRD1-S or BRD1-L or downregulating endogenous BRD1 (Fig. 3f, g, Additional file 13). Significantly more up- or downregulated genes were identified in segments where BRD1-S or BRD1-L binds in close proximity to TSSs compared to all segments across the window ( P <0.001) supporting that BRD1-S and BRD1-L are involved in regulating expression of genes when bound in close proximity to the TSSs. Within these segments, upregulation of either BRD1-S or BRD1-L resulted in more upregulated compared to downregulated genes while downregulating BRD1 resulted in more downregulated compared to upregulated genes. Furthermore, upregulation of BRD1-L proportionally generated more downregulated genes compared to BRD1-S in these segments (Fig. 3f, g). We conclude that both BRD1-S and BRD1-L primarily play a role in activation of gene expression when bound near TSSs, although they may act, and BRD1-L in particular, as a repressor for a subset of genes. Interestingly, for BRD1-L upregulated DEGs we observed the additional feature that increased expression fold changes were present when having the TSS either -50 kb upstream or 100 kb downstream the BRD1-L binding site (Fig. 3g). This suggests that BRD1-L may also participate in activating gene expression distal from its binding sites.

Fig. 3: Integration of ChIP-seq data and gene-expression data. a Examples of gene expression arrays performed. Shows (a ) upregulation of BRD1-S vs. control, (b ) upregulation of BRD1-L vs. control, (c ) siRNA downregulation of BRD1 vs. control, and (d ) two controls transfected with scrambled siRNA, control 1 vs. control 2. The red lines indicate fold change 1.5 which was used in this study to identify DEGs. The number of upregulated and downregulated DEGs is shown for each study in the upper left and lower right corner , respectively. Three siRNA knockdown experiments of BRD1 (KD1-3) were performed. More information about KD1-3 can be found in Additional file 10. eThe number of genes that overlap after upregulating BRD1-S and BRD1-L expression. After upregulating either BRD1-S or BRD1-L, the number of DEGs was identified in 10 kb segments (blue bars ) across a window of +/- 200 kb from either (f ) BRD1-S or (g ) BRD1-L binding sites, respectively. The expression fold change quartiles were obtained from segments of 10 DEGs. Gray areas indicate the upper and lower quartiles while the middle is shown as a gray line [see PDF for image]

In order to investigate BRD1 gene-regulatory potential while considering all genes on the array (and not only DEGs), we segmented and calculated the average fold change of all probes within windows +/- 10 kb and +/- 200 kb from BRD1 binding sites similarly as in the previous analysis. Overexpressing BRD1-S primarily increased the overall expression fold change of genes with BRD1-S binding approximately +/- 1 kb from their TSSs while downregulating BRD1 reversibly downregulated genes within the same range (Additional file 13). Overexpressing BRD1-L resulted in increased expression of genes having TSSs -50 kb upstream and +100 kb downstream of BRD1-L binding similar to what we observed in the previous analysis of the DEGs; however, a clear reverse pattern was not observed after downregulating BRD1 (Additional file 13).

The integrative analysis of ChIP-seq and gene expression data collectively suggests that BRD1 can activate and more rarely, repress, gene expression by binding to TSSs, as well as suggesting that specifically the BRD1-L isoform plays a role in activating gene-expression distal from its chromatin binding sites.

Pathway analysis of DEGs identified as a consequence of upregulating BRD1-S and BRD1-L

To elucidate potential functional consequences of the BRD1-mediated changes in gene expression and identify potential differences between BRD1-S and BRD1-L isoforms, the DEGs identified by upregulating either BRD1-S or BRD1-L were analyzed by IPA software to identify enrichment of specific canonical pathways. DEGs were defined as previously stated using a 1.5 threshold on the probe fold change. The top ten most significant canonical pathways can be found in Additional file 14. The genes in all the pathways were mainly upregulated. The most significant canonical pathway identified when upregulating BRD1-S or BRD1-L was Reelin Signaling in Neurons and IL-8 Signaling, respectively. Of notice, the IGF1 Signaling pathway which was identified to be the most significant pathway among the BRD1-S interaction network was also identified among the ten most significant pathways of the BRD1-S DEGs further implicating BRD1-S with IGF1 signaling. Among the identified pathways, the Reelin Signaling in Neurons pathway is of a particular interest to brain function and development [61].

Correlation of the BRD1 interaction network with spatiotemporal BRD1 expression profiles in the human brain

In order to explore the role of BRD1 in the human brain, we obtained temporal RNA-seq profiles and expression array profiles from the Brainspan atlas of the developing human brain [43] and the Human Brain Transcriptome database [62], respectively. A high BRD1 expression was observed in the fetal stages (until approximately day 180) for all brain regions whereas a high BRD1 expression in the cerebellar cortex was observed throughout the timeline of the datasets (age >40 years; Additional file 15). To investigate whether the identified BRD1 interaction network is likely to operate in the context of the human brain, we used data from the Brainspan RNA-seq database to correlate the expression of 21,906 genes with BRD1 expression across 13 developmental stages and 26 different brain regions in the human brain. Correlation values (Fig. 4a) for all genes ranged from 1 (complete correlation with BRD1 expression) to -1 (complete anti-correlation with BRD1 expression). In total, 1241 and 672 interaction partners (PTGs + PPIs) for BRD1-S and BRD1-L, respectively, were identified in the Brainspan database (Fig. 4b). We observed significantly more correlation for both of the BRD1 isoforms than expected by chance (Fig. 4c), indicating that the BRD1 interactions are likely to occur also in the human brain. The PPIs identified for BRD1-S and BRD1-L were also more correlated with BRD1 expression (P <1 x 10^-10 and P = 0.0039, respectively); in particular,SUV420H1 , DNMT1 , PBRM1 , and CHD4 were highly correlated with BRD1 (Fig. 4d).

Fig. 4: Co-expression analysis of the BRD1 interaction network in human brain. a A total of 21,906 gene-correlation values ranging from 1 to -1 were obtained across 13 developmental stages and 26 different brain regions from the Brainspan RNA-seq database. Correlation values for all genes were illustrated as a blue line and gray and white bins (#1-4) define sections (1 to 0.5, 0.5 to 0, 0 to -0.5, and -0.5 to -1) where the number of genes in each bin was shown at the top. b From the identified BRD1-S and BRD1-L PTGs + PPIs, 1241 and 672 were identified in the database, respectively. The correlation values for each gene distribute into bins #1-4 as illustrated. c The distributions of the correlation values of all genes in the Brainspan database and the correlation values of BRD1-S and BRD1-L PTGs + PPIs (BRD1-S and BRD1-L interaction networks) were summarized as relative to the number of genes and illustrated in a histogram (a.u. = arbitrary units). Mann-Whitney tests were performed for correlation values from either of the BRD1-S or BRD1-L interaction networks and all genes identified. Both tests resulted in the rejection of the null hypothesis of identical medians (P < 1 x 10^-10 ). d Correlation values for the BRD1 PPIs represent the level of co-expression with BRD1 in the developing human brain [see PDF for image]

Disease risk enrichment analysis of the BRD1 interaction network

Using available GWAS datasets obtained from the Psychiatric Genomics Consortium (PGC) we investigated whether the BRD1 interaction network is enriched for mental disorder risk applying the MAGMA program [38]. As the evidence for genetic association to the BRD1 locus is particularly strong in schizophrenia, we tested for enrichment in this disorder as our main analysis. Out of the entire BRD1 network, 1853 genes remained after filtering SNPs and removing genes in the MHC region (see Methods), 468 of these had gene-wide P values <0.05, and 44 of these were significant after adjusting for the number of genes tested (Additional file 16). In this analysis the entire BRD1 network showed significant enrichment (Fig. 5).

Fig. 5: Disease risk gene analysis of the BRD1 network. Enrichment for genetic disease risk was investigated for seven BRD1 sub-networks across 12 GWASs comprising eight brain disorders and four disorders that are not considered brain disorders. The bars indicate BRD1 sub-networks in the order: all BRD1-S and BRD1-L PTGs and PPIs (black ), BRD1-S PTGs and PPIs, BRD1-S PTGs, BRD1-S PPIs, BRD1-L PTGs and PPIs, BRD1-L PTGs, and BRD1-L PPIs. GWAS from the left are: schizophrenia ( scz ) [4, 32], attention deficit hyperactivity disorder (adhd ) [27, 32], major depressive disorder (mdd ) [28, 32], autism spectrum disorder (asd ) [32], bipolar disorder (bip ) [30, 32], psychiatric cross-disorders (cross ) [25, 32], anorexia nervosa (an ) [23, 32], Alzheimer's disease (alz ) [24, 37], Crohn's disease (crohns ) [36], rheumatoid arthritis (ra ) [29, 33], coronary artery disease (cad ) [31, 34], and type 2 diabetes (t2d ) [26, 35]. The red line indicates P = 0.05, before correction for the number of GWAS tested. The number of significant loci as well as the number of subjects in each study is noted under each study [see PDF for image]

In order to explore whether the interaction network might also be enriched for susceptibility to other mental disorders, we tested 12 other GWAS datasets, comprising eight psychiatric disorders and four disorders not considered related to brain function. In addition, to investigate whether the risk enrichment is dependent on the specific isoform, we stratified for BRD1-S and BRD1-L sub-networks. When adjusting for the number of GWASs tested, significant enrichment was observed solely for the entire BRD1 network in schizophrenia (P _adjusted = 0.03, Fig. 5). Interestingly, substantial differences in enrichment patterns between the two BRD1 isoforms were observed in several disorders. The BRD1-S sub-networks comprising PTGs showed enrichment for schizophrenia and anorexia nervosa while the PPI sub-network showed evidence of enrichment for attention deficit hyperactivity disorder (ADHD) and major depressive disorder (MDD) risk before adjusting for the number of tests (P <0.05). These results suggest that the BRD1-S network may play a larger role in mental disorder risk compared to the BRD1-L network.

Finally, we explored whether the network was enriched for de novo mutations and rare disruptive variants observed in autism and schizophrenia [40-42]. The analysis did not identify significant enrichment. However, in particular among the PPI network genes, we did observe de novo mutations in autism probands that were not identified in healthy siblings, including one missense mutation inBRD1 , one missense mutation in CHD4 , and four disruptive and two missense mutations inSUV420H1 [40, 41], as well as rare disruptive mutations seen only in schizophrenia cases and not in controls (in BRD1 , ING5 , YWHAZ , NAP1L4 , MYST2 , and HIST1H3A ) [42].

In summary, we identify the BRD1 interaction network to be enriched for schizophrenia risk thereby providing supporting evidence that BRD1 plays a role in the etiology of schizophrenia. Additionally, our analyses suggest that the BRD1-S interaction network is more enriched for mental disorder risk, including schizophrenia risk, compared to the interaction network of BRD1-L.

Exploring the BRD1 spatiotemporal interaction networks in the human brain

Taking into consideration that the BRD1 network is enriched with schizophrenia risk and that this enrichment seems to be more pronounced in the BRD1-S network, we asked whether these networks are more likely to be interacting in specific regions and developmental intervals in the human brain. In order to evaluate co-expression with BRD1 in the human brain, we applied an approach similar to the one described here [63]. We used RNA-seq data from Brainspan to build 32 spatiotemporal intervals comprising eight temporal intervals (P1-P8) and four brain regions (R1-R4) (Fig. 6). We then calculated the fraction of genes that were expressed with Spearman's correlation coefficient >0.5 compared to BRD1 for each of the 32 intervals. Previously we showed that the majority of genes in the BRD1 network are positively correlated with BRD1 expression across all samples in Brainspan, including the BRD1 PPI networks. To gain insight into which spatiotemporal intervals that may be driving this positive correlation, we investigated the fraction of BRD1 co-expressed genes in the BRD1-S and BRD1-L PPI networks and compared them to the background fraction across all genes in the Brainspan dataset. A higher fraction of co-expressed genes was observed for both networks during pcw 16-26 (P2, second trimester) in the prefrontal cortex regions (R2) and the hippocampus, striatum, and amygdala regions (R3) (Fig. 6a, top panel and Additional file 17). The BRD1-S PPI network was found to have a higher fraction of co-expressed genes across several spatiotemporal intervals including P3R3, P5R1-R3, P6R1-R4, P7R4, and P8R1-3.

Fig. 6: Identifying spatiotemporal networks of BRD1 in the human brain. a The fraction of co-expressed genes (correlation coefficient >0.5) was calculated for each spatiotemporal interval and each BRD1 sub-network and compared to the background of genes co-expressed with BRD1 in the entire dataset. The color-coding denotes the fold change of the fraction of co-expressed genes compared to the background. b The fraction of genes that co-expressed with BRD1 across each of the 32 spatiotemporal intervals was calculated for all genes (gray lines), genes in the BRD1 PTGs + PPI network (orange lines), and genes in BRD1 schizophrenia network of 468 genes (red lines). A one-sided binominal test was performed for each BRD1 sub-network and spatiotemporal interval compared to the expected fraction in the background where the asterisk (*) denote P < 0.05 after adjusting for the number of tests. Temporal intervals (P) were grouped as follows: P1 included pcw 8-13 (first trimester), P2 included 16-26 pcw (second trimester), P3 included 35-37 pcw (third trimester), P4 included 4-10 months, P5 included 1-4 years, P6 included 8-13 years, P7 included 15-19 years, and P8 included 21-40 years. Brain regions (R) were grouped as follows: R1 included the posterior inferior parietal cortex, primary auditory cortex, primary visual cortex, superior temporal cortex, and inferior temporal cortex; R2 included the primary somatosensory cortex, primary motor cortex, orbital prefrontal cortex, dorsolateral prefrontal cortex, medial prefrontal cortex, and ventrolateral prefrontal cortex; R3 included the striatum, hippocampus, and amygdala; and R4 included the mediodorsal nucleus of the thalamus, and cerebella cortex [see PDF for image]

We then expanded the analysis to include the entire BRD1 interaction network (BRD1-S and BRD1-L PTGs and PPIs) and the BRD1 schizophrenia sub-network consisting of the 468 genes that achieved a gene-wide P value <0.05. The highest fraction of genes co-expressed with BRD1 consistent across all BRD1 networks was observed for pwc 16-26 (P2) in the hippocampus, striatum, and amygdala regions (R3) (Fig. 6a). All BRD1 PTG + PPI networks showed a significantly higher fraction of genes co-expressed with BRD1 in this interval (Fig. 6b and Additional file 17). We noticed that the BRD1 schizophrenia sub-network, across several spatiotemporal intervals, was more enriched with genes that co-expressed with BRD1 compared to other networks that included PTGs (Fig. 6a, bottom panel). These intervals were largely co-occurring with those of the many times smaller BRD1-S PPI network. The BRD1 schizophrenia sub-network comprised a significantly higher fraction of genes that co-expressed with BRD1 in 12 of the total 32 spatiotemporal intervals (Fig. 6b).

Discussion

In this study we have obtained results providing novel insights into the molecular interactions of BRD1 and its relation to mental disorders.

The molecular function of BRD1

Our results provide confirming evidence that BRD1 is part of a histone acetyltransferase complex comprising ING4, ING5, MEAF6, and MYST2 [7, 49] along with supporting evidence for BRD1 and histone H3 interaction [50].

Interestingly, we identify several novel protein interactions that are important for chromatin modulation and gene regulation, including PBRM1 which is part of the SWI/SNF chromatin remodeling complex PBAF [64] and the 14-3-3 proteins YWHAE, G, H, and Z. Among several known functions, 14-3-3 proteins facilitate the recruitment of transcription factors to chromatin particularly in conjunction with phosphorylated histone H3S10 and H3S28 [51, 65, 66], e.g. the PBAF complex is recruited to chromatin via the 14-3-3 proteins and phosphorylated histone H3S10 [64]. Our results identify both BRD1 isoforms to bind 14-3-3 proteins as well as PBRM1 indicating that BRD1 plays a role in chromatin remodeling.

Although the identified PPI networks of BRD1-S and BRD1-L seem to comprise a common set of core proteins, we identify several isoform-specific PPIs that could indicate a significant functional difference between the two isoforms. For instance, the histone methyltransferase SUV420H1 was identified only for the BRD1-S isoform and QSER1, a relatively uncharacterized protein, was only identified for BRD1-L. In addition, BRD1-S and BRD1-L had a relatively limited overlap of PTGs as well as DEGs suggesting that these two isoforms, although very similar in amino acid sequence, target and regulate different gene sets. By integrating ChIP-seq and expression data, we were able to show that BRD1-S primarily regulate gene expression in proximity to TSSs while BRD1-L, in addition to regulating gene expression in proximity to TSSs, has the capacity to regulate gene expression further up- and downstream from its binding sites. Based on our observations, it is tempting to speculate that BRD1-L is involved in mediating chromatin loop status thereby regulating the expression of distal genes similar to the mechanisms identified for SATB1 [67]. Conversely, BRD1-S may primarily regulate gene expression by binding in close proximity to the TSSs of genes and by mediating acetylation of histones, in particular acetylation of histone 3 lysine 14, and chromatin remodeling through its interactions with MYST2 [7] and PBRM1, respectively.

The BRD1 interaction networks and their relation to mental disorders

Several studies have implicated 14-3-3 genes with schizophrenia and bipolar disorder [52-55, 68], YWHAE and YWHAZ interact with disrupted in schizophrenia 1 (DISC1) [69], and Ywhae ^+/- andYwhaz ^-/- mice have been shown to display neurodevelopmental and schizophrenia-associated phenotypes [70, 71]. Furthermore, the PBRM1 gene has been associated with susceptibility to schizophrenia and bipolar disorder [4, 30, 56-58, 72, 73]. Our finding that these proteins all interact with BRD1 provides a novel line of evidence implicating the network of BRD1, PBRM1, and 14-3-3 proteins with schizophrenia and bipolar disorder.

In addition to the similarities and differences of the molecular and cellular functions of BRD1-S and BRD1-L, we explored their interaction networks for enrichment of disease risk. We observed that the entire BRD1 network is enriched for schizophrenia risk and we found that the BRD1-S network was generally more enriched with schizophrenia risk and perhaps other mental disorders compared to BRD1-L. Moreover, the entire BRD1 network was more enriched with schizophrenia compared to the BRD1-S network, suggesting an enhanced effect of combining the networks from both isoforms. Our results support a recently published large-scale study of autism risk genes which demonstrates the importance of considering splice isoforms when exploring disease networks [74]. It should be noted that the identified BRD1 networks were not identified in brain-derived cells. However, our integrative analysis using spatiotemporal expression data from the human brain suggest that the identified interactions may also be occurring and operating in the human brain.

Although the BRD1 interaction network was not identified to be enriched with bipolar disorder risk genes, it is reasonable to speculate that this difference is based on power constraints in the original GWAS rather than on a less prominent role of BRD1 in the etiology of bipolar disorder. Future studies in larger bipolar disorder GWAS datasets are warranted to confirm this hypothesis. In general, the smaller GWAS sets may be underpowered which should be taken into consideration when interpreting the results.

Further integrative analyses using spatiotemporal transcriptome data from the human brain suggested that both BRD1-S and BRD1-L interaction networks play a role in brain function and schizophrenia, especially at mid-fetal stages (pcw 16-26) in the hippocampus, amygdala, and striatum. However, the BRD1-S networks seem also to play a role in brain function and schizophrenia later in life, particularly during childhood and early adulthood. These observations are supportive of the enrichment analysis showing the highest enrichment of schizophrenia risk when analyzing the entire BRD1 network as well as more enrichment of schizophrenia risk when analyzing the BRD1-S network compared to the BRD1-L network.

Conclusions

Here we expand the molecular and functional characterization of BRD1 and provide evidence that BRD1 acts as a regulatory hub in a comprehensive schizophrenia risk network and possibly risk networks for other mental disorders as well, thereby supporting previous association studies implicating the BRD1 gene with schizophrenia and bipolar disorder [1-6]. Furthermore, we identify spatiotemporal intervals in the human brain where BRD1 sub-networks are likely to play a role in brain function and schizophrenia. These results encourage further research of BRD1, for example using in vivo models.

Abbreviations: BRD1: Bromodomain containing 1; DEG: Differentially expressed gene; IPA: Ingenuity pathway analysis; PBAF: Polybromo-associated BRG1 associated factor; PPI: Protein-protein interaction; PTG: Promotor target gene; SWI/SNF: SWItch/Sucrose Non-Fermentable; TSS: Transcription start site; UTR: Untranslated region

Competing interests: Drs Mors and Barglum are co-inventors on a patent application submitted by Aarhus University entitled "Method for diagnosis and treatment of a mental disease" (EP2287340) that includes claims relating to BRD1 among other genes. The other authors declare no conflict of interest.

Authors' contributions: TF designed and performed the laboratory and bioinformatics research, analyzed and interpreted data, prepared the first draft, and wrote the manuscript. JHC, ALN, and ADB designed the research, interpreted data, and wrote the manuscript. JP(1), MEH, and MM performed bioinformatic analyses, analyzed GWAS and exome sequencing data, and wrote the manuscript. HSO performed bioinformatics analyses and wrote the manuscript. JB and TC performed laboratory research and reviewed the manuscript. JP(2) and MB performed laboratory research, analyzed data, and reviewed the manuscript. JG, DD, MN, NG, NT, and OM acquired data, interpreted data, and reviewed the manuscript. All authors read and approved the final manuscript.

Acknowledgements: We thank John Strouboulis, Frank Grosveld, and colleagues for sending us a copy of the V5 ChIP protocol. Financial support: Grants from The Lundbeck Foundation, The Danish Strategic Research Council, The Danish Council for Independent Research | Medical Sciences, John and Birthe Meyer Foundation, The Faculty of Health Sciences at Aarhus University, and The Novo Nordisk Foundation. We are grateful to The Encyclopedia of DNA Elements Consortium (ENCODE), all working groups in the Psychiatric Genomics Consortium (PGC), the DIAbetes Genetics Replication And Meta-analysis Consortium (DIAGRAM), the Coronary artery disease and myocardial infarction (CARDIOGRAM), the International Inflammatory Bowel Disease Genetics Consortium (IIBDGC), and the International Genomics of Alzheimer's Project (IGAP) for providing summary results data for these analyses. The investigators within IGAP contributed to the design and implementation of IGAP and/or provided data but did not participate in analysis or writing of this report. IGAP was made possible by the generous participation of the control subjects, the patients, and their families. The i-Select chips were funded by the French National Foundation on Alzheimer's disease and related disorders. EADI was supported by the LABEX (laboratory of excellence program investment for the future) DISTALZ grant, Inserm, Institut Pasteur de Lille, Universitñ de Lille 2, and the Lille University Hospital. GERAD was supported by the Medical Research Council (Grant no. 503480), Alzheimer's Research UK (Grant no. 503176), the Wellcome Trust (Grant no. 082604/2/07/Z) and German Federal Ministry of Education and Research (BMBF): Competence Network Dementia (CND) grant nos. 01GI0102, 01GI0711, 01GI0420. CHARGE was partly supported by the NIH/NIA grant R01 AG033193 and the NIA AG081220 and AGES contract N01-AG-12100, the NHLBI grant R01 HL105756, the Icelandic Heart Association, and the Erasmus Medical Center and Erasmus University. ADGC was supported by the NIH/NIA grants: U01 AG032984, U24 AG021886, U01 AG016976, and the Alzheimer's Association grant ADGC-10-196728.

Correspondence: Anders D. Barglum:

Author details: ¹ Department of Biomedicine, Aarhus University, Building 1242, Bartholins Allñ 6, 8000, Aarhus C, Denmark. ² iPSYCH, The Lundbeck Foundation Initiative for Integrative Psychiatric Research, 8000 Aarhus C, Denmark. ³ iSEQ, Centre for Integrative Sequencing, Aarhus University, 8000, Aarhus C, Denmark. ⁴ Research Unit for Molecular Medicine, Aarhus University Hospital, 8200, Skejby, Denmark. ⁵ Wilhelm Johannsen Centre for Functional Genome Research, Department of Cellular and Molecular Medicine, University of Copenhagen, 2200, Copenhagen N, Denmark. ⁶Bioinformatics Research Centre (BiRC, Aarhus University, 8000, Aarhus C, Denmark. ⁷ Research Department P, Aarhus University Hospital, 8240, Risskov, Denmark.

Article history: Received 22 January 2015 Accepted 15 April 2016 Published online 3 May 2016

References

1. Nyegaard M, Severinsen JE, Als TD, Hedemand A, Straarup S, Nordentoft M. Support of association between BRD1 and both schizophrenia and bipolar affective disorder. Am J Med Genet B Neuropsychiatr Genet.. 2010;153B:582.-91.
2. Purcell SM, Wray NR, Stone JL, Visscher PM, O'Donovan MC, Sullivan PF. Common polygenic variation contributes to risk of schizophrenia and bipolar disorder. Nature.. 2009;460:748.-52.
3. Severinsen JE, Bjarkam CR, Kiaer-Larsen S, Olsen IM, Nielsen MM, Blechingberg J. Evidence implicating BRD1 with brain development and susceptibility to both schizophrenia and bipolar affective disorder. Mol Psychiatry.. 2006;11:1126.-38 doi: 10.1038/sj.mp.4001885.
4. Biological insights from 108 schizophrenia-associated genetic loci. Nature. 2014;511:421.-7 doi: 10.1038/nature13595.
5. Andreassen OA, Thompson WK, Dale AM. Boosting the power of schizophrenia genetics by leveraging new statistical tools. Schizophr Bull.. 2014;40:13.-7 doi: 10.1093/schbul/sbt168.
6. Aberg KA, Liu Y, Bukszar J, McClay JL, Khachane AN, Andreassen OA. A comprehensive family-based replication study of schizophrenia genes. JAMA Psychiatry.. 2013;70:1.-9 doi: 10.1001/jamapsychiatry.2013.288.
7. Mishima Y, Miyagi S, Saraya A, Negishi M, Endoh M, Endo TA. The Hbo1-Brd1/Brpf2 complex is responsible for global acetylation of H3K14 and required for fetal liver erythropoiesis. Blood.. 2011;118:2443.-53 doi: 10.1182/blood-2011-01-331892.
8. Kueh AJ, Dixon MP, Voss AK, Thomas T. HBO1 is required for H3K14 acetylation and normal transcriptional activity during embryonic development. Mol Cell Biol.. 2011;31:845.-60 doi: 10.1128/MCB.00159-10.
9. Bjarkam CR, Corydon TJ, Olsen IM, Pallesen J, Nyegaard M, Fryland T. Further immunohistochemical characterization of BRD1 a new susceptibility gene for schizophrenia and bipolar affective disorder. Brain Struct Funct.. 2009;214:37.-47 doi: 10.1007/s00429-009-0219-3.
10. Christensen JH, Elfving B, Muller HK, Fryland T, Nyegaard M, Corydon TJ. The schizophrenia and bipolar disorder associated BRD1 gene is regulated upon chronic restraint stress. Eur Neuropsychopharmacol.. 2012;22:651.-6 doi: 10.1016/j.euroneuro.2012.01.005.
11. Fryland T, Elfving B, Christensen JH, Mors O, Wegener G, Borglum AD. Electroconvulsive seizures regulates the Brd1 gene in the frontal cortex and hippocampus of the adult rat. Neurosci Lett.. 2012;516:110.-3 doi: 10.1016/j.neulet.2012.03.069.
12. Klenova E, Chernukhin I, Inoue T, Shamsuddin S, Norton J. Immunoprecipitation techniques for the analysis of transcription factor complexes. Methods.. 2002;26:254.-9 doi: 10.1016/S1046-2023(02)00029-4.
13. Hansen J, Palmfeldt J, Vang S, Corydon TJ, Gregersen N, Bross P. Quantitative proteomics reveals cellular targets of celastrol. PLoS One.. 2011;6:e26634. doi: 10.1371/journal.pone.0026634.
14. Cox J, Mann M. MaxQuant enables high peptide identification rates, individualized p.p.b.-range mass accuracies and proteome-wide protein quantification. Nat Biotechnol. 2008;26:1367.-72 doi: 10.1038/nbt.1511.
15. Kolodziej KE, Pourfarzad F, Boer E, Krpic S, Grosveld F, Strouboulis J. Optimal use of tandem biotin and V5 tags in ChIP assays. BMC Mol Biol.. 2009;10:6. doi: 10.1186/1471-2199-10-6.
16.  Goecks J, Nekrutenko A, Taylor J. Galaxy: a comprehensive approach for supporting accessible, reproducible, and transparent computational research in the life sciences. Genome Biol.. 2010;11:R86. doi: 10.1186/gb-2010-11-8-r86.
17. Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol.. 2009;10:R25. doi: 10.1186/gb-2009-10-3-r25.
18. Zhang Y, Liu T, Meyer CA, Eeckhoute J, Johnson DS, Bernstein BE. Model-based analysis of ChIP-Seq (MACS). Genome Biol.. 2008;9:R137. doi: 10.1186/gb-2008-9-9-r137.
19. Andersen CL, Christensen LL, Thorsen K, Schepeler T, Sorensen FB, Verspaget HW. Dysregulation of the transcription factors SOX4, CBFB and SMARCC1 correlates with outcome of colorectal cancer. Br J Cancer.. 2009;100:511.-23 doi: 10.1038/sj.bjc.6604884.
20. Mansilla F, Birkenkamp-Demtroder K, Kruhoffer M, Sorensen FB, Andersen CL, Laiho P. Differential expression of DHHC9 in microsatellite stable and instable human colorectal cancer subgroups. Br J Cancer.. 2007;96:1896.-903 doi: 10.1038/sj.bjc.6603818.
21. Rozen S, Skaletsky H. Primer3 on the WWW for general users and for biologist programmers. Methods Mol Biol.. 2000;132:365.-86.
22. Identification of risk loci with shared effects on five major psychiatric disorders: a genome-wide analysis. Lancet. 2013;381:1371.-9 doi: 10.1016/S0140-6736(12)62129-1.
23. Boraska V, Franklin CS, Floyd JA, Thornton LM, Huckins LM, Southam L. A genome-wide association study of anorexia nervosa. Mol Psychiatry.. 2014;19:1085.-94 doi: 10.1038/mp.2013.187.
24. Lambert JC, Ibrahim-Verbaas CA, Harold D, Naj AC, Sims R, Bellenguez C. Meta-analysis of 74,046 individuals identifies 11 new susceptibility loci for Alzheimer's disease. Nat Genet.. 2013;45:1452.-8 doi: 10.1038/ng.2802.
25. Lee SH, Ripke S, Neale BM, Faraone SV, Purcell SM, Perlis RH. Genetic relationship between five psychiatric disorders estimated from genome-wide SNPs. Nat Genet.. 2013;45:984.-94 doi: 10.1038/ng.2711.
26. Morris AP, Voight BF, Teslovich TM, Ferreira T, Segre AV, Steinthorsdottir V. Large-scale association analysis provides insights into the genetic architecture and pathophysiology of type 2 diabetes. Nat Genet.. 2012;44:981.-90 doi: 10.1038/ng.2383.
27. Neale BM, Medland SE, Ripke S, Asherson P, Franke B, Lesch KP. Meta-analysis of genome-wide association studies of attention-deficit/hyperactivity disorder. J Am Acad Child Adolesc Psychiatry.. 2010;49:884.-97 doi: 10.1016/j.jaac.2010.06.008.
28. Ripke S, Wray NR, Lewis CM, Hamilton SP, Weissman MM, Breen G. A mega-analysis of genome-wide association studies for major depressive disorder. Mol Psychiatry.. 2013;18:497.-511 doi: 10.1038/mp.2012.21.
29. Stahl EA, Raychaudhuri S, Remmers EF, Xie G, Eyre S, Thomson BP. Genome-wide association study meta-analysis identifies seven new rheumatoid arthritis risk loci. Nat Genet.. 2010;42:508.-14 doi: 10.1038/ng.582.
30. Large-scale genome-wide association analysis of bipolar disorder identifies a new susceptibility locus near ODZ4. Nat Genet. 2011;43:977.-83 doi: 10.1038/ng.943.
31. Wain LV, Verwoert GC, O'Reilly PF, Shi G, Johnson T, Johnson AD. Genome-wide association study identifies six new loci influencing pulse pressure and mean arterial pressure. Nat Genet.. 2011;43:1005.-11 doi: 10.1038/ng.922.
32. Psychiatric Genomics Consortium (PGC). http://www.med.unc.edu/pgc.
33. Rheumatoid arthritis. http://www.broadinstitute.org/ftp/pub/rheumatoid_arthritis/Stahl_etal_2010NG/.
34. CARDIoGRAMplusC4D (Coronary ARtery DIsease Genome wide Replication and Meta-analysis (CARDIoGRAM) plus The Coronary Artery Disease (C4D) Genetics) consortium. http://www.cardiogramplusc4d.org/.
35. DIAbetes Genetics Replication and Meta-analysis Consortium. diagram-consortium.org
36. International Inflammatory Bowel Disease Genetics Consortium (IIBDGC). http://www.ibdgenetics.org/.
37. International Genomics of Alzheimer's Project. http://web.pasteurlille.fr/en/recherche/u744/igap/igap_download.php
38. Leeuw CA, Mooij JM, Heskes T, Posthuma D. MAGMA: generalized gene-set analysis of GWAS data. PLoS Comput Biol.. 2015;11:e1004219. doi: 10.1371/journal.pcbi.1004219.
39. Abecasis GR, Auton A, Brooks LD, DePristo MA, Durbin RM, Handsaker RE. An integrated map of genetic variation from 1,092 human genomes. Nature.. 2012;491:56.-65 doi: 10.1038/nature11632.
40. Iossifov I, O'Roak BJ, Sanders SJ, Ronemus M, Krumm N, Levy D. The contribution of de novo coding mutations to autism spectrum disorder. Nature.. 2014;515:216.-21 doi: 10.1038/nature13908.
41. De RS, He X, Goldberg AP, Poultney CS, Samocha K, Cicek AE. Synaptic, transcriptional and chromatin genes disrupted in autism. Nature.. 2014;515:209.-15 doi: 10.1038/nature13772.
42. Purcell SM, Moran JL, Fromer M, Ruderfer D, Solovieff N, Roussos P. A polygenic burden of rare disruptive mutations in schizophrenia. Nature.. 2014;506:185.-90 doi: 10.1038/nature12975.
43. BrainSpan: Atlas of the Developing Human Brain [Internet]. Funded by ARRA Awards 1RC2MH089921-01, 1RC2MH090047-01, and 1RC2MH089929-01. [c] 2011. Available from: http://developinghumanbrain.org.
44. Willsey AJ, Sanders SJ, Li M, Dong S, Tebbenkamp AT, Muhle RA. Coexpression networks implicate human midfetal deep cortical projection neurons in the pathogenesis of autism. Cell.. 2013;155:997.-1007 doi: 10.1016/j.cell.2013.10.020.
45. Edgar R, Domrachev M, Lash AE. Gene Expression Omnibus: NCBI gene expression and hybridization array data repository. Nucleic Acids Res.. 2002;30:207.-10 doi: 10.1093/nar/30.1.207.
46. The International Molecular Exchange Consortium, IMEx. http://www.imexconsortium.org.
47. Orchard S, Kerrien S, Abbani S, Aranda B, Bhate J, Bidwell S. Protein interaction data curation: the International Molecular Exchange (IMEx) consortium. Nat Methods.. 2012;9:345.-50 doi: 10.1038/nmeth.1931.
48. Orchard S, Ammari M, Aranda B, Breuza L, Briganti L, Broackes-Carter F. The MIntAct project--IntAct as a common curation platform for 11 molecular interaction databases. Nucleic Acids Res.. 2014;42:D358.-63 doi: 10.1093/nar/gkt1115.
49. Doyon Y, Cayrou C, Ullah M, Landry AJ, Cote V, Selleck W. ING tumor suppressor proteins are critical regulators of chromatin acetylation required for genome expression and perpetuation. Mol Cell.. 2006;21:51.-64 doi: 10.1016/j.molcel.2005.12.007.
50. Qin S, Jin L, Zhang J, Liu L, Ji P, Wu M. Recognition of unmodified histone H3 by the first PHD finger of bromodomain-PHD finger protein 2 provides insights into the regulation of histone acetyltransferases monocytic leukemic zinc-finger protein (MOZ) and MOZ-related factor (MORF). J Biol Chem.. 2011;286:36944.-55 doi: 10.1074/jbc.M111.244400.
51. Healy S, Khan DH, Davie JR. Gene expression regulation through 14-3-3 interactions with histones and HDACs. Discov Med.. 2011;11:349.-58.
52. Wong AH, Macciardi F, Klempan T, Kawczynski W, Barr CL, Lakatoo S. Identification of candidate genes for psychosis in rat models, and possible association between schizophrenia and the 14-3-3eta gene. Mol Psychiatry.. 2003;8:156.-66 doi: 10.1038/sj.mp.4001237.
53. Wong AH, Likhodi O, Trakalo J, Yusuf M, Sinha A, Pato CN. Genetic and post-mortem mRNA analysis of the 14-3-3 genes that encode phosphoserine/threonine-binding regulatory proteins in schizophrenia and bipolar disorder. Schizophr Res.. 2005;78:137.-46 doi: 10.1016/j.schres.2005.06.009.
54. Grover D, Verma R, Goes FS, Mahon PL, Gershon ES, McMahon FJ. Family-based association of YWHAH in psychotic bipolar disorder. Am J Med Genet B Neuropsychiatr Genet.. 2009;150B:977.-83 doi: 10.1002/ajmg.b.30927.
55. Pers TH, Hansen NT, Lage K, Koefoed P, Dworzynski P, Miller ML. Meta-analysis of heterogeneous data sources for genome-scale identification of risk genes in complex phenotypes. Genet Epidemiol.. 2011;35:318.-32 doi: 10.1002/gepi.20580.
56. Ripke S, O'Dushlaine C, Chambert K, Moran JL, Kahler AK, Akterin S. Genome-wide association analysis identifies 13 new risk loci for schizophrenia. Nat Genet.. 2013;45:1150.-9 doi: 10.1038/ng.2742.
57. McMahon FJ, Akula N, Schulze TG, Muglia P, Tozzi F, Detera-Wadleigh SD. Meta-analysis of genome-wide association data identifies a risk locus for major mood disorders on 3p21.1. Nat Genet. 2010;42:128.-31 doi: 10.1038/ng.523.
58. Williams HJ, Craddock N, Russo G, Hamshere ML, Moskvina V, Dwyer S. Most genome-wide significant susceptibility loci for schizophrenia and bipolar disorder reported to date cross-traditional diagnostic boundaries. Hum Mol Genet.. 2011;20:387.-91 doi: 10.1093/hmg/ddq471.
59. Luo Y, Blechingberg J, Fernandes AM, Li S, Fryland T, Borglum AD. EWS and FUS bind a subset of transcribed genes encoding proteins enriched in RNA regulatory functions. BMC Genomics.. 2015;16:929. doi: 10.1186/s12864-015-2125-9.
60. An integrated encyclopedia of DNA elements in the human genome. Nature. 2012;489:57.-74 doi: 10.1038/nature11247.
61. Rice DS, Curran T. Role of the reelin signaling pathway in central nervous system development. Annu Rev Neurosci.. 2001;24:1005.-39 doi: 10.1146/annurev.neuro.24.1.1005.
62. Kang HJ, Kawasawa YI, Cheng F, Zhu Y, Xu X, Li M. Spatio-temporal transcriptome of the human brain. Nature.. 2011;478:483.-9 doi: 10.1038/nature10523.
63. Lin GN, Corominas R, Lemmens I, Yang X, Tavernier J, Hill DE. Spatiotemporal 16p11.2 protein network implicates cortical late mid-fetal brain development and KCTD13-Cul3-RhoA pathway in psychiatric diseases. Neuron. 2015;85:742.-54 doi: 10.1016/j.neuron.2015.01.010.
64. Drobic B, Perez-Cadahia B, Yu J, Kung SK, Davie JR. Promoter chromatin remodeling of immediate-early genes is mediated through H3 phosphorylation at either serine 28 or 10 by the MSK1 multi-protein complex. Nucleic Acids Res.. 2010;38:3196.-208 doi: 10.1093/nar/gkq030.
65. Galliher-Beckley AJ, Williams JG, Cidlowski JA. Ligand-independent phosphorylation of the glucocorticoid receptor integrates cellular stress pathways with nuclear receptor signaling. Mol Cell Biol.. 2011;31:4663.-75 doi: 10.1128/MCB.05866-11.
66. Kellner WA, Ramos E, Van BK, Takenaka N, Corces VG. Genome-wide phosphoacetylation of histone H3 at Drosophila enhancers and promoters. Genome Res.. 2012;22:1081.-8 doi: 10.1101/gr.136929.111.
67. Yasui D, Miyano M, Cai S, Varga-Weisz P, Kohwi-Shigematsu T. SATB1 targets chromatin remodelling to regulate genes over long distances. Nature.. 2002;419:641.-5 doi: 10.1038/nature01084.
68. Bell R, Munro J, Russ C, Powell JF, Bruinvels A, Kerwin RW. Systematic screening of the 14-3-3 eta (eta) chain gene for polymorphic variants and case-control analysis in schizophrenia. Am J Med Genet.. 2000;96:736.-43 doi: 10.1002/1096-8628(20001204)96:63.0.CO;2-2.
69. Taya S, Shinoda T, Tsuboi D, Asaki J, Nagai K, Hikita T. DISC1 regulates the transport of the NUDEL/LIS1/14-3-3epsilon complex through kinesin-1. J Neurosci.. 2007;27:15.-26 doi: 10.1523/JNEUROSCI.3826-06.2006.
70. Cheah PS, Ramshaw HS, Thomas PQ, Toyo-Oka K, Xu X, Martin S. Neurodevelopmental and neuropsychiatric behaviour defects arise from 14-3-3zeta deficiency. Mol Psychiatry.. 2012;17:451.-66 doi: 10.1038/mp.2011.158.
71. Ikeda M, Hikita T, Taya S, Uraguchi-Asaki J, Toyo-Oka K, Wynshaw-Boris A. Identification of YWHAE, a gene encoding 14-3-3epsilon, as a possible susceptibility gene for schizophrenia. Hum Mol Genet.. 2008;17:3212.-22 doi: 10.1093/hmg/ddn217.
72. Genome-wide association study identifies five new schizophrenia loci. Nat Genet. 2011;43:969.-76 doi: 10.1038/ng.940.
73. Vassos E, Steinberg S, Cichon S, Breen G, Sigurdsson E, Andreassen OA. Replication study and meta-analysis in European samples supports association of the 3p21.1 locus with bipolar disorder. Biol Psychiatry. 2012;72:645.-50 doi: 10.1016/j.biopsych.2012.02.040.