Functional associations among G protein-coupled neurotransmitter receptors in the human brain

Background The activity of neurons is controlled by groups of neurotransmitter receptors rather than by individual receptors. Experimental studies have investigated some receptor interactions, but currently little information is available about transcriptional associations among receptors at the whole-brain level. Results A total of 4950 correlations between 100 G protein-coupled neurotransmitter receptors were examined across 169 brain regions in the human brain using expression data published in the Allen Human Brain Atlas. A large number of highly significant correlations were found, many of which have not been investigated in hypothesis-driven studies. The highest positive and negative correlations of each receptor are reported, which can facilitate the construction of receptor sets likely to be affected by altered transcription of one receptor (such sets always exist, but their members are difficult to predict). A graph analysis isolated two large receptor communities, within each of which receptor mRNA levels were strongly cross-correlated. Conclusions The presented systematic analysis shows that the mRNA levels of many G protein-coupled receptors are interdependent. This finding is not unexpected, since the brain is a highly integrated complex system. However, the analysis also revealed two novel properties of global brain structure. First, receptor correlations are described by a simple statistical distribution, which suggests that receptor interactions may be guided by qualitatively similar processes. Second, receptors appear to form two large functional communities, which might be differentially affected in brain disorders.


Background
A typical neuron receives thousands of synaptic contacts [1], and each postsynaptic site can express a number of neurotransmitter receptors. Since neurons integrate signals from all receptors on their surface, their activity is determined by receptor sets and not by individual receptors.
Functional receptor groups explain why a constitutive null-mutation of a neurotransmitter receptor often produces a mild phenotype, even when the receptor is known to be important in specific brain functions [2][3][4][5][6][7][8][9][10]. Such minor functional effects can be explained by compensatory mechanisms in the developing brain, which at least partially depend on receptors that detect other neurotransmitters [4,11].
Similarly, the absence of an entire neurotransmitter may not have a major effect on brain development and function. Serotonin (5-hydroxytryptamine, 5-HT) is an extremely abundant neurotransmitter in the brain: by some estimates, the density of serotonergic varicosities in the rat cerebral cortex is around 6·10 6 /mm 3 , with each cortical neuron receiving some 200 varicosities [12]. The density of serotonergic projections may exceed that of the brain capillary system [13,14] and must carry a significant energetic cost. However, genetic mutations and pharmacological manipulations that eliminate virtually all 5-HT in the brain produce only subtle behavioral alterations, with no gross morphological or cellular changes [15][16][17][18]. Abnormally low dopamine levels in the brain can be nearly asymptomatic until around 50-80% of the substantia nigra neurons are lost [19], and ablation of dopaminergic neurons in neonatal rats does not result in any significant motor dysfunctions [20,21]. A lack of norepinephrine due to a genetic mutation in the dopamine-β-hydroxylase gene produces an unremarkable neurological phenotype in humans [22,23]. Mammalian thalamic nuclei can function normally in the virtual absence of GABAergic interneurons, as has been shown in mice and rats [24,25]. Experimental evidence shows that the lack of one neurotransmitter can be compensated for by changes in other neurotransmitters: for example, serotonergic processes have been shown to permanently extend into brain areas previously occupied by dopaminergic terminals [26], and serotonin and/or dopamine may compensate for the lack of norepinephrine [27,28]. It should be noted that at least some adult neurons have the flexibility to switch from one neurotransmitter to another in response to the environment [29].
These findings suggest that biological information in the brain is coded not by individual neurotransmitters or their receptors, but by finely-tuned neurotransmitterreceptor sets. While this hypothesis does not require that the receptors be physically linked, it is supported by the unexpected abundance of heteromeric receptor complexes [30,31]. For example, the serotonin 5-HT 2A receptor can form complexes with the metabotropic glutamate mGluR 2 receptor [32] and the dopamine D 2 receptor [33], which may play a major role in the action of antipsychotic drugs and hallucinogens. A heterocomplex is more than the sum of its individual receptors, since heteromerization can alter receptor mobility at the neuron surface, downstream signaling, and intracellular trafficking [32,34].
Understanding neurotransmitter receptor sets will require new analytical and theoretical approaches. Functionally complete receptor sets have to be isolated and their dynamic properties investigated in specific brain regions and in individual cells [35]. While these studies pose technical challenges, they are likely to lead to major theoretical simplifications. Some researchers have already used this approach with considerable success [36][37][38]. As a further step toward this goal, the present study used the Allen Human Brain Atlas [39,40] to examine mRNA expression associations among nearly all known G protein-coupled neurotransmitter receptors in the human brain.

