Cellular and molecular signatures of in vivo GABAergic neurotransmission in the human brain

Diverse GABAergic interneuron microcircuits orchestrate information processing in the brain. Understanding the cellular and molecular composition of these microcircuits, and whether these can be imaged by available non-invasive in vivo methods is crucial for the study of GABAergic neurotransmission in health and disease. Here, we use human gene expression data and state-of-the-art imaging transcriptomics to uncover co-expression patterns between GABAA receptor subunits and interneuron subtype-specific markers, and to decode the cellular and molecular signatures of gold-standard GABA PET radiotracers, [11C]Ro15-4513 and [11C]flumazenil. We find that the interneuron marker somatostatin is co-expressed with GABAA receptor-subunit genes GABRA5 and GABRA2, and their distribution maps onto [11C]Ro15-4513 binding in vivo. In contrast, the interneuron marker parvalbumin co-expressed with more predominant GABAA receptor subunits (GABRA1, GABRB2 and GABRG2), and their distribution tracks [11C]flumazenil binding in vivo. These results have important implications for the non-invasive study of GABAergic microcircuit dysfunction in psychiatric conditions.

The GABAergic system comprises diverse interneuron subtypes, innervating different neural targets through a variety of receptors (1). This complexity poses a challenge to microcircuit investigation in the human brain in vivo that can be both selective and non-invasive. GABAergic interneurons vary in their firing threshold, spiking frequency and location of postsynaptic cell innervation, which makes them fit for various functions including the control of synaptic input into the local network and neuronal output regulation (5,12). This multitude of interneurons can be classified through the expression of specific proteins (interneuron markers) (13). While the vast majority of interneurons in the brain are positive for either parvalbumin (PVALB), somatostatin (SST), or vasoactive intestinal peptide (VIP), further specific subtypes can be identified through the expression of other markers, such as cholecystokinin (CCK) (3). This array of interneurons achieves fine-tuned inhibitory responses via the ionotropic GABAA receptor which mediates postsynaptic cell hyperpolarisation. The GABAA receptor (or GABAAR) is a pentameric chloride channel which most commonly comprises two α, two β and one  subunit (14). There are five subtypes of the α subunit and three each of the β and  subunits; moreover, β can be replaced by a  subunit, and  can be replaced by δ, ε or π (14). This generates a great variety of receptors, the biology and pharmacology of which are determined by their subunit composition. For instance, the low-affinity α1 subunit-containing GABAAR (GABAARα1) mediates phasic, or activity-dependent inhibition on the postsynaptic cell, whereas GABAARα5 present higher affinity to GABA, maintaining a more continuous inhibitory tone extrasynaptically (15)(16)(17)(18)(19).
Detailed investigation of the roles that GABAergic interneuron neurotransmission may play in brain function in health and disease requires: 1) precise knowledge of the basic principles underlying the organization of this complex system in the human brain; and 2) the ability to identify and tease apart its specific microcircuits. In this context, positron emission tomography (PET) allows in vivo GABAAR quantification in anatomically defined brain regions through the use of radiolabelled ligands, mainly [ 11 C]Ro15-4513, with high affinity to GABAARα5, and [ 11 C]flumazenil, a benzodiazepine site-specific ligand with affinity to GABAARα1-3 and α5 (20,21). Although receptor affinity for these radiotracers has been confirmed in preclinical research (22), their specificity for defined cell types is unknown. This compromises the understanding of which specific interneuron microcircuits contribute the most to inter-regional differences in signals obtained from human GABA PET measurements. Interestingly, both the distribution of GABAergic interneurons and GABA PET radiotracer binding are heterogeneous across the brain. For instance, SST and PVALB follow an anticorrelated distribution (11), as do the binding patterns of [ 11 C]Ro15-4513 and [ 11 C]flumazenil (23). Moreover, postsynaptic expression of GABAAR subunits, encoded by individual genes in target neurons, is associated with specific interneuron subtypes in rodent (19,(24)(25)(26). Brain-wide gene expression atlases such as the Allen Human Brain Atlas (AHBA) are increasingly being used to gain insight into the mechanisms linking complex brain microcircuits to measurements of human brain function in vivo (27). Thus, determining the spatial relationships between the expression of GABAAR subunits and interneuron markers may inform the basic principles that govern the spatial organization major GABAergic interneuron microcircuits in the human brain. Moreover, because these brain-wide gene expression data can be integrated with neuroimaging measures, such as binding from GABA PET tracers, this approach may help understand which interneuron microcircuits follow and contribute the most to the spatial pattern of binding of these tracers.
Here, we used state-of-the-art imaging transcriptomics to uncover patterns of co-expression between GABAAR subunits and interneuron markers in the human brain, and to decode the molecular and cellular signatures of two gold-standard GABA PET radiotracers, [ 11 C]Ro15-4513 and [ 11 C]flumazenil (28). We demonstrate that SST co-expresses with two GABAAR subunits implicated in affective function, GABRA5 and GABRA2, while PVALB strongly correlates with genes encoding subunits of the most prevalent GABAergic receptor in the brain, GABAARα1β22. While [ 11 C]Ro15-4513 signal covaries with the expression of SST, GABRA5, GABRA2 and GABRA3, [ 11 C]flumazenil signal is positively correlated with the expression of PVALB and genes of the GABAARα1β22 receptor. We also show that VIP is co-expressed with CCK, and that these two genes covary with both radiotracers. Taken together, our findings show for the first time in human that 1) the expression of markers for PVALB and SST interneurons is associated with distinct GABAA receptor complexes; and 2) that the distribution of genes from those two interneuron populations can be tracked by [ 11 C]Ro15-4513 and [ 11 C]flumazenil binding in vivo. Given the key role for PVALB and SST cells in healthy brain function, and the strong implication of PVALB and SST dysfunction in several brain disorders (8), our findings provide a detailed framework to inform future GABA PET studies in the choice of PET tracer to investigate specific GABAergic microcircuitry in health and disease.

