A functionally defined high-density NRF2 interactome reveals new conditional regulators of ARE transactivation

NRF2 (NFE2L2) is a cytoprotective transcription factor associated with >60 human diseases, adverse drug reactions and therapeutic resistance. To provide insight into the complex regulation of NRF2 responses, 1962 predicted NRF2-partner interactions were systematically tested to generate an experimentally defined high-density human NRF2 interactome. Verification and conditional stratification of 46 new NRF2 partners was achieved by co-immunoprecipitation and the novel integration of quantitative data from dual luminescence-based co-immunoprecipitation (DULIP) assays and live-cell fluorescence cross-correlation spectroscopy (FCCS). The functional impact of new partners was then assessed in genetically edited loss-of-function (NRF2−/−) and disease-related gain-of-function (NRF2T80K and KEAP1−/−) cell-lines. Of the new partners investigated >77% (17/22) modified NRF2 responses, including partners that only exhibited effects under disease-related conditions. This experimentally defined binary NRF2 interactome provides a new vision of the complex molecular networks that govern the modulation and consequence of NRF2 activity in health and disease.


Introduction
The human NFE2L2 (NRF2) gene encodes a cap'n'collar (CNC) basic leucine zipper (bZIP) transcription factor. In addition to its well established role in redox homeostasis [1,2] there is increasing evidence that NRF2 has pleotropic effects on many cellular processes, driving expression of >240 functionally diverse genes [3]. The fundamental importance of the human NRF2 response network is emphasised by its association with >60 human diseases [4]. However, despite intense research, details of conditional regulation, functional-crosstalk or feedback within the human NRF2 network remain incomplete.
Under basal conditions NRF2 is sequestered in the cytoplasm through a direct interaction with KEAP1. This interaction leads to NRF2 ubiquitination and proteasomal degradation, ensuring low basal activity. Upon oxidative or xenobiotic insult, the KEAP1/NRF2 interaction is perturbed, preventing NRF2 degradation and allowing nascent NRF2 to drive transactivation of ARE-responsive genes, including NQO1, GCLM, and HMOX1 [15][16][17][18][19]. However, there is increasing evidence that NRF2 activity can be modulated by other processes, including pro-inflammatory mitogen-activated kinase (MAPK) [20], nuclear factor kappa-B (NF-κB) [21], ER stress [5], and p53 response pathways [22]. In this context, it is conceivable that in order to maintain homeostasis, NRF2 may have evolved the ability to mediate the dynamic integration of responses from multiple stress response pathways. The diverse spectrum of NRF2 response genes and predicted interaction partners also implies a potential for complex mechanisms of conditional regulation and functional crosstalk with other biological processes.
Although induction of NRF2 counteracts the detrimental effects of inflammation [23], cardiovascular disease [24] and DNA damage [25], genetic mutations that constitutively enhance NRF2 activity are linked to therapeutic resistance and poor prognosis in several forms of cancer [26][27][28][29]. As such, there is a need to develop new strategies to selectively modulate NRF2 responses under different pathological conditions. Unfortunately, lack of experimental evidence and systems-level analyses impedes mechanistic understanding of how NRF2 responses may be modulated, or the effects that changes in NRF2 activity may have on other cellular processes. Therefore, a more detailed functionally defined map of the human NRF2 interactome is required to inform development of improved therapeutic strategies. In addition to the 37 known binary human NRF2 interaction partners, a large-scale in silico study provided an intriguing insight into the potential complexity and diversity of the human NRF2 interactome, and its connectivity to other physiological and disease-related processes [30]. As a large proportion of these putative interactions were untested, we sought to generate a new high-density experimentally defined human NRF2 interactome, including verification, stratification, and functional analysis of novel Outer shell represents edges re-confirmed in co-immunoprecipitation assays. Transcription factors (TFs) are bordered in red, bZIP TFs in pink, and a kinase in light green. As the Neh4 and Neh5 domains tended to auto-activate Y2H reporters we were unable to define partner profiles for most of these regions. Auto-activating clones are shown in red. One Neh4+5+7 prey clone did not auto-activate and was used in the screen. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.) partners. To maximise network coverage partner interaction screens were initially performed in high-stringency targeted yeast two-hybrid assays, which have the ability to detect weak yet functionally verifiable interactions [31,32]. In addition to conventional co-immunoprecipitation studies, a novel approach of combining live-cell fluorescence cross-correlation spectroscopy (FCCS) with high-sensitivity dual luminescence-based co-immunoprecipitation (DULIP) assays was used to calculate in vivo partner K d values and stratify conditional partner preferences. Finally, a predictive analysis of biological function and disease association was performed to provide new insight into potential mechanistic links to human disease phenotypes.