Results
The analysis used the mRNA expression data of 100 G protein-coupled receptors (Table 1) in 169 regions of six normal human brains presently available in the Allen Brain Atlas database ( Figure 1). The mRNA amounts of many receptor pairs were very strongly correlated ( Figure 2). The five strongest positive and negative correlations of each receptor are given in Table 2. The distribution of all 4950 correlations had a nearly symmetric shape, with a single mode shifted toward the positive values ( Figure 3A, B). This distribution failed normality tests ( Figure 3A; the Kolmogorov-Smirnov test: D = 0.024, p < 0.01; the Shapiro-Wilk test: W = 0.995, p < 0.01), but was well described by the beta distribution with the same range, mean and variance (the shape parameters α = 3.51 and β = 3.32) ( Figure 3B; D = 0.015, p = 0.209). The distributions of several functionally meaningful subsets were not significantly different from the complete set (in each pair, both receptors represent the same neurotransmitter ( Figure 3C): D = 0.076, p = 0.169); (in each pair, the receptors represent different neurotransmitters ( Figure 3D): D = 0.004, p = 1); (in each pair, both receptors have the same G protein-coupling ( Figure 3E): D = 0.015, p = 0.950); (in each pair, the receptors have different G protein-couplings ( Figure 3F): D = 0.007, p = 1)).
All statistically significant correlations among the receptors were plotted as a graph ( Figure 4A). The nociceptin receptor (ORL-1) had the highest vertex degree (connectivity), due to its significant correlations with 58 receptors ( Figure 4B). The largest clique consisted of 18 completely interconnected vertices: the glutamate receptors mGluR 2 , mGluR 4 , mGluR 7 , the adrenergic receptors α 1B , α 1D , β 1 , the serotonin receptors 5-HT 1A , 5-HT 1F , 5-HT 2A , the cholinergic receptors M 1 , M 3 , the histamine receptors H 1 , H 2 , the melanin-concentrating hormone MCH 1 , the neuropeptide Y receptors Y 1 , Y 5 , the nociceptin receptor (ORL-1), and the somatostatin receptor SST 1 . The distribution of the vertex degrees (the number of links originating in each receptor) appeared to be bimodal and did not follow the power law that is often observed in natural networks with high functional connectivity ( Figure 4A, inset). However, a recent study has shown that a bimodal degree distribution can emerge in robust networks [41].
Next, the obtained correlation information was used to examine whether some receptors groups are more tightly interlinked than other receptors or, more presicely, whether the graph ( Figure 4A) can be broken down into distinct receptor communities. Two community detection methods were used: the modularity algorithm and the clique percolation algorithm [42]. The modularity method revealed two major receptor communities ( Figure 5). As recommended by Palla et al. [42], the clique percolation method was optimized using several correlation thresholds (T = 0.5, 0.6, 0.7, 0.8, and 0.9) and k-cliques of several sizes. The best separation was achieved with T = 0.6 and k = 4, which again revealed two distinct receptor communities ( Figure 6). With the exception of one receptor (BDKRB2), the separation among the receptors was identical to the one obtained with the modularity method. Since the clique percolation method used more stringent criteria, it excluded some more weakly correlated receptors (importantly, they were not placed in the opposite community). These analyses suggest that the human brain has two functional receptor communities, within each of which the mRNA levels are strongly correlated and can potentially affect each other. The first community contains

