Abstract
GABAergic inhibitory interneurons comprise many subtypes that differ in their molecular, anatomical, and functional properties. In mouse visual cortex, they also differ in their modulation with an animal’s behavioral state, and this state modulation can be predicted from the first principal component (PC) of the gene expression matrix. Here, we ask whether this link between transcriptome and state-dependent processing generalizes across species. To this end, we analysed seven single-cell and single-nucleus RNA sequencing datasets from mouse, human, songbird, and turtle forebrains. Despite homology at the level of cell types, we found clear differences between transcriptomic PCs, with greater dissimilarities between evolutionarily distant species. These dissimilarities arise from two factors: divergence in gene expression within homologous cell types and divergence in cell-type abundance. We also compare the expression of cholinergic receptors, which are thought to causally link transcriptome and state modulation. Several cholinergic receptors predictive of state modulation in mouse interneurons are differentially expressed between species. Circuit modelling and mathematical analyses suggest conditions under which these expression differences could translate into functional differences.
Significance Statement
The brain is a complex network of many different cell types. A particularly diverse group of cells is that of inhibitory interneurons, named for their suppressive effect. These interneurons change their activity depending on an animal’s behavioral state—at least in mice. Here, we investigate if this finding generalizes to other species by comparing gene expression patterns of human, turtle, and zebra finch interneurons. Despite sharing an evolutionary past, we find that only human and mouse interneurons have similar gene expression patterns associated with state modulation. A mathematical model suggests which expression differences in individual cells translate into functional differences at the network level.
Introduction
GABAergic inhibitory interneurons are a highly diverse population of multiple cell types (Markram et al., 2004; Tremblay et al., 2016). In recent years, single-cell RNA sequencing (scRNA-seq, Tang et al., 2009) has revealed that these types can be further subdivided into tens of subtypes (Zeisel et al., 2015; Paul et al., 2017; Tasic et al., 2018) that also differ in their morphological and electrophysiological properties (Gouwens et al., 2020; Scala et al., 2021). So far, it has been difficult to understand the functional relevance of this fine-grained diversity. Bugeon et al. (2022) recently bridged this gap by revealing that interneurons show subtype-specific modulation with an animal’s behavioral state, at least in layers 1–3 of mouse primary visual cortex (VISp). Strikingly, this state modulation could be predicted from the first transcriptomic principal component (tPC1). An interneuron types’ tPC1 score also correlated with other dimensions of cellular diversity, such as electrophysiology and connectivity, hinting at an “approximate but general principle” of mouse cortical interneurons (Bugeon et al., 2022).
Intrigued by these findings, we wondered how general the principle embodied by tPC1 actually is (Fig. 1). Are transcriptomic correlates of state modulation similar across different species, or at least across mouse cortical layers and areas? If yes, then this similarity would suggest conserved principles; if no, then the difference could reveal distinct solutions to shared computational problems (Yartsev, 2017; Laurent, 2020; Jourjine and Hoekstra, 2021). The uniformity of interneurons in the mouse brain (Tasic et al., 2018) suggests that their gene expression and state modulation patterns observed in VISp might apply generally. In fact, earlier work by the authors of Bugeon et al. (2022) found that hippocampal interneurons are also organized along a single latent factor (Harris et al., 2018). Similarly, recent comparative transcriptomic analyses have emphasized the conservation of (cortical) inhibitory interneurons across mammals (Hodge et al., 2019; Bakken et al., 2021; Wei et al., 2022; Yu et al., 2023), and more distantly related species (Tosches et al., 2018; Kebschull et al., 2020; Colquitt et al., 2021). But these and other studies (Boldog et al., 2018; Krienen et al., 2020) have also discovered species-specific interneuron subtypes. Additionally, the relative proportions of interneuron types vary even across mouse cortex (Kim et al., 2017; Scala et al., 2019), as does the modulation of interneurons with brain state (Gulledge et al., 2007; Galvin et al., 2018; Bratsch-Prince et al., 2023; Khoury et al., 2023).
Schematic of our main question: could the same transcriptomic axis predict state modulation in other layers and areas of the mouse cortex and in other species? We investigate this by comparing tPCs and cholinergic (ACh) receptor expression across RNA-seq datasets (Tasic et al., 2018; Tosches et al., 2018; Hodge et al., 2019; Bakken et al., 2021; Colquitt et al., 2021; Yao et al., 2021; Bugeon et al., 2022).
We therefore investigated the transcriptomic correlates of state modulation in seven existing single-cell RNA-seq (scRNA-seq) and single-nucleus RNA-seq (snRNA-seq) datasets from the forebrains of mice (Tasic et al., 2018; Yao et al., 2021; Bugeon et al., 2022), humans (Hodge et al., 2019; Bakken et al., 2021), turtles (Tosches et al., 2018), and songbirds (Colquitt et al., 2021). These species each have homologous types of inhibitory interneurons, but their evolutionary history and brain organization vary in important ways. Reptilian cortex, for example, has a three-layered structure akin to the mammalian piriform cortex (Shepherd, 2011; Fournier et al., 2015). Yet, transcription factor expression rules out a one-to-one homology between individual reptilian and mammalian layers or projection cell types (Nomura et al., 2018; Tosches et al., 2018; Tosches and Laurent, 2019). The songbird forebrain, on the other hand, is organized in nuclei instead of layers (Reiner et al., 2004; Jarvis et al., 2005; Güntürkün and Bugnyar, 2016).
We found that transcriptomic PCs show relatively minor differences across smaller evolutionary distances (e.g., between mice and humans) but diverge over longer evolutionary time scales (e.g., mice and turtles). Between-species differences dwarf within-species differences, likely due to biological rather than technical reasons. Specifically, we trace differences in tPCs to species-specific cell type abundances and within-type gene expression patterns. We also find a combination of conservation and divergence in the expression of the cholinergic receptors correlated with state modulation in mice (Bugeon et al., 2022). Circuit modelling predicts the connectivity patterns for which differences in receptor expression translate into species-specific state modulation of interneurons and cortical information flow.
Materials and Methods
Code was written in Python and R and combined into a reproducible workflow using Snakemake (Köster and Rahmann, 2012). The code can be found at Keijser (2024).
Datasets
An overview of the analysed datasets is shown in Table 1. The datasets used for each analysis are listed in the figure captions. An overview of the dataset and preprocessing steps is shown in Table 2. Table 3 lists the relative frequency of different cell types in each dataset.
Overview of analysed datasets
Datasets and preprocessing steps for main data analysis figures
Percentage of cell types for each dataset, rounded to a single decimal place
Replication of Bugeon et al.
The starting point of our replication was the in vivo calcium imaging data and in situ transcriptomic data previously described by Bugeon et al. (2022). We preprocessed and analysed these data following the original paper unless indicated otherwise. We selected interneurons with a high-confidence assignment to a particular subtype (posterior probability at least 0.5, also see Qian et al., 2020) that belonged to a subtype with at least 3 cells. We used the previous assignment into 35 upper-layer subtypes and grouped Serpinf1 cells into the Vip class. Consistent with the original publication, this resulted in 1,065 cells, hierarchically distributed across 5 subclasses (e.g., Pvalb) and 35 subtypes (e.g., Pvalb-Reln-Itm2a). The Npy gene count of 58 cells was missing (NaN); we assumed these values were missing at random and imputed them with the subtype-specific median value. Zero-imputation gave similar results. We computed each cell’s average activity per behavioral state. Whenever a cell was recorded during multiple sessions, we used the session with the longest period of “stationary synchronized” activity since this was the least frequent state. Since 193 cells were not recorded during the stationary synchronized state, state modulation was computed for the remaining 872 cells. The expression matrix contained continuously valued estimates of gene expression instead of integer counts. We normalized these values to 10,000 “counts” per cell for consistency with the other datasets, although this slightly decreased predictive performance. Finally, we log-transformed the normalized values after adding one pseudo-count,
Other datasets
The transcriptomic datasets each consisted of raw count matrices and metadata that included cell class and subtype/cluster. For the Tasic dataset, we only considered the VISp (not the ALM) cells to allow for a direct comparison with Bugeon et al. In the mouse datasets, we merged the small number of Serpinf1 cells into the Vip cluster for consistency with the analyses from ref. Bugeon et al. (2022). From the Colquitt and Tosches datasets, we only used the zebra finch and turtle cells, respectively, since the data from other species (Bengalese finch and lizard) contained only a small number of interneurons. For both datasets, we assigned cells to putative mammalian homologues according to the correlation-based matching in the original publications (Tosches et al., 2018; Colquitt et al., 2021). For the Hodge dataset, we assigned each cell a cell type based on the original publication (Hodge et al., 2019). Cross-species comparisons were based on one-to-one orthologs, identified using the eggNOG-mapper v2 web interface (Huerta-Cepas et al., 2019; Cantalapiedra et al., 2021). To this end, we first downloaded the protein fasta files for Homo sapiens (GRCh38.p14), Mus musculus (GRCm39), Taeniopygia guttata (bTaeGut1.4.pri), and Trachemys scripta elegans (CAS_Tse
Principal component analysis (PCA)
We scaled gene expression values to 10,000 counts per cell (CP10K) to account for differences in sequencing depth across cells and log-transformed the normalized data. We then identified the top 2000 highly variable genes based on their dispersion across cells (Scanpy’s highly_variable_genes; using 3000 genes gave similar results). We computed the top 30 PCs based on these highly variable genes. For visual comparison, we made an arbitrary but consistent choice for the signs of tPC1 and tPC2.
We quantified the similarity of PCs from different datasets using principal angles (Björck and Golub, 1973; Scipy’s subspace_angles). More precisely, let WX be the gene-by-PC matrix whose columns are the PCs of dataset A. The principal angles between the PC subspaces of datasets A and B are computed from the singular value decomposition (SVD) of the PC-by-PC matrix
Subsampling gene counts
The datasets vary in their sequencing depth (the number of RNA counts per cell, see Table 1), presumably due to a combination of technical and biological differences. We aimed to control for these differences by downsampling counts to the depth of the shallower dataset as follows (Fig. 4). Let
Subsampling cells
To investigate the influence of the number of cells, we downsampled the number of cells in each dataset to that of the smallest dataset (640 cells, Table 1). Sampling cells was done without replacement for 10 different random seeds. After downsampling, we performed preprocessing PCA as described above.
Dataset integration
We used Seurat’s anchor-based integration (Stuart et al., 2019) to map datasets onto the Tasic data. To this end, we converted the AnnData objects to Seurat objects. Next, we separately log normalized each dataset as described above (this time using the equivalent Seurat function NormalizeData) and found genes that were highly variable across datasets (FindVariableFeatures, followed by SelectIntegrationFeatures, with 2000 features/genes). Next, we found mutual nearest neighbours across datasets (“anchors”) after projecting each dataset onto the other’s PCA space (reciprocal PCA). A more flexible reduction method (canonical correlation analysis) gave similar results. We then used the anchors to identify and project out dataset-specific differences. After integration, PCA was performed separately on each transformed dataset.
Correlating gene expression with state modulation
To search across all mouse genes, including cholinergic receptors, we correlated the expression from the Tasic data (Tasic et al., 2018) with the state modulation from the Bugeon data (Bugeon et al., 2022). We first averaged the expression and state modulation within each subtype and then computed the correlation coefficient across subtypes (Fig. 4b). To correlate state modulation from the Bugeon dataset with tPC1 from other datasets (Fig. 2g), we computed each tPC1 based on the 72 genes from the Bugeon et al. in situ panel.
Transcriptomic PCs capture conserved and divergent global gene expression patterns. a, Projections of gene expression data from forebrain interneurons onto each dataset’s first tPCs. Samples correspond to cells; colors indicate cell types. Variance explained by tPC1&2: 23.5% (mouse), 22.4% (human), 13.2% (turtle), and 12.6% (zebra finch). b, Schematic: principal angles measure dissimilarity between subspaces spanned by two sets of tPCs (see Materials and Methods). Larger angles indicate larger differences. c, Principal angles between human and mouse subspaces are smaller than between mouse and zebra finch or turtle subspaces. Chance level estimated by sampling random vectors. d Schematic: Variance explained in the mouse data as a measure of tPC similarity. A tPC’s length is proportional to the variance it explains. e, Variance of mouse data explained by tPCs of different datasets, normalized to the variance explained by mouse tPC1. The human, zebra finch, and turtle tPC1 explain 16.2, 4.6, and 1.1% of the variance explained by mouse tPC1. A random direction (dashed line) explains 0.5%. f, Schematic: Predicting mouse state modulation from human tPC1. The color gradient symbolizes the state modulation of mouse cells. g, State modulation of mouse interneurons can be predicted from each interneuron’s projection onto tPC1 from a reference mouse dataset and human tPC1. It can also be predicted, to some extent, from the projection onto zebra finch tPC1 but not turtle tPC1. R2: cross-validated fraction of variance explained, r: Pearson correlation. h, Subsampling procedure to control for the relative abundance of interneuron subclasses across datasets. Colors code for cell types (see panel a). i,j, Matching the relative abundance does not increase the similarity of datasets as measured using principal angles. Numbers indicate order or PCs. Data from refs. Tasic et al. (2018) (mouse), Bugeon et al. (2022) (mouse state modulation), Bakken et al. (2021) (human), Tosches et al. (2018) (turtle), Colquitt et al. (2021) (zebra finch).
Network simulations
We simulated a rate-based network of Pvalb, Sst, and Vip interneurons and excitatory pyramidal neurons. A single equation represented each cell type except for the pyramidal neurons, represented by two equations, modelling the somatic and dendritic compartments. The network state was defined by the rate vector r = (re, rd, rp, rs, rv), of somatic, dendritic, Pvalb, Sst, and Vip activity. The rate of cell type/compartment x evolved according to
The recurrent connections were chosen based on experimental (Pfeffer et al., 2013; Campagnola et al., 2022) and theoretical work (Hertäg and Sprekeler, 2019). The only difference is relatively weak mutual inhibition between Sst and Vip neurons; strong inhibition could prevent the simultaneous activation of these cell types observed in the data (Bugeon et al., 2022).
The interneurons received cholinergic modulation instead of external inputs, and their amplitudes were varied based on qualitative differences in cholinergic receptor expression (Table 4). These amplitudes were the only differences between species-specific networks: to isolate the effect of differences in ACh receptor densities, the synaptic connectivity patterns were unchanged. In the mouse network, Pvalb neurons were negatively modulated; Vip and—to a lesser extent—Sst neurons were positively modulated. This is consistent with both the activity and expression data from the mouse. For the other species, we only have expression data. In the human and turtle network, Pvalb neurons were not modulated, consistent with their weak or absent expression of, e.g., Chrm4. In the turtle network, Sst neurons were negatively modulated; Vip neurons were positively modulated, but to a smaller extent, given the under-expression of Chrm3 and Chrna5 in the turtle data. For the turtle network, we added a positive external input (amplitude 5) to the Sst equation in the absence of cholinergic modulation. A similar result could be obtained by decreasing the dendritic drive during baseline.
Network parameters: species and cell type-specific additive cholinergic modulation
The network dynamics were numerically integrated using a forward Euler scheme with a time step of 0.1 ms. Each simulation consisted of 11,000 time steps divided into a baseline period of 3,300 steps, a cholinergic modulation time of 3400 steps, and another baseline period of 4300 steps. Not shown in Figure 5 is an initial settling time of 300 timesteps. These values were chosen to let the figure highlight the effect of turning the modulation on and off.
Results
We first validated our analysis pipeline by replicating the relevant results from Bugeon et al. (2022) on their data and conducting several additional analyses. Briefly, we reproduced the systematic variation of interneuron subtypes with behavioral state (roughly, running vs. stationary) and its correlation with tPC1. This correlation seems driven by differences within and across cell types (Bomkamp et al., 2019) and is strongest within the Pvalb and Sst populations. We also found that whether interneurons form a continuum or cluster along tPC1 depends on the preprocessing of the transcriptomic data. These caveats aside, our analyses are consistent with those from Bugeon et al. and might reveal similar patterns—or the lack thereof—in other brain areas and species. A detailed description of the replication, along with supplementary figures and mathematical analyses, is available online (Keijser et al., 2024).
Conserved and divergent transcriptomic axes across species
Having validated our approach on mouse data, we next turn to the cross-species comparison using transcriptomic data from humans (Homo sapiens, Bakken et al., 2021), turtles (Trachemys scripta elegans, Tosches et al., 2018), and zebra finch (Taeniopygia guttata, Colquitt et al., 2021); see Table 1 for an overview of all analysed datasets. We compare these data with a large reference dataset from mouse VISp (Tasic et al., 2018).
We first visualized the data from each species. To this end, we preprocessed the datasets using the same analysis pipeline and applied PCA to the resulting RNA count matrices (see Materials and Methods). The projection onto the first 2 tPCs of the human, but not turtle or zebra finch data, was similar to that of the mouse data (Fig. 2a). Mouse and human interneurons clustered by developmental area (Lim et al., 2018), with medial ganglionic eminence (MGE)-born Pvalb and Sst cells occupying one side of tPC1, and caudal ganglionic eminence (CGE)-born Lamp5, Vip, and Sncg cells the other. An intermediate position was occupied by a small group of Meis2 neurons (Tasic et al., 2018), located in the white matter (Frazer et al., 2017). In contrast to the mammalian datasets, the turtle and finch data were characterized by a large population of Meis2-positive neurons (Fig. 2a, Table 3). Transcriptomic and morphological evidence suggests that these cells are likely homologous to neurons in the mammalian striatum rather than the white matter (Tosches et al., 2018; Colquitt et al., 2021).
We quantified these visual differences using the principal angles, which generalize the notion of angle between two lines in a plane (Fig. 2b; see Materials and Methods). Here, we computed the angles between subspaces spanned by each dataset’s top 30 PCs. Consistent with the impression from the first 2 PCs, the principal angles were the smallest between mouse and human subspaces (Fig. 2c). Turtle and zebra finch PCs were both dissimilar to mouse PCs. Principal angles do not require a one-to-one relationship between individual PCs but also do not consider the variance explained by these components. For instance, a pair of highly similar but low-variance dimensions will result in small principal angles—inadvertently suggesting high similarity. We therefore performed a complementary analysis by computing the variance in the mouse data explained by the PCs of other datasets (Fig. 2d). The first human PC accounted for 16% of the variance explained by the first mouse PC; the turtle and songbird tPC1 accounted for 1 and 5%, respectively (Fig. 2e). Each tPC1 explained more variance than a random direction (0.5%), consistent with the existence of some shared global structure.
We confirmed that these results were not due to technical differences in the different datasets. We first controlled for sequencing depth by downsampling gene counts to that of the lowest-depth dataset (Materials and Methods). We also controlled for sample size by downsampling cells to the number of the smallest dataset. In either case, the relatively small differences between mouse and human tPCs were maintained, suggesting this effect is not due to technical factors. We also mapped each dataset onto the mouse data using anchor-based integration (Stuart et al., 2019). This method has been widely used in cross-species analyses (e.g., Hodge et al., 2019; Colquitt et al., 2021; Hain et al., 2022; Keijser and Sprekeler, 2023). As expected, computational integration increased the similarity among the datasets, but, importantly, it preserved the larger similarity between human and mouse data.
How might the transcriptomic differences relate to state modulation? It is possible that even though the top PCs are only similar between mouse and human interneurons, mouse tPC1 does predict state modulation in other species. Because state modulation information was only available for the mouse (Bugeon et al., 2022), we projected this data onto the tPCs from other datasets to determine their predictive ability (Fig. 2f). We validated this approach by showing that tPC1 from a reference mouse dataset predicted state modulation in the Bugeon et al. data (R2 = 0.19). Next, we found that the human tPC1 predicts state modulation in the mouse approximately as well as mouse tPC1 (R2 = 0.18), but the turtle tPC1 did not (R2 = 0) (Fig. 2g). The zebra finch tPC1 showed a weak but significant ability to predict state modulation (R2 = 0.07). We conclude that human tPCs are similar to those of the mouse also on a functional level, in line with evolutionary history. We emphasise that this result relies on statistical similarities at the transcriptomic level and will need to be tested using simultaneously recorded neural and behavioral data.
What evolutionary changes underlie the differences between transcriptomic PCs? At least two non-mutually exclusive processes are possible. First, homologous subclasses could evolve in a species-dependent manner, as indicated by differences in gene expression. Second, evolution can also change the relative abundance of otherwise conserved classes (Bakken et al., 2021; Chartrand et al., 2023). We wondered if the relative abundance of cell classes was sufficient to explain the species differences. To this end, we resampled cells to equal fractions, such that the three classes (Pvalb, Sst, Vip) present in all datasets each accounted for one-third of the cells (Fig. 2h). This increased the visual similarity between the first two tPCs of the mammalian and non-mammalian datasets due to the absence of Meis2 neurons. Still, the matched-abundance datasets were as dissimilar as the original datasets (Fig. 2j). This highlights the divergence of homologous cell types as a driver of evolutionary change in the global transcriptomic landscape.
Similar transcriptomic axes across mouse datasets
The previous cross-species comparison is based on data collected with different sequencing protocols and from different brain areas. To account for these factors, we calibrated the between-species differences against within-species differences by comparing three mouse datasets (Fig. 3a): the in situ data from VISp layers (L) 1–3 (Bugeon et al., 2022), the plate-based (SMART-seq2) data from VISp L1-6 (Tasic et al., 2018), and the droplet-based (10X) data from multiple cortical and hippocampal areas (Ctx & Hpc, Yao et al., 2021).
Similar transcriptomic PCs across mouse data sets. a, Schematic of brain areas sequenced for different mouse datasets. Ctx, cortex; Hpc, hippocampal formation; VISp, primary visual cortex. The datasets were also collected using different technologies (Table 1). b, Similar tPC1&2 across mouse datasets; tPCs1&2 jointly explain 17.7, 23.1, and 29.8% of variance, from left to right. Arrows indicate a qualitative difference: the relative position of Pvalb and Sst cells along tPC1. c, Quantitative comparison with VISp L1-6 dataset based on 2,000 highly variable genes; tPC1 of the Ctx & Hc dataset explains 78.3%% of the variance explained by tPC1 of VISp L1-6. Grey lines: cross-species angles, taken from Figure 2c. d, As c, but based on the 72 genes shared by the three datasets. In this reduced space, tPC1 of the Ctx & Hc and VISp L1-3 explains 88.7 and 81.3%, respectively, of the variance explained by tPC1 of VISp L1-6. e, Relatively small change in Pvalb and Sst position after matching gene sets between L1-6 and L1-3 data. f, Larger differences due to relative cell type abundance. g, Distribution of tPC1 projection of Pvalb (blue) and Sst (orange) cells for the L1-3 data and different versions of the L1-6 data. Match subtypes: select only the Sst subtypes present in the L1-3 dataset. Match genes: select only the genes present in the L1-3 dataset. Match abundance: subsample such that Sst cells comprise only 8% of the samples, as in the L1-3 data. Projections were normalized such that the mean and variance of the respective Pvalb population were zero and one, respectively. Expression data from refs. Yao et al. (2021) (Ctx & Hpc), Tasic et al. (2018) (VISp L1-6), and Bugeon et al. (2022) (VISp L1-3).
Evolution of cholinergic receptor expression a, Schematic of the relationship between tPC1, gene expression, and state modulation. The shaded variable tPC1 is not directly observable. Arrows indicate direct dependence relationships; arrow thickness indicates the (hypothetical) strength of the relationship. b, Correlation between gene expression and state modulation in the mouse data. Gene expression and state modulation were not measured in the same cells and were therefore combined at the level of subtypes (Materials and Methods). Blue: significant correlation (p < 0.05), grey: not significant. Annotated are the cholinergic receptors that significantly predict state modulation. Chrm4 is ranked 327th (top 1.5%) with the strongest negative correlation. Chrm3 and Chrna3, 4 & 5 are among the top 1.7, 2.4, 7.1, and 8.1% with the strongest positive correlation. c, tPC projection of mouse VISp L1-6, colored by Chrna4 expression (log CP10K). Upper layer Sst types were selected before PCA to retain the cell type arrangement of Bugeon et al. (Fig. 3f) d-g, Dotplots of cell type-specific cholinergic receptor expression, z-scored across all cells. h, Schematic of RNA count subsampling to control for differences in sequencing depth. Each RNA count from the deeper dataset was sampled with a probability equal to the relative depth of the deep and the shallower dataset (Materials and Methods). i-k, Log2-fold difference in expression with mouse data after subsampling; negative and positive values indicate under—and overexpression, respectively, compared to mouse data. Each violin plot shows the distribution of 100 subsampled datasets. Differences outside of shaded areas are larger than the typical differences between different datasets of the same species. Without the subsampling procedure, a comparison of raw RNA counts would suggest that every receptor is overexpressed in the mouse by a factor of 8 or more due to the larger sequencing depth of the mouse data. Expression data from refs. Tasic et al. (2018) (mouse), Bugeon et al. (2022) (VISp L1-3), and Tasic et al. (2018) (L1-6).
Circuit model highlights functional consequences of cholinergic receptor expression divergence a, Schematic of mouse cortical circuit model. PC: pyramidal cell. Grey and red sinusoids represent inputs to pyramidal soma and dendrites, respectively; teal arrows represent cholinergic modulation of interneurons. Excitatory and inhibitory ACh receptors are shown in red and blue, respectively. b, Analytical prediction of the Pvalb and Sst ACh receptor expression values for which Pvalb and Sst cells are inhibited by ACh (see online Supplementary Material). The star indicates parameter settings used for simulation in c. At baseline, the PC rate reflects both somatic and dendritic inputs. Cholinergic activation inhibits the dendritic contribution by activating Sst cells. d Human cortical circuit model, in which Pvalb cells lack inhibitory ACh receptors. e, As (b), but with the shaded area indicating expression values qualitatively consistent with the human data. Star: parameter settings used for simulation in f. Cholinergic modulation has a qualitatively similar effect due to the indirect cholinergic inhibition of Pvalb neurons via the Sst neurons. g, Turtle circuit model with inhibitory ACh receptor expression in Sst cells. h, as (b), but with the shaded area indicating turtle expression values. Star: parameter settings used for simulation in i. In contrast to the mammalian circuit models, cholinergic modulation causes a disinhibition of dendritic inputs.
Visually, the projections onto the first tPCs were similar (Fig. 3b), with interneurons clustering by developmental area, as before. But subtle differences were also visible. For example, the L1-3 dataset lacked a Meis2 population present in both L1-6 datasets (Tasic et al., 2018; Yao et al., 2021). The tPC1 score of Sst cells also varied between datasets. In the L1-3 data, Sst cells occupied an intermediate position on tPC1 (and tPC2) compared to Pvalb cells, consistent with their weaker state modulation. In contrast, Sst and Pvalb cells occupied similar positions in the other datasets.
The three datasets were also quantitatively similar: principal angles between different mouse datasets were substantially smaller than angles between species (Fig. 3c). To compare the larger mouse datasets with the smaller dataset of Bugeon et al. (2022), we performed the same analyses after selecting the 72 genes shared by all datasets. This revealed the Ctx & Hpc data to be more similar to the VISp L1-6 data than the VISp L1-3 data (Fig. 3d), consistent with the varying relative positions of the cell types in the space of the first two PCs.
Several factors could explain the different positions of Pvalb and Sst cells along tPC1 and tPC2 (arrows in Fig. 3b). We first tested if the difference was due to layer-specific subtypes known to be transcriptomically identifiable (see, e.g., Tasic et al., 2018; Wu et al., 2023). However, selecting L1-3 subtypes from the L1-6 data followed by PCA only moved the Sst cells further along tPC1 (Fig. 3g, “match subtypes”). We next tested for the influence of gene set by performing PCA on the L1-6 data after selecting 72 genes describing the L1-3 dataset. This only modestly increased the similarity to the in situ data (Fig. 3e), reflecting the careful selection of the gene panel (Qian et al., 2020). Finally, we reasoned that the intermediate position of Sst cells in the L1-3 data could be due to their relative sparsity (8% in the L1-3 data vs 28% in the L1-6 data). After all, a given pattern of covariability explains less variance when present in a smaller number of samples. Indeed, sampling the same number of cells from the entire Sst population moved the Sst population to an intermediate tPC1 position (Fig. 3f). Therefore, the intermediate position of Sst cells in the L1-3 dataset might be due to their relative sparsity.
In summary, mouse datasets are highly similar compared to cross-species datasets despite differences in brain area and sequencing technology. Two human datasets, (Hodge et al., 2019; Bakken et al., 2021) showed equally high levels of similarity. Between-species differences, therefore, likely reflect biologically meaningful signals rather than technical artefacts.
Evolution of cholinergic receptor expression
So far, we have shown that interspecies expression differences are reflected in the first PCs. This rules out a conserved tPC1 that predicts state modulation—at least across evolutionarily distant species. However, it does not rule out that the species-specific tPCs predict state modulation. Unfortunately, this cannot be tested directly due to the lack of data on state modulation for the other species. As a proxy, we, therefore, analyzed the expression of cholinergic receptors that are known to contribute to the correlation between tPC1 and state modulation in mice (Fig. 4a, Bugeon et al., 2022).
According to our analysis (see Materials and Methods), five cholinergic receptors can predict state modulation of upper-layer subtypes in held-out data in mice. These receptors also ranked among the top genes in their correlation with state modulation (Fig. 4b). The predictive nicotinic receptors (Chrna3,4,5) showed a rough gradient along tPC1 (see, Fig. 4c). The only predictive inhibitory receptor (Chrm4), on the other hand, was expressed by Pvalb neurons (Fig. 4d), consistent with their negative state modulation.
Do the same receptors mediate state modulation in other species? If yes, one would expect differential expression across cell types, with a similar pattern as in mice. However, several receptors that predict state modulation in mice show qualitatively different patterns of expression in the other species (Fig. 4e–g). For example, Chrna4 and Chrna5 show much weaker expression in the human data than in mice (Fig. 4e). Chrm4 is overexpressed in the turtle data relative to the other species (Fig. 4f).
The general trend is that the predictive receptors are under-expressed in the other datasets. A possible explanation is a regression to the mean: predictive receptors are, by necessity, expressed in mice. But the relative expression in other datasets could also be due to technical reasons such as a lower sequencing depth (Table 1). For example, the typical mouse cell contained several orders of magnitude more detected RNA counts than the typical human nucleus (Fig. 4h). We controlled for this confound by downsampling the mouse data to the sequencing depths of the other datasets (Materials and Methods). To measure variability, we also applied this procedure to two datasets from the same species, which revealed typical log2-fold expression differences between −2 and 2. In comparison, downsampling retained larger differences between species that are qualitatively consistent with the analysis of the full datasets. In the human data, Chrm4, and to a lesser extent Chrna3 & 4, were still underexpressed after downsampling (Fig. 4i). Chrm3 and Chrna5 were underexpressed in the turtle data, whereas Chrm4 was overexpressed (Fig. 4j). In the songbird data, only Chrm4 was overexpressed (Fig. 4k).
Thus, several cholinergic receptors that might mediate state modulation in mice show species-specific expression. This suggests that homologous cell types in different species could show substantial differences in state modulation.
Robustness of state modulation to cholinergic receptor expression
How do species-specific cholinergic receptor expression patterns influence cortical information flow? Since this depends not just on the cell type-specific gene expression but also on the interplay of different interneurons, we investigated this question using a circuit model (Materials and Methods). We focused on the most salient differences in receptor expression between the three species with a cortex (mice, humans, and turtles).
The model consists of the three most common interneuron types, Pvalb, Sst, and Vip cells, whose connectivity patterns have been mapped (Pfeffer et al., 2013; Jiang et al., 2015) and are relatively conserved, at least in mice and humans (Campagnola et al., 2022). Additionally, the computational repertoire of this “canonical circuit” has previously been investigated (Litwin-Kumar et al., 2016; Yang et al., 2016; Hertäg and Sprekeler, 2019). The microcircuit connectivity in the turtle (and zebra finch) is, to our knowledge, currently unknown. Our models are therefore a proof-of-principle, to be extended when more data become available. As a parsimonious working hypothesis; here, we simulated a turtle circuit assuming similar connectivity patterns as in the mammalian cortex and mathematically analysed which changes in connectivity could lead to functional differences (see Keijser et al., 2024). To explore the effect of cholinergic modulation on excitatory activity, we also included a two-compartmental pyramidal neuron. These two compartments receive different information streams: whereas the soma receives feedforward (sensory) input, the pyramidal dendrites receive top-down input (Larkum, 2013; Ju et al., 2020; Schroeder et al., 2023). For visualization purposes, these input streams were represented by sinusoids of different frequencies (Fig. 5a).
We model cholinergic modulation as an additive input to the interneurons with a strength that is based on cell type-specific receptor expression data (Materials and Methods, Table 4). Cholinergic modulation inhibits the Pvalb population via muscarinic ACh receptors (Chrm4) while activating Vip cells and—to a lesser extent—Sst cells via nicotinic ACh receptors. The activated Sst cells suppress inputs arriving at the dendrites (Silberberg and Markram, 2007; Murayama et al., 2009; Keijser and Sprekeler, 2022), and increase the effective cholinergic inhibition of Pvalb cells. By recruiting dendritic inhibition, ACh therefore limits the influence of top-down inputs on PC rates (Parikh et al., 2007; Guillem et al., 2011; Fig. 5c) while enhancing the influence of top-down inputs (Hasselmo, 1995; Moran et al., 2013).
How might cholinergic modulation affect interneuron activity and information flow in other species? We modelled a human cortical circuit by deleting the inhibitory Chrm4 receptor from the Pvalb population, mimicking the strongest difference with the mouse VISp data (Fig. 5d). Mathematical analyses indicate that Pvalb activity is still reduced in the presence of ACh due to inhibition from Sst cells (Fig. 5e,f; see Materials and Methods and online Supplementary Material). The difference in direct cholinergic inhibition in Pvalb interneurons between humans and mice might therefore have a relatively weak functional consequence. This is consistent with experimental data (Chen et al., 2015) and the variable Chrm4 expression across mouse datasets. By contrast, we found a qualitatively different effect after changing the ACh receptor densities to mimic the turtle data (Fig. 5g): the over-expression of inhibitory Chrm4 receptors by Sst instead of Pvalb cells led to cholinergic disinhibition of dendritic inputs (Fig. 5h,i). This qualitative deviation from the mammalian cortex might affect the state-dependent processing of different input streams in reptiles (see Discussion).
Although the lack of data regarding turtle (or songbird) connectivity precludes us from making testable predictions, these findings highlight that the computation performed by the circuit can be very sensitive to certain patterns of differential expression and robust to others (Prinz et al., 2004; Marder et al., 2014).
Discussion
We have shown that the global gene expression patterns of inhibitory interneurons, as assessed by PCA, show considerable similarity between mice and humans. However, such similarity is not observed between mice and turtles or songbirds. This suggests that the first transcriptomic PC obtained from the upper layers of the mouse cortex (Bugeon et al., 2022) may serve as an organizing principle for mammalian interneurons but not for reptilian and avian interneurons. Control analyses indicate that technical factors cannot explain cross-species differences. Instead, our results suggest that the evolutionary divergence of homologous interneurons is mainly explained by changes in gene expression rather than changes in the relative abundance of cell types. Alongside the differences in global expression patterns, we have also observed differences in the expression of cholinergic receptors, suggesting that interneurons undergo species-specific modulation in their functional states.
Comparison to prior work
The large transcriptomic differences between species might be surprising given the evolutionary conservation of broad interneuron subclasses in the forebrain (Tosches et al., 2018; Bakken et al., 2021; Colquitt et al., 2021; Wei et al., 2022) and other areas (Kebschull et al., 2020; Yu et al., 2023). However, these works also found many genes with species-specific expression, suggesting cell types might be homologous across species but not preserved in their detailed properties. Moreover, fine interneuron subtypes are not necessarily conserved across larger evolutionary distances (Tosches and Laurent, 2019). Our analyses reveal global differences in gene expression between homologous cell types that rival the differences between glutamatergic cells.
Cholinergic fluctuations with arousal and other cognitive processes have been reported in several species (see, e.g., Xiang et al., 1998; Disney et al., 2007; Pinto et al., 2013; Puzerey et al., 2018; Jaffe and Brainard, 2020; Zhu et al., 2023), and might even be mediated by similar midbrain cell types (Lovett-Barron et al., 2017). However, a causal link between acetylcholine and arousal-like states has only been shown in mice (e.g., Alitto and Dan, 2013; Pinto et al., 2013), and acetylcholine seems to act via species-specific pathways (Coppola and Disney, 2018). For example, most human but not rat PV neurons express the Chrm1 receptor (Disney and Reynolds, 2014). Even within a single species, cholinergic projections and their effects vary across areas (Avendano et al., 1996; Li et al., 2018; Khoury et al., 2023) and layers (Poorthuis et al., 2013; Obermayer et al., 2017). Although the data therefore suggest a broadly conserved role for acetylcholine in state modulation, this assumption requires testing using functional experiments. In addition, technical differences between datasets naturally pose a more serious limitation when comparing the expression of individual genes than aggregate measures such as PCA. Future work will need to confirm the present results using, e.g., immunohistochemistry.
Interpretation of transcriptomic PCs
Across all datasets, one feature consistently stands out: the clustering along tPC1&2 by developmental area. A similar pattern has been previously observed based on nonlinear dimensionality reduction and clustering methods (see, e.g., Tasic et al., 2018; Keijser and Sprekeler, 2023). The structuring of top PCs by developmental origin and cell type is expected since cell types are defined by developmentally activated transcription factors that coregulate batteries of protein-coding genes (Hobert, 2008; Arendt et al., 2016). These low-dimensional patterns of gene expression are naturally picked up by a method like PCA.
A clear difference between the mouse datasets is given by two layer-specific subtypes: the deep-layer Meis2 cells (Frazer et al., 2017; Tasic et al., 2018), and long-range projecting Sst Chodl cells (Tasic et al., 2018; Melzer and Monyer, 2020). Their intermediate position along tPC1 but distinct connectivity suggests that the correlation of tPC1 with cellular properties (Bugeon et al., 2022) might not apply to these deep layer subtypes. This could be further tested using, for example, Patch-seq experiments (Cadwell et al., 2016; Lipovsek et al., 2021). A caveat is the relative sensitivity of PCA to cell-type proportions: the intermediate tPC1 scores of Meis2 cells may be caused by their scarcity in the mouse data.
Evolution of cholinergic modulation
Many genes are at least as predictive of state modulation as cholinergic receptors (Fig. 4b). Some genes might be causally related to state modulation, but others are merely co-regulated with causal genes (Fig. 4a). Thus, strong correlations—the very property that allows the reliable identification of transcriptomic PCs—also preclude the identification of causal genes based on regression analyses. Whether genes predict state modulation also depends on factors only partially under genetic control, such as synaptic connectivity. For example, our network analyses show that synaptic connectivity patterns influence which cholinergic receptor expression differences affect state modulation (Fig. 5): sufficiently strong inhibition from Sst cells causes Pvalb cells to be inhibited by ACh, even if Pvalb cells do not express cholinergic receptors. Finally, species differences in state modulation may not necessarily imply differences in function. An interesting example is turtle Sst cells’ expression of the inhibitory Chrm4 receptor, which might lead to cholinergic disinhibition of pyramidal dendrites. Since sensory inputs to turtle cortex arrive in layer 1 instead of deeper layers (Reiner, 1993; Ulinski, 2007), we speculate that acetylcholine can thus disinhibit sensory inputs, as it potentially does in mouse cortex. Alternatively, cholinergic modulation could have qualitatively different effects on the processing of somatic and dendritic inputs in turtle compared to mouse cortex. Future experiments could arbitrate between these alternatives.
Limitations
Our main findings seem robust against potential confounding factors such as sample size and sequencing depth. Yet, computational controls cannot correct all technical differences between the datasets. Some differences may therefore be due to experimental rather than biological factors. Hence, a cross-species comparison using data collected with similar protocols is an important topic for future work. We used theoretical modelling to explore the functional implications of expression differences. At this stage, a lack of experimental constraints prevented us from making specific predictions regarding species-specific effects of internal states. We expect this to change with the availability of more data from nonstandard model organisms.
Conclusion
The wide availability of transcriptomic data in different species offers new opportunities for comparative analyses. Transcriptomic data can not only predict behavioral features such as state modulation but also the electrophysiology and morphology of homologous cell types (Tripathy et al., 2017; Gouwens et al., 2020; Bernaerts et al., 2023, but see Scala et al., 2019) which are more accessible. It will be exciting to see whether these predictions generalize across species and if they correlate with high-variance transcriptomic dimensions. More generally, we expect that future cross-species experiments will complement work in the genetically accessible mouse to reveal general principles of brain function.
Footnotes
We thank Trygve Bakken, Michael Brainard, Stéphane Bugeon, Bradley Colquitt, Kenneth Harris, and Bosiljka Tasic for their feedback. We are also grateful to Sadra Sadeh for suggesting a modelling approach.
- Correspondence should be addressed to Joram Keijser at keijser{at}tu-berlin.de, joramkeijser{at}gmail.com.
