Tempo and mode of gene expression evolution in the brain across primates

  1. Katherine Rickelton  Is a corresponding author
  2. Trisha M Zintel
  3. Jason Pizzollo
  4. Emily Miller
  5. John J Ely
  6. Mary Ann Raghanti
  7. William D Hopkins
  8. Patrick R Hof
  9. Chet C Sherwood
  10. Amy L Bauernfeind
  11. Courtney C Babbitt  Is a corresponding author
  1. Department of Biology, University of Massachusetts Amherst, United States
  2. Molecular and Cellular Biology Graduate Program, University of Massachusetts Amherst, United States
  3. Department of Anthropology and Center for the Advanced Study of Human Paleobiology, The George Washington University, United States
  4. MAEBIOS Epidemiology Unit, United States
  5. Department of Anthropology, School of Biomedical Sciences, and Brain Health Research Institute, Kent State University, United States
  6. Department of Comparative Medicine, Michale E. Keeling Center for Comparative Medicine,The University of Texas M D Anderson Cancer Centre, United States
  7. New York Consortium in Evolutionary Primatology, United States
  8. Nash Family Department of Neuroscience and Friedman Brain Institute, Icahn School of Medicine at Mount Sinai, United States
  9. Department of Neuroscience, Washington University School of Medicine, United States
  10. Department of Anthropology, Washington University in St. Louis, United States

Abstract

Primate evolution has led to a remarkable diversity of behavioral specializations and pronounced brain size variation among species (Barton, 2012; DeCasien and Higham, 2019; Powell et al., 2017). Gene expression provides a promising opportunity for studying the molecular basis of brain evolution, but it has been explored in very few primate species to date (e.g. Khaitovich et al., 2005; Khrameeva et al., 2020; Ma et al., 2022; Somel et al., 2009). To understand the landscape of gene expression evolution across the primate lineage, we generated and analyzed RNA-seq data from four brain regions in an unprecedented eighteen species. Here, we show a remarkable level of variation in gene expression among hominid species, including humans and chimpanzees, despite their relatively recent divergence time from other primates. We found that individual genes display a wide range of expression dynamics across evolutionary time reflective of the diverse selection pressures acting on genes within primate brain tissue. Using our samples that represent a 190-fold difference in primate brain size, we identified genes with variation in expression most correlated with brain size. Our study extensively broadens the phylogenetic context of what is known about the molecular evolution of the brain across primates and identifies novel candidate genes for the study of genetic regulation of brain evolution.

Editor's evaluation

This is an important study that represents a significant contribution to our understanding of how gene expression in the primate brain has evolved across the extant primate phylogeny. It provides solid evidence for potential links between gene expression variation and brain size, although these are somewhat limited by the focus only on adult brains, since many key changes likely occur during development. Nevertheless, both the taxonomically broad data set and the analysis are likely to be of broad interest to the evolutionary biology, anthropology, and comparative neuroscience communities.

https://doi.org/10.7554/eLife.70276.sa0

Introduction

Primates are distinguished from other mammals by their large brains relative to body size (Boddy et al., 2012; Martin, 1981; Smaers et al., 2021). Among the diversity of primate species, there is remarkable variation in behavioral specializations, including differences in social structure, spatial, dietary and visual ecology, and locomotion (Powell et al., 2017; Barton, 2012; DeCasien and Higham, 2019). Despite the impressive array of cognitive attributes displayed by primates, molecular and cellular studies investigating various aspects of brain evolution tend to sample from a small number of species to address questions of how humans are unique. In large part, the emphasis on human brain evolution is warranted. Humans are unmatched in possessing exceptionally large brains and unparalleled cognitive abilities, such as language (Konopka and Roberts, 2016; Rilling, 2014). While valuable, the limited number of species included in prior research lacks a comprehensive perspective of the phylogenetic context in which the human brain evolved within the diversity of primates.

Although researchers have used a variety of approaches to assess whether the human brain is unique (Stout and Hecht, 2017) and how it might have evolved (Hrvoj-Mihic et al., 2013; Sousa et al., 2017a), the full potential for using gene expression to evaluate patterns of brain evolution in primates has not yet been met. Upon observing the remarkable similarity between human and chimpanzee protein sequences, King and Wilson, 1975 proposed that the basis of the physical and behavioral phenotypic differences between these two species must be found in changes within gene regulatory regions that drive expression. Previous studies have explored how changes in regulatory regions can influence gene expression but have often sampled various organs from species across broad spans of evolutionary time, such as mammals or vertebrates (Brawand et al., 2011; Breschi et al., 2016). In studies focusing on gene expression in primate brain tissues, research has mostly focused on the neocortex and cerebellum in human, chimpanzee, and rhesus macaque (Babbitt et al., 2010; Blekhman et al., 2010; Khaitovich et al., 2005; Khaitovich et al., 2004; Khrameeva et al., 2020; Konopka et al., 2012; Ma et al., 2022; Somel et al., 2009; Sousa et al., 2017b). However, new insights can be gained by sampling at greater neuroanatomical resolution from a broader array of primates. Examining gene expression of the brain from a more comprehensive landscape empowers novel inquiry in primate brain evolution, including questions pertaining to the sources of variation that drive expression differences, rates of expression change across the primate phylogenetic tree, and genes that correlate with brain size across primates.

In the current study, we sampled prefrontal cortex (PFC), primary visual cortex (V1), hippocampus (HIP), and lateral cerebellum (CBL) from 18 primate species, the broadest diversity of primates sampled in any study of gene expression in the brain to date, including species from several rarely-studied lineages. Our dataset represents 70–90 million years (Perelman et al., 2011) of primate evolution, providing a more thorough understanding of the evolution of gene expression across primates and allowing for an unprecedented view of how gene expression in the brain has changed over time across all major clades of primate phylogeny.

Results

Most variation can be explained at the species level, not by brain region

To understand how gene expression in the brain has evolved across the primate lineage, we generated and analyzed RNA-seq data from 18 primate species (including five hominoids, four cercopithecoids, four platyrrhines, and five strepsirrhines, with 1–3 biological replicates) across four brain regions, including PFC, V1, HIP, and CBL (Figure 1, Supplementary file 1). The transcriptomes and gene models were assembled de novo (Haas et al., 2013) (see Materials and methods). We quantified the expression of 15,017 orthologs within hominoids, and 3432 on-to-one orthologs across all 18 species (Supplementary files 2 and 3). Variation in interspecific mammalian gene expression has been shown to be less pronounced than that observed across samples from different organs, reflecting the diversity of underlying organ physiology (Brawand et al., 2011; Sudmant et al., 2015). Furthermore, it has been reported that the rate of gene expression divergence evolved more slowly within the cerebral cortex and cerebellum compared to other organ systems from developmentally distinct germ layers (Brawand et al., 2011; Khaitovich et al., 2005). To explore the variability of gene expression from distinct regions of the brain across our broad sampling of primates, we constructed a pairwise distance matrix of the 500 most variable protein-coding genes based on the standard deviation of expression across samples (Methods). This subset of genes was enriched with glycoproteins, signal peptides, and plasma membrane proteins, with roles in immune function, molecular trafficking, and cell signaling. Using this distance matrix, we performed a principal coordinates analysis (PCoA) on data from all brain regions. Because our samples represent disparate regions of the same organ, we expected less variation to be attributed to brain regions than primate species or taxa, reflecting the similarity in physiology of brain tissues. Unsurprisingly, the variation from our complex gene expression dataset is represented across multiple axes of the PCoA (Supplementary file 4).

Primate phylogeny showing the eighteen species sampled in this study.

