Gene regulatory effects of disease-associated variation in the NRF2 network.

Reactive oxygen species (ROS), which are both a natural byproduct of oxidative metabolism and an undesirable byproduct of many environmental stressors, can damage all classes of cellular macromolecules and promote diseases from cancer to neurodegeneration. The actions of ROS are mitigated by the transcription factor NRF2, which regulates expression of antioxidant genes via its interaction with cis-regulatory antioxidant response elements (AREs). However, despite the seemingly straightforward relationship between the opposing forces of ROS and NRF2, regulatory precision in the NRF2 network is essential. Genetic variants that alter NRF2 stability or alter ARE sequences have been linked to a range of diseases. NRF2 hyperactivating mutations are associated with tumorigenesis. On the subtler end of the spectrum, single nucleotide variants (SNVs) that alter individual ARE sequences have been linked to neurodegenerative disorders including progressive supranuclear palsy and Parkinson's disease, as well as other diseases. Although the human health implications of NRF2 dysregulation have been recognized for some time, a systems level view of this regulatory network is beginning to highlight key NRF2-targeted AREs consistently associated with disease.

1. NRF2-ARE-mediated gene regulation

Many cellular challenges – chemical, metabolic, and physical – disrupt redox homeostasis and generate excess reactive oxygen species (ROS). ROS have the potential to damage macromolecules including proteins, lipids, and DNA. The latter effect is especially harmful because it can cause DNA mutations with long-term consequences, but significant oxidative stress can also lead to apoptotic or necrotic cell death. Thus oxidative damage is a significant contributor to chronic diseases, from cancer to neurodegenerative disease [1-4]. In response to oxidative stress, cells activate a panel of cytoprotective genes, including antioxidant and detoxifying enzymes, that counteract ROS and ROS-induced damage to the cell. NRF2, a Cap-n-Collar (CNC) basic leucine zipper (bZIP) transcription factor encoded by the gene NFE2L2, is a master regulator of the transcriptional response to oxidative stress [5]. NRF2 is structurally and functionally conserved from insects to humans, and it must dimerize with one of the three small MAF (sMAF) bZIP proteins (MAFF, MAFG, MAFK) to bind ARE sequences and regulate gene expression [6]. NRF2 is widely expressed, but when ROS levels are low, nuclear NRF2 is kept low by direct interaction with the inhibitory protein KEAP1, which sequesters NRF2 in the cytoplasm and targets it for proteasomal degradation [7-10]. However, ROS modify KEAP1 and impair its ability to target NRF2 for degradation. So as ROS levels increase (oxidative stress), KEAP1’s ability to inhibit NRF2 decreases and nuclear NRF2 increases; the increase in nuclear NRF2 drives upregulation of many cytoprotective genes (Figure 1A) [10,11].

Figure 1

Regulatory effects of NRF2 network variation

(A) General schematic of NRF2 regulatory pathway. See text for details. (B) Mutations that disrupt NRF2-KEAP1 interactions lead to NRF2 hyperactivation and are associated with tumorigenesis. Select NRF2 target genes are consistently upregulated in tumors with hyperactivated NRF2. (C) Variants that disrupt individual ARE sequences create alleles with stronger ARE activity and weaker ARE activity. Such variants can alter NRF2 target gene expression and, in some cases, disease risk.

Although gene sets necessary for responding to and recovering from oxidative stress are key NRF2 targets, these are not the only genes regulated by NRF2 [12-14]. The NRF2 regulatory network also contains multiple autoregulatory loops, including a critical negative feedback loop: NRF2 transcriptionally activates its own repressor, KEAP1 [14,15]. The NRF2-KEAP1 negative feedback loop is deeply conserved, present in organisms ranging from Drosophila to human [14,16]. This ancient negative feedback loop highlights the importance of keeping theNRF2 pathway in check, and suggests precise regulation of NRF2’s nuclear concentration is paramount. The rules governing transcription factor interactions with DNA within the nucleus are complex, but are largely a function of transcription factor concentration and its binding affinity for a given target DNA sequence [17,18]. Presumably, the autoregulatory loops modulating NRF2’s nuclear concentration ensure that NRF2-ARE binding and cis-regulatory output at all NRF2 target genes is finely tuned for a wide range of stress conditions.

