Graphical abstract

   graphic file with name fx1.jpg
   [39]Open in a new tab

Highlights

     * •
       Our study provides novel insights into circadian gene expression in
       neuronal circuits
     * •
       Brain regions known for their circadian regulation show low robust
       transcript cycling
     * •
       The mammillary body’s gene oscillations challenge its traditional
       classifications
     * •
       Some related brain regions share distinct gene expression patterns
       while others do not
     __________________________________________________________________

   Neuroscience; Bioinformatics; Expression study

Introduction

   The brain, a remarkably complex organ, is organized into clusters of
   specialized structures responsible for regulating and performing
   specific functions. These structures interact and communicate with each
   other, forming interconnected networks known as neural assemblies or
   neuronal circuits.[40]^1^,[41]^2 These circuits form the basis of
   various brain functions, including sensory processing, motor control,
   and higher-order cognitive
   functions.[42]^3^,[43]^4^,[44]^5^,[45]^6^,[46]^7^,[47]^8^,[48]^9^,[49]^
   10 Neuronal circuits enable the harmonious functioning of both basic
   and higher-order brain activities, orchestrating a wide array of
   physiological processes and behaviors. Most physiological processes
   governed by neuronal circuits exhibit highly conserved circadian
   rhythms, suggesting synchronized activity and communication across
   nuclei within a particular circuit.[50]^11^,[51]^12

   Circadian rhythms are endogenous biological oscillations with a
   periodicity of approximately 24 h that regulate diverse physiological
   processes and behaviors in living organisms. These brain circuits,
   which coordinate functional processes, play a significant role in the
   timing and rhythmic activity of various physiological functions like
   thermoregulation, food intake, and energy balance, to higher-order
   brain activities such as cognition, emotion regulation, and motor
   coordination.[52]^13^,[53]^14^,[54]^15^,[55]^16^,[56]^17^,[57]^18^,[58]
   ^19^,[59]^20

   The suprachiasmatic nucleus (SCN) in the hypothalamus serves as the
   central pacemaker or “master clock” of the mammalian circadian
   system.[60]^21^,[61]^22^,[62]^23 The SCN receives information relating
   to environmental cues, primarily light, and synchronizes endogenous
   circadian clock components to the 24-h light/dark cycle. By employing a
   variety of entrainment mechanisms, the SCN coordinates phase-setting of
   circadian rhythms in other brain regions, ensuring harmonious
   physiological functioning throughout the organism.[63]^24 Beyond the
   SCN, other hypothalamic nuclei also exhibit timekeeping mechanisms and
   contribute to the regulation of circadian-sensitive physiological
   processes and behaviors, such as sleep/wake cycling and daily food
   intake.[64]^25^,[65]^26^,[66]^27 Moreover, rhythmic activity has been
   reported in multiple brain regions and structures, including the
   cortex, hippocampus, and amygdala, where neurons are believed to
   express all essential genes required for cell-autonomous circadian
   oscillations.[67]^20^,[68]^28^,[69]^29^,[70]^30^,[71]^31^,[72]^32^,[73]
   ^33

   We hypothesize that the genetic components within neuronal circuits
   display unique oscillating expression patterns, reflecting the
   synchronization of their respective molecular clocks through
   inter-neuronal communication. To investigate this hypothesis, we
   leveraged the largest available atlas of spatiotemporal gene expression
   in primates to analyze rhythmic transcriptomes across various brain
   regions. We then characterized gene expression patterns within selected
   circuits to better understand the relationship between circadian
   regulation and connectivity between anatomically defined brain regions.

   Our findings reveal that gene expression cycling within brain regions
   traditionally associated with circadian regulation is less robust
   compared to regions not typically considered circadian in nature, and
   some anatomically or functionally related brain regions share distinct
   gene expression patterns while others do not. This observation provides
   novel insights into the encoding and execution of rhythmic behaviors
   and has implications for understanding the impact of circadian rhythm
   disruption on these processes.

Results

Rhythmic gene expression across brain regions

   Our comprehensive analysis identified 15,335 transcripts displaying
   circadian oscillation patterns in at least one brain region. Over 93%
   of these transcripts cycled with a periodicity of approximately 24 h,
   signifying circadian rhythmicity in their expression ([74]Table S2).
   The prefrontal cortex (PRC) exhibited the highest number of circadian
   cycling transcripts (4,276), followed closely by the medial globus
   pallidus (MGP) (4,170) ([75]Figure 1). Conversely, the lateral
   hypothalamus (LH) and amygdala (AMY) exhibited the lowest numbers, with
   371 and 434 circadian transcripts, respectively.

Figure 1.

   [76]Figure 1
   [77]Open in a new tab

   Abundance of circadian transcripts in the primates’ brain regions

   Histogram showing the number of cycling transcripts and their
   percentages of their total transcripts in 22 brain regions. AMY:
   amygdala, ARC: arcuate nucleus, CER: cerebellum, DMH: dorsomedial
   hypothalamus, HAB: habenula, HIP: hippocampus, PRC: prefrontal cortex,
   PUT: putamen, LGP: lateral globus pallidus, LH: lateral hypothalamus,
   MMB: mammillary body, MGP: medial globus pallidus, OLB: olfactory Bulb,
   PVN: paraventricular nucleus, PRA: preoptic area, PON: pons, SUN:
   substantia nigra, SCN: suprachiasmatic nucleus, SON: supraoptic
   nucleus, THA: thalamus, VMH: ventromedial hypothalamus, VIC: visual
   cortex.

   It has been estimated that approximately 10% of genes in most tissues
   exhibit circadian expression.[78]^34^,[79]^35 Based on this, the 22
   brain regions were categorized into three groups: highly rhythmic,
   moderately rhythmic, and low rhythmic ([80]Figure 1). Highly rhythmic
   regions had over 16% rhythmic transcripts, including the MGP, PRC, and
   putamen (PUT) (about 25% each), the lateral globus pallidus (LGP) and
   MMB (24% and 23%, respectively), substantia nigra (SUN) (19%), and
   visual cortex (VIC) (17%). Moderately rhythmic regions had 10–16%
   rhythmic transcripts, including the pons (PON), ventromedial
   hypothalamus (VMH), PVN (15% each) and dorsomedial hypothalamus (DMH)
   (13%). Low rhythmic regions had less than 10% rhythmic transcripts,
   including the olfactory bulb (OLB) (5%), arcuate nucleus (ARC) (6%),
   habenula (HAB), SCN, and thalamus (THA) (9% each). The LH, AMY, and
   supraoptic nucleus (SON) (2.5% each), and hippocampus (HIP), cerebellum
   (CER), and preoptic area (PRA) (4% each) had the lowest rhythmic
   transcripts ([81]Figure 1).

