Neurotranscriptomics: The Effects of Neonatal Stimulus Deprivation on the Rat Pineal Transcriptome.

The term neurotranscriptomics is used here to describe genome-wide analysis of neural control of transcriptomes. In this report, next-generation RNA sequencing was using to analyze the effects of neonatal (5-days-of-age) surgical stimulus deprivation on the adult rat pineal transcriptome. In intact animals, more than 3000 coding genes were found to exhibit differential expression (adjusted-p < 0.001) on a night/day basis in the pineal gland (70% of these increased at night, 376 genes changed more than 4-fold in either direction). Of these, more than two thousand genes were not previously known to be differentially expressed on a night/day basis. The night/day changes in expression were almost completely eliminated by neonatal removal (SCGX) or decentralization (DCN) of the superior cervical ganglia (SCG), which innervate the pineal gland. Other than the loss of rhythmic variation, surgical stimulus deprivation had little impact on the abundance of most genes; of particular interest, expression levels of the melatonin-synthesis-related genes Tph1, Gch1, and Asmt displayed little change (less than 35%) following DCN or SCGX. However, strong and consistent changes were observed in the expression of a small number of genes including the gene encoding Serpina1, a secreted protease inhibitor that might influence extracellular architecture. Many of the genes that exhibited night/day differential expression in intact animals also exhibited similar changes following in vitro treatment with norepinephrine, a superior cervical ganglia transmitter, or with an analog of cyclic AMP, a norepinephrine second messenger in this tissue. These findings are of significance in that they establish that the pineal-defining transcriptome is established prior to the neonatal period. Further, this work expands our knowledge of the biological process under neural control in this tissue and underlines the value of RNA sequencing in revealing how neurotransmission influences cell biology.

Introduction

The rodent pineal gland is an especially attractive model for neurotranscriptomic studies because of evidence of large neurally-regulated changes in the abundance of many pineal transcripts [1-9]. The gland is innervated by sympathetic projections from the superior cervical ganglia (SCG) [10-13]. At night norepinephrine is released in response to signals originating in the hypothalamic suprachiasmatic nuclei (SCN), which house an autonomous circadian clock. The SCN are hard wired to the pineal gland by a multisynaptic neural pathway: SCN neurons project to the paraventricular nucleus, where they contact neurons that project caudally via the mesencephalic periaqueductal gray to the intermediolateral nuclei in the upper thoracic segments of the spinal cord [14-16]. From there, preganglionic neurons contact a small subpopulation of SCG cells which innervate the pineal gland via projections through the internal carotid nerve and the conarian nerves. Light acts on this system through the retina and a retinohypothalamic projection that terminates in the SCN; light entrains the SCN clock to environmental lighting and also gates stimulatory signals to the pineal gland [14, 17–20]. Experimental deprivation of neural stimulation of the pineal gland can be accomplished by removal (SCGX) or decentralization (DCN) of the SCG, thereby eliminating circadian input from the SCN [21].

Norepinephrine is the primary biogenic amine controlling the pineal gland and acts through an α1b- and β1-adrenergic receptor “AND” gate mechanism to activate adenylate cyclase and elevate cyclic AMP [22-28]. Whereas β1-adrenergic receptor activation is essential for activation of adenylate cyclase, this effect is potentiated by α1b –adrenergic receptors via Ca++ and phosphatidyl inositol (Pi) activation of protein kinase C. The adrenergic, receptor-dependent elevation of cyclic AMP activates protein kinase A, which in turn phosphorylates cyclic AMP response element binding protein, which is bound to regulatory elements, thereby triggering gene expression [27, 29–31]. In addition to releasing norepinephrine, adrenergic neurons serve to take up and sequester norepinephrine and other catecholamines; this prevents stimulation by circulating catecholamines entering the pineal perivascular space, thereby enhancing the on/off nature of neurotransmission [32]. Accordingly, whereas both SCGX and DCN block stimulation of the pineal gland, the reuptake/sequestering function of the sympathetic nerves is retained following DCN and is eliminated by SCGX [32]. In addition to neural regulation of the pineal transcriptome, peripheral clock mechanisms operating within the tissue may play a regulatory role.

Previous neurotranscriptomic studies of the pineal gland have used a variety of classical biochemical methods and cDNA microarray technologies [2, 5, 9, 33, 34]. In the current report, strand-specific next-generation RNA-Seq is used to extend this body of knowledge. This powerful technology provides unprecedented statistical power across a broad dynamic range, allowing robust detection of differential effects in more genes than ever before [35-38]. This effort has increased our knowledge of the rodent pineal gland by identifying genes not previously known to be under neural control. Examination of the effects of neonatal (5-days-of-age) SCGX and DCN on the adult rat pineal transcriptome revealed that these procedures essentially eliminate the 24-hour pattern of differential gene expression but do not markedly alter the developmental expression of the non-rhythmic component of the pineal-defining transcriptome.