Experimental screening of predicted human NRF2 binary interaction partners
A speculative binary human NRF2 network was generated by prediction of domain-motif and domain-domain interactions using the Eukaryotic Linear Motif (ELM) and PRINCESS PPI-evaluation tools, as previously described [30]; with the inclusion of: (i) additional mammalian interologs from BioGRID [33,34], (ii) interologs involving the C. elegans NRF2 ortholog SKN-1 in Wormbase [35], WormNet [36], and the CCSB Worm ORFeome [37], (iii) interologs from the Drosophila melanogaster ortholog CncC in Flybase [38] and BioGRID databases, as well as (iv) predictions from the Reactome Functional Interactome (ReactomeFI) [39]. This resulted in a predicted network containing 187 potential binary NRF2 interaction partners (Supplementary Table S1); including 132 previously predicted NRF2 partners [30], together with additional non-mammalian interologs. In total, 133 of these putative binary NRF2 partners were cloned into yeast two-hybrid (Y2H) prey (pACTBE-B/pACTBD-B) and bait (pGBAE-B/pGBAD-B) expression vectors [40]. Following removal of auto-activating clones, 87 unique bait and 124 unique prey clones were generated (Supplementary Table S2). To identify NRF2 domain-specific interactions, fragments corresponding to individual and combined NRF2 Neh domains were generated, alongside full-length NRF2 (Fig. 1A). As constructs containing NRF2 Neh4 and Neh5 had a propensity to auto-activate it was not possible to use most of these clones in Y2H screens. However, one Y2H prey clone (Neh4+5+7) did not autoactivate and was included in the screen. Using all non-auto-activating clones, a total of 1962 binary interactions were systematically tested in a series of repeated Y2H matrix screens, resulting in the identification of 46 reproducible positive NRF2 partners, which show strong evidence of partner/domain specificity (Fig. 1B). Significantly, 28% (13/46) of the binary NRF2 partners identified in this screen were predicted interologs: seven of which were conserved in Drosophila, four in C. elegans, and two in mice, thereby justifying the inclusion of more distant orthologs in this study.
Initially, 43 positive NRF2 partners were expressed as mCherrytagged constructs and tested by co-immunoprecipitation with EGFPtagged NRF2. In total, 62.8% (27/43) of partners were reproducibly co-immunoprecipitated with EGFP-tagged NRF2 (Fig. 1C, Supplementary Fig. S1). As a high number of positive interaction partners were identified with the Neh1 or CT domains, it is possible that these proteins may act as conditionally competitive partners. Therefore, insight into the rank order of partner interaction strength, under basal and derepressed conditions, may aid prediction of conditional changes in partner preference. Significantly, a large proportion of new NRF2 partners were transcription factors (TFs), eight of which (ATF3, CEBPG, CREBL2, CREBZF, DDIT3, FOSB, FOSL2, TEF) belong to the same bZIP family as NRF2 (Fig. 1C); thus expanding the spectrum of NRF2 heterodimers, which have the potential to modulate ARE transactivation and target gene expression profiles.
Two complementary assays were used to compare the relative strength of NRF2-partner interactions. Binding affinities (in vivo K d values) for KEAP1, MAFG and 23 novel NRF2 partners were first measured in HEK293T cells using fluorescence cross-correlation spectroscopy (FCCS) [41,42] (Fig. 2A). Values ranged over two orders of magnitude, with KEAP1 being the strongest binding partner, possessing an in vivo K d of 1148 nM (95% CI, 776-1520 nM). This value is higher than the published in vitro K d surface values of 167 and 580 nM for interactions between the NRF2 Neh2 domain or Neh2 peptide and KEAP1, respectively [43,44]. This difference may reflect the effects of co-complex or competitive binding events in vivo [42]. A large proportion of the partners, including KEAP1 and MAFG, have in vivo K d values of 1000-10,000 nM, however five NRF2 partners (CEBPG, FOSL2, CREBL2, EIF3J, and ELF3) displayed very low binding affinities >10, 000 nM ( Fig. 2A, Supplementary Fig. S2-3). To re-assess these trends, the high-sensitivity dual luminescence-based co-immunoprecipitation (DULIP) method [45] was used to provide an independent measure of the relative strength and range of novel partner interactions. Six known binary NRF2 partners (MAFG, MAFK, MAFF, KEAP1, ATF4, and PMF1) and five known indirect or co-complex partners (BTRC, FBXW11, APEX1, COPS7A, JUN) were also included for comparison. As FCCS K d and DULIP K d values were highly concordant, all DULIP K d values were scaled to equivalent molar concentrations via linear regression (Y = 0.213*X + 2.110, ***p ≤ 0.001) ( Supplementary Fig. S4A, Supplementary Methods). Scaled DULIP ex vivo K d values similarly spanned two orders of magnitude, with six novel partners (CEBPG, FOSB, ELF3, FOSL2, DDIT3, and MAP2K6) exhibiting weak binding to NRF2 (scaled DULIP K d > 9000 nM). Significantly, known NRF2 partners exhibited a very similar range of binding affinities, ranging from very strong (KEAP1 and sMAF proteins) to very weak (PMF1 and APEX1) (Fig. 2B). The strong correlation between DULIP and FCCS K d values are shown in Fig. 2C.
Given the dominance of KEAP1 binding under basal conditions, it is likely that other mechanisms of functional modulation will be most relevant when KEAP1/NRF2 interactions are genetically or conditionally disrupted. Therefore, the conditional hierarchy of NRF2-partners, and the potential for competitive dominance by KEAP1 was investigated by performing comparative DULIP assays in both WT and genetically edited KEAP1 − /− HEK293T cells. While the majority of interaction partners exhibit similar scaled DULIP K d values in WT and KEAP1 − /− cells ( Supplementary Fig. S4B), 13 partners exhibited a >1.25-fold difference between scaled DULIP K d values in WT and KEAP1 − /− cells, with seven partners (ARPC2, ELK1, ETV6, CREBZF, ATF4, REL, and EIF3J) exhibiting a loss of affinity and six partners (JUN, CEBPG, PMF1, FOSL2, FOSB, and BTRC) showing a gain of affinity for NRF2 in KEAP1 − /− cells ( Supplementary Fig. S4C). The resulting hierarchical binary NRF2 interaction network is shown in Fig. 2D. Considering these results, we propose that partners exhibiting lower K d values in KEAP1 − /− cells may be competitively inhibited from binding NRF2 under basal conditions. However, these partners could in principle modulate NRF2 responses either: in the nuclear compartment where KEAP1 levels are low; in several forms of cancer, where NRF2/KEAP1 interactions are genetically perturbed; or under conditions of oxidative or xenobiotic stress, when KEAP1 mediated degradation of NRF2 is inhibited. It is also important to note that the rank order of NRF2 partner binding changes under derepressed (KEAP1 − /− ) conditions. This quantitative change in partner preference may provide a mechanism by which the extent or nature of NRF2 transcriptional responses can be conditionally tuned. In contrast, an increase in partner K d values in KEAP1 − /− cells may suggest that under basal conditions, these partners could bind cooperatively with NRF2 and KEAP1 in multiprotein complexes. Quantitative data from this study also supports the current paradigm in which BTRC primarily governs nuclear degradation of NRF2 via CUL1-based proteasomal degradation under derepressed rather than basal conditions.