Regional peaking times and temporal patterning in distinct neuronal circuits

   Rhythmic transcripts exhibited region-specific expression peaks,
   reflecting the unique functional and anatomical characteristics of each
   brain region. We grouped brain regions by anatomical and functional
   proximity and then by when their transcript peaks occurred within one
   or more quarter-phase time intervals ([82]Figures 2A, 2B, and [83]S1).
   Brain structures like LGP, MGP, MMB, PON, PRC, PUT, SCN, SUN, and VMH
   showed over 50% of their rhythmic transcripts peaking within a narrow
   6-h time interval or one quarter-phase ([84]Figures 2A, 2B, and
   [85]S1). These peak times coincided with the light phase of the 24-h
   cycle. The first quarter-phase (ZT0-ZT5) included PON, PUT, SCN, and
   VMH, while the second (ZT6-ZT11) had MMB, LGP, MGP, PRC, and SUN
   ([86]Figures 2A, 2B, and [87]S1). The DMH, PVN, and VIC showed peaks in
   two distinct quarter-phase intervals ([88]Figures 2A, 2B, and [89]S1).
   Regions like AMY, ARC, CER, LH, SON, HAB, HIP, OLB, PRA, and THA peaked
   in more than two intervals. Notably, PUT, MGP, and VIC had sharp peaks
   at specific times ([90]Figures 2A, 2B, and [91]S1).

Figure 2.

   [92]Figure 2
   [93]Open in a new tab

   Brain circadian gene expression exhibit region-specific patterns

   (A) Heatmap representation of 24-h oscillation of cycling genes in
   representative regions of the basal ganglia structures, hypothalamic
   nuclei, Cortices, and other brain regions with mean expressions
   normalized between 0 (blue) and 1 (yellow) and ordered by phase
   (p < 0.05), ZT: zeitgeber time.

   (B) Radial diagram of the distribution of the peak phase of expression
   of the cycling genes in representative regions of the basal ganglia
   structures, hypothalamic nuclei, Cortices, and other brain regions. The
   radial plot displays phases (hr) on the circumference and the number of
   gene expression peaks on the radius.

   (C) Cycling genes in brain regions mediating motor activity exhibit
   sharp expression peaks opposing to cycling genes of hypothalamic
   nuclei; Histogram showing number of cycling gene expression in the
   basal ganglia, hypothalamus, and PRC regions with the x axis
   corresponding to the hour of day, ZT: zeitgeber time.

Divergent temporal patterns in basal ganglia and hypothalamic circuits

   Considering that the brain regions analyzed for cycling transcripts in
   the original study are part of distinct anatomical or functional
   circuits, we investigated the patterns of cycling transcripts in
   interacting nuclei. We determined inter-nuclei interactions based on
   relative anatomical position or functional connectivity and examined
   them for discernible patterns of overlapping cycling gene expression.

   The basal ganglia (BG) circuit, which includes the PUT, LGP, MGP, and
   SUN has connections with THA and PRC, and is crucial for motor control,
   learning, and executive functions. High levels of cycling transcripts
   were seen in these regions ([94]Figure 1), especially in PUT and MGP,
   the input and output regions of the BG, respectively. The rhythmic
   transcripts within the BG structures exhibited sharp peaking patterns,
   unfolding in a sequentially ordered manner from one region to the next.
   The PUT showed a strong circadian peak at ZT4 in the morning, whereas
   the SUN and LGP had similar peak timings, and together with MGP,
   followed PUT with sharp peaks 3–4 h later (ZT6-ZT7). Peak rhythmicity
   in the PRC occurred post-BG structures, around ZT8-ZT9 ([95]Figures 2B
   and 2C).

   In contrast to the circadian patterns observed in BG structures, the
   hypothalamic nuclei (ARC, DMH, LH, PRA, PVN, SCN, SON, and VMH) showed
   a more dispersed circadian transcript pattern, whereas MMB’s pattern
   resembled that of BG structures more than hypothalamic nuclei
   ([96]Figures 2B and 2C). Strikingly, peak times of cycling transcripts
   in the BG structures and hypothalamic nuclei were inversely related.
   The transcripts in BG structures reached their peak during the light
   phase, while those in the hypothalamic nuclei peaked during the dark
   phase ([97]Figure 2C).

The number of cyclic genes in a brain region is related to its input and
output connectivity patterns

   We observed several interesting and surprising patterns in our data,
   including many of the brain regions that showed high levels of cyclic
   genes were in the BG circuit. However, brain connectivity is more
   complicated than classical models suggest, indicating that if we want
   to assess the relationship between connectivity and gene expression
   patterns, we need a quantitative, whole-brain mesoscale connectivity
   atlas. While this does not exist in primate species, the Allen Brain
   Institute provides an open access database of connectivity in the mouse
   brain.[98]^36 By using this rodent atlas, we could approximate the
   strength of connections between brain regions. While connectivity in
   the rodent brain is certainly different than connectivity in the
   primate brain, many key brain structures are conserved between the
   species, allowing us to use this as a first-order approximation of the
   relationship between cyclic gene expression and connectivity.

   We first built a model network between the brain regions we are
   examining in this study. Using the Louvain Community Detection
   algorithm, we detected three communities of brain regions
   ([99]Figure 3A). The Louvain algorithm identifies communities—groupings
   of nodes that maximize modularity such that nodes within community
   preferentially connect with each other over those outside of their
   community. The PUT, THA, and PON were the most extensively
   interconnected of these regions, as indicated by their high betweenness
   centrality score ([100]Figure 3A). We also created a correlogram to
   detail the quantitative relationship between input/output connectivity
   for each pair of brain regions ([101]Figure 3B). As expected, regions
   in the BG cluster into a common community (Community 2), while most
   hypothalamic regions cluster into Community 1. Community 3 contains
   mostly thalamic, epithalamic (HAB), and hippocampal regions, with the
   exception of the SCN; this is likely because the SCN receives strong
   innervation from the lateral geniculate nucleus (LGN) of the thalamus,
   pulling it into a community with THA and not hypothalamic regions.