Our results provide a powerful new resource for investigators studying profiles of rhythmic and non-rhythmic expression in the pineal gland.

Results

In vivo night/day differential expression of pineal transcripts, with and without neurotransmission

More than 3000 genes showed statistically significant night/day differential expression in the Control and Sham animals at an adjusted-p-value cutoff of 0.001, with more than 1400 genes exhibiting greater than 2-fold changes in either direction (Fig 1a and 1b, Table 1). Of these, 257 genes exhibited strong differential expression in the Control group with fold change (FC) greater than 4 in either direction, 195 genes met these criteria in the Sham-operated group, and 157 genes exhibited met these criteria in both groups (Fig 2a). The RNA-Seq results replicated 255 out of the 343 genes found in previous microarray-based studies to exhibit greater than 2-fold night/day differential expression (in either direction) at the adjusted-p < 0.05 level (S1 Fig, S1 Dataset) [9]. In addition, more than a thousand genes were identified that had not been previously found to exhibit differential expression between night and day at this level.

Fig 1

MA plots for six differential expression analyses.

These plots display the mean normalized read-pair counts (x-axis) versus the estimated fold change (y-axis), on a log-log scale for four in vivo analyses and two in vitro analyses. The blue shading indicates the density of genes, and each red point represents a gene with statistically significant differential expression (adjusted-p < 0.001). Dashed horizontal lines mark 4-fold changes in both directions, dashed vertical line indicates minimum abundance threshold for the statistical tests. The four in vivo analyses compared night and day time points in adult rats for the following groups: (a) no surgery (Control); (b) neonatal sham surgery (Sham); (c) neonatal superior cervical ganglia decentralization (DCN); (d) neonatal superior cervical ganglionectomy (SCGX). The two in vitro analyses compared treated/untreated pineal glands: (e) norepinephrine-treated (NE) vs untreated and (f) dibutyryl-cyclic-AMP-treated (DBcAMP) vs untreated.

Table 1

Genes significantly differentially expressed in the pineal gland on a night/day basis in Control and/or Sham groups, in vivo (adjusted-p < 0.001).

Genes that exhibit strong differential expression in both are listed in bold (genes are classified by the least-extreme fold change across the two analyses). Note that some genes are listed twice if the fold changes of the two analyses fall in different intervals. In the SCGX and DCN groups, none of the genes displayed in this table exhibited night/day differential expression at this level (adjusted p < 0.001, FC > 4 in either direction). A complete list of all genes with fold changes, p-values, and normalized expression estimates is available in the SI (S1 Dataset).

