The Transcription Factor NRF2 Has Epigenetic Regulatory Functions Modulating HDACs, DNMTs, and miRNA Biogenesis

The epigenetic regulation of gene expression is a complex and tightly regulated process that defines cellular identity and is associated with health and disease processes. Oxidative stress is capable of inducing epigenetic modifications. The transcription factor NRF2 (nuclear factor erythroid-derived 2-like 2) is a master regulator of cellular homeostasis, regulating genes bearing antioxidant response elements (AREs) in their promoters. Here, we report the identification of ARE sequences in the promoter regions of genes encoding several epigenetic regulatory factors, such as histone deacetylases (HDACs), DNA methyltransferases (DNMTs), and proteins involved in microRNA biogenesis. In this research, we study this possibility by integrating bioinformatic, genetic, pharmacological, and molecular approaches. We found ARE sequences in the promoter regions of genes encoding several HDACs, DNMTs, and proteins involved in miRNA biogenesis. We confirmed that NRF2 regulates the production of these genes by studying NRF2-deficient cells and cells treated with dimethyl fumarate (DMF), an inducer of the NRF2 signaling pathway. In addition, we found that NRF2 could be involved in the target RNA-dependent microRNA degradation (TDMD) of miR-155-5p through its interaction with Nfe2l2 mRNA. Our data indicate that NRF2 has an epigenetic regulatory function, complementing its traditional function and expanding the regulatory dimensions that should be considered when developing NRF2-centered therapeutic strategies.


Introduction
The cells of multicellular organisms have the same genetic content, yet their functions are diverse because they express different genes. These differences are due in large part to epigenetic mechanisms that lead to heritable and stable changes in gene expression programs that occur through alterations in chromatin structure. Epigenetic gene regulation provides an adaptive layer of control of gene expression to enable the organism to adjust to a varying environment by eliciting histone modifications, DNA methylation, and gene silencing via microRNAs (miRNAs) [1,2]. Epigenetic mechanisms are not only key in the processes of cell development and differentiation, but they also play important roles in many pathologies, including neurodegenerative diseases.
Oxidative stress arises from an imbalance between reactive oxygen species and the cell's antioxidant capacity, leading to an accumulation of ROS and a disruption of the epigenetic state of the cell. In turn, oxidative damage triggers epigenetic changes in the chromatin

