Multi-Omic Epigenetic-Based Model Reveals Key Molecular Mechanisms Associated with Palmitic Acid Lipotoxicity in Human Astrocyte

Abstract

Astrocytes are critical for the metabolic, structural and functional modulatory support of the brain. Lipotoxicity or high levels of saturated fatty acid as Palmitate (PA) has been associated with neurotoxicity, the loss or change of astrocytic functionality, and the etiology and progression of neurodegenerative diseases such as Parkinson or Alzheimer. Several molecular mechanisms of PA’s effect in astrocytes have been described, yet the role of epigenetic regulation and chromatin architecture have not been fully explored. In this study, we developed a multi-omic epigenetic-based model to identify the molecular mechanisms of lipotoxic PA activity in astrocytes. We used data from nine histone modifications, location of Topological Associated Domains (TADs) and transcriptional CTCF regions, where we identified the basal astrocyte epigenetic landscape. Moreover, we integrated transcriptomic data of astrocytic cellular response to PA with the epigenetic multi-omic model to identify lipotoxic-induced molecular mechanisms. The multi-omic model showed that chromatin conformation in astrocytes treated with PA have response genes located within shared topological domains, in which most of them also showed either repressive or enhancing marks in the Chip-Seq enrichment, reinforcing the idea that epigenetic regulation has a huge impact on the lipotoxic mechanisms of PA in the brain.

1. Introduction

Obesity is referred to as the excessive accumulation of body fat. It has become a worldwide public health issue which several studies have linked hormonal impairment to other diseases like coronary pathologies, diabetes, hypertension, atherosclerosis, and certain types of cancer among others [1, 2]. Studies using insulin growth factor-1 (IGF-1) receptor, insulin receptor substrate-1 (IRS-1), insulin receptor substrate-4 (IRS-4), glial fibrillary acidic protein, as well as an increase in β-actin protein have been associated with fatty acid excess in the brain [3]. Additionally, recent evidence has linked adiposity and high fatty acid concentrations to significant brain region-specific dysfunction, atrophy, inflammation, and cognitive decline [2, 4, 5], as well as an increased risk in developing the accumulation of amyloid β and Tau associated with Alzheimer’s disease [3].

Astrocytes are the most versatile glial cells in the central nervous system (CNS) constituting from 20 to 40% of neuroglia, protecting the brain through so many signaling [6], demonstrating that these cells effectively engulf dead cells, synapses and protein aggregates of amyloid β (Aβ) and ɑ-synuclein, typical of Alzheimer’s disease (AD) and Parkinson’s disease (PD), respectively. Additionally, astrocytes have been shown to regulate K+ levels [7] and prevent excitotoxicity in Huntington’s disease (HD) [8, 9, 10]. Nonetheless, evidence suggests that elevated concentrations of fatty acids can trigger a pro-inflammatory response altering the correct functioning of astrocytes [11, 12, 13]. Recently, authors have proved that metabolic insults produced by fatty acids can trigger a pro-inflammatory response in astrocytes, due to their high recruitment and metabolic capacity. Among them is PA, a long-chain saturated fatty acid, that can trigger an increase in inflammatory cytokines [5] such as Interleukin (IL)-1B, IL-6 and tumor necrosis factor alpha (TNFα), leading to accelerated cognitive decline, decreased cell viability, increased endoplasmic reticulum stress, inhibition of autophagy, finally compromising the Blood–Brain Barrier (BBB) integrity and promoting dementia-like progression in humans and animal models [5, 7, 14, 15].

Recent evidence supports epigenetic responses in astrocytes followed by PA-lipotoxic exposure [10, 16]. Epigenetic transcriptional regulation such as chromatin accessibility by histone modifications and chromatin architecture modulate euchromatin/heterochromatin equilibrium has shown the great potential of providing groundbreaking insight into the effects of neurotoxic compounds such as PA [4, 17, 18]. Additionally, the epigenetic modulation in astrocytes produced by lipotoxic compounds like PA can trigger inflammation, neurotoxicity, astrocyte reactivity, and cell fate determination in the CNS [16, 19]. In this case, the epigenetic landscape regulatory role and its response in the PA-induced astrocyte lipotoxicity are both the key to comprehend the loss of cellular function. Furthermore, several authors have also demonstrated that multi-omic models have proved to be more efficient than conventional astrocytic models in the evaluation of non-linearity in chromatin regulation considering regulatory mechanisms such as enhancers, isolators, epigenetic marks, and non-coding RNA [20, 21, 22].