Fold Change	Gene Symbol
>32	Aanat * † , Atp7b * † , Dclk3* † , Dusp1 † , Fcer1a* † , Hcrtr1 † , Hhip † , Irs2* † , Nap1l5 † , Ribc2* † , Sik1 * † , Slc15a1 †
16–32	Cd24 † , CREM * † , Dclk3 * † , Dos † , Drd4 * † , Dusp1 † , Fcer1a * † , Hcrtr1 † , Irs2 * † , LOC100362205, MGC94891*, Nap1l5 † , nod3l † , Pde10a * † , PORF1*, Ptch1 * † , Spata2L, St8sia2*, Wnt10a †
8–16	Accn2, B3gnt8 *, Bhlha15, Brunol6, Camk1g* † , Ccrl2 † , Cd24 † , Chd5, Coq10b * † , Cyp2s1, D3ZJR7 † , Dclk1 † , Dio2 * † , Dos † , Elfn1, Errfi1 † , Etnk1* † , F1LN23, Fam160b1* † , Gadd45b, Galntl1* † , Gem * † , Gls2 † , Grm2, Il20ra, Ipcef1, Irak2 * † , Kcnq4 *, Kctd3* † , Lmo1, LOC691317*, Nr4a1 † , Nr4a3 † , Osbpl6, Padi4* † , Pde4b* † , Pde8b † , Per2 † , Pnma3, RGD1309903, RGD1310852*, RGD1560523* † , RGD1561499, RGD1565158, Rgs2 † , Rxfp2 *, Slc17a6 † , Spata2L, St8sia2 *, St8sia5 † , Syt4 † , Tinagl1 † , Zrsr1 * †
4–8	7SK*, Actn2, Acy3, Ankrd52 † , ARLY † , Arpp21 † , Asgr1, B3gnt2 † , Bex1 † , Bhlha15, Bsx *, Cacna1h* † , Calhm1 *, Calhm2* † , Camk1g * † , Ccrl2 † , Cd8a † , Ceacam1, Ceacam10, Clstn3, Csrnp1, D3ZAQ3, D3ZDX5 † , D3ZJR7 † , D3ZRS0, D3ZW88* † , D4A582, D4A5F4, D4A8J6 *, D4ACZ5 † , D4AE14 † , Dnm2 † , Drd1a, Egr1 † , Elfn1, Elfn2 † , Ell2 † , EMAL5 *, Emd † , Ephx4, Errfi1 † , Etnk1 * † , F1LNL3 * † , F1LNM2 † , F1LU85 † , F1LWE5, F1LXI5, F1LY53, F1M6I6, F1M848*, F1M8Z8, Fam160b1 * † , Fam161a*, Farp2* † , Fdx1 † , Fgd2, Fkbp5* † , Fry † , Gadd45b, Gal, Galntl1 * † , Gdf15* † , Glt8d3 * † , Gnaz † , Gpr158, Greb1l *, Grm1 † , Hspa1a † , Hspa5 † , Igf1r † , Il20ra, Ipcef1, Irf8 † , Itih1, Kcnv2, Kctd3 * † , Krt19 † , Lamb3 † , Lhx2, Lmo1, LOC100361460*, LOC100364604, LOC689856, LOC691141, LOC691317 *, Lrp8, Lrrc43, Lypd6b, Mal2, Man2a1 † , Mapk10 † , Mapk4, Mbnl2 * † , Mcam † , Mettl8 *, Mfsd2, MGC94891 *, Mmd † , Mmd2 † , Mt1a † , Myh13*, Ncald † , Neurog1*, Nos1, Nptx1 * † , Osbpl6, Padi1, Padi4 * † , Parvb † , Pcdh1 † , Pde4b * † , Pde7b, Pdzrn3 † , Per1, Plagl1, Plcxd1, Plk5* † , Pmepa1, Pnma3, PORF1 *, Ptp4a1 † , PVR † , Q7TP44 † , Rap1gap, Rcan1 * † , Rem2, RGD1309108* † , RGD1309489, RGD1309586*, RGD1309903, RGD1310852 *, RGD1560523 * † , RGD1566021 *, Rgs2 † , Rlim, Rnase1*, Rorb*, Scnn1b, Shc3, Sla2, Slc17a6 † , Slc30a2 † , Smc1b, Spag4, Tbc1d1 † , Tbc1d30 † , Tbx3, Thap4 † , Thbs4, Thrb, Tmem108 † , Tmem225, Tpcn1 † , Trank1, Trnp1, Tyro3 † , Wfikkn2*, Xpot †
1/4–1/8	Adra1a, Akr1b8, Angptl4 † , Atp10a, B4F7B7, Barhl2, Bpil3, Cacna1e, Car8, Ccdc85a, Ccl9 † , Chrna5, Cpt1a † , Crp, Crybg3, Cstl1*, D3ZV53 *, D3ZW53 *, D3ZX30, D3ZXI1, D3ZYL1, D4A0N6, D4A1L0, Dhrs9, Dscaml1, F1LUT6, F1LVH8, F1LW34, F1M0S2, F1M1I4 *, F1M7W8, Frmd4b, Fscn2, Galnt12, Galnt14, Gata1*, Gng11, Gramd2, Grip2 † , Gucy1a3 † , Hook1, Hs3st2 * † , Htra4, IGG2B, Il17re, Lad1-ps1, LOC56764, LOC690918 *, Ms4a2, Ms4a3, Ndrg1 † , Npy1r † , Opn1sw* † , Opn3, Prc1 † , Pter, RGD1308116 † , RGD1308133, RGD1310773, RGD1559884*, RGD1560925, RGD1563092, RGD1564942, Rgs20, Rnase1l2*, Slc37a2, Slc39a12, Slc6a13, Slco4c1, Sstr3 † , St8sia3 † , Sult1c3, Tfpi2 † , Ucp3, Wdr89 † , Zfp385d
1/8–1/16	Ccl9 † , D3ZW53*, F1M7W8, Fbp2, Fscn2, Hs3st2* † , Htra4, HYES, LOC100362965, LOC100365668, LOC690918*, Npy1r † , PTN21 † , RGD1304931, Rptn *
1/16–1/32	Fam111a, RGD1304963 †
<1/32	RGD1304963 †

*Gene exhibits high relative expression in the pineal gland (see S2 Table).

†Gene has previously been found to be differentially expressed between day and night conditions in microarray experiments.

Fig 2

Venn diagrams displaying the overlap between the different differential expression analyses.

Diagrams indicate the concordance in the genes found to display strong (FC > 4 in either direction) and statistically significant (adjusted-p < 0.001) differential expression in multiple separate analyses. (a) Comparison of the (in vivo) Control and Sham night/day differential expression analyses. (b) Comparison of the (in vitro) NE vs untreated and DBcAMP vs untreated differential expression analyses. (c) Comparison of the overlapping region from (a) with the results of the two in vitro analyses from (b).