The scale bar for the branch lengths represents 10 million years of evolution. The phylogenetic tree is a consensus tree of 1000 iterations produced from 10kTrees v.3 (https://10ktrees.nunn-lab.org) based on data from GenBank. The insets demonstrate the approximate locations of the four brain regions sampled on a coronal section, midsagittal view, and lateral view (displayed left to right, respectively) of a schematized adult human brain.

We plotted the first three axes and created polygons around data derived from samples sharing a common taxa (Figure 2a–c) or region (Figure 2d–f). As predicted, taxon assignment explains a large amount of variation to the dataset, with clear trends emerging independent of brain region. We find the greatest divergence in expression patterns among hominoid and strepsirrhine species, while there is more similarity observed among cercopithecoids and platyrrhines. The hominoids, displaying the greatest level of diversity of any primate phylogenetic group, demonstrate variation that is particularly apparent along Axis 1 and largely driven by human and chimpanzee expression patterns (Figure 2A, B). Strepsirrhines also exhibit a large amount of variation, especially apparent along Axis 2, which can mostly be attributed to the three species of lorises. When Axis 1 and 2 of the PCoA are plotted on the same bivariate plot (Figure 2a), the hominoids display more variation than the strepsirrhines by about 24% (Supplementary file 5). However, a large portion of the variation in the strepsirrhines is attributed to evolutionary divergence over about 63 million years (since the last common ancestor of lemurs and lorises), whereas the variation within hominoids has largely accrued over only 9 million years (since humans and chimpanzees shared an ancestor with gorillas). The hominoid and strepsirrhine samples represent similar variation in terms of sex and life stage, suggesting that these factors do not account for the variability seen in these taxa. Therefore, a remarkable finding of this analysis is how much variation is represented by hominoids, despite the fact that this lineage represents a much shorter evolutionary divergence time.

Figure 2 with 3 supplements see all
Patterns of brain gene expression across primates.

The first three axes of a principal coordinates analysis (PCoA) are plotted in both rows but have different symbols and colors to emphasize expression patterns specific to taxa (upper row, a–c) and regions (lower row, d-f). Polygons in each plot surround the data points for taxa (upper row) and regions (lower row). Axes 1, 2, and 3 represent 12.8, 10.3, and 9.4% of variance, respectively.

The cerebellum differs most significantly from the other three sampled brain regions

We observed trends in gene expression by brain region that are predominantly seen along Axis 3 (Figure 2e–f). Here, we observe that the variation attributed to CBL is distinguished from that of PFC, V1, and HIP, which are very similar in their distributions. The fact that CBL differs in its pattern of gene expression from other brain regions (Hawrylycz et al., 2015; Hawrylycz et al., 2012; Itõ, 2012) is not surprising given that it is the only sampled region that develops from a different part of the embryonic neural tube (namely, the hindbrain vs. forebrain) and exhibits a neuronal packing density (predominantly glutamatergic granule cells) that far exceeds these other brain regions (Azevedo et al., 2009; Herculano-Houzel, 2011; Itō, 2012). Enrichments for cerebellar gene expression reveal expression changes to categories such as ‘Cell Surface Receptor Signaling Pathway,’ ‘Cell Projection Organization,’ and ‘Wnt Signaling Pathway.’ This is true for humans in relation to chimpanzees as well as other species with deeper evolutionary relationships. (Supplementary files 6 and 7). Notably, genes involved in cell migration and cell surface receptor signaling Emera et al., 2016 have been shown to mediate the cell-cell interactions necessary for axon guidance (Koropouli and Kolodkin, 2014). Although all the brain regions surveyed show some enrichment for categories related to cell signaling, the cerebellum shows a unique increase in signaling activity-related terms (both in the number of enrichment categories and the degree of change associated with these categories). This is indicative of region-specific differential expression.

Gene expression profiles accurately replicate phylogenetic relationships

To determine whether gene expression profiles can reconstruct known phylogenetic relationships among primates, we built expression phenograms (Methods) for all brain regions combined (Figure 2—figure supplement 2) and each region separately (Figure 2—figure supplement 3). Samples that were derived from individuals of the same species tended to be grouped together, regardless of brain region, revealing that inter-individual differences are minor compared to other sources of variation. Gene expression profiles also replicated the phylogenetic relationships of closely related species (e.g. humans and chimpanzees; pig-tailed and rhesus macaques) when all regions were considered, but these relationships became less phylogenetically structured in the phenograms constructed using expression data from individual brain regions. All neighbor-joining phenograms accurately represent cercopithecoids and strepsirrhines as monophyletic groups; however, expression data produces paraphyletic groups of hominoid and platyrrhine species. This result potentially reflects the fact that taxa with longer periods of independent evolution (i.e. strepsirrhines) are more likely to show divergent patterns of gene expression than more closely related groups. Meanwhile, a more dense sampling of cercopithecoids (three individuals per species), permits a fairly accurate reconstruction of this taxon.

Differential expression in the context of the Ornstein-Uhlenbeck model

Previous studies have used a variety of different approaches to model gene expression changes over time (Brawand et al., 2011; Perry et al., 2012). Here, we used a recently described Ornstein–Uhlenbeck (OU) model to analyze neutral and conserved processes as determined by changing gene expression levels (Chen et al., 2019). OU processes have been proposed to model gene expression evolution as they model both drift and stabilizing selection (Rohlfs et al., 2014). Previous studies have shown how models that incorporate stabilizing selection are more accurately able to predict gene expression evolution in mammals than models that account only for neutral drift (Chen et al., 2019). Across all expressed one-to-one orthologs represented in the sampled primates, we found that ~15–20% of genes show differential expression across all species-to-species comparisons (q-value <0.05). As expected, the relative amount of differential expression increases over evolutionary time, both between species and clades (Figure 3 and Figure 3—figure supplement 1 ). However, both the species and clade-wise comparisons show larger numbers of differentially expressed (DE) genes in the comparisons when strepsirrhines are included in the contrast, a result of the more than 60 million years of independent evolution of this taxon. In the prefrontal cortex, we see slightly smaller DE in humans as compared to siamang and baboon; however, in total number, this DE is not appreciably different (DE Human-Siamang n=251, DE Human-Baboon n=255) to other comparisons (DE Human-Rhesus Macaque n=315). Upset plots allow us to further analyze these unique patterns of differential expression (Figure 3—figure supplement 2). In the cerebellum, the lemur, baboon, and human samples are particularly unique and show relatively higher levels of DE compared to other species. However, in the hippocampus, visual cortex, and prefrontal cortex, the chimpanzees appear to show higher DE than the human samples. This suggests again that the cerebellum is a particularly unique brain structure and that chimpanzee gene expression is significantly different in the other three brain regions studied, warranting further analysis.

Figure 3 with 2 supplements see all
Gene counts of differentially expressed (DE) genes between species and clades.

Each row represents one of the four brain regions examined. The size of the circle represents the number of DE genes seen at q<0.05 (5% FDR). The comparisons on the left are between exemplar species or sets of species, comparisons on the right are between clades of primates.

Human-specific enrichment for metabolic processes, neural development, and gene regulation

When comparing gene expression in the PFC of humans relative to other primates, human PFC shows an enrichment of metabolic processes, including ‘regulation of cellular metabolic process’ and ‘regulation of macromolecule metabolic process’ (Supplementary files 6 and 7). Comparing human and chimpanzee PFC reveals that categories that support neural growth and development (e.g. ‘neuron projection morphogenesis,’ ‘cell morphogenesis involved in neuron differentiation’), gene regulation, and metabolic processes are enriched as differentially expressed in human PFC relative to chimpanzees.

In addition to examining human and chimpanzee data in isolation, we also analyzed human-specific changes in expression within the context of other outgroup species. Using the Ornstein-Uhlenbeck process as a model of continuous trait evolution across our 18-species primate phylogeny, we again observe similar categories of enriched processes for the human prefrontal cortex in comparison to that of chimpanzees. These terms include ‘Nervous System Development,’ ‘Neurogenesis,’ ‘Glial Cell Differentiation,’ ‘Neuron Projection Morphogenesis,’ ‘Regulation of Gene Expression,’ and ‘Metabolism.’ To further validate our results, we also looked at the human PFC in comparison to other primate species. In analyzing differential expression between the human and siamang PFC, we note that similar trends for enrichment are also found, such as ‘Neural Growth and Development’ and ‘Metabolic Processes,’ and ‘Gene Regulation.’ Of interest, under the category of ‘Positive Regulation of Transcription by RNA polymerase II’ we find several genes that appear to be upregulated in humans compared to siamang: APP (amyloid precursor protein, related to plaque formation in Alzheimer’s disease) as well as PRKN (found to be causal in Parkinson’s disease) (Funayama et al., 2023). This supports the idea that these enrichments are human-specific, have relevance to important human neurodegenerative disease states, and are not a reflection of changes occurring within the chimpanzee lineage.

Broader species comparisons show similar trends across evolutionary time

Beyond human and chimpanzee comparisons, we also note many interesting broader temporal trends observed from the EVEE-based differential expression analysis. When examining Hominini (humans and chimpanzees) compared to other Great Apes, terms related to ‘Regulation of Metabolic Process,’ ‘Nervous System Development,’ and ‘Biosynthetic processes’ are all enriched within the hippocampus. Meanwhile, ‘Negative Regulation of Synaptic Transmission’ was enriched in the other ape species. In comparing the PFC of the Hominoid clade to that of the Strepsirrhine clade, we found that ‘Neuron Development and Differentiation’ was enriched. Overall, among various species and clade comparisons, there is a general trend of decreasing specificity in enrichment categories over increasing evolutionary time. Using PFC data, we found that the relationship between humans and chimpanzees, some of the closest relatives in our dataset, shows terms related to ‘Synapse Assembly,’ ‘Regulation of Glial Cell Differentiation,’ ‘Regulation of Astrocyte Differentiation,’ ‘Axonogenesis,’ and ‘Neuron Projection Morphogenesis.’ Looking at a more distantly related species pair, the human and rhesus macaque comparison shows enrichment for terms related to ‘Cell Growth,’ ‘Cell Development,’ ‘Biological Regulation,’ ‘Neuron Projection Development,’ ‘Regulation of Neurogenesis,’ and ‘Positive Regulation of RNA Biosynthetic Processes.’ The most distantly related species, humans and lemurs, have overall the largest number of differentially expressed genes, and with that, the broadest categories of enrichment. These include categories such as ‘Regulation of Developmental Processes,’ ‘Regulation of Nervous System Development,’ ‘Cell Development,’ and ‘Multicellular Organism Development’.

Our PCoA analyses showed that gene expression in brain regions sampled from the lorises (i.e. the slender loris, slow loris, and pygmy slow loris) diverged from other strepsirrhines, and other primates more generally (Figure 2—figure supplement 1). When strepsirrhines are compared to other primates in differential expression analyses, transcription factors, and other genes involved in gene transcription and translation and multiple biosynthetic pathways involved in cellular metabolism are among the categories of DE genes (Supplementary files 6 and 7).

Evolutionary rates of expression change across clades and brain regions

Using the OU model, we found that individual genes exhibit wide variation in expression dynamics across the primate lineage (Figure 4a). Enrichments for genes showing low variation or stabilizing selection (q=0.05) reveal categories related to transport and cellular localization (GO Biological Processes, Supplementary files 6 and 7). In contrast, genes that are less constrained or neutrally evolving (q>0.05) have a number of processes related to neuron morphogenesis, plasticity, and cell death. Yet, unlike sequence evolution, gene expression is not linear across evolutionary time but a saturation point in pairwise comparisons of gene expression is reached due to stabilizing selection pressures. Here, we find that pairwise expression differences between humans and the other species increasingly diverge with evolutionary distance in all brain regions sampled (Figure 4b); however, these pairwise comparisons do not seem to saturate with evolutionary time across the primate comparisons (Chen et al., 2019). We note that the saturation of pairwise expression differences from humans may be found at a phylogenetic node ancestral to primates (Figure 3—figure supplement 1).

Figure 4 with 1 supplement see all
Rates of change over genes and evolutionary time.

a. Exemplar genes that show constraint (left panel) and variation (right panel) across primates (colors as in Figure 2). b. Mean squared expression difference plotted by evolutionary distance to humans across all orthologs that were expressed. Shapes denote the four brain regions, and the colors represent the four major primate clades represented in our samples.

Expression of a majority of genes evolves under stabilizing selection

We utilized EVEE-tools, developed in the context of the Ornstein-Uhlenbeck process of continuous trait evolution, to classify genes as primarily under the effects of stabilizing selection vs. a model of neutral drift (Chen et al., 2019). In this analysis, we found that across tissues, on average, 64% of genes fit better under stabilizing selection (64% CBL, 59% HIP, 72% PFC, and 60% V1; FDR-corrected q-value calculated via the BH procedure to correct for multiple hypothesis testing; FDR threshold of 5% to determine significance). In the context of the OU model of continuous trait evolution, we found that the overall phylogenetic signal in brain expression divergence was slightly smaller than observed with edgeR over the entire combined dataset (EVEE: 4–6.7% across all single species comparisons using a logFC of +/-2; edgeR: 15–20%). However, the relative amount of differential expression does increase gradually with evolutionary distance (as expected based on results in edgeR). The only exception to this is the human and chimpanzee comparison, which shows considerable variation in differential expression across the four brain regions. Representing such a short period of evolutionary divergence, this increase in DE suggests direction selection within that lineage.

Evidence of potential directional selection in human and chimpanzee data

To better understand the outlier gene expression in this dataset, and overall to gain insight into genes that may be subject to directional selection pressures, we again used the EVEE-tools OU-based model to score our dataset for outlier expression. This requires the determination of the evolutionary mean and variance for each gene across our entire expression dataset, from which we then can compare to individual species expression data. Using this method, we found a small subset of genes to be defined as having an expression that deviates from the optimal OU-distribution (Z-score >2 or < –2, p-value <0.05) (Figure 4Figure 4—figure supplement 1). It is important to note that a significant FDR (<5%) was not reached in this dataset, among all species comparisons. This is expected based on previous applications of this package, in which a mammalian dataset was also unable to reach an FDR below 18% (Chen et al., 2019). This is likely a reflection of the limitations of our phylogeny, and suggests that future projects should aim to sample even more broadly.

In our outlier expression analysis, we found that patterns of outlier gene expression occur in a species and tissue-specific manner (Supplementary file 10). For example, in the PFC, the chimpanzee and marmoset samples appear to have the highest number of outlier genes deviating from average expression patterns (compared to human, siamang, baboon, rhesus macaque, and marmoset samples). We found this to be true regardless of how the dataset was normalized in order to determine the average expression for each gene (via defining a reference species). This is particularly interesting and, in combination with differential expression analysis, highlights the chimpanzee PFC as a particularly divergent structure. In contrast to the PFC, outlier analysis in the CBL reveals that the human and lemur samples have the highest number of outlier genes.

Upon gene set category enrichment analysis, we find that many of these genes that are deemed ‘outliers’ are related to functions in development, transcription, nervous system development, neurogenesis, and metabolism. For example, the chimpanzee PFC shows a significant upregulation in genes involved in energy storage and transfer, such as ETFA and NDUFS4, which are both involved in electron transport for ATP generation, as well as ANKH, implicated in phosphate transport (Szeri et al., 2020; Henriques et al., 2021; Shil et al., 2021) We also see that the chimpanzee PFC shows unique expression patterns in genes related to synaptic activity and neurotransmitter release, including significant downregulation of GABRA4 and SYN1 (Fassio et al., 2011; Fan et al., 2020 Fan et al., 2020; Fassio et al., 2011). In contrast, the human CBL and PFC both display a significant upregulation of genes related to Amyloid protein production (APP), a major component of many neurodegenerative diseases with functions in synaptic signaling (O’Brien and Wong, 2011). Unique to the human CBL we also see enrichments in genes related to neurogenesis and synaptic activity, including SDK1, FZD5, and CDH10 (Redies et al., 2012; Slater et al., 2013; Bagot et al., 2016) This is suggestive of directional selection pressures occurring in a tissue-specific manner and encourages future investigation of these outlier genes.

Implications of using humans as a reference species

We continue to use humans as a reference species in these analyses as, compared with other primates, humans have exceptionally large brain sizes and unique cognitive abilities. However, we do recognize that there are some implications for having humans as a reference, especially given our data that would suggest human gene expression as being largely different from the rest of our primate dataset. To address this, we repeated analyses using EVEE-tools by including two additional reference species: siamang and rhesus macaque. The percentage of genes that fit better under the model of stabilizing selection (in comparison to neutral drift) is not statistically different from those observed when using humans as a reference (on average, 58–64% across all four major brain regions). We additionally looked at pairwise species comparisons to determine if the general trends of directional selection and differential gene expression were comparable to the human-reference data, and again confirmed that the effects of the reference species used for normalization here are negligible in terms of analyzing differential expression across our primate tree.

Correlation of gene expression to brain size and their change over evolutionary time

Because an increase in absolute brain size is one of the most striking characteristics of humans, we asked what subset of genes is correlated with this important biological trait across the primate tree. Our dataset provides a unique opportunity to evaluate this question, since the average brain size across our study species varies by ~190 fold. Using multivariate analysis (Methods), we defined gene lists most strongly correlated with brain size in each brain region (Supplementary file 8). Results indicate that the same genes rank among those with the strongest positive correlation in PFC, V1, and HIP. CBL also shares some of these same genes but includes more variation among the genes most strongly correlated to overall brain size than the other three brain regions (Figure 5, Supplementary file 9), potentially reflecting the more recent expansion of the CBL relative to the rest of the neocortex (Miller et al., 2019; Smaers et al., 2018).

Figure 5 with 1 supplement see all
Gene correlated with brain size by region.

A clustering and heatmap of the loadings from PC2 of genes for the four regions examined (V1, HIP, PFC, and CBL).

Discussion

We performed RNA-seq on four brain regions from 18 primate species, representing the broadest sampling for any gene expression study in primate brain tissue to date. Through more representative sampling of primate species, we found substantial variation in gene expression levels within the hominoid and strepsirrhine lineages, with the diversity among hominoids particularly impressive due to the recent divergence of this taxon. Using the OU model, we found that a substantial proportion of genes showed differential expression across species. The relative amount of differential expression increased over evolutionary time, both between species and clades. Additionally, when comparing gene expression across broader species and clade comparisons, we observed trends related to brain development, nervous system regulation, and cellular metabolism.

Our findings point to human-specific enrichment for metabolic processes, neural development, and gene regulation. The considerable diversity of gene expression in human and chimpanzee brain tissue has profound implications for understanding the distinct evolutionary processes that have acted upon the brains of the ancestral species of these two lineages. We observed a wide variety of expression dynamics of individual genes in the pairwise comparisons of humans to other primate species. Enrichments for genes under stabilizing selection indicated processes related to transport and cellular localization, while less constrained genes were associated with neuron morphogenesis, plasticity, and cell death. Importantly, gene expression evolution did follow a linear pattern, but did not reach a saturation point due to stabilizing selection pressures as seen in other studies (Chen et al., 2019). Less constrained, neutrally evolving patterns appeared to be the most prevalent pattern in each brain region studied, preventing a saturation point of stabilizing pressures to be reached within Primates with increased phylogenetic distance from humans. Lastly, we identified genes that are correlated with brain size across all major primate taxa, providing candidates for further inquiry.

Our deeper analysis of gene expression has revealed evolutionary patterns that were inaccessible with a more limited sampling of primate brain tissue. We anticipate that the candidate genes and data provided by this study will serve as a resource for many other lines of inquiry into human and non-human primate brain evolution.

Materials and methods

Biological sample collection and RNA extraction

Request a detailed protocol

The sample includes brain tissue from human and nonhuman primates. All samples were obtained from adult individuals free from known neurological disease. If available, the right hemisphere was preferentially sampled. Human brain samples were obtained from the National Institute for Child Health and Human Development Brain and Tissue Bank for Developmental Disorders at the University of Maryland (Baltimore, MD). Chimpanzee brain tissue was obtained from the National Chimpanzee Brain Resource (https://www.chimpanzeebrain.org, supported by NIH grant NS092988). All other sources of brain tissue are listed in Supplementary file 1.

From each individual, we sampled four regions of the brain, including PFC, V1, HIP, and CBL. PFC was sampled from the frontal pole, corresponding to Brodmann’s area 10 in humans. In other primates, the PFC region sampled more broadly encompassed prefrontal cortical areas but was limited to the most anteriorly projecting part of the frontal pole. All V1 samples were dissected around the calcarine sulcus to include primary visual cortex (Brodmann’s area 17). Samples from both PFC and V1 contained all cortical layers and a small amount of underlying white matter (<10%). The HIP was sampled from the medial aspect of the temporal lobe and included all hippocampal subfields. The CBL was sampled from the most laterally projecting region of lateral hemisphere in all primates. In humans, the CBL region corresponded to Crus I or Crus II. CBL samples contained all layers of cerebellar cortex and a small amount of underlying white matter (<10%). Each sample was briefly homogenized using a Tissuelyzer (Qiagen), and the total RNA was isolated using an RNAeasy kit (Qiagen) with a DnaseI treatment.

Library preparation and sequencing

Request a detailed protocol

Single-end RNA-seq libraries were made using the NEBNext mRNA Library Prep Reagent Set for Illumina. Libraries were prepared in batches of 4–8 samples of randomly sampled species and brain regions. Library sizes were checked on the Bioanalyzer (Agilent). RNA-seq libraries were multiplexed on the NextSeq500 (Illumina) in the Genomics Resource Laboratory at the University of Massachusetts Amherst, also randomly distributed across NextSeq500 runs. All fastq files have been submitted to the SRA: https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA639850.

Mapping and transcriptome analysis

Request a detailed protocol

Sequencing reads were assembled into species-specific transcriptomes (containing the reads from all four brain regions) using Trinity (Grabherr et al., 2011). With Trinity assembly, fragments of the original RNA reads for each species are compiled and clustered into groups based on sequence similarity, which eventually are extended to reconstruct full-length transcripts (Grabherr et al., 2011). These transcripts were then blasted against the human blast nt database, with the alignment thresholds for the top hits from different clades listed in Supplementary file 2, similar to the approach in Perry et al., 2012. Individual libraries were then mapped to the species-specific transcriptome using bowtie (Langmead et al., 2009) in RSEM, and count tables were generated using RSEM (Li and Dewey, 2011). Orthology assignments were additionally checked using the Ensemble one-to-one orthology alignments as a guide for the subset of species with a publicly available genome. Transcriptome quality was assessed using BUSCO to determine assembly and annotation completeness (Seppey et al., 2019), Supplementary file 3. We recognize that some of the species transcriptomes show relatively lower BUSCO completeness scores (namely the Slender Loris at 35.8% complete). We hypothesize that this is likely a reflection of the limited tissue sampling in this dataset of only brain tissue. Previous studies have shown that transcriptome assemblies from single tissue regions on average have lower completeness scores than assemblies composed of reads from a variety of tissue types (Simão et al., 2015). This is likely a reflection of tissue-specific gene expression. As evidence of this, we further analyzed reads deemed as ‘missing’ by BUSCO and found that many of these showed little to no expression across the human brain (Uhlén et al., 2015) and Human Protein Atlas proteinatlas.org. We do not make any major conclusions about the Loris species in this manuscript and thus do not believe these BUSCO scores significantly affect the conclusions made in this manuscript. We only consider Loris data in concert with Lemurs, which by comparison have much more complete transcriptomes.

Distance-based data analyses (PCoA and phenograms)

We performed principal coordinates analyses (PCoA) based on a pairwise distance matrix of all 137 samples. The distance matrix was comprised of the top 500 most variably expressed protein-coding genes. Pairwise distances were calculated by leading log2 fold change, providing a symmetrical representation of the expression ratio centered around 0 (i.e. log2(2)=1 while log2(0.5) = –1) (Robinson et al., 2010). Creating the distance matrix and plotting the PCoA were performed using the plotMDS.DGElist function (based in limma) in the edgeR package in R. Although variation is represented across more than 20 axes (Supplementary file 4), the first three axes were plotted to compare patterns across primate taxa and brain region sampled (Figure 2, Figure 2—figure supplements 13). Polygons overlap the data points representing taxa or brain region. The area of each polygon was computed using functions in the sp package of R. The chull() function was used to define the points around the perimeter of each polygon, and the Polygon() function calculated the area of each. Relative areas of each polygon are listed in Supplementary file 5. Figure 3 displays the same data as the PCoA but uses an array of colors allowing the data from each individual species to be visualized.

The same log2 fold change distance matrix was then used to create phenograms representing the similarity of gene expression profiles among samples. The minimum distance neighbor-joining function in the ‘ape’ package of R created a tree based on the method proposed by Saitou and Nei, 1987. The boot.phylo function estimated the reliability of given nodes of the tree by resampling over 1000 iterations. Although our objective in this analysis was to investigate patterns of evolution across the primate order and the brain regions, our sample included multiple individuals from the same species. By treating these samples separately, our analyses represent both within- and between-species variability in gene expression over time.

Analyzing differential expression

Request a detailed protocol

Counts were filtered and normalized using edgeR (Robinson et al., 2010), with any multispecies comparisons using the GLM functionality (McCarthy et al., 2012). Gene Ontology enrichments were performed using the DAVID gene ontology tool 6.8 (Huang et al., 2009a; Huang et al., 2009b) and g:Profiler (Reimand et al., 2016). Supplementary files 6 and 7 shows results from ordered g:profiler enrichments (g:GOSt) performed on DE genes where q<0.05 (note, this is not ranked on polarity of expression, just absolute change) with all genes expressed in this study used as background.

Differential Expression was also analyzed using the package EVEE-tools to incorporate the Ornstein-Uhlenbeck model (Chen et al., 2019). Evolutionary means and variance values were calculated for each gene across the entire phylogeny for individual brain-region datasets as well as the entire RNA-seq bulk dataset in the context of the Ornstein-Uhlenbeck model. Differential expression was determined using a multivariate OU model. Regimes were defined by species/clades of interest (i.e. in human vs. chimpanzee comparisons, the human samples represent one regime in one model to compare to all other samples, represented by a second regime. In comparison, the chimpanzee samples would be represented as an additional model and a separate regime.) p-values were calculated to represent the fit of each OU model (i.e. Human or Chimp-specific expression patterns) against a Brownian motion model. These were then corrected for multiple hypothesis testing using the Benjamini-Hochberg FDR procedure. We used an FDR threshold of 0.05 to define significance. Additionally, Akaike and Bayesian Information Criterion (AIC and BIC) scores were calculated for each gene for each model, and only genes with AIC and BIC scores significant against the null were considered in further analysis. Directionality of differential expression was further determined by comparing estimated mean expression levels for each regime in each model.

We also looked at differential expression beyond pairwise comparisons, again using EVEE-tools, to see how each individual species differed across the entire dataset. With this, one species of interest was treated as a single regime while all 17 other species were grouped as a second regime. The same criteria as above were used to determine significant differential expression, and this data was utilized to construct UpSet plots using the UpsetR package (Conway et al., 2017). We included only relevant species in these UpSet plots (representing each major primate clade) to simplify the graphs. Each graph is also separated by brain region, similar to previous analyses (Figure 3—figure supplement 2).

We also analyzed sources of variation in our data to determine how significant of an effect species differences have compared to other factors, including primate families and individual variation. For this, we conducted a correlation analysis utilizing the non-parametric Mann-Whitney U (MWU) test on the differentially expressed genes between humans and chimpanzees to other primate species. We focused this analysis on only those species with three individuals per brain region (olive baboon, rhesus macaque, and lemurs). Using the set of human-chimpanzee DEGs, we calculated the Spearman rank-based correlation coefficient between each species to either human or chimpanzee expression. We determined whether or not these correlation coefficients were significantly different across species and across brain regions using the Mann-Whitney U test. We determined that there was no significant correlation between any of the three species to human or chimpanzee expression in any of the four brain regions (Figure 3—figure supplement 1; Yapar et al., 2021). This suggests that the expression profiles of these more distantly related primates are equally similar to human and chimpanzee expression patterns. We also conducted additional Analysis of Covariates (ANCOVA) to confirm that other factors, namely age and sex, were not significant sources of variation in our gene expression analyses. These analyses, along with our PCoA plots, show that taxon identity and brain region are the two most significant determinants for DE. Additionally, differences in samples from the same individual are defined by differences at the level of brain region. ANCOVA analyses showed minor residual effects that we deemed as random and were not further analyzed for the purposes of this manuscript.

Outlier expression was determined using the scoreGenes.R script from the EVEE-Tools script suite (Chen et al., 2019). We compared the total dataset-normalized mean expression of a single gene to that of a single select species, in the context of the overall mean expression and variance across the dataset. We specifically analyzed outlier expression in a brain region-specific manner, looking at the datasets subset by individual brain region. To normalize the expression of the entire dataset, a single species was selected to be used as a reference in TMM normalization. For all outlier analyses except for humans, the human samples were used as a reference. For the human outlier analysis, the rhesus samples were used as a reference for normalization. Importantly, we tested the use of different reference species for dataset normalization and did not find a significant difference in the number of outlier genes and the enrichment categories associated. Z-scores were calculated for genes whose expression patterns fit an OU model (in comparison to the null model of Brownian Motion) and whose evolutionary mean is above 5 CPM for the entire dataset (Chen et al., 2019). In order to determine significant outliers, an FDR threshold of 0.05 was again employed, however at this level limited significance was found. Supplementary file 10 shows the results of this analysis for the CBL and PFC brain regions.

Phylogenetic and evolutionary distance analysis

Request a detailed protocol

Categorical enrichments for the contrasts between species and clades are in Supplementary file 6. The phylogenetic tree for primates was downloaded from the UCSC Genome Browser (30 primate species) (Kent et al., 2002). Distances between species were extracted using the Environment for Tree Exploration Toolkit (Huerta-Cepas et al., 2010). The residuals and mean squared expression differences of all orthologs across 18 species were found using the package EVEE (Chen et al., 2019), and in all contrasts, humans were used as the reference species. We then analyzed the subset of genes showing either broadly defined conserved (low) or neutral (higher) variation across species (low <q = 0.05, high q>0.05), with categorical enrichments for these two groups in Supplementary file 7.

Comparisons of the heterogenous tissues used and single-cell gene expression data

Request a detailed protocol

In any study that derives results from homogenized tissue samples, the composition heterogeneity of the samples may drive differences in gene expression (Montgomery and Mank, 2016). To address this issue, we compared the expression of our tissue samples to recent studies that have performed single-cell RNA-seq on neurons and astrocytes. RNA-seq data from primary neurons and astrocytes were obtained from NCBI’s Gene Expression Omnibus (GEO) and processed in the same manner as the tissue samples for all human samples. These included four hippocampal astrocytes, four cortical astrocytes, and one cortical neuron from Zhang et al., 2016 (GEO accession number GSE73721) and three pyramidal neuron samples isolated from an unspecified brain region by the ENCODE project (ENCODE Project Consortium, 2012; Davis et al., 2018) (accession numbers GSM2071331, GSM2071332, and GSM2071418). Only genes with counts greater than zero in all samples and (CPM)>1 in all 23 samples were included in this analysis (n=7111). A PCoA was made from a distance matrix of the top 500 most variably expressed genes by the pairwise biological coefficient of variation (method = “bcv”) across samples (Robinson et al., 2010). Creating the distance matrix and plotting the PCoA were performed using the plotMDS.DGElistfunction in the edgeR package in R. The PCoA of our human samples in comparison to primary neurons and astrocytes suggests that our heterogeneous tissue samples are not biased to contain more neurons or astrocytes as compared to each other (i.e. one tissue is not biased within this small sample set) (Figure 5—figure supplement 1), and is consistent with other neural cell and brain tissue comparisons (Khrameeva et al., 2020).

Additionally, prior knowledge of variation across primates in cell type composition of the brain is informative in interpreting bulk RNA-seq data. It is well known that neuron densities tend to decrease as brain size increases (Sherwood et al., 2020). This suggests that larger brains accommodate a smaller number of neurons per unit volume compared to smaller brains. However, it is interesting to note that other structural elements, such as astrocytes (Munger et al., 2022), microglia (Dos Santos et al., 2020), and synapses (Sherwood et al., 2020), exhibit a relatively invariant density per unit volume across species. Despite changes in brain size, these essential components involved in neural communication and support maintain a consistent presence, emphasizing their crucial role in brain function regardless of the species' brain size. Thus, based on this previous research, in cross-species comparisons of bulk tissue, we can tentatively interpret differences in gene expression to reflect generally similar proportions of major cell types, except for neurons, which are expected to decrease in proportion with larger brain sizes.

Brain size analyses

Request a detailed protocol

Average species endocranial volumes (ECVs) were obtained primarily from Kamilar and Cooper, 2013 and Isler et al., 2008 from the mean of male and female volumes. ECVs were used since reliable brain size data does not exist for all species samples. The data for human ECV (also averaged from male and female data points) was previously published by Coqueugniot and Hublin, 2012. Isler et al., 2008 reported a minimal error when ECV was transformed to brain mass using the correction factor of 1.036 g/ml (Stephan, 1960), and we used this conversion to obtain brain mass estimates from ECVs for all species (Supplementary file 8).

Within each brain region, we performed a PCA on the species average gene expression of the 500 most variable genes by standard deviation, using the prcomp function in R. For each regional PCA, 14 PCs were required to account for about 90% of the variance in gene expression. We performed multiple regression analyses to determine which of the PCs could predict brain size. Using all 14 PCs accounted for at least 95% of the variation in each brain region. Akaike information criterion was applied. However, it was noted that for each brain region, PC2 was the most predictive of brain size by low (regional adjusted R2 values for PC2 against brain size were: PFC, 0.42; V1, 0.56; HIP, 0.50; CBL, 0.36). The 500 genes and their loadings on PC2 are listed in Supplementary file 9. Across the four sampled regions, we find general uniformity in the extent to which individual genes affect brain size, but the cerebellum displays the most unique signature of the regions sampled (Figure 5).

Data availability

Sequencing data have been deposited in the Short Read Archive: BioProject PRJNA639850.

The following data sets were generated
The following previously published data sets were used
    1. Consortium ENCODE Project
    (2016) NCBI Gene Expression Omnibus
    ID GSE78331. single cell RNA-seq from pyramidal cell (ENCLB928LID).

References

    1. Barton RA
    (2012) Embodied cognitive evolution and the cerebellum
    Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences 367:2097–2107.
    https://doi.org/10.1098/rstb.2012.0112
  1. Book
    1. Itõ M
    (2012)
    The Cerebellum: Brain for an Implicit Self
    Upper Saddle River, N.J: FT Press.
    1. Kamilar JM
    2. Cooper N
    (2013) Phylogenetic signal in primate behaviour, ecology and life history
    Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences 368:20120341.
    https://doi.org/10.1098/rstb.2012.0341
    1. Stephan H
    (1960)
    Methodische studien uber den quantitativen vergleicharchitektonischer struktureinheiten des gehirns
    Zeitschrift Fur Wissenschartliche Zoologie 164:143–172.

Decision letter

  1. Jenny Tung
    Reviewing Editor; Duke University, United States
  2. George H Perry
    Senior Editor; Pennsylvania State University, United States
  3. Mehmet Somel
    Reviewer; Middle East Technical University, Turkey

Our editorial process produces two outputs: i) public reviews designed to be posted alongside the preprint for the benefit of readers; ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Decision letter after peer review:

Thank you for submitting your article "Tempo and mode of gene expression evolution in the brain across Primates" for consideration by eLife. Your article has been reviewed by 2 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Patricia Wittkopp as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Mehmet Somel (Reviewer #2).

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential revisions:

1) Control for phylogenetic structure in analyzing gene expression changes across the phylogeny, especially in association with brain size. Strongly consider whether phylogenetically informed analysis methods provide greater insight (and reduce the potential for confounding) over pairwise comparisons, especially given advantages afforded by broad sampling across the primate order.

2) Provide key missing detail on criteria for defining orthologues; demonstrate robustness of results to decisions about cut-offs (e.g., for defining highly variable genes) and cell type compositional heterogeneity between species.

