Identification of Nrf2-responsive microRNA networks as putative mediators of myocardial reductive stress

Although recent advances in the treatment of acute coronary heart disease have reduced mortality rates, few therapeutic strategies exist to mitigate the progressive loss of cardiac function that manifests as heart failure. Nuclear factor, erythroid 2 like 2 (Nfe2l2, Nrf2) is a transcriptional regulator that is known to confer transient myocardial cytoprotection following acute ischemic insult; however, its sustained activation paradoxically causes a reductive environment characterized by excessive antioxidant activity. We previously identified a subset of 16 microRNAs (miRNA) significantly diminished in Nrf2-ablated (Nrf2−/−) mouse hearts, leading to the hypothesis that increasing levels of Nrf2 activation augments miRNA induction and post-transcriptional dysregulation. Here, we report the identification of distinct miRNA signatures (i.e. “reductomiRs”) associated with Nrf2 overexpression in a cardiac-specific and constitutively active Nrf2 transgenic (caNrf2-Tg) mice expressing low (TgL) and high (TgH) levels. We also found several Nrf2 dose-responsive miRNAs harboring proximal antioxidant response elements (AREs), implicating these “reductomiRs” as putative meditators of Nrf2-dependent post-transcriptional regulation. Analysis of mRNA-sequencing identified a complex network of miRNAs and effector mRNAs encoding known pathological hallmarks of cardiac stress-response. Altogether, these data support Nrf2 as a putative regulator of cardiac miRNA expression and provide novel candidates for future mechanistic investigation to understand the relationship between myocardial reductive stress and cardiac pathophysiology.

and extracellular matrix degradation [1][2][3] . However, the findings from clinical studies have largely discredited the efficacy of antioxidant therapies 4,5 . Furthermore, our laboratory has identified the presence of a "reductive stress" wherein aberrant induction of antioxidant response element (ARE)-dependent antioxidant genes produces pathological cardiac hypertrophy and dysfunction [6][7][8] . While the deleterious consequences of reductive stress appear highly-conserved 9,10 , the transcriptional and post-transcriptional mechanisms of the myocardial redox milieu remain unknown.
Among the mechanisms known to regulate postnatal heart function, microRNAs (miRNAs) are a class of short (~ 22 nucleotide) RNAs that post-transcriptionally regulate mRNA stability and translational efficiency, often in a tissue-specific manner 11 . While miRNAs are necessary for physiologic cardiac function and development 12,13 , many have been found to be dysregulated in the failing heart [14][15][16] . Specifically, miRNAs have been directly linked to the structural and functional deficits in cardiac function 17 . Nevertheless, it remains unclear whether ARElinked miRNAs contribute to cardiac pathogenesis.
As a transcriptional activator of cis-regulatory AREs 18 , nuclear factor erythroid 2-related factor 2 (Nfe2l2, a.k.a. Nrf2) plays a critical role in regulating cardiac redox status. Transient Nrf2 signaling is cardioprotective immediately following ischemic insult 19 , but chronic transactivation of AREs causes reductive stress and cardiac dysfunction 20,21 . We have recently shown that Nrf2 deficiency (Nrf2 −/− ) inhibits the expression of several miRNAs in the heart 22 , but the genome-wide impact of reductive stress on miRNA expression remains unknown.
In this investigation, we identify a miRNA signature for reductive stress to gain insight into potential biomarkers and/or effectors of this novel pathological phenomenon in the heart. The cardiomyocyte-specific and constitutively-active Nrf2 transgenic mouse model (CaNrf2-Tg) was used to conduct a multi-omics analysis of Nrf2-dependent and ARE-bearing miRNAs, which we term "reductomiRs". Our use of both mRNA-seq and small RNA sequencing (miRNA-seq) in caNrf2 low (TgL) and high-expressing (TgH) mouse lines reveals a distinct signature of transgenic Nrf2 dose-responsive miRNAs linked to a number of suppressed cardiac genes under pro-reductive and reductive stress conditions 23 . Collectively, this analysis uncovers several novel miRNA candidates for which future mechanistic studies will investigate the interplay between post-transcriptional regulatory responses and redox state in the myocardium of Nrf2-Tg mice.

