Abstract
Multiple schizophrenia (SCZ) risk loci may be involved in gene co-regulation mechanisms, and analysis of coexpressed gene networks may help to clarify SCZ molecular basis. We have previously identified a dopamine D2 receptor (DRD2) coexpression module enriched for SCZ risk genes and associated with cognitive and neuroimaging phenotypes of SCZ, as well as with response to treatment with antipsychotics. Here we aimed to identify regulatory factors modulating this coexpression module and their relevance to SCZ. We performed motif enrichment analysis to identify transcription factor (TF) binding sites in human promoters of genes coexpressed with DRD2. Then, we measured transcript levels of a group of these genes in primary mouse cortical neurons in basal conditions and upon overexpression and knockdown of predicted TFs. Finally, we analyzed expression levels of these TFs in dorsolateral prefrontal cortex (DLPFC) of SCZ patients. Our in silico analysis revealed enrichment for NURR1 and ERR1 binding sites. In neuronal cultures, the expression of genes either relevant to SCZ risk (Drd2, Gatad2a, Slc28a1, Cnr1) or indexing coexpression in our module (Btg4, Chit1, Osr1, Gpld1) was significantly modified by gain and loss of Nurr1 and Err1. Postmortem DLPFC expression data analysis showed decreased expression levels of NURR1 and ERR1 in patients with SCZ. For NURR1 such decreased expression is associated with treatment with antipsychotics. Our results show that NURR1 and ERR1 modulate the transcription of DRD2 coexpression partners and support the hypothesis that NURR1 is involved in the response to SCZ treatment.
SIGNIFICANCE STATEMENT In the present study, we provide in silico and experimental evidence for a role of the TFs NURR1 and ERR1 in modulating the expression pattern of genes coexpressed with DRD2 in human DLPFC. Notably, genetic variations in these genes is associated with SCZ risk and behavioral and neuroimaging phenotypes of the disease, as well as with response to treatment. Furthermore, this study presents novel findings on a possible interplay between D2 receptor-mediated dopamine signaling involved in treatment with antipsychotics and the transcriptional regulation mechanisms exerted by NURR1. Our results suggest that coexpression and co-regulation mechanisms may help to explain some of the complex biology of genetic associations with SCZ.
Introduction
Schizophrenia (SCZ) is a chronic brain disorder associated with multiple genetic risk factors (Ripke et al., 2014). Several studies showed that risk genes for SCZ are co-regulated and operate in gene networks (Pergola et al., 2017; Chen et al., 2018; Gandal et al., 2018; Li et al., 2018; Meng et al., 2018; Pergola et al., 2019). The investigation of gene networks relevant to SCZ thus allows us to identify molecular dynamics and mechanisms of co-regulation crucially involved in this brain disorder. An approach to characterize genetic networks is based on gene coexpression. Coexpression estimates the simultaneous transcription of genes and is often associated with shared gene functions and co-regulation (Serin et al., 2016; van Dam et al., 2018).
The largest genome wide association study on SCZ includes the DRD2 locus (Ripke et al., 2014), coding for the dopamine D2 receptor. Genetic variation within DRD2 is associated with working memory (WM) deficits and related prefrontal cortex (PFC) activity in SCZ (Zhang et al., 2007; Bertolino et al., 2010; Slifstein et al., 2015). Dopaminergic signaling is altered in SCZ (Abi-Dargham, 2014), and current available antipsychotic medications (APs) either hamper or modulate D2 receptors (Miller, 2009). Using a genome-wide Weighted Genes Coexpression Network Analysis approach, we have previously identified a set of coexpressed genes (module) comprising the DRD2 transcript coding for the long isoform of the D2 receptor (D2L) in postmortem dorsolateral prefrontal cortex (DLPFC) of healthy individuals (Pergola et al., 2017). This module encompassed 85 genes and was significantly enriched for SCZ risk genes (Ripke et al., 2014). Furthermore, we computed a polygenic coexpression index (PCI) related to interindividual variability of gene coexpression, which was associated with behavioral and neuroimaging phenotypes crucially associated with SCZ (i.e., WM performance and related brain activity; Bertolino et al., 2006), as well as with response to treatment with APs (Pergola et al., 2017). Nevertheless, gene networks established by statistical correlation do not provide insight into the regulatory processes underpinning coexpression (Parikshak et al., 2015). Several studies suggest that genes with strongly correlated mRNA expression profiles are more likely to be co-regulated by shared transcription factors (TFs; Allocco et al., 2004; Serin et al., 2016). Interestingly, a proportion of genetic variability associated with SCZ involves genes coding for TFs, i.e., TCF4, CTCF, EGR1, EP300, SREBF1, and STAT6 (Ripke et al., 2014). Although this evidence supports the relevance of TFs in SCZ pathogenesis, little is known as to how the coordinated action of TFs reprograms gene expression. Our previous finding that a coexpression relationship occurs among DRD2 and other genes in DLPFC suggests that correlation of expression profiles takes place at the transcriptional level, implicating that a relevant fraction of the genes associated with DRD2 coexpression are target of the same TF(s) (Pergola et al., 2017). Because dysregulation of such TF(s) in the DLPFC of patients with SCZ may alter the expression patterns of DRD2 and its coexpressed genes (Pergola et al., 2017; Selvaggi et al., 2019), we sought to investigate their role in the pathophysiology of SCZ and its pharmacological treatment.
Here, we performed motif enrichment analysis of TF binding sites located within the human promoter regions of the genes included in the DRD2 coexpression module. Then, we performed validation experiments in mouse primary cortical neurons by measuring mRNA levels of DRD2 coexpression module-selected genes upon gain [overexpression (OE)] and loss [knockdown (KD)] of function of candidate TFs. We found that two TFs, NURR1 (nuclear receptor-related 1; coded by the NR4A2 gene) and ERR1 (estrogen-related receptor 1; coded by the ESRRA gene), are likely to bind promoters of our genes of interest. In mouse cortex, gain and loss of function of Nurr1 and Err1 is associated with expression changes of genes within the DRD2 coexpression module, including Drd2 long isoform. Importantly, we found that NR4A2 and ESRRA expression in human DLPFC is associated with SCZ diagnosis, and the reduced expression observed in patients is associated to treatment with APs for NR4A2.
Materials and Methods
Experimental design and statistical analysis
To carry out the motif enrichment analysis in promoters of the DRD2 coexpression gene set, we used the Pscan software (http://159.149.160.88/pscan/; Zambelli et al., 2009). The software performed a series of z tests, associating a p value with each TF position frequency matrix included in the JASPAR database (Zambelli et al., 2009). Then, we removed from the Pscan list those TFs that were not expressed in human DLPFC, i.e., the brain area where the DRD2 coexpression gene set was detected. To determine whether TFs expression was significantly greater than zero with α = 0.05, we performed one-sample t tests in the same discovery dataset of the DRD2 coexpression gene set (Colantuoni et al., 2011). Finally, we corrected the z test p values of TFs that were expressed in human DLPFC for false discovery rate (FDR; q < 0.05).
Following the motif enrichment analysis, we performed validation experiment in mouse cortical neurons. We used the lentiviral infection to induce OE and KD of Nr4a2 and Esrra. Then, to determine the efficacy of such procedures and the differences in target genes expression upon Nr4a2 and Esrra OE/KD, we performed unpaired t tests across experimental conditions (i.e., Nr4a2 OE, Nr4a2 KD, Essra OE, Essra KD, and respective control conditions), with α = 0.05. Gene expression interquartile range (IQR) was used to detect outliers (Schwertman et al., 2004), defined as observations that fell <Q1 − 1.5 IQR or >Q3 + 1.5 IQR. Each experiment was repeated at least three times (biological replicates) and outliers removed where detected (see Fig. 4-4, and Fig. 5-1). Cortical neurons cultures were prepared weekly from mice embryos ranging from N = 5 to N = 10.
We also aimed to estimate NR4A2 and ESRRA differential expression in the DLPFC of SCZ versus healthy controls (HCs). To this purpose, we used regression analysis. The same regression model was used to assess difference in NR4A2 and ESRRA expression between SCZ who had postmortem positive toxicology testing for APs versus those who did not. Adjustments for genetic stratification (10 genomic ancestries according to Pergola et al., 2019), age, gender, RNA integrity, mitochondrial mapping rate, gene mapping rate and 9 qualitative surrogate variables, were made as suggested by Jaffe et al. (2017) as best practice. Data were analyzed by either GraphPad Prism v7 or R3.5.1.
Motif enrichment analysis
Among the five different region lengths to set for the analysis in the Pscan software web interface, we focused on the largest possible, i.e., 950 bp upstream and 50 bp downstream of the transcription start site (TSS) of genes included in the DRD2 coexpression module, as we aimed to capture also regulatory elements potentially acting as enhancers or repressors proximal to promoters. The JASPAR database was used as a reference for TF position frequency matrices determination.
We replicated the motif enrichment analysis using the Genomatix MatInspector software (Cartharius et al., 2005). To this aim, we used the MatInspector default parameters, i.e., relative scores of matrix similarity >0.80 and of TF binding sites common to at least 85% of the input sequences.
Selection of genes-of-interest
We tested a subset of eight genes from the 85 included in the DRD2 coexpression gene set (Pergola et al., 2017). In particular, we selected the genes according to the following criteria: (1) genes previously associated with SCZ, such as GATAD2A (p = 3.63E−10), SLC28A1 (p = 1.62E−11; (Ripke et al., 2014), and CNR1 (Colizzi et al., 2015; Taurisano et al., 2016); (2) genes harboring single nucleotide polymorphisms previously included in the DRD2-related PCI (Pergola et al., 2017), and showing the highest affinity scores for binding to NURR1 and ERR1 according to Pscan prediction, such as GPLD1, CHIT1, BTG4, OSR1; and (3) genes whose expression was conserved in the mouse brain. In case of multiple annotated transcripts, we designed primers annealing to exon boundaries that are common to all isoforms to detect the total gene level expression, except for the DRD2 gene, whose long and short isoforms (D2L and D2S, respectively), along with total DRD2, were detected. The rationale for this choice was to test whether the TF activity was responsible for the preferential expression of one of the two isoforms.
Primary neuronal cultures
Animal care and experimental procedures were performed in accordance with the University of Trento ethical committee and were approved by the Italian Ministry of Health. Primary cortical neurons were obtained from C57BL/6J mouse embryos at embryonic day (E)15.5 as previously described (Tripathy et al., 2017). Briefly, cortices were isolated from the whole brain, dissociated in Papain solution [20 U papain, 5 mm EDTA, and 30 mm cysteine in 1× Eagles' balanced salt solution (EBSS)] for 20 min at 37°C. After 3 min at 37°C with DNase I, the dissociated tissues were centrifuged at 1000 × g for 5 min. After removing the supernatant, the pellet was resuspended in a solution containing soybean trypsin inhibitor and bovine serum albumin in EBSS to stop the digestion. Cells were pelleted by a centrifugation step of 1000 × g for 10 min. Finally, the cells were plated on Poly-d-lysine-coated plates in neurobasal medium supplemented with B27, PenStrep (100 U), l-glutamine (2 mm), and AraC (100 μm). For cell maintenance, half of the medium was replaced with fresh medium every 3 d. For all the experiments, cells were analyzed at 8 d in vitro (DIV).
Plasmids
Mouse Lenti-ORF Nr4a2 and Esrra viral vectors were purchased from Origene. As a negative control, mCherry sequence was subcloned within the same Lentiviral vector backbone (mCherry vector). Mouse 4unique 29mer shRNA constructs in lentiviral GFP vector for Nr4a2 and Esrra, were used for KD experiments, with the relative scrambled negative control (Scramble; Origene).
Cell lines
Human embryonic kidney (HEK) 293T cells were used to package lentiviral vectors for OE and KD of Nr4a2 and Esrra. Cells were maintained in DMEM High Glucose supplemented with 10% fetal bovine serum, l-glutamine (2 mm) and PenStrep (1%).
Lentiviral vector packaging and titer
PCMV-dR8.91 (Delta 8.9) plasmid containing gag, pol, and rev genes and VSV-G envelope plasmid for lentiviral vector packaging were cotransfected along with Lenti-ORF or Lenti-shRNAs and their control vectors in HEK293T cells using the calcium phosphate method (Kingston et al., 2001). Viral particles were then quantified as described by Pizzato et al. (2009).
Nr4a2 and Esrra OE/KD
One day in vitro cortical neurons were infected with Lenti-ORF Nr4a2 and Esrra viral particles or mCherry Vector for OE experiments, and with shRNA constructs in lentiviral GFP vector or Scramble for KD experiments 4 h after plating, at a multiplicity of infection of 0.125. The next day, culture medium was entirely replaced.
MTT cell viability assay
Cortical neurons were plated on 24 wells at a concentration of 5 × 105 cells/ml. At the day of samples collection, PBS-solved MTT reagent (Sigma-Aldrich) was added to the culture medium to a final concentration of 0.5 mg/ml and incubated for 1 h. At the end of the incubation period, the medium was removed, and the converted dye solubilized with DMSO (Sigma-Aldrich). Cell viability was measured as fold-change in absorbance at 570 nm with background subtraction at 630–690 nm.
Immunoblot analysis
Cortical neurons were plated on 10 cm dishes (6 × 106 cells/dish), washed in ice-cold PBS and lysed in ice-cold RIPA buffer (150 mm NaCl, 6 mm Na2HPO4, 4 mm NaH2PO4, 2 mm EDTA pH8, 1% Na-deoxycholate, 0.5% Triton X-100, 0.1% SDS) in presence of a protease inhibitor cocktail. The lysates were collected, sonicated, and clarified by centrifugation at 4°C at 16,000 × g. Protein extracts were quantified by a BCA assay method (ThermoFisher Scientific). Lysates were denatured at 95°C in 5× sample buffer (0.5 m Tris-HCl, pH 6.8, 10% SDS, 25% glycerol, 0.5% bromophenol blue, 0.5% β-mercaptoethanol), separated on SDS- PAGE gels and transferred to nitrocellulose membranes. Membranes were blocked in 5% dry milk in 0.02 m Tris/HCl pH 7.5, 137 mm NaCl, and 0.1% (v/v) Tween 20 for 1 h at room temperature. Antibodies against Nurr1 (Abcam, ab41917), Err1 (Abcam, ab16363), and β-tubulin (Sigma-Aldrich, T7816) were diluted 1:1000, 1:500, and 1:5000, respectively. Proteins were detected using the Odyssey infrared imaging system (LI-COR Biosciences).
RNA extraction and real-time PCR
Cortical neurons were plated on 6-wells at a concentration of 5 × 105 cells/ml. At 8 DIV cells were lysed in 500 μl of TRIzol Reagent (ThermoFisher Scientific). Then, RNA was extracted following the TRIzol manufacturer's protocol with minor modifications. RNA quantification was performed by reading absorbance at 260 nm, using NanoDrop 2000 (ThermoFisher Scientific). One microgram of total RNA was reverse transcribed into cDNA using the iScript RT Supermix (Bio-Rad Laboratories).
Primers for Real Time PCR were designed using the IDT Primer Quest tool. Oligonucleotides for this experiment were synthesized by Eurofins Scientific (for complete sequences and assay efficiency, see Fig. 4-1). Standard curves for each primer pair revealed efficiency close to 100%. Real-time PCR was performed using a CFX 96 Real Time System (Bio-Rad Laboratories), in a final volume of 20 μl, using the iTaq SYBR Green Supermix (Bio-Rad Laboratories), and with 0.5 mm primers. PCR protocol was performed as follows: 95°C for 30 s, 40 cycles at 95°C for 5 s, and 60°C for 30 s. Melting curve analysis allowed to confirm specificity of the amplified products (Bustin et al., 2009). Mean Ct values of genes of interest were normalized to those of Actin Beta (ACTB) as the reference gene. These normalized values were analyzed through the comparative Ct method for the relative quantification of targets, as previously reported (Pergola et al., 2016).
NR4A2 and ESRRA mRNA expression in human postmortem brain tissue
We investigated NR4A2 and ESRRA mRNA expression in human postmortem DLPFC specimens in two large RNA-seq datasets comparing patients with SCZ to HC, both males and females. These datasets included one from the Lieber Institute for Brain Development (LIBD; “discovery” dataset, polyA+ protocol; N = 155 SCZ cases, N = 196 controls) and one from the CommonMind Consortium (CMC; “replication” dataset, RiboZero protocol; N = 159 patients, N = 172 controls). Detailed sample description has been reported previously (Pergola et al., 2019). Within the SCZ group in the LIBD sample, we tested whether NR4A2 expression levels were different between patients on APs at the moment of death compared with patients who were free from treatment with APs. AP positivity was assessed with toxicological tests performed on frozen postmortem tissue long after death and recorded as a binary yes/no variable.
Results
Promoters of the DRD2 coexpression gene set are enriched for ERR1 and NURR1 binding sites
Using the Pscan software to perform motif enrichment analysis on the upstream regions of all the genes included in the DRD2 coexpression module, we retrieved 263 position frequency matrices from the JASPAR database, each one representing a specific TF sequence binding motif. According to postmortem brain expression data collected in the 249 PFC brain tissues of the BrainCloud dataset, 32% of the matrices were referred to TFs expressed in human PFC as assessed with a one-sample t test (Fig. 1-1). Two sequence motifs were significantly enriched within the promoter regions of genes included in the DRD2 coexpression gene set. Such motifs were recognized by ERR1 (coded by the ESRRA gene) and NURR1 (coded by the NR4A2 gene) TFs (FDR-adjusted p = 0.02 and p = 0.03 for ERR1 and NURR1 matrices, respectively; Fig. 1). All predicted binding sequences identified by Pscan matched ERR1 and NURR1 binding motifs with a fitting >75% and 85%, respectively (Fig. 1-2, and Fig. 1-3). These results were confirmed using another prediction tool, the Genomatix MatInspector software (ERR1: uncorrected p = 1.2 × 10−16; NURR1: uncorrected p = 4.5 × 10−32; Fig. 1-4).
Enrichment analysis of transcription factor binding sites (TFBSs) in the promoter regions of genes belonging to the DRD2-related coexpression network. Putative TFBSs were evaluated for enrichment using the Pscan software with FDR correction for multiple testing. Shown on the x-axis are all transcription factor motifs reported by the JASPAR database and referred to TFs expressed in human PFC. Shown on the y-axis is the p value of the enrichment analysis expressed as −log10 (p value). The red line indicates FDR correction threshold. NR4A2 and ESRRA binding motifs were significantly enriched (p = 0.00076 and p = 0.00023, respectively). TFBS sequence motifs for NR4A2 and for ESRRA are illustrated with the height of each base, indicating the probability of the presence of the corresponding base at the designated position. For further details on TFs expressed in human PFC and predicted TFBSs see Figure 1-1, Figure 1-2, and Figure 1-3. For Genomatix MatInspector software confirmation results see Figure 1-4.
Figure 1-1
Figure 1-2
Figure 1-3
Figure 1-4
According to two independent prediction tools, ERR1 and NURR1 TFs were likely to regulate the expression of genes in the DRD2 coexpression module.
Gain or loss of Nr4a2 and Esrra expression do not affect the viability of primary cortical neurons
To investigate whether NURR1 or ERR1 regulate the expression of genes in the DRD2 coexpression module as predicted, we tested the effect of gain and loss of function of the two TFs in mouse primary cortical neurons. First, we verified the mRNA and protein expression levels of Nurr1 in neurons transduced with lentiviruses that either OE or KD the TF gene. Upon OE, we obtained a fourfold increase of Nr4a2 mRNA (t(6) = 3.06, p = 0.02, unpaired t test; Fig. 2A) and a threefold increase of protein levels compared with control samples (t(4) = 3.165, p = 0.03, unpaired t test; Fig. 2B; Fig. 2-1). With KD, Nr4a2 mRNA transcript levels were decreased by 75% (t(8) = 4.232, p = 0.003, unpaired t test; Fig. 2C) and protein levels by 50% compared with the control condition (t(4) = 15.02, p = 0.0001, unpaired t test; Fig. 2D; Fig. 2-1). Similarly, the transcript levels of Esrra mRNA were enhanced by 20-fold (t(6) = 2.785, p = 0.03, unpaired t test; Fig. 3A) and protein levels by 1.5-fold compared with the control condition (t(4) = 9.192, p = 0.0008, unpaired t test; Fig. 3B; Fig. 3-1). On the other hand, KD of Esrra resulted in 60% decrease in mRNA transcript levels (t(4) = 4.75, p = 0.009, unpaired t test; Fig. 3C) and 50% decrease in protein levels (t(4) = 3.123, p = 0.03, unpaired t test; Fig. 3D; Fig. 3-1). To exclude any possible toxic effect due to lentiviral infection and/or aberrant expression of the two TFs, we used a cell viability assay. Nr4a2 and Esrra OE and KD mediated by lentiviral transduction did not induce cellular death on neuronal cultures in any of our experimental conditions compared with both cells transduced with lentiviruses expressing scramble RNA or mCherry as well as non-infected neurons (WT; Fig. 2E,F; Fig. 3E,F). Together, these results show that modulation of Nr4a2 and Esrra expression in primary neurons is not associated with overt toxic effects.
Nr4a2 OE and KD efficiency and toxicity in mouse cortical neurons. Real-time PCR (A, C) and Western blotting (B, D) analyses in primary mouse cortical neurons were performed to verify the efficiency of Nr4a2 OE and Nr4a2 KD with respect to the control conditions (mCherry vector and scramble, respectively). See Figure 2-1 for full length western blots. Neuronal viability was assessed by MTT assay (E, F). All results are the mean ± SEM of three independent experiments. *p < 0.05, **p ≤ 0.01, ***p ≤ 0.001.
Figure 2-1
Esrra OE and KD efficiency and toxicity in mouse cortical neurons. Real-time PCR (A, C) and Western blotting analysis (B, D) were performed to verify the efficiency of Esrra OE and Esrra KD with respect to the control conditions (mCherry vector and scramble, respectively). See Figure 3-1 for full length western blots. Neuronal viability was also determined by MTT assay in all the experimental conditions (E, F). All results are the mean ± SEM of three independent experiments. *p < 0.05, **p ≤ 0.01, ***p ≤ 0.001.
Figure 3-1
Nr4a2 and Esrra modulation alters the expression of target genes within the Drd2 coexpression module
To assess whether Nr4a2 and Esrra OE and KD modulate expression of our genes of interest, we measured their transcript levels in basal conditions and upon OE and KD of the two TFs. We first evaluated the effect of Nr4a2 and Esrra OE in primary neurons (Fig. 4; for graph with all data points and outliers plotted, see Fig. 4-4). Increased levels of Nr4a2 were associated with a twofold upregulation of Gpld1 (t(8) = 2.302, p = 0.04, unpaired t test), as well as with a 50% downregulation of Cnr1 (t(6) = 7.787, p = 0.0002, unpaired t test) and Chit1 (t(6) = 5.49, p = 0.002, unpaired t test) expression levels compared with the control condition (Fig. 4A). Esrra OE was associated with downregulation of D2L (50%; t(6) = 8.251, p = 0.0002, unpaired t test), Slc28a1 (60%; t(4) = 3.839, p = 0.02, unpaired t test), and Osr1 (80%; t(4) = 7.505, p = 0.002, unpaired t test) compared with the control condition, and no change in the expression of Cnr1, Gatad2a, Gpld1, Btg4, and Chit1 was observed under these experimental conditions (Fig. 4B). These results indicate that gain of Nurr1 function suppressed the expression of Cnr1 and Chit1 gene and enhanced that of Gpld1. On the other hand, gain of Err1 function suppressed Slc28a1 and Osr1 gene expression, along with D2L. This analysis suggests that specific genes within the DRD2 coexpression module associated with the risk for SCZ are negatively regulated in conditions in which the function of two key TFs is aberrantly enhanced.
Nurr1 and Err1 OE induces expression changes in target genes of the DRD2-related coexpression module. Real-time PCR analysis of the indicated genes in primary mouse cortical neurons upon Nr4a2 OE (A) and Esrra OE (B) normalized to the control condition (dashed line). All results are expressed as the mean ± SEM of four independent experiments. For complete primer sequences and assay efficiency see Figure 4-1. For further analyses verifying consistency of TFs OE/KD results with the coexpression analysis of our previous work from Pergola et al., 2017 see Figure 4-2, and Figure 4-3. For graph with all data points and outliers plotted see Figure 4-4. *p < 0.05, **p ≤ 0.01, ***p ≤ 0.001.
Figure 4-1
Figure 4-2
Figure 4-3
Figure 4-4
We next assessed the effect of the two endogenous TFs by gene KD (Fig. 5A; for graph with all data points and outliers plotted, see Fig. 5-1). Nr4a2 KD was associated with downregulation of Gpld1 (50%; t(6) = 3.265, p = 0.02, unpaired t test), D2L (60%; t(6) = 3.859, p = 0.008, unpaired t test), Drd2 total gene levels (50%; t(6) = 6.886; p = 0.0004), and Gatad2a (75%; t(8) = 9.385, p = 1.36076E−05, unpaired t test), and with a twofold upregulation of Cnr1 (t(6) = 3.044, p = 0.02, unpaired t test), and a eightfold upregulation of Slc28a1 (t(4) = 10.3, p = 0.0005, unpaired t test). D2S transcript levels were not affected by Nr4a2 KD (Fig. 5A). On the other hand, Esrra KD was associated with a 1.6-fold upregulation of Drd2 total gene levels (t(4) = 37.61, p = 3E−06), and downregulation of D2S (55%; t(6) = 2.745, p = 0.03, unpaired t test), Cnr1 (55%; t(6) = 2.851, p = 0.03, unpaired t test), and Slc28a1 (50%; t(6) = 4.104, p = 0.006, unpaired t test) compared with the control condition (Fig. 5B). Altogether, six and five of the eight examined genes displayed significant expression profile changes upon modulation of expression of Nr4a2 and Esrra, respectively.
Nurr1 and Err1 KD modulates the expression of target genes of the DRD2-related co-expression module. Real-time PCR analysis of the indicated genes in primary mouse cortical neurons upon Nr4a2 KD (A) and Esrra KD (B) normalized to the control condition (dashed line). All results are expressed as the mean ± SEM of four independent experiments. For graph with all data points and outliers plotted see Figure 5-1. *p < 0.05, **p ≤ 0.01, ***p ≤ 0.001, ****p ≤ 0.0001.
Figure 5-1
We also sought to verify whether the results described above were consistent with the coexpression analysis of our previous work (Pergola et al., 2017; Fig. 4-2, and Fig. 4-3).
These results support the idea that these two TFs control the expression of the DRD2 coexpression module associated with the risk for SCZ.
NR4A2 and ESRRA expression is reduced in patients with SCZ
Next, we set out to determine the potential role of NURR1 and ERR1 for the pathophysiology of SCZ and its contribution to pharmacological treatment. Therefore, we compared NR4A2 and ESRRA expression in the DLPFC of SCZ patients and HC using RNA-seq data from postmortem DLPFC specimens of the LIBD (discovery set) and CMC (replication set) samples. We found that NR4A2 and ESRRA expression levels were significantly reduced in patients in the discovery sample (NR4A2: t = −2.8, p = 0.0048; ESRRA: t = −2.5, p = 0.013; multiple regressions). In the replication sample only NR4A2 was significantly reduced (t = −1.9, one-tailed p = 0.032, multiple regression), but not ESSRA (t = −1.0, one-tailed p = 0.32; multiple regression). To investigate whether the replicated difference between patients and control subjects was associated with treatment with APs, we analyzed the LIBD sample by comparing patients who were being treated with APs (N = 90) versus those who were not (N = 50) using a regression model. We found that patients treated with APs have reduced NR4A2 expression (t = −2.46, p = 0.015, multiple regression) compared with those who were not. Moreover, NR4A2 expression levels in controls did not differ from those detected in SCZ patients that were negative at toxicological tests for APs treatment (t = −0.9, p = 0.4, multiple regression).
Discussion
In a previous study (Pergola et al., 2017), we investigated how prefrontal gene coexpression of DRD2 is associated with clinical and system-level phenotypes, i.e., WM performance and related brain activity. Using a network approach, we identified coexpression trans-eQTLs to predict gene coexpression, and identified a coexpression gene set in which other genes associated with SCZ risk converged, such as CNR1 (Colizzi et al., 2015), GATAD2A and SLC28A1 (Ripke et al., 2014). We found that individuals with increased simultaneous expression of the DRD2 gene set manifested greater PFC activity as well as less efficient WM processing, along with greater clinical improvement in response to D2 antagonist antipsychotic drugs.
In the current study, we sought to extend our understanding of the biological mechanisms underpinning coexpression of these genes and their relevance to SCZ and its pharmacological treatment. We tested the hypothesis that the identified DRD2 coexpression gene set is associated with transcription regulation processes mediated by one or more shared TFs. According to two independent bioinformatics tools, we found that NURR1 and ERR1 are potential transcriptional regulators of the DRD2 coexpression gene set. To confirm these findings and provide further evidence of the role of NURR1 and ERR1 in the regulation of DRD2 coexpression module, we performed biological validation experiments in mouse primary neuronal cultures. We tested the effect of gain and loss of function of NURR1 and ERR1 on the expression of selected genes belonging to the identified DRD2 coexpression module. The results were consistent with the bioinformatic prediction and support a model whereby NURR1 and ERR1 are key regulators of the DRD2-related network genes.
More in detail, Nr4a2 OE resulted in reduced Cnr1 and increased Gpld1 expression, and Nr4a2 KD had the opposite effect on the expression of these two genes. Such a strong relationship between NURR1 and Cnr1 and Gpld1 levels suggests that these two genes are directly regulated by NURR1. CNR1 codes for the cannabinoid receptor 1 (CB1), which is involved in brain functions such as emotional responses, motivated behavior (Mouslech and Valla, 2009), cognitive processing and motor control (Ruiz-Contreras et al., 2013, 2014). Endocannabinoids are CB1 partial agonists. Cannabis use and abuse contribute to increase the risk for SCZ (Murray et al., 2007), and genetic variation within the CB1 receptor has been associated with cognitive and imaging endophenotypes of SCZ (Colizzi et al., 2015; Taurisano et al., 2016). Conversely, GPLD1 involvement in SCZ pathophysiology has not been explored so far. GPLD1 encodes an enzyme with phosphatidylinositol-glycan-specific phospholipase D activity, which degrades the glycosylphosphatidylinositol anchor that tethers proteins to the extracellular face of plasma membranes, thereby releasing the target protein from the plasma membrane. A recent study has listed GPLD1 among the genes differentially expressed in blood samples of patients with psychotic disorder (English et al., 2018). Nr4a2 KD was associated with altered expression of Gatad2a, Slc28a1, and Drd2 (total gene level and long isoform). Thus, expression of six of eight genes was modulated by Nr4a2 levels, supporting the view that NURR1 regulates transcriptional processes of genes in the DRD2 coexpression network acting as enhancer or repressor. Concerning Gatad2a, Slc28a1, and D2L expression, we found association with Nr4a2 KD, but not with OE. Moreover, three of eight genes were downregulated upon Esrra OE, and three upon Esrra KD, with no evident correlation between the effects of OE and KD. We obtained 4-fold and 20-fold increase in Nr4a2 and Esrra expression levels, respectively. It is possible that the levels of expression of endogenous genes are sufficient to regulate the transcription of the target genes, resulting in no further change in target gene expression upon OE. Rather, suppression of expression of the endogenous TFs resulted in major changes in the transcriptional regulation of expression of the target genes. ERR1 as well as NURR1 are known to act in molecular complexes with transcriptional coactivators or corepressors (Sacchetti et al., 2001; Huss et al., 2015). The involvement of coactivators and corepressors may induce the same effect on target genes because the TF may act in couple with a coactivator and in turn with a corepressor.
Concerning DRD2 at total gene level expression, we found downregulation and upregulation in the Nr4a2 and Esrra KD conditions, respectively. These findings are in line with those obtained for D2L upon Nr4a2 KD, whereas D2S changed neither upon KD nor upon OE. The picture is not quite as clear for Esrra, because we observed downregulation of D2L upon OE and of D2S upon KD. Our findings are consistent with the fact that the coexpression gene set included selectively the long isoform D2L and could be explained by mechanisms of TF-mediated pre-mRNA splicing. Such mechanisms are usually exerted by controlling transcription elongation rates of RNA polymerase II (Brzyżek and Świeżewski, 2015), or the recruiting of splicing factors to pre-mRNA (Auboeuf et al., 2005; Pandit et al., 2008). Consistently, the STRING database shows that NURR1 and ERR1 interact with mediator 1, a coactivator that functions as a bridge to convey information from regulatory proteins to the RNA polymerase II transcription machinery (Szklarczyk et al., 2015). In addition, NURR1 has been found to interact with splicing factor proline- and glutamine-rich (Jacobs et al., 2009). Thus, our results are consistent with a possible role of NURR1 and ERR1 in mediating pre-mRNA splicing of DRD2 gene by interacting with the transcriptional and the RNA processing machineries.
NURR1 and ERR1 are nuclear receptors with a well characterized crystal structure, even though their natural ligands are unknown (Giguère et al., 1988; Wang et al., 2003; Huss et al., 2015). NR4A2 is expressed in dopaminergic neurons of substantia nigra, ventral tegmental area, PFC, hippocampus, and other brain regions (Xiao et al., 1996; Zetterström et al., 1996a,b), where it is involved in development, differentiation, maintenance and survival processes. This TF regulates the transcription of genes crucially contributing to dopamine signaling, such as DRD2, tyrosine hydroxylase, vesicular monoamine transporter, tyrosine kinase receptor, and dopamine transporter (Wallén et al., 2001). NR4A2-deficient mice have symptoms associated with dysfunction in dopamine neurotransmission, such as hyperactivity in novel environments, abnormal retention of emotional memory and increased locomotor response to stress (Rojas et al., 2007; Vuillermot et al., 2011). Indeed, a recent study confirmed that NR4A2 has a pivotal role in neurodevelopment (Li et al., 2018) and association between NR4A2 genetic variation and diagnosis of SCZ has been reported (Liu et al., 2015). Altogether, these results are consistent with NURR1 implication in SCZ.
ESRRA is highly expressed in tissues with high-energy demands, such as brain, heart, and skeletal muscle, where it acts as energy sensor to control cellular adaptation to energy demand and stress (Heard et al., 2000). It has been implicated in the regulation of biological processes of critical relevance to neuronal function (Ranhotra, 2010), and risk for psychiatric disorders (Cui et al., 2013). Interestingly, a recent study showed that ESRRA is a target of the SCZ-associated microRNA-137 (de Sena Cortabitarte et al., 2018).
In the current study, we also analyzed NR4A2 and ESRRA expression in postmortem PFC specimens from SCZ patients and HCs. We found significant reduction of NR4A2 expression in the PFC of patients with SCZ compared with HC, whereas there is no clear relationship between ESRRA expression and SCZ. Differences between the two groups in NR4A2 expression are associated with treatment with APs in patients. Considering that toxicology was performed at the moment of death, we cannot assume that patients who were “AP-negative” had never been treated with APs for their whole life. This is a limitation of the toxicological testing approach when clinical data are lacking. However, these results are suggestive and worth further exploration, because some studies have pointed to a crosstalk occurring between NURR1 and D2 receptor, showing that D2 receptor functions as a potent NURR1 activator (Tseng et al., 2000; Kim et al., 2006). If confirmed, this evidence would suggest that APs, by blocking the D2 receptor, modulate the expression of NURR1, which in turn regulates the coexpression of the whole gene set of DRD2 partners. Therefore, our results may suggest another molecular mechanism by which APs exert their clinical effect.
A limitation of the current study is the lack of a clear evidence demonstrating NURR1 and ERR1 physical binding to regulatory sequences of the genes investigated in the present study. Moreover, we focused only on cis-acting regulatory elements, hence potentially missing trans-acting regulators, including coactivators and corepressors other than NURR1 and ERR1. This limitation calls for further studies that might help the understanding of these mechanisms of regulation.
Footnotes
This work was supported by the Italian Ministry of Health (Grant PE-2011-02347951), by the Lieber Institute for Brain Development, by a Hoffmann-La Roche Collaboration Grant awarded to G.P., Association Française contre les Myopathies (22221 to M.P. and M.B.), by the Telethon-Italy and Provincia Autonoma di Trento (TCP12013 to M.P.), and by the European Union Seventh Framework Programme for research, technological development, and demonstration under Grant 602450 (IMAGEMEND). G.P.'s position is funded by the European Union's Horizon 2020 research and innovation program under the Marie Sklodowska-Curie Grant 798181 (FLOURISH). This paper reflects only the author's views and the European Union and Research Executive Agency are not liable for any use that may be made of the information contained therein. We thank Debasmita Tripathy and Marta Tarter for assistance with brain tissue processing, and Maria J. Polanco for assistance with real time PCR data analysis.
A.B. has received consulting fees from Biogen and lecture fees from Otsuka, Janssen, Lundbeck, and A.R. has received travel fees from Lundbeck. The remaining authors declare no competing financial interests.
- Correspondence should be addressed to Alessandro Bertolino at alessandro.bertolino{at}uniba.it