3) Address methodological questions raised by the reviewers regarding multiple hypothesis testing correction and enrichment analyses.

4) Test human-specific changes in expression in relationship to outgroup species as well as chimpanzee.

Reviewer #1 (Recommendations for the authors):

1. Trinity output is extremely noisy and returns many isoforms with poor support, or virtually undistinguishable from each other except for a couple of base pairs here and there, especially when run in a de novo mode. How did the authors prioritise the many isoforms, and determine which are credible and worth analysing further? Similarly, since many of the primates lack a reference genome, how did the authors define the set of 3432 testable one-on-one isoforms, or the 15017 genes testable across hominoids? Is it simply on the basis of pairwise orthology to human using BLAST or…? Did the authors control for gene duplication (eg reciprocal BLAST)? In addition, the thresholds for orthology seem very permissive, on the basis of overall coding sequence conservation amongst mammalian species (eg see supp Figure 1B of Chen et al. 2019, where mean coding sequence identity between humans and non-primate one-to-one orthologs is ~85%).

2. Many conclusions are based on data from the top 500 most variable genes (generally defined on the basis of SD, but not always – the sanity check against scRNA from astrocytes and neurons uses CoV). Why 500? How robust are all of the results presented to this choice of threshold? Is there a particular species driving this variance (on the basis of PCoA results I dare speculate it's chimpanzees and humans)? How do these observations (PCoA, phenograms) compare to those made on the entire dataset of 3432 testable genes across the entire dataset? It seems to me that the number of genes in the whole data set is not so large as to be worth focusing only on an arbitrarily chosen threshold, and that it would be more informative to consider the entire dataset in these analyses. Similarly, when looking at positive selection, the authors only focus on the top 200 genes DE between human and chimpanzee. Why these limitations?

3. The authors test the possibility that compositional heterogenity across samples may be biasing their results (line 723 onwards), which I commend, but nonetheless find somewhat incomplete as it currently stands. The only test performed is a comparison of human bulk RNA samples to a small number of astrocyte and neuron RNA-seq samples, which they say proves their tissues are not biased towards either cell type. But combining bulk RNA and scRNA data is not trivial, and since the PCoA shows samples clustering by technology/study (I presume the outlier neuron comes from the same dataset as the astroctyes), it's unclear how to interpret it. Regardless, there's no mention of what I think is the more interesting source of heterogeneity: is there an expectation that the composition of the same tissue type would vary substantially between species? e.g., an excess of, say, glia in the PFC of humans relative to lemurs or loris that could skew results (totally hypothetical example)? I am not sure if this is possible to taste with existing datasets, but it seems to at least be worth discussing?

4. DE testing: As above, I would like to see more detail in the methods here to better contextualise the results. First, I do not think CPM is the accurate unit for comparison here, since it does not control for differences in gene length between species, which may be substantial at some of the orthology thresholds set by the authors – RPKM might be better suited. Second, how many genes were testable between each of the comparisons summarised by figure 3? Was this testing done using in a pairwise fashion, using all genes testable between the pair of species being compared (as suggested by line 190), or using a single model matrix with many different contrasts to leverage information across the entire data set? Altough the latter represents a substantial trade off in terms of testable genes, might provide additional insights in terms of polarising results and perhaps even pinpointing the emergence of expression divergence through the tree for specific genes or families. Since subsequent analyses focus primarily on genes identified as DE between human and chimp, I think it would be worth delving a bit more here into the broader temporal trends, or by comparing other interesting pairs of primates.

5. Evolutionary mode of expression levels: as mentioned above, I think the authors do themselves a disservice by not examining their data deeply, eg, by not placing results in a broader evolutionary context or clearly distinguishing between genes that appear to evolve neutrally vs those that exhibit other trends (are there any?). I think this is where the potential of the dataset most shines, and so I strongly encourage the authors to examine the EVEE output in greater detail and see if there's anything interesting hiding in there, although I am cognisant that this might fall outside the scope of the manuscript as they see it, and thus leave it up to them to decide what to do. Nonetheless, some possible questions: can the authors tease out genes evolving under positive selection or showing bursts of accelerated evolution from the overwhelming sea of neutrality? While it's obvious why the authors choose to focus on humans, is anything interesting happening in any other primate lineage?

6. Brain size and positive selection: the authors use mean expression within a species for PCA this time around, as opposed to just all available data points – but it seems to me that this might obscure some trends? Again, the reason for focusing on humans in these sections is obvious – the change in brain size between human and chimp is monumental – but I wonder if this obscures more subtle signals in the data.

Reviewer #2 (Recommendations for the authors):

Having commended the authors for compiling this remarkable dataset, I need to share a number of concerns, with regards to data analyses and also interpretation.

1) In general, the analyses and interpretation alternate between investigating general patterns of evolutionary divergence across primate brain regions, and investigating human-specific expression divergence (e.g. Figure 3 is human-based). It may be better to separate the two questions and the analyses used to address them. The authors could start by quantifying the overall phylogenetic signal in brain expression divergence, e.g. using the Chen et al. 2018 model.

