Figure3: Expression dynamics of CGs-CNs in developing human brain.

(a) Binary activity heatmap of CGs modules over cell clusters of the developing human brain, cell clusters as annotated by Nowakowski et al. for cluster description refer table S3. (b) Heatmap showing gene expression of CGs over cell clusters, which showed significant enrichment score in the co-expression analysis. Cell clusters as annotated by Nowakoski et al. for cluster description refer table S3. (c, d, e & f) Network plots showing CGs and its targets in the developing human brain, the edges are placed from highest (red) to lowest (yellow) importance score of target genes. GeneMANIA reported co-expression gene pairs are highlighted in blue color. (c) BHLHE41, (d) BHLHE40, (e) RORB and (f) NR1D2 only top15 target genes in the BHLHE40/41 module are shown, target genes with known diurnal expression are underlined. (g) Binary activity heatmap of upstream TFs co-expressing with BHLHE41 over cell clusters of the developing human brain, cell clusters as annotated by Nowakowski et al. for cluster description refer table S3. (h) Heatmap showing gene expression of upstream TFs co-expressing with BHLHE41, which showed significant enrichment score in the co-expression analysis. Cell clusters as annotated by Nowakowski et al. for cluster description refer table S3.

Knockdown of NR1D2 in adult murine neural stem cells has been shown to affect neural differentiation⁽³⁷⁾. The topmost target genes of the NR1D2 module included SLC35F2, ELF1, PON2, KPNA2 and EEF1B2 (Figure 3f). These genes are involved in vital neuronal functions like SLC35F2 in transport across blood brain barrier⁽³⁸⁾, ELF1 in neurite growth⁽³⁹⁾ and PON2 imparts neuroprotection by scavenging superoxides^(38,40). The enrichment of four CGs-CNs in diverse cell types of the developing human brain prompted us to look for upstream regulators of CGs. Our search for upstream TFs co-expressing with CGs revealed only one significant upstream module associated with BHLHE41, (Figure 3g). We found four TFs namely TAL1, IRF8, IKZF1 and ZBTB24 which are documented to have a regulatory role in microglia activation in response to injury⁽⁴¹⁾, aging⁽⁴²⁾ with TAL1 acting as a specific marker of human microglia⁽⁴³⁾.

Expression dynamics of CGs-CNs during astrocyte genesis

Previously, considered to be the support system of neurons in the brain, glial cells are being recognized as vital players in brain function. By promoting neuronal survival and modulating synapse function these cells actively participate in neuronal signal transduction. In particular, the participation of astrocytes in circadian regulation has been recognized recently⁽⁴⁴⁾. The various functions of astrocytes like modulation of intracellular calcium levels or neurotransmitter release occur rhythmically. A study⁽⁴⁵⁾ has shown that astrocyte specific genetic deletion of Bmal1 and behavioral disruption of rhythms in mice induces astrocyte activation. This effect is augmented when Bmal1 is knocked out in both neurons and astrocytes. These studies point to an important role for astrocytes in maintaining circadian function in the brain.

We analyzed the total transcriptome of human fetal and adult brain-derived astrocytes⁽⁴⁶⁾ . In contrast to neurons, astrocytes showed statistically significant upregulation of CGs namely CLOCK, CRY2, PER3, RORA, RORB, HLF, NR1D1, BHLHE40 and BHLHE41 upon maturation (Figure 4a). To strengthen this observation, we analyzed single-cell transcriptomic data of astrocytes derived from cortical spheroids⁽⁴⁷⁾ . Similar to neuronal differentiation analysis, we first constructed pseudotime trajectory of astrocyte differentiation and classified cell clusters using cellular markers of early and late stages of astrocyte development. Overall, we found an increased expression of all CGs except RORB and RORC in differentiated astrocytes (Figure 4b). This was in concordance with the trend seen in total RNA-seq of fetal and adult brain-derived astrocytes.

Figure 4: Expression dynamics of CGs-CNs during astrocyte genesis

a) Heatmap showing CGs expression in fetal and adult brain derived astrocytes profiled by total RNA-seq. TOP2A and AQP4 represent marker genes associated with early and late stages of astrocyte development, respectively. Asterisks indicate the differential expression with statistical significance calculated using DEseq2.b) Heatmap showing CGs (row) expression over pseudotime ordered cells (column) derived from astrocytic spheroid, M denoting marker genes, S and NS depicting genes showing statistically significant and non-significant pseudotime expression dynamics respectively. Top annotation denotes cell clusters assigned by TSCAN, bottom annotation denotes TSCAN calculated pseudotime in ascending order.

c) Binary activity heatmap of CGs modules over pseudotime aligned trajectory of astrocyte genesis. Top annotation denotes cell clusters assigned by TSCAN, bottom annotation denotes TSCAN calculated pseudotime in ascending order. d) Network plot of BHLHE41 module in astrocytes, the edges are placed from highest (red) to lowest (yellow) importance score of target genes. GeneMANIA reported co-expression gene pairs are highlighted in blue color, only top15 target genes are shown