NRF2, in combination with one of the sMAF proteins, regulates gene expression by binding ARE sequences in its target genes’ enhancer or promoter regions. The ARE is also referred to as the electrophile response element (EpRE) or the CNC-sMAF-binding element (CsMBE) [19,20]. Although the term ARE is more restrictive than the latter two terms – NRF2-ARE transactivation is responsive to more than just oxidative stress – it is the most commonly used term for the NRF2-sMAF-binding element. The original ARE consensus sequence was defined as GCnnnSTCAY (where S = G or C, and Y = C or T) [21,22]. Consistent with the consensus motif, genome-wide NRF2 binding sites identified using ChIP-seq (chromatin immunoprecipitation followed by high throughput sequencing) are strongly enriched for the sequence TGCTGAGTCAY [12,13,23]. Thus in vivo NRF2 DNA binding is largely driven by a higher information content version of the original ARE consensus. The sequence identified by ChIP-seq is functionally relevant, because both NRF2-sMAF DNA binding (in vitro and in vivo) and NRF2-mediated regulatory output are correlated with a target ARE’s similarity to TGCTGAGTCAY [14,23,24].

Considering the finely tuned nature of NRF2-mediated gene expression, it is likely that genetic variants disrupting either controls on NRF2 nuclear concentration (Figure 1B) or NRF2-bound ARE sequences (Figure 1C) would have significant implications for NRF2 target gene expression. Indeed, variants of both types have now been characterized. These variants can have a marked effect on NRF2-mediated gene expression, and are associated with a range of pathologies.

2. Trans-regulatory effects of NRF2 or KEAP1 mutation

NRF2 activity is cytoprotective. Loss of NRF2 is associated with genomic instability and tumorigenesis, whereas activation of NRF2 is chemopreventive and promotes longevity [25-27]. Yet NRF2 activation beyond a certain threshold can be detrimental: somatic mutations that disrupt NRF2-KEAP1 interaction promote cancer progression [28-32]. Mutations altering either the NRF2 binding domain of KEAP1 or the KEAP1 binding domain of NRF2 were first observed in lung cancer [28,29,33]. Mechanistically, disruption of protein domains at the NRF2-KEAP1 interface prevents efficient targeting of NRF2 for proteasomal degradation, which in turn leads to NRF2 accumulation and constitutive activation of the pathway. Lung tissue is particularly prone to NRF2 hyperactivating mutations, but such mutations are also found in head and neck squamous cell carcinomas, endometrial cancer, and many other solid tumor types [32,34-37]. In addition to somatic mutation, DNA hypermethylation at the KEAP1 locus is also associated with KEAP1 repression, increased NRF2 activity, and tumorigenesis [38-42].

The exact mechanisms underlying NRF2’s oncogenic properties remain unclear, but likely involve aberrant induction of NRF2 target genes (Figure 1B). Importantly, data from The Cancer Genome Atlas (TCGA) indicate that mutations altering the KEAP1 binding domain of NRF2 are associated with similar patterns of gene upregulation across diverse tumor types (Table 1) [43,44]. Because these gene expression changes are a result of genetic variation at an unlinked locus (the NFE2L2 gene), they represent the trans-effects of mutations at NFE2L2. Many of the upregulated genes are near genomic regions bound by NRF2 based on human ChIP-seq data, so these presumably represent direct NRF2 targets (Table 1).

Table 1

Genes consistently upregulated in tumors with hyperactivated NRF2

Bold genes are direct NRF2 targets based on human ChIP-seq data [13,14]. Non-bold genes are direct NRF2 targets based on additional experimental evidence [48]. “Additional genes” are those with little evidence for direct regulation by NRF2. Tumor expression data are from [43].