In general, I think that performing pairwise DE tests makes little sense with such data (except for cases where specific hypotheses are tested). It would be preferable if the authors studied human-specific expression changes also within an OU-based phylogenetic model that can also incorporate lineage-specific positive selection (e.g. https://doi.org/10.1093/sysbio/syv042).

Also, in analyses about human-specific expression changes, the authors should preferably rely on gene sets which show human-specific upregulation relative to chimpanzees and outgroup species, not just DE gene sets between human and chimpanzee, especially if they wish to interpret the results in the context of other observations (e.g. human-specific positive selection in promoters of semaphorin genes, or increased white matter connectivity).

2) Regarding the analyses on expression patterns correlated with brain size expansions, I would suggest to use some type of phylogenetic residuals analysis, because both brain size and expression will reflect phylogenetic relatedness. So, no surprise that the same genes show correlation across brain regions.

Moreover, it is not clear if the gene list presented in Table 1 and discussed in detail is indeed of statistical significance – multiple testing correction has not been applied.

3) A large number of expression changes may be driven by changes in cell type composition, especially in evolutionary time. At least there should be some discussion on this point in the main text.

Currently the only mention is the supplement, and the content is suboptimal – human brain region bulk transcriptome profiles are compared with cell type-specific transcriptomes of neurons and astrocytes, whereas it is still highly possible that DE genes across species reflect changes in composition. This could be studied explicitly: e.g. https://www.biorxiv.org/content/10.1101/010553v2.

