Immune-related transcriptomic and epigenetic reconfiguration in BV2 cells after LPS exposure: an in vitro omics integrative study

doi:10.21203/rs.3.rs-3318179/v1

Download PDF

Research Article

Immune-related transcriptomic and epigenetic reconfiguration in BV2 cells after LPS exposure: an in vitro omics integrative study

https://doi.org/10.21203/rs.3.rs-3318179/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Background

Molecular alterations affecting microglia activity have been consistently associated with the inflammatory response. These cells can have pro or anti-inflammatory activity, phenotypes thought to be heavily regulated by epigenetic mechanisms. Still, little is known about the details on how epigenetic marks regulate the expression of genes in the context of an inflammatory response.

Methods

Through CUT&RUN, we profiled four genome-wide histone marks (HM) (H3K4me1, H3K4me3, H3K27ac and H3K27me3) in lipopolysaccharide exposed-cells and compared their distributions to control cells. Transcriptomic profiles were determined through RNA-seq and differentially expressed genes were identified and contrasted with the epigenetic landscapes. Other downstream analyses were also included in this study.

Results

Our results illustrate an effectively induced M1 phenotype in microglial cells derived from LPS exposure. We observed differential bound regions associated with genes classically involved in the inflammatory response in the expected direction according to each histone modification. Consistently, our transcriptomic analysis yielded a conspicuous illustration of the LPS-induced immune activity showing the up-regulation of Nf-κB induced mRNAs (TNF-α, nfκbiz, nfκbia) and other important genes (Marco, Il-6, etc.). Furthermore, we integrated both omics profiles and identified an important reconfiguration of the genome induced by LPS. The latter was depicted by 8 different chromatin states that changed between conditions and that, associated with unique clusters of differentially expressed genes, not only represented regulatory elements, but also underlined distinct biological functions (inhibition of morphogenesis; changes in metabolism, homeostasis, and cytokine regulation; activation of the inflammatory response).

Conclusion

This study exhibits important differences in the distribution of histone modifications in treated and control BV2 cells, constituting an epigenetic reconfiguration that leads to the inflammatory response. Also, it highlights the importance of these marks' regulatory role in gene expression and provides possible targets for further studies in the context of inflammation.

BV2 cells

Neuroinflammation

Epigenetics

Chromatin

Microglia

Histone modifications

Innate immunity

The evolution of organisms' interactions led to the development of mammals' immune systems, which rely on inflammation to cope with harmful external stimuli and maintain homeostasis [1, 2]. However, while acute inflammation is crucial for appropriate immediate responses, chronic inflammation has been linked to a range of diseases, including cancers, neurodegenerative and autoimmune diseases, as well as psychiatric disorders that were previously attributed mainly to chemical imbalances in the brain [3–6]. Indeed, inflammation plays a key etiological role in the occurrence of diseases across various organs, including the central nervous system (CNS) [1, 7], but the intricacies of this process are still a matter of intense investigations.

Recent research has highlighted the pathogenic contribution of microglial cells (MCs), the predominant type of glial cells in the CNS [6, 8], to chronic neuroinflammation and the onset and progression of CNS diseases [6, 9]. MCs originate from myeloid precursors and migrate to the CNS during embryonic development, where they are responsible for brain development and homeostasis [6, 9, 10]. Considered as the macrophage-like resident cells of the CNS, MCs can adopt, among a continuum of phenotypes, two well- defined states upon activation in response to noxious stimuli: the M1 or pro-inflammatory phenotype and the M2 phenotype that, conversely, predominantly executes anti-inflammatory activities [6, 11]. M1 state can be induced by interactions with external agents such as lipopolysaccharide (LPS) or IFN-γ and contributes to the initial immune response, while the M2 phenotype can be induced through IL-4 and IL-3 to repair tissue and solve inflammatory processes [6]. Irrespectively of the complexity of the spectrum of phenotypes that can be adopted, it is well-accepted that these cellular fates stem from several mechanisms that are still poorly understood, including epigenetic mechanisms [12–14].

Epigenetic mechanisms associated with the regulation of MCs’ phenotype include DNA methylation, non-coding RNA interactions, chromatin remodeling and post-translational histone modifications (PTHMs) [6, 15]. These mechanisms are involved in different stages of the cell life and contribute to the biological adaptation to variable immune challenges [14, 16]. In particular, research on PTHMs suggests that treatments inhibiting or potentiating histone modifying enzymes can modulate inflammatory responses in MCs [17]. For instance, H3K27ac (active) and H3K27me3 (repressive) marks have been directly linked to the induction of microglial immune response and to the subsequent phenotypic polarization [15]. Specifically, LPS exposure increases histone deacetylases (HDACs) activity in MCs, hindering the gene expression by inducing DNA coiling [15]. Simultaneously, it increases EZH2 methyltransferase’s activity and concomitant H3K27me3 levels in non-inflammation related genes promoting M1 polarization in MC [15]. However, some studies have produced contradictory results [15], perhaps arising from the differences between the models used to study MCs [12, 18].

Together, this suggests that epigenetic mechanisms are directly involved in the induction and evolution of immune responses in MCs. Nonetheless, the nature of this involvement and the role of specific PTHMs on regulation of gene expression in these processes are still unknown. To address these questions, we first explored the differential binding regions of four commonly studied PTHMs, including both repressive and active marks, before and after LPS exposure. We then assessed the differential transcriptional signature of BV2 cells before and after exposure to LPS and finally assessed the functional interaction between the distribution of PTHMs and transcriptional profiles induced by LPS treatment.

BV2 cells culture

Immortalized mouse BV2 microglial cell lines were purchased (American Type Culture Collection) and subsequently grown in Dulbecco’s Modified Eagle Medium supplemented with 10% Fetal Bovine Serum, 1% sodium pyruvate and 1% penicillin G and streptomycin (Gibco, Thermo Fisher Scientific.) at 37°C in a humidified atmosphere containing 95% air and 5% CO2. The medium was changed every 2 days and the cells were split with 0.25% trypsin when they reached 80% confluence and sub-cultured for further passages.

BV2 cells stimulation and RNA extraction

BV2 cells were seeded in 6-well plates and treated with LPS (1 µg/mL) or 0.9% saline (vehicle control) for 24 h. For each of four PTHM, six samples (3 treated with LPS and 3 controls) were prepared to assess their whole-genome distribution through Cut Under Targets and Release Using Nuclease (CUT&RUN) [19]. Twelve additional samples (6 treated with LPS and 6 controls) were prepared and treated the same way for the assessment of transcriptomic profiles through RNA sequencing (RNA-Seq). After stimulation, total RNA was extracted using the miRNeasy mini kit (Qiagen) with a spin column according to the manufacturer’s protocol, which included a digestion phase with DNase I (New England Biolabs). LPS and Complete™ Mini, EDTA-free Protease Inhibitor Cocktail tablets were obtained from Millipore Sigma.

Assessment of PTHMs distribution through CUT&RUN

