Construction of a cross-species developmental brain meta-atlas.

a, Sampling strategy. Colored areas indicate the dissected structures in a PCD80 macaque brain; this regional color palette is used throughout. b, Composition of the atlas. Stacked bars show the number of samples ordered by post-conception day (PCD) per species, colored by region, see color legend in (a). c, Schematic of ANTIPODE. The encoder compresses gene counts to a latent variable Z; a classifier branch generates probability of each cluster k and ψ values; decoding is according to the equation shown with observed counts likelihood computed with the negative binomial distribution. d, Two-dimensional UMAP projection of the ANTIPODE latent space for all cells, colored by the 98 initial classes. e, UMAP projection colored by region (top, see color legend in (a)), species (middle) and by inferred cell-cycle phase (calculated by scanpy score_genes_cell_cycle) (bottom).

Taxonomy and ANTIPODE model architecture.

Heat-map summary of the 98 initial classes. For each class (rows) we display (left-to-right): regional abundance across broad areas, total cell count, Spearman correlation of gene expression across pairs of species, mean absolute inter-species log fold change (“ANTIPODE divergence”). Row dendrogram shows Kendall correlation hierarchical clustering based on ANTIPODE latent space.

Benchmarking ANTIPODE against existing cross-species integration methods.

a, Schematic detailing the benchmarking strategy for evaluating ANTIPODE as an integration method. b–e, Bubble plots summarizing integration metrics calculated by scib for: (b) mammalian retinal ganglion cells from Hahn et al.,, (c,d) cortical cross-areal (XA) and cross-species (XS) clusters from combined Jorstad et al. and Jorstad et al.. (f) whole developing brain initial classes from the current study f–h, Aggregated integration space k-nearest neighbor classification entropy for each method’s latent space calculated for (f) mammalian retinal ganglion cells from Hahn et al.,(g) combined cortical cross-areal (XA) and cross-species (XS) clusters from combined Jorstad et al. and Jorstad et al., (h) whole developing brain initial classes from the current study. Boxplots show median and inter-quartile range.

Modes and landscape of evolutionary differential expression.

a, Cumulative distribution showing proportion of genes by their per species log fold change. 71.6% exceed |log2FC > 1 in at least one pairwise species comparison (based on raw log expression, not ANTIPODE fit). b, Conceptual schematic of four categories of differential expression learned by ANTIPODE: DA (differential-by-all states intercept), DM (by-module), DC (by-co-expression) and DI (by-identity intercept). c, Boxplot showing median and quartile divergence (overall DE) for each ANTIPODE cluster, grouped by division. d, ANTIPODE UMAP colored by each cluster’s summary divergence (mean |LFC for all categories combined).

Gene expression evolution in context.

a, Gene size vs DA for all genes. Line of best fit calculated by scipy linregress. b, Best fit line slopes between divergence categories vs 0-1 scaled log10 gene size and 0-1 scaled log10 intergenic distance. c, Divergence of DA, DI, DM and DC sets split by gene class (Shared Neu TF refers to TFs with shared expression in development and adult in neurons, Shared NN TF refers to TFs with shared expression in development and adult in nonneurons, TFs refers to TFs not shared between development and adult, Non-TF refers to all other genes not in one of the other categories). P values are Bonferroni corrected two sided Mann-Whitney test, from python’s statannotations package * p<0.05, ** p< 0.01, *** p<0.001, **** p<0.0001. d, Divergence of DA, DI, DM and DC sets split by conserved synteny. P value annotations are the same as c. e, Summary divergence versus specificity (tau, 0-1 score where 1 represents gene is expressed only by one type) colored by log gene mean expression value. f, Schematic depicting evolutionary implications of paired neuropeptide and receptor gene expression divergence in sender and recipient initial classes, respectively. Neuropeptides are represented by keys and receptors are represented by locks. g–h, Histogram showing observed ligand–receptor correlations (blue) are more concordant in direction (g) and magnitude (h) than 10 000 shuffled pairs (orange); vertical lines indicate mean values.

