Functional predication of differentially expressed circRNAs/lncRNAs in the prefrontal cortex of Nrf2-knockout mice

In the central nervous system, nuclear factor erythroid-2-related factor 2 (Nrf2) protects neurons from oxidant injury, thereby ameliorating neurodegeneration. We explored the key circular RNAs (circRNAs) and long non-coding RNAs (lncRNAs) involved in Nrf2-induced neuroprotection. We used microarrays to examine the circRNAs (DEcircRNAs), lncRNAs (DElncRNAs) and mRNAs (DEmRNAs) differentially expressed between Nrf2 (+/+) and Nrf2 (-/-) mice and identified DEcircRNA/DElncRNA-miRNA-DEmRNA interaction networks. In total, 197 DEcircRNAs, 685 DElncRNAs and 356 DEmRNAs were identified in prefrontal cortical tissues from Nrf2 (-/-) mice. The expression patterns of selected DEcircRNAs (except for mmu_circ_0003404) and DElncRNAs in qRT-PCR analyses were generally consistent with the microarray analysis results. Functional annotation of the DEmRNAs in the DEcircRNA/DElncRNA-miRNA-DEmRNA networks indicated that five non-coding RNAs (mmu_circ_0000233, ENSMUST00000204847, NONMMUT024778, NONMMUT132160 and NONMMUT132168) may contribute to Nrf2 activity, with the help of mmu_circ_0015035 and NONMMUT127961. The results also revealed that four non-coding RNAs (cicRNA.20127, mmu_circ_0012936, ENSMUST00000194077 and NONMMUT109267) may influence glutathione metabolism. Additionally, 44 DEcircRNAs and 7 DElncRNAs were found to possess coding potential. These findings provide clues to the molecular pathways through which Nrf2 protects neurons.

Non-coding RNAs (ncRNAs), including circular RNAs (circRNAs), long non-coding RNAs (lncRNAs) and microRNAs (miRNAs), are the dominant products of eukaryotic transcription, and function directly at the In previous studies, we examined the effects of Nrf2related circRNAs and lncRNAs on the brain by detecting the altered expression of circRNAs, lncRNAs and mRNAs in different brain areas in Nrf2-null (Nrf2 (-/-)) mice, including substantia nigra, corpus striatum and hippocampus [23][24][25]. The prefrontal cortex has often been called the command and control center of the brain, and is mainly responsible for cognitive control, especially the ability to orchestrate thought and action in accordance with internal goals [26]. Thus, in this study, we selected the prefrontal cortex to continue this series of studies in Nrf2 (-/-) mice. First, we performed microarray analyses to identify differentially expressed circRNAs (DEcircRNAs), lncRNAs (DElncRNAs) and mRNAs (DEmRNAs) between prefrontal cortex tissues from Nrf2 (-/-) and Nrf2 (+/+) mice. Then, we constructed DEcircRNA/DElncRNA-miRNA-DEmRNA networks and performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses to annotate the functions of the DEmRNAs in the networks. Additionally, we evaluated the coding potential of the DEcircRNAs/DElncRNAs, the expression correlations between the DEcircRNAs/DElncRNAs and their host mRNAs, and the intersections between the DEcircRNAs/DElncRNAs/DEmRNAs in different brain areas. Our study has provided clues into the molecular pathways whereby Nrf2 protects neurons.

Quantitative real-time polymerase chain reaction (qRT-PCR) validation of DEcircRNA and DElncRNA expression
Considering that DEcircRNAs and DElncRNAs with higher FCs and lower p-values were more likely to be regulated by Nrf2, and that higher expression may enable more accurate microarray detection, we selected DEcircRNAs and DElncRNAs with relatively high FCs, high expression values and low p-value to validate the accuracy of the microarray data using qRT-PCR. We used different screening criteria to screen different RNAs. For DElncRNAs, those with p-values < 0.015, FCs > 2.8 or < 0.15, and Standardized Expression Levels (SELs) > 1 were selected. For DEcircRNAs included in "circBase", those with p-values < 0.05, FCs >1.7 or < 0.55, and SELs > 1 were selected. For DEcircRNAs not included in "circBase", those with p-values < 0.01, FCs >2.5 or < 0.25, and SELs > 1 were selected.
Subsequently, qRT-PCR was used to screen primers, and only six DEcircRNAs and eight DElncRNAs were found to have suitable primers. Thus, three up-regulated DEcircRNAs (mmu_circ_0000233, mmu_circ_0015035 and mmu_circ_0003404), three down-regulated DE circRNAs (mmu_circ_0012936, mmu_circ_0008393 and cicRNA.20127), six up-regulated DElncRNAs (NON MMUT132168, NONMMUT127961, NONMMU T132160, NR_045161, ENSMUST00000204847 and NONMMUT024778) and two down-regulated DElnc RNAs (NONMMUT109267 and ENSMUST000 00194077) with relatively high FCs, high expression values and low p-values were selected to validate the accuracy of the microarray data using qRT-PCR. As shown in Figure 2, the qRT-PCR results for all the DElncRNAs and DEcircRNAs except for mmu_circ_0003404 were consistent with the microarray results.