Phosphomimetic mutation of NRF2 Y576 inhibits ARE transactivation and alters bZIP binding preferences
Although the Neh1 domain of NRF2 is known to contain the bZIP heterodimerization domain, quantitative DULIP assays show that both the Neh1 and Neh3 domains of NRF2 are essential for binding to MAFG ( Supplementary Fig. S5A). This raises the possibility that changes within the Neh3 domain may have broader effects on the binding of other bZIP proteins in the neighbouring Neh1 region (Fig. 3A). A preliminary mutagenesis study of six predicted NRF2 phosphorylation sites S40, T80, S215, S408, T559 & Y576 [46] showed that only the phosphomimetic Y576E mutation significantly inhibited ARE transactivation (Supplementary Fig. S5B and Fig. 3B). As bZIP proteins are known to either facilitate or modulate NRF2 transcription, we performed a series of DULIP assays to investigate the effects of the Y576E mutation on bZIP protein recruitment. Interestingly, this mutation significantly inhibited interactions between NRF2/MAFK and NRF2/MAFF, while binding to MAFG was unaffected (Fig. 3C). As eight additional NRF2 bZIP-partners were identified in this study, the effects of the Y576E mutation on these interactions was also tested. While interactions with CREBZF, ATF3, CEBPG, and FOSB all remained unchanged, binding to TEF, CREBL2, and DDIT3 were significantly increased by the Y567E mutation (Fig. 3D).
These quantitative results provide the first evidence that the Neh3 domain of NRF2 is essential for heterodimerization with sMAF proteins. While the Neh3 domain may be required to maintain conformational integrity of the Neh1 region, our data indicates a possible new mechanism by which Y576 phosphorylation in Neh3 could change the rank order of bZIP partner binding, thereby providing a means by which NRF2 activity and target gene expression may be conditionally tuned.