Methods
Animals. Our method for establishing the cardiac-specific constitutively active Nrf2 transgenic mouse model (caNrf2-Tg) has been described previously 23 . Briefly, cDNA encoding a truncated Nrf2 protein lacking the Neh2 domain was ligated into an α myosin heavy chain (αMHC) expression vector, the plasmid backbone was digested, and the αMHC-caNrf2 insert was used for pronuclear injection. Transgenic low (TgL) and transgenic high (TgH) founders were determined using caNrf2 primer sets in real-time qPCR which compared relative transgene expression to endogenous Nrf2 mRNA, and transgenic mice were back-crossed onto the C57BL/6J background for six generations. For expression analyses, sex-matched male and female TgL, TgH and nontransgenic (NTg) littermate heart tissue was harvested at 6 months of age. All the animal studies including noninvasive cardiac imaging were conducted in accordance with the Guide for Care and Use of Laboratory Animals developed by the National Research Council at the National Institutes of Health (NIH). The Institutional Animal Care and Use Committee (IACUC#14-10160) at the University of Alabama at Birmingham has approved the study, which was carried out in compliance with the Guidelines for Animal Research Reporting of In Vivo Experiments (ARRIVE).

GSH-NEM immunofluorescence assay.
To measure the redox status of Tg & NTg mouse hearts, we performed immunofluorescence using anti-GSH staining 24 . The mouse was anesthetized (2% isoflurane) and the heart was perfused in situ with 10 mM PBS-N-Ethylamaleimide (NEM) before sacrificing. The heart tissue sections (5 µM) were incubated in ethanol containing 10 mM NEM for 30 min and washed three times with PBS prior to blocking. Slides were incubated at 4 °C overnight with the primary antibody mouse anti-Glutathione:N-Ethylamaleimide adduct (1:500 v/v; EMD Millipore Corp., Billerica, MA, US). Following the incubation with the secondary Ab conjugated to Alexa-Flour 488 and washing, imaging was performed using a Fluorescence microscope (Olympus BX43). The intensity of green fluorescence representing glutathione (GSH) was calculated using Image-J (NIH) in 6 images from 3 hearts per group.
Non-invasive cardiac imaging. Two-dimensional trans-thoracic echocardiography was performed in NTg and Tg mice at 6-7 months of age using the Vevo2100 system (Fujifilm Visual Sonics Inc., Ontario, Canada). In brief, animals were prepared for imaging as previously described 25 . Using a 38-MHz mechanical transducer, strain analysis in parasternal long-axis were acquired under isoflurane (1.5-2.0%) anesthesia. Structural and functional (systolic and diastolic) measurements were obtained using Vevo lab 3.1 software.
RNA and miRNA sequencing analysis. Details of the R coding scripts and other bioinformatics tools used in the current study are available as online Supplemental Methods and GitHub data repository: https:// github. com/ mepep in/ Reduc tomiRs. High-throughput RNA-and miRNA-sequencing and analysis were performed at the University of Utah. Alignment of reads to the mm10 genome was accomplished using NovoAlign (Singapore, Malaysia), raw counts generated using Samtools 26 , and differential gene expression performed using DESeq2 27 (1.18.1) within the R (3.4.2) statistical computing environment as described previously 28 . Due to limited sample sizes (n = 3-4), dispersion estimates were first determined via maximum-likelihood, assuming that genes of similar average expression strength possess similar dispersion, as previously described 27 . Gene-wise dispersion estimates were then shrunken according to the empirical Bayes approach, providing normalized count data for genes proportional to both the dispersion and sample size. Differential expression was determined from www.nature.com/scientificreports/ normalized read counts via Log 2 (fold-change) using the Wald test followed by Bonferroni-adjusted P-value (i.e. Q-value) for each aligned and annotated gene. Statistical significance was assessed with unpaired two-tailed Bonferroni-adjusted P-value (Q-value) < 0.05 used to identify differentially-expressed candidate genes, with a lower-stringency statistical threshold used for pathway enrichment (P < 0.01). Differentially expressed genes (DEGs) were generated with biological significance assumed when |Log 2 FoldChange|> 1 with normalized count sum > 1. Data visualization. Functional and network gene set enrichment analysis (GSEA), along with curated literature-supported candidate upstream regulators, were performed using Enrichr, an interactive web-based tool for compiling multiple bioinformatics databases 29 . Within this software, Wikipathways 30 and KEGG 31 enrichment as well as ChIP Enrichment Analysis (ChEA) were done both on the individual datasets and as a combined comparative analyses to determine overlapping enriched pathways and associated transcriptional regulators, respectively. Heatmap generation and hierarchical clustering were performed using pheatmap package (1.0.8) within R, and VennPlex 32 was used to create the Venn diagrams and determine overlapping gene lists.

RNA isolation and real
Statistics. All data are represented as mean ± standard error of the mean (SEM) unless otherwise indicated.
Statistical significance was determined using unpaired Student's t-tests or, where appropriate, one-way (ANOVA) with Tukey post-test for multiple comparisons. Statistical analyses and data visualization were completed using GraphPad Prism (GraphPad Software, San Diego, CA) and R software, version 3.4.2 (R Foundation for Statistical Computing, Vienna, Austria). Statistical significance was assigned at the P < 0.05 level for two-way comparisons.

Consent for publication. All authors verified the content and approved the final version for submission
and publication.

Results
Evidence of reductive stress and pathologic cardiac remodeling in CaNrf2-transgenic mice. The concept myocardial reductive stress was first described in a mouse model of human cardiac disease 8 . However, the underlying molecular mechanisms remained unknown until the recent development of unique mouse models to study this phenomenon 25,33 . We recently described that constitutive activation of the transcription factor, Nrf2/NFE2L2, resulting in a dose-dependent increase of glutathione and other antioxidants, establishes a reductive and hyper-reductive redox condition in the heart 33 . Furthermore, a chronic hyperreductive condition (i.e. "reductive stress") was shown to cause pathological cardiac remodeling and diastolic dysfunction 25 . Here, we sought to identify a redox-responsive miRNA signature associated with reductive-and hyper-reductive-conditions in caNrf2 transgenic low (TgL) and high (TgL) hearts. Quantitative genotyping in TgL and TgH mice indicated dose-dependent transgene expression (respective Ct values) when compared to NTg, which did not amplify a PCR product for the primers that recognize only the truncated transgene (caNrf2) (Fig. 1A). Interestingly, qPCR using a primer that recognizes endogenous mouse Nrf2, but not caNrf2, also showed a dose-wise expression pattern consistent with reports of Nrf2 transcriptional auto-up-regulation ( Fig. 1B) 34 . As Nrf2 is well-characterized to promote glutathione biosynthesis 35 , we performed immunofluorescent staining of glutathione-N-ethylmaleimide (GSH-NEM) adducts in TgL and TgH to confirm the biochemical effects of dose-dependent transgene expression (Fig. 1C). These results confirm our previous studies employing kinetic-based spectrophotometric measurements for glutathione in CaNrf2-Tg hearts 33 . To determine the consequences of Nrf2 expression on cardiac function, echocardiographic assessments (strain analysis) were performed (n = 10/group), which provided clear evidence of diastolic dysfunction in TgL and TgH relative to NTg mice (Fig. 1D). Taken together, these results suggest that Nrf2 expression leads to both reductive conditions and pathological cardiac remodeling.

Nrf2-dependent cardiac mRNA expression.
In a recent study, we analyzed mRNA-sequencing data from caNrf2 Tg mice to identify differences associated with the physiologic state of myocardial reductive stress 36 . However, this report lacked a comparative analysis of pooled vs. dose-responsive DEGs and the mRNA coregulatory relationships with miRNAs. Thus, we analyzed differential expression for caNrf2-Tg vs. NTg (Supplemental Table S9) and TgL vs. TgH (Supplemental Table S10) to relate mRNA dysregulation with miRNA expression patterns. Principal Component Analysis (PCA) (Fig. 3A) and correlation heatmap (Fig. 3B) of DEGs revealed a striking pattern of mRNA expression of TgH compared to TgL mice. Consistent with our functional characterization (Fig. 1) and those of our earlier work 37 , GSEA of "dose-only" mRNA revealed pathways enriched for cardiac pathophysiology while "transgene-only" DEGs implicated noncanonical -possibly indirecteffects of transgenic Nrf2 expression including "longevity" and "FoxO signaling" (Fig. 3C). Notably, upregulated "transgene-dose" effect DEGs reflected quintessential functions of Nrf2 signaling such as "Glutathione" and "Xenobiotic Metabolism, " thereby suggesting that this subset of genes represent putative downstream targets of caNrf2 (Fig. 3C). Owing to its role as a positive transcriptional regulator, we were surprised to find that CaNrf2 decreased the expression of 206 "transgene-dose" DEGs (Fig. 3C). Because miRNAs silence gene expression at the post-transcriptional level, we hypothesized that this distinct DEG subset could represent miRNA targets induced downstream of Nrf2 transactivation. These DEGs thus became the focus of our study as putative targets of negative regulation via Nrf2-responsive miRNA, or cardiac "reductomiRs".  Fig. 4A. To identify DEmiRs containing upstream AREs ("ReductomiRs"), we used the position-weighted matrix (PWM) developed by Malhotra et al. 38 within PWMscan 39 to identify genome-wide (Mus musculus, mm10) Nrf2 consensus elements. The resulting genomic coordinates were then used to locate a total of 142 AREs within the promoter most-proximal to the coding sequence (+ 5 kB) of 40 upregulated Nrf2-dependent DEmiRs (Fig. 4B). From this, approximately half (22/40) of the DEmiRs were found to exhibit a transgene dose-responsive pattern of expression in caNrf2 TgL and TgH hearts. To determine whether ARE bearing doseresponsive miRNAs accounted for downregulated mRNAs within the transgene-dose mRNA subset, we used multiMiR 40 as a computational tool to systematically compile predicted mRNA targets from the following 8 algorithms: DIANA-microT, ElMMo, MicroCosm, miRanda, miRDB, PicTar, PITA, and TargetScan. Targets were then compiled using cross-algorithm validation to predict miRNA-mRNA interactions 41 . From this analysis, 19 miRNAs exhibited sequence complementarity to seed sequences within 61 down-regulated DEGs. Gene-set enrichment of these DEGs using KEGG pathways disproportionately represented "hypoxia-inducible factor-1α signaling" and "cardiac muscle contraction" (Fig. 4C). Enrichment analysis using the ARCHS4 42 gene expression atlas revealed a cardiac-specific signature among these DEGs (Fig. 4D). Taken together, these preliminary findings implicate "reductomiRs" as possible Nrf2-responsive mediators of myocardial reductive stress. www.nature.com/scientificreports/ Integrative Nrf2-responsive miRNA-mRNA functional network. It became apparent from our in silico analysis that multiple miRNAs shared one or more putative downstream mRNA targets. To understand these interactions further, we generated a bipartite network of the 19 "dose-responsive" reductomiRs (P < 0.01) found to possess a downstream seed sequence among the suppressed mRNAs (Fig. 5A). The resulting network of interactions suggested that Nrf2-mediated miRNAs display a varied capacity to influence mRNA transcript abundance in caNrf2-Tg hearts. Whereas miR-301a-5p and miR-31-3p were predicted to target relatively few mRNAs, the seed sequences of miR-449a-5p, miR-96-5p, miR-298-5p, miR-491-5p, and miR-671-5p matched a larger number of cardiac-specific mRNAs within this network. Once putative miRNA-mRNA interactions were developed, the curated database from Ingenuity Pathway Analysis (IPA) was used to incorporate known miRNA-mRNA interactions into a caNrf2-dependent reducto-miR-mRNA network, which we then compiled to identify nodal DEGs (Fig. 5B). This approach identified several key mRNAs, which together comprise established mediators of pathological cardiac remodeling: sarcomeric genes myomesin (Myom2), myosin light chain kinase 4 (Mylk4), voltage gated calcium and chloride channels subunits (Cacnad2d2, Clcn1), the striated muscle isoform subunit VIa of cytochrome C oxidase (Cox6a2), and adrenergic alpha1a receptor (Adra1a). Furthermore, reductomiR targets of caNrf2 were predicted to mediate the negative regulation of SET and MYND domain containing 1 (Smyd1) and myocyte enhancer factor 2A (Mef2a), two transcriptional regulators essential for cardiac development (Fig. 5B) [43][44][45] . Although future work should verify direct post-transcriptional regulatory interaction and resultant changes in proteins encoded by these mRNAs, these data implicate miRNAs as putative mediators of Nrf2 signaling and, ultimately, myocardial reductive stress effectors.

Discussion
Since their initial discovery in 1993 46 , two and a half decades of research have shown miRNAs to be central determinants of human biology and disease. As such, miRNA therapies are continuously being patented, and a few have entered clinical trials 47 . The heart's susceptibility to miRNA dysregulation is well-known, with aberrant miRNA induction found to directly cause pathological cardiac remodeling 48,49 . Nevertheless, endogenous miRNA activity is essential for cardiac morphogenesis 12,13 and postnatal functions 50,51 . While efforts to identify candidate miRNAs in human patients [14][15][16] have led to the identification of key miRNAs governing calcium handling, extracellular matrix remodeling, and vascularization 52 , their role in the regulation of myocardial reductive potential has remained elusive. In this context, our recent discovery of suppressed miRNAs in Nrf2 knockout (Nrf2 −/− ) hearts supports the notion that Nrf2-dependent miRNA expression offers a mechanistic link between redox disequilibrium and cardiac pathophysiology 22 .
In the absence of oxidative conditions, Nrf2 activation is generally inhibited through its interaction with the cytosolic repressor protein Keap1 53-55 . Thus, despite its cardioprotective roles in response to pathological oxidative stress stimuli, such as pressure overload 56 and ischemia-reperfusion injury 19 , the adaptive phase of Nrf2 signaling is generally transient in an unstressed adult myocardium. Our laboratory has shown that this basal restriction of Nrf2 is required, since its sustained activation leads to myocardial reductive stress characterized by accumulation of reducing equivalents and attenuation of physiological oxidant signaling 8,20,21,37 . In addition to our observations in mouse models, others have reported that human germline gain-of-function mutations in Nrf2 produce early cardiac pathologies 57 . Sariam and colleagues have reported that a subset of human heart failure patients exhibit a hyper-reductive redox state 9 . This suggests that divergent patient populations may contribute to the insufficiency and/or adverse effects of exogenous antioxidant supplementation in cardiac patients 4,5 .
In the present investigation, we observed that constitutive Nrf2 activation led to profound and dose-responsive miRNA changes in the mouse heart. We specifically found multiple patterns of Nrf2-dependent miRNA expression, with 165 "transgene-only" miRNAs, 22 "transgene-dose" miRNAs, and 40 "dose-only" miRNAs were differentially expressed in caNrf2-Tg mice. Our in silico search for ARE-containing caNrf2-dependent miRNAs, or "reductomiRs", identified 40 candidate miRNAs likely to be trans-activated by cardiac Nrf2, 22 of which also exhibited dose-responsive induction in TgH and TgL. Lastly, our mRNA target prediction identified 61 putative downstream transcriptional targets for 19 (86%) reductomiRs. While a causal role of these Nrf2-mediated reductomiRs remains untested, several miRNAs have already been implicated in the failing human myocardium (e.g. miR-326-3p and miR-214-5p) and plasma (e.g. miR-671-5p) 58,59 . Therefore, this "reductomiR" signature may provide important future insights into redox disequilibrium during cardiac pathophysiology.
Owing to the bioinformatic tools used in this study, mechanistic studies are still needed to confirm our analyses. Since RNA was isolated from heart homogenates, we cannot exclude the possibility that miRNA could originate from non-myocyte cells. Furthermore, future work should differentiate between indirect and effector www.nature.com/scientificreports/ reductomiRs in myocardial reductive stress as the current investigation does not provide causal evidence that ARE-harboring miRNAs contribute to cardiac pathology or altered redox state. Nevertheless, our data strongly implicate Nrf2 as regulator of the myocardial transcriptome, from which direct and indirect regulation of numerous mRNAs and miRNAs are likely to occur. While Nrf2-independent changes in the caNrf2-Tg model remain unexplained, we demonstrate a strong link between Nrf2-mediated miRNA activity and post-transcriptional regulation of genes involved in pathological cardiac remodeling. The www.nature.com/scientificreports/ current study therefore identifies a novel network of Nrf2-responsive miRNAs that amy represent mediators of myocardial reductive stress.