Bioinformatics Analysis
The script used in this study, with some modifications, was previously described [17]. Briefly, it uses, as input, a Browser Extensible Data (BED) file containing the chromatin immunoprecipitation sequencing (ChIP-seq) peaks for the transcription factors of analysis, a text file containing a list of RefSeq transcript accession numbers, and a position frequency matrix (PFM) file from the JASPAR database containing the consensus transcription factorbinding sites to be computed. Additionally, it makes use of the BED file at the UCSC Genome Browser, Table Browser resource (https://genome.ucsc.edu/cgi-bin/hgTables) (accessed on 3 February 2023) containing the locations of every transcript and its RefSeq accession number in the genome. To identify regulatory elements, a combined segmentation BED file was generated by concatenating Combined Segmentations [18] at the UCSC Genome Browser for the hg19 human genome, using BEDTools, or the ENCODE Candidate Cis-Regulatory Elements (cCREs) [19] combined from all cell type tracks was used for the mm10 mouse genome. The script retrieves the genomic coordinates for the desired transcripts, extends them 5000 bp upstream of the transcription start site, and intersects them with ChIP-seq peaks downloaded from all experiments in ChIP-Atlas [20] for the given transcription factors using the wrapper of BEDTools for Python, pybedtools [21,22]. In this analysis, all the available binding sites for NFE2L2, MAFF, MAFK, MAFG, and BACH1 or their mouse orthologs were downloaded and intersected with the extended transcripts of HDAC1, HDAC2, HDAC3, SIRT1, DNMT1, DNMT3A, DNMT3B, DROSHA, DGCR8, DICER1, and TARBP2 genes or their mouse orthologs. Then, the sequences of the ChIP-seq peaks were extracted using pybedtools from the FASTA file of the hg19 human genome (Supplementary Table S1) or the mm10 mouse genome (Supplementary Table S2). The profile for NFE2L2 was downloaded from the JASPAR database [23] in PFM format from the entry MA0150.1. Absolute frequencies were turned into a PSSM (position specificscoring matrix), containing scores through the log2(odds-ratio) (odds ratio: observed frequency/expected frequency). One unit was added as a pseudo-count to each absolute frequency to avoid log(0). The scoring of each site followed a similar procedure as we have previously described [5]. Briefly, a sliding window of a width dependent on the profile to be used was passed over the extracted sequences. Each nucleotide in the sliding window received a score according to the PSSM, and then, the score from each nucleotide was added up in order to provide an absolute score for the site. The relative score, the maximal and minimal scores, were obtained with a given PSSM and computed as (absolute score + |minScore|)/(|maxScore| + |minScore|). Sites with a relative score below 0.8 were discarded, and the remaining ones were provided as a BED file. In order to detect active regions, the script makes use of pybedtools to intersect the segmentation file with the regions described in the regulatory element bed file.

Analysis of mRNA Levels via Quantitative Real-Time PCR
Total RNA extraction, reverse transcription, and quantitative polymerase chain reaction (qPCR) were performed as detailed in previous articles [27]. Primer sequences are shown in Supplementary Table S3. Data were analyzed using the 2 −∆∆CT method, with normalization of the raw data based on the geometric mean of Actb and Gapdh (Applied Biosystems), encoding housekeeping proteins. All PCR amplifications were performed in triplicate.

Luciferase Assays
Transient transfections of MEFs or HT22 cells and luciferase assays were performed as described in [26]. pTK-Renilla was used as an internal control vector (Promega).

Antisense Oligonucleotide (ASO) Pull-Down Assay
To identify Nfe2l2 mRNA-associated miRNAs, ASO pull-down was performed using non-overlapping biotinylated ASOs recognizing LacZ (four ASOs) and Nfe2l2 (eleven ASOs). Incubation of whole-cell lysates with biotinylated ASOs was followed by incubation with Streptavidin-coupled Dynabeads™ (Invitrogen by Thermo Fisher Scientific, Waltham, MA, USA). RNAs were isolated from the pull-down materials, and qPCR was performed. Briefly, whole-cell lysates were prepared using RIPA buffer with a cocktail of protease and phosphatase inhibitors (Thermo Fischer Scientific, Waltham, MA, USA). One milligram of whole-cell extract was incubated with 1 µg of either Nfe2l2 ASO or control LacZ ASO for 16 h at 4 • C, whereupon RNA complexes were isolated with M-280 Streptavidin Dynabeads (Invitrogen) for 2 h at 25 • C. Total RNA was isolated with Trizol-Choloroform and cDNAs were synthesized with the Maxima Reverse Transcription kit following the manufacturer's protocols (Thermo Fisher Scientific, Waltham, MA, USA); real-time quantitative (q)PCR analysis was then performed using a PCR kit according to manufacturer's instructions (KAPA Biosystems, Wilmington, MA, USA), and the mRNA abundance was calculated via the 2 −∆∆CT analysis method, using Gapdh mRNA levels as the control transcript for normalization. miRs were reverse-transcribed using the MiR-X kit (Takara Bio, Shiga, Japan) and quantified via qPCR analysis using U6 as the normalization control RNA.

Statistical Analyses
Data are presented as the mean ± SEM. To determine the statistical test to be used, we employed GraphPad Instat 3, which includes the analysis of the data to a normal distribution via the Kolmogorov-Smirnov test. In addition, statistical assessments of differences between groups were analyzed (GraphPad Prism 8 by Dotmatics, San Diego, CA, USA) by performing an unpaired Student's t-tests. A one-way ANOVA with post-hoc Newman-Keuls test was used.

Results
The epigenome is comprised of modifications to chromatin, including histone modifications and DNA methylation. One of the main histone modifications is acetylation, which is often a necessary precursor to other modifications, such as phosphorylation, methylation, and ubiquitylation. Acetylation is controlled by the opposing functions of two families of enzymes: histone acetyltransferases (HATs) and histone deacetylases (HDACs). HDACs are involved in key biological functions, such as transcription, cell cycle, autophagy, DNA damage repair, stress responses, and senescence [28,29]. HDACs are classified according to their sequence similarities with yeast HDACs into class I, class II (IIa and IIb), class III, and class IV. HDACs have been found to regulate NRF2 signaling by directly modulating NRF2 acetylation [30]. In this study, we focused on the impact of NRF2 on the expression of HDAC1, HDAC2, and HDAC3. We also tested SIRT1 (a type-IV HDAC) given that its crosstalk with NRF2 has been previously described [31]. DNA methylation, another major epigenetic change, results from the transfer of a methyl group onto the cytosine to form 5-methylcytosine. DNA methylation regulates gene expression by recruiting proteins involved in gene repression or by inhibiting the binding of transcription factor(s) to DNA [32]. DNMTs include DNMT1, with a maintenance function to copy DNA methylation patterns from the parental DNA strand onto the newly synthesized daughter strand during DNA replication, and DNMT3a and DNMT3b, with de novo functions to establish new methylation patterns on unmodified DNA. DNA methylation is crucial for regulating tissue-specific gene expression, genomic imprinting, and X chromosome inactivation. Given that oxidative stress modulates DNA methylation levels in cancer [33], cardiovascular diseases, and type 2 diabetes [34], we sought to study if the expression of HDACs and DNMTs is dependent on the transcription factor NRF2.

Identification of Putative ARE Sequences in HDCAs and DNMTs
First, to define comprehensively the role of NRF2 in the transcriptional regulation of HDACs and DNMTs, we searched for putative ARE sequences in the ChIP-Atlas database, an integrative database covering almost all public data archived in the Sequence Read Archive of NCBI, EBI, and DDBJ, using ChIP-seq data [20] of the human (Table 1) or mouse ( Table 2) genomes for HDAC1, HDAC2, HDAC3, SIRT1, DNMT1, DNMT3a, and DNMT3b. Table 1. Table showing the putative ARE sequences identified in the human genome with a relative score over 80%. Columns 1 and 2 describe the genes containing the transcription factor-binding sites and their specific locations according to the GRCh37 (hg19) human genome assembly. The Motif column contains the specific sequence that was identified as a putative ARE sequence. The TGA motif in nucleotides 2, 3, and 4 of the motifs and the GC motif in nucleotides 9 and 10 of the motif are in bold, if present. Relative scores were calculated against the consensus binding sequence according to the position frequency matrix of the JASPAR database. The strand column indicates in which strand of the DNA, relative to the transcript sense is the putative ARE sequence: sense (+) or antisense (-). The Regulatory Element column provides the segmentation annotation from the Combined Segmentation Track at the UCSC Genome Browser. The TFs column indicates the transcription factor for which the ChIP-seq site belonged in the ChIP-Atlas database. * TSS, transcription start site; T, transcribed region; R; repressive or low activity region; PF, promoter flanking region; E; enhancer; WE, weak enhancer; CTCF, CTCF-enriched region. The ChIP-Atlas database includes experimental data from chromatin immunoprecipitation (ChIP) analysis of the ARE-binding transcription factors MAFG, MAFF, MAFK, BACH1, and NFE2L2. We used Python-based bioinformatic analysis to scan this binding region for the consensus ARE, as established in the JASPAR database [17,35]. Depending on the gene, we detected zero, one, or several putative ARE sequences with a relative score higher than 80%, a commonly used threshold for transcription factor binding-site analysis [36,37]. These putative ARE sequences in the promotor region have a high degree of similarity with the consensus ARE sequence (NTGACNNNGCN) described by [38]. As shown in Tables 1 and 3, bioinformatic analyses suggested that there are ARE sequences in HDACs, especially in DNMTs. We then performed functional assays to test the putative function of these sequences.  Table showing the putative ARE sequences identified in the mouse genome with a relative score over 80%. Columns 1 and 2 describe the genes containing the transcription factor-binding sites and their specific locations according to the GRCm38 (mm10) mouse genome assembly. The Motif column contains the specific sequence that was identified as a putative ARE sequence. The TGA motif in nucleotides 2, 3, and 4 of the motifs and the GC motif in nucleotides 9 and 10 of the motif are in bold, if present. Relative scores were calculated against the consensus binding sequence according to the position frequency matrix of the JASPAR database. The strand column indicates in which strand of the DNA, relative to the transcript sense is the putative ARE sequence: sense (+) or antisense (-). The Regulatory Element column provides the segmentation annotation from the ENCODE Candidate Cis-Regulatory Elements (cCREs) combined from all cell types of Track at the UCSC Genome Browser. The TFs column indicates the transcription factor for which the ChIP-seq site belonged in the ChIP-Atlas database. *, PLS, promoter-like signature; dELS; distal enhancer-like signature.

Gene
Coordinates

HDCAs Are NRF2-Dependent Genes
To test the bioinformatic predictions, we analyzed the levels of HDAC mRNAs in mouse embryonic fibroblasts (MEFs) derived from Nfe2l2 +/+ and Nfe2l2 −/− mice (Supplementary Figure S1A-verification of MEF genotype). As shown, the ablation of NRF2 led to a significant decrease in mRNA levels ( Figure 1A). These reductions were reflected at the protein level ( Figure 1B), further supporting the view that the production of HDAC1, HDAC2, HDAC3, and SIRT1 was regulated by NRF2. We then studied whether inducing NRF2 levels might enhance HDAC expression in hippocampal HT22 neurons (Supplementary Figure S1B,C). We treated cells with dimethyl fumarate (DMF), an activator of the NRF2 pathway through both the KEAP1 and GSK-3 pathways, as previously described [39]. Treatment with DMF induced a time-dependent increase in the levels of all HDAC mRNAs ( Figure 1C). These results were corroborated at the protein level, as we found that DMF treatment significantly induced protein levels of HDAC1, HDAC2, and SIRT1 ( Figure 1D). Taken together, these data indicate that NRF2 is able to modulate the expression of HDAC1, HDAC2, HDAC3, and SIRT1 in different cell types, suggesting a novel function with an epigenetic impact on this transcription factor. Newman-Keuls posterior test was used to evaluate differences in significance between groups: * p < 0.05, ** p < 0.01, *** p < 0.001 and **** p < 0.0001 compared to basal levels. The one-way ANOVA test with a Newman-Keuls posterior test was used to evaluate differences in significance between groups: * p < 0.05, ** p < 0.01, *** p < 0.001 and **** p < 0.0001 compared to basal levels.

NRF2 Is a Modulator of DNMTs Expression
After bioinformatic analysis of ARE elements in the promoters of DNMTs, we followed the same strategy for the analysis of HDACs. First, we analyzed the levels of Dnmt1, Dnmt3a, and Dnmt3b mRNAs in NRF2-deficient MEFs. We observed that the levels of Dnmt1 and Dnmt3b mRNAs were significantly reduced, while the changes in Dnmt3a mRNA were more modest (p = 0.07) (Figure 2A). These changes were mirrored at the protein level for DNMT1; DNMT3a levels were moderately reduced, while DNMT3b was not detectable (see Figure 2B). We then analyzed the effects of inducing NRF2 by treating HT22 cells with DMF. As shown, DMF treatment induced DNMT1 and DNMT3b at both mRNA and protein levels, validating the effects observed in MEFs and further supporting the notion that they are regulated by NRF2. For DNTM3a, there was a slight induction of mRNA levels but not at the protein level ( Figure 2C,D); further analysis is needed to determine whether the expression of DNMT3a is regulated by NRF2.

NRF2 Is a Modulator of DNMTs Expression
After bioinformatic analysis of ARE elements in the promoters of DNMTs, we followed the same strategy for the analysis of HDACs. First, we analyzed the levels of Dnmt1, Dnmt3a, and Dnmt3b mRNAs in NRF2-deficient MEFs. We observed that the levels of Dnmt1 and Dnmt3b mRNAs were significantly reduced, while the changes in Dnmt3a mRNA were more modest (p = 0.07) (Figure 2A). These changes were mirrored at the protein level for DNMT1; DNMT3a levels were moderately reduced, while DNMT3b was not detectable (see Figure 2B). We then analyzed the effects of inducing NRF2 by treating HT22 cells with DMF. As shown, DMF treatment induced DNMT1 and DNMT3b at both mRNA and protein levels, validating the effects observed in MEFs and further supporting the notion that they are regulated by NRF2. For DNTM3a, there was a slight induction of mRNA levels but not at the protein level ( Figure 2C,2D); further analysis is needed to determine whether the expression of DNMT3a is regulated by NRF2.

NRF2 Regulates the Expression of Proteins Implicated in miRNA Biogenesis
Besides histone acetylation and DNA methylation, miRNAs (small non-coding RNAs) reduce the stability and translation of target mRNAs and are also epigenetic regulators. ROS can modulate miRNA biogenesis at many levels, and several enzymes and components of the miRNA processing machinery can be affected by ROS [40]. Therefore, we examined the possibility that expression of the proteins involved in microRNA biogenesis (DGCR8, DROSHA, DICER, and TARBP2) might be regulated by NRF2.
Therefore, as before, we investigated bioinformatically whether there were ARE sequences in the promoter regions of the DGCR8, DROSHA, DICER1, and TARBP2 genes. We observed that, especially in the human genome (Table 3), there were several possible ARE sequences within all genes analyzed. In the mouse genome, we only found possible ARE sequences in the regulatory regions of the Dicer1 and Drosha genes ( Table 4). In MEFs from Nfe2l2 +/+ and Nfe2l2 -/mice, we observed significant decreases in Dgcr8, Drosha, Dicer1, and Tarbp2 mRNAs ( Figure 3A). Similarly, the expression levels of DGCR8, DROSHA, and DICER1 proteins decreased greatly in the absence of NRF2 ( Figure 3B), while the levels of TARDBP2 protein did not change. On the other hand, treatment of HT22 cells with DMF significantly induced the levels of Dgcr8, Drosha, Dicer1, and Tarbp2 mRNAs ( Figure 3C), as well as the proteins DGCR8 and DICER, while TARBP2 was less induced and DROSHA was unchanged ( Figure 3D). These data suggest that many proteins involved in miRNA biogenesis are regulated by NRF2.  cells described in A), the levels of DGCR8, DROSHA, DICER, and TARBP2 proteins were analyzed via immunoblotting and quantification based on densitometry using the same samples as in A). n = 3−4 samples per experimental group ± SEM. Asterisks denote significant differences: ** p < 0.01 and **** p < 0.0001, comparing the indicated group with the wild-type mice according to a Student's t-test. Time-course analysis of treatment of hippocampal HT22 cells with DMF (20 μM). C) qPCR analysis of the levels of Dgcr8, Drosha, Dicer1, and Tarbp2 mRNAs. n= 4−5 samples ± SEM. D) In the cells described in C), the levels of DGCR8, DROSHA, DICER, and TARBP2 were analyzed via immunoblotting and respective quantification based on densitometry. n = 3−4 samples per experimental group ± SEM. The one-way ANOVA test with a Newman-Keuls posterior test was used to evaluate differences in significance between groups: *p < 0.05 and **p < 0.01 compared to basal levels. (D) In the cells described in (C), the levels of DGCR8, DROSHA, DICER, and TARBP2 were analyzed via immunoblotting and respective quantification based on densitometry. n = 3−4 samples per experimental group ± SEM. The one-way ANOVA test with a Newman-Keuls posterior test was used to evaluate differences in significance between groups: * p < 0.05 and ** p < 0.01 compared to basal levels.

