Genome-Wide Analysis Identifies NURR1-Controlled Network of New Synapse Formation and Cell Cycle Arrest in Human Neural Stem Cells.

Nuclear receptor-related 1 (Nurr1) protein has been identified as an obligatory transcription factor in midbrain dopaminergic neurogenesis, but the global set of human NURR1 target genes remains unexplored. Here, we identified direct gene targets of NURR1 by analyzing genome-wide differential expression of NURR1 together with NURR1 consensus sites in three human neural stem cell (hNSC) lines. Microarray data were validated by quantitative PCR in hNSCs and mouse embryonic brains and through comparison to published human data, including genome-wide association study hits and the BioGPS gene expression atlas. Our analysis identified ~40 NURR1 direct target genes, many of them involved in essential protein modules such as synapse formation, neuronal cell migration during brain development, and cell cycle progression and DNA replication. Specifically, expression of genes related to synapse formation and neuronal cell migration correlated tightly with NURR1 expression, whereas cell cycle progression correlated negatively with it, precisely recapitulating midbrain dopaminergic development. Overall, this systematic examination of NURR1-controlled regulatory networks provides important insights into this protein's biological functions in dopamine-based neurogenesis.

INTRODUCTION

The cell bodies of midbrain dopaminergic (mDA) neurons are localized in the substantia nigra and ventral tegmental area, where these neurons play critical roles in the central regulation of motor and motivational behaviors (Elsworth and Roth, 1997; Wallen and Perlmann, 2003). Progressive degeneration of dopaminergic (DA) cells leads to the cardinal symptoms of Parkinson’s disease (Shastry, 2001; Wallen and Perlmann, 2003). Similar to the development of other neuron types in the central nervous system (CNS), maturation of mDA neurons may require neuron migration, cell positioning, and axon wiring, all processes for establishing neural networks in different neuronal systems.

A previous report described a series of events that are important for mesencephalic DA cell differentiation in mice (Andersson et al., 2013; Blakely et al., 2011; Deng et al., 2011; Di Salvio et al., 2010; Inestrosa and Arenas, 2010; Panman et al., 2014; Prakash and Wurst, 2006; Smits et al., 2006; Theofilopoulos et al., 2013; Van den Heuvel and Pasterkamp, 2008; Wallen et al., 1999; Yang et al., 2013; Zhang et al., 2013). Initially, the mDA progenitor appears in the ventricular-most cell layers within the ventral mesencephalon at embryonic day (E)10.5, undergoes neurogenesis, and gives rise to mDA neuroblasts. The mDA precursor cells then migrate from the ventricular zone to a medial location in the mantle zone and exit the cell cycle (become post-mitotic) between E11.5 and E13.5. At this stage, expression of nuclear receptor-related 1 (Nurr1) protein, a transcription factor, is strongly upregulated.

This transcription factor is a member of a family of nuclear receptors that are critical for the development and survival of DA neurons. During the development of DA neurons, Nurr1 expression is induced in early postmitotic DA progenitors and then maintained during differentiation and in the adult brain (Perlmann and Wallen-Mackenzie, 2004). Of interest, in mice lacking the Nurr1 gene, DA neurons fail to differentiate. In addition, in contrast to wild-type animals, in the knockout mice, DA progenitor cells show no lateral migration in the midbrain, fail to innervate the striatal target area, and become apoptotic (Castillo et al., 1998; Saucedo-Cardenas et al., 1998; Zetterstrom et al., 1997). In adult DA neurons, Nurr1 also maintains fiber integrity, and Nurr1 ablation results in a progressive pathology associated with reduced striatal DA, impaired motor behaviors, and dystrophic axons and dendrites (Kadkhodaei et al., 2013). The early development of mDA neurons has received much attention in recent years, particularly with regard to morphogenesis, progenitor specification, mDA differentiation, and neuritogenesis (Andersson et al., 2013; Blakely et al., 2011; Deng et al., 2011; Di Salvio et al., 2010; Inestrosa and Arenas, 2010; Prakash and Wurst, 2006; Smits et al., 2006; Theofilopoulos et al., 2013; Van den Heuvel and Pasterkamp, 2008). Little is known, however, about the downstream targets of Nurr1 that are involved in this pathway.

Few reports are available that characterize the gene expression profile of Nurr1 overexpression (Jacobs et al., 2009a; Sousa et al., 2007), and no reports focused on human development are available. The dramatic differences between mice and humans in brain development—including size, events at specific stages, and structural proportions—likely reflect meaningful differences in developmental gene expression (Bohland et al., 2010; Hawrylycz et al., 2012; La Manno et al., 2016; Lein et al., 2007; Myers et al., 2015; Ng et al., 2009). In the absence of human developmental data, gaps in understanding of gene expression profiling are likely, given the limitations of mouse-derived data. Furthermore, the two available reports describe Nurr1 as regulating many genes in the murine mesencephalon MN9D cell line (Jacobs et al., 2009b) and in meso-diencephalic dopamine neurons (Sousa et al., 2007). Neither group, however, confirmed whether the affected genes are direct targets of Nurr1 regulation. For this reason, using human neural stem cells (hNSCs), we sought to address this gap and identify direct targets of NURR1. Here, we show that NURR1 modulates sets of genes implicated in cell migration, synapse wiring, and postmitotic events.

MATERIALS AND METHODS

Cell cultures

Immortalized hNSC lines (HB1.F3, HB1.F5, and HB1.A4) (Kim, 2004) and the human embryonic kidney cell line HEK293 were maintained and passaged on uncoated culture dishes in Dulbecco’s modified Eagle medium (Gibco – ThermoFisher Scientific, USA) with 10% fetal bovine serum (Hyclone), and 10 µg/ml penicillin-streptomycin (Gibco, USA). All cells transduced with the plasmids pLPCX or pLPC-NURR1 were maintained in the same culture condition. Cells from the mouse amphotropic retrovirus packaging cell line PA317 were cultured in RPMI 1640 medium with 10% fetal bovine serum. All cells were incubated at 37°C with 5% CO2.

Retrovirus-mediated gene transfer

We used an amphotropic replication-incompetent retroviral vector to infect target cells, as previously described (Kim et al., 2013). A vector encoding the human NURR1 gene was generated using pLPCX to infect the PA317 retrovirus packaging cells. Cells then were infected with a retrovirus encoding human NURR1 (pLPC-NURR1) or with the control vector (pLPCX).

Reverse transcription polymerase chain reaction (RT-PCR) and Northern blot analysis

Total RNA was extracted from all cells using the RNeasy Mini Kit (Qiagen, Germany). Total RNA (2 μg) was subjected to reverse transcription with Superscript II reverse transcriptase (Invitrogen, USA) according to the manufacturer’s instructions. Amplification reactions were performed with 1/10 volume of the reverse-transcribed product in a final volume of 25 ml using recombinant Taq DNA polymerase (Invitrogen), as described previously (Kim et al., 2013).

For Northern blotting, 20 mg of total RNA was separated on a 1% formaldehyde-agarose gel and transferred to polyvinylidene fluoride membranes (Millipore, USA). Blots were hybridized with a 1.8-kb BamHI–MluI fragment of the NURR1 gene labeled with [α-32P]-dCTP (3000 Ci/mmol; Amersham, USA). X-ray film was exposed for 2 days at –80°C and developed.

DNA microarray analysis