Functional annotation of DEmRNAs that were targets of DEcircRNAs and DElncRNAs
GO and KEGG pathway analyses were employed to annotate the functions of the DEmRNAs in the DEcircRNA/DElncRNA-miRNA-DEmRNA network.

Coding potential of DEcircRNAs and DElncRNAs
Among the 197 DEcircRNAs and 685 DElncRNAs, we found that 44 DEcircRNAs and 7 DElncRNAs had scores > 0 in the Coding-Non-Coding Index [27], suggesting that they may have coding potential (Supplementary Table 1).

Expression correlation analysis between DEcircRNAs/DElncRNAs and their host mRNAs
We next evaluated whether the levels of the DEcircRNAs/DElncRNAs correlated with those of their host mRNAs, in order to determine whether the DEcircRNAs/DElncRNAs altered the expression of their host genes. The levels of 38% of the DEcircRNAs (74 of 197) correlated with the levels of their host mRNAs, with 70% correlating positively and 30% correlating negatively ( Figure 5A and Supplementary Table 2). After we removed 527 intergenic DElncRNAs without maternal genes, we found that the levels of 35% of the DElncRNAs (55 of 158) correlated with the levels of their host mRNAs, with 67% correlating positively and 33% correlating negatively ( Figure 5B and Supplementary Table 2).

Comparison of DEcircRNAs, DElncRNAs and DEmRNAs among the prefrontal cortex, substantia nigra, corpus striatum and hippocampus
Next, we determined the shared DEcircRNAs, DElncRNAs and DEmRNAs among the prefrontal cortex, substantia nigra, corpus striatum and hippocampus, in order to identify common regulatory mechanisms of Nrf2 in these areas. Surprisingly, few circRNAs/lncRNAs were consistently regulated by Nrf2 in all four brain regions. Based on their differential expression between Nrf2 (-/-) and Nrf2 (+/+) mice, no circRNAs/lncRNAs were consistently regulated by Nrf2 in the prefrontal cortex and substantia nigra; only one down-regulated ncRNA (mmu_circ_0010156) was consistently regulated by Nrf2 in the prefrontal cortex, corpus striatum and hippocampus; AGING only one up-regulated ncRNA (ENSMUST00000144818) was consistently regulated by Nrf2 in the prefrontal cortex and corpus striatum; and only two up-regulated ncRNAs (cicRNA.4014 and NR_027831) and one down-regulated ncRNA (ENSMUST00000145790) were consistently regulated by Nrf2 in the prefrontal cortex and hippocampus ( Figure 6A, 6B).