The co-expression analysis of astrocyte scRNA-Seq revealed the BHLHE41 associated module with 97genes to be the prominent module that showed enrichment in the later stages of astrocyte genesis (Figure 4c). BHLHE41 acts as a transcriptional repressor in the circadian regulation of gene expression. Previous studies have reported the association of mutations in BHLHE41 with short sleep timing in both humans and experimental models^(48,49). Also, it shows high abundance in astrocytes as compared to neurons⁽⁵⁰⁾. Few of the topmost target genes in the BHLHE41 module included CLU, PSAP, BCL6, CRYAB, TPP1, TUBA4A and ADD3 (Figure 4d). The co-expression of a few genes like ITPKB, SSPN and CRYAB have been previously reported in the GeneMANIA. Additionally, we report various novel genes, specifically GFAP, a well-known marker of astrocytes and also shown to be associated with the etiology of Alexander disease⁽⁵¹⁾. We also observed CD44 co-expression with BHLHE41, CD44 acts as a vital adhesion protein during astrocyte development⁽⁵²⁾ and an increase in CD44 positive astrocytes have been shown in Alzheimer’s brain⁽⁵³⁾. Our approach highlights the potential regulation of these vital genes involved in astrocyte function by BHLHE41.

Discussion

We combined pseudotime/co-expression analysis and explored potential co-regulatory modules of CGs in the developing brain. Our workflow reveals CGs-CNs with cell type and cell state identity, a feat not possible from bulk RNA-seq based co-expression analysis. Also, high gene drop-out rates pertaining to scRNA-seq greatly impact the pseudotime analysis of low expressing genes like CGs. Therefore, we looked at co-expression modules that can accurately reflect variations in gene expression trends related to underlying biological processes (neurogenesis).

Our analysis of the murine and human embryonic brain scRNA-seq datasets reveals CGs-CNs of RORB, whose expression and role in cortical development is well studied in murine models. Further, many of the identified co-expression pairs are also reported in GeneMANIA with few of the target genes showing oscillatory expression in the adult brain⁽⁵⁴⁾. These validations highlight the biological relevance of the CGs-CNs captured through our analysis.

Many recent studies reveal the regulatory role of non-neuronal cell types in brain physiology. Our analysis also identified four CGs-CNs in the non-neuronal population of the developing human brain. Among the four modules, BHLHE41 strongly co-express with genes having a crucial function in microglia such as P2RY12 and CLEC7A. Further, our co-expression analysis for TFs upstream of CGs reveal BHLHE41 as a target of microglia associated TFs -TAL1, IRF8, IKZF1 and ZBTB24. In addition, astrocytes derived from cortical spheroids show high expression of the BHLHE41 module in mature astrocytes. Above results clearly favour BHLHE41 as a prioritized module for future biological validation in the context of glial biology.

Our meta-analysis uncovers CGs-CNs in rare cell types of the developing brain- radial glia, microglia, and astrocytes. Similar analysis workflow on scRNA-seq datasets on pre-sorted rare neuronal populations, with high sequencing depth in future can reveal more CGs-CNs in the brain. Further, a rapid increase in public availability of human scRNA-Seq dataset can be utilized to construct CGs-CN map of all cell types in human body. Additively, integrating emerging single cell chromatin accessibility maps with scRNA-seq will improve confidence of co-expression analysis in general. We also want to underlie that through this analysis we only prioritized CGs based on pseudotime and co-expression analysis in developing brain cells, understanding their functional consequence need experimental validations. Thus, studies involving genetic perturbation of these CGs in organoid models of human brain development is needed to fully appreciate the intricate role of circadian transcription factors in human brain development. Considering the known dysregulation of clock genes in psychiatric disorders it will also be interesting to extend similar analysis into scRNA-seq datasets derived from diseased conditions.

Methods

Pseudotime analysis