Figure 3.

   [102]Figure 3
   [103]Open in a new tab

   Relationship between cyclic gene expression and network connectivity

   (A) Network connectivity matrix between the brain regions analyzed in
   this dataset. The Louvain Community Detection Algorithm was used to
   assign regions into Communities based on their connectivity patterns.

   (B) Correlogram indicating the strength of input-output connections
   between all pairs of regions in the dataset. The color bar indicates
   the log-normalized projection density between any two regions.

   (C) Linear regression plot of the relationship between the number of
   cyclic genes expressed in each region and the betweenness score of that
   region. No significant association was detected (r^2 = 0.03, p = 0.43).

   (D) Visualization of the network connectivity matrix between the brain
   regions analyzed in this dataset, with the node size reflecting the
   number of cyclic genes expressed in that region.

   (E) Bar graph of the number of cyclic genes in each region, binned by
   community. The number of cyclic genes was significantly different
   between these three communities (One-way ANOVA, p = 0.05), with regions
   in Community 2 expressing higher numbers of cyclic genes relative to
   Community 1 (t-test, p = 0.06) and Community 3 (p = 0.13). Communities
   1 and 3 expressed similar numbers of cyclic genes (p = 1.0). T-tests
   include corrections for multiple comparisons.

   (F) Linear regression plot of the relationship between the numbers of
   cyclic genes shared between two brain regions, and the log normalized
   projection density of the connection between them. A highly significant
   association was detected (r^2 = 0.04, p = 0.004). The pairs of brain
   regions with < 200 shared cyclic genes were excluded from this
   analysis.

   Based on this connectivity matrix, we first wanted to explore whether
   the number of connections each brain region received was related to the
   number of cyclic genes that the region expressed. We did not find a
   significant correlation between these metrics (p = 0.43; standard error
   of slope = 7,397, Y-intercept = 355; [104]Figure 3C), indicating that
   whether or not a region serves as a hub gene has no effect on the
   number of cyclic genes it expresses. Next, we aimed to explore the
   relationship between the innerconnectivity of any two regions and the
   numbers of cyclic genes that are shared between these regions. We first
   visualized the connectivity network, using the same edges as in
   [105]Figure 3A, and visualized each node by correlating the size with
   the total number of cyclic genes expressed in that brain region
   ([106]Figure 3D). Using this approach, we find that brain regions in
   Community 2, which contains regions in the BG, have higher levels of
   cyclic genes than regions in Communities 1 or 3 ([107]Figure 3D). We
   quantified the number of cyclic genes in each region within each
   community, and found a strong trend toward a significant difference
   between groups (One-way ANOVA, p = 0.053; [108]Figure 3E). We followed
   up on this finding by building a mixed effects model, modeling
   community membership as a random effect on number of cyclic genes.
   Community membership accounted for 26% of variance in number of cyclic
   genes due to random effects, suggesting community membership is related
   to the number of cyclic genes expressed in a given brain region.
   Notably, we excluded the MMB from the previous analysis, given that its
   cyclic gene expression pattern more closely aligns with patterns
   observed in BG structures (as described later in this manuscript).

   Lastly, to more quantitatively explore the relationship between the
   number of shared cyclic genes and the strength of projection between
   them, we performed a linear regression analysis of the numbers of
   shared cyclic genes between pairs of regions relative to the
   log-normalized projection density between them ([109]Figure 3F). Region
   pairs with 200 or fewer shared cyclic genes were discarded as they were
   overrepresented in our data. We found a statistically significant
   correlation, indicating that the number of shared cyclic genes between
   two regions is strongly related to how directly interconnected those
   two regions are (p = 0.004; standard error of slope = 214.5,
   Y-intercept = 143). Data failed a test for heteroscedasticity using the
   Breush-Pagan test (p = 0.006), so we tested the linear coefficient p
   value for bias by performing scrambled linear regression tests. We
   scrambled the relationship of number of shared genes with the
   projection strength and performed a linear regression test 1000 times
   and recorded all the linear coefficient p values. The p value for the
   actual data were less than 99.5 percent of all scrambled p values,
   suggesting that it is still statistically significant. Even so, we also
   tested these data for monotonicity with a Spearman correlation test,
   and found a positive relationship between them (p = 0.022), providing
   further evidence of a relationship between connectivity and the number
   of shared cyclic genes between brain sites.

Overlapping expression of circadian transcripts across brain regions

   We found varying overlaps in circadian transcript expression across
   brain regions, with two (Metazoa_SRP and U6) shared globally
   ([110]Table S3). Half the regions shared 106 transcripts, while 5220
   transcripts were unique to individual regions.

   There were 231 regional overlaps of the circadian transcriptome, with
   the number of overlapping transcripts ranging from 15 (between AMY and
   LH) to 2368 transcripts (between LGP and MGP, [111]Figure 4A). MGP had
   the highest number of cyclic transcripts overlapping with other regions
   (3969). Besides the MGP, the LGP, SUN, MMB, PRC, PUT, and SCN shared
   over 90% of their cycling transcripts ([112]Figure 4B). The LH had the
   least overlap (259 transcripts, 70% interregional) ([113]Figure 4B).
   Most regions showed the greatest number of cyclic transcript overlaps
   with the MGP, PRC, LGP, and MMB ([114]Figure S2). BG structures
   exhibited a higher overlap of circadian transcripts among each other
   than with other brain regions, with each nucleus sharing at least 45%
   of their cycling transcripts with the others ([115]Figures 5A and
   [116]S2; [117]Table S3). They shared 1083 circadian transcripts among
   themselves, and 871 with the PRC ([118]Figures 5A and [119]S2;
   [120]Table S3). These shared transcripts showed sharp peaks at specific
   times: 4ZT for PUT, 6-7ZT for SUN and LGP, 7ZT for MGP, and 8-9ZT for
   PRC ([121]Figure 5B). Interestingly, the order of these peaks aligned
   with the sequence of anatomical connectivity in controlling movement.

