TempShift Reveals the Sequential Development of Human Neocortex and Skewed Developmental Timing of Down Syndrome Brains

Development is a complex process involving precise regulation. Developmental regulation may vary in tissues and individuals, and is often altered in disorders. Currently, the regulation of developmental timing across neocortical areas and developmental changes in Down syndrome (DS) brains remain unclear. The changes in regulation are often accompanied by changes in the gene expression trajectories, which can be divided into two scenarios: (1) changes of gene expression trajectory shape that reflect changes in cell type composition or altered molecular machinery; (2) temporal shift of gene expression trajectories that indicate different regulation of developmental timing. Therefore, we developed an R package TempShift to separates these two scenarios and demonstrated that TempShift can distinguish temporal shift from different shape (DiffShape) of expression trajectories, and can accurately estimate the time difference between multiple trajectories. We applied TempShift to identify sequential gene expression across 11 neocortical areas, which suggested sequential occurrence of synapse formation and axon guidance, as well as reconstructed interneuron migration pathways within neocortex. Comparison between healthy and DS brains revealed increased microglia, shortened neuronal migration process, and delayed synaptogenesis and myelination in DS. These applications also demonstrate the potential of TempShift in understanding gene expression temporal dynamics during different biological processes.


Introduction
Gene expression has been observed to be dynamically regulated across development in many organisms. The temporal dynamics are often associated with the occurrence of developmental processes, such as generating new cells, responding to internal or external signals, and so on. While some of the developmental processes are universal, some others are distinct to species, tissues or disorders, which are accompanied by different temporal gene expression patterns. The differences during development can be categorized into two scenarios: (1) different shape (DiffShape) of gene expression trajectories may reflect increase or decrease in certain cell types, and disrupted or altered molecular machinery of developmental processes; (2) temporal shift of gene expression trajectories reflect different regulation of developmental timing.
Development of the cerebral neocortex is an example. The neocortex is organized into structurally similar but functionally distinct areas. Development of neocortical areas involves in common processes, such as neurogenesis, neuron differentiation, synaptogenesis, myelination, and so on, but the maturation rate of distinct neocortical areas has been found different. Furthermore, these neurodevelopmental disorders have been observed impaired in neurodevelopmental disorders, such as Down syndrome (DS). DS, also known as trisomy 21, is a genetic disorder caused by the presence of all or part of a third copy of chromosome 21. As the most common neurodevelopmental disorder, the estimated prevalence of DS is as high as 13.65 per 10,000 live births [1]. DS is typically characterized by delays of physical growth and intellectual disabilities. Morphometric and cellular studies in human and functional studies in mouse models have indicated defects in neurogenesis, neuronal differentiation, synaptic plasticity and myelination [2][3][4]. However, understanding the regulation of developmental timing across neocortical areas and defects in DS remains challenging.
The time-series design of transcriptome analyses provides a unique opportunity to detect the two scenarios of developmental changes in a high-throughput way. To date, multiple approaches have been developed for differential expression analysis, clustering and alignment of time-series data [5]. A Gaussian process-based model, TempShift (Appendix A) has been proposed to distinguish the above two scenarios in our previous study on spatiotemporal transcriptomic divergence across human and macaque brain development [6]. TempShift divided the genes into three categories: DiffShape genes-these genes follow different shapes of expression trajectories under different conditions, shift genes-these genes express at different times without changing the expression trajectories over time under different conditions, and no-shift genes-these genes express at the same time and follow the same expression trajectories under different conditions ( Figure 1). TempShift used Gaussian process [7] to model the gene expression trajectories so that it has the following advantages [8,9]: (1) it explicitly addresses the dependencies between consecutive measurements and thus can deal with non-matched time-point sampling in different groups; (2) it can handle an arbitrary number of replicates; (3) it provides a statistical framework to distinguish DiffShape genes from shift genes and no-shift genes, and to infer the time shift between groups; (4) it can be used to analyze two or more groups. Here, we programmed TempShift into an R package. The performance of TempShift was tested using different types of simulated data first. Then, we applied TempShift to human brain transcriptome data including 11 neocortical areas, revealing the gradual development of the human neocortex. Further application to a DS dataset demonstrated the ability of TempShift in handling data with unmatched age and revealed increased microglia and skewed developmental timing of neuronal and oligodendrocyte development in DS.