It has been demonstrated that epigenetic data such as Chip-Seq and Hi-C with transcriptomics allows the detailed identification of specific molecular mechanisms associated with impairment conditions. In the present study, we report a multi-omic model to describe the epigenetic baseline of astrocytes as well as the astrocytic response to PA-lipotoxicity over specific astrocytic processes such as inflammation, autophagy, endoplasmic reticulum stress, energetic metabolism, mitochondrial dysfunction, and astrocyte-neuron interaction pathways, herein described here as astrocytic PA response (APAR) mechanisms.

2. Materials and methods

2.1 Hi-C, ChiP-Seq and transcriptomics datasets acquisition

Hi-C has been adopted as a method to measure pairwise contacts between pairs of genomic loci and allows a mapping of the three-dimensional conformation of chromatin within a population of cells, as well as to detect the structural variation and corrects assembly of chromosomal missed junctions [23]. Chip-Seq data also allows the analysis of histone marks interaction with DNA in an activation/repression mechanism. In this study, we analyzed nine treatments which were controls, H3K36me3, H3K27me3, H3K9ac, H3K9me3, H3K79me2, H4K20me1, H3K4me1 and CTCF (entry: GSM733678, GSM733751, GSM733729, GSM1003534, GSM1003491, GSM1003490, GSM1003525, GSM733710 and GSM733765 respectively). Tissue-specific datasets for astrocyte Hi-C from cerebellum and spinal cord were downloaded from the Encyclopedia of DNA Elements (ENCODE), as part of the ENCODE project consortium with ID numbers 200105194 and 200105957, respectively [24, 25, 26]. From ChiP-Seq data we obtained nine astrocyte datasets from NHA cells culture from the ENCODE database. Moreover, the whole human genome GRC version hg19 was obtained from ENSEMBL (https://www.ensembl.org/index.html) to map and enrich all the datasets. Transcriptomic data was experimentally obtained in the laboratory of Experimental and Computational Biochemistry of the Pontificia Universidad Javeriana, Bogotá D.C, Colombia.

2.2 Transcriptomic data

We used Normal Human Astrocytes (NHA; Lonza, CC-2565) divided in three different batches (#0000612736, #00005656712, #0000514417), which were cultured according to the manufacturer’s specifications. All batches were cultured in a supplemented medium with SingleQuots supplements. In order to induce PA toxicity, NHA cells were seeded in 48, 24, 12, and 6-well plates and incubated in a humidified incubator for 12 days at 37°C and 5% CO2. Then the NHA cells were washed with PBS 1x and starved in medium with serum-free DMEM without phenol red, L- and supplements (Lonza, Basel, Switzerland) for 6 h.

RNA extraction was performed using mini kit RNeasy (Qiagen, USA). The RNA quantification of the preparations was performed with NanoDrop (Thermo Fisher Scientific, Waltham, Ma, 174 USA). To remove possible DNA contamination, RNA was treated with DNase I. The RNA integrity (RIN) and 28S/18S ratio were determined with the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA). Transcriptomic datasets were obtained for NHA astrocytes treated in DMEM medium. RNA-seq libraries were prepared using the TruSeq Stranded mRNA library prep kit following the manufacturer’s protocol (Illumina, Cat# RS-122-2101) [27, 28]. The RNA-seq libraries were sequenced in HiSeq platform (Illumina) using protocol 2x150bp paired-end configuration, single index per lane. Scores and nucleotide composition were assessed with FastQ to evaluate accuracy using the Nextflow (V18.10.1) pipeline QUARS (https://github.com/TainVelasco-Luquez/QUARS).

Salmon package was used for mapping and quantifying the expression level (https://combine-lab.github.io/salmon/, V0.13.1) on the genome assembly GRCh38 (patch 12) from the NCBI without ALT regions using Gencode [26, 29]. NOISeq was used to import data into R (V3.6.1) and assessed sequence plot quality diagnosis [5]. Gene level Ensembl IDs were used with tximport function to create the count matrix. We used DESeq2 for modeling the average expression in function of the treatments correcting for sex (design formula: ~ sex + condition) [30]. Moreover, DESeq2 was used for normalization by size, variance shrinkage, outliers filtering, and hypothesis testing. The Wald test was used for assessing genes differentially expressed above |LFC| > = 0.5 with an alpha cutoff at 0.05.

The overlap analysis was performed using Fisher’s exact test using alpha at 0.05 implemented in the package GeneOverlap. To correct for multiple testing, p-values were adjusted using the Benjamini-Hochberg method [31]. Additionally, we leveraged the GeneOverlap’s Odds Ratio and Jaccard Index as measures of the strength of association and the similarity, respectively. Odds ratio equal to or less than 1 means no associations and greater than 1 represents strong associations. Jaccard Index is a measure of similarity that varies between 0, no similarity, and 1, completely identical lists. The package WGCNA from the platform iDEP (V0.90) was used to perform the co-expression analysis, employing the 1000 more variable genes across all samples, with a soft threshold of 16 and minimum module size of 20 [32]. The Pearson’s correlation was used on the count matrix, normalized and regularized using the log transformation of the DESeq2 library.

2.3 Epigenetics data analysis

To create the multi-layer model, first we obtained the ChIP-Seq data of an astrocyte in homeostatic conditions from the ENCODE database, then re-analyzed the BED/BAM files using ChIP-Seq model-based analysis implemented in MACS2. We integrated H3K4me3, H3K27ac, H3K27me3, H3K9ac, H3K4me1, H4K20me1, H3K36me3, and H3K79me2 to the human genome, in order to identify the core active regulatory and repressed genes in astrocyte, H3K27me3, H4K20me1 and H3K9me3. The active genes were also identified by the integration of H3K4me3, H3K4me1, H3K79me, 2H3K9ac, H3K27me3, and H3K27ac [33] to the same genome (hg19) and those who shared both repressing and enhancing modifications were identified as bivalent genes. All the individual samples, the core activation, repression and bivalent samples were enriched using the cutoff p-value set at 0.01 for both molecular function and biological process excluding redundant gene ontology (GO) terms.

Hi-C data was obtained from the ENCODE database. In this database, it can be found up to 80% of the annotated genome, in which for the interest of the investigation cerebellum and spinal cord data with ID numbers 200105194 and 200105957 respectively [24]. Hi-C data of spinal cord and cerebellum were compared to identify the potential tissue-specific differences in astrocyte functioning. All the individual samples, the core activation, repression, and bivalent samples were enriched using the cutoff p value set at 0.01. Molecular enrichment was performed using ShinyGO [34].

2.4 Data integration

In order to identify euchromatin and heterochromatin regions in astrocyte, we overlapped all activation/repression specific histone modifications covering regions. With this approach, it was possible to identify activation, repression and bivalent core genomic regions. Thus, both the omic and epigenetic integrations were performed through the adjudication of the data described above into a multi-omic model. The model consisted in three different layers of the Chip-seq and Hi-C of an astrocyte under homeostatic conditions and the transcriptomic data of an astrocyte under the lipotoxic effects of PA, where these three layers were used to make an inference about the possible epigenetic effects of PA in an astrocyte. Considering that there are no Hi-C or Chip-Seq data for astrocytes under the effects of PA, the integration of the transcriptomic data allows the identification and analysis of genes associated with PA response. Accordingly, the identification of changes in the epigenetic regulation of genes was performed as follows: first the differentially expressed genes in the transcriptome presented within the TADs of euchromatin and with the groups of genes were identified. Then, the proteomic data was sought to identify whether gene expression patterns are correlated with histone modification data sets [35, 36, 37]. Also, both the core regions and the specific histone marks with the Hi-C were overlapped with the TADs of an astrocyte in order to identify patterns between chromatin architecture and modulation of histone expression [37, 38]. Later, the Chip-Seq core regions integrated with TADs were compared with gene expression and proteomic data.

To identify the role of chromatin conformation in the expression and the effect of TADs in PA activity, gene expression and gene localization were associated with each other. This approach allowed us to identify additional epigenetic regulatory events related with TAD genes [39]. Subsequently non-coding regions such as enhancers and promoter regions were identified to be able to explain the patterns of expression of PA lipotoxicity. All data analysis was developed using R-Bioconductor suite (https://www.bioconductor.org/) as well as publicly available databases to ensure reproducibility and robustness. All the resulting sets of genes were enriched for molecular processes and biological functions.

3. Results

3.1 Chromatin/Histone expression regulation

Considering the functional importance of histone modifications in the cellular behavior [40], we identified a set of 34852 genes with known activation roles across the seven ChiP-Seq samples. Histone marks H3K4me3, H3K27ac, H3K9ac, H3K4me1, H4K20me1, H3K36me3, and H3K79me2 and the number of genes per activation mark were 2432, 3034, 3773, 5133, 8228, 6766 and 5486, respectively, with a non-homogenous pattern. We obtained a set of 11214 genes enhanced between the 7 studied active histone modifications with at least each gene included in two or more of the samples. Moreover, samples H3K27me3 and H3K9me3 showed 6276 and 6747 repressed genes, respectively. We identified a set of 9796 genes repressed in astrocytes based on H3K27me3 and H3K9me3 data with the condition that each gene should be present in at least one of the ChiP samples. Considering that a bivalent region is due to the presence of a repressor and an enhancer in the modifications of histones H3K4me3 and H3K27me3, we identified 7608 genes in bivalent sites. Moreover, shared genes for both activation and repression were identified as coding regions present in at least one of the mark datasets in each group for the specific markers [33].

In general, we performed a functional enrichment of the dataset where it was possible to identify the biological process and molecular functions associated with GO terms. As a result, 30 biological processes relevant to APAR biological mechanisms were presented (Table 1). Among these biological processes identified for the activation of ChIP-Seq datasets, glutamate-cysteine ligase activity, CD4 receptor binding, ion channel binding, and extracellular matrix binding were the top-enriched functions. Besides, functional groups associated with the homeostatic astrocytic activity were identified as highlighting transferase activity, carbohydrate derivative binding, hydrolase activity, DNA-binding transcription factor activity, molecular function regulator, transporter activity, oxidoreductase activity, enzyme regulator activity, transmembrane transporter activity, signaling receptor activity, extracellular matrix structural constituent, lipid binding, extracellular matrix binding, electron transfer activity, and lyase/ligase activity. The genes that encode the molecular functions mentioned above, are associated with metabolic support in astrocytic activity and neural functionality present in euchromatin regions.

Category	Process	Gene number
Function	Regulation of biological quality	4319
	Homeostatic processes	2004
	Ion homeostasis	836
	Response to nutrients	730
	Regulation of membrane potential	450
	Apoptotic mitochondrial changes	130
	Membrane depolarization	88
	Regulation of membrane depolarization	44
	Mitochondrial depolarization	24
	Regulation of mitochondrial depolarization	21
	Non-ribosomal peptide biosynthetic process	19
	Glutathione biosynthetic process	17
Group	Immune response	45
	Regulation of biological quality	32
	Response to stress	23
	Regulation of response to stimulus	20
	Regulation of molecular function	20
	Response to external stimulus	17
	Regulation of signaling	15
	Cell adhesion	11
	Catabolic processes	11

Table 1.

Top biological processes associated with the core activation dataset from the histone markers H3K4me3, H3K27ac, H3K9ac, H3K4me1, H4K20me1, H3K36me3, and H3K79me2. Accordingly, biological processes were separated into functions and groups, in terms associated with GO for a more detailed analysis.

It was also possible to identify active coding regions tightly regulated for cellular ion maintenance and response to stimuli that are essential for astrocytes well-functioning [7]. Additionally, we were able to associate the presence of constant euchromatin regions with genes that encode for metabolic and cellular exchange mechanisms necessary for astrocyte function [41]. Therefore, the model demonstrated that the presence of genes in regular euchromatin regions are often associated with many regulatory elements such as promoters, enhancers, insulators and silencers, all related with cell adhesion, support and exchange processes [42, 43].

In the case of bivalent expression regions, 30 biological processes were identified corresponding with APAR mechanisms, highlighting transcription regulation, sequence-specific DNA binding, RNA polymerase II/III distal enhancer, regulatory region, and proximal promoter sequence-specific DNA binding, gamma-aminobutyric acid transmembrane transporter activity, nerve growth factor binding, ubiquitin-protein transferase activator activity, mannosyl-transferase activity and cofactor, corepressor and coactivator transcription binding (Table 2). Additionally, these same genes that were also associated with specific functional groups such as macromolecule binding (i.e., carbohydrates, sulfurates, lipids, amides), DNA-binding transcription factor activity, cofactor binding, extracellular matrix structural constituent, structural constituent of ribosome, structural constituent of cytoskeleton, extracellular matrix binding, neurotransmitter binding and structural constituent of myelin sheath and activity of hydrolase, transferase, peroxidase, oxidoreductase, isomerase, lyase/ligase signaling, transmembrane transporters, enzyme regulation, antioxidant, electron transfer, neurotransmitter transport, cytochrome-c oxidase, MAPKK, and glutathione dehydrogenase functionality.

Category	Process	Gene number
Function	Positive regulation of biosynthetic process	4091
	Positive regulation of nucleic acid-templated transcription	3293
	Positive regulation of cellular biosynthetic process	2082
	Positive regulation of RNA metabolic process	1789
	Positive regulation of RNA biosynthetic process	1695
	Extracellular matrix assembly	28
	Negative regulation of autophagy of mitochondrion	9
	Interleukin-23-mediated signaling pathway	9
	Positive regulation of axon extension involved in axon guidance	7
Group	Regulation of response to stimulus	1409
	Regulation of biological quality	1360
	Response to stress	1312
	Regulation of signaling	1242
	Immune system process	973
	Catabolic process	863
	Response to external stimulus	793
	Cell proliferation	708

Table 2.

Biological processes associated with histone markers H3K4me3 and H3K27me3 were separated between functions and groups. All the terms were associated with GO terms for further analysis [33].

Further, we identified all the specific genes associated with the APAR mechanisms in astrocytes under homeostatic conditions (Figure 1). It was possible to clarify the epigenetic basal response of astrocytes and identify the gene activation/expression profiling under homeostatic conditions of these cells to elucidate the potential response to PA detrimental conditions. Both activation and bivalent processes and functions made it possible to understand the expression and regulatory mechanisms associated with an epigenetic chromatin landscape in astrocytes [4]. For both activation and bivalent datasets, the biological processes and functions were consistent with the basal astrocytic activity (Tables 1 and 2).

Figure 1.

Graphical representation of activation, repression, and bivalent genes present in every APAR categories. Specifically, AN interaction has 42%, 25% and 31%; autophagy has 42%, 20% and 37%; ER Stress has 33%, 23% and 43%; inflammation has 39%, 28% and 31% all for activation, bivalent and repression genes respectively. Note: the colors correspond to Green-Activation, Orange-Bivalent, and Purple-Repression.

3.2 Chromatin architecture involvement in APAR

Overexpressed genes in NHA astrocytes exposed to lipotoxic PA concentrations were fully integrated to ChIP-Seq data, APAR-related data and the chromatin conformation in order to identify co-regulated genomic regions in astrocytes [44]. We identified 328 molecular processes that were found overexpressed when astrocytes are exposed to PA, of which 27 molecular APAR associated processes were selected based on the role and significance (Table 3).

Molecular process	Adjusted p-value	Negative log10 adjusted p-value	Gene number
Inflammatory response	1.0487E-16	15.979344	53
Response to external stimulus	9.8873E-13	12.0049202	101
Response to lipid	1.8892E-11	10.7237215	51
Response to stress	8.5056E-11	10.070294	121
Response to stimulus	3.438E-09	8.46369645	205
Regulation of lipid metabolic process	4.1334E-08	7.38369637	29
Cellular response to stimulus	1.0314E-05	4.98659017	169
Regulation of biological quality	1.8704E-05	4.72806228	107
Regulation of immune response	0.0001047	3.98004731	43
Regulation of cell activation	0.00018219	3.73947057	30
Positive regulation of inflammatory response	0.00029935	3.52381953	14
Regulation of inflammatory response	0.00057697	3.23884686	22
Regulation of response to stress	0.00213314	2.67098006	50
Cellular lipid metabolic process	0.00528208	2.27719528	38
Positive regulation of biological process	0.00848563	2.07131588	136
Neuroinflammatory response	0.0114169	1.94245175	9
Regulation of ERK1 and ERK2 cascade	0.01342102	1.8722145	17
Glial cell activation	0.01563567	1.8058835	8
Interleukin-1 secretion	0.01762805	1.75379577	8
Positive regulation of metabolic process	0.01849939	1.7328425	89
Regulation of response to stimulus	0.02555143	1.59258475	101
ERK1 and ERK2 cascade	0.0278352	1.55540561	17
Positive regulation of immune response	0.03214889	1.49283399	32
Hippocampal neuron apoptotic process	0.0381717	1.41825855	3
Regulation of hippocampal neuron apoptotic process	0.0381717	1.41825855	3
Negative regulation of transport	0.04292557	1.36728396	23
Synapse pruning	0.04548992	1.34208479	4

Table 3.

Differentially expressed biological processes associated with the APAR mechanisms of PA-lipotoxicity in astrocytes. All considered processes have p > 0.05 as a threshold value of significance.

Although comparing specific TAD regions from Hi-C experiments and differentially expressed genes associated with APAR mechanisms in astrocytes, we identified clusters sharing the same TAD (Table A1). In this regard, we identified 3039 and 3048 TAD regions for spinal cord and cerebellum astrocytes, respectively. We focused on differentially expressed genes present among the APAR gene sets, located in the corresponding TAD regions to identify co-regulated genes or regulatory profiles. Moreover, due to the CTCF role in the conformation of chromatin folding architecture, TAD regions were overlapped with CTCF in order to identify true TAD regions. In this sense, these results elucidated some of the potential role of epigenetic modulation in the APAR molecular mechanisms in astrocytes in response to PA-lipotoxicity [45, 46, 47].

4. Discussion

In the present study, our model showed that astrocytes regulate enzymatic and protein activity from the genomic to the protein level, considering the protein functional modulation at different molecular levels. Additionally, during normal conditions, the enzymatic activity of hydrolase, transferase, peroxidase, oxidoreductase, isomerase and lyase/ligase were identified as constantly regulated due to elevated metabolic rates and plasticity in astrocytes [48, 49]. In this case, metabolic maintenance and support are not permanently regulated by epigenetic processes due to a dynamic environmental-dependent mechanism in astrocytes. Nevertheless, the presence of metabolic processes in bivalent regions implies the presence of highly active metabolic processes that change across time due to the fact that a genomic region can present both marks and become active or repressed [43, 50]. In terms of astrocyte-neuron interaction, we identified the presence of antioxidant activity associated with glutathione biosynthetic processes, reductase activity, as well as the activity of structural constituents of myelin sheath [51, 52]. Relationship between astrocytes and neurons in the context of antioxidant defense to ensuring neuronal well-being during pathological conditions play a significant role in metabolic support by neuroprotective capacity from oxidative stress, supply of glutathione to neurons, modulation of the extracellular matrix assembly, among others [52, 53].

To examine the molecular response to PA or APAR mechanism in astrocytes, we integrated the epigenetic data with the transcriptomic data from NHA to elucidate the potential damaging conditions by the PA activity in the brain. The shared TAD regions from both cerebellum and spinal cord astrocyte Hi-C data were compared to each other in order to establish the differences and possible considerations associated with tissue-specific stimuli. Thus, our multi-omic model showed that during PA lipotoxicity in astrocytes, inflammatory and stress responses are overexpressed. Our results also indicated that lipid droplets are epigenetically regulated in order to respond to free fatty acid concentrations in homeostatic conditions by the presence of apolipoprotein-E (APOE) gene in euchromatin regions [4, 40]. For instance, recent evidence has shown that maintenance of the homeostasis between astrocytes and neurons mitigate the lipotoxic effects of fatty acids as well as modulating APOE-lipid particles becomes of vital importance [54].

The presence of PA is associated with the overexpression of biological processes such as response to cellular lipid metabolism, which can lead to disease [5, 55]. Moreover, high concentrations of PA induce the expression of markers involved in pro-inflammatory response where the secretion of IL-1 activates endothelial cells and astrocytes to propagate the inflammatory signals in CNS [56, 57]. Overall, IL-1 is a typical biomarker associated with lipotoxicity and inflammation in astrocytes, as LC3-II, p62, or TLR2 have been directly linked to the astrocytic response to PA [5, 11, 58]. Likewise, IL-1 supports mechanisms as extracellular matrix binding modulation and regulation obtained in experimental studies that are essential for the response to mechanical stimuli in astrocytes [41]. In this sense, our results support the involvement of epigenetic regulation over cellular functional determinants in astrocytes during neurodegeneration but are necessary to develop more precise algorithms associated with gene screening [4].

Moreover, our model shows and support evidence from experimental studies, highlighting the expression and regulation of transporters such as the glutamate and lactate shuttle, redox stress reduction, transfer mitochondrial, among others, which are associated with the APAR mechanisms. Many of these biological functions associated with the response of astrocytes seem to be regulated by some of the tested histone modifications. Also, the response to external stimulus can be associated with the presence of neurotransmitter receptors, evidencing the neuron-astrocyte interaction beyond the metabolic support. Interestingly, we also report the presence of genes involved in the biosynthetic process of glutathione in the euchromatin regions, meaning a recurrent antioxidant activity process in astrocytes. Glutathione biosynthesis and release have been associated as a strategy for the balance and detoxify of the neural activity mediated by mitochondrial reactive oxygen species (ROS) in neurons linked to neurotransmission, neuroinflammation, neural disease etiology and progression [59, 60]. Glutathione biosynthesis is related to astrocytes antioxidant defense activity during pathological and non-pathological conditions.

Transcriptomic data, epigenetic landscape of TDAs, and histone modification regions data allowed the identification of APAR genes in the transcriptomic dataset and their localization (bivalent activation) [61]. TNFRSF1B, IL1R2, IL18RAP, IL1A, IL5RA, CXCL10, IL5, PIK3CG, IL10RA, and CCL8 genes were identified and associated with APAR mechanisms in astrocytes. Recently it has been demonstrated that during non-stimulating conditions, astrocytes secrete cytokines such as GM-CSF, CXCL1, CCL2, CXCL8, IL-6, and IL-8, all of those displayed at different levels [22, 37]. Moreover, administration of IL-1B and TNF activates astrocytes response with the production of cytokines IL-1B, IL-1RA, TNFA, CXCL10, CCL3, CCL5 and IL6 [62, 63, 64], being IL-6 response more efficient at higher concentration [65, 66].

Chromatin conformation in astrocytes has shown that PA response genes were located within shared TADs. During inflammation interleukin-1 receptor type II (IL1R2) has been described as a key receptor of which the expression reduces IL1A and IL1B activity [9]. On the other hand, the interleukin 18 receptor accessory protein (IL18RAP) that is associated with the pro-inflammatory response of IL18 by intracellular signaling was located in the same TAD region, suggesting that they share the same regulatory response when inflammatory processes occur in astrocyte [67]. Additionally, this TAD region also contains IL1R1 which is a key molecular mechanism associated with astrocytic response to inflammation by interaction with IL1A, IL1B and IL1R-agonists. Likewise, the TAD contains IL1RL2, and IL18R1, both interleukin receptors related to inflammatory cellular processes [5, 68, 69].

The coregulation of certain gene groups can also be associated with either master regulatory regions in TADs or architecture proximity regulation in the nucleus [70]. It is plausible that PA-lipotoxic responses in regulation of astrocytes by activating TAD regions depends upon extracellular signaling. This is possible because of the proximity of TAD to nucleus for cooperative organized regulation of genomic regions [44, 71]. Our results finally suggest that epigenetic modulation has an important role in the regulation of APAR mechanisms, yet further experiments are necessary to explore the TAD proximity involved in APAR regulation.

5. Conclusions

We present the first comprehensive data integration of epigenetic involvement in the astrocytic response to PA through the analysis from Hi-C, ChIP-Seq, and transcriptomic data in a multi-omic level. We described the role of epigenetics as a key mechanism of astrocytic PA response within which we found histones markers with bivalent capacity associated with repression of genomic activity (H3K4me3 and H3K27me3). This finding determines the adaptability and response to environmental stress, provided through complex astrocyte metabolic plasticity networks. In addition, our results showed that markers as H3K4me3, H3K27ac, H3K9ac, H3K4me1, H4K20me1, H3K36me3 and H3K79me2 have regions associated with homeostatic processes linked to exchange processes, regulation of the extracellular matrix, protein maintenance and ion channels regulation. These processes were found in euchromatin regions, highlighting that it is associated with essential basal functions in astrocytes. Likewise, signaling pathways modulation (i.e., PI3K/AKT), antioxidant activity (a recurrent mechanism in astrocytes), among others, were associated with glutathione biosynthesis processes, glutamate transport and glutamatergic neuronal support, identified as active basal coding regions.

APAR mechanisms proved to be highly regulated by histone modifications along the genome which is essential for the response to PA. Additionally, our results revealed the presence of highly regulatory regions in the TADs associated with IL1R2 and IL18RAP. Moreover, the location of genes encoding to interleukins in the genome and chromatin conformation revealed the putative epigenetic regulation of the inflammatory response. In this sense, our results support the involvement of APAR mechanisms on the lipotoxic effect of PA in astrocytes. While integrating transcriptomics with epigenetics data was possible to identify associated genes with APAR mechanisms and genes in response to PA located inside the topologically associated, genes found in the TAD region that shared the regulator responses linked to inflammatory processes were likely modulated by lipotoxicity actions. Additionally, it is possible to suggest that additional epigenetic mechanisms such as lncRNA, miRNA and extracellular signaling could be involved in the astrocytic response to PA. Considering that deterministic mechanisms of expression are still unknown for astrocytes in lipotoxic conditions, we suggest that epigenetic modulation is essential for an efficient and dynamic cellular response. This work is a novel approach that involves epigenetic regulation in the cellular response to PA-lipotoxicity in astrocytes. Therefore, it should be emphasized that it is recommended the development of new methodologies and algorithms for more accurate analysis associated especially to genetic encryption. Finally, an accurate investigation of this new multi-omic epigenetic-based model by integrating multiple underlying data sources about the cellular mechanisms of the response to PA-lipotoxicity in astrocytes, might help in the future to detect shared genetic patterns found in the TAD region among the neurodegenerative diseases, identifying biomarkers for differentiating disease states and thereby facilitating the decision-making process and treatment management.

Acknowledgments

We would like to acknowledge the members of the Computational and Experimental Biochemistry Group, whose help was indispensable for envisioning the research. RNA sequencing, and analysis were funded by Minciencias grants 8845, contract 886–2019, and 20304 awarded to the Pontificia Universidad Javeriana.

Conflict of interest

The authors declare no conflict of interest.

Gene (Regulation)	Location	TAD region	Genes in TAD
TNFRSF1B	chr1 12,166,949-12,209,228 [+Strand]	11,920,000-14,360,002	KIAA2013, PRAMEF13, KAZN, PLOD1, HNRNPCL2, AADACL3, PRAMEF5, MFN2, DHRS3, PRDM2, C1orf158, PRAMEF17, PRAMEF12, PRAMEF20, PRAMEF1, PRAMEF14, TNFRSF8, LRRC38, PRAMEF11, PDPN, HNRNPCL1, PRAMEF2, PRAMEF4, PRAMEF10, PRAMEF6, VPS13D, PRAMEF7, AADACL4, PRAMEF18, PRAMEF27, HNRNPCL3, PRAMEF25, PRAMEF26, HNRNPCL4, PRAMEF9, PRAMEF8, PRAMEF33, PRAMEF15
IL1R2	chr2 101,991,816-102,028,544 [+Strand]	101,880,002-102,560,000	SLC9A4, IL1R1, IL1RL1, IL1RL2, IL18RAP, IL18R1
IL18RAP	chr2 102,418,558-102,452,568 [+Strand]	101,880,002-102,560,000	SLC9A4, IL1R1, IL1RL1, IL1RL2, IL1R2, IL18R1
IL1A	chr2 112,773,915-112,785,394 [-Strand]	112,600,002-113,400,000	POLR1B, IL36G, PSD4, IL37, IL1F10, CHCHD5, IL36A, SLC20A1, IL36B, NT5DC4, IL36RN, CKAP2L, IL1B, IL1RN, PAX-AS1, PAX8
IL5RA	chr3 3,066,324-3,126,613 [-Strand]	2,360,002-3,160,000	CNTN4, TRNT1, CRBN
CXCL10	chr4 76,021,116-76,023,536 [-Strand]	75,760,002-76,440,000	USO1, NAAA, SCARB2, NUP54STBD1, PPEF2, CXCL11, SDAD1, FAM47E, CXCL9, FAM47E-STBD1, ART3, CCDC158
IL5	chr5 132,541,444-132,556,890 [-Strand]	132,440,002-133,360,000	IRF1-AS1, IL13, SOWAHA, CCNI2, KIF3AUQCRQ, FSTL4, SHROOM1, HSPA4, IRF1, GDF9, IL4, RAD50, SEPTIN8, LEAP2, AFF4, ZCCHC10
PIK3CG	chr7 106,865,278-106,908,980 [+Strand]	106,840,002-107,600,000	PRKAR2B, HBP1, COG5, GPR22, BCAP29, DUS4L
IL10RA	chr11 117,857,063-117,872,198 [+Strand]	117,240,002-118,320,000	DSCAML1, FXYD2, FYD6, CEP164, SMIM35, RNF214, PCSK7, PAFAH1B2, SIDT2, TAGLN, BACE1, TMPRSS13, TMPRSS4, SCN4B, SCN2B, JAML, MPZL3
CCL8	chr17 34,319,047-34,321,402 [+Strand]	33,520,002-35,720,000	ASIC2, CCL2, TMEM132E, FNDC8, CCL7, CCL11, CCT6B, PEX12, LIG3, CCL13, CCL1, ZNF830, NLE1, C17orf102, RFFL, AP2B1, RAD51D, UNC45B, SLC35G3, SLFN12L, SLFN5, SLFN13, SLFN11, SLFN14, SLFN12

Table A1.

Topological architecture of the differentially expressed genes associated with APAR mechanisms in astrocytes. All gene regions have been obtained from ENSEMBL (GRCh37/hg19). All the TAD described contained promoters, enhancers and promoter flanks.

Abbreviations

PA	palmitic acid
TADs	topological associated domains
BBB	blood–brain barrier
ROS	reactive oxygen species
APAR	astrocytic PA response