4) Humans are used as reference species in a large number of analyses, but the motivation and rationale behind this is not explained. In fact for questions about general expression divergence in the brain, humans may not be a good reference.

5) Many methodological details are lacking, including crucial information on RNA-seq data quality, how different gene sets were defined, and motivations behind different cutoff choices used, etc. Overall rewriting of methods would be helpful.

[Editors’ note: further revisions were suggested prior to acceptance, as described below.]

Thank you for resubmitting your work entitled "Tempo and mode of gene expression evolution in the brain across Primates" for further consideration by eLife. Your revised article has been evaluated by George Perry (Senior Editor) and a Reviewing Editor.

The manuscript has been improved in the revision process, and both reviewers note their appreciation for your responsiveness to the previous round of reviews. I am therefore returning it so that you can address some final issues regarding methodological clarity identified by Reviewer 1, and to give you the opportunity to fix some errors identified by Reviewer 2. In addition, I noticed that even though you have removed the analysis of positive selection on regulatory regions from the manuscript, it is still referenced in your abstract--so please do a thorough read-through for consistency throughout.

Reviewer #1 (Recommendations for the authors):

I thank the authors for their reply to my comments and the accompanying revisions and additional details – it's good to see this paper again, as I remain impressed by the dataset and by the scope of the questions the authors are hoping to address.

My responses focus primarily on the points they have addressed in this revision, to avoid dragging this process on forever; I have tried hard to ask only for clarifications, or the absolute minimum, when I have asked for something new. As last time, most of my questions center on methods, so I've brought those to the top, but in all cases, they're requests for additional detail, not for things to be repeated or rerun, so despite the length of the comments I hope they're not too onerous.

1. Trinity output: I would still like a bit more detail on this. Did the authors simply take the best scoring match to a given human gene from the blastnt search? Were there any controls for length or similar? The BUSCO scores are a welcome addition, but some of them are quite low, so additional clarity would be good to help understand what was done. Please note that I am not saying 'redo everything with different thresholds,' but I think more detail would be valuable in parsing the surprisingly low number of genes with data – if 50-70% of the transcriptome is evolving under stabilising selection, but that only covers the 3000 testable genes with one-to-one orthology across the clade, that's actually not a lot of genes at all…

2. PCoA etc: I thank the authors for the additional detail, but I confess I am a bit confused as to how they did what they have done. Line 476 (and 570 for the scRNA data) states that the distance metric is based on log2 fold change distances and that it was generated with plotMDS, but that is not what plotMDS does (plotMDS, by the way, is a limma function, not an edgeR function; limma is loaded in the background when edgeR is loaded). As per the Limma manual, "The distance between each pair of samples (columns) is the root-mean-square deviation (Euclidean distance) for the top genes. Distances on the plot can be interpreted as leading log2-fold-change, meaning the typical (root-mean-square) log2-fold-change between the samples for the genes that distinguish those samples."

If this is not what was done and the authors used a different approach to calculating the distance matrix, why was this done, and how should it be interpreted, especially in light of figure 2? Are the trends reported in this figure (and associated supplements) contingent on the distance metric choice?

3. DE testing methods: It's not clear from the text whether any additional covariates (eg known covariates such as sex or age of the individuals, batches even though they were randomised,) were included in the model or whether samples from the same individual were treated as random effects. As above, I am not asking for a redo at this point, but I think it is important to state these details clearly to ensure reproducibility.

4. Outlier detection: It's not clear to me how this works, from the text (especially as it suggests the opposite of that is described in the response to my original comment). The authors first calculate a mean and sd, and then choose humans as a reference species against which to make comparisons. Does this mean that an outlier gene between humans and chimps is an outlier in the chimpanzee lineage? How can it be an outlier in the human lineage if the human is the reference? I do really appreciate the authors for discussing (line 351) the implications of using humans as the reference, given the general evolutionary shifts we naively expect in the human brain, but in that case, can the data please be included as a supplement?

Other questions:

1. Line 182: Since hominoid expression levels are driving so much of the variation in the PCoA analyses, why do the authors think that this isn't showing up in the phenograms?

2. Line 224: Why are the differences between human and siamang treated as human-specific, rather than informative about all great apes? I do not see how this can be disambiguated with the current setup…

3. Figure 3: I think an upset plot or similar here could be an effective way of visualising results, in terms of exploring how far down the phylogeny some signals are shared, but this is optional.

4. Discussion: I think the sentences beginning in 404 and 407 are fundamentally contradictory, and only 407 seems to truly reflect the results above. Is a reference or a word missing from 404?

https://doi.org/10.7554/eLife.70276.sa1

Author response

Essential revisions:

1) Control for phylogenetic structure in analyzing gene expression changes across the phylogeny, especially in association with brain size. Strongly consider whether phylogenetically informed analysis methods provide greater insight (and reduce the potential for confounding) over pairwise comparisons, especially given advantages afforded by broad sampling across the primate order.

We have recalculated DEGs across brain regions in multiple pairwise comparisons using the OU model and EVEE package (Chen 2018). We added these new analyses to the text (starting at line 189) and have updated Figure 3 with these numbers. We also added information in the methods section on the package used and the specific methods.

2) Provide key missing detail on criteria for defining orthologues; demonstrate robustness of results to decisions about cut-offs (e.g., for defining highly variable genes) and cell type compositional heterogeneity between species.

We have put in a BUSCO table and much more extensive methods sections concerning transcriptome builds and ortholog cutoffs.

The reviewer raises an important point regarding the interpretation of bulk tissue RNA-Seq data in the context of cell type compositional heterogeneity across species. We agree that cell type differences can introduce confounding factors and potentially impact the interpretation of gene expression patterns. However, as we also say to the reviewers that brought up this concern:

Prior knowledge of variation across primates in cell type composition of the brain is informative in interpreting bulk RNA-seq data. It is well known that neuron densities tend to decrease as brain size increases (Sherwood et al., 2020). This suggests that larger brains accommodate a smaller number of neurons per unit volume compared to smaller brains. However, it is interesting to note that other structural elements, such as astrocytes (Munger et al., 2022), microglia (Dos Santos et al., 2020), and synapses (Sherwood et al., 2020), exhibit a relatively invariant density per unit volume across species. Despite changes in brain size, these essential components involved in neural communication and support maintain a consistent presence, emphasizing their crucial role in brain function regardless of the species' brain size. Thus, based on this previous research, in cross-species comparisons of bulk tissue, we can tentatively interpret differences in gene expression to reflect generally similar proportions of major cell types, except for neurons, which are expected to decrease in proportion with larger brain sizes.