Functional effects of new NRF2 partners on ARE transactivation
Due to the dominant nature of the KEAP1/NRF2 interaction, the ability of new NRF2 partners to modulate ARE transactivation was analysed under both basal and derepressed conditions, using a series of CRISPR/Cas9-edited NRF2 gain-of-function (NRF2 T80K and KEAP1 − /− ) and loss-of-function (NRF2 − /− ) cell-lines. Evaluation of NRF2 levels following treatment with the NRF2 activator sulforaphane (SFN) or the proteasomal inhibitor MG132 confirmed that expression of NRF2 was abolished in NRF2 − /− cell-lines but enhanced in NRF2 T80K/T80K , NRF2 T80K/and KEAP1 − /− cells (Fig. 4A-B). Interestingly, comparative qPCR analysis of several NRF2 target genes under unstimulated conditions revealed differences in expression of HMOX1, FTH1, and AKR1C1, but not GCLM or NQO1, between KEAP1 binding-impaired (NRF2 T80K ) and KEAP1 − /− cell-lines (Fig. 4C-G). This shows that KEAP1 can influence ARE transactivation, and by extension NRF2 target gene expression, by some mechanism that is independent of its ability to bind NRF2. As KEAP1 is an E3 substrate adaptor these effects could result from changes in the degradation of other target proteins, which directly or indirectly affect the nature of NRF2 transcriptional responses. Having confirmed the functional phenotypes of each cell-line, ARE-luciferase reporter assays were performed to compare changes in ARE transactivation in response to ectopic expression of NRF2 interaction partners in HEK293T WT, NRF2 T80K/T80K , NRF2 T80K/-, KEAP1 − /− and two NRF2 − /− cell-lines. In total, the effects of 26 partners were examined in each of the six cell-lines. The resulting data was fitted using a linear mixed model; with the presumption that the effects of partners dependent on derepression of NRF2 would be significantly enhanced in the NRF2 T80K and KEAP1 − /− cell-lines, while the effects of partners specifically dependent on NRF2 would be dampened in the NRF2 − /− cell-lines.
Nine binary NRF2 partners (MAFK, MAFF, MAFG, ETV4, FOSL2, CEBPG, ARPC2, EIF3J, and CREBZF) significantly enhanced AREluciferase activity in HEK293T WT cells, while six partners (ATF3, TBP, NFAT5, TEF, RELA, and ELF3) caused inhibition compared to empty vector controls (Fig. 4H). ETV4 and CEBPG both exhibited a greater upregulation of NRF2 activity under derepressed conditions while TBP and NFAT5 caused greater downregulation in both NRF2 T80K and KEAP1 − /− cells. Significantly, RELA impaired ARE-luciferase transactivation in all three derepressed NRF2 cell-lines, demonstrating a striking conditional dominance under disease-like conditions. While ATF3 downregulated ARE-luciferase transactivation in both WT and KEAP1 − /− cells, this effect was absent in NRF2 T80K cells. This again emphasises the functional differences between KEAP1 ablation and  selective genetic derepression of NRF2. Analysis of ARE perturbation in NRF2 − /− cell-lines shows that the effects of ATF3, TBP, TEF, and ELF3 are all NRF2-dependent (Fig. 4I). Interestingly, six partners (KPNA4, MAP2K6, TADA2A, KPNA2, DDIT3, and KEAP1) that did not perturb ARE-luciferase activity in WT cells, did induce significant changes in derepressed (NRF2 T80K and/or KEAP1 − /− ) cells. Of these, KPNA2, KPNA4, MAP2K6, and TADA2A all enhanced ARE-luciferase activity. Despite being one of the weakest NRF2 binding partners MAP2K6 displaying a strong dependence on disease-related NRF2 derepression; with effects being seen in all three derepressed NRF2 cell-lines. Significantly, ectopic expression of KEAP1 increased ARE-luciferase transactivation in both T80K cell-lines, confirming that KEAP1 can indirectly enhance ARE transactivation, via an NRF2-independent mechanism. These functional screens demonstrate that the new high-density binary NRF2 interactome contains novel activators and inhibitors of ARE transactivation, as well as a subset of partners that specifically exhibit NRF2-dependent effects under disease-related conditions (Fig. 4J). Finally, it is important to note that several of the partners that modulate ARE-luciferase transactivation correspond to very weak NRF2 binding partners.