Discussion
The Allen Human Brain Atlas is a relatively new database [39,40] that continues to be updated and refined. In the absence of generally accepted standards for how gene expression data should be normalized and presented, the published expression values should be treated with caution. The expression of many receptors is highly consistent across individuals, but some receptors show a high degree of variability (Table 1). Among them is the dopamine receptor 4 (coded by the DRD4 gene), which has been extensively studied because of a functionally important polymorphism in its exon 3 [43][44][45][46][47]. The origin of its inconsistent distribution across individuals is not clear and may be due to either unreliable detection or true expression differences in the population. Notably, no relationship has been established between the DRD4 alleles and their mRNA levels [48]. Caution should also be exercised in the interpretation of mRNA levels that show a consistent pattern across the subjects. For many receptors, the relationship between the mRNA and protein quantities is often poorly understood, and a change in one of these measures may not indicate a change in the other. A recent large-scale study has shown that, on average, mRNA levels explain around 40% of the variability in protein levels, and that the abundance of a protein is predominantly controlled by translation [49]. In addition, many neurotransmitter receptors operate in two different signaling modes: at the membrane surface through G-proteins and, when internalized, in an arrestin-dependent fashion [50]. Neurotransmitter receptor genes can produce several mRNA splice variants, some of which may be constitutively active (ligand-independent) [51], or have different internalization properties [52]. Also, protein molecules can be phosphorylated, glycosylated, and undergo other modifications [50,[53][54][55]. These post-translational processes place severe limitations on functional inferences from mRNA data. On the other hand, mRNA quantification allows a high degree of specificity, which remains difficult to achieve in protein analyses. Post-translational modifications of protein molecules and the absence of specific antibodies for a number of neurotransmitter receptors (contrary to the claims of manufacturers) currently do not allow large scale quantifications of proteins in the entire brain. Even though the protein data remain limited, Table 1 provides information about the possible inter-individual variability of the mRNA levels of nearly all neurotransmitter receptors and will facilitate the interpretation of completed and future studies.
The analyzed mRNA levels in brain structures reflect the cumulative gene expression in several types of neuronal and glial cells, with a possible contribution from endothelial and ependymal cells. This lack of spatial precision makes the obtained results too "coarse" for the modeling of local neural circuits. However, an association between the abundance of two receptors over many brain regions is functionally meaningful, just as biologically The inter-subject correlations were obtained from 6 subjects (15 cross-correlations).
meaningful information can be obtained from the correlation between the population numbers of two species across geographic areas (even if the species do not directly interact). Hypofunction or hyperfunction of a receptor in a class of cells is likely to affect the activity of their local neuroecosystems, which may induce changes in the expression of receptors in other cells [56]. Therefore, estimates of the most likely associations among receptors are important for the interpretation of receptor knockout models, as well as for the prediction of changes in other receptors associated with pharmacological targeting of a specific neurotransmitter receptor [57,58]. At present, no comprehensive quantitative analysis exists to facilitate these theoretical considerations, and published data are likely to be biased by hypothesis-driven approaches, funding agency priorities, and attracting nodes in researcher networks. While the current analysis is a step forward, major theoretical breakthroughs can be expected when technical capabilities become sufficiently advanced to dynamically monitor entire receptor sets in single neurons and glial cells [35]. Currently little information is available about the correlation between the mRNA levels of receptors that form heteromeric complexes. Among them, the complex between the metabotropic glutamate receptor mGluR 2 and the serotonin receptor 5-HT 2A has been particularly well studied, partly because of its potential importance in schizophrenia and other related brain disorders [32,59,60]. It has been recently reported that the disruption of 5-HT 2A receptor-dependent signaling can suppress mGluR 2 transcription through epigenetic modifications in the mGluR 2 gene promoter [61]. The present analysis found a highly significant positive correlation between the mRNA levels of these two receptors (0.49). However, mGluR 2 had the highest positive correlations with the nociceptin receptor (ORL-1), the adrenergic α 1D and β 2 receptors, the histamine H 2 receptor, and the purine P2Y 6 receptor; and the 5-HT 2A receptor had the highest positive correlations with the histamine H 1 receptor (Figure 1), the serotonin 5-HT 1F receptor, the muscarinic cholinergic M 3 receptor, and the melanin-concentrating hormone receptors MCH 1 and MCH 2 ( Table 2).
Two receptor communities were extracted from the data ( Figures 5 and 6). There are many neural circuits where these receptors interact, but it remains unclear whether the entire communities can be assigned a biologically Figure 1 The mRNA expression profiles of six neurotransmitter receptors. The horizontal axis represents the 169 brain regions and the vertical axis represents the mRNA amounts (averaged across the probes and subjects). Since the numerical mRNA values are normalized and relative, they are not indicated on the vertical axis (for all six genes, it ranges from −2 to 2). Note high similarity between some of the profiles (e.g., mGluR 1 and CRF 1 , 5-HT 2A and H 1 ). meaningful role. It should be noted that many of the receptors in the two "minimal" communities ( Figure 6) control global brain functions, such as wakefulness and sleep [62]. The two communities can be differentially affected in some brain disorders. For example, the 5-HT 1A , 5-HT 2A and 5-HT 4 receptors belong to the same community ( Figure 6) and are expressed by neurons in the medial prefrontal cortex (mPFC) that project to the dorsal raphe nucleus and control serotonin release [63,64]. Altered activity of these neurons has been implicated in mood  Figure 2 The scatter plots of the six most strongly correlated receptor pairs (out of the 4950 pairs). Each point represents the mRNA amounts (averaged across the probes and subjects) of the two receptors in one of the 169 brain regions. Considerable clustering is apparent, which indicates that in most brain regions the mRNA levels of the two receptors are either low or high simultaneously, with few regions in between. All correlations are significant at p < 10 -15 .     65,66]. The exact structure of the communities is likely to become more refined as more data become available in the Allen Brain Atlas. In general, receptor network analyses hold great promise for understanding the brain in health and disease, as has been demonstrated by recent research [36,[67][68][69].
The relatively simple distribution of correlations ( Figure 3) raises interesting questions. Theoretically, such a distribution can be obtained from a single dynamical interaction. Depending on the numerical values of its coefficients, the same process can produce uncorrelated or highly correlated equilibrium values, even if the two receptors are strongly dynamically coupled [70]. Since theoretically all receptors can be expressed by all brain cells and they can only differ in their equilibrium levels (some mRNA amounts may be undetectably small), such qualitative uniformity remains an intriguing possibility.

Conclusions
Progress in neuroscience requires both accurate factual observations and complexity reduction. Since information processing in the brain depends on thousands of unique neurotransmitter-receptor interactions, understanding how these neurotransmitter-receptor pairs operate in functional groups is not only a theoretical imperative, but also a practical necessity. The obtained results suggest that the apparent complexity of neurotransmitter signaling has an underlying global structure, which is not readily detectable if receptor interactions are studied in isolation.

Methods
The human brain expression data of one hundred G protein-coupled neurotransmitter receptors (Table 1) were downloaded from the Allen Brain Atlas data portal (http:// human.brain-map.org; the data release of March 7, 2013). Technical details about the brain donors, tissue preparation, specificity controls, and data normalization (including normalization across brains) are described in the Allen Human Brain Atlas Technical White Papers (Case Qualification and Donor Profiles, Microarray Survey, Microarray Data Normalization).
The normalized mRNA amounts in 169 brain regions were obtained from six available subjects (three Caucasian males (31, 55, and 57 years of age), two African-American males (24 and 39 years of age), and one Hispanic female (49 years of age)).
Of the analyzed regions, the 14 regions from the myelencephalon were the central glial substance, the arcuate nucleus of the medulla, the inferior olivary complex, the gracile nucleus, the cuneate nucleus, the raphe nuclei of the medulla, the central medullary reticular group, the lateral medullary reticular group, the gigantocellular group, the hypoglossal nucleus, the dorsal motor nucleus of the vagus, the spinal trigeminal nucleus, the vestibular nuclei, and the cochlear nuclei. The 12 regions from the pontine tegmentum were the pontine nuclei, the superior olivary complex, the central gray of the pons, the paramedian pontine reticular formation, the locus ceruleus, the nucleus subceruleus, the pontine raphe nucleus, the medial parabrachial nucleus, the lateral parabrachial nucleus, the facial motor nucleus, the abducens nucleus, and the trigeminal nuclei. The 27 regions from the cerebellum were 12 vermal areas (I-II, III, IV, V, VI, VIIAf, VIIAt, VIIB, VIIIA, VIIIB, IX, X), 11 lobules (III, IV, V, VI, Crus I, Crus II, VIIB, VIIIA, VIIIB, IX, X), and the four deep cerebellar nuclei (the fastigial nucleus, the globose nucleus, the emboliform nucleus, and the dentate The receptor numbers (Table 1) are given in the left column and the correlations in the right column (e.g., mGluR1 has the strongest positive correlation with CRF1 (0.80)). Correlations with an absolute value of less than 0.34 are not significant after the Bonferroni correction (at α = 0.05) and should be interpreted with caution. The actual correlations may be somewhat stronger due the imperfect reliability of the expression data (the statistical "attenuation" phenomenon). For compactness, the following abbreviations were used: GABA1 = GABABR1, GABA2 = GABABR2.  nucleus). The 14 regions from the mesencephalon were the ventral tegmental area, the substantia nigra, the red nucleus, the central gray of the midbrain, the midbrain raphe nuclei, the midbrain reticular formation, the trochlear nucleus, the oculomotor nuclear complex, the Edinger-Westphal nucleus, the inferior colliculus, the superior colliculus, the pretectal region, the interstitial nucleus of Cajal, and the nucleus of Darkschewitsch. The 11 regions from the thalamic area were the subthalamus, the ventral thalamus, the posterior group of nuclei, the medial geniculate complex, the dorsal lateral geniculate nucleus, the dorsal division of the lateral group of nuclei, the ventral division of the lateral group of nuclei, the anterior group of nuclei, the medial group of nuclei, the caudal group of intralaminar nuclei, and the rostral group of intralaminar nuclei. The 19 regions from the epithalamus and the hypothalamus were the pineal, the habenular nuclei, the paraventricular thalamic nuclei, the posterior hypothalamic area, the lateral hypothalamic area, the mammillary region of the lateral hypothalamic area, the mammillary body, the supramammillary nucleus, the tuberomammillary nucleus, the tuberal region of the lateral hypothalamic area, the lateral tuberal nucleus, the perifornical nucleus, the ventromedial hypothalamic nucleus, the dorsomedial hypothalamic nucleus, the anterior hypothalamic area, the arcuate nucleus of the hypothalamus, the preoptic region, the paraventricular nucleus of the hypothalamus, and the supraoptic nucleus. The 9 regions from the striatum, pallidum, and septum were the head of the caudate nucleus, the body of the caudate nucleus, the tail of the caudate nucleus, the nucleus accumbens, the putamen, the external segment of the globus pallidus, the internal segment of the globus pallidus, the substantia innominata, and the septal nuclei. The 6  Figure 5 The two communities of G protein-coupled neurotransmitter receptors detected with the modularity method.
regions from the amygdala were the central nucleus, the basomedial nucleus, the cortico-medial area, the basolateral nucleus, the lateral nucleus, and the amygdaloid transition zone. Two other regions from the lateral pallium were the claustrum and the piriform cortex. The 7 regions from the hippocampal formation were the parahippocampal gyrus, the dentate gyrus, the CA1 field, the CA2 field, the CA3 field, the CA4 field, and the subiculum. The cingulate gyrus was subdivided into the frontal, parietal, and retrosplenial parts, and the insula was subdivided into the short and long gyri. The 9 regions from the temporal lobe were the temporal pole, planum polare, the transverse gyri, Heschl's gyrus, the planum temporale, the superior temporal gyrus, the middle temporal gyrus, the inferior temporal gyrus, and the fusiform gyrus. The 6 regions from the occipital lobe were the occipital pole, the cuneus, the lingual gyrus, the superior occipital gyrus, the inferior occipital gyrus, and the occipito-temporal gyrus. The 6 regions from the parietal lobe were the precuneus, the posterior paracentral lobule, superior parietal lobule, the angular gyrus, the supramarginal gyrus, and the postcentral gyrus. The 19 regions from the frontal lobe were the paraterminal gyrus, the subcallosal gyrus, the parolfactory gyri, the anterior paracentral lobule, the precentral gyrus, the superior frontal gyrus, the middle frontal gyrus, the frontal operculum, the opercular part of the inferior frontal gyrus, the triangular part of the inferior frontal gyrus, the orbital part of the inferior frontal gyrus, the lateral orbital gyrus, the medial orbital gyrus, the posterior orbital gyrus, the anterior orbital gyrus, the gyrus rectus, the superior rostral gyrus, the inferior rostral gyrus, and the frontal pole. In addition, the expression data from the cingulum, the corpus callosum, and the choroid plexus of the lateral ventricle were used.
Expression data were available from four or more subjects in 86% of the 169 brain regions, and 53% of the 169 brain regions were represented by all six subjects. Only four brain regions (2%) were represented by a single subject. The median number of mRNA probes per gene was 3. One gene was analyzed with only one probe (ADORA2B) and one gene was analyzed with 89 probes (CNR1). Probes were located on different exons as much Figure 6 The two communities of G protein-coupled neurotransmitter receptors detected with the clique percolation method. The absolute values of correlations were thresholded at 0.6 (with the corresponding p < 10 -15 ) and the clique number was set at k = 4.