TempShift: A Statistical Model to Detect Global Temporal Shift between Multiple Time Series
TempShift is a statistical model designed for two goals: (1) distinguishing DiffShape, shift and no-shift genes, and (2) for shift genes, estimating the time shift between groups. To achieve the above goals, it builds three temporal expression models based on Gaussian process regression: independence model, shift model and no-shift model (Figure 1; Appendix A). The DiffShape genes are first separated from genes with the same shape of trajectories (shift and no-shift genes) by comparing the independence model with the shift model, and the shift genes are then selected by comparing the shift model with the no-shift model. These comparisons are based on the log likelihood ratio of the shift model versus the independence model (LLR shape ) and the log likelihood ratio of the shift model versus the no-shift model (LLR shi f t ), respectively (Appendix A). Genes with LLR shape smaller than a threshold Λ shape are selected as DiffShape genes; genes with LLR shape greater than Λ shape and LLR shi f t greater and Λ shi f t are selected as shift genes, and the rest genes are no-shift genes. The time shift (∆t) with respect to the reference group is estimated in the shift model.

Gaussian Process
For group i , the time vectors x i consist of 100 points randomly sampled from the uniform distribution between 5 and 15. For each model, we generated 100 genes with rmvnorm function in the mvtnorm R package [10], using their corresponding covariance matrix. The length of each gene expression vector is 100, corresponding to 100 time points of expression. The following parameters were used in the Gaussian kernel: the amplitude σ f = 5 and the length scale l = 3. The σ f determines the average distance of the function away from its mean. For two-group data, ∆t between two groups were fixed to 2. For three-group data, the first group was considered as the reference group, and the time shifts of the other two groups with respect to the reference group were randomly sampled from Gaussian distribution N (2, 1) or N (−2, 1).

Periodic Data
The sine function sin(·) is a well-known periodic function, so we can simulate periodic data based on the sine function. For group i, the time vectors x i consists of 100 points randomly sampled from the uniform distribution between 5 and 15. For both the noshift model and the shift model, we simulated 100 genes each. The length of each gene expression vector is 100, corresponding to 100 time points of expression. Two-group and three-group data were simulated. For the no-shift model, the gene expression of each group was modeled as y i = sin 2πx i 10 + i , with each element of the error vector i following independent and identically distributed N (0, 0.3). For the shift model, the first group (i = 1) was considered as the reference group, and was modeled by the same function as the no-shift model. For group i (i = 1), the gene expression was modeled as is the time shift between group i and the reference group; → 1 n i is a vector of all ones with length of n i , the sample size of group i; i is the error vector with all elements following independent and identically distributed Gaussian distribution N (0, 0.3).

Polynomial Data
For each group, the time vectors x i were generated in the same way as in the Gaussian process and periodic data. Two-group and three-group data were simulated. For the noshift model, gene expression was modeled as y i = a + b × (x i − 10) + c × (x i − 10) 2 + ε i , for group i. The elements of i follow independent and identically distributed N (0, 3). For the shift model, gene expression of the reference group was simulated by the same formula as in the no-shift model, and gene expression of other groups (i = 1) was generated by shift between group i and the reference group, and i is the error vector with all elements following independent and identically distributed Gaussian distribution N (0, 3).