Predictive analysis of the expanded human NRF2 interactome
Functional annotation of the expanded binary interactome was performed using Gene Ontology Biological Processes (GO BP) with biological processes being grouped according to functional class, using the Cytoscape ClueGO plugin [47] (Fig. 5A). In total, 76 GO BP terms (grouped into 18 functional classes) were significantly enriched; including four previously identified themes (ER stress, cellular response to external stimuli, cellular transcription, and cell cycle/DNA damage response) and two additional themes (nuclear import, and cell differentiation/development) identified in this study. Considering novel partners that significantly alter ARE transactivation, ATF3 and DDIT3 (both NRF2 repressors) were over-represented in GO BP terms involved in ER stress, thus highlighting a novel intersection between oxidative and ER stress responses (Fig. 5A).
An analysis of disease association was also performed on the new expanded human NRF2 interactome. In total 334 gene-disease associations were identified, including 130 (39.9%) from 22 (51.2%) newly identified NRF2 partners (Fig. 5B, Supplementary Table S3). Six of the new NRF2 partners associated with known NRF2-disease phenotypes ( Table 1, Supplementary Fig. S6A), and 13 disease phenotypes involved four or more binary NRF2 partners (Table 2). Given the potential for functional co-operativity between NRF2 partners, the DisGeNET network was merged with a PPI network constructed from five publicly available human PPI databases (BioGRID [33,34], IntAct [48], HPRD [49], HuRI [50], and InnateDB [51]). Disease-centric sub-networks involving two or more NRF2 partners were then extracted, excluding partner pairs with known physical association ( Supplementary  Fig. S6B). Interestingly, STAT3 and ATM were jointly associated with five different NRF2 disease phenotypes (pancreatic neoplasm, squamous cell carcinoma, stomach neoplasms, prostatic neoplasms, and mammary neoplasms), suggesting possible mechanistic similarities in the pathophysiology of these diseases.

Conclusions
Although NRF2 is known to play a major role in redox homeostasis, the complex mechanisms by which NRF2 responses are regulated cannot be explained by current knowledge of well characterised pathway components. This study represents the first large-scale quantitative analysis of the human NRF2 protein interaction network. It was designed to provide a broader insight into the diverse range of cellular components that interact directly with the different functional domains of NRF2. Using a combination of high stringency protein interaction assays we identified 46 novel binary NRF2 partners (Fig. 6). Significantly, >77% of the novel partners tested were shown to effect NRF2 driven ARE-transactivation. Quantitative analysis of protein interactions enabled the rank ordering and strength of partner interactions to be defined, both in vitro and in live cells. These data confirm the dominance of KEAP1 and sMAF proteins under basal conditions, but reveal an emerging model of regulation in which changes that disrupt binding of dominant basal partners (KEAP1 and sMAF proteins) allow weaker binding partners to interact with NRF2, thereby facilitating changes in transcriptional activity. In addition, we provide the first evidence that the NRF2 Neh3 domain is essential for recruitment of MAFG, and demonstrate that a phosphomimetic Y576E mutation not only inhibits ARE transactivation, but also selectively inhibits interactions with MAFK and MAFF, while enhancing interactions with other bZIP partners. Finally, we present quantitative data to show that KEAP1 can influence ARE transactivation in an NRF2 independent manner.
These data significantly expand the range of candidate proteins that have the potential to modulate NRF2 responses. Together they support an emerging paradigm that NRF2 acts as a dynamic sensor to facilitate communication and balance between multiple conditional stressresponses. This new high-density binary NRF2 interactome provides a resource to inform future investigation into the complex conditional regulation of NRF2 responses, and the development of better rationally designed strategies for therapeutic intervention or drug sensitisation.

Lead contact and materials availability
Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Christopher M Sanderson (, cmsand@liverpool.ac.uk).

Generation of bait and prey yeast two-hybrid clones
Generic primers containing Gateway flanking sequences were designed to amplify any ORF from pDONR201/221/223 vectors by proof-reading PCR to facilitate in vivo gap repair cloning into Gateway compatible BamHI linearised bait (pGBAD/E-B) and prey (pACTBD/E-B) Y2H vectors [40]. This approach allowed all clones where possible to be expressed as both bait and prey fusions. Gap repair reactions were performed as previously described [52]. Y2H assays utilised the bait PJ69-4a yeast strain and prey PJ69-4α strain.

