RNA-seq of Isolated Chromaffin Cells Highlights the Role of Sex-Linked and Imprinted Genes in Adrenal Medulla Development

Adrenal chromaffin cells and sympathetic neurons synthesize and release catecholamines, and both cell types are derived from neural crest precursors. However, they have different developmental histories, with sympathetic neurons derived directly from neural crest precursors while adrenal chromaffin cells arise from neural crest-derived cells that express Schwann cell markers. We have sought to identify the genes, including imprinted genes, which regulate the development of the two cell types in mice. We developed a method of separating the two cell types as early as E12.5, using differences in expression of enhanced yellow fluorescent protein driven from the tyrosine hydroxylase gene, and then used RNA sequencing to confirm the characteristic molecular signatures of the two cell types. We identified genes differentially expressed by adrenal chromaffin cells and sympathetic neurons. Deletion of a gene highly expressed by adrenal chromaffin cells, NIK-related kinase, a gene on the X-chromosome, results in reduced expression of adrenaline-synthesizing enzyme, phenyl-N-methyl transferase, by adrenal chromaffin cells and changes in cell cycle dynamics. Finally, many imprinted genes are up-regulated in chromaffin cells and may play key roles in their development.

their phenotypic similarity and low numbers. In our previous study 17 , TH and CART were revealed as potential makers for the early separation of sympathetic neuroblasts and adrenal chromaffin cells in E12.5 mice because of their differential expression levels; sympathetic neuroblasts cells are relatively low in TH-immunoreactivity while the adrenal chromaffin cells have significantly higher TH-immunoreactivity, reflecting differences in Th RNA expression 4 . In addition, only sympathetic neuroblasts are immunoreactive for the neuropeptide, CART (Cocaine and Amphetamine Regulated Transcript) from E12.5 to E13.5. Therefore, in the present study we used TH-Cre activation of enhanced yellow fluorescent protein (EYFP) expression in transgenic mice coupled with fluorescence-activated cell sorting (FACS) to isolate and collect sufficient number of sympathetic neuroblasts and adrenal chromaffin cells at E12.5 for transcriptomic analysis by RNA sequencing. This allowed the assessment of all differentially expressed genes, and the identification of potentially important transcription and cell signaling genes. Subsequent studies tested the leading candidate gene for a role in chromaffin cell development along with assessing the expression of imprinted genes.

Differential EYFP Expression in Sympathetic Neuroblasts and Adrenal Chromaffin Cells.
We have shown that TH immunoreactivity in developing chromaffin cells is significantly higher than in sympathetic neuroblasts 17 . We sought to separate developing chromaffin cells from sympathetic neuroblasts based on this difference using TH-Cre::R26R-EYFP reporter mice. In E13.5 mice (Fig. 1A-E), where developing chromaffin cells and sympathetic neuroblasts were anatomically distinct, surprisingly the native EYFP signal (and EYFP immunoreactivity seen using a green fluorescent protein antiserum) in the adrenal gland anlagen was weaker than in the suprarenal and other prevertebral ganglia (Fig. 1E), the inverse of the staining intensity difference seen with antisera to TH 17 . In E12.5 TH-Cre::R26R-EYFP mice ( Fig. 1F-J), where anatomical boundaries between developing chromaffin cells and sympathetic neuroblasts were much less distinct, there was also heterogeneity in both native EYFP and EYFP immunoreactivity. EYFP+ cells with both high and low levels of expression were usually intermingled without clear anatomical boundaries.
We then examined whether cells expressing high levels of EYFP from the TH transgene (EYFP+ Hi ) corresponded to sympathetic neuroblasts while cells expressing low levels of EYFP (EYFP+ Lo ) corresponded to developing chromaffin cells. We quantified and plotted the relative fluorescence intensity for TH-IR against EYFP-IR for each cell in the abdominal region of E13.5 and E12.5 TH-Cre::R26R-EYFP mice (Fig. 2). The distributions of both TH and EYFP-IR immunofluorescence at E13.5 ( Fig. 2A) appeared largely bimodal, with adrenal chromaffin cells clustering in the lower right region of the graph (TH-IR Hi /EYFP-IR Lo cells) and sympathetic neuroblasts from the suprarenal ganglia clustered in the upper left region of the graph (TH-IR Lo /EYFP-IR Hi cells). Chan et al. 17 showed that immunoreactivity to CART was specific to sympathetic neuroblasts and not present in developing chromaffin cells. CART-immunoreactivity was present almost exclusively in TH-IR Lo /EYFP-IR Hi cells ( Fig. 2A), confirming that the levels of EYFP immunofluorescence can be used to separate the two cell types. A handful of cells, previously identified as extra-adrenal chromaffin or paraganglionic cells 17 , expressed CART-IR in addition to high levels of EYFP-IR and TH-IR ( Fig. 1A-D).
At E12.5, where the anatomy cannot be reliably used to separate the developing chromaffin cells from sympathetic neuroblasts, there were again two distinct types of EYFP+ cells; one with high intensity of EYFP-IR (EYFP-IR Hi ) and the other low intensity of EYFP-IR (EYFP-IR Lo ), intermixed around the dorsal aorta ( Fig. 1F-J). When EYFP-IR was plotted against TH-IR for each cell, they segregated into two populations, as seen in E13.5 embryos (Fig. 2B). Again, when CART-IR (a marker of sympathetic neuroblasts) was considered, it was expressed predominantly by cells with a TH-IR Lo /EYFP-IR Hi phenotype that lay in the upper left part of the distribution (Fig. 2B). These data show that separation of developing adrenal chromaffin cells and sympathetic neuroblasts is possible in TH-Cre::R26R-EYFP mice as early as E12.5 on the basis of their differential EYFP expression. This was confirmed below, using RNA sequencing analysis after FACS isolation.

