A single-cell multi-omic atlas spanning the adult rhesus macaque brain

Abstract

Cataloging the diverse cellular architecture of the primate brain is crucial for understanding cognition, behavior, and disease in humans. Here, we generated a brain-wide single-cell multimodal molecular atlas of the rhesus macaque brain. Together, we profiled 2.58 M transcriptomes and 1.59 M epigenomes from single nuclei sampled from 30 regions across the adult brain. Cell composition differed extensively across the brain, revealing cellular signatures of region-specific functions. We also identified 1.19 M candidate regulatory elements, many previously unidentified, allowing us to explore the landscape of cis-regulatory grammar and neurological disease risk in a cell type-specific manner. Altogether, this multi-omic atlas provides an open resource for investigating the evolution of the human brain and identifying novel targets for disease interventions.

Fig. 1.. Experimental setup and summary of the Macaque Brain Atlas snRNA-seq dataset.

(A) Schematic of biopsied brain regions for sci-RNA-seq3 experiment. A full list of sampled regions is provided in table S1. arc, arcuate sulcus; cgs, cingulate sulcus; cs, central sulcus; ecal, external calcarine sulcus; iarc, inferior arcuate sulcus; ic, internal capsule; ips, intraparietal sulcus; ls, lateral sulcus; lv, lateral ventricle; p, principal sulcus; rf, rhinal fissure; sarc, superior arcuate sulcus. (B) UMAP visualization of all snRNA-seq profiled cells colored by cell type [with color code shown in (C)]. (C) Barplots showing the log₂ -transformed cell counts (left), regional specificity score (middle), and regional composition [right, with color code shown in (E)] of each cell type. (D) UMAP visualization of all snRNA-seq cells colored by cell type [with color code shown in (E)]. (E) Barplots showing the cell type composition [left, with color code shown in (C)], log₂-transformed ratio of glutamatergic neurons and GABAergic neurons (middle), and log₂-transformed ratio of neurons and glial cells (right) of each region. Regions are organized by the regional subclass to which they belong.

Fig. 2.. Cell subtype distribution and variation across the brain.

(A) Barplots showing the region specificity score (i.e., Jensen-Shannon divergence statistic) and composition for cell subtypes (with color code shown in Fig. 1E). (B) Heatmap showing scaled log₂ ratios of GABAergic neuron and astrocyte subtype compositions within cortical region, compared to the average across all regions. Cell subtypes with at least 100 cells profiled are shown in the order of abundance (x axis, left to right) in the cortical regions organized by region subclasses (y axis). The color and direction of each pie correspond to relative enrichment (blue, clockwise) and depletion (red, anticlockwise) of a cell subtype in a region. Log₂ ratios were capped at positive and negative 2 before scaling. (C) UMAP visualizations of GABAergic neurons colored by cell subtype (left) and regional subclass (right). (D) UMAP visualizations of GABAergic neurons colored by cell subtype marker gene expression. (E) UMAP visualization of astrocytes colored by the region with the highest lochNESS, indicating enrichment of a region subclass in the cell’s transcriptional vicinity. LochNESS distribution in a few example regions (occipital lobe, basal ganglia, brainstem, thalamus, and frontal lobe) is highlighted in separate panels as examples. (F) UMAP visualizations of astrocytes colored by lochNESS-derived region-related marker genes.

Fig. 3.. Generation of the Macaque Brain Atlas sci-ATAC-seq dataset and identification of cell classes.

(A) UMAP visualization of all snATAC-seq cells colored by brain region [with color code shown in (C)]. (B) UMAP visualizations of promoter accessibility scores of cell markers (GAD2: GABAergic neurons, ENTPD1: microglia, SLC1A2: astrocytes, ATP10A: vascular cells) reveal high specificity. (C) Barplots showing nuclei counts by brain region of the snRNA-seq, snATAC-seq, and integrated datasets. (D) UMAP visualizations of integrated multimodal data, with cell classes colored separately for (left) snRNA-seq and (right) snATAC-seq nuclei [with color code shown in (F)]. (E) Spearman’s rank correlation coefficients showing the correlation between cell class proportions in the snRNA-seq and snATAC-seq datasets within each region (representing data generated from the same homogenized sample). (F) Scatterplot showing the correlation between cell class proportions in the overall snRNA-seq and snATAC-seq datasets (combined across brain regions). (G) Integration-derived cell class annotations visualized over the same snATAC-seq UMAP visualization shown in (A) [with color code shown in (F)].

Fig. 4.. Enrichment of TF binding site motifs in candidate regulatory elements.

(A) Barplots showing summary statistics for peak sets called separately on reads derived from cells assigned to each of 11 cell classes. (B) Heatmap showing enrichment (log₂ OR) of TF binding motifs among cell classes. The top five most-enriched nonredundant TF motifs (all P_adj < 0.05) are shown per cell class, ordered from left to right by increasing P_adj. Log₂ OR color ranges are capped at ±1.5. (C) Position weight matrices of the most-enriched TF motifs for six example cell classes. ORs are shown in parentheses. (D) Scatterplots showing correlation between snATAC-seq accessibility of TF binding motifs and snRNA-seq gene expression of corresponding TF genes within cell classes in regional classes for four example TFs.

Fig. 5.. The landscape of cis-regulatory interactions in the Macaque Brain Atlas.

(A) Schematic outlining criteria for identification of cCREs. Squares represent peak:gene pairs, darker colors symbolize stronger evidence for a given measure, and solid borders represent statistically significant measures. (B) Distribution of gene-peak GLUE regulatory scores binned according to the minimum signed distance (left: upstream, right: downstream) between peaks and gene TSSs. Distributions are shown separately according to whether the gene:peak pair also exhibited a significant association (P_adj < 0.05) based on the metacell-based logistic regression analysis. The ratio between significant and not-significant gene-peak pairs for proximity bins according to the logistic regression model is shown in the top margin, while the distribution of GLUE regulatory scores is shown in the right margin. (C) Candidate regulatory elements are shown in relation to, from top to bottom, (i) the strength of and direction of inferred regulatory links connecting peaks to MBP expression (based on metacell logistic regression analysis). The height of links represents the strength (−log₁₀ P values) of evidence for regulatory connections and the color symbolizes the direction of the relationship; (ii) the differential accessibility (−log₂ FC) of peaks in oligodendrocytes relative to all other cell classes; (iii) the distribution of normalized snATAC-seq reads by cell class; (iv) gene and transcript boundaries of MBP and its known isoforms in the rhesus macaque genome, with exons shown in blue; (v) the distribution of normalized snRNA-seq reads by cell class. Oligodendrocyte reads are shown in relation to all other cell classes on the upper portion of the plot. On the bottom portion, the y axis is magnified ×60 and cropped to highlight more subtle differences among cell classes.

Fig. 6.. Enrichment of heritable disease-relevant sites among candidate regulatory elements.

The heatmap displays heritability enrichment (log₂ OR) of diseases among cell class snATAC-seq peaks for tested diseases, syndromes, and phenotypes. Only results passing a threshold of P_adj < 0.05 are shown. The log₂ OR color range is capped at 3.0.