Stimulus deprivation resulting from neonatal SCGX or DCN eliminated nearly all night/day differential changes observed in the Control- and Sham-operated groups (see Fig 1c and 1d). Principal component analysis reveals that the SCGX and DCN samples for both day and night cluster tightly with the Sham-day and Control-day samples (S2e Fig). Furthermore, when the SCGX or DCN groups were directly compared with the Sham group, the vast majority of the differences were found at night in rhythmic genes (S2a–S2d Fig). All of these observations are consistent with the hypothesis that the changes resulting from SCGX and DCN are primarily due to the loss of SCG-mediated central neural stimulation.

A small set of genes were identified in which the expression differed significantly (adjusted-p < 0.001), substantially (FC > 2 in either direction), and consistently (FC in the same direction at both day and night) between the DCN and Sham groups (Table 2, S2 Dataset). Of particular note: at both day and night the gene Serpina1 exhibited a >6-fold reduction in expression in the DCN group compared with the Sham group. Weaker but still significant effects were observed in Serpina1 when comparing the SCGX group to the Sham group.

Table 2

Genes with consistent differential expression in the DCN group relative to the Sham group.

These genes displayed significant (adjusted-p < 0.001) differentials across 2 analyses: DCN-day vs Sham-day and DCN-night vs Sham-night. Genes are classified by the least-extreme fold change across both analyses. The complete results for these analyses are available in the SI (S2 Dataset).

Fold Change	Gene Symbol
>8	-
4–8	Rcvrn, LOC681316, RGD1304731
2–4	Gabra1, F1LZ56, RGD1306717, Kcnab1 † , E9PSL1, March4, Cldn7, Rnf182, Ppp1r3c, Tuft1, Kcnj11, ABCC8, Rnf222*, D4A4G4, RGD1563692, RGD1563070, Gata5, LOC688006
1/2–1/4	Fry † ‡ , Lnx1 † , Parm1 § , Wfikkn2* ‡ , F1LVN4, Pde4b* † ‡ § , RGS6, Gngt1, F1LWK7, Scg2 † , D3ZEA8, Ptprn*, Prr16, F1MAG9*, Sult2b1 † , LOC691300 § , Rdh2, RGD1563903*, Dclk3* † ‡ § , CAD23, Gpr123, Glod5, Cpe
1/4–1/8	Sftpa1, Rln3* § , Pde10a* † ‡ § , Il20ra ‡ , Gda, Tyrp1, Serpina1, Ribc2* † ‡ §
<1/8	-

*Gene exhibits high relative expression in the pineal gland (see S2 Table).

†Gene has previously been found to be differentially expressed between day and night conditions in microarray experiments (adjusted-p < 0.05).

‡Gene has strong (FC > 4 or FC < 0.25) and statistically significant (adjusted-p < 0.001) differential expression in the Control and/or Sham night/day in vivo analyses (Table 1).

§Gene has strong and statistically significant differential expression in either or both of the treated vs untreated in vitro analyses (Table 3).

Table 3

Genes significantly differentially expressed in pineal gland tissues in response to in vitro treatment with NE or DBcAMP (adjusted-p < 0.001).

Genes that exhibit strong differential expression in both the NE and DBcAMP analyses are listed in bold (genes are classified by the least-extreme fold change across the two analyses). A complete list of all genes with fold changes, p-values, and normalized expression estimates is available in the SI (S3 Dataset).