Results
GABAergic interneuron markers co-express with specific GABAAR subunits Our first aim was to identify co-expression patterns between interneuron cell-type markers and GABAAR subunits prior to their integration with the PET imaging data. Interneuron markers of interest comprised: the GABA-synthesising enzymes GAD67 (GAD1) and GAD65 (GAD2) (29) Data on all available GABAAR subunits passing quality threshold were included: (α1-5 (GABRA1-5), β1- , ε (GABRE) and δ (GABRD)). Thus, the subunits α6 (GABRA6), π (GABRP),  (GABRQ) and ρ1-3 (GABRR1-3) were not used in further analyses as they did not show levels of expression above background. We performed weighted gene co-expression network analysis (WGCNA) (35) on gene expression data from the AHBA (36). This dataset contained microarray data on 15,633 genes from six post-mortem samples across the left (n=6 healthy donors) and right hemispheres (n=2 healthy donors), which were resampled into 83 brain regions of the Desikan-Killiany atlas (37). WGCNA is a data-driven approach that allows to identify clusters (modules) of highly correlated genes across the whole transcriptome (35).
WGCNA identified 52 co-expression clusters, 13 of which included genes encoding interneuron markers and GABAAR subunits of interest. We selected these 13 clusters to investigate which genes shared cluster allocation ( Figure 1A). SST was located in the same cluster as GABRA5, GABRA2 and GABRB1. PVALB had its own cluster (i.e., it was not located in the same cluster as any other gene of interest). VIP was found in the same cluster as CCK and no other genes of interest. As those three individual clusters included three of the main non-overlapping interneuron markers, labelling the majority of GABAergic inhibitory cells in the mammalian brain (2), we investigated their enrichment in genes co-expressed in specific cell types defined by previous single-cell transcriptomic analysis (38) using the WEB-based GEne SeT AnaLysis Toolkit (39). These analyses revealed inhibitory interneuron cell-type enrichment in the SST and PVALB clusters, and excitatory cell-type enrichment in the VIP cluster (for full results, see Supplement). This indicated that the subsequent analyses including the SST and PVALB clusters were related to GABAergic interneuron cell-types. Other separate clusters of genes included GAD1, GABRA1, GABRB2, GABRG2 and GABRG3; CALB1, CALB2, GABRG1 and GABRE; GABRA4, NPY and TAC3; and GAD2 and NOS1. Finally, GABRB1, GABRB3, GABRA3, RELN, TAC1, GABRD and TAC4 were all found individually in separate clusters that did not share assignment with any other gene of interest.
While we used WGCNA to identify clusters of co-expressed genes, we sought to complement those findings with a pairwise correlation analysis. This served both as a validation step and as a method to investigate co-expression patterns between genes of interest that might not pertain to a discrete co-    (Table 1; for full PLS regression analysis  results, see Supplementary table 1). Interestingly, PVALB expression had a significant negative PLS1 weight (Z = -2.46, pFDR = 0.0255), which suggested an anti-correlation between PVALB and [ 11 C]Ro15-4513 binding.
The radiotracer binding and the distribution of weights resulting from both PLS analyses (cluster-wise and gene-wise) followed and antero-posterior distribution gradient in the brain (Figure 2), consistent with the analogous gradient of SST expression shown previously (11). We then followed up these results with a cell-type enrichment analysis, accounting for the weights associated with each gene     Radiotracer binding, as well as the distribution of weights resulting from both PLS analyses, followed and postero-anterior distribution gradient in the brain (   Statistically significant results (pFDR<0.05) shown only. PLS weight and pFDR shown to third significant figure. PLS, partial least squares regression analysis.

Discussion
Integrating transcriptional and molecular neuroimaging data in humans, we demonstrate that the spatial pattern of expression of specific GABAergic interneuron markers covaries with that of different GABAAR subunits, and that these co-expression patterns explain a substantial portion of the variation in GABA PET radiotracer binding, a measure of in vivo brain neurotransmission. Our main finding is that  (23), resembling a largely developmentally preserved gradient of SST to PVALB distribution in the human brain, as recently reported (11). Using cutting-edge imaging transcriptomics, we show that these findings may not be coincidental, and that specific GABAergic microcircuits may be investigated in humans in vivo with existing neuroimaging methods.
Our approach relies on indirect spatial associations between PET radiotracer binding and gene expression across brain regions. been shown to present 10-15-fold higher affinity to human cloned GABAARα5 than to GABAARα1-3 (23) and the co-expression of GABRA2 and GABRA5 was previously shown by immunohistochemistry in the rat brain (16), which we observed in our analyses as well. Interestingly, we also found covariance between the expression of CCK and GABRA3 and the distribution of [ 11 C]Ro15-4513 signal.
GABAARα2/3 expression in rodents has been most consistently found in the post-synapse of principal cells targeting CCK basket cells (3,19,24,41). It is plausible that the association between [ 11 C]Ro15-4513 binding and CCK cell-related microcircuits we found is circumstantial, if enough spatial overlap between the two microcircuits exists. Alternatively, the discrepancy could be a result of secondary affinity of the tracer to GABAARα2/3, as the latter constitute less than 5% of all GABAARs in the brain (41) and PET radiotracers are administered systemically with a bolus injection. Future research using precise methods such as immunocytochemistry or autoradiography with pharmacological blocking is warranted to determine whether [ 11 C]Ro15-4513 binds to GABAARα2/3 in addition to GABAARα5.
The spatial pattern of [ 11 C]flumazenil binding covaried most strongly with the cluster containing GAD1, GABRA1, GABRB2, GABRG2 and GABRG3. Prior evidence that [ 18 ]flumazenil accumulation across the mouse brain after mutations in α2, α3 and α5 subunits but not in α1 remained similar to that in wildtype mice suggests that flumazenil binding to GABAARα1 predominates in the mammalian brain (42).
Our finding is also in agreement with the notion that the abundance of GABAARα1 in the brain may be reflected in greater GABAARα1 flumazenil binding (14,23). Indeed, GABAARα1β2γ2 is the most widely expressed GABAAR in brain (8,20) and the co-expression of GABRA1, GABRB2 and GABRG2 is supported by analogous observation in preclinical immunohistochemistry studies (43,44). Interestingly, [ 11 C]flumazenil signal was also associated with PVALB expression, consistent with our observation of a high association between PVALB expression with GABRA1, GABRB2 and GABRG2, and in line with preclinical evidence that GABAARα1 is located in synapses formed by PVALB interneurons onto each other and onto principal cells (19,(24)(25)(26)41). Finally, covariance of both GABRA4 and GABRD expression with [ 11 C]flumazenil binding is consistent with the finding that those subunits are commonly co-expressed in the forebrain and that the δ subunit associates extrasynaptically with α1, primarily in the cerebellum but also in the cerebrum (45,46) where [ 11 C]flumazenil uptake is higher than that of Other noteworthy patterns of gene expression covariance and cell-type enrichment with radiotracer signal patterns emerged. For instance, the lack of PVALB assignment to a cluster containing any other genes of interest may reflect the widespread expression of this protein in the human brain.
Additionally, the binding of both radiotracers was associated with CCK and VIP expression. The latter observation may reflect a co-localisation of several cell types within regions showing radiotracer binding, since VIP interneurons primarily innervate SST interneurons (3,26) and VIP/CCK cells were shown to target principal cells and PVALB basket cells (47). On the other hand, the VIP co-expression cluster was enriched in excitatory cell type markers. These observations warrant further investigation with more precise methods such as immunocytochemistry, chemogenetics or pharmacological manipulations.
Our findings have important implications for future studies investigating GABAergic dysfunction with PET. Abnormalities in both PVALB interneuron/GABAARα1 and SST interneuron/GABAARα2,3,5 pathways have been hypothesised to contribute to the pathophysiology of several brain disorders, including those with affective abnormalities (8). For instance, GABAARα2 and GABAARα3 agonism is implicated in benzodiazepine-mediated anxiolysis (14,19 (14,19), [ 11 C]Ro15-4513 may be used in drug discovery paradigms for new anxiolytic medications with more GABAARα2/3 selectivity (51). SST interneuron dysfunction has also been implicated in the aetiology of depressive disorders (41,52), and a preclinical model of SST cell disinhibition produced an anxiolytic-and antidepressant-like effect akin to that of benzodiazepines or ketamine (53). As α5 subunit increases have been found post-mortem in patients with depression (9) Interestingly, people at clinical high-risk for psychosis show increased hippocampal perfusion (59) which is correlated with prefrontal GABA levels, particularly in those individuals who subsequently developed psychosis (60). Taken together, these findings may suggest that SST interneuron dysfunction in the hippocampus and PVALB interneuron abnormalities in cortical regions may play a role in the onset of psychosis. Future PET studies in early psychosis may address this hypothesis by examining PVALB interneuron-mediated networks with [ 11 C]flumazenil, and SST interneuronassociated inhibition with [ 11 C]Ro15-4513.
Our study has some limitations worth acknowledging. First, we relied on indirect spatial associations between gene expression and radiotracer binding, which alone does not directly imply co-expression in the same cell or in interconnected neurons, nor direct radiotracer binding. However, we note that our findings are broadly supported by the preclinical literature using more fine-grained methods, which lends support to the plausibility of our imaging transcriptomics findings in humans. Our findings were generated though several hypotheses, which may guide focused molecular studies in the postmortem human brain, and which can be easily extrapolated to other neurochemical systems.
Moreover, our work follows a similar approach to previous studies investigating relationships between interneuron marker expression and resting-state activity, (11) and benzodiazepine receptor availability (61). Second, the AHBA includes data from six donors only. Samples from the right hemisphere were only collected for two donors, which led us to restrict our analyses to the left hemisphere. Although not a specific limitation of this study, this raises questions about whether this small sample can capture well the principles of organisation of the canonical architecture of gene expression in the human brain and generalise well. Finally, because we applied an intensity threshold to the microarray dataset to minimise inclusion of unreliable measures of gene expression, we were not able to investigate some genes of interest, including GABRA6, GABRP, GABRQ and GABRR1-3, due to the low intensity of signal for these genes in the AHBA dataset as compared to background. Future studies using high sensitivity methods to measures expression of these genes across the whole brain will help to complement our findings in this respect.
In summary, we provide evidence of spatial alignment between the expression of: 1) SST, GABRA5 and GABRA2; 2) PVALB and GABRA1, GABRB2 and GABRG2; and 3) VIP and CCK in human. These findings expand our understanding of the canonical transcriptomic architecture of different GABAergic interneuron microcircuits in the human brain. Furthermore, we provide first evidence that these separate interneuron subtype-specific microcircuits covary with

