- Research article
- Open Access
Evolution and transition of expression trajectory during human brain development
BMC Evolutionary Biology volume 20, Article number: 72 (2020)
The remarkable abilities of the human brain are distinctive features that set us apart from other animals. However, our understanding of how the brain has changed in the human lineage remains incomplete, but is essential for understanding cognition, behavior, and brain disorders in humans. Here, we compared the expression trajectory in brain development between humans and rhesus macaques (Macaca mulatta) to explore their divergent transcriptome profiles.
Results showed that brain development could be divided into two stages, with a demarcation date in a range between 25 and 26 postconception weeks (PCW) for humans and 17-23PCWfor rhesus macaques, rather than birth time that have been widely used as a uniform demarcation time of neurodevelopment across species. Dynamic network biomarker (DNB) analysis revealed that the two demarcation dates were transition phases during brain development, after which the brain transcriptome profiles underwent critical transitions characterized by highly fluctuating DNB molecules. We also found that changes between early and later brain developmental stages (as defined by the demarcation points) were substantially greater in the human brain than in the macaque brain. To explore the molecular mechanism underlying prolonged timing during early human brain development, we carried out expression heterochrony tests. Results demonstrated that compared to macaques, more heterochronic genes exhibited neoteny during early human brain development, consistent with the delayed demarcation time in the human lineage, and proving that neoteny in human brain development could be traced to the prenatal period. We further constructed transcriptional networks to explore the profile of early human brain development and identified the hub gene RBFOX1 as playing an important role in regulating early brain development. We also found RBFOX1 evolved rapidly in its non-coding regions, indicating that this gene played an important role in human brain evolution. Our findings provide evidence that RBFOX1 is a likely key hub gene in early human brain development and evolution.
By comparing gene expression profiles between humans and macaques, we found divergent expression trajectories between the two species, which deepens our understanding of the evolution of the human brain.
Our highly developed and distinctive brains, which set humans apart from other mammals, are the product of evolution [1, 2], the mechanism of which has fascinated people for centuries . Based on compelling differences in cognitive and behavioral capacities, but relatively close phylogenetic relationship between humans and non-human primates (NHPs) [4,5,6], recent comparative analyses have provided a novel strategy to study human-specific neurodevelopment [7,8,9]. Increasingly persuasive evidence suggests that brain development is not static but is a continuous process of molecular changes throughout life, including changes in gene expression, glucose metabolism, and synaptic density [10,11,12,13,14]. Previous comparative analyses between humans and NHPs have only offered a snapshot in time [15, 16]. However, it is necessary to compare the whole process of brain development to provide a more objective and comprehensive understanding of human brain evolution.
Earlier research noted that neurodevelopmental timing is impacted by different developmental rates and life history strategies . These differences in neurodevelopmental timing among species, also called heterochrony, have long been considered a crucial impetus for evolution [17,18,19]. Humans have an unusually extended childhood and slow rate of neurodevelopment (known as neoteny) relative to other animals, which is considered a possible mechanism for human brain evolution [18, 20]. While previous studies have primarily focused on heterochronic gene expression during postnatal brain development , the macroscopic layout of the brain is nearly complete at the time of birth . Thus, extending comparative analysis to the prenatal stages is necessary for exploring the features of neurodevelopment.
Currently, it is widely accepted that changes in spatiotemporal gene expression play a critical role in the emergence of the sophisticated human brain, and several attempts have been made to estimate divergent gene expression patterns between humans and NHPs [15, 22, 23]. However, our understanding of how gene expression patterns have changed in the human lineage remains incomplete. With increasing high-quality brain transcriptome data [24,25,26], an unprecedented opportunity to investigate gene expression trajectory in multiple brain regions and different developmental stages among primates has become possible [22, 27]. In this study, we collected large-scale gene expression data from humans and macaques to systematically investigate and compare their divergent gene expression trajectory. We aimed to identify critical states during brain development as well as novel molecular mechanisms underlying human brain evolution.
Figure 1 highlights the strategy used to investigate evolution of gene expression trajectory in humans, including hierarchical clustering and dynamic network analyses to identify demarcation times of brain development in humans and macaques, expression heterochrony analysis to explore the mechanism of neurodevelopmental timing in humans, and differential expression and gene co-expression network analyses to identify key genes in human brain development and evolution.
Human brain transcriptome data were used in this study, including RNA-sequencing (RNA-seq) and microarray data across multiple brain regions downloaded from the Allen Brain Atlas [24, 25] (Table 1; Additional file 1: Table S1 ~ Table S2). These data cover 14 brain regions spanning 27 different developmental ages from 8 postconception weeks (PCW) to postnatal 40 years old (Table 2). In total, 30,881 and 17,280 genes had detectable expression signals in the RNA-seq and microarray data, respectively. We also used microarray data from macaques, which contained five brain regions (22 brain subregions) corresponding to brain regions in humans (Table 1; Additional file 1: Table S3 ~ Table S4) and spanning 8 PCW to postnatal 48 months (Table 2). In total, 15,381 genes exhibited detectable expression signals.
Different developmental trajectories and demarcation times in humans and macaques
We performed hierarchical clustering analysis based on gene expression levels to determine whether transformation exists during brain development. In humans, the RNA-seq and microarray data supported the division of brain development into two stages, with a demarcation time of 25–26 PCW (Fig. 2a; Additional file 2: Figure S1). In macaques, gene expression levels from the microarray data in most brain regions were also clustered into two groups, with a demarcation time of 17 PCW ~ birth (17 -23PCW)(Fig. 2a; Additional file 2: Figure S2). These results suggest that brain development in both humans and macaques could be divided into two major stages, separated by species-specific demarcation points that occurred prior to birth (25–26 PCW in humans and 17–23 PCW in macaques), rather than birth times, which have been widely used as a uniform demarcation time of neurodevelopment across species [28, 29].
Transitional state and critical transitions during brain development based on dynamic network biomarker (DNB) analysis
Having identified the demarcation points of brain development in humans and macaques, we further applied DNB analysis to verify whether the above demarcation points were also transitional states of brain development. Based on nonlinear dynamic theory, biological processes, such as brain development, are not always smooth, but can exhibit dramatic transitions from one state to another [30, 31]. When the process occurs near a critical transition phase, a dominant group of genes/molecules, i.e., DNBs, can drive transition of the dynamic process [30, 32, 33]. We thus performed a genome-wide DNB analysis to identify the transitional states and DNB genes in both humans and macaques.
The DNB results demonstrated that the transitional states of brain development in humans and macaques occurred at around 26 PCW and around 17PCW, respectively (Fig. 2b-c; Additional file 2: Figure S3). After the transitional state, gene expression patterns changed markedly. The transitional states identified by DNB analysis were largely consistent with the demarcation points above (Fig. 2a), thus supporting the robustness of our results. We also obtained the corresponding DNBs, which included 369 DNB genes in humans and 34 DNB genes in macaques (Additional file 1: Table S5). The greater number of DNB genes in humans suggests more dramatic changes between the early and later stages of brain development in humans relative to that in macaques.
Transcriptional profile and cell fate change from early to later stages of brain development
We used differential expression analysis to compare the degree of change in early and later gene expression between humans and macaques. We selected five brain regions (i.e., hippocampus (HIP), striatum (STR), anterior cingulate cortex (ACC), amygdala (AMY), and primary visual cortex (V1C)) that coexist and contain similar sample sizes in the two species. Based on microarray probes, which matched between humans and macaques, we found a larger number of differentially expressed genes (DEGs) (Benjamini-Hochberg FDR < 0.05, fold change [FC] > 1.5) between early and later stages in humans than in macaques (Fig. 3a; Additional file 1: Table S6). These results also suggest more dramatic changes between early and later stages of brain development in humans relative to macaques.
To better understand brain developmental processes in the human lineage, we conducted functional enrichment analysis of the DEGs between early and later stages in humans. Results showed that up-regulated genes in the early stage were mainly involved in cell cycle, DNA packaging, and meiosis (Fig. 3b), whereas, up-regulated genes in the later stage were enriched in synaptic signaling, myelination, and axon establishment (Fig. 3c). Remarkably, the DEG patterns well matched the reported properties of the neurodevelopmental timeline in humans .
Increasing evidence suggests that a cell fate switch from neurons to glial cells is operational in prenatal brains and represents a key process in brain development [36,37,38]. We thus considered whether the demarcation points in human correspond to the transient time of known neuron to glial cell fate switch. We found up-regulated genes in the early stage were predominantly enriched in neuronal genes (Fig. 3d), whereas up-regulated genes in the later stage represented a diversity of cell-type-associated genes, including astrocytes, oligodendrocytes, and neurons (Fig. 3d), which likely reflects the transformation of cell fate switch from neurons to glial cells during the demarcation point. This conclusion is confirmed by the previous single-cell transcription analysis , which reported that neurons developed from neural progenitor cells in early gestational weeks (GW8, GW9, GW10, GW12, GW13, GW16, GW19, GW23), whereas glial cells (oligodendrocyte progenitor cells and astrocytes) differentiated from neural progenitor cells in later weeks (GW26).
Earlier studies have also reported several pathways that govern the neuron to astrocyte cell switch, including the gp130/JAK/STAT and MEK/MAPK pathways … [38, 40]. We found that DEGs between early and later stages were significantly enriched in the MEK/MAPK pathway (P = 4.6e-04, Fisher’s exact test) (Additional file 1: Table S7). Although DNBs were not significantly enriched in the MEK/MAPK pathway (P = 0.15, Fisher’s exact test), 15 DNBs were still involved (Additional file 1: Table S7). This suggests that the MEK/MAPK pathway likely plays an important role in the transformation ratio of cell type from neurons to glial cells during the human demarcation point (25–26 PCW).
Molecular mechanism underlying protracted timing of early human brain development
The different demarcation points (25–26 PCW in humans and 17-23PCWin macaques) identified here reflect prolonged timing of early brain development in the human lineage. Thus, we performed heterochronic analysis to explore the molecular mechanism underlying different timing of early neurodevelopment between humans and macaques.
We again used five brain regions (i.e., HIP, STR, ACC, AMY, and V1C) that coexist in humans and macaques to test heterochronic gene expression. After rigorous quality control (see Methods), we retained 9758 genes with microarray probes well matched between humans and macaques. Among these genes, we selected several that showed both age-related and species-specific differences for each brain region (see Methods; Additional file 1: Table S8). We then sorted these genes into two categories: (i) human acceleration genes, whose expression change was significantly faster during human brain development than that during macaque brain development (Fig. 4a); and, (ii) human neoteny genes, whose expression change was significantly delayed during human brain development than that during macaque brain development (Fig. 4b), as defined in previous study . Compared to macaques, more genes displayed a neotenic pattern in all five human brain regions (Fig. 4c; Additional file 1: Table 9), consistent with the above delayed demarcation point in the human lineage. In addition, the result suggest that neoteny of human brain development could be traced to prenatal period.
Interestingly, the neotenic genes from the five brain regions significantly overlapped (Additional file 1: Table 10), suggesting that neotenic mechanisms among different brain regions are largely convergent. The functions of these neotenic genes were mainly involved in neurodevelopment-related pathways (Fig. 4d), suggesting that more neurodevelopmental genes exhibited neotenic features, which eventually help humans develop a more complex brain.
Co-expression analysis identifies gene network during early brain development in humans
Extended timing of early neurodevelopment in humans is important for brain evolution [2, 18]. We applied weighted gene co-expression network analysis (WGCNA) to further explore the transcriptional profile of early neurodevelopment in humans [41,42,43]. A total of 38 modules related to early human neurodevelopment were identified (see Methods; Fig. 5a; Additional file 1: Table S11). To quantify network reorganization across early and later brain development, we applied modular differential connectivity (MDC), which is the ratio of the average connectivity for any pair of modules sharing genes in the early stage compared to that in the later stage. Among the 38 early modules, five (GCM1–GCM5) showed gain of connectivity compared to later development, with co-regulation enhanced between genes in these modules. In contrast, 21 modules (LCM1–LCM21) showed loss of connectivity and 12 modules (NCM1–NCM12) (31.5%) showed no change in connectivity and were conserved during development (Additional file 2: Figure S4A-S4B).
For modules with gains or losses of connectivity, we ranked them according to their degree of DEG enrichment across brain regions (Fig. 5a) and MDC scores (Additional file 1: Table S11). Module GCM1, which showed a gain of connectivity in the early stage, was identified as the most highly ranked module. The genes in this module were enriched in neurogenesis, neuron projection morphogenesis, and axon development (Fig. 5d). Additionally, GCM1 showed highly significant enrichment for known autism susceptibility markers (P = 5.49E-07; Fisher’s exact test) , and the expression levels of genes in this module were significantly higher during early brain development (Additional file 2: Figure S4C). These results suggest that GCM1 likely plays an important role in early development of the human brain.
Remarkably, most genes in the GCM1 module were located in human-accelerated conserved non-coding sequences (HACNSs) (P = 3.19E-16; Fisher’s exact test)  or in human DNA sequence accelerated regions (HARs) (P = 7.49E-12; Fisher’s exact test)  (Fig. 5c), suggesting that genes in the GCM1 module also likely played an important role in human brain evolution.
We next mapped the genes in GCM1 to single-cell expression data derived from 20,262 prenatal human prefrontal cortex cells that ranged in age from 8 to 26 gestational weeks  (Additional file 2: Figure S5), which represent a broad diversity of cell types including neural progenitor cells, interneurons, astrocytes, oligodendrocyte progenitor cell, microglia and excitatory neurons. The expression patterns of the GCM1 genes closely matched the cell type specific of excitatory neurons (see Methods) (Fig. 5e), confirming that genes in the GCM1 module function through excitatory neurons.
We further reconstructed the network structure of genes within the GCM1 module based on their connectivity and identified 53 hub genes, 36 of which were early stage-specific hub genes  (see Methods, Fig. 6a). Among these genes, RBFOX1(RNA Binding Fox-1 Homolog 1) was of particular interest. RBFOX1 is a highly conserved splicing regulator that displays higher expression in the brain than in other tissues  (Fig. 6b). RBFOX1 is implicated in autism, epilepsy syndromes, and Alzheimer’s disease [49,50,51] and plays an important role in mammalian brain development . Interestingly, RBFOX1 is also a DNB gene, and therefore considered to play an important role in critical transition during brain development .
We next used evolutionary analysis to test if RBFOX1 experienced positive selection in the human lineage. Although the protein coding sequence of RBFOX1 has not changed in humans compared to other primates (Additional file 2: Table S12; see Methods), six HARs were found in the non-coding regions for this gene (Fig. 6c). HARs are non-coding regions conserved across mammals, and which have acquired many sequence changes in humans since their divergence from chimpanzees . Only seven human RefSeq genes from the entire human genome (based on hg18) contain six or more HARs , suggesting that strongly human-specific accelerated evolution occurred recently in the non-coding region of the human RBFOX1 gene. Our findings provide evidence that RBFOX1 is a likely key hub gene in early human brain development and evolution. In addition, RBFOX1 also showed cell specificity to excitatory neurons by single cell transcriptome analysis (Fig. 6d), suggesting that RBFOX1 functions through excitatory neurons.
The remarkable abilities of the human brain set us apart from NHPs. With the advent of large-scale genomic, transcriptomic, and epigenomic data, many genetic underpinnings of the rapid evolution of the human brain have been revealed [54,55,56,57,58]. However, our understanding of how the brain has changed in the human lineage remains incomplete . Based on large-scale transcriptomic and genomic data, the results of the current study provide new insight into the evolution and transition of gene expression trajectory in the human brain.
Firstly, we found that brain development could be divided into two stages in both humans and macaques; more specifically, demarcation times of 25–26 PCW and 17–23 PCW in humans and rhesus macaques, respectively. Further DNB analysis indicated that the demarcation points were nearly the same as the critical transitional states during brain development in humans and macaques. Previous studies on brain development have primarily used birth as the default boundary [28, 29, 59]. However, we suggest that the demarcation points identified here should be considered in the future to minimize biases in studies of brain development.
Secondly, we also found that neoteny of human brain development could be traced to the prenatal period. Previous studies have primarily focused on heterochronic gene expression during postnatal brain development . The macroscopic layout of the brain is nearly complete at the time of birth . Thus, extending comparative analysis to the prenatal stages is necessary to explore the features of neurodevelopment, which is lacking in prior studies. Thus, in this paper, we performed heterochronic analysis across prenatal samples between humans and macaques, and found that more genes displayed a neotenic pattern in humans than in macaques, consistent with the delayed demarcation time in the human lineage, and proving that neoteny in human brain development could be traced to the prenatal period.
Thirdly, we used gene co-expression network analysis to identify transcription profiles in early human brain development and identified that the RBFOX1 gene likely plays an important role in early human brain development and displays positive selection in its non-coding region [50,51,52]. Therefore, we speculated thatRBFOX1 is a likely key hub gene in early human early brain development and evolution. As such, we propose that RBFOX1 should be considered in further neurodevelopmental research.
Finally, we highlighted the importance of excitatory neurons in human brain development and evolution. Over the past few decades, the comparisons of excitatory neurons between humans and NHPs have mainly focused on their differences in morphology and abundance [60, 61]. Further molecular biology research on excitatory neurons is limited. In this paper, we found that the GCM1 module and RBFOX1 gene were related to early human brain development and evolution and were enriched in excitatory neurons. Therefore, studies on excitatory neurons would be promising for exploring human brain evolution.
We note that the study presented here is far from being comprehensive: Firstly, Based on currently available transcriptome data, we identified a demarcation line time-frame of brain development in humans and macaques. The precise demarcation point could not be concluded from existing data but should be explored in future studies.
Secondly, due to the relatively small sample sizes used in the current study, as well as the sparse distribution of samples across ages, we cannot rule out certain important changes in transcriptional profiles during neurodevelopment that may have occurred beyond the sampling range used in this study. For instance, previous studies have reported that during juvenile development in humans (1–8 postnatal years), the cerebral cortex consumes nearly twice the amount of glucose as observed during adulthood and is accompanied by dramatic changes in synaptic density during that developmental window [11, 12]. Thus, the transition state we identified is not absolute, with more saturated samples across different ages required to confirm our conclusions.
Thirdly, comparative analysis between humans and macaques was based on microarray data only, which rely on pre-existing knowledge of RNA sequences; as such, some important genes may be missed.
Finally, due to the lack of prenatal transcriptome data on brain development in hominoids, it is difficult to compare hominoids with humans, which would be valuable when exploring human brain evolution.
Further analyses, including expression data analysis across additional development stages, comparative analysis of RNA-seq data between humans and NHPs, as well as analysis incorporating more hominoids, are needed to expand our results.
In this study, we integrated transcriptomic analysis to reveal the evolution and transition of expression trajectories during human brain development. By comparing gene expression profiles between humans and rhesus macaques, our results provide new insights into the gene expression trajectory of human brain development, which will deepen our understanding on evolution of the human brain.
The normalized gene expression human and rhesus macaque datasets were obtained from the Allen Brain Atlas (http://www.brain-map.org) (Table 1; Table 2) [24, 25]. We used two datasets for humans, including RNA-seq and microarray data, which contained 14 brain regions and 27 developmental stages. The RNA-seq data were summarized to Gencode 10 gene-level reads per kilobase million mapped reads (RPKM), whereas the microarray data were based on the Affymetrix GeneChip Human Exon 1.0 ST Array platform. Several quality control measures were implemented to reduce errors due to spatial artifacts on the chips, technical differences between chips in probe saturation, or other unaccounted for batch effects. Detailed information can be found in the Supplementary Materials of Kang et al. (2011) . For rhesus macaques, we used the microarray dataset based on GeneChip Rhesus Macaque Genome Arrays from Affymetrix. From the 52,865 probe sets included in the macaque microarray data, 12,441 high-confidence probe sets were retained after quality control filtering. Detailed description of the macaque data can be found in Bakken et al. (2016) . The macaque microarray dataset contained five brain regions (22 brain subregions) and nine developmental stages (Table 1; Table 2).
Clustering of genes in each tissue
The human microarray and RNA-seq datasets and rhesus macaque microarray dataset were used to cluster genes for each brain region according to their expression levels. To reduce the influence of technical noise, only genes with expression values of more than 0 in 80% of the available samples were considered expressed. Before clustering, we log2 transformed and then z-transformed the expression levels (normalized the mean to 0 and variance to 1). Using agglomerative hierarchical clustering with the average and complete methods in the flashClust R package , the RNA-seq and microarray data from most human tissues were clustered into two groups, separated with a demarcation point of 25 PCW. The microarray data from the rhesus macaques were also clustered into two groups, with 17 PCW as the demarcation point.
Dynamic network biomarker (DNB) analysis
Based on DNB theoretical analysis [30, 32, 33], we can prove that when a system is near the critical state/transition phase, a dominant group of genes/molecules, i.e., DNBs, can drive transition of the dynamic process. These molecules must satisfy the following three criteria :
Deviation (or fluctuation) for each molecule inside the dominant group (SDd: standard deviation) drastically increases.
Correlation between molecules inside the dominant group (PCCd: Pearson correlation coefficients in absolute values) rapidly increases.
Correlation between molecules inside and outside the dominant group (PCCo: Pearson correlation coefficients in absolute values) rapidly decreases.
The dominant group is considered a DNB and plays an important role in phase transition. A quantification index (CI) considering all three criteria can then be used as the numerical signal of the critical state or transition phase and also for the identification of DNB members/molecule, with the following equation:
where, PCCd is the average Pearson’s correlation coefficient (PCC) between the genes in the dominant group (or DNB) of the same time stage in absolute value; PCCo is the average PCC between the dominant group (or DNB) and others of the same time stage in absolute value; and, SDd is the average standard deviation (SD) of genes in the dominant group (or DNB). The three criteria together construct the composite index (CI). The CI is expected to peak or increase sharply during the measured stages when the system approaches the critical state, thus indicating imminent transition or transition phase of the biological process .
We applied this DNB method to detect the critical points and DNB members in humans and macaques based on the transcriptome data of the primary visual cortex. In each sampling stage, there were 1–4 samples with gene expression profiles. To increase the reliability of the DNB results, the slide window method was incorporated into the DNB model to process data .
Differentially expressed genes (DEGs) between early and later stages of brain development
To remove the potential effect of different high-throughput platforms on gene expression values, we only used microarray data for DEG analysis for both species. Pairwise differential expression was investigated using the edgeR R package . To determine the DEGs between the two developmental times for humans and macaques, the demarcation times were set 25 PCW and 17 PCW, respectively. A nominal significance threshold of Benjamini-Hochberg FDR < 0.05 and fold change [FC] > 1.5 was used to identify DEGs.
For DEGs in humans, we applied g:Profiler (https://biit.cs.ut.ee/gprofiler/)  for functional annotation analysis (GO and KEGG). To assess cell-type specificity in the 14 brain regions of humans, we used genes expressed at least five-fold higher in one cell type than all other cell types (neuron, microglia, astrocyte, oligodendrocyte, endothelial) from brain-based RNA expression data as the cell marker .
Heterochrony analyses with dynamic time warping algorithm (DTW-S)
We combined data from microarray probes of humans and macaques to study heterochronic gene expression in five brain regions (i.e., hippocampus (HIP), striatum (STR), anterior cingulate cortex (ACC), amygdala (AMY) and primary visual cortex (V1C)). We used the “sva” R package  to remove batch effects between microarray datasets of humans and rhesus macaques.
To choose age-related genes, we first used a log2 transformed age scale to ensure a more linear relationship between age and phenotype . We tested the effect of age on the expression levels using polynomial regression models, as described previously . We next tested each gene for expression divergence between humans and rhesus macaques using analysis of covariance  (F-test P < 0.05). Identification of age-related genes and species-specific genes were based on the adjusted r2 criterion. The identification methods have been described previously . Consequently, we selected 892 genes in AMY, 2431 genes in HIP, 1961 genes in STR, 1899 genes in ACC, and 2416 genes in V1C as the test gene set for DTW-S, satisfying the following criteria: (i) significant expression change with age and (ii) significant expression difference between humans and macaques.
The DTW-S algorithm was then used to analyze the data for heterochrony . We defined genes showing significant heterochrony into two categories: (i) human acceleration genes, whose expression changes were significantly faster during human brain development than that during macaque brain development; and, (ii) human neoteny genes, whose expression changes were significantly delayed during human brain development compared with that during macaque brain development. Using those genes that showed significant age-related and species-specific differences, as defined above, we aligned the macaque and human expression trajectories and estimated the time-shift (heterochrony) between humans and macaques, with simulations conducted to estimate the significance of the shifts. We considered genes as ‘significantly heterochronic’ if they showed a shift at P < 0.05. A detailed description of the DTW-S algorithm can be found in Yuan et al. 2011 .
Construction of gene co-expression modules for early human brain development
We constructed multi-tissue co-expression networks that simultaneously captured intra- and inter-tissue gene-gene interactions using the human RNA-seq expression data [42, 70]. Before identifying co-expressed gene modules, we used linear regression to correct sex and brain region covariates in the expression data. To quantify the differences in the transcript network organization between the early and late stages, we employed the modular differential connectivity (MDC) metric . In brief, MDC represents the connectivity ratios of all gene pairs in a module from the early stage to the same gene pairs from the later stage: MDC > 0 indicates a gain of connectivity or enhanced co-regulation between genes in the early stage, whereas MDC < 0 indicates a loss of connectivity or reduced co-regulation between genes in the early stage.
To identify key regulator (driver) genes of the GCM1 module, we applied key driver analysis to the module-based unweighted co-expression networks derived from ARACNE . ARACNE first identified significant interactions between genes in the module based on their mutual information and then removed indirect interactions through data processing inequality (DPI). For each ARACNE-derived unweighted network, we further identified the key regulators by examining the number of N-hop neighborhood nodes for each gene.
Identification of cell type and subtype from single cell data
Single-cell RNA-seq data (accession number GSE104276) were reported in previous study . Transcript counts for each cell were normalized to transcripts per million (TPM), with the TPM values then normalized by log ((TPM/10) + 1) for subsequent analysis . The Seurat package  v1.2.1 implemented in R was applied to identify major cell types among the 2394 single cells from the prefrontal cortex. Only cells that expressed more than 1000 genes were considered, and only genes with normalized expression levels greater than 1 and expressed in at least three single cells were included, which left 20,262 genes across 2344 samples for clustering analysis. After initial clustering, PAX6, NEUROD2, GAD1, PDGFRA, AQP4, and PTPRC were used as markers to identify the major cell types in the brain: i.e., neural progenitor cells, excitatory neurons, interneurons, oligodendrocyte progenitor cells, astrocytes, and microglia, respectively.
Coding sequence evolutionary analysis of RBFOX1
To analyze the evolution of the coding regions of RBFOX1, we obtained the human, chimpanzee, rhesus macaque, marmoset, mouse lemur, mouse, rat, cow, dog, and opossum coding sequences for this gene from Ensembl . The coding sequences were aligned using Prank . Gblocks v0.91b was used to remove poorly aligned regions in the resulting nucleotide sequence alignments . We then used the modified branch-site model A from the PAML package v4.9 to test positive selection of RBFOX1 in the human and primate lineages, respectively . The null hypothesis of the branch test was that all lineages shared the same dN/dS ratio. The alternative hypothesis was that human or primate lineages had a different dN/dS ratio from other lineages, with w0, w1, and w2 representing codons under negative, null, and positive selection, respectively. The Chi-square test was used to calculate the P value P_adjust< 0.05 was considered as significant.
Availability of data and materials
The RNA-seq and microarray dataset of human analysed during the current study are available in the http://www.brainspan.org/static/download.html repository . The microarray dataset of macaque analysed during the current study are available in the http://blueprintnhpatlas.org/static/download repository .
- DNB :
Dynamic network biomarker
- NHPs :
- RNA-seq :
- HIP :
- ACC :
Anterior cingulate cortex
- AMY :
Primary visual cortex
- DEGs :
Differentially expressed genes
- WGCNA :
Gene co-expression network analysis
Modular differential connectivity
Modules that show gain of connectivity
- LCM :
Modules that show loss of connectivity
- NCM :
Modules that show no change of connectivity
- RBFOX1 :
RNA Binding Fox-1 Homolog 1
- HACNSs :
Human-accelerated conserved non-coding sequences
Human DNA sequence accelerated regions
- TPM :
Transcripts per million
Sousa AMM, Meyer KA, Santpere G, Gulden FO, Sestan N. Evolution of the human nervous system function, structure, and development. Cell. 2017;170(2):226–47.
Workman AD, Charvet CJ, Clancy B, Darlington RB, Finlay BL. Modeling transformations of neurodevelopmental sequences across mammalian species. J Neurosci. 2013;33(17):7368.
Bae B-I, Jayaraman D, Walsh CA. Genetic changes shaping the human brain. Dev Cell. 2015;32(4):423–34.
Lieberman P. The evolution of language and thought. J Anthropol Sci. 2016;94:127–46.
Penn DC, Holyoak KJ, Povinelli DJ. Darwin's mistake: explaining the discontinuity between human and nonhuman minds. Behav Brain Sci. 2008;31(2):109–30 discussion 130-178.
Zhang M-L. M-LLAOARWMD-DWYS: Conserved sequences identify the closest living relatives of primates. Zool Res. 2019;40(6):532–40.
Enard W. The molecular basis of human Brain evolution. Curr Biol. 2016;26(20):R1109–17.
Franchini LS, Pollard K. Genomic approaches to studying human-specific developmental traits. Development. 2015;142:3100–12.
Silver DL. Genomic divergence and brain evolution: how regulatory DNA influences development of the cerebral cortex. BioEssays. 2016;38(2):162–71.
Chugani HT, Phelps ME, Mazziotta JC. Positron emission tomography study of human brain functional development. Ann Neurol. 1987;22(4):487–97.
Peter RH. Synaptic density in human frontal cortex — developmental changes and effects of aging. Brain Res. 1979;163(2):195–205.
Sterner K, Weckle AT, Chugani H, Tarca A, Sherwood C, Hof P, Kuzawa C, Boddy A, Abbas A, Raaum R, et al. Dynamic Gene Expression in the Human Cerebral Cortex Distinguishes Children from Adults. PloS one. 2012;7:e37714.
Jacobs B, Chugani HT, Allada V, Chen S, Phelps ME, Pollack DB, Raleigh MJ. Developmental changes in Brain metabolism in sedated rhesus macaques and Vervet monkeys revealed by positron emission tomography. Cereb Cortex. 1995;5(3):222–33.
Ye L-Q, Zhao H, Zhou H-J, Ren X-D, Liu L-L, Otecko NO, Wang Z-B, Yang M-M, Zeng L, Hu X-T, et al. The RNA editome of Macaca mulatta and functional characterization of RNA editing in mitochondria. Sci Bull. 2017;62(12):820–30.
Cáceres M, Lachuer J, Zapala MA, Redmond JC, Kudo L, Geschwind DH, Lockhart DJ, Preuss TM, Barlow C. Elevated gene expression levels distinguish human from non-human primate brains. Proc Natl Acad Sci. 2003;100(22):13030.
Gu J, Gu X. Induced gene expression in human brain after the split from chimpanzee. Trends Genet. 2003;19(2):63–5.
Kim J, Kerr JQ, Min G-S. Molecular heterochrony in the early development of <em>Drosophila</em>. Proc Natl Acad Sci. 2000;97(1):212.
Langer J. The Heterochronic evolution of primate cognitive development. Biol Theory. 2006;1(1):41–3.
Moss EG. Heterochronic genes and the nature of developmental time. Curr Biol. 2007;17(11):R425–34.
Somel M, Franz H, Yan Z, Lorenc A, Guo S, Giger T, Kelso J, Nickel B, Dannemann M, Bahn S, et al. Transcriptional neoteny in the human brain. Proc Natl Acad Sci U S A. 2009;106(14):5743–8.
Stiles J, Jernigan TL. The basics of brain development. Neuropsychol Rev. 2010;20(4):327–48.
Zhu Y, Sousa AMM, Gao T, Skarica M, Li M, Santpere G, Esteller-Cucala P, Juan D, Ferrández-Peral L, Gulden FO, et al. Spatiotemporal transcriptomic divergence across human and macaque brain development. Science. 2018;362(6420):eaat8077.
Khaitovich P, Muetzel B, She X, Lachmann M, Hellmann I, Dietzsch J, Steigele S, Do H-H, Weiss G, Enard W, et al. Regional patterns of gene expression in human and chimpanzee brains. Genome Res. 2004;14(8):1462–73.
Kang HJ, Kawasawa YI, Cheng F, Zhu Y, Xu X, Li M, Sousa AM, Pletikos M, Meyer KA, Sedmak G, et al. Spatio-temporal transcriptome of the human brain. Nature. 2011;478(7370):483–9.
Bakken TE, Miller JA, Ding S-L, Sunkin SM, Smith KA, Ng L, Szafer A, Dalley RA, Royall JJ, Lemon T, et al. A comprehensive transcriptional map of primate brain development. Nature. 2016;535:367.
Sunkin SM, Ng L, Lau C, Dolbeare T, Gilbert TL, Thompson CL, Hawrylycz M, Dang C. Allen Brain atlas: an integrated spatio-temporal portal for exploring the central nervous system. Nucleic Acids Res. 2013;41(Database issue):D996–D1008.
Li M-L, Wu S-H, Zhang J-J, Tian H-Y, Shao Y, Wang Z-B, Irwin DM, Li J-L, Hu X-T, Wu D-D. 547 transcriptomes from 44 brain areas reveal features of the aging brain in non-human primates. Genome Biol. 2019;20(1):258.
Raznahan A, Greenstein D, Lee NR, Clasen LS, Giedd JN. Prenatal growth in humans and postnatal brain maturation into late adolescence. Proc Natl Acad Sci. 2012;109(28):11366.
Somel M, Guo S, Fu N, Yan Z, Hu HY, Xu Y, Yuan Y, Ning Z, Hu Y, Menzel C, et al. MicroRNA, mRNA, and protein expression link development and aging in human and macaque brain. Genome Res. 2010;20(9):1207–18.
Chen L, Liu R, Liu ZP, Li M, Aihara K. Detecting early-warning signals for sudden deterioration of complex diseases by dynamical network biomarkers. Sci Rep. 2012;2:342.
Richard A, Boullu L, Herbach U, Bonnafoux A, Morin V, Vallin E, Guillemin A, Papili Gao N, Gunawan R, Cosette J, et al. Single-cell-based analysis highlights a surge in cell-to-cell molecular variability preceding irreversible commitment in a differentiation process. PLoS Biol. 2016;14(12):e1002585.
Liu X, Chang X, Liu R, Yu X, Chen L, Aihara K. Quantifying critical states of complex diseases using single-sample dynamic network biomarkers. PLoS Comput Biol. 2017;13(7):e1005633.
Liu R, Wang X, Aihara K, Chen L. Early diagnosis of complex diseases by molecular biomarkers, network biomarkers, and dynamical network biomarkers. Med Res Rev. 2014;34(3):455–78.
Cahoy JD, Emery B, Kaushal A, Foo LC, Zamanian JL, Christopherson KS, Xing Y, Lubischer JL, Krieg PA, Krupenko SA, et al. A transcriptome database for astrocytes, neurons, and oligodendrocytes: a new resource for understanding brain development and function. J Neurosci. 2008;28(1):264–78.
Silbereis JC, Pochareddy S, Zhu Y, Li M, Sestan N. The cellular and molecular landscapes of the developing human central nervous system. Neuron. 2016;89(2):248–68.
Shen Q, Wang Y, Dimos JT, Fasano CA, Phoenix TN, Lemischka IR, Ivanova NB, Stifani S, Morrisey EE, Temple S. The timing of cortical neurogenesis is encoded within lineages of individual progenitor cells. Nat Neurosci. 2006;9(6):743–51.
Qian X, Shen Q, Goderie SK, He W, Capela A, Davis AA, Temple S. Timing of CNS cell generation: a programmed sequence of neuron and glial cell production from isolated murine cortical stem cells. Neuron. 2000;28(1):69–80.
Miller FD, Gauthier AS. Timing is everything: making neurons versus glia in the developing cortex. Neuron. 2007;54(3):357–69.
Zhong S, Zhang S, Fan X, Wu Q, Yan L, Dong J, Zhang H, Li L, Sun L, Pan N, et al. A single-cell RNA-seq survey of the developmental landscape of the human prefrontal cortex. Nature. 2018;555(7697):524–8.
Li X, Newbern JM, Wu Y, Morgan-Smith M, Zhong J, Charron J, Snider WD. MEK is a key regulator of Gliogenesis in the developing Brain. Neuron. 2012;75(6):1035–50.
Bagot RC, Cates HM, Purushothaman I, Lorsch ZS, Walker DM, Wang J, Huang X, Schluter OM, Maze I, Pena CJ, et al. Circuit-wide transcriptional profiling reveals Brain region-specific gene networks regulating depression susceptibility. Neuron. 2016;90(5):969–83.
Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC bioinformatics. 2008;9:559.
Parikshak NN, Luo R, Zhang A, Won H, Lowe JK, Chandran V, Horvath S, Geschwind DH. Integrative functional genomic analyses implicate specific molecular pathways and circuits in autism. Cell. 2013;155(5):1008–21.
Meyer KA, Marques-Bonet T, Sestan N. Differential gene expression in the human Brain is associated with conserved, but not accelerated, noncoding sequences. Mol Biol Evol. 2017;34(5):1217–29.
Lindblad-Toh K, Garber M, Zuk O, Lin MF, Parker BJ, Washietl S, Kheradpour P, Ernst J, Jordan G, Mauceli E, et al. A high-resolution map of human evolutionary constraint using 29 mammals. Nature. 2011;478(7370):476–82.
Willsey AJ, Sanders SJ, Li M, Dong S, Tebbenkamp AT, Muhle RA, Reilly SK, Lin L, Fertuzinhos S, Miller JA, et al. Coexpression networks implicate human midfetal deep cortical projection neurons in the pathogenesis of autism. Cell. 2013;155(5):997–1007.
Margolin AA, Nemenman I, Basso K, Wiggins C, Stolovitzky G, Dalla Favera R, Califano A. ARACNE: an algorithm for the reconstruction of gene regulatory networks in a mammalian cellular context. BMC bioinformatics. 2006;7(Suppl 1):S7.
Zerbino DR, Achuthan P, Akanni W, Amode MR, Barrell D, Bhai J, Billis K, Cummins C, Gall A, Girón CG, et al. Ensembl 2018. Nucleic Acids Res. 2018;46(D1):D754–61.
Alkallas R, Fish L, Goodarzi H, Najafabadi HS. Inference of RNA decay rate from transcriptional profiling highlights the regulatory programs of Alzheimer’s disease. Nat Commun. 2017;8(1):909.
Lal D, Reinthaler EM, Altmuller J, Toliat MR, Thiele H, Nurnberg P, Lerche H, Hahn A, Moller RS, Muhle H, et al. RBFOX1 and RBFOX3 mutations in rolandic epilepsy. PLoS One. 2013;8(9):e73323.
Lee JA, Damianov A, Lin CH, Fontes M, Parikshak NN, Anderson ES, Geschwind DH, Black DL, Martin KC. Cytoplasmic Rbfox1 regulates the expression of synaptic and autism-related genes. Neuron. 2016;89(1):113–28.
Gehman LT, Stoilov P, Maguire J, Damianov A, Lin CH, Shiue L, Ares M Jr, Mody I, Black DL. The splicing regulator Rbfox1 (A2BP1) controls neuronal excitation in the mammalian brain. Nat Genet. 2011;43(7):706–11.
Kamm GB, Pisciottano F, Kliger R, Franchini LF. The developmental brain gene NPAS3 contains the largest number of accelerated regulatory sequences in the human genome. Mol Biol Evol. 2013;30(5):1088–102.
He Z, Han D, Efimova O, Guijarro P, Yu Q, Oleksiak A, Jiang S, Anokhin K, Velichkovsky B, Grünewald S, et al. Comprehensive transcriptome analysis of neocortical layers in humans, chimpanzees and macaques. Nat Neurosci. 2017;20:886.
Kronenberg ZN, Fiddes IT, Gordon D, Murali S, Cantsilieris S, Meyerson OS, Underwood JG, Nelson BJ, Chaisson MJP, Dougherty ML, et al. High-resolution comparative analysis of great ape genomes. Science. 2018;360:eaar6343.
Sousa AMM, Zhu Y, Raghanti MA, Kitchen RR, Onorati M, Tebbenkamp ATN, Stutz B, Meyer KA, Li M, Kawasawa YI, et al. Molecular and cellular reorganization of neural circuits in the human lineage. Science. 2017;358(6366):1027–32.
Vermunt MW, Tan SC, Castelijns B, Geeven G, Reinink P, de Bruijn E, Kondova I, Persengiev S, Netherlands Brain B, Bontrop R, et al. Epigenomic annotation of gene regulatory alterations during evolution of the primate brain. Nat Neurosci. 2016;19:494.
Xu C, Li Q, Efimova O, He L, Tatsumoto S, Stepanova V, Oishi T, Udono T, Yamaguchi K, Shigenobu S, et al. Human-specific features of spatial gene expression and regulation in eight brain regions. Genome Res. 2018;28(8):1097–110.
Bakken TE, Miller JA, Luo R, Bernard A, Bennett JL, Lee C-K, Bertagnolli D, Parikshak NN, Smith KA, Sunkin SM, et al. Spatiotemporal dynamics of the postnatal developing primate brain transcriptome. Hum Mol Genet. 2015;24(15):4327–39.
Bianchi S, Stimpson CD, Bauernfeind AL, Schapiro SJ, Baze WB, McArthur MJ, Bronson E, Hopkins WD, Semendeferi K, Jacobs B, et al. Dendritic morphology of pyramidal neurons in the chimpanzee neocortex: regional specializations and comparison to humans. Cereb Cortex. 2013;23(10):2429–36.
Elston G, Benavides-Piccione R, Elston A, Manger P, Defelipe J. Pyramidal cells in prefrontal cortex of Primates: marked differences in neuronal structure among species. Front Neuroanat. 2011;5:2.
Langfelder P, Horvath S. Fast R functions for robust correlations and hierarchical clustering. J Stat Softw. 2012;46(11):i11.
Guttman M, Garber M, Levin JZ, Donaghey J, Robinson J, Adiconis X, Fan L, Koziol MJ, Gnirke A, Nusbaum C, et al. Ab initio reconstruction of cell type-specific transcriptomes in mouse reveals the conserved multi-exonic structure of lincRNAs. Nat Biotechnol. 2010;28(5):503–10.
Reimand J, Arak T, Adler P, Kolberg L, Reisberg S, Peterson H, Vilo J. g:profiler-a web server for functional interpretation of gene lists (2016 update). Nucleic Acids Res. 2016;44(W1):W83–9.
Zhang Y, Chen K, Sloan SA, Bennett ML, Scholze AR, O'Keeffe S, Phatnani HP, Guarnieri P, Caneda C, Ruderisch N, et al. An RNA-sequencing transcriptome and splicing database of glia, neurons, and vascular cells of the cerebral cortex. J Neurosci. 2014;34(36):11929–47.
Leek JT, Johnson WE, Parker HS, Jaffe AE, Storey JD. The sva package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics. 2012;28(6):882–3.
Clancy B, Darlington RB, Finlay BL. Translating developmental time across mammalian species. Neuroscience. 2001;105(1):7–17.
Yuan Y, Chen Y-PP, Ni S, Xu AG, Tang L, Vingron M, Somel M, Khaitovich P. Development and application of a modified dynamic time warping algorithm (DTW-S) to analyses of primate brain expression time series. BMC bioinformatics. 2011;12:347.
Faraway JJ. Practical regression and ANOVA using R; 2002.
Zhang B, Gaiteri C, Bodea LG, Wang Z, McElwee J, Podtelezhnikov AA, Zhang C, Xie T, Tran L, Dobrin R, et al. Integrated systems approach identifies genetic nodes and networks in late-onset Alzheimer's disease. Cell. 2013;153(3):707–20.
McKenzie AT, Katsyv I, Song WM, Wang M, Zhang B. DGCA: a comprehensive R package for differential gene correlation analysis. BMC Syst Biol. 2016;10(1):106.
Butler A, Hoffman P, Smibert P, Papalexi E, Satija R. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat Biotechnol. 2018;36(5):411–20.
Loytynoja A. Phylogeny-aware alignment with PRANK. Methods Mol Biol. 2014;1079:155–70.
Castresana J. Selection of conserved blocks from multiple alignments for their use in phylogenetic analysis. Mol Biol Evol. 2000;17(4):540–52.
Yang Z, dos Reis M. Statistical properties of the branch-site test of positive selection. Mol Biol Evol. 2011;28(3):1217–28.
This work was supported by the Animal Branch of the Germplasm Bank of Wild Species, Chinese Academy of Sciences (the Large Research Infrastructure Funding).
This work was supported by the Strategic Priority Research Program of the Chinese Academy of Sciences (XDB13000000, XDB38040400), National Key R&D Program of China (2017YFA0505500), National Natural Science Foundation of China (31671325, 31822048, 31771476, 31930022), and the Bureau of Science and Technology of Yunnan Province (2019FI010). The funding bodies played no role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript. Y.L. was supported by the Young Academic and Technical Leader Raising Foundation of Yunnan Province.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
About this article
Cite this article
Li, ML., Tang, H., Shao, Y. et al. Evolution and transition of expression trajectory during human brain development. BMC Evol Biol 20, 72 (2020). https://doi.org/10.1186/s12862-020-01633-4
- Expression trajectory
- Brain evolution