Figure 4.

   [122]Figure 4
   [123]Open in a new tab

   Overlapping expression of circadian transcripts across brain regions

   (A) Heatmap of region-by-region overlap showing overlapping expression
   of every region’s circadian transcriptome with all other brain regions
   contributing a total of 231 regional overlaps.

   (B) Percentage (y axis) of overlapped genes in each brain region (x
   axis). MGP showed highest rate of overlapping circadian genes with
   other regions (95% of its circadian gene) while LH showed the lowest
   overlapping of cycling genes with other regions (<70% of its circadian
   genes).

Figure 5.

   [124]Figure 5
   [125]Open in a new tab

   Rhythmic patterns of the shared cycling transcripts in the basal
   ganglia-cortex and hypothalamic circuits

   (A) Venn diagram showing the shared cyclic transcripts across the
   structures of the basal ganglia circuit and prefrontal cortex.

   (B) Radial plot of the peak phase of expression of the shared circadian
   genes across the structures of the basal ganglia and the prefrontal
   cortex.

   (C) Venn diagram showing the overlap of regions of the basal ganglia
   structure with MMB.

   (D) Radial plot of the distribution of the peak phase expression of
   shared circadian genes in the basal ganglia and MMB.

   (E) Radial plot showing peak expression phase of the clock genes with
   circadian pattern in basal ganglia and hypothalamic structures.

   Circadian transcripts in hypothalamic structures showed minimal overlap
   with each other and other brain regions ([126]Figure S3). Given the low
   degree of interregional overlap across hypothalamic structures
   (excluding the MMB), we analyzed shared transcripts across different
   hypothalamic nuclei, assuming that hypothalamic circuits encompassed at
   least three nuclei. This analysis produced 84 sets of three
   hypothalamic nuclei ([127]Table S3). The overlap of circadian
   transcripts across these sets ranged from five transcripts (in
   ARC/LH/SON and MMB/LH/SON sets) to 348 transcripts (in DMH/MMB/PVN set)
   ([128]Table S3).

   Though traditionally classified as a hypothalamic nucleus due to its
   location, MMB displayed greater gene overlap with BG structures (PUT,
   LGP, MGP, and SUN) than other hypothalamic nuclei ([129]Figures 5C and
   [130]S3). Over 1000 circadian transcripts were shared between MMB and
   any two BG structures, contrasting with only 106–348 shared with any
   two hypothalamic nuclei ([131]Figures 5C and [132]S3). Shared circadian
   genes across MMB and BG showed a pronounced peak at 4ZT for PUT and at
   7ZT for LGP, MGP, and MMB ([133]Figure 5D).

   Finally, clock genes (Dbp, Per1, Per2, Per3, Cry1, Cry2, Clock, Bmal1,
   Bmal2, Npas2, Ror2, Nr1d1, and Nr1d2) exhibited circadian rhythms in at
   least one BG and hypothalamic structures, except the LH
   ([134]Figure 5E). Notably, the BG structures PUT and MGP expressed the
   highest number of circadian clock genes (seven), whereas the
   hypothalamic MMB expressed the most circadian clock genes (eight)
   ([135]Figure 5E).