Fold Change	Gene Symbol
>32	Aanat * † , RGD1305627, Slc15a1 †
16–32	Dclk3* † , Dusp1 † , Fcer1a* † , Il11, Nr4a1 † , Rgs2 †
8–16	Actn2, Atp7b* † , Bhlha15, Calca, Ccrn4l, Cd24 † , Chd5, Chst2 † , CREM* † , D4A599, Dclk3 * † , Dct, F1LWE5, Fam148b, Fcer1a * † , Giot1, Gpd1, Grm2, Has1, Hcrtr1 † , Il6, Ipcef1, Irs2 * † , Kcnk10, Kcnq4*, LOC100360880, LOC100361460*, LOC691300, Mstn, Mustn1, Nap1l5 † , Nefh, Nlrp10, Nmu, Nptx1 * † , Nr4a1 † , Nr4a3 † , Pde10a * † , Pdp2 † , Plagl1, Ptch1 * † , Rgs2 † , Ribc2 * † , Rln3*, Rrad, Sik1 * † , Slc17a6 † , Sstr2, Tinagl1 † , Vgf
4–8	Adamts15, Ankrd13b, Atp7b * † , B3gnt2 † , B3gnt8 *, Bhlha15, Bsx*, Camk1g * † , Cbln1, Ccrn4l, Chst2 † , Cited1*, Clic6, Clstn3, Coq10b * † , Cpd, CREM * † , Csf3, Cxcl3, Cxcl6, D3ZD52, D3ZHM0, D3ZVN6, D3ZW88 * † , D4AA43, Dos † , Drd4* † , Errfi1 † , Etnk1 * † , F1LU85†, F1LXI5, F1LZR9, F1M1S6, F1M578, F1M8Z8, Fam107b, Fam160b1* † , Fam46b, Fam71e2, Fdx1 † , Fgd2, Fos, Gadd45b, Gchfr, Gem * † , Ggcx *, Giot1, Gjc2, Glt8d3* † , Gpd1, Gpr83, Gprc5a, Greb1l*, Grm2, Hcrtr1 † , Herc4, Hhip † , Id1 † , Il1b, Il1rn, Il22, Inhba, Irak2 * † , Kcnq4 *, Klf4, Lamb3 † , LOC100360507, LOC100360880, LOC100361460 *, LOC681749, LOC691300, Lrrn3 † , Lxn † , Mal2, Map3k8, Masp1, Mfsd2, MGC94891 *, Mmd2 † , Mstn, Musk*, Myo16, Nap1l5 † , Nmu, Nod2, nod3l † , Nos1, Nos2, Npas4, Nppc, Nr4a2, Nt5dc3, Nt5e, Osbpl6, Padi4 * † , Parm1, Pde4b * † , Pde8b † , Pdxk, Plk3, Prickle1, Prokr2, Prrg4, Ptp4a1 † , Q6MG91, Q7TP44 † , Rasd1, Rassf8, Rcan1* † , Rem2, RGD1560523 * † , RGD1564664, RGD1566021*, RGD1566035, Rhobtb3, Rlim, Rln3 *, Rrad, Rxfp1, Rxfp2*, Scnn1b, Scnn1g, Sgk493, Slc16a6, Slc17a6 † , Slc23a2 † , Slc25a25, Slc30a2 † , Slc30a3, Sostdc1, Spt1, Sstr2, St8sia2 *, Syt13, Syt4 † , Tbc1d1 † , Thbd, Tinagl1 † , Tkt, Top1, Wnt10a † , Xpot † , Yod1, Zrsr1 * †
1/4–1/8	F1M1I4*, Klkb1, Lgi2, LOC690918 *, PTN21 † , RGD1308133, Slitrk6, Wnt2b*
<1/8	-

*Gene exhibits high relative expression in the pineal gland (see S2 Table).

†Gene has previously been found to be differentially expressed between day and night conditions in microarray experiments.

Three genes dedicated to melatonin synthesis: Asmt, Tph1, and Gch1, were not found to exhibit strong changes in expression as a result of SCGX or DCN. In comparisons between the SCGX/DCN groups with the Sham group at day and night, these three appeared to consistently differ by less than 35%. While some of these differentials were moderately statistically significant (adjusted-p < 0.05), the very small effect size implies that neural stimulation is not required to maintain relatively-high levels of expression.

In vitro differential expression of pineal transcripts resulting from treatment with norepinephrine or dibutyryl cyclic AMP

As indicated in the introduction, neural control of pineal transcription involves release of norepinephrine and elevation of the second messenger cyclic AMP. This was extended using pineal organ culture. Glands were incubated for 48 hours under control conditions and then treated with norepinephrine or dibutyryl cyclic AMP. Subsequent RNA-Seq analysis found that 97 genes displayed strong and statistically significant differential expression in the norepinephrine-treated vs untreated analysis, and 170 in the dibutyryl-cyclic-AMP-treated vs untreated analysis (adjusted-p < 0.001, FC > 4 in either direction, see Fig 1e and 1f). A set of 78 genes met these criteria in both analyses (see Table 3, Fig 2b, S3 Dataset).

Interestingly, even though the Serpina1 gene displays no apparent night/day differential expression in vivo, it does display a moderate increase in expression in vitro when treated with NE (adjusted-p = 0.00001, FC = 2.61) or dibutyryl cyclic AMP (adjusted-p = 0.030, FC = 1.73) (S3 Dataset).

Comparison of the in vivo and in vitro experiments revealed strong concordance between the in vivo Control and Sham-operated day/night experiments and both in vitro treated/untreated experiments, with 49 genes exhibiting strong (FC > 4 in either direction) and highly-significant (adjusted-p < 0.001) differential expression in all four analyses (see Fig 2c; S1 Table).

Identification of pineal marker genes

Comparison of the abundance of transcripts in the pineal gland (at both day and night) with other tissues identified a set of 349 strongly-expressed genes (normalized read-pair count > 30) with a maximum-likelihood fold change greater than 8 in the pineal gland at either day or night (S3 Fig, S2 Table, S4 Dataset). Due to the unavailability of biological replicates in this comparison, hypothesis tests could not be performed on these comparisons. Thus, in the absence of other corroborative evidence (such as data from other studies [9]), genes identified as being highly expressed in the pineal gland by this analysis alone are best viewed as potential candidates for further investigation using other experimental approaches.

Discussion