TSCAN was used for pseudotime analysis (https://github.com/zji90/TSCAN). Each single-cell gene expression matrix was subjected to TSCAN default workflow (Figure 1a), the workflow included preprocessing, dimensionality reduction followed by clustering. The genes with zero count were removed and scRNA-seq dropouts were handled by clustering genes having similar expression patterns using euclidean distance and complete linkage. The above pre-processing method was applied to all the datasets. Next, PCA was used for dimension reduction followed by clustering, yielding a minimum spanning tree. The identification of clusters was performed using marker genes. p-value and FDR (False discover rate) were calculated on the TSCAN ordered cells. TSCAN assigned cell clusters, were inspected for a panel of marker genes from the publications of the respective datasets, while a single marker gene was used for visualization purpose. If required, TSCAN cluster ordering was further manually corrected to place neural stem cell clusters as the starting point of the pseudotime. The TSCAN employs generalized additive model (GAM) for fitting the individual gene expression dynamics over pseudotime. GAM fitting was compared to null-model by assuming gene expression was constant over pseudotime. P-value was computed over the above comparison using likelihood ratio test and then converted to false discovery rate (FDR) using Benjamini-Hochberg Procedure. An FDR of <0.05, was used as a cut-off to score for genes with statistically significant expression dynamics over calculated pseudotime. The accession number of the datasets are provided in the supplementary table 2.

Co-expression analysis

Two steps in the pySCENIC (V0.10.3)⁽¹¹⁾and SCENIC (Single-Cell rEgulatory Network Inference and Clustering)⁽¹²⁾ pipeline were implemented in our co-expression analysis. First the co-expressing gene modules were identified using GRNBoost (https://github.com/aertslab/GRNBoost). Input for GRNBoost included a list of annotated transcription factors (https://github.com/aertslab/pySCENIC/tree/master/resources) and scRNA-seq expression matrix. In return GRNBoost outputs adjacency matrix with TFs, target genes and importance score. All the CGs (TFs) and its associated targets genes with importance scores greater than 99th percentile were subsetted as CGs-CNs.

Next AUCell (V1.11.0) was used for scoring the activity of above CGs-CNs. AUCell uses “Area under curve” to score the activity of a given geneset (CGs-CNs) within expressed genes of scRNA-seq dataset. AUCell was implemented using R (4.0.2) using default parameters. The CGs-CNs which crossed the default threshold of AUCell geneset activity parameters are shown as binary CGs-CNs activity heatmap.

ScRNA-seq dataset description

Mouse Cortex

Yuzwa et al. (GSE107122) performed drop-seq of CD1 mouse embryos at E11.5, 13.5, 15.5 and 17.5 stage. For E15.5 embryonic stage, 13 embryos were used and their brains dissected and dissociated for Drop-seq droplet collection. Following droplet collection single-cell libraries were prepared and sequenced on Illumina NextSeq500 platform.

 Mouse Subventricular zone (SVZ)

Babodilla et al (GSE67833) performed microdissection of SVZ from C57BL/6 adult mice. GLAST⁺ Prominin1⁺ and PSA-NCAM⁺ were used as marker to sort neural stem cells and neuroblasts from SVZ region. scRNA-seq libraries were prepared using Smart-seq2 and sequenced on HiSeq2000 platform. Trimmed reads were mapped to the mouse genome (ENSEMBL Release 78) using STAR_2.4.0g followed by read quantification as FPKM.

Mouse Subgranular zone (SGZ)

Shin et al. (GSE71485) used transgenic mice expressing nuclear localized CFP under the control of nestin (Nes-CFP^nuc). The CFP containing cells marking quiescent-NSCs and their progenies were sorted under a fluorescent microscope using glass pipette. The single-cell libraries were prepared by SMART protocol and sequenced by HiSeq2500 followed by read mapping and quantification to give TPM matrix.

Developing human fetal cortex

Nowakowski et al. performed stereotaxic based microdissection of brain regions like prefrontal cortex (PFC), primary visual cortex (V1) and MGE, the brain region aged from 5-37pcw. The single-cell suspension was prepared using FluidigmC1 followed by sequencing of single-cell libraries on the Illumina Hiseq2500 platform. The sequencing reads were aligned to GRCh38 and quantified to reveal transcript count as CPM.

Astrocyte spheroid

Astrocyte lineage cells were purified from cortical spheroids derived from human iPSC line by Sloan et al (GSE99951). The immunopanning method was used to isolate astrocytes using HepaCAM as an astrocyte marker. The sequencing libraries were prepared from astrocytes derived using immunopan method at 100, 130, 175 and 450 days of in vitro differentiation. Single-cell libraries were prepared using SMART-seq protocol and sequenced on NextSeq500 sequencing platform. The raw reads were aligned to hg19 human genome followed by transcript quantification as FPKM values.

Total RNA-seq of human astrocytes

Zhang et al. (GSE73721) performed RNA-sequencing of astrocytes derived from juvenile (8-18years old) and adult (21-63 years old) human brain. The astrocytes from the brain tissue were purified in culture using immunopanning. Briefly, single-cell suspension of donor tissues was subjected to dishes coated with antibodies against various cell types, Hepa-CAM was used as a marker for astrocytes. The cDNA libraries from the two samples were sequenced on illumina NextSeq sequencer as 150bp paired end reads. The RNA-seq reads were mapped using TopHat2 to hg19 as the reference genome and expression level was estimated as FPKM.

Acknowledgements The authors acknowledge Dr. Beena Pillai and Dr. Souvik Maiti for their help in writing the manuscript and discussions. We are grateful to them for allowing us to communicate the manuscript. We also acknowledge the high-performance-computing (HPC) facility of CSIR-IGIB.

Author contributionsSR and SS contributed equally. SR and SS conceptualized the work, SR and AHA performed the analysis, SR and SS prepared the manuscript.

Supplementary tables

Table S1. List of core clock genes analyzed in this study.

Table S2. List of datasets used in the analysis.

TableS3. Identities of cell clusters in developing human brain (derived from Nowakowski et al. 2017).

References