Identifying Developmental Sequences in the Human Neocortex
Normalized human brain microarray data were downloaded from GSE25219 [11]. Samples from 11 neocortical areas were used for analysis. The samples used in the analysis were listed in Supplementary Table S1. Genes with temporal dynamics were pre-selected by cubic regression gene expression ∼ 1 + age + age 2 + age 3 , where age = log 2 (postconceptual days). A total of 4153 genes with variance explained by the cubic model R 2 > 0.5 were selected for further analysis by TempShift. The R 2 of cubic model measures the relationship between gene expression and developmental time. The larger the R 2 , the greater the change in gene expression with developmental time. Therefore, we selected the genes with R 2 > 0.5 as input temporal dynamics genes to the downstream analysis with TempShift. Next, temporally shifted genes were selected with LLR threshold Λ shape > 50 and Λ shi f t > 50. Enrichment analysis was performed by PANTHER Overrepresentation Test (release 20160321) (http://pantherdb.org/, accessed on 18 April 2016 ) [12].

Gene Selection
Microarray data of DS samples were downloaded from GSE59630 [13]. DFC samples were used and combined with the human DFC samples from GSE25219 in the above section. No obvious batch effect was observed between two data sets, as the same platform and normalization methods were applied. All 17,542 genes were analyzed by TempShift. DiffShape genes were selected with LLR shape below two standard deviations below the mean and shift genes were selected with LLR shape above mean and LLR shi f t above 10.

Cell Type Enrichment
Single cell RNA-seq of human neocortex were downloaded from GSE67835 [14]. Log2-RPM (reads per million mapped reads) were used to measure expression. To assign cell types, we first reduced the dimension of the data by tSNE and then performed k-means clustering. The cell type of a cluster of cells is determined by the expression of cell-type markers. We identified two clusters of prenatal single cells, representing progenitors and neurons, respectively, and six clusters of adult single cells, including pyramidal neuron, interneurons (two clusters), oligodendrocytes, astrocytes and microglia. For each gene, cell type enrichment was calculated by one-way ANOVA followed by post hoc Tukey's honest significant difference (HSD) test. A gene is enriched in cell type A, if ANOVA p < 0.01 and Tukey's HSD test p < 0.01 for comparison between cell type A and at least 6 other cell types. The enrichment of cell type-enriched in DiffShape genes and shift genes was tested by Fisher's exact test.

Testing TempShift by Simulation
To test the performance of TempShift [6] (https://github.com/YingZhuLab/TempShift, accessed on 31 May 2023) implemented in R, we applied it to different types of simulated data with known time shift, including Gaussian process data, periodic data, and polynomial data. The specific principles of data simulation are described in Section 2.

Gaussian Process Data
We first tested the performance of our model on two-group and three-group simulated data generated from Gaussian process. We generated 100 genes following the independence model (  In both the two-group and the three-group data, we observed clear separation between genes following the independence model and genes following the no-shift or shift model by LLR shape and further separation between shift genes and no-shift gene by LLR shi f t (Figure 2d,f). Therefore, based on LLR shape and LLR shi f t , we can achieve the purpose of screening out DiffShape, no-shift and shift genes. Furthermore, our model accurately predicted the time shift in the shift model, with small mean squared prediction error (MSPE; two-group MSPE = 0.0061; three group MSPE = 0.0039 for group 2 and MSPE = 0.0027 for group 3; Figure 2e,g). The mean squared prediction error measures the expected squared distance between our estimated time shift and the true time shift.
Next, we analyze the robustness of the TempShift to noise. Residual variance refers to the variance in a model that cannot be explained by the variables in the model. The higher the residual variance of a model, the higher the noise level. In the above two-group data with σ = 0.5, the residual variance constituted 0.3-65.2% total variance in the noshift model (6.4 ± 9.8%) and 0.4-30.5% (5.7 ± 5.5%) in the shift model, where TempShift successfully distinguished all shift genes from no-shift genes under this noise level. To further explore the effects of noise level, we then increased the residual variance by setting σ = 1, which resulted in the residual variance constituting 2.1-67.2% total variance in the no-shift model (16.7 ± 14.8%) and 1.7-76.4% (11.9 ± 10.8%) in the shift model. In this case, TempShift is also able to distinguish most temporally-shifted genes from no-shift genes, except one gene with 76.4% noise (Figure 2h,i ). In general, LLR shi f t decreases and the prediction error increases with the level of noise (Figure 2j,k), and TempShift is robust to noise.

Periodic Data
We then assessed whether TempShift could be applied to data generated from other models. TempShift was first applied to simulated periodic data, the expression changes associated with periodic processes, such as cell cycle or circadian rhythm. The residual variance constituted approximately 15% of the total variance in the simulated datasets. This percentage of noise is selected according to previous estimation in the microarray experiments [15]. For two-group data, we generated 100 no-shift genes (∆t = 0) and 100 shift genes (∆t = 2; Figure 3a), and for three-group data, we generated 100 no-shift genes (∆t = 0) and 100 shift genes with the randomly sampled time shift (∆t) (Appendix A; Figure 2d). In both settings, TempShift successfully distinguished the shift genes from no-shift genes by LLR shi f t , and accurately estimated ∆t (two-group data MSPE = 0.0095; three-group data MSPE = 0.0085 for ∆t 1 and 0.0097 for ∆t 2 ) (Figure 3b,c,e,f).

Polynomial Data
The non-periodic data representing gene expression changes across development were simulated with a residual variance accounting for approximately 10%. Then, TempShift was applied to the simulated non-periodic data with two-group and three-group, respectively. Again, TempShift successfully distinguished the shift genes from no-shift genes with LLR shi f t and accurately estimated ∆t using both two-group (MSPE = 0.012) and threegroup data sets (MSPE = 0.0052 for ∆t 1 and MSPE = 0.0065 for ∆t 2 ) (Figure 3h,i,k,l).
In summary, the simulation results demonstrated that TempShift is able to identify temporally shifted genes and accurately estimate time shift between groups even with high noise.

Identifying Developmental Sequences in the Human Neocortex
We then applied TempShift to investigate the developmental sequences in the human neocortex. The developmental sequences represent the order in which changes in structure or function occur during the process of development of an organism. The cerebral neocortex consists of functionally distinct sensory, motor and association areas. Previous studies have suggested differential gene expression and maturation rates of different neocortical areas.
To explore the developmental sequences in the human neocortex, we applied TempShift to a previously published microarray data set, including samples from 11 neocortical areas (Supplementary Table S1) across prenatal and postnatal development and adulthood [11]. Temporal dynamic genes were pre-selected by fitting a cubic regression model combining all regions (Details see Section 2). The 4153 selected genes were then used as input to TempShift. The dorsolateral prefrontal cortex (DFC) was used as the reference area (∆t = 0). Since the pre-selection step is prone to select genes with the same shape, the majority of the genes (4111 genes) passed the shape selection Λ shape = 50. Among these genes, 366 genes were selected as multi-area temporally shifted genes with threshold Λ shi f t = 50 (Figure 4a; Supplementary Table S2). On average, the selected temporally shifted genes showed delayed development in the prefrontal cortex (MFC, OFC, DFC and VFC) (Figure 4b), inferring the hierarchical development within neocortex [16,17]. The development of the primary visual cortex (V1C), the only area located in the occipital lobe in our data, is also delayed, but with higher variation. Gene enrichment analysis showed that the temporally shifted genes are enriched in synaptic transmission and neuron differentiation (Figure 4c).

Sequential Expression of Neurotransmitter Receptors and Axon Guidance Molecules
In the synaptic transmission group, we found multiple glutamate receptor genes (GRIA4, GRIK1 and GRIK4) (Figures 4d and A1a) and GABA receptor genes (GABRA1, GABRA4, GABRB2 and GABRG1) (Figures 4e and A1b). GRIA4 gene encodes GluR4 subunit of AMPA receptor, which is the main subunit expressed in the late postnatal and adult. This gene demonstrated delayed expression in the prefrontal cortex and the primary visual cortex (V1C). On the other hand, GRIK1 and GRIK4 encode Kainate receptor GluR5 and KA-1 subunits, respectively. The expression of GRIK1 started from the temporal lobe (ITC, A1C, and STC), while GRIK4 started its expression sporadically from ITC, VFC and MFC to other areas. The GABA receptor genes we identified encode subunits of GABAA receptors. Among them, GABARA1, GABRB2, and GABRG1 exhibited a developmental sequence that started from ITC and then spread through the temporal lobe to other neocortical areas. In contrast, GABRA4 started from the prefrontal cortex (OFC, VFC, and DFC) and was delayed in ITC. In the neuron differentiation group, we found four genes involved in the PANTHER pathway: axon guidance mediated by Slit/Robo [12]. One way that Slit/Robo signaling mediates repulsion from the midline is by silencing the receptor of the attractive guidance cue netrin-1, netrin-2 and DCC [18]. While all four genes exhibit delayed expression in the prefrontal cortex, ROBO1 displays opposite medial-lateral temporal patterns of NTNG1, NTNT2 and DCC. ROBO1 demonstrates the medial-to-lateral pattern, with earlier expression in V1C, S1C, IPC, M1C and A1C, while NTNG1, NTNG2 and DCC demonstrate the lateral-to-medial expression (Figures 4f and A1c).

Expression of Marker Genes Predicts the Interneuron Migration Pathway
In the above analysis, we also identified GAD1, a synthetic enzyme of interneuron neurotransmitter GABA [19], as a temporally shifted gene. The expression of GAD1, like GABA receptor genes GABRA1, GABRB2, and GABRG1, started from ITC and then spread from the temporal lobe to other neocortical areas. Unlike cortical projection neurons that derive from the dorsal telencephalon and migrate radically to the cortical plates [20], despite the existence of an additional subset of neocortex-originated neocortical GABAergic interneurons in primates, cortical interneurons mostly arise from the ventral telencephalon and migrate tangentially to the neocortex [21][22][23]. Therefore, we investigated the expression trajectories of markers of GABAergic interneurons (GAD1, GAD2) [19] and their progenitors (DLX1, DLX2) [24] to explore whether the migratory pathway of interneurons is reflected by a temporal shift in the gene expression. DLX1 and DLX2, transcription factors required for interneuron differentiation and migration [24], peaked in mid-fetal brains and were expressed until infancy (Figure 5a). On the other hand, GAD1 and GAD2, markers of pan GABA interneurons, reached highest expression in infancy following the peak of DLX1 and DLX2 expression in all neocortical areas (Figure 5a). Application of TempShift suggested that DLX1, DLX2, GAD1 and GAD2 exhibit the same temporal shift pattern of expression trajectories across neocortical areas. All four genes demonstrated sequential expression from ventrolateral to dorsomedial areas (Figure 5b,c), indicating the migratory streams of interneurons from ganglionic eminence to neocortex in humans (Figure 5d).

Identifying Developmental Changes in the Down Syndrome
To further explore the changes in neurodevelopmental processes and related genes in DS, we applied TempShift to a published data set including complementary DS samples [13]. The original study combined this data set with the human brain development set in the previous section and only used the age-matched controls for differential expression analysis. In this section, we analyzed samples from dorsal prefrontal cortex (DFC; Supplementary Table S5) and demonstrated that TempShift is able to include all control samples by comparing gene expression trajectories inferred from age-unmatched samples and to identify DiffShape and shift genes that indicate disrupted and time-shifted biological processes.

Increased Microglia Gene Expression in DS
We first selected DiffShape genes based on LLR shape and 260 genes were selected as DiffShape genes with LLR shape below two standard deviations below the mean ( Figure A2a and Supplementary Table S3). By integrating with single-cell RNA-seq [14], we found that the DiffShape genes are enriched in microglia (Fisher's exact test p < 0.05; Figure 6a). All 15 microglia-enriched genes show higher expression in DS than controls ( Figure A2c) and this difference increases with development ( Figure 6b).

Delayed Expression of Oligodendrocyte Genes
We identified 66 shift genes with criteria LLR shape above mean and LLR shi f t greater than 10 ( Figure A2b and Supplementary Table S4). Shift genes identified are enriched in oligodendrocytes (Figher's exact test p < 0.05; Figure 6a). All three genes enriched in oligodendrocytes (Figure 6d) demonstrated delayed expression in DS (Figure 6c,d,g). The detected shift genes include myelin basic protein (MBP; ∆t = 0.43), transmembrane protein 144 (TMEM144; ∆t = 0.7), and solute carrier family 5 member 11 (SLC5A11; ∆t = 2.5). The original study selected candidate genes based on paired t-test and gene co-expression network analysis. Only MBP were identified by gene co-expression network analysis, while no statistics were available for quantification. None of the genes were detected as differentially expressed based on paired t-test combining all samples.
Another cell type found to be enriched with shift genes is interneuron 2 (Figure 6a). This cluster of interneurons is enriched with genes involved in synapses. The shift genes enriched in this cell type include vesicle-associated membrane protein 1 (VAMP1; ∆t = 1.3) and regulator of calcineurin 2 (RCAN2; ∆t = 0.4), both of which exhibited obvious delayed expression in DS (Figure 6e,g). VAMP1, also known as synaptobrevin 1, is a member of the synaptobrevin family and is involved in the docking and/or fusion of synaptic vesicles with the presynaptic membrane. RCAN2, also known as DS Candidate Region 1-Like 1 (DSCR1L1), binds to the atalytic domain of calcineurin A and has been previously associated with DS. Another shift gene enriched in interneuron 2 is DMKN; however, since the estimated expression trajectory of this gene is linear, it is indistinguishable whether this gene is upregulated or shifted leftward in DS (Figure 6e).
Furthermore, additional shift genes playing critical roles in cortical development further implied delayed development in oligodendrocytes generation, neuronal migration, neurite growth and synaptogenesis (Figure 6c,g). CNTN6, also named NB-3, encodes the neural cell adhesion molecule contactin-6 and has been found to be delayed (∆t = 0.43) in expression in DS. Contactin 6 has been implicated as an autism risk gene. It interacts with NOTCH1 and promotes oligodendrocyte generation [25]. CNTN6 also interacts with cell adhesion molecule L1-like (CHL1) to regulate oriented growth of apical dendrites in the mouse neocortex [26]. Similarly, cerebellin 2 (CBLN2), a gene found to be involved in synaptogenesis induced by neurexin-neuroligin signaling [27], and synaptotagmin-2 (SYT2), a gene functioning as a Ca 2+ sensor for fast neurotransmitter release, are both delayed in DS, suggesting delayed synaptogenesis in the DS. In addition, ASTN1 encodes astrotactin, a neuron-glial adhesion molecule that mediates glial-guided neuronal migration [27]. The expression of ASTN1 decreases earlier in the DS (∆t = −4.2), indicating a shortened neuronal migration process in the DS, which may explain the reduced brain size and altered cortical lamination in DS [28].

Conclusions
TempShift is a framework that provides flexible modeling of global temporal shift of time-series data which can deal with an arbitrary number of replicates, does not require matched time points across conditions, and can be used for two or more conditions. For ease of use, we implemented TempShift into an open-source R package (https://github. com/YingZhuLab/TempShift, accessed on 31 May 2023) in this study. At first, we validated that TempShift works well for both periodic and non-periodic data and is robust to noise through simulation. We adopted the Gaussian kernel to fit the gene expression trajectories in this study, but it can be replaced with other kernels according to the properties of the data. As we maximize likelihood using a Quasi-Newton method, which finds local maxima, we would suggest using a larger Initial value of the length scale l of the Gaussian kernel to avoid overfitting. Otherwise, multiple Initial values can be tried to get the best fit. In summary, TempShift provides a framework that can be applied broadly to study temporal differences across different conditions, such as different tissue types, disease status and species. It can be applied not only to gene expression data, but also to other time-series measurements.
The implementation of TempShift to human brain transcriptome data demonstrated the capability of TempShift in identifying shift genes and estimating temporal shift between as many as 11 neocortical areas at the same time. In addition to comparing multiple groups, Tempshift is able to detect developmental sequences of multiple biological processes in a high-throughput way. In the above application, we found that shift genes are enriched in synaptic transmission and neuron differentiation, and reconstructed the migratory streams of interneurons in the human neocortex.
Using the DS data, we demonstrated the application of TempShift to analyze groups of time-series data with unmatched ages. Using TempShift, we selected DiffShape genes and shift genes, each of which respectively suggested increased microglia, and altered regulation of developmental timing in the DS, including delayed development of oligodendrocytes, neurite outgrowth, synaptogenesis and shortened period of neuronal migration. Anatomical changes observed in DS include reduced brain size, altered cortical lamination, reduced dendritic ramifications, diminished synaptic formation, and delayed myelination [28]. TempShift successfully detected changes in these processes and provided a list of candidate genes associated with these changes in the developmental processes for future functional studies.
The TempShift detects shape difference and global shift of gene expression trajectories currently. For further studies, other models can be developed to identify more transformation of trajectories, or to refine the time interval, during which the trajectories are different in shape or temporally shifted.
In summary, we believe that not only can TempShift be used for transcriptome data analysis of the human brain, but that it also has great potential for understanding the temporal dynamics of gene expression in other biological processes.
Supplementary Materials: The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/brainsci13071070/s1, Table S1: Human NCX samples.csv; Table S2: Shift genes in the human neocortex.xlsx; Table S3: DiffShape genes between DS and control.xlsx; Table S4: Shift genes in between DS and control.xlsx; Table S5: Down Syndrome samples.csv.  Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The normalized human brain microarray data used in this study are openly available in GSE25219 [11]. The microarray data of DS samples used in this study are openly available in GSE59630 [13]. The single cell RNA-seq of human neocortex used in this study are openly available in GSE67835 [14].

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A. The Mathematical Principle of Tempshift
In our previous study on spatiotemporal transcriptomic divergence across human and macaque brain development [6], we proposed TempShift to distinguish DiffShape genes, shift genes and no-shift genes across groups and to estimate the corresponding time difference for shift genes. This is accomplished through the following procedure. We consider one gene at a time in our analysis. Assume we have a total of m groups, and there are n i samples in the ith group. The total sample size of all m groups is n = ∑ m i=0 n i . In group i, the time vector of length n i is denoted as x i , and the corresponding gene expression vector is denoted as y i . The combined time vector across all the samples is where "T" denotes transpose. We consider the following three models: independence model (Equation (A1)), no-shift model (Equation (A2)) and shift model (Equation (A3)), where y follows a multivariate Gaussian distribution with different covariance matrices for these three models: Independence model: No -shift model: Shift model: where K(x i , x j ) is a covariance function or kernel between groups i and j, with n i rows and n j columns. In this study, we adopted Gaussian kernel K x i , where l denotes the length scale and σ 2 f denotes the amplitude, which determines the average distance of the function away from its mean. Other kernel functions can be incorporated in the framework if they provide better fit to the data.
For the independence model, the between-group covariance matrix is set to a zero matrix 0 (Equation (A1)). Each position of the covariance matrix represents the correlation between the horizontal and vertical elements, and if the covariance is a diagonal matrix, then each component of the multivariate Gaussian distribution is independent of each other. For the no-shift model, all groups are combined into one group (Equation (A2)). For the shift model, the first group is set as the reference group, and for group i (i = 1), the time shift is ∆t i and the shifted time vector is x i − ∆t i → 1 n i ; → 1 n i is a vector of all ones of length n i . The covariance matrix of the shift model is defined by combining the shifted time vectors into one group (Equation (A3)). The σ 2 is the variance of noise, and I n is the n × n identity matrix. Therefore, the no-shift model can be considered as a special case of the shift model when ∆t = 0 for all groups.
We evaluated whether a gene shows temporal shifts across groups in two steps: (1) we compared the shift model with the independence model to identify genes following the same trajectory across groups based on log-likelihood ratio LLR shape (Equation (A4)); we then compared the shift model with the no-shift model to identify genes with temporal shift by calculating log-likelihood ratio LLR shi f t (Equation (A5)).
LLR shape = log( max L ∆t, θ shi f t x, y max L θ independence x, y ), LLR shi f t = log( max L ∆t, θ shi f t x, y L ∆t = 0, θ shi f t x, y ), where ∆t is the vector of time shift from the reference group, θ shi f t and θ independence are the other parameters in the shift and independence models, respectively. L(·) is the likelihood function. The parameters were estimated using optim function in R.
The likelihood function L(·) measures the likelihood that the time-series data formed by the time vector x and the corresponding gene expression vector y follow independent model or shift model. Thus, we are able to filter Diffshape genes by comparing the values of LLR shape and and select shift genes by comparing the value of LLR shi f t .