NRF2 Modulates miRNA Expression
The sequences of 3 -untranslated regions (3'UTRs) of messenger RNAs (mRNAs) govern their stability, localization, and expression [41]. miRNAs typically bind to the 3 UTRs of target mRNAs with which they share partial complementarity and reduce their stability and translation. Therefore, we evaluated whether miRNAs interacting with the 3'UTR of NFE2L2 mRNA suppressed NRF2 expression and whether NRF2 was able to modulate this loop. To test this hypothesis, the ability of miRNAs to regulate the 3'UTR of Nfe2l2 mRNA was evaluated using luciferase reporters. MEFs from Nfe2l2 +/+ and Nfe2l2 −/− mice were transfected with the pSGG luciferase vectors bearing the NRF2 3'UTR, along with a control renilla vector. As shown in Figure 4A, the absence of NRF2 led to increased expression of miRNAs that negatively modulate the expression of the 3'UTR. Similarly, these data were confirmed in HT22 cells, where transfection of a dominant negative NRF2 (DN-NRF2) negatively regulated the expression of the 3'UTR of NRF2 ( Figure 4C). In contrast, an increase in NRF2 levels, either in KEAP1-deficient MEFs (where there is an increase in NRF2 levels) ( Figure 4B) or via transfection of the NRF2-∆ETGE-V5 ( Figure 4C) that lacks four residues (ETGE) essential for recognition by the E3 ligase complex Cul3/KEAP1, led to a decrease in the levels of miRNAs that inhibit the 3'UTR of NRF2 and thus an increase in its expression. These data were verified using the NRF2 inducer DMF at different doses ( Figure 4D). The results indicate a dose-response effect in which the higher the increase in NRF2 levels the greater the de-repression of the reporter 3'UTR. Although in these experiments we cannot rule out the impact of other factors, such as RNA-binding proteins, other non-coding RNAs, or RNA modifications on this 3'UTR, in addition to the action of miRNAs, the data are consistent with a model whereby NRF2 reduces the levels of specific microRNAs capable of binding the Nfe2l2 mRNA and reducing Nfe2l2 mRNA and NRF2 protein.  (DN), or pcDNA3.1/V5HisB-mNRF2∆ETGE, respectively. Asterisks denote significant differences: ** p < 0.01, *** p < 0.001 and **** p < 0.0001, comparing the indicated group with the control group according to a Student's t-test. (D) HT22 cells were transfected with pSGG-NRF2-3'UTR and treated with DMF at 6 µM, 20 µM, or 60 µM for 16h. All luciferase activities were measured 24 h after transfection. The luciferase activities were normalized to the renilla luciferase activities from the co-transfected reporter. The relative luciferase activities were calculated by normalizing them to those of controls. n = 3 samples per experimental group ± SEM. The one-way ANOVA test with a Newman-Keuls posterior test was used to evaluate differences in significance between groups: ** p < 0.01 and *** p < 0.001 compared to basal levels.
As mentioned above, miRNAs control target gene expression by inhibiting translation and degrading target RNAs. In addition, there is evidence that mRNAs can affect microRNA activity in two manners [42]. First, some mRNAs can function as competing endogenous RNAs (ceRNAs), as seen when two endogenous targets compete with each other for binding to a shared miRNA [43]. In this case, if one of the endogenous mRNA targets changes its expression, the activity of miRNAs upon the other targets will change accordingly. This mechanism appears to have been ruled out, since NRF2 overexpression itself reduced the levels of the miRNAs that bind to it. The second mechanism is targetdirected miRNA degradation (TDMD) [44]. In TDMD, the target mRNA promotes the degradation of the miRNA that binds to it [45]. To test whether this might be what was happening with NRF2, we treated HT22 cells with DMF (conditions similar to Figure 4D, 20 µM) and analyzed the expression levels of miR-27a-3p, miR-128-3p, and miR-155-5p. As observed in Figure 5C, NRF2 induction promoted a significant decrease in miR-155-5p levels, while other miRNAs were unchanged. Finally, we sought to determine in which biological processes miR-155-5p is involved, mainly at the neuronal level, in order to determine the role that its modulation through NRF2 may have. To do this, we determined its targets through the TargetScan database, and using the ShinyGo 0.76.2 platform and its pathway database, we performed a GO biological process analysis. Although miR-155-5p participates in many processes, it was found to be involved in the regulation of 45 mRNAs encoding proteins that participate in the development of the central nervous system (CNS) and~64 mRNAs encoding proteins implicated in neurogenesis (Figure 6), respectively. Thus, a more exhaustive analysis of the involvement of miR-155-5p and the possible regulation by NRF2 in the development of the central nervous system is required.   Using the TargetScan platform, the list of genes that can be modulated by miR-155-5p was extracted. These genes were analyzed using the ShinyGO 0.76.2 pathway database GO biological process platform.

Discussion
Although the expression of the transcription factor NRF2 at the epigenetic level has been studied extensively, the reciprocal process, i.e., the impact of NRF2 on epigenetic gene regulation, has not been explored. In this study, we show for the first time that NRF2 is able to regulate the expression of type-I HDACs (HDAC1, HDAC2, HDAC3) and SIRT1, DNMTs, and proteins involved in miRNA biogenesis, DROSHA, DGCR8, DICER1, and TARBP2. Our data further suggest that NRF2 may be involved in the regulation of the levels of certain miRNAs through TDMD. These new data are of great potential relevance to the pharmacological regulation of NRF2 as a therapeutic strategy for various pathologies since epigenetic changes triggered by NRF2 will also have to be weighed.
NRF2 is a pleiotropic transcription factor capable of regulating the expression of genes involved in different processes, including xenobiotic, redox, and carbohydrate metabolism, inflammation, and proteostasis [5] (Figure 7). Therefore, its dysregulation has been described in a multitude of pathologies, in many cases along with epigenetic changes that cause aberrant gene expression programs and the loss of homeostasis. These data suggest that NRF2 affects many cellular functions underlying disease, although to-date, we only understand how epigenetic modifications affect the expression and function of the NRF2 pathway. The fact that NRF2 can promote the expression of type-I HDAcs (Figure 1), DNMTs (Figure 2), and proteins involved in miRNA biogenesis (Figure 3) opens Neuronal-related functions are highlighted (blue). Using the TargetScan platform, the list of genes that can be modulated by miR-155-5p was extracted. These genes were analyzed using the ShinyGO 0.76.2 pathway database GO biological process platform.

Discussion
Although the expression of the transcription factor NRF2 at the epigenetic level has been studied extensively, the reciprocal process, i.e., the impact of NRF2 on epigenetic gene regulation, has not been explored. In this study, we show for the first time that NRF2 is able to regulate the expression of type-I HDACs (HDAC1, HDAC2, HDAC3) and SIRT1, DNMTs, and proteins involved in miRNA biogenesis, DROSHA, DGCR8, DICER1, and TARBP2. Our data further suggest that NRF2 may be involved in the regulation of the levels of certain miRNAs through TDMD. These new data are of great potential relevance to the pharmacological regulation of NRF2 as a therapeutic strategy for various pathologies since epigenetic changes triggered by NRF2 will also have to be weighed.
NRF2 is a pleiotropic transcription factor capable of regulating the expression of genes involved in different processes, including xenobiotic, redox, and carbohydrate metabolism, inflammation, and proteostasis [5] (Figure 7). Therefore, its dysregulation has been described in a multitude of pathologies, in many cases along with epigenetic changes that cause aberrant gene expression programs and the loss of homeostasis. These data suggest that NRF2 affects many cellular functions underlying disease, although to-date, we only understand how epigenetic modifications affect the expression and function of the NRF2 pathway. The fact that NRF2 can promote the expression of type-I HDAcs (Figure 1), DNMTs (Figure 2), and proteins involved in miRNA biogenesis (Figure 3) opens new perspectives on the spectrum of actions of NRF2, its epigenetic influence, and its implications in disease. In this study, we have focused on the study of type-I and SIRT1 HDACs, but in further studies, it will be interesting to determine the involvement of NRF2 related to the other HDACs and to determine whether the results of the study presented here can be extrapolated to all types of HDACs or is specific to type-I HDACs and SIRT1.
new perspectives on the spectrum of actions of NRF2, its epigenetic influence, and its implications in disease. In this study, we have focused on the study of type-I and SIRT1 HDACs, but in further studies, it will be interesting to determine the involvement of NRF2 related to the other HDACs and to determine whether the results of the study presented here can be extrapolated to all types of HDACs or is specific to type-I HDACs and SIRT1.  Our studies suggest that when NRF2 levels are decreased in pathologies, this reduction may lead to decreased expression of Type-I HDACs and increased expression of genes that were repressed under physiologic conditions. Similarly, NRF2 deficiency would also lower DNMT levels, in turn inducing the expression of genes that were previously repressed. On the other hand, NRF2 inducers, such as sulforaphane, were found to induce NRF2 regulation at the epigenetic level, mainly associated with DNA methylation [47]. Epigenetic changes in the mechanism of action of DMF have also been described [48][49][50]. In contrast with the several agents that function as NRF2 inducers, very few molecular components have been recognized as NRF2 inhibitors. Brusatol, luteolin, trigonelline, and retinoic acid are several compounds that have been described as having inhibitory effects on NRF2 signaling [51], but their effects at the epigenetic modulation level have not yet been described. In future experiments, it will be interesting to analyze the effect of NRF2 inhibitors on the levels of HDACs, DNMTs, and miRNAs, to determine the possible impact of NRF2 inhibitor treatments at the epigenetic level. Our data support the fact that the pharmacological regulation of NRF2 involves different downstream effectors (Figure 7), including HDACs, DNMTs, and miRNAs, and thus broadening the spectrum of action of this transcription factor.
The results presented in this study can be potentially relevant to a wide range of pathologies where epigenetic mechanisms are of particular importance, such as cancer [52,53], allergies [54], and neurodegenerative diseases [55,56]. In cancer, specific methylation and other alterations of the NFE2L2 promoter have been documented [12,57,58], which can alter the expression levels of NRF2, linked to carcinogenesis through metabolic reprogramming, tumor promotion, inflammation, and resistance to therapy. NRF2 silencing or pharmacological inhibition by brusatol reduced the proliferation and migration of breast cancer (BC) cells, inhibited proliferation, activated apoptosis, sensitized BC cells to cisplatin in vitro, and slowed tumor cell growth in vivo [59]. Along these lines, esophageal adenocarcinoma (EAC) displayed increased NRF2, and both the knockdown of NRF2 and pharmacological inhibition with brusatol inhibited tumor growth by inducing ferroptosis and apoptosis [60]. In the case of metastatic Ewing sarcoma and osteosarcoma, it has been described that oxidative stress attenuates metastasis; here, treatment with the class-I HDAC inhibitor MS-275 inhibited the deacetylation of YB-1 (Y-box binding protein 1), reduced its binding to the 5'UTR of NFE2L2, reduced the translation of NRF2, and increased the levels of intracellular ROS [61]. However, other studies showed that HDAC inhibitors increased NRF2-signaling in tumor cells [62]. Therefore, further studies will be needed to look at the crosstalk between HDACs and NRF2, taking into consideration that NRF2 has pleiotropic roles in cancer cells [51].
Outside of cancer, type I and II HDAC inhibition mediated by Trichostatin A (TSA) activated transcription factor NRF2 and protected against cerebral ischemic damage. On the other hand, SIRT1 activation was also found to induce the NRF2 signaling pathway, with beneficial effects in focal cerebral ischemia [63], indicating that further studies are needed to unravel the crosstalk between HDAC and NRF2 in this pathology [64]. In neurodegenerative diseases, HDAC inhibitors were found to improve the redox balance and attenuate neuronal degeneration in Huntington's disease [65] and amyotrophic lateral sclerosis [66] mouse models and in Alzheimer's disease-like pathological changes in SH-SY5Y neuroblastoma cells [67]. One mechanism of action described in this regard is that HDACs remove acetyl groups in histones associated with the KEAP1 promoter region, inducing an increase in KEAP1 transcription, and therefore, the inhibition of HDAC might have opposite effects [65]. Additionally, NRF2 levels were found to be elevated in many neurodegenerative diseases [27,68,69], and thus, further studies are warranted to fully understand the interaction between NRF2 and HDACs and the therapeutic value of interventions directed at NRF2 and/or HDACs. Similar to HDACs, DNMTs are also potential therapeutic targets, since alterations in their activity have also been associated with various pathologies. Therefore, as with HDACs, more studies are needed to establish the value of therapeutic strategies that modulate NRF2 and DNMTs pharmacologically [70][71][72].
A multitude of miRNAs are predicted to repress NRF2 production [73], and many miRNAs are involved in oxidative stress processes associated with physiological and pathological conditions [74]. The fact that oxidative stress is capable of regulating the expression of proteins involved in the biogenesis of miRNAs [75][76][77] led us to evaluate the potential involvement of NRF2 in this mechanism. Our data indicate that NRF2 promotes the production of proteins DROSHA, DGCR8, DICER1, and TARBP2 involved in miRNA biosynthesis (Figure 3), expanding the function of this transcription factor into the post-transcriptional space. Although miRNAs typically modulate the stability and translation of their target mRNAs, we have shown that the absence of NRF2 leads to the increased expression of miRNAs that negatively modulate the expression of the 3'UTR-Nfe2l2 ( Figure 4A,C) and vice versa, underscoring the fact that NRF2 levels can also repress the actions of some miRNAs. More detailed studies are necessary to elucidate the specific miRNAs involved in this mechanism. As mentioned above, we cannot rule out the impact of other factors, such as RNA-binding proteins, other non-coding RNAs, or RNA modifications on this 3'UTR, in addition to the action of miRNAs. Further studies of additional regulatory factors should also be considered.
Oxidative stress is one of the main inducers of the NRF2 pathway, and thus, its activation is linked to the induction of oxidative stress-associated miRNAs, the so-called "redoximiRs", such as miR-27a-3p and miR-155-5p [75,78]. In HT22 hippocampal cells, our data indicate that of all miRNAs analyzed ( Figure 5A), only miR-27a-3p, miR-27b-3p, miR-128-3p, and miR-155-5p associate with Nfe2l2 mRNA ( Figure 5B). There is previous evidence that miR-27a-3p and miR-27b-3p are redox-sensitive miRNAs [79][80][81][82] and modulate NRF2 levels, in agreement with our results. For example, in maternal diabetes-induced oxidative stress, miR-27a-3p levels rise, in turn suppressing NRF2 production [81]. Furthermore, an analysis of miRNA signatures in transgenic mice expressing a constitutively active, cardiacspecific NRF2 (caNrf2-Tg) [15] revealed that increasing the levels of NRF2 leads to reduced miR-155-5p levels. These results support our data that the induction of NRF2 decreased miR-155-5p levels ( Figure 5B). Here, we focused on "redoximiRs" and found that NFE2L2 mRNA might drive the degradation of specific miRNAs mediated by TDMD and thereby reduce the levels of specific miRNAs. Further experiments are necessary to determine the exact mechanism by which NRF2 causes the degradation of other microRNAs, such as miR-155-5p, given its implication in neuroinflammation and other pathologies, and is the main miRNA induced by LPS treatment in microglia [83,84]. Furthermore, miR-155 alters the expression of genes that regulate axon growth [85], supporting the bioinformatic prediction that miR-155 can regulate the expression of genes involved in CNS development and neurogenesis ( Figure 6). Therefore, the modulation of this miR-155-5p could be of great relevance in relation to neurodegenerative diseases, such as Alzheimer's disease. The precise regulation of miR-155 and other microRNAs by NRF2 awaits future study. It would also be interesting to determine, in future experiments, whether these miRNAs are specific to a specific neuronal type (e.g., hippocampal neurons) or have a more general functions, as well as whether there are differences between the various brain cell types.

Conclusions
This study demonstrates for the first time that the molecules involved in epigenetic gene regulation, namely HDACs, DNMTs, and miRNA biogenesis factors, are under the control by NRF2, expanding the field of action of this transcription factor.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/antiox12030641/s1, Figure S1: Verification of MEF genotype and DMF treatment, Suppl. Table S1-List of transcripts scanned for putative ARE sequences in the human genome (hg19); Suppl. Table S2-List of transcripts scanned for putative ARE sequences in the mouse genome (mm10); Suppl. Table S3. List of primers used in this study; Suppl. Table S4. List of antibodies used in this study