Furthermore, it is important to note that while bulk tissue RNA-seq has limitations in capturing cell type-specific gene expression patterns, it can still provide valuable insights into the overall gene expression landscape. By comparing bulk transcriptomes across species, we can identify conserved or divergent expression patterns that may highlight important biological processes or evolutionary trends. Additionally, our study serves as a starting point, providing a foundation for future investigations that can employ more advanced techniques, such as single-cell RNA-seq, to delve deeper into cell type-specific expression patterns.

3) Address methodological questions raised by the reviewers regarding multiple hypothesis testing correction and enrichment analyses.

We have chosen to remove the section looking at positive selection in putative promoter regions from the manuscript and the corresponding table.

4) Test human-specific changes in expression in relationship to outgroup species as well as chimpanzee.

We thank the reviewers for this suggestion, and agree that analyses in the context of outgroup species are vital to understanding what is truly unique to humans. To explore this issue, we took Human and Chimpanzee PFC data from EVEE-tools and ran an enrichment analysis. In the EVEE-tools package, based on the Ornstein-Uhlenbeck model of trait evolution, we can look at expression in human and chimpanzee tissue relative to the rest of the dataset, meaning that all other species in the phylogeny are used as an outgroup. With this, the human PFC was enriched for categories related to “neural growth and development”, “metabolic processes”, and “gene regulation” compared to the chimpanzee PFC in the context of the OU model (as opposed to a null model of random drift). We find that both the number of DEGs as well as the specific categorical enrichments are similar using both EVEE and EdgeR for DE analyses, which ultimately suggests that these changes are relevant in the context of an outgroup species. With this, we are confident that these human-specific changes are representative of the phylogeny presented.

We also looked at Human and Siamang PFC data to compare with an outgroup sample, and the same general trends of enrichment are true in this comparison as well, with enrichments for “Neural Growth/Development”, “Gene Regulation”, and “Metabolic Processes”. This supports the idea that this is human-specific and not just a reflection of changes on the chimpanzee branch.

Reviewer #1 (Recommendations for the authors):

1. Trinity output is extremely noisy and returns many isoforms with poor support, or virtually undistinguishable from each other except for a couple of base pairs here and there, especially when run in a de novo mode. How did the authors prioritise the many isoforms, and determine which are credible and worth analysing further? Similarly, since many of the primates lack a reference genome, how did the authors define the set of 3432 testable one-on-one isoforms, or the 15017 genes testable across hominoids? Is it simply on the basis of pairwise orthology to human using BLAST or…? Did the authors control for gene duplication (eg reciprocal BLAST)? In addition, the thresholds for orthology seem very permissive, on the basis of overall coding sequence conservation amongst mammalian species (eg see supp Figure 1B of Chen et al. 2019, where mean coding sequence identity between humans and non-primate one-to-one orthologs is ~85%).

For this and other reviewer questions we have substantially added more details to our Methods section (starting on pg. 12) and moved up more details from the supplement. The testable one-to-one orthologs were taken from the Ensembl database of one-to-one orthologs. For the many species that have published genomes, we also mapped directly to the genome and compared our results to the ones presented here (they were very similar, but more genes are included when mapping directly to the genome, as expected). The thresholds were modeled on another primate transcriptome study (Perry et al. 2012). We have also added BUSCO scores in the Materials and methods (Table S3). The Saki Monkey is the least complete, but we don’t have any results or discussion about that specific lineage and it wasn’t an outlier in any downstream analyses.

2. Many conclusions are based on data from the top 500 most variable genes (generally defined on the basis of SD, but not always – the sanity check against scRNA from astrocytes and neurons uses CoV). Why 500? How robust are all of the results presented to this choice of threshold? Is there a particular species driving this variance (on the basis of PCoA results I dare speculate it's chimpanzees and humans)? How do these observations (PCoA, phenograms) compare to those made on the entire dataset of 3432 testable genes across the entire dataset? It seems to me that the number of genes in the whole data set is not so large as to be worth focusing only on an arbitrarily chosen threshold, and that it would be more informative to consider the entire dataset in these analyses. Similarly, when looking at positive selection, the authors only focus on the top 200 genes DE between human and chimpanzee. Why these limitations?

The top 500 genes is the default state for edgeR and is typical of PCA figures using that package. We also ran the analysis with all of the genes and saw the same patterns. We just show the default for simplicity. The positive selection section has been removed.

3. The authors test the possibility that compositional heterogenity across samples may be biasing their results (line 723 onwards), which I commend, but nonetheless find somewhat incomplete as it currently stands. The only test performed is a comparison of human bulk RNA samples to a small number of astrocyte and neuron RNA-seq samples, which they say proves their tissues are not biased towards either cell type. But combining bulk RNA and scRNA data is not trivial, and since the PCoA shows samples clustering by technology/study (I presume the outlier neuron comes from the same dataset as the astroctyes), it's unclear how to interpret it. Regardless, there's no mention of what I think is the more interesting source of heterogeneity: is there an expectation that the composition of the same tissue type would vary substantially between species? e.g., an excess of, say, glia in the PFC of humans relative to lemurs or loris that could skew results (totally hypothetical example)? I am not sure if this is possible to taste with existing datasets, but it seems to at least be worth discussing?

We agree that mixing the bulk and scRNA is tricky and that there can be confounders. We think that the best approach is to look at the extensive histological literature in primates from the brain.

Prior knowledge of variation across primates in cell type composition of the brain is informative in interpreting bulk RNA-seq data. It is well known that neuron densities tend to decrease as brain size increases (Sherwood et al., 2020). This suggests that larger brains accommodate a smaller number of neurons per unit volume compared to smaller brains. However, it is interesting to note that other structural elements, such as astrocytes (Munger et al., 2022), microglia (Dos Santos et al., 2020), and synapses (Sherwood et al., 2020), exhibit a relatively invariant density per unit volume across species. Despite changes in brain size, these essential components involved in neural communication and support maintain a consistent presence, emphasizing their crucial role in brain function regardless of the species' brain size. Thus, based on this previous research, in cross-species comparisons of bulk tissue, we can tentatively interpret differences in gene expression to reflect generally similar proportions of major cell types, except for neurons, which are expected to decrease in proportion with larger brain sizes.

4. DE testing: As above, I would like to see more detail in the methods here to better contextualise the results. First, I do not think CPM is the accurate unit for comparison here, since it does not control for differences in gene length between species, which may be substantial at some of the orthology thresholds set by the authors – RPKM might be better suited. Second, how many genes were testable between each of the comparisons summarised by figure 3? Was this testing done using in a pairwise fashion, using all genes testable between the pair of species being compared (as suggested by line 190), or using a single model matrix with many different contrasts to leverage information across the entire data set? Altough the latter represents a substantial trade off in terms of testable genes, might provide additional insights in terms of polarising results and perhaps even pinpointing the emergence of expression divergence through the tree for specific genes or families. Since subsequent analyses focus primarily on genes identified as DE between human and chimp, I think it would be worth delving a bit more here into the broader temporal trends, or by comparing other interesting pairs of primates.

Figure 3 was computed using a single model matrix with many different pairwise contrasts. Due to this, the figure does represent fewer testable genes (only 3400 total), but as mentioned above we believe this provides us with the most interesting and relevant points of comparison within the space of this dataset.

We thank the reviewer for the suggestions to look into broader temporal trends, as our expansive dataset allows us the unique ability to better understand some of the more evolutionarily distant primate relationships. In addition to the hominoid comparisons that are the primary focus of the text, we also look at additional species and clade comparisons. Many of these comparisons are shown in Figure 3 and details on the enriched genes in these broader comparisons have been added to the text (lines 250-283).

Significantly, we did not compare all of the species in our dataset in a pairwise manner. This would mean 153 separate differential expression dataset analyses for the whole brain data, and 612 total if we were to analyze on a tissue-specific level. This type of analysis is not within the scope of this paper, and as a result we limited our analyses to only a select few comparisons of interest. Namely, those referenced in the paper for gene categorical enrichment analyses were selected based on the interesting number of differentially expressed genes produced from the EVEE output. However, these other species comparisons can be pursued in future projects and publications.

5. Evolutionary mode of expression levels: as mentioned above, I think the authors do themselves a disservice by not examining their data deeply, eg, by not placing results in a broader evolutionary context or clearly distinguishing between genes that appear to evolve neutrally vs those that exhibit other trends (are there any?). I think this is where the potential of the dataset most shines, and so I strongly encourage the authors to examine the EVEE output in greater detail and see if there's anything interesting hiding in there, although I am cognisant that this might fall outside the scope of the manuscript as they see it, and thus leave it up to them to decide what to do. Nonetheless, some possible questions: can the authors tease out genes evolving under positive selection or showing bursts of accelerated evolution from the overwhelming sea of neutrality? While it's obvious why the authors choose to focus on humans, is anything interesting happening in any other primate lineage?

We thank the reviewer for this comment. We analyzed our dataset for evidence of directional selection and accelerated evolution through a combination of methods that have been added to the manuscript.

Through the EVEE package from Chen 2018, we analyzed the evolutionary history of this dataset through the Ornstein-Uhlenbeck process of continuous trait evolution. With this, we were able to distinguish between gene expression evolving under neutral, stabilizing, and directional selection. We note the percentage of genes that appear to be under stabilizing selection (a majority of the genes in our dataset) (lines 312-326). In contrast to that, we also analyze gene sets that appear to be under directional selection, and thus may be of particular interest to our understanding of primate brain evolution. With this, we utilized two aspects of the EVEE toolset: the ability to score genes for outlier expression (scoreGenes.R) and the ability to identify ‘Differential Expression’ across the phylogenetic tree (ouRegimes.R).

Although the outlier expression script was originally proposed as a method to analyze clinical data, we utilized it here in order to determine the genes that appeared to be the most divergent in terms of expression within our dataset. With this, we expect that genes having significant divergence from the evolutionary mean level of expression and variance (calculated through the EVEE package) would be those that are also under directional selection pressures. We compared our reference sample (human) expression data to that of a test sample, which, in this case, was the mean expression for each individual species comparison. Again, we did not analyze outlier data in all species comparisons across all tissue regions. Rather, we limited this analysis to specific species comparisons of interest. These methods are explained in the manuscript from lines 546-554.

Identification of Differential Expression using the EVEE model employs a multivariate OU model in order to account for multiple regimes of selection (Chen 2018). The methods of this particular model and script are detailed in the manuscript (lines 529-544). With this, we considered “differentially expressed” genes that showed consistent directional shifts to be those likely under directional selection pressures, and thus, evolving under interesting trends (compared to neutral drift, which represented the null hypothesis for this analysis).

The results from the outlier analysis are mentioned in lines 328 to 362, and the results from the differential expression analysis using multiple regimes of selection are detailed from lines 203-283 as well as in Figure 3.

6. Brain size and positive selection: the authors use mean expression within a species for PCA this time around, as opposed to just all available data points – but it seems to me that this might obscure some trends? Again, the reason for focusing on humans in these sections is obvious – the change in brain size between human and chimp is monumental – but I wonder if this obscures more subtle signals in the data.

We kind of have to use ECVs when dealing with this wide of a range of primates since brain size data doesn’t exist that we know of. We’ve added to the methods that we added a human ECV (also averaged of male and female) from data published by (Coqueugniot and Hublin, 2012; https://doi.org/10.1002/ajpa.21655). Isler et al. (2008) reported minimal error when ECV was transformed to brain mass using the correction factor of 1.036g/cc (Stephan, 1960). WE transformed all the ECV to brain sizes using this conversion. This way we can truly say that we were looking for correlations with brain size and not ECV.

Reviewer #2 (Recommendations for the authors):

Having commended the authors for compiling this remarkable dataset, I need to share a number of concerns, with regards to data analyses and also interpretation.

1) In general, the analyses and interpretation alternate between investigating general patterns of evolutionary divergence across primate brain regions, and investigating human-specific expression divergence (e.g. Figure 3 is human-based). It may be better to separate the two questions and the analyses used to address them. The authors could start by quantifying the overall phylogenetic signal in brain expression divergence, e.g. using the Chen et al. 2018 model.

Thank you very much for the suggestion of the model from Chen 2018.

We used the EVEE toolset as a method to look at both tissue-specific expression divergences, as well as for quantification of the overall phylogenetic signal in the brain in regards to this comment.

In the context of the OU model, the overall phylogenetic signal in brain expression divergence was slightly smaller than was observed with EdgeR. This is likely a reflection of the fact that this model is looking at constrained trait expression evolution in the context of phylogenetic relationships, which represents a more stringent type of analysis.

In general, I think that performing pairwise DE tests makes little sense with such data (except for cases where specific hypotheses are tested). It would be preferable if the authors studied human-specific expression changes also within an OU-based phylogenetic model that can also incorporate lineage-specific positive selection (e.g. https://doi.org/10.1093/sysbio/syv042).

As we have noted in other reviewer comments we did DE in the context of the OU model and have included those numbers in the text and have remade Figure 3 with those numbers.

Also, in analyses about human-specific expression changes, the authors should preferably rely on gene sets which show human-specific upregulation relative to chimpanzees and outgroup species, not just DE gene sets between human and chimpanzee, especially if they wish to interpret the results in the context of other observations (e.g. human-specific positive selection in promoters of semaphorin genes, or increased white matter connectivity).

We agree that human-specific DE trends should be analyzed in reference to multiple outgroup species, and we have addressed this point above in a previous reviewer comment. We have also revised statements made about selection in semaphorin genes to look at more broad species comparisons, and these are now in the manuscript.

2) Regarding the analyses on expression patterns correlated with brain size expansions, I would suggest to use some type of phylogenetic residuals analysis, because both brain size and expression will reflect phylogenetic relatedness. So, no surprise that the same genes show correlation across brain regions.

Moreover, it is not clear if the gene list presented in Table 1 and discussed in detail is indeed of statistical significance – multiple testing correction has not been applied.

Based on this concern and those listed below, we have chosen to remove this second section of the manuscript and the table. We address the first part where it’s repeated below.

3) A large number of expression changes may be driven by changes in cell type composition, especially in evolutionary time. At least there should be some discussion on this point in the main text.