A Bayesian model of progenitor-state progression across species and regions.

a, Progenitor timing model schematic. Progenitor abundances are fitted with asymmetric split-Gaussian curves with parameters: μ, σₐ, σᵦ, Δt, Δh; right inset is a cartoon depicting that larger structures tend to be later-developing. b, UMAP of progenitor classes used for timing model, colored by initial class. Abbreviations used throughout: (HB/RE: hindbrain/rhombencephalon, MB/ME: Midbrain/Mesencephalon, DE: Diencephalon, FB: Forebrain i.e. non-pallial, non-diencephalic prosencephalon, Ctx: Cortex i.e. pallium) c, Estimated centre-of-mass (COM) lateness for distinct progenitor state across regions (points colored by species) with NPC classes on the left and Astro-OPC progenitors on the right (these are the progenitor groupings that are analogous across all regions). Boxplots show median and quartile values. Boxplot including all values can be found in Supplementary Figure 9. d, Scatter of overall mean |logFC DE versus stereotypical COM for each progenitor state; color encodes the initial class. e, Schematic summarizing the known acceleration rostral-caudal spatiotemporal neurogenic gradient. f, Relative acceleration of cortical progenitor progression fit by progenitor state model using only cortical regions along the rostral-caudal axis for human and macaque.

Gene programs linked to early versus late neurogenesis.

a, Illustration of three lateness metrics: global, stereotypical and differential. b, Histogram of gene expression correlation scores with each progenitor lateness fit values c, Heatmap showing the correlation between the magnitude of each score category with the magnitude of each differential expression parcel’s overall divergence. d, Dot-plot of top genes correlated with each lateness metric (dot size = fraction expressing; color = log counts/count). e, Enrichment networks for genes associated with stereotypical lateness (left) and global lateness (right). Terms represent the top 15 positive and negative terms with resampling pvalue less than FDR<0.05. Edge width represents proportion of genes overlapping in gene sets, node color represents the mean lateness score for the genes in associated with that term. f, Gene set enrichment analysis (using gseapy prerank) of stereotypical lateness genes shows enrichment for neuromorphological disorder gene sets; curves plot directional enrichment score against ranked gene list.

Dataset composition, batch structure and basic QC.

a–c, ANTIPODE UMAP colored by post-conception day (PCD) for human (a), macaque (b) and mouse (c) reveal continuous developmental trajectories that align across species after integration. d, UMAP colored by the 15 10X scRNA-seq datasets showing that biological structure, not dataset of origin, dominates the embedding. e, The same UMAP colored by individual library-preparation batches. f, Stacked-bar plot of the 1,854,767 high-quality cells grouped by species, structure (colors) and PCD demonstrates broad temporal and anatomical coverage. g–i, Violin plots of log10 UMI counts per cell across datasets for human (g), macaque (h) and mouse (i).

Goodness-of-fit for ANTIPODE gene-specific parameters.

a, Reconstructed (posterior) mean log-expression (log counts/count) for every gene-cluster pair versus its empirical mean (black points; red density contours). b, Empirical zero probability (p(count = 0)) versus empirical log-mean expression suggest a small degree of zero-inflation (left-right shift of the curve). c, Full heat-map summary of the 98 initial classes from figure 2. For each class (rows) we display (left-to-right): regional abundance across broad areas, total cell count, Spearman correlation of gene expression across pairs of species, mean absolute inter-species log fold change (“ANTIPODE divergence”), membership in 200 ANTIPODE gene modules. The bottom section shows divergence by DC and DM. Gene set enrichments are colored by gseapy prerank normalized enrichment scores of each ANTIPODE component, with color alpha equal to 1 – significance q value. Row and column dendrogram shows Kendall correlation hierarchical clustering by ANTIPODE latent space.

Conserved marker expression across initial classes.