Gene symbol	Gene name	Tumor type
Upregulated in 4/4 NRF2 hyperactivated tumor types:
ABCB6	ATP binding cassette subfamily B member 6	B, H, L, U
ALDH3A1	Aldehyde dehydrogenase 3 family member A1	B, H, L, U
FECH	Ferrochelatase	B, H, L, U
GCLM	Glutamate-cysteine ligase, modifier subunit	B, H, L, U
ME1	Malic enzyme 1	B, H, L, U
NQO1	NAD(P)H quinone dehydrogenase 1	B, H, L, U
PGD	Phosphogluconate dehydrogenase	B, H, L, U
SRXN1	Sulfiredoxin 1	B, H, L, U
TALDO1	Transaldolase 1	B, H, L, U
TKT	Transketolase	B, H, L, U
TXNRD1	Thioredoxin reductase 1	B, H, L, U
OSGIN1	Oxidative stress induced growth inhibitor 1	B, H, L, U
AKR1C1	Aldo-keto reductase family 1 member C1	B, H, L, U
EPHX1	Epoxide hydrolase 1	B, H, L, U
G6PD	Glucose-6-phosphate dehydrogenase	B, H, L, U
Additional genes: AKR1B10, AKR1C2, AKR1C3, CABYR, CES1, CYP4F3, JAKMIP3, PANX2, TRIM16L
Upregulated in 3/4 NRF2 hyperactivated tumor types:
GCLC	Glutamate-cysteine ligase catalytic subunit	B, H, U
GSR	Glutathione reductase	B, H, L
PRDX1	Peroxiredoxin 1	B, H, L
SLC7A11	Solute carrier family 7 member 11	B, H, L
TXN	Thioredoxin	B, H, L
PTGR1	Prostaglandin reductase 1	B, H, L
Additional genes: AGPAT9, CBR3, CES4, CLDN8, CYP4F11, FTHL3, KIAA0319, MAP2, MDGA1, NAMPT, RBM19, RIT1, SAMD5, SPP1, TDP2,    TSPAN7, WNT5A
Upregulated in 2/4 NRF2 hyperactivated tumor types:
ABHD4	Abhydrolase domain containing 4	H, L
ADAM17	ADAM metallopeptidase domain 17	B, L
COA6	Cytochrome c oxidase assembly factor 6	B, L
FTH1	Ferritin heavy chain 1	H, U
KEAP1	Kelch like ECH associated protein 1	B, L
MAF-G	MAF bZIP transcription factor G	B, L
PIR	Pirin	B, U
SLC12A8	Solute carrier family 12 member 8	H, L
SLC3A2	Solute carrier family 3 member 2	B, L
TLK1	Tousled like kinase 1	H, L
TMTC3	Transmembrane and tetratricopeptide repeat containing 3	B, L
ZNF746	Zinc finger protein 746	B, L
ABCC1	ATP binding cassette subfamily C member 1	B, H
ABCC3	ATP binding cassette subfamily C member 3	H, L
ADH7	Alcohol dehydrogenase 7	B, L
CBR1	Carbonyl reductase 1	B, H
GPX2	Glutathione peroxidase 2	H, L
IDH1	Isocitrate dehydrogenase 1	H, L
PRDX6	Peroxiredoxin 6	H, L
Additional genes: ABCA4, ADAM23, AKR1B15, ANXA10, ASF1A, ASPH, C14orf149, CREG1, DNAJB4, EPS8, EPT1, ETFB, FAM190A, FBXO30,    GLI2, GSTM3, GSTM4, HHIPL2, LOC729082, MAP1B, MEGF9, NECAB2, NSUN3, PHEX, PHKB, RAP1GAP, RNF217, SLC47A2, SLC9A3R1,    SMOC2, SOST, TNPO1, TPD52L1, TRIM16, TSKU, UGT1A7

