香京julia种子在线播放

All experimental procedures were approved by the local ethical committee (Procedure Code: ISA-11-1358). Wild type mice (C57BL/6J) and Tg(Thy1-YFP) [strain B6.Cg-Tg(Thy1-YFPH)2Jrs/J No. 003782; The Jackson Laboratories] were kept and bred according to local (Catalan law 5/1995 and Decrees 214/97, 32/2007) and European regulations (EU directives 86/609 and 2001-486).

After weaning (21 days of age), mice were randomly reared under either non-enriched (NE) or enriched (EE) conditions for 30 days. In NE conditions, animals were reared in conventional Plexiglas cages (20 × 12 × 12 cm height) in groups of two to three animals. The EE group was reared in spacious (55 × 80 × 50 cm height) Plexiglas cages with toys, small houses, tunnels, and platforms. The arrangement was changed every 3 days to maintain the novelty of the environment. To stimulate social interactions, six to eight mice were housed in each EE cage. All groups of animals were maintained under the same 12 h (8:00 to 20:00) light-dark cycle in controlled environmental conditions of humidity (60%) and temperature (22°C), with free access to food and water. The experiments were conducted using only females, since male mice showed hierarchical behavior similar to that observed previously that may affect the outcome of EE (Martínez-Cué et al., 2002). To stimulate social interactions, six to eight mice were housed in each EE cage. The use of female mice allowed to mix different mice in the same cage from different litters.

Cortical data analysis derived from the same samples and EE protocol that was used by Pons-Espinal et al. (2013) and PhD thesis “Role of DYRK1A in hippocampal neuroplasticity: Implications for Down syndrome¹,” where behavioral studies (Morris water maze testing) showed successful EE treatment. At 5 (P28) or 8 weeks of age (P51), mice were euthanized by carbon dioxide, and the cerebral cortex was dissected within 1 min of death. The tissue was immediately flash frozen in liquid nitrogen. In the case of Tg(Thy1-YFP) animals, the tissue was immersed in Hank’s Balanced Salt Solution (HBSS 1X, Gibco 14065-049) before proceeding with the sample preparation for the FAC-sorting (see section below). For each condition and replicate, we pooled the cortices of five mice with some exceptions: mice for insitu-HiC where single replicas as well as Thy-YFP mice for ATACseq (N_cortex = 20 mice in 4 bioreplicates; N_sorted = 30 mice in 6 bioreplicates, N_{NeuN+_HiC} = 4 mice in 4 bioreplicates, N_Thy+ = 4 mice in 4 bioreplicates). The frozen cortices for pooled animals were ground together in a frozen mortar containing liquid nitrogen, to obtain a fine powder of pooled cortex tissue. The powder was aliquoted and stored at −80°C.

DNA was extracted using Phenol:chroloform:iso-amyl alcohol (25:24:1) according to manufacturer guidelines (Sigma 77617-500ml). RNA was extracted using Qiazol total RNA (Qiagen Cat No:79306) kit according to the manufacturer’s instructions. The RNA was quantified by Qubit^® 2.0 Fluorometer (Life Technologies) and the quality was assessed using a Nanodrop 2000c (Thermo Fisher Scientific) and a 2100 Bioanalyzer (Agilent Technologies, CA, United States).