This effort has expanded our knowledge of the impact that neural stimulation has on the pineal transcriptome by identifying regulated genes that are expressed differentially on a 24-hour schedule. Perhaps the most remarkable element is the requirement of sympathetic neural stimulation for essentially all of the observed 24-hour changes in thousands of pineal transcripts. This adds to evidence that the transcriptome profile can be dynamically regulated by neurotransmission, while providing little-to-no support for the view that peripheral mechanisms operating within the rat pineal gland play an important role in the night/day changes in gene expression in this tissue, a topic that deserves further study.

In the case of the pineal gland, it is likely that many of the observed changes occur in pinealocytes, which constitute more than 95% of the cells in the rodent pineal gland [39]. However, it is also possible that some of the observed changes occur in other resident cell types including interstitial cells, endothelial cells, glia cells, or in transient cells. Establishing the precise localization of the observed changes in transcript abundance is important in understanding the functional relevance of the observed changes.

Previous analyses of the daily changes in the pineal transcriptome have made it clear that essentially all aspects of cell physiology are subject to neural control [9]. It would appear that this reflects the dedication of the pineal gland to the daily production of melatonin. Although one of the enzymes in this pathway is under neural control, such control appears to involve processes with no obvious link to melatonin production, including many unrelated functions linked to lipid biology, small molecule transport, extracellular matrix processing, cell-cell contact and other processes.

Transsynaptic control of the pineal transcriptome

It is well-established from previous studies that production of melatonin is neurally controlled via the norepinephrine/cyclic AMP mechanism [9]. In this study we extended the observation that many of the changes in transcript abundance that occur at night are reproduced in vitro by treatment with either norepinephrine or dibutyryl cyclic AMP. This included the transcription factor encoding genes Nr4a1 and Nr4a3 which were found in previous studies to increase at night and in response to norepinephrine [9, 40]. Thus, this adds to evidence that these transcription factors are likely to mediate cyclic AMP control of the pineal transcriptome. All of these observations provide strong support for the hypothesis that neural stimulation via norepinephrine drives the bulk of the observed changes in the pineal transcriptome. However, the results also provide reason to consider that some neurally induced changes seen in vivo on a night/day basis are not simply a reflection of the presence of norepinephrine acting through cyclic AMP. Many genes with observed in vivo night/day differentials did not exhibit strong and statistically significant differential between untreated and treated samples, in vitro. There are several possible hypothetical explanations for this. First of all, the statistical methods used here do not guarantee uniform power between experiments and such apparent discrepancies can occur by chance. Secondly, it is possible that the pineal gland is less responsive in organ culture because in vitro conditions do not optimally reproduce in vivo conditions, and (undetermined) essential factors may be absent from the defined culture medium. This phenomenon has been observed in the pineal gland with the adrenergic stimulation of Drd4 expression, which is known from previous studies to require the addition of thyroid hormone [41]. Similar dependencies may exist for other genes. Thirdly, it is possible that these apparent inconsistencies between night/day and treated/untreated differentials may be due to the involvement of other transmitters. One candidate for this role is acetylcholine, which might be released via either central or peripheral cholinergic innervation of the pineal gland and act through highly-expressed nicotinic acetylcholine receptors [42]. Similarly, neuropeptide Y might play a role; it is present in sympathetic nerves innervating the pineal gland and could act through the neuropeptide Y receptor 1 [43, 44].

Furthermore, the observation that many of the genes induced by norepinephrine are also induced by dibutyryl cyclic AMP treatment supports the view that the effects of norepinephrine are mediated by cyclic AMP. This observation does not, however, preclude the possibility that second messengers other than cyclic AMP might be required for norepinephrine to control the abundance of some genes. Candidates for this role include Ca++, cyclic GMP and phospholipids, all of which are controlled by norepinephrine and could modulate cyclic AMP induction of genes [25, 45–47].

The effects of SCGX and DCN: Similarities and Differences

The ganglionectomy and decentralization surgical procedures were expected to produce somewhat similar phenotypes, as both eliminate or at least greatly reduce the rhythmic stimulation of the pineal gland. By comparing and contrasting these groups with one another and with the sham-surgery group we were able to identify previously unidentified subtleties in the regulation of the pineal transcriptome.

Our analyses detected more statistically significant night/day differentially expressed genes in the SCGX animals than in the DCN animals. This apparent difference could be a coincidence, or could be due to the effects of circulating catecholamines, which gain access to pinealocytes in the absence of sympathetic nerves in the pineal gland of the SCGX animals. The sympathetic innervation of the pineal gland, responsible for uptake of extracellular catecholamines, is eliminated by SCGX, but remains intact in DCN animals, thereby preventing stimulation by these circulating catecholamines [32].