Tumor type abbreviations: B = Bladder Urothelial Carcinoma; H = Head–Neck Squamous Cell Carcinoma; L = Lung Squamous Cell Carcinoma; U = Uterine Corpus Endometrial Carcinoma.

Oncogenic NRF2 provides a selective advantage to cells across diverse tissue environments, and the TCGA data demonstrate that it also casts a transcriptional shadow that is relatively consistent across these environments. Thus it is likely that subsets of the consistently upregulated genes are responsible for the metabolic, proliferative, and chemoresistance advantages afforded to cells with constitutive NRF2 activity. Multiple genes in Table 1 represent high priority candidates. For example, NRF2 targets necessary for generating the antioxidants glutathione and thioredoxin are consistently overexpressed in NRF2 gain-of-function tumors, and both antioxidants are essential for cancer initiation and progression [45]. Synthesis and regeneration of glutathione is largely controlled by the catalytic and modifier subunits of the glutamate-cysteine ligase enzyme (encoded by GCLC and GCLM, respectively), the cysteine transporter subunit encoded by SLC7A11, and glutathione reductase (encoded by GSR). Thioredoxin activity is induced by expression of TXN, which codes for this antioxidant, and regenerated by thioredoxin reductase (encoded by TXNRD1). But antioxidant genes are not the only group of NRF2 targets with potential cancer relevance. Upregulation of enzymes in the pentose phosphate pathway, including 6-phosphogluconate dehydrogenase (PGD), glucose-6-phosphate dehydrogenase (G6PD), transaldolase (TALDO1), and transketolase (TKT), drives metabolic reprogramming and NRF2-dependent proliferation in lung adenocarcinoma cells (A549) and other cell lines [46,47]. Additionally, a number of the NRF2 targets that are repeatedly induced in cancer likely promote chemoresistance by increasing drug metabolism (ALDH3A1, NQO1, AKR1C1, EPHX1) or transport (ABCB6) [48]. Combined, the genes that are consistently responsive to NRF2 activity in diverse tissue contexts can explain much of NRF2’s oncogenic potential.

The molecular profile of oncogenic NRF2 hyperactivation, at least with regard to gene expression signatures, is largely independent of a tumor’s tissue of origin. The upregulated genes in NRF2-driven cancer are responsible for many of NRF2’s core cytoprotective functions in non-pathological situations, so it makes sense that they can be activated in many cell types. Nevertheless, what mechanistically differentiates this significant subset of NRF2-targeted AREs from its other targets remains unclear. At the cancer-induced AREs, NRF2 binding and transactivation is ostensibly not inhibited by other tissue-specific regulatory networks (transcriptional repressors, chromatin environment, etc.), and this lack of constraint can have dire consequences.

3. Cis-regulatory effects of ARE variation

Genetic variants that disrupt individual NRF2-targeted ARE sequences are expected to have more specific effects than variants that alter overall levels of nuclear NRF2. Whereas disruption of overall NRF2 levels impacts many genes in the network (Table 1), a cis-regulatory ARE variant will primarily alter expression of the gene regulated by the ARE (Figure 1C). Although the effects of cis-regulatory variation are more precise, they are not without biological consequence. Most disease-associated variants identified by genome-wide association studies (GWAS) fall in non-protein coding DNA. Most of the non-coding disease-associated variants are innocuous – they simply reach significance because they are co-inherited with a functional variant – but those with functional relevance often disrupt transcription factor binding sites [49-51]. This is indeed the case for NRF2: polymorphic ARE motifs have recently been linked to allele-specific enhancer activity, gene expression, and disease risk [24,52].