Pathway enrichment analysis of cycling transcripts

   We performed a pathway enrichment analysis on cycling transcripts in
   each brain region, calculating a Log2 fold enrichment (Log2(FE), with a
   cut-off of a false discovery rate (FDR) < 0.05. Brain regions
   displaying a Log2(FE) greater than 3 are presented in [136]Figure 6A,
   and the rest are presented in [137]Figure S3. We classified brain
   regions by the density of associated pathways into low, moderate, or
   highly enriched groups, based on the number of enriched pathways. Low
   enrichment regions were AMY (enriched in transmembrane receptor protein
   tyrosine kinase signaling pathway and circadian rhythm), LH
   (neuropeptide signaling pathway), OLB (extracellular matrix
   organization), PRA (cell communication regulation), SCN
   (posttranslational gene expression regulation), and SUN (base excision
   repair) ([138]Figures 6A and [139]S4). Moderate enrichment pathways’
   numbers were seen in ARC, DMH, MMB, PON, PVN, and THA, whereas HAB,
   LGP, MGP, PRC, PUT, SON, VIC, and VMH exhibited numerous enrichment
   pathways ([140]Figures 6A and [141]S4). Notably, HIP and CER showed no
   pathway enrichment, and the number of cycling transcripts in certain
   regions did not correlate with enriched pathways.

Figure 6.

   [142]Figure 6
   [143]Open in a new tab

   Enrichment of the biological pathways of cycling transcripts in
   different brain regions

   (A) Enriched pathways in different brain regions displaying a Log2(FE)
   greater than 3.

   (B) Enriched pathways in shared circadian transcripts across two-third
   (15) brain regions and (right) 50% of all the brain regions (11
   regions).

   (C) Enrichment pathways in shared circadian transcripts in the basal
   ganglia.

   Circadian rhythm and related processes were enriched pathways in
   cycling transcripts of many regions, particularly the AMY, MMB, HAB,
   PRC, SON, THA, VIC, and VMH. Pathways regulating proteasomal
   ubiquitin-dependent protein processes were enriched in cycling
   transcripts of the DMH, LGP, MGP, PRC, PUT, PVN, SUN, VIC, and VMH.
   Microtubule-based process pathways were enriched in the ARC, DMH, MGP,
   PRC, PVN, and VIC. Specifically, regulatory pathways involving
   microtubule-based organelles, such as cilia and centrosomes, were
   enriched in the ARC and DMH, respectively ([144]Figures 6A and
   [145]S4).

   Pathways linked to RNA splicing and post-transcriptional regulation
   were enriched in the HAB, PON, SCN, THA, VMH, and VIC. Mitochondrial
   processes were significant in MMB, PUT, PON, VMH, and VIC. MGP, HAB,
   and PRC showed enrichment in axonal myelination pathways, whereas LGP,
   PUT, and VIC favored the ERAD (endoplasmic reticulum to Golgi
   vesicle-mediated transport) pathway. Lipid response and metabolism
   pathways were enriched in MGP, DMH, HAB, PRC, PVN, and MMB, whereas DNA
   transcription and damage response pathways were enriched in DMH, HAB,
   MGP, MMB, LGP, PON, PVN, PUT, and VIC. Telomere maintenance pathways
   were active in LGP, PON, and VMH, while HAB, PON, VMH, and MGP were
   enriched in histone modification pathways. Synapse development pathways
   were enriched in SON and VMH, with SON, VMH, and PRC. The regulatory
   pathway of calcium-dependent exocytosis was enriched for the PRC, PVN,
   and AMY, and protein localization to the nucleus was enriched for the
   LGP and PUT. Response to insulin and regulation of cell glucose
   pathways were enriched for the HAB and VMH. The neurogenesis pathway
   was enriched for the SON, HAB, and PRC, and the WNT pathway was
   enriched for the LGP and PUT ([146]Figures 6B and 6C).

   We conducted a gene analysis across 15 brain regions, representing
   two-third (68%) of all areas studied. The results revealed shared
   biological pathways related to circadian rhythm, positive triglyceride
   biosynthesis regulation, negative peptidyl-threonine phosphorylation
   regulation, ncRNA transcription regulation, base-excision repair,
   histone deacetylation, and JUN kinase activity regulation
   ([147]Figure 6B).

   Shared cycling genes across the BG structures showed enrichment in
   pathways for positive cholesterol efflux regulation, chondroitin
   sulfate proteoglycan metabolism, selective autophagy, snRNA metabolism,
   endoplasmic reticulum stress response regulation, fatty acid
   catabolism, NADH dehydrogenase complex assembly, mitochondrial
   respiratory chain complex I assembly, and positive regulation of
   proteasomal ubiquitin-dependent protein catabolism ([148]Figure 6C).
   For hypothalamic circuits, we analyzed shared transcripts across 84
   sets of three hypothalamic nuclei. This analysis, however, did not
   reveal any significantly enriched pathways.

Discussion

   In this study, we investigated whether genetic components of neuronal
   circuits display distinct oscillating expression patterns based on
   synchronized molecular clocks in interneuronal communication. We also
   tested if brain regions with more rhythmic physiological functions
   express more circadian genes than those traditionally less linked with
   circadian function. We analyzed publicly available transcriptomic data
   from a non-human primate’s spatiotemporal gene expression atlas.
   Unexpectedly, we found that nuclei involved in rhythmic processes show
   moderate to low transcript rhythmicity, and some anatomically or
   functionally related brain regions share distinct gene expression
   patterns while others do not.

Methodological consideration

   While numerous studies have analyzed oscillatory expression patterns in
   peripheral tissues of humans and animals, including primates and
   rodents, only a few have focused specifically on brain regions and
   assessed the data derived from the underlying circuits. The dataset
   used here encompassed an atlas of transcriptomic data from 64 tissues,
   a third of which were brain regions.[149]^34 Our analysis revealed
   significant differences in the number of circadian transcripts
   identified by BIO_CYCLE, used in our study, compared to MetaCycle, used
   in the original study. Except for the paraventricular nucleus (PVN),
   BIO_CYCLE identified more circadian transcripts in all brain regions
   than MetaCycle.[150]^37 BIO_CYCLE, being a Bayesian algorithm, provides
   a statistically robust framework for detecting cyclical patterns in
   time-series data. Its in-built features, including the ability to
   handle missing data, the use of phase-adjusted errors, and the
   assessment of periodicity across multiple cycles, make it particularly
   suitable for our analysis. Furthermore, BIO_CYCLE’s Bayesian framework
   enables more effective handling of uncertainty, which is beneficial in
   the context of high-throughput transcriptomic data, where intrinsic
   variability is common.[151]^38^,[152]^39^,[153]^40^,[154]^41^,[155]^42

   As noted in the results section, we used a mouse mesoscale brain atlas
   to model connectivity, and related this to the cyclic gene expression
   patterns observed in the primate brain. Of course, the brain-wide
   connectivity patterns between rodents and primates are not wholly
   conserved. In addition, the mouse is a nocturnal species, while the
   baboon is diurnal; there is thus some possibility that the
   relationships between cyclic gene expression and anatomical
   connectivity may differ between the two species. Thus, a degree of
   caution should be taken in directly extrapolating these results
   directly into the relationship between cyclic gene expression and
   connectivity in baboons. Nonetheless, despite these key differences in
   animal behavior, the interconnectivity between key brain structures
   such as the BG is largely conserved across species. This certainly is
   the case between those evolutionally relatively closely related, such
   as the mouse and primate, but also in species that are less closely
   related, such as mice and birds and
   reptiles.[156]^43^,[157]^44^,[158]^45 In all, the results relating
   connectivity to cyclic gene expression are best viewed as a first-order
   approximation of the relationship between these factors, rather than a
   rigorous quantitative assessment. Nonetheless, we believe that the key
   finding of these analyses—that the number of shared cyclic genes
   between two brain regions is related to the interconnectivity between
   them—is supported by the evidence presented in our manuscript. This
   finding also is intuitively logical, and suggests that brain regions
   that function within a network have similar shared transcripts, and
   suggests that the cyclic nature of these genes extends beyond one brain
   region, but likely impacts the functional output of entire circuits.
   How entire circuits may be impacted by the shared cyclic nature of
   these genes, however, remains to be explored.

Brain structures typically associated with regulating circadian processes do
not display robust cycling patterns in their gene expression

   It is reasonable to presume that brain structures mediating more
   rhythmic functions, such as hypothalamic nuclei, hippocampus, and
   amygdala, would express more circadian genes. This presumption, indeed,
   influenced prior research, evident in the original transcriptomic atlas
   focusing on circadian genes from nine hypothalamic subregions, whereas
   other regions, such as the striatum with its distinct nuclei, were
   underrepresented or neglected, focusing only on the putamen. Vital
   structures like the ventral tegmental area and SUN pars reticulata (a
   subdivision of the SUN) were also omitted.[159]^34

   Contrary to expectations, our analysis revealed that hypothalamic
   nuclei, which govern various rhythmic physiological processes,
   exhibited only moderate to low gene rhythmicity, with the MMB as an
   exception.[160]^28^,[161]^46 Notably, the LH, a key player in circadian
   functions, displayed the lowest rhythmicity. Our enrichment analysis
   revealed a limited number of defined circadian pathways, although the
   enriched pathways did involve neuropeptides associated with circadian
   functions like feeding behavior and sleep. Similar to the LH and other
   hypothalamic structures, the HIP, known for its circadian activity and
   synaptic plasticity
   modulation,[162]^47^,[163]^48^,[164]^49^,[165]^50^,[166]^51^,[167]^52
   also showed very low levels of rhythmic transcripts. Our findings
   revealed that the PRC, MMB, and BG nuclei (PUT, LGP, MGP, and SUN) have
   the highest number of circadian genes, a surprising result given the
   traditional association of BG nuclei with motor control, learning, and
   executive functions rather than circadian systems. The underexplored
   MGP, showing the highest proportion of cycling transcripts across all
   brain regions, and LGP present an unexpected aspect of rhythmic
   process, aligning with recent findings of a high proportion of
   circadian genes in human PUT.[168]^53^,[169]^54^,[170]^55 Remarkably,
   the peak times of cycling transcripts in the BG structures and
   hypothalamic nuclei were inversely related. The cyclic transcripts in
   BG structures reached their peak during the light phase, while those in
   the hypothalamic nuclei peaked during the dark phase. This contrasting
   expression is interesting in light of the distinct functional roles of
   the BG in motor functions and hypothalamic circuits in sleep.

   The MMB, currently classified as a hypothalamic structure, demonstrated
   high rhythmicity akin to the BG nuclei and PRC rather than its
   hypothalamic counterparts. Given MMB’s involvement in memory and
   movement control and its strong gene cycling correlation with BG, we
   suggest reconsidering its categorization. This could also have
   implications for disorders like Wernicke-Korsakoff syndrome, affecting
   the MMB, as it could enhance our understanding of motor-related aspects
   and refine treatment strategies.[171]^56^,[172]^57

Pathway enrichment analysis of cyclic transcripts reveals discrete roles of
different brain structures

   Our pathway enrichment analysis of cycling transcripts demonstrated
   distinct roles of various brain structures. Notably, proteasomal
   ubiquitin-dependent protein degradation pathways were enriched in 9 out
   of 22 brain regions. These pathways are responsible for breaking down
   proteins and recycling amino acids to create new proteins. Given that
   organisms require diverse proteins for time-dependent and
   region-specific functions and protein degradation rates fluctuate
   widely, robust cycling behavior in gene expression for these pathways
   is anticipated to match daily cellular activity
   variations.[173]^58^,[174]^59^,[175]^60^,[176]^61^,[177]^62^,[178]^63^,
   [179]^64^,[180]^65 The malfunction of this pathway due to deficient
   Ubiquitin ligase (UbL) expression has been linked to disrupted
   circadian systems and impaired neuropeptide transport in animal UbL
   knockout
   models.[181]^58^,[182]^60^,[183]^62^,[184]^63^,[185]^64^,[186]^65
   Microtubule-based processes, facilitating intracellular protein
   transport, were found enriched in most of the same brain regions as
   this protein regulatory pathway.[187]^59 Additionally, an overlap was
   observed between areas where ERAD, a vesicle-mediated transport
   pathway, and Ub regulatory pathways were enriched (i.e., LGP, PUT, and
   VIC). This implies a deeper functional interrelationship among these
   pathways, contributing to similar pathway enrichments across multiple
   brain regions.

   Despite its role in circadian-sensitive memory processes, our analysis
   found that the HIP did not show pathway enrichment for any cycling
   gene. The HIP clock, governed by the transcription factor BMAL1,
   influences time-dependent memory deficits, and reduced BMAL1 levels
   impact neurogenesis in the HIP.[188]^66^,[189]^67^,[190]^68 It is
   plausible that the HIP function could be indirectly affected by the
   cyclical gene expression patterns in other areas, like the HAB, a
   region interconnecting the BG and the limbic system, thereby
   influencing the limbic memory circuit that relies on the HIP.[191]^69

Functional and anatomical circuit-driven transcriptomic circadian rhythm

   We sought to determine whether cycling genes in anatomically and
   functionally defined brain circuits, such as the BG and hypothalamic
   structures, exhibit coordinated oscillatory patterns. These patterns
   could shed light on the information flow within neuronal circuits,
   giving rise to complex behaviors.

Circadian oscillator network within the hypothalamus

   The hypothalamus, composed of interconnected nuclei, and housing
   central “master” clock in the SCN, acts as a network of oscillators
   controlling numerous physiological and behavioral aspects. Circadian
   clocks, which are sensitive to temperature changes, integrate cues such
   as body temperature fluctuations, informing thermoregulatory
   mechanisms. This coordinated activity between SCN and other nuclei is
   vital for rhythmic temperature control. Moreover, SCN-independent
   hypothalamic clocks help synchronize diurnal energy expenditure and
   thermogenesis.[192]^70^,[193]^71

   Key thermogenic circuits involve the ARC, PRA, DMH, and LH, which help
   regulate body temperature and metabolism.[194]^72 The thermoregulatory
   processes impact metabolic activities like energy expenditure and
   adaptive thermogenesis, which responds to environmental changes like
   temperature fluctuations and food intake. Studies have shown body
   temperature fluctuations coincide with metabolic heat production in
   constant darkness and temperature-controlled environments.[195]^73
   Energy homeostasis, achieved by balancing energy expenditure with food
   intake, involves hypothalamic nuclei like the ARC, DMH, PVN, and LH,
   which cooperate with the SCN to generate feeding
   rhythms.[196]^74^,[197]^75^,[198]^76^,[199]^77^,[200]^78 Neural
   networks connecting these hypothalamic nuclei express signaling
   molecules affecting appetite, such as proopiomelanocortin (POMC),
   neuropeptide Y (NPY), and agouti-related peptide (AgRP).[201]^79 The
   sleep/wake cycle, another circadian-related behavior, is controlled by
   neural connections from the DMH and LH to the
   PRA.[202]^80^,[203]^81^,[204]^82 The DMH, receiving direct input from
   the SCN, is instrumental in maintaining the sleep/wake cycle’s
   circadian rhythm.[205]^83^,[206]^84 Its neuronal projections to the POA
   and LH support rhythmic changes, enabling transitions between
   wakefulness and sleep.[207]^84 Despite the known interactions among
   hypothalamic nuclei, our findings surprisingly reveal no significant
   overlap among cyclic transcripts or any metabolic or functional
   pathways across these nuclei. This may suggest more parallel than
   hierarchical or synchronized action. We recognize that our analysis of
   bulk RNAseq data have a significant limitation regarding its ability to
   reflect the cellular heterogeneity within each brain nucleus. Notably,
   structures such as the hypothalamic nuclei are comprised diverse
   neuronal populations, each of which could exhibit distinct rhythmic
   gene expression patterns. This complexity might be masked when
   analyzing bulk tissue transcriptomic data, potentially leading to an
   underrepresentation or overlook of individual population-specific
   rhythmicities. This potential confounder is particularly relevant when
   considering that distinct neuronal populations within the same nucleus
   might have gene expression peaks at varying times of the day. Such
   complexities could result in an apparent lack of rhythmicity at the
   tissue level, even as a variety of rhythmic patterns occur at the
   cellular level.

Circadian oscillator network within the basal ganglia circuit

   The BG circuit, comprising the PUT, MGP, LGP, and SUN, controls
   movement, learning, executive functions, and emotional processing. Our
   study identified rhythmic dopamine D2 receptor expression in the PUT,
   with most rhythmic genes in the BG nuclei clustering in one or two
   narrow time windows. The order of peak times for shared rhythmic genes
   was PUT > SUN > LGP/MGP, occurring within 1–2 h. However, the precise
   nature of these circadian oscillations is unclear due to the
   uncertainty regarding SUN tissues analyzed, which could involve
   dopaminergic (SUNc) or GABAergic (SUNr) neurons. The sequential peak
   times might be linked to changes in locomotor activity, given the
   variation in activity levels and molecular factors across BG nuclei.
   During rest periods, striatal neurons are generally inactive, with
   information processing and relay occurring during movement
   initiation.[208]^85^,[209]^86 The complex architecture of BG circuitry
   might underpin the sequential peak times, with structures coordinating
   their activities after periods of sleep-related inactivity. While
   locomotor activity exhibits circadian rhythms, traditional views do not
   consider BG structures involved in motor activity as part of the
   circadian system. Nonetheless, our findings indicate high circadian
   transcript numbers in these regions, possibly due to motor activity
   suppression during sleep and the need for rapid activation upon waking.

   Sensory decoupling and sleep-related muscle relaxation prevent
   execution of motor plans in response to external or internal
   triggers.[210]^87 It is plausible that circadian rhythm may influence
   mechanisms underlying these processes, with genes showing oscillating
   expression patterns in relevant brain regions corresponding to
   sleep/wake cycling. Moreover, BG neurons’ spatiotemporal expression of
   rhythmic genes follows sequential order of functional circuitry
   controlling movement. Our findings suggest that drugs targeting BG
   circuits (e.g., via D2 receptors) may be influenced by rhythmic
   cycling. As BG abnormalities are implicated in schizophrenia, Tourette
   syndrome, Parkinson’s disease, and obsessive-compulsive disorder, our
   results may have clinical implications for therapeutic administration
   timing. Further research on circadian rhythmicity and
   pharmacotherapeutics’ dosing schedules could enhance medication
   efficacy by leveraging natural cycling of circadian processes.

Concluding remarks

   Our investigation into spatiotemporal gene expression patterns in brain
   regions with related anatomical and functional properties has revealed
   interesting insights into circadian processes. We found that brain
   structures traditionally associated with circadian functions, like
   hypothalamic nuclei, do not always exhibit strong gene cycling.
   Conversely, regions like the BG, not typically linked to circadian
   physiology, display dynamic cycling patterns. The MMB, often considered
   hypothalamic nuclei, exhibited gene oscillations similar to BG nuclei,
   potentially impacting their classification. Further investigation into
   the interregional effects of pathway activation could enhance our
   understanding of molecular activity’s correlation with circuit
   functions of different brain structures.

Limitations of the study

   We recognize the limitations of our study associated with using
   publicly available data and the intrinsic limitations tied to
   transcriptomic analysis. The extrapolation of findings from a non-human
   primate model to human physiology should be done with caution. Our
   focus on brain structures excludes potential connections with key
   peripheral organs such as the liver, suggesting an avenue for future
   research. Additionally, our largely computational approach lacks
   experimental validation from model organisms, restricting the
   establishment of causal relationships between gene oscillations and
   their roles across brain regions. Despite these caveats, our work is
   expected to provide a solid foundation for future studies exploring the
   complexities of circadian rhythms.

STAR★Methods

Key resources table

   REAGENT or RESOURCE SOURCE IDENTIFIER
   Deposited data
     __________________________________________________________________

   Transcriptome of baboon brain (22 regions) (Mure et al.[211]^34)
   [212]GSE98965
   Allen Mouse Brain Connectivity Atlas (Oh et al.[213]^36)
   [214]https://doi.org/10.1038/nature13186
   Chimpanzee gene database PANTHER 17.0
   [215]https://doi.org/10.5281/zenodo.7942786
     __________________________________________________________________

   Software and algorithms
     __________________________________________________________________

   BIO_CYCLE CircadiOmics [216]https://circadiomics.igb.uci.edu/
   NetworkX (Hagberg,[217]^89) [218]https://pypi.org/project/networkx/
   Plotly (Sievert,[219]^90) [220]https://github.com/plotly/plotly.py
   GOGrapher (Muller et al.[221]^91)
   [222]https://github.com/fjukstad/gographer
   SciPy (Virtanen et al.[223]^92) [224]https://scipy.org/
   lme4 (Bates et al.[225]^93) [226]https://github.com/lme4/lme4/
   [227]Open in a new tab

Resource availability

Lead contact

   Further information and requests for resources and reagents should be
   directed to and will be fulfilled by the lead contact, Amal Alachkar
   ([228]aalachka@uci.edu).

Materials availability

   This study did not generate new unique reagents.

Method details

Transcriptomic data

   Transcriptomic time series in 22 brain regions in baboon were retrieved
   from the publicly available database.[229]^34 The original study
   consisted of twelve male baboons, aged between 5 to 6 years and
   weighing between 7 to 11 kilograms. A total of 284 brain sections were
   collected, amounting to 22 to 25 per animal. Samples were collected
   from one animal every 2 hours over a 24-hour period. To ensure the
   validity of the gene expression data, only genes with an average FPKM
   value greater than 0.1 across the 12 time points were deemed as
   expressed.[230]^34 All tissues were immediately frozen on dry ice and
   stored at -80°C, preserving the integrity of RNA and ensuring
   high-quality sequencing data.[231]^34 Brain regions included: amygdala
   (AMY), arcuate nucleus (ARC), cerebellum (CER), dorsomedial
   hypothalamus (DMH), habenula (HAB), hippocampus (HIP), prefrontal
   cortex (PRC), putamen (PUT), lateral globus pallidus (LGP), lateral
   hypothalamus (LH), mammillary body (MMB), medial globus pallidus (MGP),
   olfactory bulb (OLB), paraventricular nucleus (PVN), preoptic area
   (PRA), pons (PON), substantia nigra (SUN), suprachiasmatic nucleus
   (SCN), supraoptic nucleus (SON), thalamus (THA), ventromedial
   hypothalamus (VMH), and visual cortex (VIC).

Circadian analysis

   The transcriptomic time series were extracted and analyzed through
   CircadiOmics,[232]^39^,[233]^42^,[234]^88 the largest repository and
   web portal for circadian -omic time series datasets. BIO_CYCLE, a
   deep-learning-based software developed to analyze periodicity in -omic
   time series, was used to identify statistically significant circadian
   transcripts and to provide the corresponding amplitudes and
   phases.[235]^37 BIO_CYCLE is trained and calibrated using both
   synthetic and real-world biological time series datasets containing
   labels for periodic and aperiodic signals. It uses a classification
   deep neural network (DNN) to discriminate periodic from aperiodic
   signals and a regression DNN to estimate the signal's period, phase,
   and amplitude. BIO_CYCLE computes p-values, and we used a cutoff of
   p=0.05 in order to consider a gene's transcript as being circadian.
   BIO_CYCLE is publicly available on the CircadiOmics web portal at:
   [236]http://circadiomics.igb.uci.edu.

   Cycling transcripts were clustered into four discrete time phases
   associated with four 6-hour intervals (quarter phases) depending on
   their expression peaks. Each phase is defined as follows: first
   quarter-phase (ZT0-ZT5), second quarter-phase (ZT6-ZT11), third
   quarter-phase (ZT12-ZT17), and fourth quarter-phase (ZT18-23). The
   overlapping expression of rhythmic transcripts across different brain
   regions was determined through gene-view and region-view overlap. The
   distribution of circadian transcripts across functional brain circuits
   was examined and grouped based on regional localization of expression.

Network Construction

   The structural network was constructed using data from the Allen Mouse
   Brain Connectivity Atlas,[237]^36 as there currently is no equivalent
   mesoscale baboon or primate connectivity atlas. All experiments with
   injections in our regions of interest were included. One exception was
   the thalamus. Due to the vast number of experiments with injections
   that included the thalamus, we restricted the thalamus experiments to
   just those that primarily targeted it and manually removed a couple
   experiments that we found to have labeling that extended too far out of
   the region. The mapping between the baboon regions of interest in this
   manuscript and the mouse regions of interest used from the Allen
   Institute are in [238]Table S1 ([239]Table S1). Projections from source
   regions to the other regions of interest were normalized across
   experiments by dividing the projection density by the injection volume.
   For each pair of source regions and target regions, the normalized
   projection density was averaged across all experiments to produce a
   connectivity matrix. These values were then all log normalized to
   produce a better distribution to use as edge weights in a network.
   Network analysis was completed in Python using NetworkX to represent
   the networks.[240]^89 Plotly was used to create network
   visualizations.[241]^90 Community detection was performed using the
   Louvain algorithm as implemented in scikit-network with a resolution
   parameter of 1.[242]^91 The code used for our analysis can be found at
   the following
   link:([243]https://github.com/pderdeyn/neural-cyclic-gene-networks).

Network statistical analysis

   The numbers of cyclic genes per node were compared across communities
   using a one-way ANOVA and then individual t-tests, both as implemented
   in scipy.[244]^92 A mixed effects model was also built to test the
   effect of community membership on number of cyclic genes, using the
   lme4 library in R.[245]^93 Community membership was modeled as a random
   effect on number of cyclic genes. Linear correlation models were fit
   using ordinary least squares, implemented in Python in statsmodels.
   p-values of the linear coefficient were reported.[246]^94
   Heteroscedasticity of data was tested using the Breusch-Pagan test in
   statsmodels. Linear correlation coefficient p-values were tested for
   bias by scrambling the data and fitting a new linear regression model
   1000 times. Monotonicity of data was tested using the Spearman
   correlation in scipy.

Functional enrichment analysis

   Pathway enrichment analysis of cycling transcripts was performed using
   the PANTHER platform [247]http://www.pantherdb.org/. Since the baboon
   gene database is not available in the Panther platform, the transcripts
   in the chimpanzee (Pan Troglodytes) database were selected as the
   reference list. Both chimpanzees (Pan troglodytes) and baboons (Papio
   anubis) share a common ancestor with humans. This shared lineage
   results in a significant degree of genetic homology among these
   species.[248]^95^,[249]^96 Genome sequencing of chimpanzees and bonobos
   revealed a high degree of identity (approximately 98.7%) with human
   orthologous sequences.[250]^95^,[251]^97 Despite their differences in
   cognitive abilities and social behaviors, the genetic similarities
   between baboons and chimpanzees have allowed to utilize the extensively
   annotated chimpanzee gene database for analyzing baboon gene expression
   data.[252]^95^,[253]^97 Data from the PANTHER pathways annotation
   dataset were compared using Fisher's exact test, setting log-fold
   enrichment value of (2) as the threshold (FDR<0.05).

Statistical analysis

   For circadian rhythm insights, transcriptomic time series data was
   analyzed using CircadiOmics and BIO_CYCLE, with a significance cut-off
   at p=0.05. Network analysis employed one-way ANOVA and t-tests in
   scipy, mixed-effects models in R's lme4, and linear models in Python's
   statsmodels. Functional enrichment of cycling transcripts was conducted
   using the PANTHER platform, leveraging the chimpanzee database due to
   its genetic closeness to baboons, with significance determined through
   Fisher's exact test at an FDR<0.05.

Acknowledgments