Our comparisons of the DCN and SCGX animals to the sham-surgery animals suggest that the neonatal interruption of neural stimulation does not broadly alter the pineal transcriptome other than preventing rhythmic changes in gene expression. The principal component analysis and the cross-group analyses demonstrated that any differences in the pineal transcriptome resulting from developmental changes in the SCGX or DCN animals are (at worst) minor in comparison to the profound changes that occur between night and day conditions in intact rats (S2 Fig, S2 Dataset). In addition, we found only small changes of three pineal marker genes, Tph1, Asmt, and Gch1, indicating that neural stimulation is not required to maintain these outstanding features of the pineal transcriptome during post-neonatal development.

An interesting observation was made in this study regarding Serpina1, a member of the Serpin family of peptidase inhibitors, which has been well studied biochemically and in human disease [48, 49]. This gene was downregulated following SCGX and DCN, and was upregulated following treatment with norepinephrine or dibutyryl cyclic AMP. This, together with the knowledge that it is a secreted protease inhibitor, provides reason to speculate that the encoded protein mediates neurotransmission-dependent control of the extracellular matrix by inhibiting extracellular proteases, as is the case for another serpin family member, neuroserpin [49, 50]. Serpina1 might mediate neural regulation of cell-cell attachments, sympathetic synaptic architecture, features of vascularization, or some combination thereof.

Impact on future investigations

The results presented here provide a number of valuable resources for future investigators. Tables in the SI provide the abundances and fold changes for all annotated genes and for all analyses. This includes genes exhibiting strong and statistically significant in vivo night/day differential expression and those exhibiting strong and statistically significant in vitro treated/untreated differential expression. Additionally, the SI provides data comparing the abundance of each transcript in the pineal gland to average expression in other tissues. These resources could be useful in identifying potentially productive targets for further study.

Methods

In vivo night/day differential expression of pineal transcripts

To identify and characterize rhythmic transcripts and to study the effects of surgical stimulus deprivation in vivo, pineal gland samples were taken from adult rats in each of four groups and each of two time points. The four surgical groups were: no surgical intervention (Control); neonatal sham surgery (Sham); neonatal bilateral SCG decentralization (DCN); neonatal bilateral superior cervical ganglionectomy (SCGX) [21]. Each surgical group had 18 to 20 animals, which were then divided into two sub groups which were euthanized at mid-day (ZT7) or midnight (ZT19). Samples were pooled into three biological replicates for each of the eight biological conditions (24 pooled samples, total).

In vitro differential expression of pineal transcripts resulting from treatment with norepinephrine or cyclic AMP

To further elucidate the relationship between innervation and transcript regulation, a second experiment was run using pineal gland organ culture[51]. Three treatment conditions were compared: untreated control (Untreated), norepinephrine-treated (NE), and dibutyryl-cyclic-AMP-treated (DBcAMP). For each of the three treatment conditions, samples were pooled into three biological replicates (9 pooled samples, total).

To identify transcripts that are highly expressed in the pineal gland relative to other tissues, pineal gland samples were taken at mid-day (ZT7) or midnight (ZT19). Additional samples were taken from the following tissues at ZT7: cortex, cerebellum, midbrain, hypothalamus, hindbrain, spinal cord, retina, pituitary, heart, liver, lung, kidney, skeletal muscle, small intestine, and adrenal gland. Total RNA was extracted from each tissue, and then equal amounts of each of the 15 tissues were combined for the final pooled “mixed-tissue” sample. Three pooled samples were created: one each for pineal-day, pineal-night, and mixed-tissue (3 pooled samples, total).

Animals

Sprague Dawley rats (Taconic Farms, Germantown, NY) were used for these studies and were maintained on a 14:10 light:dark (L:D) cycle.

For the in vivo differential expression experiments, timed pregnant animals were used; on postnatal day 5 male and female pups were divided into the four surgical groups and subjected to the respective surgical interventions. Animals were then maintained in a 14:10 L:D cycle until day 40. For the in vitro experiments, female rats were used that had been maintained for at least two weeks in a 14:10 L:D cycle before being euthanized (2 to 3 months old).

Animals were euthanized by CO2 asphyxiation and decapitated. Animals sacrificed at ZT19 were decapitated in dim red light. Pineal glands were immediately placed in microtubes on solid CO2 (three glands per tube); samples of other tissues (approximately 10 mg) were rapidly prepared and placed in individual tubes. Tissues were stored at -80°C until use.

Animal use and care protocols were approved by the NIH Institutional Animal Care and Use Committee and followed the guidelines of the National Research Council's Guide for Care and Use of Laboratory Animals (Vol. 8) and the Animal Research: Reporting In vivo Experiments (ARRIVE) guidelines.

Surgical procedures

The neonatal DCN and SCGX were based on modifications to previously published procedures for adult rats (See SI) [21].

Organ culture

For the in vitro experiments, pineal glands were cultured as described [51] for a period of 48 hours under control conditions, then treated with norepinephrine (NE; 1 μM) or dibutyryl cyclic AMP (DBcAMP; 500 μM), or left untreated (CN) for 6 hours before being harvested into microtubes placed on solid CO2 (three glands per tube).