The GeneChip Human Genome U133 Plus 2.0 Array (Affymetrix, USA) was used to analyze differential gene expression profiles modulated by NURR1 protein expression. Total RNA was extracted according to the manufacturer’s protocol (Affymetrix). For extraction of RNA from hNSCs overexpressing NURR1, we harvested the cells within one week of puromycin selection after transduction to avoid shutdown of NURR1 expression. Hybridization to the U133A DNA microarray, washing, and scanning were performed according to the manufacturer’s protocol, and expression patterns were compared between samples.

To normalize signals, we divided each probe by the average value of the chip to avoid differences between different chips and experiments. NCBI UniGene cluster nomenclature was used to describe uncharacterized sequences. All of the raw microarray data are available from Gene Expression Omnibus (GEO) (series record No. GSE58475). For hierarchical clustering of the list of differentially expressed genes using Hamming distance correlation, we used software made publicly available by the lab of Michael Eisen at the University of California, Berkeley.

Quantitative real-time PCR

Quantitative real-time RT-PCR assays were carried out using gene-specific double-labeled fluorescent probes and sets of specific primers in an ABI PRISM 7700 Sequence detection system (PE Applied Biosystems, USA). The primers and probes were obtained for 12 selected genes and the glyceraldehyde phosphate dehydrogenase (GAPDH) gene control from Assays-on-Demand system, which are ready-made and validated assays for the human genome (Applied Biosystems). To compensate for differences in cell number and/or RNA recovery, the copy number of hNSC mRNA was determined relative to human GAPDH mRNA simultaneously using the comparative ΔCt method.

RNA interference knockdown

Four sets of human NURR1 small interfering (si) RNA duplexes (19 nucleotides) were purchased from Dharmacon RNA Technologies (USA), and used as a cocktail. As a negative control, the commercially available control siRNA (Non-specific Control Duplex IX; Dharmacon RNA Technologies, USA) was used. Transfection of the siRNAs was performed according to the manufacturer’s instructions. Briefly, cells were transfected with 100 nM NURR1 siRNA by Lipofectamine Plus (Invitrogen). Forty-eight hours after transfection, cells were harvested for real-time PCR analysis.

Cloning of the 5'-flanking region of human NMU and G0S2 genes and mutagenesis

We cloned the 5'-flanking region of the human NMU (neuromedin U) and G0S2 (lymphocyte G0/G1 switch) genes using a PCR-based method of genomic DNA from SH-SY5Y cells. The following sets of primers were designed: forward, NMUp (5'--3') and reverse NMUp (5'-gggggtaccACATTCCAAGCAGCCTGGTTCA-3'); and forward, G0S2p (5'-gggaagcttCTCGGCGGCTGCGGGAGGCT-3') and reverse, G0S2p (5'gggggtacccatttctgccagttatcagaggt-3'). To facilitate subcloning, the enzyme site (indicated with lowercase letters) was added to the 5'-end of each primer. The PCR conditions were 25 cycles at 94°C for 1 min, at 60°C for 45 s, and at 72°C for 3 min. PCR products were directly cloned into a pGEM-T Easy vector (Promega, USA) and confirmed by sequence analysis. To generate luciferase reporter constructs, each promoter fragment (human NMU, 2102 bp; human G0S2, 2082 bp) was subcloned upstream of the firefly luciferase reporter gene of pGL3-Basic vector (Promega) employing KpnI/HindIII and KpnI/XhoI sites.

The QuickChange Mutagenesis kit (Stratagene, USA) was used to perform site-directed mutagenesis of the NurRE element within the human NMU promoter. The human NMU-luc construct was used as a template. Oligonucleotides with 34 to 36 nucleotides containing the desired point mutations were created as follows, with the mutated residues underlined: mNurRE, 5'-GTTCCCTGTTTTTCAAAAACAGGTCAAATA-3. The site-directed mutations were confirmed by sequence analysis.

Electrophoretic mobility shift assay (EMSA)

Nuclear extracts from HB1.F3 cells were prepared as described earlier (Kim et al., 2013). Sense and antisense oligonucleotides were annealed and then end-labeled with [γ-32P]ATP (Amersham) and T4 polynucleotide kinase. Labeled probes were purified on 19% non-denaturing polyacrylamide gels. The DNA–protein binding reaction was performed in a final volume of 20 μl reaction buffer containing 10 mM Tris (pH 7.6), 50 mM NaCl, 0.5 mM dithiothreitol, 0.5 mM EDTA, 1 mM MgCl2, 5% glycerol, and 250 μg of poly(dI–dC) per milliliter. Nuclear extract (20 μg of protein) was added to the reaction buffer in the absence or presence of unlabeled competitor DNA and pre-incubated for 15 min on ice. Radioisotope-labeled probes (50,000 cpm) were added, and the mixture was incubated for a further 30 min at room temperature. To resolve DNA–protein complexes, electrophoresis was performed on a 5% non-denaturing polyacrylamide gel. Gels were fixed, dried, and visualized by autoradiography. The oligonucleotides employed were as follows (only sense strands presented):

NMU-NurRE, 5'-GTTCCTCACCTTTCAAAGGGAGGTCAAATA-3';

NMU-mtNurRE, 5'-GTTCCCTGTTTTTCAAAAACAGGTCAAATA-3'; and

G0S2-NBRE2, 5'-CATCACTGACCTTTGCAATT-3'.

Chromatin immunoprecipitation (ChIP)

Soluble chromatin from HB.F3 cells was prepared as described previously (Kim et al., 2013) and immunoprecipitated with antibody against NURR1(TransCruz, sc-991X). Specific primer pairs were designed to amplify the promoter region of human NMU (forward, 5'-GGCTTAATTCCTAGCTGACGGGT-3'; and reverse, 5'-TCTAACTGGTCTGGAGCGCTGGT-3') and human G0S2 (forward, 5'-AGTGGGACCTTCGCGTGCAC-3'; and reverse, 5'-CCAGCTAATCTTGGGGAGGGCTC-3') from human genomic DNA. PCR conditions were optimized to allow for semiquantitative measurement, and PCR products were visualized on 1% agarose-Tris-acetate EDTA gels

Transfection and luciferase reporter assays

Transfections were performed using Lipofectamine PLUS reagent (Invitrogen) according to the manufacturer’s instructions. A total of 3 mg of DNA was used per transfection, as previously described (Kim et al., 2013), and briefly as follows: 1 to 2 mg of pLPCX or pLPCX-NURR1, 1 mg of reporter vector, and 0.5 μg of the internal control plasmid, pSV-β-galactosidase (Promega). All transfections were carried out three times in triplicate. Promoter activities were determined using the Single-Luciferase Reporter Assay System (Promega). Luciferase activities were normalized based on β-galactosidase activity in each well.

BioGPS gene expression analyses