CUT&RUN was performed following EpiCypher® CUTANA™ CUT&RUN Protocol (v1.6) (https://www.epicypher.com/resources/protocols/cutana-cut-and-run-protocol/). Briefly, 500,000 cells were harvested by scraping 24 h post treatment (LPS stimulation). Cells were washed with Wash Buffer (20 mM HEPES pH 7.5, 150 mM NaCl, 0.5 mM spermidine and one complete protein inhibitor tablet per 50 mL), attached to BioMag®Plus Concanavalin A-coated magnetic beads (Bangs Laboratories, Inc.) and incubated with 1 µg primary antibody targeting each PTHM (H3K4me1, H3K4me3, H3K27Ac and H3K27me3) diluted in Wash Buffer containing 0.01% digitonin (Digitonin Buffer) overnight at 4°C. In this step we used rabbit H3K4me1 (C15410194) and H3K27Ac (C15410196) antibodies (Diagenode) and rabbit H3K4me3 (9751S), H3K27me3 (9733S) and control IgG (66362S) (Cell Signaling Technology). The BV2 cells were then washed and incubated with 2.5 µg protein A-MNase (pAG-MNase) (EpiCypher) for 10 min at room temperature. The material was washed again and placed on an ice-cold block and incubated with Digitonin Buffer containing 2 mM CaCl2 to activate pA-MNase digestion. After 2h digestion at 4°C, one volume of 2X stop buffer (340 mM NaCl, 20 mM EDTA, 4 mM EGTA, 0.02% Digitonin, 0.05 mg/mL glycogen and 50 µg/mL RNase A) was added to stop the reaction and fragments were released by 30-min incubation at 37°C. Supernatant was recovered by placing the tube on a magnet stand until the mix was cleared. DNA fragments were extracted from the supernatant and purified using the MinElute PCR Purification Kit (Qiagen) per manufacturer’s instructions.

Library preparation and sequencing

CUT&RUN. Sequencing libraries were prepared from DNA fragments following the NEBNext UltraTM II DNA Library Prep Kit for Illumina protocol (New England Biolabs) with some modifications. Briefly, after adapter ligation, 1.1X volume of Axygen select beads (Fisher Scientific) was added to the reaction to ensure high recovery efficiency of short fragments. After 14 cycles of PCR amplification, the reaction was cleaned up with 1.1X volume of Axygen beads. The concentration and molarity of each of the 24 barcoded libraries were quantified using the Agilent High Sensitivity DNA Kit on an Agilent 2100 Bioanalyzer (Agilent Technologies) according to the manufacturer’s instructions. The libraries were sequenced at paired end with read length of 100 bp on the NovaSeq 6000 platform (Illumina). Samples were pooled 32 per lane and sequenced at a depth of 22 million reads per sample. Sequencing service was provided by the Plateforme de Séquencage de nouvelle génération in the Centre de recherche du CHU de Québec-Université Laval.

RNA-Seq. Quality and concentration of the RNA retrieved with the miRNeasy mini kit was assessed using Agilent RNA 6000 Nano Kit on an Agilent 2100 Bioanalyzer (Agilent Technologies) according to the manufacturer’s instructions. RNAseq libraries were constructed using the NEBNext Ultra Directional RNA Library Prep Kit for Illumina (New England Biolabs) preceded by ribosomal RNA depletion starting with 100ng RNA. Libraries were barcoded and sequenced in multiplex at pair-end with read length of 100 bp on NovaSeq 6000 platform (Illumina). Samples were pooled 12 per lane and sequenced at a depth of 22 million reads per sample. Sequencing service was provided by the Plateforme de Séquencage de nouvelle génération in the Centre de recherche du CHU de Québec-Université Laval.

Bioinformatic analyses

PTHM assessment. We first assessed the quality of the paired-end fastq files with the open-source softwares fastQC (Available online at: http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and MultiQC [20]. The main quality control outputs are shown in Additional File 1: Figure S1. Reads were trimmed following a two-step protocol starting with Trimmomatic [21] and ending with Kseq tool [22] as suggested by Zhu and collaborators [22, 23]. Samples were then aligned to the mm10 UCSC genome using Bowtie2 (https://bowtie-bio.sourceforge.net/bowtie2/index.shtml) [24, 25] as described before by Zheng and collaborators (dx.doi.org/10.17504/protocols.io.bjk2kkye). Fragment length distribution and between-sample Pearson correlations were assessed through a script adapted from [26] (Additional File 1: Figures S2 and 2A). To evaluate the strength of the relationships between samples we built two multiple regressions where the outcome was the pair-wise Pearson correlation coefficients and the independent variables were categorical variables designating the PTHM, the replicate, the type of activity of the mark and the condition (LPS vs saline) (See Additional file 1: Extended Methods for further details). Through several packages, files were converted to different formats (Bed, BEDPE, BigWig, bedgraph) to allow downstream manipulations (https://www.encodeproject.org/software/bedgraphtobigwig/) [28, 29].

Differential binding analysis (DBA) was performed using diffReps with a window size of 1000 bp, sliding step of 100 bp, and fragment size of 200 bp [30]. A FDR < 10% and P value < 0.0001 for the negative binomial tests were established as initial significance cutoffs. Using a custom R script, we performed Over representation analysis (ORA) of Gene Ontology (GO) terms associated to the differentially bound regions (DBR) identified by diffReps per PTHM using ChiPeakAnno (default parameters) and Genomic ranges packages of the Bioconductor library [31, 32]. Additionally, we explored the different intersections between DBR through manipulation of the “genomic ranges” objects (intersect function) and plotted the data with the ChiPeakAnno, ComplexUpset and UpSetR packages [32–34]. To see if the co-occurrence of two or more PTHMs was observed significantly more than expected by chance we simulated 1000 datasets of random regions with the same mean width and standard deviation that the DBR observed. We determined the pooled proportion of overlapping PTHMs across the simulated datasets and performed Fisher exact tests (FET) adjusted by BH to see if the proportion observed was significantly greater.

Subsequently, we used the ChromHMM software [35, 36] to evaluate the joint distribution of the different PTHMs. Briefly, the genome was partitioned into 200-bp bins for which each PTHM was assigned as present or absent based on the mean read count in the bin and a local Poisson distribution (p-value < 0.004). This was achieved with the BinarizeBam function and all the CUT&RUN samples. Furthermore, an 8-state model was chosen considering biological coherence with previous annotated states in similar cells and the redundancy of the states found in alternative 10, 15 and 16-states models. The LearnModel function was used with default parameters changing the number of states with the -printposterior flag to obtain posterior probabilities. We selected the bins for which the latter was > 0.5 (~ 99% of the genome), determined the differences of the genome-wide state distribution between LPS and control and annotated the different states based on previous findings in the literature [35–37]. We next identified the regions that most likely endured a state conversion (e.g., from quiescent state to active TSS) across conditions and performed ORA. Finally, we explored known motifs enrichment (ME) in the promoter regions from the genome that overlapped with the regions that changed from quiescent state to TSS and from intergenic regions that overlapped with regions that changed from poised enhancer to enhancer using the SEA module from the meme-suite [38]. The visualization of the reads across the Genome was done with IGV browser with merged RPKM normalized bigwig per PTHM in both conditions [39].

RNA-Seq Analysis. Quality of the paired-end RNAseq data was assessed with fastQC and MultiQC [20]. We used Trimmomatic [21] to remove the universal adapter sequences. Reads were aligned with kallisto [40] on the GRCm38/mm10 version of the mouse genome retrieved from GENCODE. A gene was considered the union of all its exons in any known isoforms, based on GENCODE annotation. Any reads that clustered in different transcripts were excluded from the analysis. The threshold for filtering out genes expressed at low levels was set at 5 reads minimum in ~ 90% (n = 10 out of 11) of the samples. As all the samples were sequenced simultaneously, we only accounted for the variance explained by the RNA integrity number independently of a specific threshold as suggested by Gallego Romero et al. [41] and for surrogate variables representing unknown sources of variation [42, 43]. Specifically, we found a significant linear association between the second principal component (PC, prcomp function), which explained more than 15% of the variance, and RIN value (p < 0.03). Therefore, RIN was added as a covariate in the generalized linear model (GLM). Additionally, regularized log transformed counts were evaluated through principal component analysis (PCA) and visual means which led to the exclusion of one of the samples (Additional file 1: Figure S3).

The differential expression analysis was performed with Deseq2 [44] (median ratio normalization, apelgm algorithm, Negative binomial GLM per gene, Wald test and Benjamini Hochberg (BH) adjustment). An individual gene was called differentially expressed if the adjusted P-value (BH) of its Wald test (t-statistic) was ≤ 0.05. Globally, this approach allowed the identification of genes differentially expressed in BV2 cells after LPS treatment while controlling for baseline variations in gene expression and for the effects of technical covariates across biological replicates using the sva package [42]. ORA of GO terms was carried out with cluster profiler [45, 46] on differentially expressed genes (DEGs) with significant enrichment fixed at p-value adjusted by BH < 0.05. Figures were generated with ggplot2 [47], EnhancedVolcano [48] and enrichplot [49] packages.

Integrative analysis. To evaluate the relationship between the transcriptomic and the epigenomic profile we first sought to assess the enrichment of individual PTHMs according to the level of expression. We ranked the genes from RNA-seq based on their basal expression across samples, identified their promoter regions (TSS from GENCODE ± 2000Kb), exons (as annotated by GENCODE), introns (CDS after subtracting exonic regions) and known microglial enhancers (retrieved from the Enhancer Atlas 2.0 under “mus musculus”, http://www.enhanceratlas.org/data/download/enhancer/mm/Microglia.bed [50]) and determined their corresponding binding score from normalized RPKM tracks for each PTHM. Afterwards, we divided genes into highly and lowly expressed according to the median value of expression across all samples. The enrichment of highly and lowly expressed genes in specific genomic features was evaluated through the ChromHMM’s OverlapEnrichment function using the previously 8 states learned model. Similarly, we explore the overlap between DEGs and the 8 model states. In both cases, saline and LPS conditions were considered when the assessment was made. We also identified the regions that were differentially annotated between control and LPS and considered them as “transitions” induced by the immunological challenge. We calculated the frequency of these transitions and evaluated if their observed occurrence was significantly greater than that expected by chance if all transitions had the same probability of occurrence (112 possible transitions) based on the binomial distribution. Finally, we explored the biological function of DEGs associated with specific transitions.

Statistical analyses

Although sample size calculation was not performed, the sample size in this study is justified based on several previously published reports using similar or even smaller sample sizes and showing enough power to detect significant statistical differences. In this study, we included three replicates per condition per PTHM (CUT&RUN) and six replicates per condition for RNA-seq. Alignment proportions, number of reads and other quality checks can be seen in the supplementary information (Additional file 1: Figure S2). Normality assumptions (for linear modelling and Pearson correlation) were evaluated through visual means although given the sample size we could fairly assume a normal sampling distribution according to the Central Limit Theorem. As aforementioned, DBR analysis was performed using the diffreps software. However, before running the analysis the distribution of counts was assessed through histograms and qqplots to see if the distribution diverged significantly from a negative binomial distribution with mean estimated across samples and different dispersion values. Diffreps, includes a normalization step that is clearly described by Shen et al. Briefly, the package calculates numeric factors to linearly scale each sample based on the geometric mean and the median of the factors [30]. Deseq2, used in this case for RNA-seq, uses the median ratio normalization which is very similar to the method used by diffreps but does not include a window segmentation [30, 44]. The two methods are similar to the relative log expression normalization proposed by Anders and Huber [43, 44, 51]. The level of significance α was established at 0.05 and the adjustment method for multiple comparisons was BH except stated otherwise. All the analyses were carried out integrating bash and R 4.2.2 built-in scripts and ran in an Ubuntu 20.04 local machine. Only the meme-suite was used directly online.

The main objective of this study was to define the epigenetic patterns underlying the acute establishment of the active immune phenotype in microglial cells. To do so, we used a model of immune activation in BV2 cells induced by the administration of LPS. We first established the epigenetic profiles of four chromatin marks and evaluated their dynamics following LPS treatment. We paired these experiments with the analysis of the transcriptional signatures induced by LPS. Finally, by combining results from both analyses, we characterized the dynamic remodeling of epigenetic profiles underlying the transcriptional reorganization required for the transition from a naive to an activated immune state in immortalized MCs.

Epigenetic reprogramming after LPS exposure: differential binding analysis.

We started by assessing whether LPS treatment interferes with the epigenetic signatures of four active (H3K4me3; H3K4me1, H3K27ac) and repressive (H3K27me3) histone marks across the genome. To do so, we first mapped the genome-wide distribution of each histone mark and compared control with LPS-treated BV2 cells’ profiles through differential binding analysis (DBA). Our analysis revealed a large number of sites with differential occupancy for every chromatin mark across the genome (Fig. 1a, left and additional file 2: Tables S1 and S2). The greatest number of differentially bound regions were identified for the functionally active histone marks H3K4me3 (12860; 82% increase vs 18% loss), H3K27ac (8953; 41% increase vs 59% loss) and H3K4me1 (5808; 61% increase vs 39% loss) while the repressive mark H3K27me3 (5531; 59% increase vs 41% loss) exhibited the smaller number of regions affected by LPS treatment. These sites of differential occupancy were found across the genome, within intergenic, regulatory, and coding regions, with important differences between individual marks (Fig. 1c, 1e and additional file 1: Figure S4-6). Globally, these results are consistent with previous findings in the same or similar cellular types [52–55] and suggest that the activation of MCs by LPS is underlined by a considerable reorganization of active and repressive chromatin marks across the genome.

Although important differences were identified at the single-PTHM level, remodeling of the epigenetic landscape involves simultaneous changes of several marks. Accordingly, our analysis suggests that the probability of differential co-occupancy by the combinations of two active PTHMs is more than it would be expected by chance considering the genome-wide occupancy of each mark (H3K4me3 and H3K4me1: p_adj ≈1.6e-155; H3K4me3 and H3K27ac: p_adj ≈1.8e-187; H3K27ac and H3K4me1: p_adj ≈ 4e-73, Fig. 1a, right and additional file 1: table S3). In contrast, the probability of finding differential enrichment of H3K27me3 along with any of the other active marks is low (H3K4me3: p_adj ≈0.99 and H4K4me1: p_adj ≈0.55), with the exception of H3K27ac (H3K27ac and H3K27me3: p_adj=6.4e-9), findings consistent with previously reported results [54, 56, 57]. Figure 1d highlights a representative example of an epigenetic reorganization involving every chromatin mark at the Apol9a gene. Apol9a encodes the apolipoprotein L 9a protein, a phosphatidylethanolamine-binding protein that has previously been associated with anti-inflammatory responses induced by IFN- α [58]. Our results show that LPS treatment increases the occupancy of all three active chromatin marks in the gene’s distal and proximal promoter region while concomitantly depleting H3k27me3 in the same region (Fig. 1d). Together, this further supports the distinct contribution of these marks on the global molecular reorganization underlying the transition toward an active immune state in BV2 cells.

We next computed the nearest bidirectional genes for every site of differential occupancy and performed Over representation analysis (ORA) of GO terms to reveal the functional mechanisms associated with either enrichment or depletion of each chromatin mark following LPS treatment (Fig. 1d, 1f and additional file 1: Figures S4-6). Our results on single active marks highlighted terms enriched with genes relevant to immunity related to interferon and NF-κβ signaling (H3K4me3: p_adj< 4.1e-6), response to external stimulus (H3K4me1: p_adj< 8.3e-5), cytokine production (H3K27ac: p_adj< 1.1e-5), regulation of cytokine production and myeloid leukocyte activation (H3K27ac: p_adj< 2.2e-4) and intracellular signal transport (H3K4me3: p_adj< 3.7e-7). In contrast, our analysis of H3K27me3 revealed no term related to immunity but rather associated with nervous system development (p_adj< 5.1e-8), neurogenesis (p_adj< 1.3e-6) and organ morphogenesis (p_adj< 9.2e-6) (Additional file 5). Further analysis of the regions enriched with co-occurring differentially represented chromatin marks revealed functional terms related to immune response, although these terms were more specific (Additional file 1: Figures S4-6 and additional files 3,4,5 and 6), supporting the idea that specific combinations of PTHMs underly distinct functions required throughout the acute immune response and the subsequent anti-inflammatory phase. Together, this suggests that BV2 cells’ transition into an activated immune state requires a coordinated epigenetic remodeling of active and repressive histone marks promoting, in one hand, the molecular and cellular cascades promoting immune activation while, in parallel, repressing pathways related to non-immune functions.

Epigenetic reprogramming after LPS exposure: Epigenetic states’ transition.

Histone marks follow a specific pattern of co-occurrence across the genome which all together contribute to the establishment of the epigenetic landscape [59, 60]. The dynamic remodeling of the epigenetic landscape allows DNA-based task required for cell function and adaptation to environmental challenges [61]. We next assessed the level of occurrence and co-occurrence of active (H3k4me3, H3k4me1, H3K27ac) and repressive (H3K27me3) chromatin marks in BV2 cells in a genome-wide manner. Our analysis highlights marked patterns of co-occurrence, revealing strong positive mean correlation within (r ≈ 0.95; p_adj< 0.001) and amongst (r ≈ 0.8; p_adj< 0.001) active marks, with H3K4me3 and H3K27ac showing the highest level of co-occurrence (r ≈ 0.96; p_adj< 0.001, Fig. 2a and additional file 1: Tables S1 and S2). On the contrary, we found H3K27me3 to be inversely correlated with each of the three active marks (r≈-0.2; p_adj< 0.001, Fig. 2a). Together, these results are consistent with findings from our intersection analysis (Fig. 1A) and suggest that different chromatin states across the genome control gene expression and associated cellular functions [35, 53, 62].

We next used ChromHMM to estimate the identity and prevalence of chromatin states across the genome of BV2 cells and test how these states change following LPS treatment. We selected a model describing eight different states of co-occurrence throughout the genome of BV2 cells (Fig. 2b). For instance, the active enhancer state, defined by high levels of H3K4me1 and H3K27ac, low H3K4me3 and absence of H3K27me3, is found mainly at TES, introns, exons and gene bodies, while the active TSS state, enriched with H3K4me3, H3K27ac, H3K4me1 and depleted for H3K27me3, is found mainly at TSS enriched with CpG islands (Fig. 2b). In contrast, the repressed by polycomb state, enriched with the repressive mark H3K27me3 only, is depleted in TSS and CpG islands but enriched in intergenic regions, pattern similar to the quiescent state but more intense for the latter (Fig. 2b). Our data show that LPS induces a significant reorganization of these states across the genome. This results into a reduction of the quiescent state (control: 67.5% vs LPS: 66.5%) at the expense of other states such as active promoter (control: 3.8% vs LPS: 4.0%), flanking active promoter (control: 0.8% vs LPS: 1.6%), bivalent TSS (control: 0.8% vs LPS: 1.2%) and repressed by polycomb (control: 12.5% vs LPS: 14.1%; Fig. 2b). Interestingly, the transitions amongst enhancer-related states seem to be of bigger magnitude than those of TSS-related states. This suggests that LPS induces a coordinated reorganization of the epigenetic landscape resulting in the transition from quiescent inactive sites towards active promoter and enhancer sites allowing the expression of gene programs required for immune response in immortalized microglial cells as well as the tuning of the response in accordance with the type of harmful stimuli.

We next performed a gene ontology analysis using the nearest bidirectional genes associated with a change of state following LPS treatment. Interestingly, our results show that genomic regions transitioning from a quiescent state to an active state (n = 371) or from repressed by polycomb to an active TSS (n = 30) state after LPS treatment showed strong enrichment for genes involved in the immediate immune response terms such as response to interferon β (p_adj<0.0001), defense response to symbiont (p_adj<0.005), responses to cytokine (p_adj<0.007) and regulation of response to stress (p_adj<0.001) and external stimulus (p_adj<0.001; Fig. 2c and additional file 7). A representative example of such a transition is found in the promoter of Gbp6 (Fig. 2d). One of the variants of this gene is mainly identified as quiescent in the control cells, but an increase of both H3K4me3 and H4K27Ac occupancies drives an acquisition of an active TSS state in LPS-induced samples (LFC ≈ 3.92 and p_adj≈2.2e-66; Fig. 2d and 3c). Also, under control condition, the Rtp4 gene is maintained under a repressive state with a small bivalent TSS state in the promoter of one of its alternative variants. Remarkably, LPS treatment induced a strong enrichment of each active mark (H3K4me3, H4K27Ac and H3K4me1) while depleting the repressive H3K27me3 mark in the promoter region of this alternative variants, switching a bivalent promoter state into an active TSS state and increasing the expression of this variant of the Rtp4 gene (LFC ≈ 6.5 and p_adj≈0.0001; Fig. 2d and 3c).

Alternatively, transition from poised enhancer to the enhancer state after the LPS challenge (n = 30 486) was ten times more frequent than the abovementioned transitions and was mostly associated with terms related to intracellular cascade and immune response maintenance such as intracellular signal transduction (p_adj< 8.6e-12), metabolic processes (p_adj< 3.7e-8), cytokine production (p_adj< 3.0e-7), protein phosphorylation (p_adj< 3.8e-7) and biosynthetic processes (p_adj< 5.0e-7; Fig. 2c and additional file 7). For example, our analysis at baseline revealed a broad poised enhancer state, defined by the presence of H3K4me1 with low levels of any of the other marks, in the vicinity of the Ccl5 gene. Interestingly, LPS treatment increased the presence of H3K27ac and H3K4me3 in the region upstream Ccl5 which finally resulted in an active enhancer state, here again associated with elevated expression of Ccl5 (LFC ≈ 6.5 and p_adj≈3.0e-48; Fig. 2d and 3c).

In contrast, regions transitioning from an active TSS to a repressed by polycomb state after LPS treatment (n = 108) were mostly enriched with genes relevant to lipid catabolic processes (p_adj<0.005), neuron differentiation (p_adj<0.004), fatty acid homeostasis (p_adj<0.01) and histone phosphorylation (p_adj<0.01; Fig. 2c and additional file 7). Finally, we performed motif analysis (MA) in the centered regions that transitioned from quiescent to TSS active in promoter regions (n = 204) and those that transitioned from poised enhancer to enhancer in intergenic regions (n = 138 747). Interestingly, for the former we found the Interferon Receptor Factor 1 (p_adj<0.008) and 2 (p_adj<0.02) and other transcription factors like the Retinoid Acid Receptor Alpha (RARA, p_adj<0.02) (additional file 1: Fig S7A and additional file 8). In the second case, we found the Myc-associated zinc finger protein (Maz, p_adj<5.91e-995), the Retinoid X Receptor Alpha (RXRA, p_adj<2.48e-820) and the Retinoid X Receptor Gamma (RXRG, p_adj<6.12e-756) to be enriched (additional file 1: Fig S7B and additional file 8). Together, this suggests that the activation of BV2 cells after LPS follows a specific epigenetic reorganization involving a coordinated transition of chromatin states simultaneously promoting and repressing immune and non-related functions, respectively, via the recruitment of state-specific transcription factors.

Transcriptional reorganization induced by LPS exposure.

We assessed the extent to which LPS treatment would induce a consistent reorganization of transcriptional profiles in BV2 cells using RNAseq. Our analyses show that LPS and saline treated samples clustered very strongly (Fig. 3a). Further analysis confirmed that LPS treatment could explain roughly 71% of the overall variance (Fig. 3b) observed in our samples suggesting that LPS induces a strong a consistent impact on expression profiles in BV2 cells. We next used differential expression analysis to identify genes significantly up or downregulated by LPS. Overall, our results show that 1546 genes were significantly up regulated and 1419 downregulated (p_adj<0.1) following LPS treatment (Fig. 3c and additional file 9). Noticeably, Acod1, Lcn2 and Tnfrsf1b were three of the most significantly upregulated genes. This is consistent with previous reports showing that Acod1 is key to switch from an oxidative based metabolism to glycolysis during the immune response [63, 64] and Lcn2 is a key amplifier of the M1 polarization of MCs [65, 66]. Tnfrsf1b also participates in the TL4 cascade embedded in the immune response induced by LPS [67]. On the other hand, Arg1, Lpl and Fosb showed the strongest level of downregulation following LPS treatment. Consistently as well, Arg1 is known to be a hallmark signature of the anti-inflammatory phenotype in MCs [68–70]. ORA of GO terms confirmed that the functional terms most significantly associated with DEG are related to immune functions such as regulation of defense response (p_adj< 1.2e-18), cellular responses to cytokine stimulus (p_adj< 1.3e-22), cytokine production (p_adj< 6.6e-23) and response to lipids (p_adj<1.1e-21; Fig. 3d and additional file 10). Finally, we visually inspected the relationship between the 8 most enriched GO terms and DEGs. We can observe many crucial immune response related genes belonging to multiple enriched terms (Tnf, Oas2, Nfkb1, Irak3, etc.). Overall, the assessment of the transcriptomic changes induced by LPS suggests that we captured BV2 cells strongly polarizing into a pro-inflammatory phenotype (M1) with a few characteristics of the recently described M2a phenotype (for instance, differentially expressed IL-10 receptors) [71]. Our results also show that certain changes in the transcriptome only partially overlap with the epigenomic changes, which could correspond to the transcriptomic-epigenomic delay.

Integration of epigenetic modifications and expression

We studied the relationship between both epigenomics and transcriptomic profiles and evaluated the DEGs in accordance with changes in chromatin states after LPS exposure. Initially, we explored the abundance of each PTHM according to expression in different locations in the genome. In general, our results show that H3K27me3 abundance increased as expression diminished independently of the genomic feature, supporting the idea that its repressive effect is achieved through occurrence in very broad genomic regions rather than through peaks in specific genomic features [69]. In comparison, high abundance of H3K4me3 and H3K27Ac was found in highly expressed genes through genomic features but mainly in TSS and exonic regions. Interestingly, H3K4me1 mimicked the behavior of H3K4me3 and H3K27Ac in intronic and enhancer regions of genes according to their expression. In contrast, the level of H3K27me3 occupancy was highly correlated to lowly expressed genes in both conditions and in all the genomic features (Fig. 4a). There were no important differences regarding the relationship between expression and occupancy of the marks between control and LPS-induced cells. The little effect exhibited in Fig. 4a is likely due to the property of invariance for which RPKM transformation does not account for across samples [72]. Overall, this is consistent with previous reports on the role of PTHM on gene expression [35, 73] and supports the idea of different marks co-occurring in a synchronized way independently of the condition.

We next tested whether the eight chromatin states identified above were enriched for genes differentially expressed following LPS treatment (Fig. 4b). Although the changes were overall subtle (in part due to most of the genome staying the same), in TSS regions of the highly expressed genes, the relative occupancy of the TSS active state increased at the expense of the bivalent TSS, repressed by polycomb and poised enhancer states. Also, we saw that down-regulated DEGs diminished their relative enrichment of TSS bivalent whereas up-regulated DEGs showed a higher relative enrichment for poised enhancer. Then, we explored the proportion of different transitions among the 191 000 genomic regions that were differentially annotated between saline and LPS. Interestingly, among these dynamic genomic regions, the proportions of active TSS (8 to 10%), flanking active TSS (6 to 17%) and repressed by polycomb (7 to 13%) increased particularly at the expense of the enhancer (25 to 14%) and enhancer in gene (25 to 18%) states (Fig. 4c). Amongst these transitions, some occurred more than it would have produced solely by chance. For instance, the transitions from quiescent to represssed polycomb (8.3%, p_adj< 2.2e-16) and from enhancer to poised enhancer (7.6%, p_adj< 2.2e-16) were the most frequently observed transitions. In general terms, the significant transitions were related to modifications in the regions annotated with enhancers and poised enhancers (Fig. 4c and additional file 11). Finally, we checked the functions of the DEGs that overlapped the top 5 more frequently observed transitions (Fig. 4d, 4e, additional file: Fig S8 and additional file 12). For the transition from quiescent to represssed polycomb state, as expected, most of the functions were related to non-immune activities such as heart contraction (p_adj< 3.4e-05), blood circulation (p_adj< 3.4e-05) or circulatory system process (p_adj< 5.0e-05) (Fig. 4d and additional file 12). However, interestingly, other top significant term, regulation of system process(p_adj< 3.4e-05), included genes such as Il6, Pde4b and KDM4A, all of them upregulated after LPS exposure and involved in M1 polarization. This suggests that, not only non-necessary functions are repressed through chromatin states establishment. Also, given the differences in the evolution of the immune response in time, controlling the immune response after the initial stimulus is washed out is also achieved through PTHMs coordinated modifications. The second most significant transition, from enhancer to poised enhancer, overlapped with DEGs associated with cellular adaptation embedded in an inflammatory response endocytosis (p_adj< 9.7e-12), regulation of cell adhesion (p_adj< 1.8e-10) and cytokine production (p_adj< 1.8e-07). This reinforces the idea that the tunning and posterior anti-inflammatory control of the initially broad and unspecific response is executed, at least partly, through PTHMs modifications. Overall, our analysis suggests that the coordinated changes in chromatin states are associated with changes in expression. However, not only in a dichotomous expressed or non-expressed fashion but also in the functional role that a gene might play during acute immune response and depending on the stage of the response.

Monocytes and tissue-resident macrophages play important roles in development, metabolism, and immunity, thereby contributing to the maintenance of homeostasis [54]. In the brain, MCs are the first line of response to external irritants. Depending on the stimuli, MCs can adopt different phenotypes that are context-dependent [17]; however, an aberrant chronic polarization towards phenotypes related to inflammatory responses has been associated with several types of neural diseases [9]. Although the reprogramming of this phenotype is a complex process that results from an intricate relationship between several elements, epigenetic mechanisms are thought to play a major regulatory role [15, 17, 74, 75]. In this study, we assessed the relationship between epigenomic, and transcriptomic modifications endured by BV2 cells during an acute response induced by an external bacterial stimulus (LPS) exposure. Our analysis suggests that immunological cells activated by LPS undergo a reconfiguration of the epigenetic landscape (here partially depicted by the presence of four PTHMs) that is implemented in a coordinated fashion and that will have an important effect on the transcriptome. Additionally, our findings suggest a dynamic acquisition of the phenotype that shows mainly pro-inflammatory (M1) characteristics at the transcriptome level, but that exhibits epigenomic anti-inflammatory characteristics, suggesting a reversible reconfiguration after the initial stimulus.

Our results show important functional associations with changes in the distribution of PTHMs after LPS exposure. DBRs corresponding to each PTHM were annotated with distinct biological functions. The active marks changed abundance near immunological-related genes whereas the repressive mark was modified in regions related to organ development and vascular organization. Moreover, each one of the active PTHMs seemed to have distinct general functions in the inflammatory response. For instance, H3K4me3 seemed to increase in the vicinity of genes embedded in a broad innate immune response, whereas H3K27ac seemed to be more involved in the activation of, not only immune-related genes, but also in the regulation and tuning of specific pathways related to cytokine production, metabolism, and differentiation. This is consistent with other work showing that H3K27ac is particularly important in the early acquisition of phenotype [52, 73] in immune cells and that it correlates with gene expression (in enhancers and promoters) immediately, but also several hours after an initial LPS challenge [13, 55, 73]. On the other hand, H3k4me3 is correlated with expression 1 to 4 hours after the initial challenge, but decreases in promoters [55] and distal elements [13, 73] after that time window [55, 73]. Just as others [35, 55, 62], we observed that combinations of PTHMs do not occur in a random way. This suggests, as the “histone code” theory states [76], that combinations of PTHMs have unique functional repercussions on transcription. As previously demonstrated, we found a high co-occurrence of H3K4me3 and H3K27ac [73], H3K27ac and H3K4me1 [55, 77, 78], and H3K4me1 and H3K4me3 [35]. This led us to investigate how these combinatorial states changed after LPS induction or, in other words, how an initial immunological challenge reconfigures the epigenetic landscape in a coordinated fashion. We described 8 different chromatin states built from the combinations of the 4 marks studied. These states were similar to those described in similar cellular models [36, 79–81], which reinforces a pervasive “histone code” organization [76]. Our findings exhibit a reconfiguration affecting approximately 10% of the genome, which is in line with previous findings in monocytes/macrophages [73], especially at the expense of transitioning from states with low abundance of PTHMs to states rich in active marks [52, 57, 73]. This activation was previously described to be due mainly to changes in distal elements involving the active marks and particularly H3K27ac [52, 61, 73, 82]. To deepen our analysis of the immune induced reconfiguration, we sought to evaluate the transcriptomic profile of the challenged cells and couple it with the epigenetic changes. After LPS exposure, just like in other studies, a transition towards the M1 phenotype occurred [12, 17, 83]. Important immunological genes such as Il-6, Marco, Tnf-α were significantly up-regulated compared to control. Finally, we compared both transcriptomic and epigenomic profiles, showing that the different states and, moreover, the transitions from state to state induced by LPS, are associated with particular biological processes.

Interestingly, although the transcriptome clearly depicted an M1 phenotype, the epigenomic profile did not completely overlap it. Most observations were consistent between mRNA expression and PTHMs. Highly expressed genes were accompanied by high abundance of active marks in their promoters and up-regulated genes were associated with PTHM abundance gain in nearby regions. Conversely, genes lowly expressed were greatly bound by H3K27me3. Most of these paired changes illustrated a pro-inflammatory activity. Nonetheless, some regions near anti-inflammatory genes were annotated with active transcription chromatin states. We considered thus that we captured a BV2 cell classically activated (M1 phenotype) reverting towards a naïve-state epigenetic landscape. In this sense, only the innate immune memory component affected by LPS exposure would remain epigenetically modified [12, 74]. Indeed other studies have shown that BV2 cells exhibit transcriptomic changes within the first hours (1-4h), but once the challenge is removed the expression of those genes is dampen [74]. Moreover, a similar model showed a peak of immunological activity 12 hours after exposure followed by a complete reversion 48 hours after exposure [82]. Hence, we consider that, although, from a transcriptomic point of view, we captured a classical M1 phenotype, the epigenomic landscape depicted a phenotype of stable and transient changes induced by LPS but also characteristics of an anti-inflammatory regulatory M2 phenotype.

Given that in this study we evaluated MCs in a unique timepoint after exposure, it is crucial to consider the different combinations of PTHMs (a.k.a, the chromatin states) and their transitions between conditions to more accurately describe the true epigenomic dynamics. Our results show that this categorical evaluation of states (in specific timepoints) represents non-random functional units that have an impact on MCs. Although these states are cross-sectional measurements of a continuous and dynamic activity of each PTHM [52, 55, 73], they seem to accurately characterize transient and definitive modifications undergone by histones in a specific time point (24 hours in our case). Specifically, we found that the most frequent transitions induced by LPS reflected changes in clusters of DEGs involved in different functions. Independently of the label attributed to each state, the transitions seemed to represent activation/inactivation of different programs (inhibition of morphogenesis; changes in metabolism, homeostasis, and cytokine regulation; activation of the inflammatory response, etc.). This has already been observed in development [84], but to our knowledge this is the first time that such an analysis is performed in MCs in the context of inflammation. Indeed, our data show some unexpected changes in abundance of unique PTHM that, when studied in the context of states seem to gain consistency with the expected phenotypic dynamics of the model. The different states, which are illustrations of the previous activity of each PTHM can bring a better understanding of the impact of the epigenomic landscape in cellular behavior.

Besides the importance of our findings, there are important points to be raised in this study. For instance, we discuss many microglial behavior in the context of inflammation, but the suitability of BV2 cells as a surrogate of mouse MCs has been previously questioned [12, 85]. Thus, the generalizability of our results must be made with caution. Also, we studied 4 out of a myriad of possible PTHMs, therefore, these might change if more information is taking into account. Lastly, in this study we evaluated a single timepoint and estimated the impact between two groups, but a longitudinal study is warranted to truly determine the effect of LPS in MCs and the states that participate therein.

In this study we described the epigenomic adaptation that MCs undergo in an acute inflammatory context. Through the depiction of a phenotypic adaptation that is regulated by coordinated transitions taking place in different elements of the genome, we show the relevance of studying chromatin states rather than PTHMs in the context of neuroinflammation. Assessing states may reduce the historical ambiguity of the functional annotation attributed to each PTHMs (e.g., H3K4me1 as a memory mark [86] or H3K27ac as an activator [57]). Furthermore, we exhibit a basis to evaluate not only acute, but chronic effects of noxious external stimuli in cellular models. This is of particular importance in neurons and microglia, since deleterious changes in these cells derived from detrimental environments related to stress seem to be at the very basis of most neurological and psychiatric disorders [6, 87].

BH: Benjamini Hochberg

bp: Base pairs

CNS: Central nervous system

CUT&RUN: Cut Under Targets and Release Using Nuclease

DBA: Differential binding analysis

DBR: differentially bound regions

DEGs: Differentially expressed genes

FDR: False discovery rate

GLM: Generalized linear model

GO: Gene Ontology

HDACs: Histone deacetylases

IGV: Integrative Genomics Viewer

LPS: Lipopolysaccharide

MA: Motif analysis

ME: Motifs enrichment

MCs: Microglial cells

ORA: Over representation analysis

PCA: Principal component analysis

PCR: Polymerase chain reaction

PTHMs: Post-translational histone modifications

RIN: RNA integrity number

RNA-Seq: RNA sequencing

RPKM: Reads per kilo base per million mapped reads

TSS: Transcription start site

TF: Transcription factor

UCSC: University of California Santa Cruz

Ethics approval and consent to participate:

Not applicable

Consent for publication:

Not applicable

Availability of data and materials:

The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.

Competing interests

The authors declare that they have no competing interests.

Funding

BL holds a Sentinelle Nord Research Chair, is supported by the Canadian Institutes of Health

Research (Grant No. PJT-451728 and PJT-451858), and the Natural Science and Engineering

Research Council of Canada (Grant No. RGPIN-2019-06496) and receives Fonds de Recherche

en Santé du Québec (FRQS) Junior-2 salary support. AMR receives the doctoral scholarship from FRQS as well.

Authors' contributions

A.M.R, B.L, and C.Z designed the study. J.B and C.Z performed and optimized the experiments. A.M.R, A.M.P and S.M analyzed the data, created the figures and interpreted the results with support from B.L and A.D. A.M.R, A.D and B.L drafted the manuscript. All authors reviewed and approved the final manuscript.

Acknowledgements

The authors thank the funding sources for supporting this work.

Authors' information

¹CERVO Brain Research Center, Québec City, QC, Canada

²Department of Psychiatry and Neurosciences, Faculty of Medicine, Université Laval, Québec City, QC, Canada

³ Department of Social and Preventive Medicine, Faculty of Medicine, Université Laval, Québec City, QC, Canada

⁴ Genomics Center, Centre Hospitalier Universitaire de Québec—Université Laval Research Center, Québec City, QC, Canada

⁵ Department of Molecular Medicine, Faculty of Medicine, Université Laval, Québec City, QC, Canada

Chen L, Deng H, Cui H, Fang J, Zuo Z, Deng J, et al. Inflammatory responses and inflammation-associated diseases in organs [Internet]. Oncotarget. Impact Journals, LLC; 2018 [cited 2022 Jul 28]. p. 7204–18. Available from: /pmc/articles/PMC5805548/
Chaplin DD. Overview of the immune response. J Allergy Clin Immunol [Internet]. NIH Public Access; 2010 [cited 2022 Jul 28];125:S3. Available from: /pmc/articles/PMC2923430/
Bauer ME, Teixeira AL. Neuroinflammation in Mood Disorders: Role of Regulatory Immune Cells [Internet]. Neuroimmunomodulation. 2021 [cited 2022 Jul 28]. p. 99–107. Available from: www.karger.com/nim
Dunn GA, Loftis JM, Sullivan EL. Neuroinflammation in psychiatric disorders: An introductory primer. Pharmacol Biochem Behav. 2020;196.
Moncrieff J, Cooper RE, Stockmann T, Amendola S, Hengartner MP, Horowitz MA. The serotonin theory of depression: a systematic umbrella review of the evidence [Internet]. Mol. Psychiatry. 2022 [cited 2022 Jul 29]. Available from: http://www.ncbi.nlm.nih.gov/pubmed/35854107
Wang L, Yu C-C, Liu X-Y, Deng X-N, Tian Q, Du Y-J. Epigenetic Modulation of Microglia Function and Phenotypes in Neurodegenerative Diseases. Neural Plast [Internet]. 2021 [cited 2022 Jul 28]; Available from: https://doi.org/10.1155/2021/9912686
Stewart AG, Beart PM. Inflammation: Maladies, models, mechanisms and molecules [Internet]. Br. J. Pharmacol. Wiley-Blackwell; 2016 [cited 2022 Jul 29]. p. 631–4. Available from: /pmc/articles/PMC4742297/
Ransohoff RM, Schafer D, Vincent A, Blachère NE. Neuroinflammation: Ways in Which the Immune System Affects the Brain. Neurotherapeutics. 2015;12:896–909.
Neal M, Richardson JR. Epigenetic regulation of astrocyte function in neuroinflammation and neurodegeneration. Biochim. Biophys. Acta - Mol. Basis Dis. 2018. p. 432–43.
Yeh H, Ikezu T. Transcriptional and Epigenetic Regulation of Microglia in Health and Disease. Trends Mol. Med. 2019. p. 96–111.
Almeida PGC, Nani JV, Oses JP, Brietzke E, Hayashi MAF. Neuroinflammation and glial cell activation in mental disorders. Brain, Behav Immun - Heal. Elsevier BV; 2020;2:100034.
Das A, Hwa Kim S, Arifuzzaman S, Yoon T, Choul Chai J, Seek Lee Y, et al. Transcriptome sequencing reveals that LPS-triggered transcriptional responses in established microglia BV2 cell lines are poorly representative of primary microglia. 2016 [cited 2022 Aug 7]; Available from: http://genome.ucsc.edu/
Zhang X, Kracht L, Lerario AM, Dubbelaar ML, Brouwer N, Wesseling EM, et al. Epigenetic regulation of innate immune memory in microglia. J Neuroinflammation 2022 191 [Internet]. BioMed Central; 2022 [cited 2022 Aug 7];19:1–19. Available from: https://jneuroinflammation.biomedcentral.com/articles/10.1186/s12974-022-02463-5
Wu H, Zheng J, Xu S, Fang Y, Wu Y, Zeng J, et al. Mer regulates microglial/macrophage M1/M2 polarization and alleviates neuroinflammation following traumatic brain injury. J Neuroinflammation [Internet]. BioMed Central Ltd; 2021 [cited 2022 Aug 7];18:1–20. Available from: https://jneuroinflammation.biomedcentral.com/articles/10.1186/s12974-020-02041-7
Qiu M, Xu E, Zhan L. Epigenetic Regulations of Microglia/Macrophage Polarization in Ischemic Stroke. Front Mol Neurosci. Frontiers Media S.A.; 2021;14:199.
Yao B, Christian KM, He C, Jin P, Ming GL, Song H. Epigenetic mechanisms in neurogenesis [Internet]. Nat. Rev. Neurosci. 2016 [cited 2022 Aug 7]. p. 537–49. Available from: www.nature.com/nrn
Cheray M, Joseph B. Epigenetics control microglia plasticity. Front Cell Neurosci. Frontiers Media S.A.; 2018;12:243.
He Y, Yao X, Taylor N, Bai Y, Lovenberg T, Bhattacharya A. RNA sequencing analysis reveals quiescent microglia isolation methods from postnatal mouse brains and limitations of BV2 cells. J Neuroinflammation [Internet]. BioMed Central; 2018 [cited 2023 Mar 1];15. Available from: /pmc/articles/PMC5964710/
Skene PJ, Henikoff S. An efficient targeted nuclease strategy for high-resolution mapping of DNA binding sites. Elife. 2017;6:1–35.
Ewels P, Magnusson M, Lundin S, Käller M. MultiQC: Summarize analysis results for multiple tools and samples in a single report. Bioinformatics [Internet]. Oxford Academic; 2016 [cited 2023 Mar 2];32:3047–8. Available from: https://academic.oup.com/bioinformatics/article/32/19/3047/2196507
Bolger AM, Lohse M, Usadel B. Trimmomatic: A flexible trimmer for Illumina sequence data. Bioinformatics. Oxford University Press; 2014;30:2114–20.
Zhu Q, Liu N, Orkin SH, Yuan GC. CUT and RUNTools: A flexible pipeline for CUT and RUN processing and footprint analysis. Genome Biol [Internet]. BioMed Central Ltd.; 2019 [cited 2022 Aug 8];20:1–12. Available from: https://genomebiology.biomedcentral.com/articles/10.1186/s13059-019-1802-4
Yu F, Sankaran VG, Yuan GC. CUT&RUNTools 2.0: a pipeline for single-cell and bulk-level CUT&RUN and CUT&Tag data analysis. Bioinformatics [Internet]. Oxford University Press; 2022 [cited 2023 Mar 30];38:252. Available from: /pmc/articles/PMC8696090/
Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods 2012 94 [Internet]. Nature Publishing Group; 2012 [cited 2022 Aug 8];9:357–9. Available from: https://www.nature.com/articles/nmeth.1923
Langmead B, Wilks C, Antonescu V, Charles R. Scaling read aligners to hundreds of threads on general-purpose processors. Bioinformatics [Internet]. Oxford Academic; 2019 [cited 2022 Aug 8];35:421–32. Available from: https://academic.oup.com/bioinformatics/article/35/3/421/5055585
Zheng Y, Ahmad K, Henikoff S. CUT&Tag Data Processing and Analysis Tutorial [Internet]. protocols.io. 2020. Available from: dx.doi.org/10.17504/protocols.io.bjk2kkye
Suzuki R, Shimodaira H. Pvclust: An R package for assessing the uncertainty in hierarchical clustering. Bioinformatics [Internet]. 2006 [cited 2023 Mar 30];22:1540–2. Available from: http://www.is.
Kent WJ, Zweig AS, Barber G, Hinrichs AS, Karolchik D. BigWig and BigBed: enabling browsing of large distributed datasets. Bioinformatics [Internet]. Oxford University Press; 2010 [cited 2023 Mar 1];26:2204. Available from: /pmc/articles/PMC2922891/
Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics [Internet]. Oxford University Press; 2010 [cited 2023 Mar 1];26:841. Available from: /pmc/articles/PMC2832824/
Shen L, Shao N-Y, Liu X, Maze I, Feng J. diffReps: Detecting Differential Chromatin Modification Sites from ChIP-seq Data with Biological Replicates. PLoS One [Internet]. 2013;8:65598. Available from: www.plosone.org
Lawrence M, Huber W, Pagès H, Aboyoun P, Carlson M, Gentleman R, et al. Software for Computing and Annotating Genomic Ranges. PLOS Comput Biol [Internet]. Public Library of Science; 2013 [cited 2023 Apr 24];9:e1003118. Available from: https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1003118
Zhu LJ, Gazin C, Lawson ND, Pagès H, Lin SM, Lapointe DS, et al. ChIPpeakAnno: A Bioconductor package to annotate ChIP-seq and ChIP-chip data. BMC Bioinformatics [Internet]. BioMed Central; 2010 [cited 2023 Apr 24];11:1–10. Available from: https://bmcbioinformatics.biomedcentral.com/articles/10.1186/1471-2105-11-237
Conway JR, Lex A, Gehlenborg N. UpSetR: an R package for the visualization of intersecting sets and their properties. Bioinformatics [Internet]. Oxford Academic; 2017 [cited 2023 Apr 24];33:2938–40. Available from: https://academic.oup.com/bioinformatics/article/33/18/2938/3884387
Krassowski M. “ComplexUpset.” 2020 [cited 2023 May 24]; Available from: https://zenodo.org/record/7314197
Ernst J, Kellis M. Chromatin-state discovery and genome annotation with ChromHMM. Nat Protoc [Internet]. Nature Publishing Group; 2017 [cited 2022 Dec 23];12:2478–92. Available from: https://www-nature-com.acces.bibl.ulaval.ca/articles/nprot.2017.124
van der Velde A, Fan K, Tsuji J, Moore JE, Purcaro MJ, Pratt HE, et al. Annotation of chromatin states in 66 complete mouse epigenomes during development. Commun Biol 2021 41 [Internet]. Nature Publishing Group; 2021 [cited 2023 Mar 3];4:1–15. Available from: https://www.nature.com/articles/s42003-021-01756-4
Lutz PE, Chay MA, Pacis A, Chen GG, Aouabed Z, Maffioletti E, et al. Non-CG methylation and multiple histone profiles associate child abuse with immune and small GTPase dysregulation. Nat Commun 2021 121 [Internet]. Nature Publishing Group; 2021 [cited 2023 Mar 21];12:1–16. Available from: https://www.nature.com/articles/s41467-021-21365-3
Bailey TL, Boden M, Buske FA, Frith M, Grant CE, Clementi L, et al. MEME Suite: Tools for motif discovery and searching. Nucleic Acids Res. 2009;37.
Robinson JT, Thorvaldsdóttir H, Winckler W, Guttman M, Lander ES, Getz G, et al. Integrative genomics viewer. Nat Biotechnol [Internet]. 2011 [cited 2023 Apr 25];29:24–6. Available from: http://maq.sourceforge.net/
Bray NL, Pimentel H, Melsted P, Pachter L. Near-optimal probabilistic RNA-seq quantification. Nat Biotechnol 2016 345 [Internet]. Nature Publishing Group; 2016 [cited 2023 Apr 26];34:525–7. Available from: https://www.nature.com/articles/nbt.3519
Gallego Romero I, Pai AA, Tung J, Gilad Y. RNA-seq: Impact of RNA degradation on transcript quantification. BMC Biol [Internet]. BioMed Central Ltd.; 2014 [cited 2023 Apr 26];12:1–13. Available from: https://bmcbiol.biomedcentral.com/articles/10.1186/1741-7007-12-42
Leek JT. svaseq: removing batch effects and other unwanted noise from sequencing data. Nucleic Acids Res [Internet]. Oxford University Press; 2014 [cited 2023 Apr 26];42:e161. Available from: /pmc/articles/PMC4245966/
Yu L, Qu H, Jia Q, Wang X, Jia Z. Transformation, Normalization, and Batch Effect Removal. BIO-PROTOCOL. Bio-Protocol, LLC; 2022;12.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol [Internet]. BioMed Central Ltd.; 2014 [cited 2022 Oct 6];15:1–21. Available from: https://genomebiology.biomedcentral.com/articles/10.1186/s13059-014-0550-8
Yu G, Wang LG, Han Y, He QY. ClusterProfiler: An R package for comparing biological themes among gene clusters. Omi A J Integr Biol. 2012;16:284–7.
Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation [Internet]. Elsevier; 2021;2:100141. Available from: http://dx.doi.org/10.1016/j.xinn.2021.100141
Wickham H. ggplot2 Elegant Graphics for Data Analysis. Use R! Ser. Springer; 2016;211.
Blighe K, Rana S, Lewis M. EnhancedVolcano version 1.10.0: Publication-ready volcano plots with enhanced colouring and labeling [Internet]. R-Package. 2021 [cited 2023 Apr 26]. Available from: https://bioconductor.org/packages/release/bioc/html/EnhancedVolcano.html
Yu G. enrichplot: Visualization of Functional Enrichment Result. 2023 [cited 2023 Apr 26]; Available from: https://bioconductor.org/packages/release/bioc/html/enrichplot.html
Gao T, Qian J. EnhancerAtlas 2.0: an updated resource with enhancer annotation in 586 tissue/cell types across nine species. Nucleic Acids Res [Internet]. Oxford Academic; 2020 [cited 2023 May 24];48:D58–64. Available from: https://academic.oup.com/nar/article/48/D1/D58/5628925
Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010;11:R106.
Novakovic B, Habibi E, Wang SY, Arts RJW, Davar R, Megchelenbrink W, et al. β-Glucan Reverses the Epigenetic State of LPS-Induced Immunological Tolerance. Cell. Cell Press; 2016;167:1354-1368.e14.
Naler LB, Hsieh YP, Geng S, Zhou Z, Li L, Lu C. Epigenomic and transcriptomic analyses reveal differences between low-grade inflammation and severe exhaustion in LPS-challenged murine monocytes. Commun Biol 2022 51 [Internet]. Nature Publishing Group; 2022 [cited 2023 Apr 24];5:1–17. Available from: https://www.nature.com/articles/s42003-022-03035-2
Zhang X, Kracht L, Lerario AM, Dubbelaar ML, Brouwer N, Wesseling EM, et al. Epigenetic regulation of innate immune memory in microglia. J Neuroinflammation. BioMed Central Ltd; 2022;19.
Vandenbon A, Kumagai Y, Lin M, Suzuki Y, Nakai K. Waves of chromatin modifications in mouse dendritic cells in response to LPS stimulation. Genome Biol [Internet]. BioMed Central Ltd.; 2018 [cited 2023 May 28];19:1–15. Available from: https://genomebiology.biomedcentral.com/articles/10.1186/s13059-018-1524-z
Kitazawa T, Machlab D, Joshi O, Maiorano N, Kohler H, Ducret S, et al. A Unique Bipartite Polycomb Signature Regulates Stimulus-Response Transcription during Development. Nat Genet [Internet]. Europe PMC Funders; 2021 [cited 2023 May 25];53:379. Available from: /pmc/articles/PMC7610396/
Fujita Y, Yamashita T. Alterations in Chromatin Structure and Function in the Microglia. Front. Cell Dev. Biol. Frontiers Media S.A.; 2021. p. 1884.
Kreit M, Vertommen D, Gillet L, Michiels T. The Interferon-Inducible Mouse Apolipoprotein L9 and Prohibitins Cooperate to Restrict Theiler’s Virus Replication. PLoS One [Internet]. PLOS; 2015 [cited 2023 Mar 4];10. Available from: /pmc/articles/PMC4510265/
Ernst J, Kheradpour P, Mikkelsen TS, Shoresh N, Ward LD, Epstein CB, et al. Mapping and analysis of chromatin state dynamics in nine human cell types. 2011;
Seisen T, Roupret M, Gomez F, Malouf GG, Shariat SF, Peyronnet B, et al. A comprehensive review of genomic landscape, biomarkers and treatment sequencing in castration-resistant prostate cancer. Cancer Treat Rev [Internet]. Netherlands; 2016;48:25–33. Available from: http://ovidsp.ovid.com/ovidweb.cgi?T=JS&PAGE=reference&D=med12&NEWS=N&AN=27327958
Gosselin D, Skola D, Coufal NG, Holtman IR, Schlachetzki JCM, Sajti E, et al. An environment-dependent transcriptional network specifies human microglia identity. Science (80- ) [Internet]. NIH Public Access; 2017 [cited 2023 Apr 24];356:1248–59. Available from: /pmc/articles/PMC5858585/
Vu H, Ernst J. Universal annotation of the human genome through integration of over a thousand epigenomic datasets. Genome Biol [Internet]. BioMed Central Ltd; 2022 [cited 2023 Apr 24];23:1–37. Available from: https://genomebiology.biomedcentral.com/articles/10.1186/s13059-021-02572-z
Wu R, Chen F, Wang N, Tang D, Kang R. ACOD1 in immunometabolism and disease. Cell Mol Immunol [Internet]. Nature Publishing Group; 2020 [cited 2023 Mar 4];17:822. Available from: /pmc/articles/PMC7395145/
Pires-Afonso Y, Muller A, Grzyb K, Oudin A, Yabo YA, Sousa C, et al. Elucidating tumour‐associated microglia/macrophage diversity along glioblastoma progression and under ACOD1 deficiency. Mol Oncol [Internet]. Wiley-Blackwell; 2022 [cited 2023 Jun 5];16:3167. Available from: /pmc/articles/PMC9441003/
Wei L, Du Y, Xie Y, Yu X, Chen H, Qiu Y. Lipocalin-2 Regulates Hippocampal Microglial Activation in Poststroke Depression. Front Aging Neurosci [Internet]. Frontiers Media S.A.; 2021 [cited 2023 Jun 5];13. Available from: /pmc/articles/PMC8710735/
Jang E, Lee S, Kim JH, Kim JH, Seo JW, Lee WH, et al. Secreted protein lipocalin-2 promotes microglial M1 polarization. FASEB J [Internet]. 2013 [cited 2023 Jun 5];27:1176–90. Available from: https://faseb.onlinelibrary.wiley.com/doi/10.1096/fj.12-222257
Gao H, Danzi MC, Choi CS, Taherian M, Dalby-Hansen C, Ellman DG, et al. Opposing functions of microglial and macrophagic TNFR2 in the pathogenesis of experimental autoimmune encephalomyelitis. Cell Rep [Internet]. NIH Public Access; 2017 [cited 2023 Jun 5];18:198. Available from: /pmc/articles/PMC5218601/
Jurga AM, Paleczna M, Kuter KZ. Overview of General and Discriminating Markers of Differential Microglia Phenotypes. Front Cell Neurosci. Frontiers Media S.A.; 2020;14:198.
Cai Y, Zhang Y, Loh YP, Tng JQ, Lim MC, Cao Z, et al. H3K27me3-rich genomic regions can function as silencers to repress gene expression via chromatin interactions. Nat Commun [Internet]. Nature Publishing Group; 2021 [cited 2023 Jun 8];12. Available from: /pmc/articles/PMC7846766/
Ishii M, Wen H, Corsa CAS, Liu T, Coelho AL, Allen RM, et al. Epigenetic regulation of the alternatively activated macrophage phenotype. Blood [Internet]. American Society of Hematology; 2009 [cited 2023 Mar 5];114:3244–54. Available from: https://ashpublications.org/blood/article/114/15/3244/125102/Epigenetic-regulation-of-the-alternatively
Marrocco A, Ortiz LA. Role of metabolic reprogramming in pro-inflammatory cytokine secretion from LPS or silica-activated macrophages. Front. Immunol. Frontiers Media S.A.; 2022. p. 936167.
Zhao S, Ye Z, Stanton R. Misuse of RPKM or TPM normalization when comparing across samples and sequencing protocols. RNA [Internet]. Cold Spring Harbor Laboratory Press; 2020 [cited 2023 Jul 24];26:903–9. Available from: /pmc/articles/PMC7373998/
Saeed S, Quintin J, Kerstens HHD, Rao NA, Aghajanirefah A, Matarese F, et al. Epigenetic programming of monocyte-to-macrophage differentiation and trained innate immunity. Science (80- ) [Internet]. NIH Public Access; 2014 [cited 2023 Mar 5];345:1251086. Available from: /pmc/articles/PMC4242194/
Das A, Chai JC, Kim SH, Park KS, Lee YS, Jung KH, et al. Dual RNA sequencing reveals the expression of unique transcriptomic signatures in lipopolysaccharide-induced BV-2 microglial cells. PLoS One. 2015;10.
Jiang W, Zhang L-X, Tan X-Y, Yu P, Dong M. Inflammation and histone modification in chronic pain. Front Immunol [Internet]. Frontiers Media SA; 2023 [cited 2023 Feb 28];13. Available from: /pmc/articles/PMC9880030/
Bonitto K, Sarathy K, Atai K, Mitra M, Coller HA. Is There a Histone Code for Cellular Quiescence? Front Cell Dev Biol [Internet]. Frontiers Media SA; 2021 [cited 2023 May 28];9:739780. Available from: /pmc/articles/PMC8586460/
Creyghton MP, Cheng AW, Welstead GG, Kooistra T, Carey BW, Steine EJ, et al. Histone H3K27ac separates active from poised enhancers and predicts developmental state. Proc Natl Acad Sci U S A [Internet]. 2010 [cited 2023 Apr 25];107:21931–6. Available from: www.pnas.org/cgi/doi/10.1073/pnas.1016071107
Beacon TH, Delcuve GP, López C, Nardocci G, Kovalchuk I, van Wijnen AJ, et al. The dynamic broad epigenetic (H3K4me3, H3K27ac) domain as a mark of essential genes [Internet]. Clin. Epigenetics. BioMed Central; 2021 [cited 2023 May 24]. p. 1–17. Available from: https://clinicalepigeneticsjournal.biomedcentral.com/articles/10.1186/s13148-021-01126-1
Forouzanfar MH, Afshin A, Alexander LT, Anderson HR, Bhutta ZA, Biryukov S, et al. Global, regional, and national comparative risk assessment of 79 behavioural, environmental and occupational, and metabolic risks or clusters of risks, 1990–2015: a systematic analysis for the Global Burden of Disease Study 2015. Lancet. 2016;388:1659–724.
Ernst J, Kellis M. ChromHMM: automating chromatin state discovery and characterization. [cited 2023 Apr 24]; Available from: http://compbio.mit.edu/ChromHMM.
Hoffman MM, Ernst J, Wilder SP, Kundaje A, Harris RS, Libbrecht M, et al. Integrative annotation of chromatin elements from ENCODE data. Nucleic Acids Res [Internet]. Oxford Academic; 2013 [cited 2023 Apr 28];41:827–41. Available from: https://academic.oup.com/nar/article/41/2/827/1071531
Huang M, Malovic E, Ealy A, Jin H, Anantharam V, Kanthasamy A, et al. Microglial immune regulation by epigenetic reprogramming through histone H3K27 acetylation in neuroinflammation. Front Immunol. NLM (Medline); 2023;14:1052925.
Lang GP, Li C, Han YY. Rutin pretreatment promotes microglial M1 to M2 phenotype polarization. Neural Regen Res [Internet]. Wolters Kluwer -- Medknow Publications; 2021 [cited 2023 Aug 7];16:2499. Available from: /pmc/articles/PMC8374565/
van der Velde A, Fan K, Tsuji J, Moore JE, Purcaro MJ, Pratt HE, et al. Annotation of chromatin states in 66 complete mouse epigenomes during development. Commun Biol 2021 41 [Internet]. Nature Publishing Group; 2021 [cited 2023 Apr 24];4:1–15. Available from: https://www.nature.com/articles/s42003-021-01756-4
Henn A, Lund S, Hedtjärn M, Schrattenholz A, Pörzgen P, Leist M. The suitability of BV2 cells as alternative model system for primary microglia cultures or for animal experiments examining brain inflammation. ALTEX. 2009. p. 83–94.
Ostuni R, Piccolo V, Barozzi I, Polletti S, Termanini A, Bonifacio S, et al. Latent enhancers activated by stimulation in differentiated cells. Cell [Internet]. Elsevier B.V.; 2013 [cited 2023 Aug 6];152:157–71. Available from: http://www.cell.com/article/S0092867412015401/fulltext
Guo S, Wang H, Yin Y. Microglia Polarization From M1 to M2 in Neurodegenerative Diseases [Internet]. Front. Aging Neurosci. Frontiers Media S.A.; 2022 [cited 2022 Aug 7]. Available from: /pmc/articles/PMC8888930/

No competing interests reported.

Download PDF

Version 1

posted

You are reading this latest preprint version

Immune-related transcriptomic and epigenetic reconfiguration in BV2 cells after LPS exposure: an in vitro omics integrative study

Status:

Version 1

Abstract

Background

Methods

Results

Conclusion

Figures

Background

Materials and methods

BV2 cells culture

BV2 cells stimulation and RNA extraction

Assessment of PTHMs distribution through CUT&RUN

Library preparation and sequencing

Bioinformatic analyses

Statistical analyses

Results

Integration of epigenetic modifications and expression

Discussion

Conclusion

Abbreviations

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1