Currently the only mention is the supplement, and the content is suboptimal – human brain region bulk transcriptome profiles are compared with cell type-specific transcriptomes of neurons and astrocytes, whereas it is still highly possible that DE genes across species reflect changes in composition. This could be studied explicitly: e.g. https://www.biorxiv.org/content/10.1101/010553v2.

We’ve added a large section describing what was already known about this potential issue to that part of the methods section.

4) Humans are used as reference species in a large number of analyses, but the motivation and rationale behind this is not explained. In fact for questions about general expression divergence in the brain, humans may not be a good reference.

We thank the reviewer for this perspective. In order to investigate this we repeated all of our EVEE-based analyses using two other reference species, siamang and rhesus macaque. In this type of analysis, the reference species refers to a select subset of our entire dataset with which we compare expression differences to. This is used in order to normalize counts in our dataset for downstream analysis.

In the manuscript, we report the percentage of genes that fit better under the model of stabilizing selection as compared to neutral drift. These numbers are not significantly different from the human comparison data and ultimately, a majority of genes still fit better under the OU model. This information has been added at lines 364-387.

We also looked at select pairwise species comparisons with these alternate reference species in order to determine if the general trends of directional selection and DGE were also similar to the human reference dataset. We again saw similar numbers of DEGs, and the enrichment categories of those significant terms also followed the same trends no matter what reference species was used in the analysis. The only terms that would differ would be those that had the highest p-values and thus were the least statistically significant. We have added this to the manuscript at lines 364-387.

5) Many methodological details are lacking, including crucial information on RNA-seq data quality, how different gene sets were defined, and motivations behind different cutoff choices used, etc. Overall rewriting of methods would be helpful.

We agree and have now extensively rewritten and expanded the methods section (starting on pg 12).

[Editors’ note: what follows is the authors’ response to the second round of review.]

The manuscript has been improved in the revision process, and both reviewers note their appreciation for your responsiveness to the previous round of reviews. I am therefore returning it so that you can address some final issues regarding methodological clarity identified by Reviewer 1. In addition, I noticed that even though you have removed the analysis of positive selection on regulatory regions from the manuscript, it is still referenced in your abstract--so please do a thorough read-through for consistency throughout.

We have addressed the reviewers’ comments and have added a optional figure (an Upset Plot now Figure 3 – supplemental figure 2) and a table looking at the effects of which species is chosen as the outlier (Supplemental Table 10). We’re grateful for these thoughtful reviewer comments. Thank you for the catch on the abstract!

Reviewer #1 (Recommendations for the authors):

I thank the authors for their reply to my comments and the accompanying revisions and additional details – it's good to see this paper again, as I remain impressed by the dataset and by the scope of the questions the authors are hoping to address.

My responses focus primarily on the points they have addressed in this revision, to avoid dragging this process on forever; I have tried hard to ask only for clarifications, or the absolute minimum, when I have asked for something new. As last time, most of my questions center on methods, so I've brought those to the top, but in all cases, they're requests for additional detail, not for things to be repeated or rerun, so despite the length of the comments I hope they're not too onerous.

1. Trinity output: I would still like a bit more detail on this. Did the authors simply take the best scoring match to a given human gene from the blastnt search? Were there any controls for length or similar? The BUSCO scores are a welcome addition, but some of them are quite low, so additional clarity would be good to help understand what was done. Please note that I am not saying 'redo everything with different thresholds,' but I think more detail would be valuable in parsing the surprisingly low number of genes with data – if 50-70% of the transcriptome is evolving under stabilising selection, but that only covers the 3000 testable genes with one-to-one orthology across the clade, that's actually not a lot of genes at all…

We thank the reviewers for this comment and are happy to provide more explanation on our Trinity assemblies (as added to the methods section, lines 494-497). With Trinity assembly, fragments of the original RNA reads for each species are compiled and clustered into groups based on sequence similarity, which eventually are extended to reconstruct full-length transcripts. Following Trinity assembly, we mapped our individual samples to each assembled transcriptome using Bowtie2, and from there we further process samples using HTSeq to assign read counts to annotated genomic features (using a GTF file obtained from Ensembl). With this, the BLAST database is never used to directly map to the human genome. Additionally, there are no specific controls for length of transcripts.

Regarding the lower range of BUSCO scores, we would like to highlight the fact that the dataset only consists of brain tissue. Previous studies, including the original BUSCO publication, show that in general transcriptome assemblies from single tissue regions on average have lower completeness scores than assemblies composed of reads from a variety of tissue types (Simao 2015). This is likely a reflection of tissue-specific gene expression. As evidence of this, we further analyzed reads deemed as “missing” by BUSCO and found that many of these showed little to no expression across the human brain (The Human Protein Atlas; Uhlen 2015). This information has been added to the Methods section (lines 505-511).

We would also like to note that the species with the lowest BUSCO scores (and therefore the least complete transcriptomes) are the Slow Loris and Slender Loris. We do not make any major conclusions about the Loris species in this manuscript and thus do not believe these BUSCO scores will be majorly affecting the conclusions made in this manuscript. We only consider Loris data in concert with Lemurs, which by comparison have much more complete transcriptomes (comparable to other primate brain datasets; Simao 2015). This information has been added to the Methods section (lines 511-516).

We are limited to roughly 3000 genes for analysis in order to compare one-to-one orthologs across 18 primates and 70-90 million years of evolutionary history. We recognize that this is a significant limitation and that we are missing a large number of genes in this dataset that may be important to brain evolution on a smaller timescale. However, in order to better characterize long-term primate evolution in these particular brain regions, we believe that the addition of more evolutionary divergent species enhances the significance of the data presented.

2. PCoA etc: I thank the authors for the additional detail, but I confess I am a bit confused as to how they did what they have done. Line 476 (and 570 for the scRNA data) states that the distance metric is based on log2 fold change distances and that it was generated with plotMDS, but that is not what plotMDS does (plotMDS, by the way, is a limma function, not an edgeR function; limma is loaded in the background when edgeR is loaded). As per the Limma manual, "The distance between each pair of samples (columns) is the root-mean-square deviation (Euclidean distance) for the top genes. Distances on the plot can be interpreted as leading log2-fold-change, meaning the typical (root-mean-square) log2-fold-change between the samples for the genes that distinguish those samples."

If this is not what was done and the authors used a different approach to calculating the distance matrix, why was this done, and how should it be interpreted, especially in light of figure 2? Are the trends reported in this figure (and associated supplements) contingent on the distance metric choice?

We thank the reviewer for this method clarification and have added this information about plotMDS to the text (lines 521-525). To build the PCoAs, we are utilizing EdgeR’s plotMDS.DGElist function that is designed for gene expression datasets. Within this (as described in the EdgeR manual), gene counts from a DGElist are converted to log-counts-per-million and then passed on to limma where distances are calculated based on leading log fold changes. To define this further, the leading log fold chance is the “average (root-mean-square) of the largest absolute log-fold-changes between each pair of samples”. (Robinson, McCarthy, and Smyth, 2010)

3. DE testing methods: It's not clear from the text whether any additional covariates (eg known covariates such as sex or age of the individuals, batches even though they were randomised,) were included in the model or whether samples from the same individual were treated as random effects. As above, I am not asking for a redo at this point, but I think it is important to state these details clearly to ensure reproducibility.

We agree with the comment above that information on other sources of variation is important to ensure the conclusions made from this dataset are reproducible. To that end, no additional covariates were included in either the EdgeR DE analysis or the EVEE-tools OU-model based DE analysis. However, we did conduct additional ANCOVA covariate testing to confirm that other factors, including age and sex, were not significant sources of variation in gene expression (The median contribution of age and sex to the gene expression data was about 0.3% for both factors). These analyses, along with our PCoA plots, show that species identity (12% contribution) and brain region (6.6% contribution) are the two most significant determinants for DE. Additionally, differences in samples from the same individual were assumed to be mostly linked to brain region, as both region and individual show similar levels of contribution in the ANCOVA data (6.6% and 5% respectively). These ANCOVA analyses did show minor residual effects that we deemed asrandom and were not further analyzed for the purposes of this manuscript.

This information is added to the methods section (lines 590-596).

Author response table 1 is the data from the ANCOVA analysis (a multi-way ANOVA) performed using the lm function from the stats package in EdgeR.

Author response table 1
termminq25medianmeanq75maxsd
Residuals0.012589560.234783660.311341620.314098940.387743080.752044130.11557519
age8.2272E-100.000491120.0033550.010219050.012965460.191665290.01660386
family0.000189170.1518150.27370520.304220550.435945550.874221170.18970219
individual0.000425910.027881990.053294020.067739060.094400770.377807840.05345294
region0.000215860.034128930.066446260.100444760.128916810.881634520.10373119
sex5.0908E-100.000527580.003136190.009088440.011747990.15012280.01403884
species2.0257E-060.041365670.120981790.19418920.30788350.80053670.18885516