Yeast two-hybrid matrix screens
NRF2 clones were systematically mated against the array of predicted interaction partners in both bait/prey orientations where possible. The Neh4 and Neh5 domains of NRF2 showed a propensity to auto-activate in Y2H assays, however this was not the case for the combined Neh4+5+7 domains, when expressed as a Y2H prey construct. Haploid yeast were initially mated on YPAD media for 24 h prior to replication onto SD-Leu/-Trp diploid selection media and grown for 48 h. Protein-protein interaction was detected following the replication of diploid yeast onto the above media also lacking adenine or histidine (supplemented with 2.5 mM 3-AT). Activation of reporter genes ADE2 and HIS3 and thus growth, indicative of protein-protein interaction, was scored over a period of 14 days. Only interactions reproducibly observed in two independent assays were scored as

Co-immunoprecipitation
HEK293T cells were co-transfected with 1 μg each of EGFP-tagged NRF2 and mCherry-tagged binary partner vectors for 24 h. Cells were lysed in NP-40 lysis buffer (0.5% NP-40 substitute, 25 mM Tris-HCl pH 7.5, 150 mM NaCl, 2 mM MgCl 2 , 1 mM EGTA, protease and phosphatase inhibitor cocktail), following which clarified protein lysates were incubated with GFP-Trap agarose beads (ChromoTek) for 2 h at 4 • C, with rotation. Beads were collected by centrifugation and washed thrice in NP-40 lysis buffer, then once in IP wash buffer (10 mM Tris pH 7.5, 2 mM MgCl 2 ). Proteins were eluted directly with NuPAGE LDS Sample Buffer (Invitrogen) for 10 min at 80 • C before being resolved by SDS-PAGE and Western blotting.

FCCS quantification
FCCS data collection was performed as previously described [41], using a 63x/1.4 Plan-Apochromat objective on a Zeiss LSM 780. Briefly, HEK293T cells were seeded overnight in four-compartment glass bottom dishes (2.5 × 10 4 cells/compartment) and transfected with 25 ng of EGFP-tagged NRF2 and mCherry-tagged binary partner for 24 h. EGFP fluorescence was excited using a 488 nm Argon laser and emission collected between 500 and 530 nm, while mCherry fluorescence was excited using a 561 nm DPSS diode and emission collected between 590 and 640 nm. Data from individual cells were collected via 5 × 2 s measurements either in the cytoplasm or nucleus dependent upon partner localisation using 0.15-0.3% laser power to minimise bleaching and a suitable count rate of approximately 1 kHz count rate per molecule (CPM). The autocorrelation G auto (τ) and cross-correlation G cross (τ) functions were calculated using the ZEN 2012 software. The calculated autocorrelation and cross-correlation values were fitted to a two-component diffusion model with or without triplet state correction respectively using the Levenberg-Marquardt algorithm in MATLAB optimisation toolbox. The in vivo K d values were obtained from nonlinear fitting by plotting the fraction of bound EGFP-tagged and mCherry-tagged proteins as a function of free mCherry-tagged or EGFP-tagged proteins respectively, and averaging the fitted values.

ARE-luciferase activity screens
HEK293T cells were co-transfected for 24 h using 100 ng of 8XARE-FL, 5 ng of Renilla luciferase (pRL-SV40), as well as 100 ng of an expression vector containing the binary partner in a 96-well plate. Cells were lysed in Passive Lysis Buffer (Promega) and transferred to a 96-well plate before firefly and Renilla luminescence were measured sequentially on a GloMax Multi-Detection System. Log2-normalised luciferase fold change was fitted to a linear mixed model with 'condition' and 'cellline' as fixed effects and 'plate' as a random effect using the 'nlme' [53] and 'lmerTest' [54] R packages. p-values were calculated by Satterthwaite's method for unequal variance and corrected for multiple testing using Bonferroni's method.

Network construction
Protein interaction data was extracted and merged from BioGRID

Data and code availability
MATLAB custom analysis code is available from the lead contact upon reasonable request. NRF2 protein-protein interaction data have been submitted to the IMEx (http://www.imexconsortium.org) consortium through IntAct [48] (IMEx ID: IM-27648).

Author contributions
JP, AHP, JB contributed experimental data. CMS, AHP, JP, JB, DM, NH contributed to experimental design. US, EW, JDW contributed to clone generation and development of assay reagents. All authors contributed to manuscript preparation.

Declaration of competing interest
The authors declare no competing interest.