Sample preparation and sequencing

RNA was prepared from pools of pineal glands and from samples of “mixed-tissue” RNA as detailed in the SI.

To minimize sequencing batch effects all samples from each experiment were barcoded and pooled into combined libraries. The strand-specific pooled libraries were sequenced on a HiSeq2000 (Illumina, San Diego, CA) using version 3 chemistry. The resultant reads were aligned with the RNA-STAR aligner, using the rn4/RGSC3.4 genome assembly and Ensembl transcript annotation, release 69 [52, 53]. Quality control metrics were calculated and visualized using the QoRTs software package, and no major artifacts or abnormalities were found [54]. Additional details concerning the sequencing, library preparation, and alignment are described in the SI.

Differential expression analysis

Gene read counts were provided by QoRTs, using the same algorithm defined by the HTSeq package [55]. Differential expression analysis was performed using DESeq2 [56]. In each analysis, genes with a mean normalized read-pair count of less than 5 were not tested.

The pineal marker gene analysis followed a slightly different methodology. Two separate comparisons were run: pineal-day vs mixed-tissue and pineal-night vs mixed-tissue. These comparisons were run without biological replicates, and thus the biological dispersion could not be rigorously estimated and no statistical hypothesis tests could be performed. Instead, the maximum likelihood fold changes were calculated between the pineal and mixed tissue samples using DESeq2. Gene lists were generated using simple thresholds on the mean read-pair count and the maximum likelihood fold changes between the pineal and mixed-tissue samples. We selected a normalized read-pair count threshold of 30 and a fold change threshold of 8.

Supplementary methods.

(DOCX)

Click here for additional data file.

Complete results tables for genes differentially expressed in the pineal gland on a night/day basis, in vivo.

(XLSX)

Click here for additional data file.

Complete results tables for the SCGX/DCN vs Sham analyses.

(XLSX)

Click here for additional data file.

Complete results tables for genes differentially expressed in pineal gland in response to in vitro treatment with norepinephrine or dibutyryl cyclic AMP.

(XLSX)

Click here for additional data file.

Complete results tables for mixed-tissue vs pineal-gland comparisons.

(XLSX)

Click here for additional data file.

Comparison of gene lists detected by the Control day-vs-night RNA-Seq analysis and a previous microarray-based analysis.

Gene lists were compiled using an adjusted-p-value threshold of 0.05 and a fold change threshold of 2 (in either direction). *Note: 22 gene names that were found in the microarray results could not be matched with any current or former gene name of any known gene name in the Ensembl (release 69) or RGD (version 6) databases. Additionally: 155 microarray probe loci met the significance/fold-change criteria, but were not annotated with any known gene. These genes/probes were not counted in the figure above.

(TIF)

Click here for additional data file.

Summary plots for the cross-group analyses, in which SCGX and DCN groups were compared with the Sham group at both day and night.

(a-d) MA plots for all four comparisons. Statistically significant genes (adjusted-p < 0.001) are marked with dots, colored based on their night/day differential status from the Control night/day analysis: blue indicates that the gene is upregulated at night, red indicates that it is downregulated at night, black indicates that the gene is not significantly differentially expressed between night and day (adjusted-p > 0.05). The displayed counts indicate the number of dots of each color above or below the FC > 2 or FC < 0.5 thresholds. A complete list of all genes with fold changes, p-values, and normalized expression estimates is available in the SI (S3 Dataset). (e) Plot of the first two principal components from a principal component analysis across all 24 in vivo samples. Note that based on the first principal component, the Control-Night and Sham-Night samples cluster strongly to the right, while all other samples (including the SCGX and DCN samples for both day and night) cluster strongly to the left.

(TIF)

Click here for additional data file.

MA plot of the comparisons between mixed tissue with pineal glands at both day and night.

Genes were marked in red if the mean normalized read-pair count exceeded 30 and if the (maximum-likelihood) fold change was greater than 8. Note that due to the lack of replicates no significance p-values could be calculated.

(TIF)

Click here for additional data file.

Genes exhibiting differential expression in both in vivo and in vitro experiments.

The genes listed exhibit statistically significant (adjusted-p < 0.001) differentials in all of the following four comparisons: Control night vs day, Sham night vs day, untreated vs NE-treated, and untreated vs DBcAMP-treated.

(DOCX)

Click here for additional data file.

Genes with relatively high expression in the pineal gland.

Gene symbols (when available) of all genes with > 32-fold enrichment in the pineal gland (day and/or night) and mixed non-pineal tissues. Genes that exhibit high relative expression in both the day and night analyses are listed in bold (based on the minimum of the two fold changes). A more complete list is available in the SI (S4 Dataset).

(DOCX)

Click here for additional data file.