To obtain fresh nuclei, ground frozen tissue was resuspended in tissue lysis buffer (1x PBS containing 0.1% Triton X-100, 1x Complete Protease inhibitor cocktail tablette (cOMPLETE mini EDTA-free, Roche Cat No. 11836170001) and 1 mg/ml AEBSF (Pefabloc, Roche Cat No. 11585916001) and dissociated by 60–90 strokes in a glass douncer (7 ml tissue grinder Tenbroek, Wheaton Cat No. 357424). Nuclei were counted using a hemacytometer and constantly checked under a microscope (Leica DM-IL).

Two different procedure were performed: (1) Sorting total neurons using NeuN (Rbfox3) marker and (2) Sorting pyramidal neurons in Tg(Thy1-YFP) mice. Briefly, after the nucleus preparation, nuclei were resuspended in 1 ml of PBS-PI 1X [1X-PBS, 1x Complete Protease inhibitor cocktail (cOMPLETE mini EDTA-free, Roche Cat No. 11836170001), 1 mg/ml AEBSF Pefabloc (Pefabloc, Roche Cat No. 11585916001) and 0.1 mg/ml of BSA (BSA, Molecular Biology Grade NEB, Cat No. B9000S)]. 1.5 μl of anti-NeuN, clone A60, Alexa Fluor^® 555 Conjugate (Merk Millipore Cat No. MAB377A5) was added to the solution and incubated at 4°C for 1.5 h protected from light. The sample was centrifuged for 10 min, 500 g at 4°C and washed with 1 ml of PBS-PI 1X. Next, 1μl of 4′,6-Diamidine-2′-phenylindole dihydrochloride was added (DAPI, Roche Cat. No. 10236276001) and the sample was given immediately to the FACS-sorting Facility. Samples were sorted at 12PSI in cold condition in an INFLUX sorter (BD Biosciences INFLUX^TM). The sorted samples were centrifuged for 40 min at 700 g at 4°C to collect the nuclei before proceeding with the desired technique. For sorting pyramidal neurons Tg(Thy1-YFP), animals were dissected and tissue was immediately submerged in Hank’s Balanced Salt Solution (HBSS 1X, Gibco 14065-049). Brain samples were dissociated using the Neural Dissociation Kit (P) (MACS Milteny Biotec Cat. No. 130-092-628; LS Columns Cat. No 130-042-401; Myelin removal beads Cat No. 130-096-731; MidiMACS^TM separator Cat. No. 130-042-302), according to manufacturers’ instructions. Cells were sorted using an INFLUX^TM sorter (BD Biosciences INFLUX^TM). After sorting, samples were centrifuged for 40 min at 700 g at 4°C to collect the nuclei before proceeding with the desired technique.

WGBS was performed by CNAG Genome Facility on two independent sets of biological replicates. Briefly, genomic DNA (1–2 μg) was spiked with unmethylated λ DNA (5 ng of λ DNA per microgram of genomic DNA; Promega). DNA was sheared by sonication to 50–500 bp in size using a Covaris E220 sonicator, and fragments of 150–300 bp were selected using AMPure XP beads (Agencourt Bioscience). Genomic DNA libraries were constructed using the Illumina TruSeq Sample Preparation kit following Illumina’s standard protocol. DNA was treated with sodium bisulfite after adaptor ligation, using the EpiTexy Bisulfite kit (Qiagen), following the manufacturer’s instructions for formalin-fixed, paraffin-embedded tissue samples. Two rounds of bisulfite conversion were performed to ensure a conversion rate of >99%. Enrichment for adaptor-ligated DNA was carried out through seven PCR cycles using PfuTurboCx Hot-Start DNA polymerase (Stratagene). Library quality was monitored using the Agilent 2100 Bioanalyzer, and the concentration was determined by quantitative PCR with the library quantification kit from Kapa Biosystems. Paired-end DNA sequencing (2 × 100 bp) was then performed using the Illumina HiSeq 2000 platform.

Open chromatin studies were performed by ATAC-seq and SONO-seq procedures. Briefly, ATAC-seq was performed with minor modifications from Buenrostro et al. (2013). 100,000 nuclei were treated with 2.5 μl Tn5 at 37°C for 30 min, followed by cleanup on a Qiagen Minelute column. Fragments >1 kb in size were removed using AmpureXP beads (Agencourt AMPure XP, Beckman Coulter Cat. No. A63881). DNA fragments were amplified by 11 cycles of PCR with custom adapter primers from Buenrostro et al. (2013). PCR reactions were cleaned up with AmpureXP beads, quantified by Qubit and quality controlled by Bioanalyzer (Agilent Technologies). SONO-seq consists of isolating and fragmenting crosslinked chromatin, before reversing crosslinks, purifying the DNA and preparing it for sequencing (Auerbach et al., 2009). Chromatin was fragmented by sonication using the same COVARIS specifications as for ChIP-seq (see above) to obtain a median fragment size of 200 bp. Sequencing libraries were prepared using the NEBNext Ultra (New England Cat. No. E7370L) kit according to the manufacturer’s protocol. Both techniques were sequenced on a HiSeq2000 sequencer (Illumina).

Nuclei obtained in section “Nucleus Isolation” were cross- linked with 0.5% formaldehyde (Sigma F8775-25ml) for 5 min at room temperature (RT). Residual formaldehyde was quenched by addition of glycine (MAGnify^TM Glycine P/N 100006373) to a final concentration of 0.125 M and incubation for 5 min at RT. Nuclei were pelleted by centrifugation at 500 g during 10 min at 4°C and resuspended in 300 μl lysis buffer (MAGnify^TM Chromatin Immunoprecipitation System, Cat no. 49-2024). Chromatin was fragmented by sonication in a Covaris S2 [Duty Cycle: 20, Intensity: 8, Cycles per Burst:200, for 15 min (histone marks), for 25 min (FACS-sorted nuclei)] to a median size of 200 bp, aliquoted and stored at −80°C until further use. For non-histonic proteins such CTCF, no nuclei preparation was performed. Homogenized tissue was crosslinked with 0.5% formaldehyde for 10 min at RT, quenched and fragmented as above (Duty Cycle: 5, Intensity: 2, Cycles per Burst: 200, Time: 25 min). Chromatin immunoprecipitation was performed using antibodies against histone modifications (H3K27ac, H3K4me3, H3K4me1, H3K79me2, H3K36me3, H3K27me3, H3K9me2, and CTCF) with the MAGnify^TM Chromatin Immunoprecipitation System (Invitrogen Cat no. 49-2024), according to manufacturer’s instructions. For whole cortex and NeuN histonic ChIP-seq a total amount of 50 k nuclei was used per ChIP (∼330 ng), 700 k nuclei (∼4 μg) for non-histonic ChIP-seq. Recovered ChIP DNA was used to construct sequencing libraries, using the NEBNext Ultra (New England Cat. No. E7370L) kit according to the manufacturer’s protocol, and sequenced on a HiSeq2000 sequencer (Illumina). The quality of the ChIP-seq was determined by qPCR, using positive and negative primers to detect the regions where the histones should be placed in the genome (Supplementary Table 1). Primers were diluted to a final concentration of 300 ng in Power SYBR Green PCR MM (Applied Biosystems Cat. No 4367659). Samples were run in a Applied Biosystem qPCR system (7900 HT Fast Real-Time PCR System) as follows: 50°C/2 min, 95°C/10 min, 40 cycles of 95°C/15 s, 60°C/1 min, 95°C/15 s, 60°C-15 s, and 95–15 s.

Transcriptome study involved both poly-A RNA, directional RNA and small RNA libraries. Poly-A RNA sequencing libraries were prepared from total RNA using the TruSeq^TM RNA sample preparation kit (Illumina Inc., Cat. No. RS-122-2001) according to the manufacturer’s protocol. Directional RNA libraries were prepared using the ScriptSeqTM Complete Gold Kit (Human/Mouse/Rat) (Epicentre Biotechnologies), according to the manufacturer’s protocol. Briefly, 3 μg of total RNA were depleted of both cytoplasmic and mitochondrial rRNAs using the Ribo-Zero^TM Gold rRNA Removal Reagents. The total rRNA depletion of the samples was confirmed on a 2100 Bioanalyzer RNA 6000 Pico Chip. For the library preparation we used 50 ng of Ribo-Zero-treated RNA as starting material for the ScriptSeqTM v2 RNA-Seq Library Preparation Kit, followed by amplification by 12 cycles of PCR, using the FailSafeTM PCR Enzyme Mix (Epicentre Biotechnologies) before purification with AMPure XP beads (Agencourt, Beckman Coulter Cat. No. A63881). Both the directional mRNA and the Poly-A RNA libraries were sequenced in paired end mode with read length 2 × 101 bp on a HiSeq2000 (Illumina, Inc.) following the manufacturer’s protocol. After computational analysis, we validated 21 differential expressed genes in a new batch of biological replicates following the method of Schmittgen and Livak 2008 (see Supplementary Table 4). To validate the data generated in the RNA-seq analysis, 1 μg of the sequenced RNAs were used to prepare cDNA with SuperScriptIII (Invitrogen) according to manufacturer’s instructions, cDNAs were normalized to 100 ng/μl. RT-PCRs for the alternative splicing events were performed using oligos annealing to the adjacent constitutive exons and performed under standard conditions; 2% agarose gels were used to resolve the different bands. Image J software was used for quantification of the observed bands and determination of the PSIs for each event. Small RNA libraries were generated using the TruSeq (Illumina, Cat. No. RS-122-2001) kit according to the manufacturer’s protocol. The resulting 22 bp insert libraries were sequenced on a HiSeq2000 sequencer (Illumina), yielding 15–20 million reads per sample. After the analysis, we validated 6 out of 10 miRNA in a second independent group of animals, following Chen et al. protocol with minor modifications: instead of using Taqman probes we designed our own RT-miRNA oligonucleotides and performed qPCRs with Power SYBR Green PCR MM (Applied Biosystems Cat. No. 4367659, see Supplementary Table 5; Chen et al., 2005).

Samples were minced with RIPA-M buffer (1% NP40, 1% Sodium deoxycholate, 0.15 M NaCl, 0.001 M EDTA, 0.05 TrisHCl pH = 7.5, 1X cOMPLETE Mini EDTA free, 0.01 M NaF, 0.01 M Sodium pyrophosphate, 0.005 M β-glycerolphophate), sonicated with a Diagenode Bioruptor (cycles of 0.5 min ON 0.5 min OFF, medium intensity during 5 min). Samples were centrifuged during 10 min 16,000 rpm at 4°C and precipitated with acetone at -20°C for 1 hour. Samples were pelleted by centrifugation during 10 min 16,000 rpm at 4°C, dried and resuspended in Urea/200 mM ABC, sonicated again during 10 min (cycles of 0.5 min ON 0.5 min OFF, medium intensity) and quantified prior to mass spectrophotometry injection following isobaric tags for relative and absolute quantitation (iTRAQ) or Liquid Chromatography/Mass-Spectophotometry (LC/MS).

Cerebral cortex samples from individual C57BL/6J mice (2 bio-replicates per EE and CTL conditions) were sorted using NeuN⁺ (Rbfox3+) as described above. After sorting, approximately 1 million of nuclei were used for in situ HiC following previous specifications (Rao et al., 2014). Libraries were sequenced on a HiSeq2000 (PE × 125 bp) yielding approximately 300 M of reads per sample.

Active enhancers for EE and CTL cortex were predicted and annotated using a machine learning approach called Generalized Enhancer Predictor (GEP)² (Jhanwar et al., 2018). The method performed classification of epigenetic patterns coming from cortical ATAC-seq, H3K4me1, H3K4me3, H3K27ac, H3K36me3, and H3K27me3 using Random Forest (RF) and Support Vector Machine (SVM) classifiers to build predictive models for identification of active enhancers. The total amount of enhancers used in the present study corresponded to the consensus of GEP prediction for EE, control and ENCODE data, resulting in 347112 enhancers (Supplementary Table 1).

Chromatin interactions were predicted from chromatin modifications using Epitensor with minor modifications to adapt the script to the mouse genome (Zhu et al., 2016). We used epigenetic data from following tissues from the mouse ENCODE project (forebrain, heart, hindbrain, kidney, liver, lung, midbrain, stomach) (Yue et al., 2014). As well as in-house generated data from cortex of animals with and without environmental enrichment. We used the following epigenetic marks: H3K9me3, H3K4me3, H3K4me1, H3K36me3, H3K27me3, H3K27ac, CTCF as well as RNA-seq coverage data. Chromatin accessibility was measured by DNase-seq for ENCODE tissues and ATAC-seq for in-house cortex samples. Promoter regions were defined using the TxDb.Mmusculus.UCSC.mm10.knownGene package in BioConductor as 1,500 up- and 500 down-stream of any possible TSS. Active enhancers in each of the input tissues were defined by training a machine learning classifier on a list of validated enhancers using core epigenetic mark intensities as features (DNase, H3K27ac, H3K27me3, H3K36me3, H4K4me1, H3K4me3 as well as intensity ratio between the last two marks). The trained classifier was then applied to the epigenetic data from ENCODE tissues and local cortex datasets. To limit computation load, possible interactions were limited to intra-TAD interactions, which were based on the set of mouse cortex TADs³ (Dixon et al., 2012). Differential activity in enhancers was annotated using annotate.enhancers.with.genes.sh utility⁴.

Methylated CpGs were called from the raw reads using Bismark (Krueger and Andrews, 2011) following the user guide. Differential methylation analysis was carried out using bsseq (Hansen et al., 2012). Briefly, CpG methylation values were locally smoothed using the BSmooth function and CpGs that had a coverage of less than 5 reads in any of the samples were removed. We calculated the t statistic for the smoothed CpG values using the BSmooth.tstat function with paired design in local correction mode and the dmrFinder function was used to identify differentially methylated regions (DMRs) that contained multiple CpGs with a t statistic below -4.5 or above 4.5. DMRs located inside or within 5 kb of a gene or enhancer were annotated accordingly. Additionally, DMRs that overlapped enhancers, were annotated with the enhancer’s target gene according to the Epitensor predictions.

ATAC-seq libraries were aligned to mm10 using bwa-mem with predefined parameters (Li and Durbin, 2010). Duplicate read pairs were marked using the MarkDuplicates command from the Picard software suite⁵. We used a peak-independent approach on the one hand using our predicted enhancers and promoters as defined above, and on the other hand a peak-dependent method using F-seq or MACS2 with default parameters with and without duplicates (Boyle et al., 2008; Zhang et al., 2008). We also used MACS2 with the shifted strategy with the following parameters: “–nomodel –shift -75 –extsize 150 –broad –keep-dup all.” Together with enhancer and promoter regions, each annotation was loaded into DiffBind (Ross-Innes et al., 2012) providing as background input SONO-seq chromatin (Ross-Innes et al., 2012). For SONO-seq libraries, the pipeline followed the same steps as the ChIP-seq (see below).

Footprinting analysis was done using the Centipede software using the core transcription factor binding motifs from the Jaspar database (version 2016-03-02) (Pique-Regi et al., 2011; Khan et al., 2018). Instances of each motif in the genome were kept if the PWM score was superior or equal to 13. We created and used a mm10 mappability file to filter out instances that are located in regions of the genome with low mappability (gem-mappability-retriever; Marco-Sola et al., 2012). A motif was considered bound if the P-value-Z-score-combined statistics was inferior to 0.05. A motif was considered differentially bound if the ANOVA p-value was inferior to 0.05. Individual instances of a motif were considered bound when their posterior probabilities were superior to 0.99.

Samples were mapped to mouse mm10 (peak-independent) using Bowtie (Langmead et al., 2009) with “–quiet –sam –best –strata -m 1” parameters. Sam files were converted to bam files allowing only reads with mapping quality greater than 30 (“-q 30”). Files were visualized using bamCoverage utility from deepTools (Ramírez et al., 2016). We initially used a peak-dependent strategy as a benchmarking method and also as a tool to define our peak-independent approach (data not shown). The quality of ChIP-seq was assessed by the irreproducible discovery rate (IDR) (Li et al., 2011). The peak-independent method used in the present study consisted in the calculation of the coverage at defined regions as follows. The broad marks H3K36me3 and H3K79me2 were measured along the full gene body, defined as the distance from the TSS to the TES. Enhancer regions were provided by the GEP enhancer predictor in combination with ENCODE enhancers (Jhanwar et al., 2018; Supplementary Table 1). Promoter regions were defined as the interval from 1,500 bp upstream to 500bp downstream of the TSS. Reads were quantified by featurecounts (Liao et al., 2014) with the following parameters: “–ignoreDup –minReadOverlap 25 -Q 1 –O.” Counts tables were supplied to the batch effect corrector RUVseq for further differential analysis using edgeR (Robinson et al., 2010; Risso et al., 2014). A binned approach was used to determine differential changes in 1Mb bins to be associated with inter-chromosomal interactions hubs. This differential analysis was performed with Diffreps (Shen et al., 2013). Using the following parameters: “-pval 0.001 -frag 150 -window 1000.” Batch effects are a major issue in sequencing studies. In this study, we acknowledge the lack of bio-replicates in the cortical tissue by having a wide set of different techniques as well as a pooling strategy of animals per bio-replicate (5 individuals per sample). To minimize batch effects, we randomized the samples, applied standardized procedures and parallelized the experiments as much as possible. However, FAC-sorting experiments could not always be parallelized for different reasons (i.e., availability of the sorter). To remove batch effects, we decided to use the strategy of Russo et al. that was also used in data different from RNAseq (Risso et al., 2014; Koberstein et al., 2018). A principal component analysis was the main criterion to evaluate each dataset and the requirement for batch correction. If a PCA on uncorrected data could clearly separate the conditions, we ruled that there was no need for batch correction. However, if required, we then selected the method(s) (among RUVg, RUVs, and/or RUVr) that were able to separate conditions and intersected the results with a FDR threshold of 0.1 to obtain a conservative final gene list of changes induced by EE. Coverage plots were produced using the function bamCoverage of deepTools to produce bigwig files (Ramírez et al., 2016) (v2.0 with Python 2.7.14). These files were supplied to the function computeMatrix using the “scale-regions” parameter that allows to plot the coverage profile using the plotProfile function along regions of interest. The cell-specificity of sorted populations was assessed by the tool MakerGeneProfile that estimates cell specificity using a curated single-cell mouse brain RNAseq database⁶. We transformed gene-body assigned reads of H3K79me2 NeuN⁺ and NeuN^– into RPKM values and ran MakerGeneProfile to compare neuronal and non-neuronal enrichments to oligodendrocytes, astrocytes and pyramidal neurons specific markers. MGP states for MakerGeneProfile and it is described as a correlate value of cell type proportions. Specifically, MGP scores consist in the summarized expression profiles into the first principal component to estimate cell ratios (Mancarci et al., 2017). Gene ontology term enrichment analysis was performed using the Cytoscape tools clueGO and cluepedia, Metascape and SynGO (Bindea et al., 2009; Koopmans et al., 2019; Zhou et al., 2019). We performed the analysis individually for each histone mark or together, supplying upregulated and downregulated DBS separately. We used Bonferroni step down or Benjamini-Hochberg with p-value thresholds of 0.05 or 0.01 depending of the amount of data provided. In general, larger datasets needed more stringency in the statistical correction (Benjamini and Hochberg, 1995).

mRNA reads were aligned to mm10 using STAR with standard parameters (Dobin et al., 2013). Reads were counted using featureCounts by “gene_name” (Liao et al., 2014), batch corrected by RUV-seq before differential analysis in edgeR (Robinson et al., 2010; Risso et al., 2014). We also provide a splicing analysis results compiled in Supplementary Table 4. To identify and quantify alternative splicing, we used vast-tools v1.0.0 (Tapial et al., 2017). To increase effective read coverage at splice junctions, we pooled all replicates for polyA, respectively ribominus, samples using vast-tools merge. Differential splicing analysis was performed using vast-tools compare, comparing EE with Ctl samples in a paired manner, using the parameters –min_dPSI 10 –min_range 5 –p_IR –noVLOW. This resulted in, respectively, 40 and 53, cassette exons being up- or downregulated with EE (including 2 and 1 microexons; Irimia et al., 2014), 49 and 43 retained introns, 27 and 22 alternative 3’ss choices and 28 and 23 alternative 5’ss splice choices.

Regarding the lncRNA analysis, total RNA reads were aligned to mm10 using Subjunc splice-aware aligner with default settings (Liao et al., 2013), and reads overlapping exons were summarized at the gene-level to the corresponding genes using featureCounts (Liao et al., 2014). Read assignment to exons was carried out in a strand-aware manner, only fragments with both mates correctly aligned were considered and genomic regions with multiple overlapping exons on the same strand were disregarded. The count matrix was further filtered to retain only GENCODE long non-coding RNA (GRCm38.p5_M15). Reads with an average log2 CPM < 0.9 across all samples were filtered out yielding 4,472 genes. Normalization factors were calculated using the TMM method from edgeR package (Robinson et al., 2010). Observational-level weights were calculated using voom (Law et al., 2014) and used to fit gene-wise linear models (Smyth, 2004). Multiple testing p-value adjustment was performed using the Benjamini-Hochberg method (Benjamini and Hochberg, 1995).

Small RNA reads with homo-polymer and low PHRED scores were removed using FASTQ-Toolkit and a custom script. SeqBuster was used to remove adapters and align using miraligner.jar with mouse miRbase v18 annotation (Pantano et al., 2010). The small RNA dataset presented a strong batch effect which could not be corrected by RUVseq (Risso et al., 2014). We therefore devised an alternative strategy which consisted of filtering microRNAs for low read coverage (<50 counts), normalizing the libraries by read counts per million (rpm), before contrasting EE against CTL batch-wise and considering exclusively reproducible direction of change. Additionally, we calculated z-scores validating partially previous strategy (Supplementary Table 5).

The analysis of iTRAQ data consisted in sorting the discrepancy (δ = φ/β) between both biological contrasts (φ = EE1/CTL1 and β = EE2/CTL2), following the premise that consistent results will be indicated by:

Discrepant results then will be considered as values far from δ = 1. Based on this consideration we select all the proteins that follow the condition and we established a threshold of ±0.1. The protein expression values from the LC/MS were log2-transformed and loess normalized using the normalize.ExpressionSet.loess function from the BioConductor package AffyPML. Differential expression analysis was conducted with a standard limma (BioConductor) pipeline by calculating sample weights, fitting a linear model for each gene across all samples and calculating moderated F-statistics. Unadjusted p-values were used to rank the proteins (Smyth, 2004).

For quality check, sequencing reads were mapped to the mouse reference genome assembly (mm10), artifacts were filtered, and library was ICED normalized using the Hi-C-Pro (Servant et al., 2015) (v2.9.0) (Supplementary Table 7). For visualization we used the hicpro2juicebox.sh utility to convert the previous into a ^∗.hic format to visualize heatmap interaction matrices in Juicebox (Durand et al., 2016a). Juicer_tools were used to calculate A/B compartments using eigenvector utility (Durand et al., 2016b). We called chromatin loops with HICCUP at 5 and 10 kb of resolution (Durand et al., 2016a). For differential analysis, resulting interactions from Hi-C-Pro were splitted into intra and interchromsomal interactions, sex chromosomes and self-interacting bins were removed. Then piped into the edgeR wrapper RUVseq were RUVr for intra and RUVg for inter were used to normalized batch effects.

Required EE and CTL gtrack file to run chrom3D (Paulsen et al., 2017) was produced using the chrom3D wrapper automat_chrom3D utility⁷. Domains were called using Arrowhead (Juicer tools 1.7.6; Durand et al., 2016b). “–ignore_sparsity” parameter was used, and calls could be only produced at not lower than 10 kb. For the present study, we finally selected 5 M iterations including the parameter “–nucleus” to force the beads to remain confined inside the designed radius: “-r 3.0.” Domain coloring was produced by automat_color⁸ that allows to color any region of interest in the model. Both in silico models are available Supplementary Table 7.

We used Metascape to intersect the totality of the data. Some of the results were clumped as the maximum amount of sets allowed is 30 (Zhou et al., 2019). Transcriptomic and proteomic changes percentages explain by the rest of the data were calculated out of the evidences table reported as Metascape result (Supplementary Table 8).

Linear dependency test was performed by using Pearson and Spearman correlations to test how enhancer and promoter activity of different marks influence differential changes observed in the transcriptome and proteome. Briefly, normalized counts by RPKM mapping into enhancers and promoters were first averaged by target gene name they interact with using EpiTensor (Zhu et al., 2016). Then a fold change was calculated of EE vs. CTL samples. These values then are correlated with the corresponding fold changes found in differential RNA-seq and proteomic analysis. We used Spearman for proteomic changes as the scale of the fold changes were considerably different coming directly from iTRAQ data and the method is less sensitive for these outliers.

In order to test the association of differential regions induced by EE with inter-chromosomal interactions we run a permutation analysis using the package regioneR. Both permutations shown in the study were performed with a total number of 100,000 iterations. GWAS brain trait studies are collected in detail in Supplementary Table 8.