A position weight matrix scanning approach identified over two million ARE sequences in the human genome, but NRF2 and sMAF ChIP-seq data suggest that less than 2% are functional AREs [52]. This is consistent with numbers observed for other human transcription factors [18,53]. Approximately 7.5% (2689 out of 35,659) of the AREs consistently bound by NRF2 and/or sMAFs contain a potential ARE-altering single nucleotide variant (SNV), and 14 of the variable AREs are in linkage disequilibrium (i.e., co-inherited) with disease-associated variants identified by GWAS [52]. Disease-associated, variant AREs represent instances where cis-regulatory variation in the NRF2 network might have a phenotypic impact (disease risk); a subset of these ARE-altering SNVs are outlined in Table 2. A connection to cancer is still evident: two polymorphic AREs are associated with testicular germ cell tumors. However variant AREs are also linked to disease beyond cancer. Interestingly, the list includes hits for neurodegenerative disorders including progressive supranuclear palsy (PSP), Parkinson disease (PD), and corticobasal degeneration, as well as gastrointestinal disorders including celiac disease and colitis. These disease associations are based on common germline SNVs (minor allele frequency greater than 1%), but rare inherited variants that disrupt AREs could also be important [24].

Table 2

Disease-associated, ARE-disrupting single nucleotide variants

Summary of significant SNVs identified in [52]. SNVs represented are those falling within 2 base pairs of an ARE containing the GCnnnnTCA core sequence, and with a position weight matrix (PWM) match score >10. For the ARE Sequences column, SNVs are underlined and highlighted in red/blue – the variant that generates a stronger PWM match is highlighted red, and the weaker PWM match is highlighted blue. Allele frequency data are from the 1000 Genomes Project.

SNP ID	ARE sequence	Allele frequency	Nearest gene	Disease association(s)
rs242561		0.90	MAPT	Progressive Supranuclear Palsy;Parkinson’s Disease; Corticobasal
		0.10		Degeneration; Interstitial Lung Disease
rs241032		0.57	CRHR1-IT1	Parkinson’s Disease
		0.43
rs6426833		0.59	RNF186	Ulcerative Colitis
		0.41
rs17035378		0.52	PLEK	Celiac Disease
		0.48
rs369184		0.85	TEX14	Testicular Germ Cell Tumor
		0.15
rs4818832		0.64	YBEY	Testicular Germ Cell Tumor
		0.36
Additional high priority polymorphic AREs: rs6426519, rs9603754, rs9884209, rs12638492, rs13067040, rs16857611, rs62033400, rs62094906

One ARE-altering SNV, rs242561, falls within an NRF2-bound ARE at the MAPT locus [52]. The major allele of rs242561 creates a mismatch ARE (CGCTGAGTCAC – variant sequence is underlined) and the minor allele creates a perfect ARE (TGCTGAGTCAC). Thus, most people carry one or two copies of a mismatch NRF2-targeted ARE in a MAPT enhancer region. A smaller subset of the population carries a perfect-match ARE at this MAPT enhancer. MAPT encodes the protein tau, which plays a central role in multiple neurodegenerative diseases, and rs242561 is in linkage disequilibrium with variants associated PSP, PD, and corticobasal degeneration [54-56]. Importantly, the minor allele of rs242561 (perfect-match ARE) is associated with decreased risk of all three aforementioned neurodegenerative disorders, acts as a hypermorphic ARE in reporter assays, and is associated with increased expression of a protective isoform of MAPT [52,57-59]. Although it is possible that additional non-coding variants affect MAPT expression and neurodegenerative disease risk, these data suggest that the ARE impacted by rs242561 plays a significant functional role at this locus.

The above example suggests that inherited variation in specific NRF2-targeted ARE sequences can influence gene expression and, ultimately, disease risk. However, somatic variation in ARE sequences might also have an impact on disease. Cancer genome sequencing data from TCGA indicate that somatic variants disrupting ARE-like sequences are under positive selection in cancer cells [60]. This intriguing finding places the ARE among a small subset of transcription factor binding motifs commonly mutated in cancer. Thus, ARE mutations and aberrant expression of select NRF2 target genes are also likely to play an important functional role in cancer.