Methods
The Allen Human Brain Atlas (AHBA) dataset The AHBA dataset includes microarray data of gene expression in post-mortem brain samples from six healthy donors (one female, mean age +/-SD 42.5 +/-13.38, range 24-57) (36). The brain (cerebrum including the brainstem) was sampled systematically across the left hemisphere in all six donors, and the right hemisphere in two of the donors. Manual macrodissection was performed on the cerebral and cerebellar cortex, as well as subcortical nuclei, in 50-200mg increments. Subcortical areas and cerebellar nuclei were sampled with laser microdissection in 36mm 2 increments. RNA was then isolated from these dissections and gene expression quantified with microarray. More information about donor characteristics and dataset generation can be found in the Allen Human Brain Atlas website (https://human.brain-map.org/).

Gene expression data: preprocessing and spatial mapping
Human gene expression microarray data was extracted from the AHBA with the abagen toolbox (https://github.com/netneurolab/abagen) (62) in JupyterLab Notebook through anaconda3 in Python 3.8.5. We mapped AHBA samples to the parcels of the Desikan-Killiany, including 83 brain regions across both brain hemispheres (34 cortical and 7 subcortical regions per brain hemisphere, plus brainstem). Genetic probes were reannotated using information provided by Arnatkeviciute et al., 2019 (63) instead of the default probe information from the AHBA dataset to exclude probes that cannot be reliably matched to genes. According to the existing guidelines for probe-to-gene mappings and intensity-based filtering (63), the reannotated probes were filtered based on their intensity relative to background noise level; probes with intensity less than background in ≥50% of samples were discarded. A single probe with the highest differential stability, ΔS(p), was selected for each gene, where differential stability was calculated as (64): , where ρ is Spearman's rank correlation of the expression of a single probe p across regions in two donor brains, Bi and Bj, and N is the total number of donor brains. This procedure retained 15,633 probes, each representing a unique gene.
Next, tissue samples were assigned to brain regions using their corrected MNI coordinates (https://github.com/chrisfilo/alleninf) by finding the nearest region within a radius of 2 mm. To reduce the potential for misassignment, sample-to-region matching was constrained by hemisphere and cortical/subcortical divisions. If a brain region was not assigned to any sample based on the above procedure, the sample closest to the centroid of that region was selected in order to ensure that all brain regions were assigned a value. Samples assigned to the same brain region were averaged separately for each donor. Gene expression values were then normalized separately for each donor across regions using a robust sigmoid function and rescaled to the unit interval. Scaled expression profiles were finally averaged across donors, resulting in a single matrix with rows corresponding to brain regions and columns corresponding to the retained 15,633 genes.
The genes of interest list included data on all available GABAAR subunits and interneuron markers defined according to the Petilla terminology (13)  indicating co-expression (35). Gene expression correlation matrix was transformed into an adjacency matrix using the soft threshold power of 14. This power value was chosen as it was the first value at which the network satisfied the free-scale topology criterion at R 2 > 0.8, therefore maximising mean network node connectivity (see Supplementary Figure S5). The adjacency matrix was then transformed into a dissimilarity measure matrix, representing both the expression correlation between pairs of genes as well as the number of the genes they both highly correlated with positively ('neighbours') (35). Finally, average-linkage hierarchical clustering using the dissimilarity measure was performed.
Individual modules were identified through the classic 'tree' dendrogram branch cut (65). hemisphere. The analysis was then repeated using the 52 WGCNA module eigengenes as the predictor variables. Prior to each PLS analysis, both predictor and response matrices were Z-scored.

Bivariate correlation of genes of interest
The first PLS component (PLS1) is the linear combination of the weighted gene expression scores that have a brain expression map that covaries the most with the map of tracer binding. As the components are calculated to explain the maximum covariance between the dependent and independent variables, the first component does not necessarily need to explain the maximum variance in the dependent variable. However, as the number of components calculated increases, they progressively tend to explain less variance in the dependent variable. Here, we tested across a range of components (between 1 and 15) and quantified the relative variance explained by each component. The statistical significance of the variance explained by each component was tested by permuting the response variables 1,000 times, while accounting for spatial autocorrelation using a combination of spin rotations for the cortical parcels and random shuffling for the subcortical ones. We decided to focus on the component explaining the largest amount of variance, which in our case was always the first component (PLS1). The error in estimating each gene's PLS1 weight was assessed by bootstrapping, and the ratio of the weight of each gene to its bootstrap standard error was used to calculate the Z scores and, hence, rank the genes according to their contribution to PLS1. The code used to implement these analyses can be found in https://github.com/SarahMorgan/Morphometric_Similarity_SZ.

Funding
This research did not receive any grant from funding agencies in the commercial or not-for-profit sectors. PBL is in receipt of a PhD studentship funded by the National Institute for Health Research

Declaration of interest
The authors declare no competing interests.