FACS Isolation of Live Sympathetic Neuroblasts and Adrenal Chromaffin Cells.
Cells from E11.5, 12.5, and 13.5 TH-Cre::R26R-EYFP mice were subjected to FACS based on native EYFP fluorescence intensity and 7-AAD staining to exclude dead or damaged cells (Fig. 3A) and then on side-scattered light (SSC, Fig. 3B-D). Consistent with our hypothesis that differences in EYFP intensity can separate developing adrenal chromaffin cells from sympathetic neuroblasts in embryonic TH-Cre::R26R-EYFP mice, E12.5 was the earliest stage at which two different populations of EYFP+ cells could be separated (Fig. 3B,C). Careful gating at E12.5 and E13.5 in the FACS eliminated the small number of cells that lay intermediate between the two populations. Gating also eliminated cells showing very high intensity in both EYFP and SSC, as they likely represent the presumptive paraganglionic cells. At E11.5, EYFP+ cells formed a homogeneous population without obvious segregation into EYFP+ Hi and EYFP+ Lo populations (Fig. 3C). Our approach enables separation of progenitor cells that give rise to chromaffin cells and sympathetic neurons as early as E12.5 and was used to identify candidate genes involved in the development of the two cell types.
RNA sequencing analysis. The transcriptomic profiles of sympathetic neuroblasts and adrenal chromaffin cells were determined by RNA sequencing in four biological replicates of FACS-isolated adrenal chromaffin precursors and sympathetic neuroblasts from E12.5 TH-Cre::R26R-EYFP mice. Sequencing yielded more than 30 million 50 base pair single-end informative reads per sample and 17,169 annotated genes were detected. Principal component analysis indicates the differences in transcriptomic profiles were mainly due to cell types, confirming a good segregation of adrenal chromaffin cells and sympathetic neuroblasts (Fig. 4A). When the differential transcriptome between adrenal chromaffin precursor cells and sympathetic neuroblasts was generated (Fig. 4B), 4,786 genes were differentially expressed with log2 fold change >1, and both false discovery rate (FDR) and Immunostaining of transverse sections through the adrenal region of TH-Cre::R26R-EYFP mouse embryos at E13.5 (A-E) and E12.5 (F-J). A shows the native EYFP (yellow) signal after fixation of TH-Cre::R26R-EYFP mouse embryos at E13.5, the prevertebral suprarenal ganglion (solid line) and the adrenal medulla (dashed line) marked. EYFP-immunoreactivity for the same section is shown in (B), THimmunoreactivity in (C) and CART-immunoreactivity in (D). (E) Is a merge of images (B,C). Note that TH immunoreactivity shows the reverse pattern of intensity to both native EYFP and EYFP-immunoreactivity. (F-J) is an equivalent region from an E12.5 embryo as (A-E). The dorsal aorta (a) is indicated. Note that differential expression of TH-driven EYFP was observed in that some TH-expressing cells were brighter in EYFP than the others, but there was no obvious anatomical segregation of cells differentially expressing EYFP. www.nature.com/scientificreports www.nature.com/scientificreports/ p-value < 0.05. Fold change is calculated by the gene expression level in adrenal chromaffin cells divided by the gene expression level in sympathetic neuroblasts.
Droplet digital PCR (ddPCR) was used to confirm the RNA sequencing results. Thirteen genes that covered a range from high to low expression levels and fold changes were tested. For all of the genes, mRNA expression patterns by ddPCR were highly consistent with RNA sequencing results (Fig. S1) in that the fold changes of expression in adrenal chromaffin cells relative to sympathetic neuroblasts revealed by both methods were strongly correlated in linear regression analysis (r = 0.989).
The RNA sequencing data confirmed the cell type-specific gene expression of markers of mature sympathetic neurons and adrenal chromaffin cells (Fig. 4C,D). Messenger RNA for known adrenal chromaffin cell markers, expressed between 6 and 24-fold greater in chromaffin cells (Fig. 4C), included Dbh (dopamine beta hydroxylase) 18,19 , Chga (chromagranin A) and Chgb (chromagranin B) 20 and Vmat1 or Slc18a (vesicle monoamine transporter 1 or solute carrier family 18 member A1) 10,21,22 . Known markers of sympathetic neurons expressed at between 2.5 and 4-fold difference in sympathetic neuroblasts (Fig. 4D) included Nefl (neurofilament, light) 9,19,23 , Nefm The upper left quadrant of each plot contained TH-IR Lo /EYFP-IR Hi cells while lower right quadrant contained TH-IR Hi /EYFP-IR Lo cells. For E13.5, cells are identified and categorized into sympathetic neuroblasts (○), paraganglia (△) and adrenal chromaffin cells (◊) based on anatomical segregation into adrenal medulla or sympathetic ganglia. On E12.5, it was not possible to discriminate cell type by anatomical organisation and the cells are unclassified. However, the pattern observed in the plot is identical to that seen at E13.5, leading to the conclusion that cells that are TH-IR Lo /EYFP-IR Hi and which also express CART-IR are likely to be sympathetic neuroblasts while TH-IR Hi /EYFP-IR Lo cells lacking CART-IR are likely to be adrenal chromaffin cells. For E13.5, n = 3 embryos, 521 cells and for E12.5, n = 2 embryos, 505 cells.
Gene Ontology and Pathway Analysis. Gene ontology (GO) analysis classified the transcriptome of E12.5 adrenal chromaffin cell and sympathetic neuroblasts into 17 and 28 GO functional terms respectively (Table 1). For the adrenal chromaffin cells, genes were abundant in biological processes that included; pattern specification processes, cell-cell signaling, system development, cell-cell adhesion and synaptic vesicle exocytosis. Among these, "pattern specification -dorsal/ventral axis specification" was the most highly over-represented term with 4.54-fold enrichment. For the transcriptome of sympathetic neuroblasts, biological processes such as chromatin organization, cell cycle, DNA metabolic process and regulation of gene expression -epigenetic were overrepresented. "Organelle organization -chromatin assembly" was the most enriched term with 8.46-fold. Genes that underlie adrenal chromaffin cell and sympathetic neuroblast development. In addition to the GO enrichment and pathway analysis that consider mainly the fold difference, our analysis also took into account of the transcript abundance for prioritizing genes of interest. For the 4,786 genes expressed at more than 2-fold difference across the two cell types and with FDR cutoff at 0.05, RNA sequencing provided a comparison of the relative numbers of copies present (as counts per million transcripts -CPM) and of the relative expression levels (fold difference) between adrenal chromaffin cells and sympathetic neuroblasts. Although it is difficult to predict from the number of copies alone whether a gene actively plays a role in development, genes that are strongly differentially expressed by one of the two cell types may be important genes in the development of that cell type. We therefore ranked genes by fold difference and in our initial analysis examined only genes with a CPM of 40 or more (33.55% of total genes detected). To confirm whether the expression of any gene selected was peaking at the adrenal chromaffin cell stage, rather than in bridge cells or in Schwann cell precursors, we examined the raw data (Accession No. GEO99933 at NCBI) by Furlan et al. 4 . Finally, attention was paid to transcription factors and genes forming signaling pathways as they are most likely to direct differentiation. www.nature.com/scientificreports www.nature.com/scientificreports/ In adrenal chromaffin cells, of the 2,938 genes expressed at more than 2-fold difference and with FDR cutoff at 0.05, we identified five transcription factors (Elf3, Elf4, Foxq1, Bhlhe40 and Trim16) and 13 genes likely to be involved in cell signaling (3 G-protein receptors; Gipr, Casr and Gpr139, 5 receptors; Acvrl1, Baiap3, Dll4, Epha4 www.nature.com/scientificreports www.nature.com/scientificreports/ and Procr, a protein kinase; Nrk, a PI3k pathway member; Dapp1, a connexin; Gjb5, a secreted protein; Wnt9a and a regulator of G-protein signaling; Rgs5) ( Table 2) as meeting these selection criteria. All of these genes were differentially highly expressed (130 to 9-fold) and, based on the raw data of Furlan et al. 4 showed peak expression in developing adrenal chromaffin cells rather than bridge cells. None has been associated previously with chromaffin cell development 9,10 .
Among the 1,848 genes that were more highly expressed in sympathetic neuroblasts (Table 3) were five transcription factors (Ebf1, Aff3, Hoxc8, E2f7 and Foxm1), and 13 genes likely to be involved in cell signaling (1 receptor; Islr2, 3 signaling molecules; Vip, Fgf13 and Bmp3, 5 genes involved in intracellular signaling; Shcbp1, Crabbp1, Socs2, Lrr1 and Fam83d and 4 genes potentially involved in Rho GTPase signaling; Arhgap22, Iqgap3, Depdc1b, Racgap1). Many of the remaining genes highly expressed in sympathetic neuroblasts are involved in regulation of the cell cycle or the cytoskeleton. Potential Regulators are Transiently and Differentially Expressed During the Key Stage for Sympathoadrenal Development. We then further investigated the temporal expression patterns of four genes by ddPCR gene expression analysis (Fig. 5). We chose Dlk1, which is expressed at very high levels (9,286 CPM) by developing chromaffin cells. The gene has been identified previously as being associated with developing chromaffin cells 9,10 but it is not known whether it is developmentally regulated. From the novel genes in Table 2 we chose Elf3 and Foxq1 as the transcription factors with the highest fold difference, and Nrk as a novel and highly expressed kinase in the JNK pathway 32 . All of the genes were transiently expressed during chromaffin cell development with highest expression on E12.5 and declined thereafter, although Dlk1 also had a second peak of expression on E17.5 before declining in expression for a second time. In every case, the expression of all of the genes was very low in sympathetic neuroblasts at all ages examined.