DISCUSSION
The prefrontal cortex is responsible for cognition, emotion, pain and behavior management, and its function is disrupted in various neurodegenerative diseases [28,29]. There is abundant evidence that insufficient Nrf2 signaling is a major cause of nervous system disorders, particularly neurodegenerative diseases [8,[30][31][32]. We hypothesized that Nrf2 could exert its neuroprotective effects by altering the expression of circRNAs and lncRNAs. Although circRNAs and lncRNAs differ significantly in shape and size, both are expressed in tissue-and spatiotemporal-specific patterns and regulate their target genes by similar mechanisms [33][34][35][36]. Since the ceRNA hypothesis was proposed [37], both circRNAs and lncRNAs have been reported to function as miRNA sponges that alter the stability or translation of mRNAs [21,38]. Thus, we used microarray and bioinformatics analyses to detect DEcircRNAs and DElncRNAs between the prefrontal cortex tissues of Nrf2 (-/-) and Nrf2 (+/+) mice.
In addition, four ncRNAs (cicRNA.20127, mmu_circ_0012936, ENSMUST00000194077 and NONMMUT109267) were found to be ceRNAs for GSTM1, GSTM2 and MGST1. These three DEmRNAs were significantly down-regulated in the prefrontal cortex in Nrf2 (-/-) mice and were enriched in the pathway of Glutathione metabolism (KEGG: mmu00480). GSTM1 and GSTM2 encode glutathione S-transferases, which mainly participate in the phase II detoxification of endogenous oxidative stress products and exogenous electrophilic chemical compounds [39]. Reduced GSTs activity has been observed in multiple brain regions and ventricular cerebrospinal fluid from Alzheimer's disease patients shortly after death [40]. Furthermore, the cytoprotective effects of Nrf2 are thought to be due to its strong induction of enzymes involved in glutathione synthesis and recycling [41]. MGST1 is a member of the Membrane-Associated Proteins in Eicosanoid and Glutathione metabolism family, an ancient and diverse family thought to predate the cytosolic family of GSTs [42]. MGST1 detoxifies reactive intermediates such as metabolic electrophile intermediates and lipophilic hydroperoxides through its glutathione-dependent transferase and peroxidase activities [43]. Thus, cicRNA.20127, mmu_circ_ 0012936, ENSMUST00000194077 and NONMMU T109267 may facilitate the effects of Nrf2 on cognitive function by regulating glutathione metabolism.
Interestingly, the results of the GO analysis included immune-related function such as Lymphocyte chemotaxis (GO: 0048247). The central nervous system lacks a classical lymphatic drainage system; however, the discovery of the lymphatic system in the central nervous system in 2015 suggested that central nervous system diseases may be associated with immune system disorders [44][45][46]. Our results support the notion that the lymphatic system is present in the brain, and indicate that its activity may be influenced by the Nrf2 pathway. However, the mechanisms whereby immune cells enter and exit the central nervous system remain poorly understood.
Although circRNAs and lncRNAs commonly alter gene expression as ceRNAs, when they are present at low levels they may not impact their target miRNAs [47]. CircRNAs and lncRNAs may function through other mechanisms, such as binding with proteins [48,49] or being translated into proteins/peptides [50,51]. Interestingly, we identified 44 DEcircRNAs and 7 DElncRNAs with coding potential, in keeping with the view that some ncRNAs may also be mRNAs.

AGING
In addition, some reports have indicated that circRNAs and lncRNAs can alter the expression of their host mRNAs [33]. Cis-acting circRNAs and lncRNAs regulate gene expression in a manner dependent on the location of their own transcription sites, at varying distances from their targets in the linear genome [52,53]. Thus, we analyzed the expression correlations between DEcircRNAs/DElncRNAs and their corresponding host mRNAs. Our results suggested that only a small percentage of DEcircRNAs/DElncRNAs (38% of DEcircRNAs and 35% of DElncRNAs) altered the expression of their host mRNAs. These results were consistent with previous reports indicating that the expression of many circRNAs and lncRNAs is independent of the expression of their related mRNA isoforms [54,55]. However, other studies have demonstrated that circRNAs are typically generated at the cost of their canonical mRNA isoforms [56,57].
Among the DEcircRNAs that were correlated with their host mRNAs, only 30% exhibited negative correlations, suggesting that many DEcircRNAs activate transcription. CircRNA, especially exon-intron circRNAs, have been reported to bind specific ally with U1 small nuclear ribonucleoprotein RNA and promote the transcription of their parental genes in cis, although they may also influence other loci in trans [58,59]. Likewise, the nonrandom positioning of lncRNAs implies that they are functionally associated with their nearby genes [60]. Cisacting lncRNAs have been demonstrated to activate, repress or otherwise alter the expression of their target genes through various mechanisms [52]. In the present study, among the DElncRNAs that were associated with their host mRNAs, 67% exhibited positive correlations and 33% negative correlations.
We also determined the shared DEcircRNAs, DElncRNAs and DEmRNAs in four brain regions: prefrontal cortex, substantia nigra, corpus striatum and hippocampus. Interestingly, few circRNAs/lncRNAs were consistently regulated by Nrf2 in these four brain regions, indicating that Nrf2 may not have the same function in different brain areas. Recent studies have demonstrated that alters the expression of not only antioxidant genes, but also genes involved in autophagy, intermediate metabolism, stem cell quiescence and unfolded protein reactions [9]. Therefore, the function of Nrf2 is complex. Moreover, our research suggested that Nrf2 activity may be tissue-specific, which has not been reported before and is worth studying in the future.
In conclusion, our study demonstrated that four DEcircRNAs (mmu_circ_0000233, mmu_circ_00150 35, cicRNA.20127 and mmu_circ_0012936) and seven DElncRNAs (ENSMUST00000204847, NONMMUT02 4778, NONMMUT132160, NONMMUT132168, NO NMMUT127961, ENSMUST00000194077 and NON MMUT109267) may contribute to the effects of Nrf2 on cognitive function by regulating protein synthesis, cell degeneration, neuronal development, intersynaptic signal transduction and glutathione metabolism. Our study has laid the foundation for future investigations into the molecular pathways through which Nrf2 protects neurons against oxidative stress and exerts other functions.

Nrf2-knockout mice and ethics
Adult male Nrf2 (+/+) and Nrf2 (-/-) mice (25-30 g, 3-4 months, n = 3) on an ICR background were kindly provided by Professor Chun-Yan Li (Department of Neurology, Second Hospital of Hebei Medical University, Shijiazhuang, China) for this study. PCR amplification of genomic DNA from tails was used to determine the genotypes of the mice. All the mice were sacrificed using an overdose of an isoflurane/oxygen mixture (Huazhong Haiwei Gene Technology Co, Ltd, Beijing, China, 021400). Prefrontal cortex tissues were surgically obtained from each mouse, and immediately homogenized for the extraction of total RNA (TRIzol, Invitrogen, Carlsbad, CA, USA, 15596026).
The animal experiments in this study complied with the regulations of the Animal Welfare Act and the National Institutes of Health Guide for the Care and Use of Laboratory Animals (NIH Publication No. 85-23, revised 1996) and were approved by the ethics committee of Hebei Medical University (IACUC-Hebmu-Glp-2016017).

Microarray analysis
Six mice (Nrf2 (+/+) group, n = 3; Nrf2 (-/-) group, n = 3) were anesthetized, and their prefrontal cortex tissues were collected and rapidly stored in liquid nitrogen. Total RNA was extracted and purified using an RNeasy Mini Kit The raw data were normalized with the Quantile algorithm and LIMMA packages in RStudio (http:// www.rstudio.com/). DEcircRNAs, DElncRNAs and DEmRNAs with p-values < 0.05 and FCs > 1.5 were selected through scatter plots and volcano diagrams using the 'pheatmap' package in R language. The data were also visualized through hierarchical clustering analyses.

qRT-PCR validation
Total RNA was isolated from Nrf2 (-/-) and Nrf2 (+/+) prefrontal cortex tissues using Trizol reagent (Thermo Fisher Scientific, Wilmington, DE, USA) according to the instruction manual, and the RNA quality was examined using 1% agarose gel electrophoresis. Then, the total RNA (2 μg) was reverse-transcribed to cDNA using an M-MLV First Strand Kit (Thermo Fisher Scientific Inc., 00341186) with random primers. The levels of six DEcircRNAs and eight lncRNAs were assessed via qRT-PCR using SuperReal PreMix Plus (SYBR Green) (TIANGEN, China, FP205) on an Illumina Eco PCR machine (Illumina, San Diego, CA, USA). The primers are displayed in Table 3. The qRT-PCR reaction conditions were as follows: an initial denaturation step of 15 min at 95° C, followed by 40 cycles of 15 s at 95° C, 20 s at 55° C, and 20 s at 72° C, and a final step of 5 min at 72° C. For normalization, β-actin was selected as an internal control. All experiments were performed in triplicate. The comparative cycle threshold (2 -ΔΔCT ) method was used to calculate the relative FCs.

Construction of ceRNA networks based on DEcircRNA/DElncRNA-miRNA-DEmRNA interactions
DEcircRNAs and DElncRNAs verified by qRT-PCR were selected for the construction of interaction networks by Shanghai Biotechnology Corporation. First, a DEcircRNA/DElncRNA-DEmRNA co-expression analysis was performed based on the Pearson correlation coefficients between DEmRNA and DEcircRNA/ DElncRNA levels in the SBC Mouse (4*180K) ceRNA Array analysis. Genes with Pearson correlation coefficients > 0.81 and p-values < 0.05 were recommended for further analysis. Then, DEcircRNA /DElncRNA-miRNA-DEmRNA interaction networks were constructed using miRanda, and binding sites with relatively high scores (≥ 140) were identified. Finally, the networks were visualized using Cytoscape_v3.7.0 (http://www.cytoscape.org/).

Functional annotation of DEmRNAs that were targets of DEcircRNAs and DElncRNAs
GO and KEGG pathway enrichment analyses were conducted for the target DEmRNAs of miRNAs sponged by the DEcircRNAs/DElncRNAs in the prefrontal cortex tissues. The analyses were performed using the clusterProfiler of R/bioconductor (http://enrich.shbio. com/index/ga.asp) to determine the potential functions of the DEcircRNAs and DElncRNAs.

Coding potential of DEcircRNAs and DElncRNAs
The Coding-Non-Coding Index (https://github.com/ www-bioinfo-org/CNCI) [27], which distinguishes between coding and noncoding transcripts based on their triplet base composition, was used to predict the protein coding potential of the DEcircRNAs and DElncRNAs (The full-length sequences of some of the DEcircRNAs have yet to be confirmed, so the sequences within 300 base pairs around the cut site were used to predict their coding potential).

Expression correlation analysis between DEcirc RNAs/DElncRNAs and their host mRNAs
Pearson correlation coefficients were used to determine the expression correlations between DEcircRNAs/ DElncRNAs and mRNAs encoded by the same parental genes. DEcircRNA/DElncRNA-host mRNA pairs with p-values < 0.05 were retained.

Statistical analysis
All data were analyzed using SPSS 22.0 software (IBM). The means and standard deviations were determined, and independent-samples t-test were performed. P-values < 0.05 was considered significant.