We used the publicly available gene expression atlas at the BioGPS website (http://biogps.org/) (Wu et al., 2016) to determine whether the identified genes showed tissue-specific expression. The embedded GeneAtlas U133A gcrma gene expression activity chart was used to interrogate gene expression data from 176 human cell types and tissues. In accordance with previously established stringency thresholds (Pennings et al., 2011), we considered genes to be tissue-specific if (1) they corresponded to a single organ or tissue with an expression value > 30 multiples of the median and (2) there was no unrelated tissue with an expression value greater than one-third of the maximum expression value.

Disease-associated single nucleotide polymorphisms (SNPs) from genome-wide association study (GWAS) data

SNPs associated with disease susceptibility at the loci of genes that are direct targets of NURR1 were searched in July 2016 from the previously identified variants of GWAS datasets from different chip platforms. Because the GWAS database is a mixture of disease- and trait-associated SNPs from many different sources, we selected SNPs linked to conditions such as schizophrenia, bipolar disorder, depression, and Parkinson’s disease. The URL for data presented here is GWASdb (GWASdb v2; http://jjwanglab.org/gwasdb) (Li et al., 2016b).

Statistics

Data analysis was performed using GraphPad Prism5 (GraphPad Software, USA). Comparisons between 2 groups were analyzed using a two-tailed unpaired Student’s t-test. Comparisons between multiple groups were analyzed using one-way ANOVA with Benjamini & Hochberg’s method. Statistical significance was defined as a P value of less than 0.05. Levels of significance were indicated as *P < 0.05, **P < 0.01, and ***P < 0.001.

RESULTS

Profiling NURR1-regulated gene expression by microarray in hNSCs

Loss of Nurr1 resulted in ventral mDA agenesis (Zetterstrom et al., 1997). Despite the many approaches that have been used to identify Nurr1 target genes in rodents (Jacobs et al., 2009a; Sousa et al., 2007), the full range of Nurr1 target genes is unlikely to have been identified. In particular, NURR1 targets have not been systematically studied in DA neurogenesis through human embryo stages, despite the requirement for Nurr1 during this time for specification of ventral mDA neurons. Thus, to identify additional NURR1 targets in human systems, we performed microarray expression profiling of three hNSC lines overexpressing NURR1 (Fig. 1A).

Fig. 1

Genome-wide screen for direct NURR1 targets.

(A) Schematics for identification of direct NURR1 target genes. (B) Semi-quantitative RT-PCR analysis of the early CNS, mesencephalic, neuronal, astrocyte, and DA markers in the three hNSC lines used in this study. (C) Comparison of NURR1-modulated transcript expression patterns among cell lines using the Eisen lab hierarchical cluster analysis software. The cluster contained 180 genes that were up- or downregulated in the three cell lines and contained NURR1-binding site(s) in their promoter. The dendrogram on the left shows different clusters of genes segregated according to the pattern of regulation in the three cell lines. (D) Functional classification of direct NURR1 target genes, obtained using software from the Eisen lab.

These cell lines—HB1.F3, HB1.A4, and HB1.F5—represent separately isolated clones from primary human mesencephalon tissues immortalized with the v-myc gene (Kim, 2004; Kim et al., 2002) and express several markers, including nestin (Fig. 1B). In addition, microarray data showed that these cell lines highly express genes such as ALDH2, TNC, Nestin, and FOXA2 that are associated with DA progenitors, involved in DA differentiation, or characterize mature phenotypes (Supplementary Table S1). Based on these data, we determined to use the human NSC lines to investigate NURR1-related gene expression that is relevant for DA development. To trigger NURR1 overexpression in the hNSC lines, we transduced them with pseudovirus harboring plasmids overexpressing NURR1 cDNA (pLPCX or pLPC-NURR1). To avoid detecting secondary gene expression changes resulting from downstream effects of a direct NURR1 target, we treated cells with puromycin for only one week to select transduced cells and checked the surviving cells for NURR1 expression (Supplementary Fig. S1A). We used Affymetrix U133 Plus 2.0 oligonucleotide microarrays to analyze mRNA gene expression in RNA extracted from these cells (Supplementary Fig. S1).

Software made available by the lab of Michael Eisen at the University of California, Berkeley, was used to filter the microarray data for each cell line. Of the > 40,000 gene queries, we found 2,437 that were upregulated by NURR1 overexpression (617 in F3, 1,580 in A4, and 240 in F5 cells) and 1,399 that were downregulated (214 in F3, 865 in A4, and 320 in F5 cells). The expression profiling data of three hNSC lines with NURR1 overexpression have been submitted to the GEO (accession No. GSE58475).

In comparing the three cell lines, we found that genes showing regulation in a coordinated manner (i.e., in all three lines and in the same direction) were not highly represented. For example, we found no genes among any of the three cell lines that were within 30 orders of the greatest magnitude of change by NURR1. We inferred from this outcome that NURR1 may modulate many genes indirectly rather than directly (Supplementary Fig. S1B).

Direct targets of NURR1 in hNSC lines

For an understanding of the logic of NURR1-dependent neurogenesis, the newly identified downstream genes must be placed within an underlying regulatory hierarchy. The transcription factor Nurr1 is characterized by binding as a monomer to the NBRE sequence motif, as a heterodimer with retinoid X receptor (RXR) to DR5, or as a dimer to NurRE (Perlmann and Jansson, 1995). To filter genes that NURR1 targets directly, based on the presence of NURR1 consensus binding sites, we searched promoter sequences (retrieved from the University of California, Santa Cruz at http://genome.ucsc.edu/cgi-bin/hgTables) within 3 kb of the transcriptional start site (TSS) for NURR1 binding elements with no more than one base pair difference compared to the consensus sites (Fig. 1A). Of the putative NURR1 direct targets, using the Eisen lab software, we detected 187 sequences representing 180 genes, based on their NCBI UniGene designations, that were predicted to modulated by NURR1 overexpression (Supplementary Table S2). A hierarchical cluster diagram is shown in Fig. 1C, and Venn diagrams depicting the number of genes commonly and differentially up- or downregulated in the three hNSC lines are shown in Figs. 2A and 2B.

Fig. 2

NURR1 overexpression in hNSCs identifies novel NURR1 direct target genes.

(A and B) Venn diagram showing the shared and unique genes upregulated (A) and downregulated (B) by NURR1 in three hNSC lines. (C) Heatmaps showing subclasses of genes coordinately regulated by NURR1 in three hNSC lines (derived from Fig. 1C); upper panel, genes upregulated by NURR1 in all three cell lines (23 genes); lower panel, genes downregulated by NURR1 in all three cell lines (15 genes). (D) The top 20 genes exhibiting the greatest magnitude of change by NURR1 in each cell line. All of these genes contained NURR1-binding site(s) in the promoter. Four genes are marked that were commonly ranked within the top 20 in the three cell lines: ● = C9orf16; ▶ = CNN1; ■ = NMU; ∇ = G0S2. All genes or expressed sequence tags were identified by the significance analysis of microarrays algorithm as being statistically significant. For all cell lines, the false discovery rate was set at 10%.

Based on the functions of the genes assigned by the Eisen lab software, NURR1 direct target genes were classified into categories (Fig. 1D). The classes of genes represented included nervous system development, cytoskeletal components, cell adhesion, DNA replication, regulation of progression through the cell cycle, transcription, and protein folding and secretion. In addition, we derived several important observations from this analysis. A significant proportion of the 180 differentially expressed genes were regulated in the same (coordinated) manner among the three hNSC lines, indicating transcriptional modulation that was independent of cell context (Fig. 1C, groups iv and v; Fig. 2C). Our microarray data also showed that several uncharacterized NURR1 targets, including chromosome 9 open reading frame 16 (C9orf16; highlighted in blue in Fig. 2D), are modulated in hNSC lines. Thus, NURR1 appears to modulate a number of characterized and uncharacterized genes, acting both as an activator and a repressor.

Genes regulated in a coordinated manner in all 3 lines

Next, to obtain a more restricted list of unambiguously regulated genes, we applied two stronger filters based on coordinated regulation of expression among the three cell lines (Fig. 2). First, we considered for further validation and quantitative analysis only those genes regulated in the same direction in all three NURR1-overexpressing hNSC lines (shown in Fig. 2C, Supplementary Tables S3 and S4). Of this group, 23 genes, including H2A.W histone (HIST3H2A), were upregulated, and 15 genes, including UDP-glucose ceramide glucosyltransferase (UGCG), were downregulated in all three lines. Almost all genes that were downregulated or upregulated in this coordinated manner showed a strong association with neuronal function, indicating a high stringency of our microarray data (Tables 1-3).

As a second step, we identified genes showing the greatest degree of change with NURR1 overexpression in each cell line (Fig. 2D). Upregulated genes that were within the 20th order of the greatest degree of change in all three lines were NMU, calponin 1 (CNN1), and C9orf16. The downregulated gene that met this criterion was G0S2.

Gene ontology (GO) and pathways of coordinately regulated genes

Because CNS development is still being investigated and novel functions and genes are always being added in this area, the GO analysis does not include the most recently published results. To identify in detail the function of each gene showing coordinated expression, we searched the possible functions of each candidate using UniProt (http://www.uniprot.org/), BioGPS, the Human Protein Atlas (http://www.proteinatlas.org/), Eukaryotic Linear Motif analysis, and PubMed. The results point to some striking features. In Tables 1-3, we present a detailed listing of these findings according to gene groups and pathways.

Overall, the analyses identified distinct categories of genes modulated by NURR1 in hNSCs, although they are not the well-known DA neuronal markers such as tyrosine hydroxylase (TH) and dopamine transporter (DAT). Among the gene sets showing coordinated upregulation by NURR1 in all three cell lines, we found a marked enrichment of synapse formation-associated dendritic morphogenesis categories as well as morphogenic movements of neuronal precursors during brain development, suggesting the importance of NURR1 in new synaptic formation and brain morphological development. Of the set of 38 genes included in this analysis, 13 of 23 upregulated genes (~56%) showed this enrichment, as did 1 of 15 downregulated genes (Table 1). Of this group, 10 genes, including sushi repeat containing protein X-linked 2 (SRPX2), localize in the synapse, primarily in postsynaptic elements (Table 1) (Sia et al., 2013; Soteros et al., 2018). In addition, fibronectin leucine-rich transmembrane protein 2 (FLRT2) was the most upregulated gene (3.5-fold) in the microarray data set for both HB1.F3 and HB1.F5 cells (Fig. 2D). Four genes, including FLRT2, regulate radial migration of pyramidal neurons and their tangential spread during embryonic brain development (Table 1) (Seiradake et al., 2014; Yamagishi et al., 2011). Of interest, C9orf16, the only uncharacterized gene, may function in synapse formation through trafficking of neurotransmitter-containing vesicles, based on its interaction with proteins in high-throughput proteomics analyses (UniProt site) (Huttlin et al., 2015; Marley and von Zastrow, 2010).

Table 1

Genes associated with synapse formation, neuronal morphogenesis and cell migration during brain development

Gene symbol	GenBank ID	Cellular localization	GO function	Fold change			P value		Reference
Synapse formation
SRPX2	NM_014467	Secreted, synapse, focal adhesion	Positive regulation of synapse assembly	0.5	0.7	0.6	0.000692	Excitatory synapse	(Sia et al., 2013; Soteros et al., 2018)
CYFIP2	AL161999	Ddendrites, synaptosome	Interaction with FMRP, a role in local protein translation at neuronal dendrites and in dendritic spine maturation, to initiate actin polymerization and branching	0.5	1.6	0.6	0.000244		(Han et al., 2015)
LTBP3	NM_021070	Extracellular regionMidbrain	Secretion, correct folding, and matrix deposition of TGF-b2Regulation of extracellular activation of TGF-b2Induction of midbrain DA neurons	0.6	2.2	1.0	0.001709	Inbibiroty synapse formation (TGFb1, CaMKII)	(Dobolyi & Palkovits, 2008)
BAALC	NM_024812	Postsynaptic lipid rafts	Synaptogenesis at postsynaptic lipid rafts by interacting with CAMK2A(Isoforms 2, 6: brain development)	1.3	1.0	1.1	0.000692	Inbibiroty synapse formation (TGFb1, CaMKII)	(Wang et al., 2005)
CAMK2N1	NM_018584	Sypaptosome/postsynaptic cell membrane (nervous system)	Potent and specific inhibitor of CAMK2	1.0	0.9	1.2	0.001994	Excitatory synapse	(Lucchesi et al., 2011; Saha et al., 2006)
C9orf16	NM_024112	Lysosome, endosome	Synpase formation, transmission with interaction with Dysbindin-D2Ra Post-endocytic sorting into lysosomea	0.5	2.3	1.6	0.000244	Coiled coil, uncharacterized protein families	(Marley & von Zastrow, 2010)
Neuronal morphogenesis
CNN1	NM_001299	Post-synaptic elements (cytoskeleton or focal adhesion)Growth cones in developing brain	Modulation of morphological structure of cytoskeleton in postsynpatic areasNeuronal and glial plasticity	0.6	2.4	1.6	0.001221	Interaction: part of cGMP kinase signaling complex at least composed of ACTA2, CNN1, PLN, PRKG1 and ITPR1	(Blunk et al., 2014; Liu & Jin, 2016; Rami et al., 2006)
ACTA2	NM_001613	Actin cytoskeleton	Cytoskeleton organization	0.5	3.2	0.4	0.000244	Interaction with CamkII	(D'Arco et al., 2018)
MYLK	NM_005965	Lamellipodium, stress fiber, focal adhesion	Actin filament organizatioPositive regulation of cell migrationGrowth initiation of astrocytic processes and in transmitter release at synapses	1.2	0.7	0.8	0.000244		(Li et al., 2016a)
BMCC1/PRUNE2	AB002365	CytoplasmCNS	BMCC1s (brain-isoform): cell morphologyNeurite outgrowth	0.8	3.2	0.8	0.000529		(Arama et al., 2012)
Cell migration and positioning during embryonic development
DPYSL3	NM_001387	F-actin, cell projection, growth cone, lamellipodiumStrong expression in fetal brain & spinal cord	For signaling by class 3 semaphorins and subsequent remodeling of the cytoskeletonAxon extension and guidancePositive regulator of convergent extension movement and correct positioning of certain type neurons during zebrafish neurulation	0.5	1.0	0.9	0.000244	Convergent extension movement and correct positioning	(Tan et al., 2015)
IVNS1ABP	NM_006469	Cytoplasm, cytoskeleton, nucleoplasm	Negative regulator of cell migrationNegative regulator of convergent extension in neurulation during embryonic development	–0.3	–0.7	–0.5	0.000244	Convergent extension movement and correct positioning	(Perconti et al., 2007)
FLRT2	NM_013231	Integral component of plasma membrane	Axon guidance/cell adhesion/regulation of neuron migraionRegulation of radial migration of neurons as well as their tangential spread	1.9	0.3	3.5	0.00293	Radial migration of pyramidal neurons as well as their tangential spread	(Seiradake et al., 2014; Yamagishi et al., 2011)
NMU	NM_006681	Terminal boutonSecreted	Positive regulator of cell migration, invasivenessPositive regulation of synaptic transmission	1.2	1.5	1.6	0.001831	Cell migration	(Kasper et al., 2016)

Predicted from interaction binding partners.

Moreover, nine genes were implicated in DNA replication and cell cycle categories. Eight of these were from a group of 15 downregulated genes (~53%) that were overrepresented among the coordinated NURR1-downregulated gene sets in the three cell lines. These results implicate NURR1 in cell cycle exit and in establishing post-mitotic DA neuroblasts during brain development (Table 2). In addition, nuclear and mitochondrially encoded genes, G0S2 and mitochondrially encoded NADH dehydrogenase 3 (MTND3), are modulated by NURR1 to increase oxidative phosphorylation (Lee et al., 2015; Yamada et al., 2013). mDA neurons are highly energy-demanding and have high requirements for ATP-generating oxidative phosphorylation. Our findings suggest that NURR1 is an important transcription factor for sustaining this distinguishing trait of DA neurons, consistent with previous results (Kadkhodaei et al., 2013). In this set of coordinated downregulated genes, we also detected leucine rich repeat containing 17 (LRRC17), UGCG, desmoplakin (DSP), and QKI, KH domain containing RNA binding (QKI), all of which occur in non-neuronal lineage differentiation categories, suggesting that NURR1 inhibits differentiation into non-neuronal lineages (Table 3).

Table 2

Genes associated with DNA replication, cell cycle progression, apoptosis and mitochondrion functions

Gene symbol	GenBank ID	Description	Cellular localization	GO function	Fold change			P value	Reference
MCM2	NM_004526	Minichromosome maintenance complex component 2	Nucleus, cytoplasm	Component of MCM2-7 complex which is replicative helicase essential for 'once per cell cycle' DNA replication initiation and elongation in eukaryotic cells	–0.4	–0.3	–0.2	0.000244	(Im et al., 2009)
MCM7	AF279900	Minichromosome maintenance complex component 7	Nucleus, cytoplasm		–0.3	–0.5	–0.3	0.000732	(Im et al., 2009)
MCM10	AB042719	Minichromosome maintenance complex component 10	Nucleus, cytoplasm	Acts as a replication initiation factor that brings together MCM2-7 helicase and DNA polymerasea/primase complex in order to initiate DNA replication	–0.1	–1.0	–0.4	0.000732	(Im et al., 2009)
PRC1	NM_003981	Protein regulator of cytokinesis 1	Nucleus, cytoplasm, spindle	Controlling spatiotemporal midzone formation and successful cytokinesis	–0.4	–0.6	–0.4	0.000244	(Polak et al., 2017)
ASF1B	NM_018154	Anti-silencing function 1B histone chaperone	Nucleus, nuclear chromatin	Cooperates with chromatin assembly factor 1 to promote replication dependent-chromatin assembly	–0.3	–0.9	–0.4	0.001221	(Jiangqiao et al., 2019)
PDAP1	NM_014891	PDGFA associated protein 1		Enhances PDGFA-stimulated cell growth	–0.9	–0.3	–0.3	0.000244	(Sharma et al., 2016)
DTL	AK001261	Denticleless E3 ubiquitin protein ligase homolog	Nuclear matrix	Decrese G2-arrest, p53 and p21 induction, and enhance cell proliferation	–0.5	–1.3	–0.4	0.000244	(Cui et al., 2019)
PFDN5	NM_002624	Prefoldin5/MM-1	Nucleus, cytoplasm Prefoldin complex	Inhibition of cell growth through repression of c-Myc activitInhibition of huntingtin/a-synuclein aggregationNecessary for proper folding of nascent proteins, in particular, tubulin and actin	0.4	0.6	0.3	0.000244	(Kadoyama et al., 2019)
G0S2	NM_015714	G0/G1 switch 2	Mitochondrion	Negative regulator of oxidative phosphorylationCell cycle switch from G0 to G1 phasePromotes apoptosis by binding to BCL2	–1.1	–0.9	–1.2	0.002889	(Lee et al., 2015; Yamada et al., 2013)
MTND3	NM_173710	Mitochondrially encoded NADH dehydrogenase 3	Mitochondrial inner membrane	Mitochondrial electron transport	0.3	0.9	0.4	0.000244	(Fu et al., 2019)

Table 3

Genes associated with differentiation into non-neuronal lineages

Gene symbol	GenBank ID	Cellular localization	GO function	Fold change			P value	Reference
UGCG	NM_003358	Membrane of Golgi apparatus	Myelin formation	–0.9	–0.3	-0.5	0.000244	(Watanabe et al., 2010)
DSP	NM_004415	Desmosome/adherens junction	Epidermis development	–1.4	–1.1	–0.6	0.000244	(Bharathan & Dickinson, 2019)
QKI	AF142421	Cytoplasm, nucleus	Regulator of oligodendrocyte differentiation and maturation in the brain	–0.5	–0.4	–0.7	0.000244	(Irie et al., 2016)
LRRC17	NM_005824	Extacellular space	Negative regulation of osteoclast differentiation	0.8	1.6	0.8	0.000244	(Kim et al., 2009)

Validation of NURR1 direct targets identified by microarray analysis

Of the 38 genes that were modulated in the same direction in the three hNSC lines (Fig. 2B, Supplementary Tables S3 and S4), 10 genes were chosen for validation by real time RT-PCR (Fig. 3). We selected eight genes involved in nervous system development: BAALC binder of MAP3K1 and KLF4 (BAALC) (Wang et al., 2005), calcium/calmodulin dependent protein kinase II inhibitor 1 (CAMK2N1) (Lucchesi et al., 2011; Saha et al., 2006), myosin light chain kinase (MYLK) (Li et al., 2016a), NMU (Kasper et al., 2016), C9orf16 (Huttlin et al., 2015; Marley and von Zastrow, 2010), prune homolog 2 with BCH domain (BMCC1) (Arama et al., 2012), FLRT2 (Seiradake et al., 2014; Yamagishi et al., 2011), and CNN1 (Blunk et al., 2014; Liu and Jin, 2016; Rami et al., 2006). The other two included genes, protein regulator of cytokinesis 1 (PRC1) (Polak et al., 2017) and G0S2 (Lee et al., 2015; Yamada et al., 2013), are involved in regulation of the cell cycle and in oxidative phosphorylation. Among the genes on this list, NMU, CNN1, C9orf16, and G0S2 were in the top 20 showing the greatest degree of change in all three lines. In all cases, we used real-time RT-PCR to confirm the direction of change in expression indicated by microarray (Figs. 3A-3C). In some cases, the array data suggested that NURR1 could modulate the same gene in opposite directions in different hNSC lines. Of these genes, we chose three for further analysis. Figure 3D shows real-time PCR data confirming this cell line–based directionality for decorin (DCN), pleiotrophin (PTN), and galanin and GMAP prepropeptide (GAL). In 30 of 39 (~77%) comparisons of gene expression using microarray and RT-PCR, results with the two methods showed agreement regarding the direction of change by cell line (Fig. 3).

Fig. 3

Real-time RT-PCR validation of selected NURR1 gene targets.

(A-C) Confirmation of NURR1-regulated genes in HB.F3, HB.FA4, and HB.F5 cells. RNA was isolated from 3 hNSC lines transduced with empty vector and NURR1-overexpressing hNSC lines, and assayed by real-time RT-PCR to validate the expression profile of NURR1-regulated transcripts. Gene expression levels were calculated by normalizing PCR data from vector-transduced and NURR1-overexpressed samples to the GAPDH control gene. Fold up- or downregulation was calculated using the ratio of normalized PCR values with the – ∆∆Ct method and the ratios of NURR1 overexpression/vector only signals from the microarray data. The left histograms show relative normalized expression from the microarray, and the right histograms show the relative normalized expression levels per RT-PCR. (D) Quantitative verification of NURR1 target genes with differential expression among different hNSC lines. Expression of three genes was analyzed using real-time RT-PCR, and fold changes in expression were similar to those from the array analyses (A-C). Asterisks (*) indicate genes showing agreement in the direction of expression change between two cell lines and between microarray and RT-PCR. Asterisks (**) indicate genes showing directional agreement among three cell lines and between microarray and RT-PCR.

NURR1 modulates target gene expression in a hNSC-specific manner

To discern which NURR1 target genes identified in this microarray analysis were involved only in hNSCs, we analyzed NURR1-associated gene expression in another cell type, human embryonic kidney HEK293 cells (Fig. 4A). In this cell line, expression of most target genes was not affected by increased levels of NURR1, except for NMU and BAALC, which showed 0.71- and –0.85-fold changes, respectively. This result contrasts with the direction of NURR1-associated gene expression change in hNSCs and suggests differential NURR1 functions in hNSCs versus HEK293 cells.

Fig. 4

NURR1 regulates target genes in a hNSC-specific manner.

NURR1 target expression in a human embryonic kidney cell line, HEK293 (HEK). Target gene levels were measured by real-time PCR in HEK293 cells overexpressing NURR1. Inset: RT-PCR analysis of HEK293 cells transfected with pLPCX (–) or pLPCX-NURR1 (+) plasmids. (B) Converse changes of NURR1 target genes in NURR1-knockdown hNSC lines. hNSCs were transfected with pLPCX-NURR1 (NURR1) or NURR1 siRNA (siNurr1). Target gene levels were assessed by real-time RT-PCR, and relative gene expression is given as the ratio of GAPDH-normalized values.

To functionally characterize the role of NURR1 in gene expression, we also inhibited NURR1 using siRNA (Fig. 4B). For this purpose, we focused on five genes (NMU, BAALC, CAMK2N1, FLRT2, and CNN1) that showed a consistent direction of change in expression levels on microarray and with RT-PCR (Fig. 3). Figure 4B demonstrates the successful inhibition of NURR1 using targeted siRNA compared to cells transfected with control siRNAs. The expression direction of the five genes was reversed in hNSCs expressing diminished levels of NURR1 compared to those overexpressing NURR1 (Fig. 4B). These data offer further evidence supporting the results of the real-time PCR analysis.

Transactivation of the human NMU and G0S2 promoters by NURR1

We next evaluated whether NURR1 may be involved in transcriptional regulation of human NMU and G0S2 genes. If NURR1 does regulate these genes during embryonic brain development, NURR1 binding elements may be evolutionarily conserved among species. We used the mVista (Frazer et al., 2004) program to identify evolutionarily conserved regions in the upstream sequences of the TSSs of human NMU and G0S2 that were upregulated or downregulated in the same direction in three cell lines (Fig. 5A). As noted above, NURR1 can bind as a monomer, a heterodimer, or dimer depending on the target. Two potential NURR1-binding elements (NBEs) were identified in each of the promoters: one NBRE and one NurRE in the NMU promoter, and two NBREs in the G0S2 promoter (Fig. 5A). Alignment of the human, rat, and mouse promoters revealed that NurRE in the human NMU promoter and NBRE2 in the human G0S2 promoter are perfectly conserved across all three species.

Fig. 5

Transcriptional regulation of human NMU and G0S2 genes by NURR1.

(A) mVista plot of 2-kb regions upstream of the TSS of human (h) NMU and G0S2 genes. The human genome was used as the base sequence for comparison with mouse and rat genomes. Peaks in the mVista plot represent conserved regions. Conserved and non-conserved binding sites for NURR1 are indicated in red and blue, respectively. (B) Sequence-specific binding activities of NURR1 to the conserved NurRE and NBRE2. EMSA was performed with radiolabeled oligonucleotides containing the NurRE of the human NMU promoter or the NBRE2 of the human G0S2 promoter using nuclear extracts from HB1.F3 cells overexpressing NURR1. Slow-migrating complexes (arrowheads) were detected in both promoters. Competitors (a 40- or 80-fold molar excess) or antibody against NURR1 was added to confirm specific NURR1 complexes. The DNA- binding activities were competed out by addition of unlabeled wild-type oligonucleotide (W), but not with addition of mutated NurRE (M). (C) ChIP analysis of NURR1 binding to human NMU or G0S2 promoter sequences. ChIP was performed with a NURR1 antibody and immunoglobulin G (IgG) as control in HB.F3 cells, and quantitative PCR analysis of the promoter sequences, as indicated. IP, input; N, no antibody added. (D and E) Reporter gene analysis of the human NMU-2102 promoter. Luciferase (LUC) activity assays were performed in hNSCs transfected with the constructs, as indicated. NURR1 binding sites are shown as gray boxes; gray boxes marked with X represent a mutation of the NurRE site that abolished the NURR1 binding. To compare NURR1 transactivation activities directly, we set the luciferase activity of the reporter construct with an empty vector (pLPCX) to 1. (F) Reporter gene analysis of the human G0S2 gene promoter from–2082 to +1 in 3 hNSC lines. Luciferase assays were performed three times with similar results. The figures represent the mean ± SD (bars), and data were normalized as in Figs. 5D and 5E. *P < 0.05, **P < 0.01, and ***P < 0.001.

We performed EMSA and a ChIP assay to verify that NURR1 is recruited to these sites in in vitro and in vivo (Figs. 5B and 5C). EMSA showed that NURR1 could bind evolutionarily conserved NBEs in both promoters in vitro (Fig. 5B). In addition, competition experiments and supershift assay confirmed that NURR1 specifically bound NurRE in the NMU promoter and NBRE2 elements in the G0S2 promoter (Fig. 5B). For the ChIP assay, HB1.F3 cells overexpressing NURR1 were cross-linked to chromatin and immunoprecipitated with an antibody against NURR1, followed by PCR using primers corresponding to NurRE or NBRE2. As shown in Fig. 5C, NURR1 was recruited to the NBEs in both promoters.

To verify whether the predicted binding sites for NURR1 are functional, we cloned a 2102-bp fragment of the human genomic sequence upstream of the NMU translation start site into a luciferase reporter plasmid. Transient transfection showed that NURR1 overexpression increased transcription from this promoter by 7.1-, 13.0-, and 3.5-fold in F3, A4, and F5 cells, respectively (Fig. 5D). To investigate whether the evolutionarily conserved NurRE is important for NURR1-dependent activation of the human NMU promoter, we used site-directed mutagenesis to introduce point mutations in the promoter that disrupt these sequences. These mutations clearly reduced the ability of NURR1 to activate the promoter (approximately 66.6%, 56.8%, and 31.3% decreases in F3, A4, and F5 cells, respectively) (Fig. 5E). The NurRE element thus is essential for NURR1-dependent NMU expression in hNSCs. We also evaluated the transcriptional regulation of the G0S2 promoter by NURR1 (Fig. 5F). A 2082-bp fragment of the human G0S2 promoter was cloned into a luciferase reporter plasmid. As shown in Fig. 5F, NURR1 comprehensively repressed human G0S2 promoter activity in all three hNSC lines in a dose-dependent manner. Collectively, the effect of NURR1 on the promoter-driven reporter activity in the NMU and G0S2 genes paralleled changes in mRNA levels detected by microarray and RT-PCR.

Expression of Nurr1 direct target genes from the ventral midbrain of embryonic mice during in vivo mDA differentiation

To confirm the in vivo expression pattern of these gene sets during midbrain development, we performed quantitative PCR analysis of mRNA extracted from the mouse embryonic ventral midbrain (Fig. 6). Because Nurr1 expression becomes evident at E10.5 during mDA cell fate specification (Fig. 6A) (Perlmann and Wallen-Mackenzie, 2004), we prepared embryonic midbrain tissues from E9.5 to E13.5 for RT-PCR analysis (Figs. 6B and 6C). We detected an increase in expression of Nmu, C9orf16, dihydropyrimidinase-like 3 (Dpysl3), G0s2, and Baalc at successive intervals from E9.5 to E13.5 (Figs. 6B and 6C). In addition, Camk2n1 and Cnn1 expression became evident at E10.5 and was maintained thereafter (Figs. 6B and 6C). Figure 6D shows a summary of the timing of target gene expression. These results validate our microarray data and indicate that these genes are truly in vivo targets of Nurr1.

Fig. 6

Expression of Nurr1 target genes in embryonic midbrain tissues.

(A-a) Schematic illustrating mouse brains ranging from E9.5 to E13.5. The red line marks the plane of section shown in immunocytochemistry. (A-b) Specification of the mDA neuronal identity within the ventral midline (VM). (A-c) Immunohistochemical analysis of Nurr1 and TH expression in midbrain sections from E12.5 mice. (B and C) Representative RT-PCR analysis for Nurr1 target genes from embryonic midbrain tissues during mDA neuron development. DAPI, 4’,6-diamidino-2-phenylindole. (D) The timing of expression of Nurr1 direct target genes in the midbrain of mice at E9.5 to E13.5. Tel, telencephalon; Di, diencephalon; Rhomb, rhombencephalon.

BioGPS analysis: CNS-specific expression of NURR1 direct target genes

To discern which Nurr1 target genes identified in our microarray analysis are biologically relevant, we searched each gene on the BioGPS website (http://biogps.org), using the GeneAtlas U133A gene expression (Wu et al., 2016). The most highly represented systems were CNS tissues, and other systems did not show a significant percentage of tissue-specific genes. As indicated in Table 4, 12 transcripts (31.6%) exhibited ≥ 1.5 fold increased expression in the CNS or fetal brain as compared to the mean expression values of all human tissues in the datasets. The gene symbols and CNS and fetal brain enrichment for these 12 genes are provided in Table 4, and additional details on the microarray probe ID and CNS versus average human expression levels are given in Supplementary Table S5. Of these 12 genes, three (BAALC, CAMKIIN1, and QKI) were highly expressed, with an expression value >30 multiples of the median in the brain/CNS, in accordance with previously established stringency thresholds (Pennings et al., 2011). The other dominant genes were DPYSL3 and FLRT2, both of which were specifically expressed in fetal brains.

Table 4

BioGPS gene expression analysis of Nurr1 direct target genes

Gene name	Gene symbol	CNS-enriched (fold change)	FB-enriched (fold change)	CNS expression
Upregulated genes
Calcium/calmodulin dependent protein kinase II inhibitor I	CAMK2N1	69.63a	56.92	Prefrontal cortex, amygdala, whole brain
Brain and acute leukemia, cytoplasmic	BAALC	48.25a	40.16	Spinal cord, hypothalamus, amygdala, prefrontal cortex, whole brain
Dihydropyrimidinase-related protein 3	DPYSL3	18.65	156.10b	Fetal brain
Protein prune homolog 2	Prune2/BMCC1	22.01	3.03c	Prefrontal cortex, hypothalamus
Fibronectin leucine rich transmembrane protein 2	FLRT2	7.32	25.27b	Fetal brain, prefrontal cortex
Cytoplasmic FMR1 interacting protein	CYFIP2	10.02	5.20	Amygdala, prefrontal cortex, pineal
Chromosome 9 open reading frame 16	C9orf16	2.03	2.91	Whole brain (6.5×), temporal lobe (5×), amygdala, prefrontal cortex, fetal brain (3-4×)
Latent transforming growth factor beta binding protein 3	LTBP3	3.67	0.99	Retina, pineal
Prefoldin 5	PFDN5	1.20	1.63	Pineal, fetal brain
Histone 3, H2a	HIST3H2A	1.25	2.05	Cerebellum peduncles, cerebellum, fetal brain
Downregulated genes
Quaking homolog, KH domain RNA binding (mouse)	QKI	38.96a	15.84c	Spinal cord, hypothalamus, prefrontal cortex
Influenza virus NS1A binding protein	IVNS1ABP	2.00	9.12b	Fetal brain

Genes are associated with CNS-specific expression > 30MoMs in BioGPS.

Genes are expressed with the b highest or the c lowest level in the fetal brain of CNS tissues.

GWAS of NURR1 direct target genes

NURR1 is a Parkinson’s disease susceptibility gene because it is a key factor in the neurodevelopmental control of DA neurons in the substantia nigra pars compacta and ventral tegmental area. In addition, NURR1 mutations are associated with both schizophrenia and bipolar disorder (Buervenich et al., 2000; Vuillermot et al., 2011; Xing et al., 2006). Considering a possible role for NURR1 direct target genes in DA neurogenesis, we hypothesized that genetic alterations within the loci of these genes (including non-coding DNA sequences and regulatory elements) may be associated with susceptibility to schizophrenia, bipolar disorder, and depression, as well as Parkinson’s disease. For this reason, we searched for SNPs associated with susceptibilities to these conditions, focusing on the loci of NURR1 direct target genes from the previously identified variants of GWAS datasets derived using different chip platforms (Li et al., 2012; 2016b).

About half of the up- and downregulated targets showed associations with these diseases (Supplementary Table S6). Of these targets, we identified 11 SNPs with Parkinson’s disease susceptibilities at seven independent loci: rs1385331 (chr 3) for phosphatidylinositol-glycan biosynthesis class X (PIGX); rs10056132 (chr 5) and rs10041339 (chr 5) for DPYSL3; rs4874150 (chr 8) for tissue specific transplantation antigen P35B (TSTA3); rs2453998 (chr 8) for BAALC; rs3003602 (chr 9) and rs2502731 (chr 9) for C9orf16; rs1011711 (chr 9) for PRUNE; and rs5966709 (chr X), rs4828037 (chr X), and rs932437 (chr X) for SRPX2 (Table 5). All of the variations were found in non-coding DNA within ~50 kb of NURR1 target genes. Because a non-coding DNA sequence at the loci of direct target genes may encode information guiding the proper level, time, and place of gene expression during DA neurogenesis (Lee et al., 2020), we searched the genomic region for each variant associated with Parkinson’s disease, and annotated this region for sequence conservation and candidate enhancers based on H3K4me1 and H3K27ac histone modifications and DNase hypersensitive sites (Fig. 7). Of 11 SNPs with Parkinson’s disease susceptibilities, two variations in the C9orf16 locus (rs3003602 and rs2502731), two in DPYSL3 (rs10056132 and rs10041339), one in BAALC (rs2453998), and one in PIGX (rs2453998) directly overlapped with a H3K4me1 and a H3K27ac marker of regulatory element sequences (candidate enhancer). These variants were not evolutionarily conserved except in primates, and we did not detect any overlap between the other five polymorphisms and any regulatory feature. These GWAS data strongly support the inference that these polymorphisms are likely to disrupt an enhancer of C9orf16, DPYSL3, BAALC, and PIGX genes. In addition, there were many SNPs associated with schizophrenia, bipolar disorder, and depression, which are implicated in DA neuronal activity in the CNS (Supplementary Table S6).

Table 5

Genetic Associations for Parkinson’s Disease at Nurr1 direct target genes

SNP	Chr: position (hg19)	Locus	Risk/non-risk allele	P valuea	SNP type	PubMed ID (No. of independent studies)
rs1385331	3: 196351215	PIGX	T/C	7.78 × 10–5	Upstream	23936387, 17052657 (2)
rs10056132	5: 146831553	DPYSL3	G/A	3.00 × 10–6	Intron	21876681, 24665060 (2)
rs10041339	5: 146949107	DPYSL3	G/C	4.00 × 10–6	Intergenic	21876681 (1)
rs2453998	8: 104089618	BAALC	A/C	2.98 × 10–4	Upstream	24023788, 17052657, pha002887 (3)
rs4874150	8: 144636272	TSTA3	A/G	7.68 × 10–4	Intron	17052657 (1)
rs1011711	9: 79550279	Prune2	T/C	5.39 × 10–4	Upstream	20877124, 16252231 (2)
rs3003602	9: 130981064	C9orf16	T/C	1.00 × 10–5	Intron	19772629 (1)
rs2502731	9: 130976557	C9orf16	T/C	2.00 × 10–6	Intron	18839057, 19772629, 17255346 (3)
rs5966709	X: 99844506	SRPX2	T/G	1.55 × 10–4	Intron	17052657 (1)
rs4828037	X: 99845684	SRPX2	C/T	1.59 × 10–4	Intron	17052657 (1)
rs932437	X: 99837874	SRPX2	C/G	2.23 × 10–4	nearGene-5	17052657 (1)

For each locus, the SNP, genomic location (chromosome, position [hg19]), risk/non-risk allele, and association P value are shown.

The most significant P value.

Fig. 7

Genetic associations for Parkinson’s disease at the DPYSL3, C9orf16, BAALC, and PIGX loci.

(A) The genomic region around the DPYSL3 gene showing locations of all known common (> 10%) SNPs, sequence conservation, and enhancer marks from human cell lines (DNase I hypersensitive sites, H3K4me1, H3K27ac). All SNPs studied are plotted against their hg38 genomic location. The Parkinson’s disease-associated SNPs rs10056132 and rs10041339 are highlighted in red. DPYSL3 and additional genes near the peak SNPs rs10056132 and rs10041339 are depicted at the bottom of the window with arrows indicating the direction of transcription. (B) The peak SNPs rs3003602 and rs2502731 near the C9orf16 gene. (C) The peak SNP rs2453998 near the BAALC gene. (D) The peak SNP rs1385331 near the PIGX gene.

DISCUSSION

Human brain development is very different from mouse brain development (Bernard et al., 2012; Bohland et al., 2010; La Manno et al., 2016; Lein et al., 2007; Myers et al., 2015; Ng et al., 2009). Previous reports have described gene expression profiling by the DA transcription factor Nurr1 in a rodent system (Jacobs et al., 2009a; Sousa et al., 2007) but did not confirm whether the modulated genes are direct targets of Nurr1. We addressed this question by identifying NURR1 direct target genes in human neural progenitors.

Our analyses have yielded several important findings. First, several of the 38 genes differentially and directly affected by NURR1 were regulated in the same (coordinated) manner among the three hNSC lines, indicating that NURR1 transcriptional regulation occurs independent of cell context. The results also suggest that NURR1 modulates many genes, acting both as an activator and a repressor. These findings contrast with earlier reports (Kadkhodaei et al., 2013; Panman et al., 2014) describing Nurr1 as functioning mainly as a transcriptional activator to regulate a battery of genes expressed in DA neurons. This discrepancy may trace to the fact that these authors analyzed Nurr1 ablation of already mature and postmitotic DA neurons.

Moreover, genes associated with synapse formation (6 genes), neuronal morphogenesis (4 genes), and cell migration during brain development (4 genes) were highly upregulated in the current study, whereas genes involved in cell cycle progression and DNA replication (9 genes) were downregulated. In the E10.5-E11.5 mouse embryo, a series of events is important for mesencephalic DA cell differentiation (Wallen et al., 1999). Specifically, DA progenitor cells migrate from the ventricular zone to the mantle zone, exit the cell cycle (become post-mitotic), and start to innervate the striatal target area. Simultaneously, Nurr1 expression is strongly upregulated. Our microarray data clearly fill in gaps from previous reports showing that in contrast to wild-type animals, in Nurr1–/– mutant mice, DA progenitor cells in the midbrain cannot migrate laterally, fail to innervate the striatal target area, and become apoptotic (Wallen et al., 1999). These authors did not directly demonstrate that the reason for these failures is a lapse in regulation of genes involved in migration, synaptogenesis of DA progenitor cells, and cell cycle exit. To date, a direct relationship of Nurr1 expression with other DA differentiation–related events has not been demonstrated. Our findings fill in these gaps, providing direct evidence that NURR1 selectively contributes to the expression of genes involved in the correct positioning of migrating neurons and in synapse formation with correct partners, as well as ensuring postmitotic status during mDA neuron development.

We also found that genes encoding proteins involved in regulating the DA neurotransmitter phenotype and expressed in the final stage of DA differentiation, including TH, vesicular monoamine transporter, dopa decarboxylase, and DAT, were not upregulated by NURR1 in our hNSC lines (Supplementary Table S2). This result is consistent with our earlier work (Kim et al., 2013) and with another previous report (Kadkhodaei et al., 2013). An explanation may lie in the timing of expression relative to the window of development we examined here. Our hNSC lines show gene expression profiling similar to what is seen from E9.5 to E10.5 (Supplementary Table S1), with in vivo Nurr1 expression becoming evident in the midbrain and DA precursors shifting to a postmitotic DA neuroblast state. These other well-known markers of DA neurons may not have shown expression changes because the timing of their action is outside of this window. Nevertheless, we did identify a large number of neurogenesis-related genes as direct targets of NURR1, all candidates that, to the best of our knowledge, have not yet been described previously as downstream of NURR1. The restricted modulation pattern of these genes reveals potential roles for NURR1 in brain development, signaling, metabolism, and disease. Taken together, our findings are consistent with published reports that in adult DA neurons, Nurr1 ablation is associated with a progressive pathology involving reduced striatal DA, impaired motor behaviors, and dystrophic axons and dendrites, as well as downregulation of oxidative phosphorylation–related genes in Parkinson’s disease (Kadkhodaei et al., 2013).

Supplemental Materials

Note: Supplementary information is available on the Molecules and Cells website (www.molcells.org).