Nrk-deficiency in Mice Impaired the Acquisition of an Adrenergic Phenotype by Chromaffin Cells and altered Cell Proliferation in Both Cell types.
Among the potential regulators of chromaffin cell development, Nrk is expressed at the highest level in E12.5 chromaffin cells (Table 2), more than 54-fold higher than in sympathetic neuroblasts. To examine the function of Nrk in chromaffin cell development, Nrk-deficient mice 33 were examined at E18.75, when development of the adrenal medulla and sympathetic ganglia are mostly complete. E18.75 Nrk mutants were generated by crossing heterozygous Nrk mutant female mice (Nrk-het) to Nrk-null mutant male mice. As a result, hemizygous male (Nrk +/Y ) offspring are wild-type; heterozygous mutant females (Nrk +/− ) contain approximately 50% of Nrk-null cells due to the random effect of X-inactivation in the wild-type Nrk allele; homozygous mutant females (Nrk −/− ) and hemizygous mutant males (Nrk −/Y ) are Nrk-null. Four embryos of each genotype across two litters were examined.
In Nrk mutant mice, the size, appearance and anatomical position of the adrenal glands were not noticeably different from wild-type mice. The numbers of adrenal chromaffin cells and sympathetic neuroblasts ( Fig. 6A-I) were examined in transverse sections of the upper abdomen at the adrenal level. Both cell types expressed TH, but only adrenal chromaffin cells were immunoreactive for phenylethanolamine N-methyltransferase (PNMT), the enzyme that generates adrenaline from noradrenaline. Nrk-null and Nrk-het mutant mice expressed TH in sympathetic neurons of the suprarenal ganglion and adrenal medullary chromaffin cells as in wild-type mice. Subpopulations of the TH-IR+ adrenal medullary chromaffin cells also co-expressed PNMT in both wild-type and Nrk mutant mice. The number of PNMT-IR+ cells as a proportion of TH-IR+ adrenal chromaffin cells was compared for each genotype (Fig. 6J). In wild-type mice, 72.2% of adrenal chromaffin cells were PNMT-IR+ (adrenergic), while in Nrk-null mice, only 41.1% of adrenal chromaffin cells was PNMT-IR+ (Fig. 6J). In Nrk-het mice, the proportion of adrenergic chromaffin cells was intermediate between the wild-type and Nrk-null with 65.4% of TH-IR+ cells in the adrenal medulla being PNMT-IR+. The apparent intensity of PNMT-IR in adrenal chromaffin cells varied in the same way, being highest in wild-type mice. The difference in proportions of PNMT-IR+ cells was significant (Chi-squared test, X 2 , 762, p = 0.05, df = 2, N = 10833, critical value = 5.99), confirming that the proportions of PNMT+ chromaffin cells among the wild-type, Nrk-het and Nrk-null mice are not equal. In order to test where any significant differences lay, the 2 × 3 contingency table was subdivided 34 into 3, 2 × 2 tables (wild-type versus Nrk-null, wild-type versus Nrk-het and Nrk-het versus Nrk-null), which revealed that Nrk-null mice were significantly different from both wild-type and Nrk-het mice (Fig. 6J).
These results demonstrate that disruption of Nrk caused developmental defects in the acquisition of an adrenergic phenotype in adrenal chromaffin cells. To examine whether the reduction of adrenergic chromaffin cells was associated with a defect in cell proliferation, the proportion of adrenergic (TH+/PNMT+) chromaffin cells, noradrenergic (TH+/PNMT−) chromaffin cells and sympathetic neuroblasts that were cycling (Ki67 immunoreactive) was quantified. Analysis of the growth fraction (proportion of cycling to non-cycling cells) showed loss of Nrk caused significant changes in proliferative behavior in all cell types ( Fig. 6K; Chi-squared test, adrenergic chromaffin cells X 2 , 20.6, p = 0.05, df = 2, N = 4680; noradrenergic chromaffin cells, X 2 , 58.96, p = 0.05, df = 2, N = 4212; sympathetic neuroblasts, X 2 , 17.27, p = 0.05, df = 2, N = 3103, in each case critical value = 5.99). As before, each 2 × 3 contingency table was broken down to 3, 2 × 2 tables and re-tested to determine which pairwise comparisons were significantly different. In all three cell types, Nrk-null were significantly different from wild-type and from Nrk-het. Nrk-hets were also significantly different from wild type animals for PNMT− chromaffin cells but not for PNMT+ chromaffin cells or from sympathetic neurons. Details of the statistical analyses are shown in Table S1.
Loss of Nrk in null mutant mice led to a significant (40-50%) decrease in growth fraction so that only about 7.4% and 13.1% of adrenergic and noradrenergic chromaffin cells respectively were dividing, compared to 11.7% and 23.4% respectively in wild-type mice (Fig. 6K). As the proliferation of both adrenergic and noradrenergic www.nature.com/scientificreports www.nature.com/scientificreports/ chromaffin cells changed to a similar extent, the reduction in the proportion of chromaffin cells with an adrenergic phenotype is unlikely to be due solely, or even mainly, to defects in proliferation.
In contrast to the situation in developing adrenal chromaffin cells, loss of Nrk in sympathetic neuroblasts resulted in a significantly increased growth faction from 26.3% to 32.5% (Fig. 6K). As was the case for the proportion of PNMT+ cells, Nrk-het mice showed effects on proliferation intermediate between wild-type mice and Nrk-null mice in all cell types (Fig. 6K). The changes in the proportion of proliferating cells for neuroblasts in Nrk-null mice occurred despite this cell type failing to express Nrk. Imprinted genes. The Imprinting Resource at Mousebook (http://www.mousebook.org/mousebook-catalogs/ imprinting-resource) lists 151 imprinted genes in the mouse. Of these, our RNA sequencing survey detected 83, largely somatic genes with a few large non-coding RNA genes. Many of the genes not detected were microRNA or small nucleolar RNA, both of which would not have been retained by our RNA purification and detection methods. Of the 83 detected genes, 50 were up-regulated more than 2-fold in adrenal chromaffin cells. Only 9 genes were up-regulated more than 2-fold in sympathetic neuroblasts (Table 4). While more genes overall were upregulated in adrenal chromaffin cells than in sympathetic neuroblasts, imprinted genes represented 1.24% of all genes upregulated more than two fold in adrenal chromaffin cells but only 0.38% of genes upregulated more than two fold in sympathetic neuroblasts.

Discussion
In this study, we describe a technique for separating adrenal chromaffin and sympathetic neuron precursors in embryonic mice before they segregate anatomically, which allowed us to generate a resource of RNA sequencing data on genes differentially expressed by the two cell types, identify a role for Nrk in adrenal chromaffin cell development and produce evidence that many imprinted genes are upregulated during adrenal chromaffin cell development.

Isolating live adrenal chromaffin cells and sympathetic neuroblasts from as early as E12.5.
The separation of developing adrenal chromaffin cells from sympathetic neuroblasts in TH-Cre::R26R-EYFP mice by FACS depended on the differing intensity of the EYFP reporter in the two cell types. Consistent with our earlier study 17 , we showed that, compared to adrenal chromaffin cells, sympathetic neuroblasts show relatively low TH immunoreactivity; unexpectedly, we also showed that in TH-Cre::R26R-EYFP mice, sympathetic neuroblasts show relatively high intensity of the EYFP reporter compared to chromaffin cells. Why the inverse levels of TH immunostaining and EYFP intensity occur in the two cell types is unclear. It could be due to differences in the activity of the promoter in the endogenous TH gene versus the same promoter after recombination in the EYFP transgene. It is also possible that the handling of the foreign protein, EYFP, differs between the two cells or that the intrinsic cytoplasmic microenvironment of the two cell types differs sufficiently to alter the fluorescence of the EYFP reporter. Finally, the timing of the activation of the transgene may be different in the two cells and levels of EYFP may not yet be maximal in the cell type activating the gene last. It should be noted that the difference in intensity level of EYFP by the two cell types only lasts until E14.5, after which they are the same, which favors the  www.nature.com/scientificreports www.nature.com/scientificreports/ latter explanation above. Whatever the explanation for the difference, it is sufficient to separate efficiently the two cell types at E12.5 and E13.5, so that the separated cells carry the correct molecular signature as previously determined by single cell PCR 4 . One advantage of isolating the cells using FACs based on expression of EYFP driven from the TH promoter is that this approach yields living embryonic cells that can cultured for further study.

NRK gene.
Our study showed that Nrk is transiently and highly expressed (903 CPM) in the developing mouse adrenal chromaffin cells on E12.5, but is expressed at only low levels (16.5 CPM) in sympathetic neuroblasts. Nrk is an X chromosome-linked gene that encodes for a Ser/Thr kinase 35 . The protein, NRK, also known as NESK, belongs to the Group I germinal center kinase (GCK) subfamily 32 . Nrk mRNA is expressed by the myotome during embryogenesis as well as in the spongiotrophoblast layer of the placenta, but is not expressed in any adult tissues with exception of the mammary gland of pregnant female mice 33,35,36 .
Like other group I GCKs, NRK appears to activate JNK signaling in model in vitro systems 32 . In this process it sits downstream of TNF receptor associated factor 2 (TRAF2) and upstream of MEKK1 and MKK4 and mediates the effect of TNF alpha 32 . Our data show that Traf2, Mekk1 (as Map3k1) and Mkk4 (as Map2k4) are all expressed at moderate levels in both adrenal chromaffin cells and sympathetic neuroblasts. The roles of NRK in vivo are poorly understood. A previous study in Nrk mutant mouse found that NRK is important for placental development and labor induction, while absence of NRK led to overgrowth of the spongiotrophoblast in the placenta 33 . It was suggested that the phenotype of spongiotrophoblast overgrowth was similar to that seen when another X-linked homeobox gene, Esx1, was knocked out, raising the possibility that NRK can act in the same pathway as Esx1 33 , but Esx1 is not expressed in either adrenal chromaffin cells or sympathetic neuroblasts.
In our study, we found Nrk-deficient mice displayed a significant reduction in the proportion of adrenergic adrenal chromaffin cells in neonates. Thus, one possible action of NRK is to promote the differentiation of adrenergic adrenal chromaffin cells. On the limited evidence available, the differentiation and phenotype of noradrenergic chromaffin cells seemed unaffected by loss of NRK, except there were both less PNMT+ and more PNMT− chromaffin cells per section in the Nrk-null animals, which is consistent with the idea that PNMT− (noradrenergic) chromaffin cells were failing to adopt a PNMT+ (adrenergic) phenotype in the absence of Nrk. Furthermore, disruption of Nrk also lead to an increase in proliferation of sympathetic neurons but a decrease in proliferation of both adrenergic and noradrenergic chromaffin cells in E18.75 Nrk-null mice.
One puzzle is how loss of NRK exerts an action on sympathetic neuroblast proliferation. Levels of mRNA for Nrk were much higher in adrenal chromaffin cells (903 CPM) at E12.5 than in sympathetic neuroblasts (16 CPM) at the same age. Furthermore, expression of Nrk was only elevated in adrenal chromaffin cells at E12.5 and E13.5 and was at low levels in neuroblasts from E11.5 to P0. It is possible that the effect observed at E18.75 in adrenal chromaffin cells could have been exerted by the absence of Nrk in adrenal chromaffin cells at E12.5 and E13.5.  Table 3. Potential regulators of differentiation in sympathetic neuroblasts. Fold changes and mean expression were obtained using the Bioconductor edgeR package 58 with variance between samples by trended estimate. www.nature.com/scientificreports www.nature.com/scientificreports/ Imprinted genes. In our study we showed that 50 imprinted genes were up-regulated more than 2-fold in adrenal chromaffin cells, but only 9 imprinted genes were up-regulated more than 2-fold in sympathetic neuroblasts. Imprinted genes have not previously been specifically associated with the developmental of neural crest-derived cells. Imprinted genes have been suggested to co-express as a gene network in a range of circumstances, including embryonic development, postnatal growth and stem cell regulation [37][38][39] , although there is no evidence of a dominant function(s) for imprinted genes 37 . Varrault et al. 38 suggested that around 15 imprinted genes were commonly expressed together to regulate embryonic growth. Subsequent studies identified many of the same genes as also being involved in regulating postnatal growth and in the regulation of adult stem cells [39][40][41] . A largely similar, but not identical list of imprinted genes comprised the network in each of these studies [38][39][40][41] and overlapping genes included H19, Igf2, Ndn, Peg3, Cdkn1c, Mest, Dlk1, Gtl2 and Grb10. All are also highly differentially expressed in developing chromaffin cells in our data, suggesting that the same imprinted gene network is active. Recently, the putative imprinted gene network has been expanded to include a large proportion of known imprinted genes and it has been suggested that many non-imprinted somatic genes are also co-activated with the imprinted gene network 37 .

Gene symbol Gene name
Among imprinted genes are a number of non-coding RNAs 42,43 including the long non-coding RNA, Meg3 (maternally expressed gene 3), and many microRNAs. Loss of Meg3 and the microRNAs it controls leads to reduced neural differentiation in human embryonic stem cells 43 . Recently, microRNAs have been implicated in adrenal chromaffin cell fate determination 24 , as deletion of the RNAse, Dicer, which processes miRNAs, leads to a reversion of developing chromaffin cells to a more neuronal phenotype 24 .
Many of the same imprinted genes are also highly expressed in adrenal chromaffin cells in the raw data of Furlan et al. 4 . Analysis of the levels of imprinted gene expression across the four stages of adrenal chromaffin cell differentiation (Schwann cell precursor, early bridge cell, late bridge cell, adrenal chromaffin cell) in that data set suggests that, for most imprinted genes, expression is at very low or modest levels in Schwann cell precursors and at successively higher levels as differentiation proceeds (Fig. 7). A few imprinted genes were exceptions to this pattern; for instance, Cdkn1c was highest in Schwann cell precursors and decreased thereafter.
Overall, the importance of imprinted genes in development of adrenal chromaffin cells remains unclear. The silencing occurs by DNA-methylation and/or histone modification 44,45 and is important for maintaining normal embryonic development and growth 46 . Deregulation of genomic imprinting gives rise to developmental disorders such as Beckwith-Wiedemann syndrome as well as tumorigenesis [47][48][49] . Most imprinted genes exist in chromosomal clusters called imprinted loci where the genes in each locus share an imprinting control region within the differentially methylated region 50,51 . The expression of imprinted genes may indicate that the development of www.nature.com/scientificreports www.nature.com/scientificreports/ adrenal chromaffin cells is a site of maternal vs paternal genome resource conflict, with different evolutionary pressures being exerted on the male and female genomes in this stage of development. Alternatively these genes are known to play developmental roles in other tissues and stages and imprinting may be the consequence of the role of these genes in the other tissues and limit proliferation. Expression of the network of imprinted genes often appears to correlate with withdrawal from the cell cycle and differentiation, perhaps to support subsequent tissue growth. More specific roles in embryonic development cannot be ruled out. One important differentially methylated region (DMR) includes the imprinted genes Dlk1 and Gtl2 (Meg3) and changes in expression or imprinting at this DMR is associated with a range of cancers, including neuroblastoma 52 . Changes in the expression of the imprinted microRNAs mir134a and mir335 have also been reported in many cases of neuroblastoma 5 . Given the upregulation of so many imprinted genes in developing adrenal chromaffin cells and the potential role of epigenetic dysregulation in neuroblastoma, imprinted genes represent novel targets for study.

Animals. All animal experiments were approved by the University of Melbourne Animal Experimentation
Ethics committee and complied with National Health and Medical Research Council of Australia (NHMRC) guidelines. TH-IRES-Cre;ROSA26-EYFP (TH-Cre::R26R-EYFP) mouse embryos were generated by mating heterozygous TH-IRES-Cre mice with homozygous ROSA26-EYFP mice. The TH-IRES-Cre (nomenclature: B6.129X1-Th tm1(cre)Te /Kieg) mice were a kind gift from Prof A. Allen (University of Melbourne) that were originally obtained from Dr. T. Ebendal's laboratory via the European Mutant Mouse Archive repository and maintained on a C57BL/6 background 53 . ROSA26-EYFP mice (nomenclature: B6.129X1-Gt(ROSA)26Sor tm1(EYFP)Cos /J) were obtained from Jackson Laboratories and were also maintained on a C57BL/6 background. Nrk mutant mice (nomenclature: B6.129S4-Nrk tm1Mkom ) were generated by crossing Nrk-null males with Nrk-het females and maintained on a C57BL/6 background. Wild-type (WT) C57BL/6 mice were used for validation experiments. All mice were time plug-mated and the morning of detection of the plug was counted as embryonic day (E) 0.5. Pregnant dams at 11.5-14.5 days post-fertilization were killed by cervical dislocation and the embryos were collected and decapitated. See Supplementary data for genotyping and phenotyping methods.
Tissue Dissection. Embryos were dissected in ice-cold "dissecting medium" (DMEM/F-12, Gibco, cat.  www.nature.com/scientificreports www.nature.com/scientificreports/ urogenital ridges and the tissues in between, including the dorsal aorta were then detached from the dorsal body wall by inserting the tips of fine forceps underneath, and pulled away by holding the dorsal aorta. In E11.5 tissues, the urogenital ridges were removed from the tissue of interest. In E12.5, the tissue of interested was further sub-dissected by removing the urogenital ridges, the kidneys and the paravertebral sympathetic chains, if attached. For E13.5 and E14.5 embryos, the urogenital tract was dissected according to Buehler et al. 54 , with modifications. Briefly, the viscera from the liver down to the genital tubercle were pulled away from the dorsal body wall by inserting fine forceps underneath the dorsal aorta at the level of the heart. In E13.5 tissues, this tissue mass was further microdissected to retain just the region containing the adrenal anlagen and the adjacent tissue surrounding the dorsal aorta. By E14.5, the adrenal glands and sympathetic ganglia were present with clearer anatomical segregation. Therefore, in E14.5-P0 tissues, the adrenal glands and the surrounding sympathetic ganglia were separately sub-dissected into different samples. Comparable tissues were also dissected from EYFP negative (−) embryos for compensation controls in the FACS. Tissues from 5-7 embryos were cut into smaller pieces and pooled into a 15 mL centrifuge tube (Corning) with 4 mL of ice-cold dissecting medium.
Cell Dissociation and FACS Isolation. The tissues were then dissociated in 200 µL embryo −1 Accumax (Innovative Cell Technologies) followed by incubation in a 37 °C incubator for 45 min. 1 mL of quench medium (dissecting medium supplemented with 10% FBS, Life Technologies; 100 units mL −1 penicillin and 100 µg mL −1 streptomycin, Gibco; 2 mM Glutamax, Thermo Fisher; 37.5 µg mL −1 DNase, Sigma-Aldrich) was added to each tube to stop the enzymatic activity and the tissues were triturated gently until completely dissociated. The dissociated cells were filtered through a 70 µm nylon Cell Strainer (Falcon). The tube and filter was then washed with 1 mL of quench 1:5 medium (quench medium supplemented with 7.5 µg mL −1 DNase) to release any remaining cells and the cells in the filtrate were pelleted by centrifugation at 250 g, 4 °C for 5 min. The pellet was re-suspended in 500 µL quench 1:5 medium and viability dye, 7-Aminoactinomycin D (7-AAD, BD Pharmingen) was added at 1:1000 (v/v) dilutions to label dead cell for FACS. Compensation controls were also prepared with 200 µL EYFP− cell suspension (double negative), 200 µL EYFP− cell suspension labelled with 0.2 µL 7-AAD (7-AAD only) and 200 µL 1:50 (v/v) dilution EYFP+ cell suspension (EYFP only).
Cells were filtered through a 35 µm cell strainer cap again into a 5 mL polystyrene round bottom tube (Falcon) prior to FACS. Flow cytometry analysis and cell sorting were performed on a BD Influx or BD FACSAria Fusion cell sorter with a 100 µm nozzle at 20 psi (BD Biosciences). The compensation controls were evaluated to optimize the voltages and gating for sorting. Live EYFP+ single cells were first isolated from dead cells on the basis of 7-AAD (excitation: 488 nm, emission 692/40 nm) and EYFP (excitation: 488 nm, emission 530/40 nm) intensity. The populations from E12.3 and E13.5 samples were then further analysed on the basis of side-scattered light (SSC) and EYFP intensity and sorted into EYFP+ Lo and EYFP+ Hi populations with the elimination of cells lying in between. For the live EYFP+ cells from E11.5 and the separated adrenal gland and sympathetic ganglia samples in E14.5-P0, a single homogenous population was observed so that no further FACS gating was needed. Dot-plots were generated by BD FACS TM Sortware software. The cells were sorted directly into 750 µL ice-cold TRIzol ® LS reagent (Ambion) in a 1.5 mL nuclease-free LoBind microcentrifuge tube (Eppendorf) for RNA extraction.

RNA Extraction and Droplet Digital PCR.
Total RNA from isolated EYFP+ Hi and EYFP+ Lo cells were extracted using by TRIzol/chloroform extraction followed by RNeasy Micro Kit (Qiagen). RNA pellets with ~20,000 cells from more than 10 embryos with a sex ratio within a range of 1:1.3 (male:female or female:male) www.nature.com/scientificreports www.nature.com/scientificreports/ were pooled. Total RNA quality and quantity from each sample were analysed on a Bioanalyzer 2100 (Agilent) and the RIN numbers of all samples were >8 with average RIN = 9.7. For Droplet digital PCR, cDNA was synthesized using the iScript Advanced cDNA Synthesis Kit (Bio-Rad Technologies). Droplet digital PCR was performed in biological triplicates and technical duplicates on a QX100 Droplet Digital PCR system (Bio-Rad Technologies) with TaqMan Assays primer/probe mixture (Thermo Fisher Scientific). See Supplementary data for detail method. RNA Sequencing. Four biological replicates of each EYFP+ Hi and EYFP+ Lo paired RNA sample were generated for transcriptomic analysis. cDNA library was prepared by random priming using the SMART-seq v4 Ultra Low Input RNA Kit (Clontech, Cat. 634889) followed by PCR amplification with the addition of Illumina adaptors for multiplexing experiment using Nextera XT Kit (Illumina, Cat. FC-131-1096). Sequencing was performed with 50 base pair single-end reads on an Illumina HiSeq 2500 platform.
Bioinformatics Analysis. Samples from RNA-seq were demultiplexed and raw reads were generated in fastq format by HiSeq Analysis Viewer software (Illumina). Around 29 million reads were detected per sample. Quality assurance and quality control for all samples were evaluated using FastQC analysis. The effect size of biological versus technical variances were analyzed by principal component analysis. Reads were aligned with HISAT2 55 to the Genome Reference Consortium mouse genome, GRCm38 with gencode consortium gene annotation M8 56 . Reads were assigned to genes using HTSeq 57 . Differential expression between adrenal chromaffin precursor cells and sympathetic neuroblasts was obtained using the Bioconductor edgeR package 58 . Fold changes were calculated by the gene expression level in adrenal chromaffin cells divided by the gene expression level in sympathetic neuroblasts. Genes with fold change greater than 2 and p-value of less than 0.05 after Benjamini-Hochberg false discovery rate correction were considered significantly differentially expressed. Gene ontology analysis for genes higher in chromaffin precursor cells and sympathetic neuroblasts of the differential transcriptomes was performed by PANTHER v.11 59 statistical overrepresentation test with Mus musculus reference list. Ingenuity Pathways Analysis v.01-06 tool was employed to the differential expressed genes in chromaffin cells to sympathetic neuroblasts for gene set and signaling pathways enrichment analysis.
Immunohistochemistry and Data Analysis. Immunohistochemistry followed the methods described in 17 .
See Supplementary data for antisera used. For comparison of TH, CART and EYFP immunoreactivity, confocal images were taken on a Zeiss Meta 501 scanning confocal microscope. For comparing TH, Ki67 and PNMT immunoreactivity in Nrk mutant mice, confocal images were taken on a Zeiss LSM 800 confocal microscope. All images were processed by Zeiss Image Browser (v4.0.0241, Carl Zeiss Microimaging). Graphs were prepared by using Numbers software (version 3.6.2, Apple Inc.) or Prism software (version 7.0a, GraphPad). Statistically analysis was performed with Prism 7 software using two-way ANOVA for gene expression pattern, linear regression for correlation analysis, and chi-squared test for Nrk loss-of function analysis.