When considering the effects of ARE variation, one must consider that our current models of NRF2 DNA binding and regulatory output are incomplete. Identification of ARE sequences are dependent on model used (position weight matrix or other models). Some models focus on the core ARE sequence described above, while others include significant stretches of flanking sequence – both may be functionally relevant. In addition, there are many variant AREs associated with gene expression changes that have not yet been linked to a disease phenotype (SL and MS, not shown). It is likely that many of these changes are not strong enough to have a phenotypic impact; however some may be revealed as important in future disease association studies.

It is also important to recognize that ARE activity is context-dependent. That is, ARE activity generally increases under various stress conditions. Therefore the cis-regulatory effects of some ARE variants might only be evident under conditions where AREs are active. Oxidative stress, from both exogenous and endogenous sources, is an early and ongoing contributor to many diseases, so it makes sense that some variant AREs will appear significant in standard disease association studies. This explanation holds for rs242561, as oxidative stress plays a significant role in the pathology of neurodegenerative diseases such as PD and PSP [61-64]. However, certain ARE variants may only display significant disease associations in the presence of other genetic variants or environmental stressors that disrupt redox homeostasis. GWAS are become increasingly expansive, allowing for exploration of gene-gene and gene-environment interactions, so it is possible that additional disease-associated ARE variants will be identified in future studies.

4. A role for additional ARE binding proteins?

Precise regulation of NRF2-ARE binding is clearly important. Nuclear NRF2 concentration and ARE quality (i.e., similarity to sequence TGCTGAGTCAY) are significant contributors to this protein-DNA interaction, however additional variables must be considered. NRF2 is part of a family of CNC transcription factors, all with very similar DNA binding properties. Mammalian genomes contain six CNC proteins: four transcriptional activators (NFE2, NRF1, NRF2, NRF3), and two repressors (BACH1, BACH2) [5]. With the exception of NFE2, which regulates developmental transitions in the hematopoietic system [65-67], all CNC factors have been implicated as regulators of stress responsive genes [5]. And like NRF2, the other CNC proteins all dimerize with sMAF proteins, and all bind ARE sequences to regulate gene expression [6].

Compared to NRF2, less is known about the other four stress responsive CNC factors. All are expressed in a variety of cell types except for BACH2, which is most prevalent in the brain and B cells [5]. Current models of CNC mediated gene expression posit that regulatory output at an ARE is the result of competition between activator and repressor CNC factors, with activator CNC factors dominating in stress conditions [68-72]. However, this competition model is based largely on the opposing actions of NRF2 and BACH1 at AREs associated with two canonical antioxidant genes (NQO1 and HMOX1), and may not apply equally to all CNC target genes [68,71,73]. Further complicating matters, non-CNC proteins can also modulate ARE activity. The ARE, like all cis-regulatory sequences, can be repressed by nucleosomes, which hinder transcription factor access to a binding motif. ARE-like motifs can also be directly repressed by MAF homodimers [74] (see for a [75] comprehensive review). Additionally, stress responsive AP-1 protein complexes, which consist of heterodimers of proteins from the FOS, JUN, and ATF families, bind a target sequence very similar to the ARE [76].

A complete understanding of NRF2-mediated gene regulation must take additional ARE-binding factors into account (Figure 2). Models that integrate cooperation and competition for ARE binding among the activating and repressive factors will certainly further our understanding of the normal, homeostatic functions of NRF2. More comprehensive models might also explain why only a subset of AREs is consistently misregulated in cancer, or why some ARE-disrupting SNVs affect disease risk while others do not. The genomics era has yielded tremendous insights into NRF2 biology. Approaches that view NRF2 in the context of additional ARE-binding factors, and integrate with disease related functional genomics data (e.g., TCGA and GWAS data), will provide a comprehensive view of the regulatory mechanisms at play in this network in both physiological and pathological contexts.

Figure 2

Potential ARE-binding transcription factor complexes

Current models suggest ARE regulatory output is driven by competition between activator and repressor CNC-MAF proteins. The impact of additional ARE binding proteins and nucleosomes, and whether this model extends equally to a wide range of functional AREs, remains unclear.