Dot plots show mean normalised expression (dot size) and mean inter-species log fold change (color scale; blue = human-up, orange = mouse-up, green = macaque-up) for representative markers of a brain-endogenous non-neuronal divisions, b brain-adjacent initial classes a.k.a neighbors c neuroepithelial / neural progenitor cells (NPC), d brain-exogenous non-neurons a.k.a Exo NN, e telencephalic (TE) neurons and f diencephalic/mesencephalic/rhombencephalic (DMR) neurons.

Mapping developing clusters to adult mouse subclasses.

a, UMAP of developmental clusters colored by their most-correlated adult subclass from Yao et al. 2023. b, Heat-map of Pearson correlations (using Yao identity TF set) between 380 developmental clusters (rows) and 337 adult subclasses (columns). c, Distribution of correlation coefficients for maximal matches (orange) versus all other comparisons (blue). d, Jaccard index of gene-marker overlap for the same developing cluster–adult subclass pairs (because matches are called by cluster all jaccard indices are 0 or 1).

Shared transcription-factor programs between development and adulthood.

a, Heat-map of binary transcription-factor (TF) expression (presence > 0.1 fraction of cells) across 337 adult cortical subclasses (rows); clustered blocks highlight subclass-specific TF sets. b, Equivalent heat-map to a for the 380 developmental clusters. c, Histogram of TF “bimodality” scores (1 = all-or-none expression) shows higher discreteness in adult (orange) than in development (blue) cell types (dashed lines: medians). d, Four example subclass/cluster pairs illustrating TF sharing: each scatter plots mean TF expression in the developmental state (x-axis) versus its matched adult subclass (y-axis); dashed lines mark the 0.33 threshold used to call shared TFs.

Qualitative comparison of integration methods on three benchmark datasets.

UMAPs colored by species (left of each pair) or cell identity (right) for methods used for analysis: ANTIPODE, scVI, Harmony, Scanorama and PCA. a, Adult cross-areal (XA) + cross-species (XS) primate cortical dataset integrations. b, Adult multi-species mammalian retina. c, Developmental whole-brain atlas used in this study.

Additional evolutionary gene expression divergence.

a, Spearman correlation across all species pairs for 380 developmental clusters (named by initial class + cluster). b, Mean correlation across species pairs shown in (a) by each cluster, plotted in UMAP space c, ECDF plots of raw divergence category parameters. d, Scatterplot of the correlation of DM, DC, and DI divergence categories pairwise. Note the consistently high correlation of DC and DM as these are multiplied together in the model (point color). e, Boxplot showing the number of genes with mean expression in clusters greater than 20 counts per million (cpm), grouped by division. f, Boxplot showing the effect of each divergence category for each cluster, grouped by division. Note that divisions that contain many clusters are biased towards divergence being categorized as DM+DC, while those consisting of fewer clusters seem somewhat biased towards DI.

Neuropeptide ligand–receptor expression landscape across progenitor and neuronal states.

a, Tile-map of neuropeptide-receptor interaction scores (multiplied 0-1 scaled expression). Ligands in neurons (rows) vs receptors in non-neurons (columns). Colored by ligand in main (central heatmap). Along left and bottom shows raw 0-1 expression of ligands/receptors in each species, and tricolor species divergence tile-map. b, Tile-map showing species divergence, with 0-1 scaled tri-species means colored according triangular legend.

Additional analyses of the progenitor timing model.

a, Model-free trapezoid sum-based COM calculations based only on raw data for species and regions (colors). b, Heat-map of annotated finer regional dissection vs de novo inferred regions. c, Heat-map of annotated broad regional dissection vs de novo inferred regions. d, Species-level time-warp parameters (scale α, shift Δ) estimated from the model; mouse shows a pronounced left-shift (earlier development). e, Global comparison of σₐ versus σᵦ for all progenitor states. f, g, Heat-map of ANTIPODE cluster (f) or batch_name (g) vs de novo inferred regions h, i, Box plot of fit COM values for each progenitor region + cluster, colored by species in arbitrary warped time (i) or real time (j) and ordered by mean time across species j, Stacked barplots showing proportion of total cells included of each progenitor state at each timepoint, for each species, grouped by simplified regions.