4. Outlier detection: It's not clear to me how this works, from the text (especially as it suggests the opposite of that is described in the response to my original comment). The authors first calculate a mean and sd, and then choose humans as a reference species against which to make comparisons. Does this mean that an outlier gene between humans and chimps is an outlier in the chimpanzee lineage? How can it be an outlier in the human lineage if the human is the reference? I do really appreciate the authors for discussing (line 351) the implications of using humans as the reference, given the general evolutionary shifts we naively expect in the human brain, but in that case, can the data please be included as a supplement?

We thank the reviewers for this comment and have added a more detailed explanation of the methods used for outlier analysis to the text (lines 598-614). First, upon further analysis, we agree that we cannot define outliers for species being used as the reference. With this, we have adjusted the conclusions made in this section of the paper to reflect different methods of testing: for conclusions on “human outliers” we are utilizing rhesus macaque as a reference (lines 343-380).

With this revised methodology, we see outliers in the human prefrontal cortex that include disease-relevant genes, such as Amyloid Precursor Protein (APP) (O'Brien and Wong 2011) and Parkin (PRKN) (Funayama et al. 2023). With this, we also see a select number of metabolism relevant genes, WNT signaling regulators, and members of the neurexin family. We also find outliers in the chimpanzee PFC which include many genes implicated in metabolic pathways, such as ANKH (Szeri et al. 2020) and ETFA (Henriques et al. 2021) (ETFA is also found to be an outlier in humans which suggests that its upregulation is shared within the hominini tribe). This new data has been added to the text (lines 343-380), with citations to support the functional information provided.

We'd also like to clarify the text and better define the term “reference” as the species with which expression values across the entire dataset are TMM-normalized to. In this analysis, we tested two different reference species (human and rhesus) depending on which species we were analyzing for outlier expression. We’ve added in the data for both of these analyses to Supplementary Table 10 as well as further clarified the methods (lines 598-614) for this analysis.

Other questions:

1. Line 182: Since hominoid expression levels are driving so much of the variation in the PCoA analyses, why do the authors think that this isn't showing up in the phenograms?

We think that’s due to the way the phenograms were constructed. We used neighbor-joining for a simple look at the way gene expression was clustering the samples. Since neighbor-joining looks for a minimum sum of branch lengths, and is somewhat weighted by what pairs are grouped in the initial analysis (Felsenstein 2004), we used it more as a simple way to look at how the different brain regions were clustering and less of a metric than the PCoA.

2. Line 224: Why are the differences between human and siamang treated as human-specific, rather than informative about all great apes? I do not see how this can be disambiguated with the current setup…

We refer to these differences as “human specific” due to the nature of the analysis. With the OU model, we are defining our particular regimes for multivariate analysis (Chen et al., 2018). In this case, our first regime of interest is human and the second regime is siamang. We directly are comparing the gene expression patterns for these species, while all other 16 species in the dataset (including the great apes) are included as “background” and used to infer whether stabilizing selection is the most significant evolutionary model. While it is likely that these differences are informative about all great apes relative to the siamang, we did not conduct an individual analysis of each great ape species to the siamang to make this conclusion.

3. Figure 3: I think an upset plot or similar here could be an effective way of visualising results, in terms of exploring how far down the phylogeny some signals are shared, but this is optional.

We thank the reviewer for this suggestion and have added upset plots, separated out by brain region, as Figure 3 – supplementary figure 2. We also reference these in our analyses to show how in general, the number of differentially expressed genes increases over evolutionary time (lines 200-216; methods lines 568-575). We would like to note, however, the way that these upset plots are constructed (using the UpsetR package in RStudio) does not actually allow us to analyze shared signals as the plot only shows the total number of differentially expressed genes, rather than the identity and functional categories of these genes (I.e. 200 DE genes in one species may not be the same 200 in another species). Rather, we find that the gene ontology enrichment analysis conducted using GProfiler is more useful to look for shared signals, and those analyses are included in this manuscript along with the data presented in Figure 3.

4. Discussion: I think the sentences beginning in 404 and 407 are fundamentally contradictory, and only 407 seems to truly reflect the results above. Is a reference or a word missing from 404?

Thank you for noting that! It has been changed and a reference added to make the contrast with another study looking at all mammals. It is now lines 436-438.

https://doi.org/10.7554/eLife.70276.sa2

Article and author information

Author details

  1. Katherine Rickelton

    1. Department of Biology, University of Massachusetts Amherst, Amherst, United States
    2. Molecular and Cellular Biology Graduate Program, University of Massachusetts Amherst, Amherst, United States
    Contribution
    Formal analysis, Visualization, Methodology, Writing – review and editing
    For correspondence
    krickelton@umass.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-6059-7505
  2. Trisha M Zintel

    1. Department of Biology, University of Massachusetts Amherst, Amherst, United States
    2. Molecular and Cellular Biology Graduate Program, University of Massachusetts Amherst, Amherst, United States
    Contribution
    Data curation, Formal analysis, Writing – review and editing
    Competing interests
    No competing interests declared
  3. Jason Pizzollo

    1. Department of Biology, University of Massachusetts Amherst, Amherst, United States
    2. Molecular and Cellular Biology Graduate Program, University of Massachusetts Amherst, Amherst, United States
    Contribution
    Formal analysis, Investigation, Methodology
    Competing interests
    No competing interests declared
  4. Emily Miller

    Department of Biology, University of Massachusetts Amherst, Amherst, United States
    Contribution
    Formal analysis, Visualization
    Competing interests
    No competing interests declared
  5. John J Ely

    1. Department of Anthropology and Center for the Advanced Study of Human Paleobiology, The George Washington University, Washington, United States
    2. MAEBIOS Epidemiology Unit, Alamogordo, United States
    Contribution
    Resources, Writing – review and editing
    Competing interests
    No competing interests declared
  6. Mary Ann Raghanti

    Department of Anthropology, School of Biomedical Sciences, and Brain Health Research Institute, Kent State University, Kent, United States
    Contribution
    Resources, Writing – review and editing
    Competing interests
    No competing interests declared
  7. William D Hopkins

    Department of Comparative Medicine, Michale E. Keeling Center for Comparative Medicine,The University of Texas M D Anderson Cancer Centre, Bastrop, United States
    Contribution
    Resources, Writing – review and editing
    Competing interests
    No competing interests declared
  8. Patrick R Hof

    1. New York Consortium in Evolutionary Primatology, New York, United States
    2. Nash Family Department of Neuroscience and Friedman Brain Institute, Icahn School of Medicine at Mount Sinai, New York, United States
    Contribution
    Resources, Writing – review and editing
    Competing interests
    No competing interests declared
  9. Chet C Sherwood

    Department of Anthropology and Center for the Advanced Study of Human Paleobiology, The George Washington University, Washington, United States
    Contribution
    Conceptualization, Resources, Investigation, Writing – review and editing
    Competing interests
    No competing interests declared
  10. Amy L Bauernfeind

    1. Department of Neuroscience, Washington University School of Medicine, St. Louis, United States
    2. Department of Anthropology, Washington University in St. Louis, St. Louis, United States
    Contribution
    Conceptualization, Investigation, Writing – original draft, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-8518-3819
  11. Courtney C Babbitt

    Department of Biology, University of Massachusetts Amherst, Amherst, United States
    Contribution
    Conceptualization, Formal analysis, Investigation, Writing – original draft, Writing – review and editing
    For correspondence
    cbabbitt@bio.umass.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-8793-4364

Funding

National Science Foundation (BCS-1750377)

  • Courtney C Babbitt

National Institutes of Health (T32 GM135096)

  • Katherine Rickelton

James S. McDonnell Foundation (https://doi.org/10.37717/220020293)

  • Chet C Sherwood

National Institutes of Health (NS-092988)

  • Chet C Sherwood

National Science Foundation (SMA-1542848)

  • Chet C Sherwood

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We would also like to acknowledge our funding from NSF BCS-1750377, the Wenner Gren Foundation, James S McDonnell Foundation (220020293), NSF INSPIRE (SMA-1542848), NIH (NS-092988), and a fellowship to KR from NIH T32 GM135096. https://doi.org/10.37717/220020293.

Senior Editor

  1. George H Perry, Pennsylvania State University, United States

Reviewing Editor

  1. Jenny Tung, Duke University, United States

Reviewer

  1. Mehmet Somel, Middle East Technical University, Turkey

Version history

  1. Preprint posted: April 21, 2021 (view preprint)
  2. Received: May 12, 2021
  3. Accepted: January 25, 2024
  4. Accepted Manuscript published: January 26, 2024 (version 1)
  5. Version of Record published: February 19, 2024 (version 2)

Copyright

© 2024, Rickelton et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 488
    Page views
  • 101
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Katherine Rickelton
  2. Trisha M Zintel
  3. Jason Pizzollo
  4. Emily Miller
  5. John J Ely
  6. Mary Ann Raghanti
  7. William D Hopkins
  8. Patrick R Hof
  9. Chet C Sherwood
  10. Amy L Bauernfeind
  11. Courtney C Babbitt
(2024)
Tempo and mode of gene expression evolution in the brain across primates
eLife 13:e70276.
https://doi.org/10.7554/eLife.70276

Share this article

https://doi.org/10.7554/eLife.70276

Further reading

    1. Chromosomes and Gene Expression
    2. Genetics and Genomics
    Li Sun, Libo Liu ... Quan-wen Jin
    Research Article Updated

    Eukaryotic cells are constantly exposed to various environmental stimuli. It remains largely unexplored how environmental cues bring about epigenetic fluctuations and affect heterochromatin stability. In the fission yeast Schizosaccharomyces pombe, heterochromatic silencing is quite stable at pericentromeres but unstable at the mating-type (mat) locus under chronic heat stress, although both loci are within the major constitutive heterochromatin regions. Here, we found that the compromised gene silencing at the mat locus at elevated temperature is linked to the phosphorylation status of Atf1, a member of the ATF/CREB superfamily. Constitutive activation of mitogen-activated protein kinase (MAPK) signaling disrupts epigenetic maintenance of heterochromatin at the mat locus even under normal temperature. Mechanistically, phosphorylation of Atf1 impairs its interaction with heterochromatin protein Swi6HP1, resulting in lower site-specific Swi6HP1 enrichment. Expression of non-phosphorylatable Atf1, tethering Swi6HP1 to the mat3M-flanking site or absence of the anti-silencing factor Epe1 can largely or partially rescue heat stress-induced defective heterochromatic maintenance at the mat locus.

    1. Chromosomes and Gene Expression
    Qinyu Hao, Minxue Liu ... Kannanganattu V Prasanth
    Research Article Updated

    Out of the several hundred copies of rRNA genes arranged in the nucleolar organizing regions (NOR) of the five human acrocentric chromosomes, ~50% remain transcriptionally inactive. NOR-associated sequences and epigenetic modifications contribute to the differential expression of rRNAs. However, the mechanism(s) controlling the dosage of active versus inactive rRNA genes within each NOR in mammals is yet to be determined. We have discovered a family of ncRNAs, SNULs (Single NUcleolus Localized RNA), which form constrained sub-nucleolar territories on individual NORs and influence rRNA expression. Individual members of the SNULs monoallelically associate with specific NOR-containing chromosomes. SNULs share sequence similarity to pre-rRNA and localize in the sub-nucleolar compartment with pre-rRNA. Finally, SNULs control rRNA expression by influencing pre-rRNA sorting to the DFC compartment and pre-rRNA processing. Our study discovered a novel class of ncRNAs influencing rRNA expression by forming constrained nucleolar territories on individual NORs.