Complex Control of GABA(A) Receptor Subunit mRNA Expression: Variation, Covariation, and Genetic Regulation

GABA type-A receptors are essential for fast inhibitory neurotransmission and are critical in brain function. Surprisingly, expression of receptor subunits is highly variable among individuals, but the cause and impact of this fluctuation remains unknown. We have studied sources of variation for all 19 receptor subunits using massive expression data sets collected across multiple brain regions and platforms in mice and humans. Expression of Gabra1, Gabra2, Gabrb2, Gabrb3, and Gabrg2 is highly variable and heritable among the large cohort of BXD strains derived from crosses of fully sequenced parents—C57BL/6J and DBA/2J. Genetic control of these subunits is complex and highly dependent on tissue and mRNA region. Remarkably, this high variation is generally not linked to phenotypic differences. The single exception is Gabrb3, a locus that is linked to anxiety. We identified upstream genetic loci that influence subunit expression, including three unlinked regions of chromosome 5 that modulate the expression of nine subunits in hippocampus, and that are also associated with multiple phenotypes. Candidate genes within these loci include, Naaa, Nos1, and Zkscan1. We confirmed a high level of coexpression for subunits comprising the major channel—Gabra1, Gabrb2, and Gabrg2—and identified conserved members of this expression network in mice and humans. Gucy1a3, Gucy1b3, and Lis1 are novel and conserved associates of multiple subunits that are involved in inhibitory signaling. Finally, proximal and distal regions of the 3′ UTRs of single subunits have remarkably independent expression patterns in both species. However, corresponding regions of different subunits often show congruent genetic control and coexpression (proximal-to-proximal or distal-to-distal), even in the absence of sequence homology. Our findings identify novel sources of variation that modulate subunit expression and highlight the extraordinary capacity of biological networks to buffer 4–100 fold differences in mRNA levels.


Introduction
There is an extraordinarily high level of variation in expression of messenger RNAs of key inhibitory GABA type-A receptors (GABA(A)R) in human brain. For example, GABRA1-a key constituent of the great majority of functional receptors-varies over a 100-fold range in the neocortex of a sample of 187 healthy but elderly humans [1]. In comparison, in the same sample, the GABRG2 and GABRB1 receptors ''only'' vary 10-20 fold. This is a remarkable range that exceeds that which is often achieved in knock-in and knock-down experiments in genetically engineered lines of mice.
This variation in expression is doubly remarkable because dysregulation of GABA(A)R have been linked to a wide range of abnormalities and neurological diseases, including epilepsy, autism, impulsivity, substance abuse disorders, mood, psychiatric disease, and chronic pain. This raises an important question as to the causes and consequences of the high level of endogenous variation among normal humans. Is it a technical artifact of arraybased methods? Is it due to difficulties in obtaining high quality RNA from human brain? Or does it reflect flexible use of GABA(A)R subunits to assemble pentameric receptors?
To answer these questions we need to weigh the relative importance of genetic, environmental, and technical sources of variation. What fraction of variation is heritable and what fraction is due to environmental or technical error? What is the functional relevance of these differences at the RNA level? In the large sample of humans studied by Webster and colleagues (2009), each individual has a unique genotype and it is not practical to resample the same genotype many times to estimate or eliminate technical errors. However, these questions can be addressed efficiently using diverse sets of fully inbred strains of mice that model human populations.
We have exploited a large set of inbred strains of mice-the BXD family-to study the expression of GABA(A)R subunits in the brain. This family is composed of both parent strains (C57BL/ 6J and DBA/2J) and ,160 lines. While each line is fully inbred, the entire collection is highly diverse and members are segregating for ,5 million sequence variants. Both parents have been fully sequenced and the progeny have been genotyped at over 7,000 genetic markers. Furthermore, members of this family have been well phenotyped for ,40 years. This remarkably dense data allows us to define genetic and phenotypic differences between parents and among their progeny. As in humans we detect large differences in subunit expression among members of this family. For example, a2, a4, and b3 subunits vary more than 3-fold even after averaging multiple samples per strain. Here we address four key questions about the molecular genetics of GABA(A)R subunits: N What are the sources of variation in subunit mRNA expression?
N What genomic regions (QTLs), genes, or even sequence variants contribute to subunit variation? N To what extent do subunits interact with each other or with other genes critical in synaptic function at the mRNA level?
N What is the functional impact of this variation in terms of behavior and disease susceptibility?

Expression, variation, and heritability
We examined expression of GABA(A)R subunits across six brain regions, two species, and multiple expression platforms using a large number of probes and probe sets, and RNA sequencing (RNA-seq). To simplify the presentation much of the analysis is based on expression in the hippocampus and on data generated using the Affymetrix M430 and Exon 1.0 ST platforms (Table S1). Expression values are expressed on a log2 intensity scale. Values of 8 to 9 units correspond to moderate levels of expression (2 to 4 pM), whereas values of 12 to 13 correspond to high levels of expression (,40 to 80 pM).
A majority of the GABA(A)R subunits are well expressed in the hippocampus. The most abundant subunits include Gabrb1, Gabra5, Gabra1, Gabrg2, and Gabrb2. Other subunits are more moderately expressed (Gabra2, Gabra3, Gabra4, Gabrb3, Gabrg1, and Gabrd) or expressed at very low levels (Gabrg3). Several subunits, including Gabra6, Gabrr1, Gabrr2, Gabrr3, Gabre, Gabrp, and Gabrq are not expressed above background in the hippocampus and are not considered in detail further. Extensive data is available for all of these genes and their corresponding probe sets in GeneNetwork (www.genenetwork.org). As shown in Table 1 and Figure S1, there is a good correlation between mean subunit expression on the arrays, gene level expression based on RNA-seq of whole brain, and the intensity of in situ expression in the Allen Brain Atlas [2].
The expression of a subset of subunits is markedly different among the BXD strains. For example, the expression level of Gabra1, Gabra2, Gabra4, Gabrb2, and Gabrb3 differs by more than 2.5-fold among strains (Table 1). This high variation in mRNA levels is attributable to a combination of technical, environmental, and genetic factors. Some of the variation could be due to hybridization artifacts caused by polymorphisms that overlap probe target sequences [3,4,5]. To eliminate signals resulting from this artifact we removed probes overlapping SNPs and small indels (Methods , Table S1). We are confident that the high level of variation shown in Table 1 is not an artifact. To estimate the influence of genetic versus environmental modulation we calculated the heritability of mRNA expression using the large hippocampal data set (Table S2). For many subunits, heritability ranged from 0.40 to 0.95-demonstrating that genetic factors are a major contributor to variation in expression. This heritable variation enables us to identify genetic loci or so-called expression quantitative trait loci (eQTLs) that control subunit expression using genomic extensions of classical gene mapping methods introduced by Damerval and colleagues [6]. These eQTLs can be classified broadly into two types: cis eQTLs and trans eQTLs. A cis eQTL is caused by sequence variants located within or very close to the gene from which mRNA is transcribed. In contrast, a trans eQTL is caused by sequence variants located on a different chromosome or more than 10 Mb distant from the cognate gene on the same chromosome.

Modulation by cis eQTLs
In agreement with their high heritability, Gabra1, Gabra2, Gabrb2, Gabrg2, and Gabrb3 are strongly modulated by cis eQTLs that have likelihood ratio statistic (LRS) scores above 20 ( Figure 1). (An LRS score can be converted to a logarithm of the odds (LOD) score by dividing by 4.61; where LOD is approximately equal to 2log(P) of linkage.) Three subunits-Gabra1, Gabrb2, and Gabrg2are located close together on chromosome (Chr) 11 between 33 and 42 megabases (Mb) and the remaining two subunits, Gabra2 and Gabrb3, are located within GABA(A)R subunit clusters on Chr 5 (near 71 Mb) and Chr 7 (near 65 Mb), respectively. Higher expression of all three genes on Chr 11 is generally associated with the B allele inherited from the maternal C57BL/6J strain whereas higher expression of Gabra2 and Gabrb3 is always associated with the D allele inherited from the DBA/2J paternal strain.
We tested the consistency of the cis eQTLs initially detected in hippocampus from the BXD population at four different levels; (1) brain tissue, (2) expression platform, (3) mRNA region, and (4) mouse population.
Comparison of subunit expression across multiple CNS tissues (whole brain, cerebellum, striatum, cortex, hypothalamus, and amygdala) and several array platforms in the BXD family reveals both common as well as unique genetic control of expression of Gabra1, Gabrb2, Gabrb3, Gabrg2, and Gabra2 ( Figure 2). Both Gabra1 and Gabra2 have highly conserved cis-regulation and high LRS values are always associated with either the B allele (Gabra1) or D allele (Gabra2), regardless of brain region or array platform. In contrast, allelic control of Gabrb2, Gabrb3, and Gabrg2 is more complex and varies across tissue and platform ( Figure 2). Whole brain RNA-seq reveals significant expression differences based on genotype for Gabra2 (fold change = 1.9, p = 0.05) and Gabra1 (fold change = 19.6, p = 0.0002). Higher expression is associated with the D allele for both subunits. Gabra6, a gene expressed almost exclusively in the cerebellum, is also cis modulated (higher expression associated with the D allele) (fold change = 4.5, p = 0.0004). We detected no significant difference in expression of Gabrb2, Gabrg2, or Gabrb3 in our RNA-seq analysis. Cis-regulation of Gabrg2 by local polymorphisms has been confirmed by Ciobanu and colleagues using traditional allele-specific expression assays [4]. For Gabrg2, higher expression is associated with the D allele (fold change = 1.5, p = 4 E-05). Allele-specific expression assays for Gabra1 failed to confirm a significant expression difference (fold change = 1.03, p = 0.09, higher D allele expression). Quantitative real time PCR validated the difference in expression for both Gabra1 Figure 1. Summary of significant genetic regulation in the hippocampus. Probe sets represent columns and ascending genomic location for each chromosome represent rows. Only chromosomes with a significant QTL (LRS.15) are shown. Several GABA(A)R subunits are regulated by both cis and trans eQTL (Gabra1, Gabra2, Gabrb2, and Gabrb3) whereas some are modified exclusively in trans (Gabra5, Gabrg1). Hue intensity provides an indication of the strength of the association between gene expression and genomic location. Red and blue alleles indicate that the D or B allele increases trait values, respectively. Arrowheads show the location of each cognate gene. PC1 = first principal component, Chr = Chromosome, Mb = Megabase, Prox = proximal, mid = middle, dist = distal. doi:10.1371/journal.pone.0034586.g001 Figure 2. Consistency of cis modulation across tissues. Genetic regulation of the 5 subunits regulated by strong cis eQTLs in the hippocampus was compared across three different platforms Affymetrix M430 (WB = whole brain, CB = cerebellum), Affymetrix MoGene 1.0 ST (AMY = amgydala and HYP = hypothalamus), and Illumina (CTX = neocortex, STR = striatum). Each probe set was required to have a mean expression of at least 8 in each tissue database. Blue and red indicate an association between the B6 and D2 allele at each locus and higher probe set expression, respectively. The more intense hue indicates a greater association and a corresponding increase in significance. For all data except the HYP and AMY, probes overlapping SNPs were filtered out and a principal component analysis was performed on remaining probes to capture a component (PC1) explaining the majority of the variation in probe set expression. Percent of variance explained by each PC1 ranged from (,45% to 95%). Because of the large number of probes included for measurement in the HYP and AMY, a probe level analysis was not possible. Some probe sets (indicated by *) contained 1 to 3 probes overlapping a variant. It is hard to predict the overall effect of these variants but, in general, those probe sets with SNPs that have higher expression associated with the B6 allele (blue) have a higher chance of being false due to SNP artifacts. Cis modulation of expression across tissues is most consistent for Gabra2 and Gabra1. doi:10.1371/journal.pone.0034586.g002 (higher expression in the B6 strain) and Gabra2 (higher expression in the D2 strain) [7].
We analyzed cis-regulation of expression at the level of mRNA region using probe sets designed to target distinct transcript regions-59 UTRs, coding exons, and 39 UTRs. Gabra2 has robust and consistent cis-modulation across all transcript regions and higher expression is always associated with the D allele ( Figure 3). In contrast, cis-modulation of Gabra1 and Gabrg2 expression tends to be associated with high B allele expression but this regulation is not congruent across all parts of the transcript. The genetic modulation of Gabrb2 and Gabrb3 expression is even more complex. Proximal 59 features of the Gabrb2 transcript show opposite allelic regulation compared to distal 39 features, and only one exon shows strong genetic cis modulation for Gabrb3. Coexpression between exon probe sets for individual subunits hints at alternative exon usage for Gabra1 and may provide evidence of multiple isoforms of Gabrb2 and Gabrb3 (Table S3) in the hippocampus.
As shown in Table 2, Gabra1, Gabrb2, and Gabrg2 are strongly cis-modulated in both the BXD family and one other genetic cross. Gabrb3 and Gabra2 are regulated by a strong cis eQTL in the BXD population and two other crosses.
Downstream consequences of genetic variation on higher order phenotypes GABA(A)R subunits that contain strong polymorphisms controlling their own expression may have multiple downstream effects on the expression of other genes or on higher order functional and Figure 3. Summary of genetic regulation of GABA(A) R subunits across transcript features. Genetic regulation of exon, intron, and 39 UTR features for each of the five subunits regulated by strong cis eQTLs was measured in the hippocampus using the Affymetrix exon platform (database accession number GN206). Each probe set was required to have a mean expression of at least 10 and all probe sets overlapping SNPs were excluded. Blue and red indicate an association between the B6 and D2 allele at each locus and higher probe set expression, respectively. The more intense hue indicates a greater association and a corresponding increase in significance. Only Gabra2 shows strong consistency of genetic regulation across mRNA regions. doi:10.1371/journal.pone.0034586.g003 behavioral traits. These effects should in turn map back to the location of the gene and the causal sequence variant. For example Li, Mulligan, and colleagues used this approach to extract a large number of phenotypes controlled by a B vs D polymorphism in catechol-O-methyltransferase (Comt) [8]. We find that anxietyassociated phenotypes are associated with variants near Gabrb3 (LRS.12, B allele) on Chr 7 (GeneNetwork BXD Trait IDs: 11389, 11390, 11385, 11388) [9]. Despite high endogenous variation in expression, no phenotypes map back to the locations of Gabra1, Gabrb2, and Gabrg2 on Chr 11 or Gabra2 on Chr 5.

Trans modulation of individual subunits
We identified five subunits with hippocampal expression variation associated with unique trans eQTLs that have LRS scores above 15 ( Figure 1 and Table S2). Within each of these five trans eQTL we nominated candidate genes that may control subunit expression based on five criteria: (1) correlation with GABA(A)R subunit expression in the hippocampus, (2) overlapping expression patterns in the hippocampal formation (Allen Brain Atlas), (3) true cis modulation of candidate regulatory genes in the hippocampus, (4) partial correlation analysis, and (5) significant correlation in normal aged human brain (database accession number GN314) [10]. These criteria were selected based on previous studies, principals of QTL mapping, and translational relevance (see methods).
Gabra1. Collectively, Gabra1 probe sets target the last coding exon and proximal 39 UTR (1421280_at), the proximal to middle portions of the 39 UTR (1421281_at), and the distal segment of the 39 UTR (and 1436889_at). Probe sets 1421280_at and 1421281_at are modulated by trans eQTLs on Chr 16 at 92.79 Mb (LRS = 18.4, higher expression associated with the B allele) and 97 Mb (LRS = 17.7, higher expression associated with the B allele), respectively. No candidates were identified near 97 Mb. The best candidate within the trans QTL region at 92.79 Mb is Son cell proliferation protein (Son). Two probe sets of Son (1420951_a_at and 1439074_a_at) are correlated with Gabra1 probe set 1421280_at (r.0.39) and retain the highest partial correlation. Son is polymorphic between B and D haplotypes and is highly expressed in hippocampus. In human brain, SON (201086_x_at) is highly correlated (r = 0.56) with expression of the middle and distal region of GABRA1 (244118_at) but not the last two exons and proximal 39 UTR. Chr = Chromosome, Mb = Megabase, N/S = Not significant or suggestive. Shading = Possible SNP artifact as the probe overlaps a SNP between parental strains (rs28224563), * = Probes overlapping SNPs removed and mapping performed with the first principal component of the remaining probes, ** = SNP between parental strains (rs28218735) overlapping probe but the allele effect (Cast high) indicates no effect on probe binding. NOTE: ILM4150164 is overlapping rs37397520 which is a SNP between B6 and PWK/PhJ but is probably not segregating in the LXS strains which are derived from A, AKR, BALB/c, C3H/2, C57BL, DBA/2, IS/Bi, and RIII. However the genotype of IS/Bi and RIII is unknown and IS/Bi is now extinct. For Gabra1, a probe from probe set rs28224563 overlaps a SNP between parental strains. doi:10.1371/journal.pone.0034586.t002 Gabra2. A region on Chr 10 centered on ,86 Mb contains genes and variants that influence Gabra2 expression in the hippocampus (LRS.15, higher expression associated with the D allele). Within this chromosomal region there are several candidates whose expression is strongly cis-regulated, correlated with Gabra2 expression, and that share the pattern of expression of Gabra2 in the hippocampus. Candidates include host cell factor C2 ( Expression of Nt5dc3, Syn3, and Hcfc2 is negatively correlated with Gabra2 (r = 20.52, 20.35, and 20.43, respectively) while Stab2 expression is positively correlated (r = 0.39). Nt5dc3 and Hcfc2 have moderate to high expression in hippocampal regions CA1 through CA3, as well as in the dentate gyrus. Syn3 is relatively abundant in dentate gyrus, whereas Stab2 is not abundantly expressed. The highest partial correlation after correction was with Nt5dc3. Of the three candidates, Hcfc2 (235264_at) has the highest correlation with GABRA2 (1554308_s_at) expression in human brain (r = 0.63).
Gabra5. Expression of Gabra5 is significantly regulated by sequence variants on Chr 19 around 16 Mb (LRS = 17, higher expression associated with the D allele). The strongest candidate is riboflavin kinase (Rfk, 17.4 Mb). Expression of Rfk (1415737_at) varies 2-fold across BXD strains, is positively correlated (r = 0.431) with Gabra5 probe set expression, and is well expressed in the hippocampus. Rfk is also highly polymorphic and expression is associated with a strong cis eQTL (LRS.50) with ,50% higher expression from the D allele. A ,3-fold variation between the B6 and D2 strains is also confirmed by RNA-seq. In humans, expression of RFK (203225_s_at) is strongly correlated (r = 0.56) with GABRA5 (206456_at) expression in several brain tissues from normal aged brain.
Gabrb2. Gabrb2 probe sets collectively target different mRNA features. Expression of probe sets targeting the middle to distal segment of the 39 UTR are modulated in trans by a region on Chr 8 at ,14.5 Mb. Modulation is significant for some probe sets (1429685_at, LRS = 16.3) and suggestive for others (1428205_ x_at and 1428204_at, LRS.10) and always associated with higher expression from the D allele. The best candidate gene is kelch repeat and BTB (POZ) domain containing 11 (Kbtbd11, 1430073_at, 15.0 Mb). Kbtbd11 is polymorphic, has high expression in the hippocampus, and is correlated with the expression of Gabrb2 In situ binding intensity in the hippocampus is low for Atp8a2, moderate for Tnfrsf19, and high for Ebpl. All three remain correlated with Gabrg1 after control for linkage disequilibrium, but Atp8a2 and Tnfrsf19 retain the highest partial correlation. All three are polymorphic between B and D haplotypes and there is evidence of differential isoform expression in RNA-seq data sets for both Atp8a2 and Tnfrsf19 (see ucscbrowser.genenetwork.org). Only EBPL (223306_at) is modestly correlated (r = 0.43) with GABRG1 (241805_at) expression in human brain, although the function of this gene is not known.

Subunit coexpression
Genetic variation results in profound and highly variable patterns of transcript expression among individuals. We can leverage this covariation to identify potential partnerships among subunits and among members of larger GABA(A)R-associated expression networks. We examined subunit coexpression signatures from BXD hippocampus (database accession number GN110) to detect possible subunit combinations and joint expression of functional GABA(A)R subtypes ( Figure 4 and Table S4). As expected, subunits known to form the major heterotrimeric receptor [11,12]-Gabra1, Gabrb2, and Gabrg2-are generally well and positively coexpressed-higher or lower expression of one isoform in a particular strain being associated with matched variation of functionally associated isoforms.

Discordant expression
A more intriguing and unexpected finding is that probe sets representing different regions of the same subunit (e.g, last coding exons, proximal 39 UTRs, and distal 39 UTRs) are often not correlated and have non-overlapping sets of covariates ( Figure 4 and Table S4). Gabra4, Gabrb2, Gabrb3, and Gabrg2 each have two or more probe sets targeting different parts of the same mRNA that are uncorrelated. These discordant probe sets have excellent sequence specificity, are all well expressed, and have biologically relevant expression correlates. The inconsistency can be resolved by comparing expression patterns between similar regions of different subunits. For example, probe sets of Gabra1, Gabrb2, and Gabrg2 that target 39 coding exons and the proximal 39 UTR are highly correlated, as are probe sets targeting the middle and distal 39 UTR. Within these two regions probe sets have an average correlation of 0.50 and 0.65, respectively. Many highly correlated (r.0.40) subunit pairs also follow this trend of coexpression within mRNA region: Gabra2 and Gabrb1 (last coding exons and proximal 39 UTR), Gabra4 and Gabrg1 (more distal 39 UTR), Gabra4 and Gabrg2 (more distal 39 UTR), Gabra5 and Gabrb3 (more distal 39 UTR), and Gabra4 and Gabra5 (more distal 39 UTR).
To evaluate the consistency and biological relevance of this finding, we examined coexpression between probe sets targeting different mRNA regions of the same subunit across brain tissues in both mouse and human. We were able to replicate the striking correlation between probe sets targeting different regions of the same mRNA in whole brain and cerebellum (database accession numbers GN123 and GN56, respectively) from the BXD family for Gabra1, Gabra3, Gabrb1, Gabrb2, Gabrg2, and Gabrb3 (Table S5). We were also able to replicate the complex correlation pattern between coding exons and different regions of the 39 UTR for Gabra5, Gabrb2, Gabrb3, and Gabrg3 in normal adult human brain tissue (database accession number GN314) ( Figure 2). This pattern of mRNA regional covariation across probe sets in human varies depending on cortical region.

RNA-sequencing
To explore GABA(A)R subunit coexpression networks on a more global level, we calculated correlations based on gene-level expression data from whole brain RNA-seq ( Figure S3). This analysis corroborates some features of the array data, including the strong coexpression of Gabra1, Gabrb2, and Gabrg2, as well as a strong positive correlation for most of the subunits. Subunits with negative correlations in the RNA-seq coexpression network relative to most other isoforms include Gabra2. Similar to the complex correlation pattern between different mRNA regions of GABA(A)R subunits on the M430 array, we also see strong negative correlations between each predicted isoform of Gabrb3 and Gabra6.

Coordinated genetic regulation of subunit expression and impact on phenotype
Are sets of different subunits and their coexpression networks in the hippocampus controlled jointly by common trans eQTLs? We identified three unlinked regions, all on Chr 5 with peaks at ,90-92 Mb, 117-120 Mb, and 137-138 Mb ( Figure S4) that modulate expression of multiple subunits. The most proximal region of Chr 5 controls expression of the last coding exon and the proximal and middle part of the 39 UTR of Gabrb3 and the distal 39 UTR of Gabra5. The middle QTL controls the expression of middle to distal parts of the 39 UTR of four subunits-Gabra1, Gabrb2, Gabrb3, and Gabrg2. The last QTL modulates the expression of the last coding exon and the proximal 39 UTR of six subunits-Gabra1, Gabra2, Gabra3, Gabrb2, Gabrb3, and Gabrg2. The last coding exons and proximal 39 UTR of four of these subunits (Gabra1, Gabrb2, Gabrb3, and Gabrg2) are jointly-albeit unequally-regulated by the middle and distal Chr 5 loci ( Figure S4). Higher expression is associated with the D allele at these three loci with the exception of the distal 39 UTR of Gabra5 and Gabrb3, and the middle 39 UTR of Gabrg1. We identified candidate genes in each interval that may regulate the network covariation based on criteria described previously: (1) correlation with GABA(A)R subunit expression in the hippocampus (accession number GN110), (2) overlapping expression patterns in the hippocampal formation (Allen Brain Atlas), (3) true cis modulation of candidate regulatory genes in the hippocampus, (4) partial correlation analysis.  PCA was used to extract PC1 for a single subunit when more than two probe sets were strongly correlated (r.0.8) and the average (AVG) was used when two probe sets were strongly correlated. prox = proximal, mid = middle. The highest positive correlations are between subunits targeting the same mRNA region. This is especially true for members of the most common heterotrimeric receptor-Gabra1, Gabrb2, and Gabrg2. doi:10.1371/journal.pone.0034586.g004 tion is between Naaa (N-acylethanolamine acid amidase) and both Gabra5 and Gabrb3. N The distal QTL, identified previously as Trans5a, modulates the expression of multiple subunits and also regulates the expression of hundreds of other transcripts in the hippocampus [13].
Overall and colleagues nominated Zkscan1 as a candidate trans regulatory gene in this region. We re-examined Trans5a with specific reference to GABA(A)R subunits. This region contains three genes associated with strong cis eQTLs that contain multiple sequence polymorphisms and are involved in transcriptional processes-Taf6, Zipro1, and Zkscan1. Both Taf6 and Zkscan1 are abundantly expressed in the hippocampal formation but Taf6 is not as well correlated with subunit expression. Both Zipro and Zkscan1 are correlated (|r|.0.5) with the expression of one or more GABA(A)R subunits but Zipro1 expression is most correlated with Gabrg1 expression, only. In contrast, Zkscan1 (1447944_at) expression is highly correlated (r.0.5) with that of multiple subunits (last exon and proximal 39 UTR probe sets of Gabrb3, Gabrg2, Gabra1, and Gabrb2). Zkscan1 probe set 1447944_at also retained the highest partial correlation with GABA(A)R subunits, indicating a residual biological signature that preserves the structure of the network.
Gene variants within the middle and distal regions (118-120.5 Mb, and 137-138 Mb) modulate the expression of multiple GABA(A)R subunits and have a large impact on the expression of synaptic genes in the hippocampus. Variation within this region would be expected to have downstream effects on biological function. Multiple phenotypes are modestly to highly associated with genetic variation on Chr 5 including, doublecortin staining in the hippocampus (a measure of adult neurogenesis), behavioral response to cocaine injection, locomotor response to paraoxon injection, climbing scores after methamphetamine injection, locomotor response after saline injection, and multiple anxiety measures (Table S6).

Conserved coexpression in mice and humans
A survey in GeneNetwork across human and mouse data sets generated using different tissues and array platforms (see methods) identified a novel subset of genes that covary (r.0.30) with six or more GABA(A)R subunits. These include Gucy1a3 and Gucy1b3 (Table S7) that interact functionally to form a soluble guanylate cyclase heterodimer involved in GDP/GTP balance and nitric oxide-mediated signal transduction [14]. Another example is platelet-activating factor acetylhydrolase 1B sububit alpha (Pa-fah1b1). Also known as Lis1, this gene is a key regulator of neuronal signal transduction and the cytoskeleton, specifically through the rho GTPase Cdc42 (Table S8) [15].

Expression variation and local modulation of subunit expression
Many GABA(A)R subunits, including, Gabra1, Gabra2, Gabra5, Gabrb2, Gabrb3, and Gabrg2, have remarkably high levels of variation in expression among both inbred strains of mice and normal human populations [1]. Sources of variation include genetic, environmental/epigenetic, and technical factors. We have controlled for several major technical sources of variation, including hybridization artifacts generated by SNPs and small insertions and deletions [4]. Environmental variation within colonies of laboratory mice is modest. After the reduction of many non-genetic sources of variance, there are still large differences that are clearly due to downstream effects of sequence variants-either cis effects in or near GABA(A)R genes, or trans effects mediated by distant variants and gene loci.
Compared to genetic modulation of expression by trans eQTLs, modulation by cis eQTLs is generally much stronger and more consistent across platforms, tissues, and mRNA regions. In general recombination that occurs within a gene is exceedingly rare and most genes and gene clusters are inherited as a block from either the B6 or D2 parent. Highly consistent expression would then be expected across mRNA region and tissues when a single mRNA product is produced from a single gene loci inherited in this fashion. For example, most regions of Gabra2 mRNA are consistently cis modulated with low expression associated with the B allele regardless of platform, tissue, or cross ( Figure 1, Figure 2, and Figure 3). Remarkably, the widely used C57BL/6J parent strain expresses by far the lowest level of Gabra2 mRNA among 99 inbred strains we have studied [13]. This locus is largely identical by descent among the parental strains and the only B vs. D sequence variants near Gabra2-eight small indels and nine SNPs-fall within introns or intergenic regions (ucscbrowser.genenetwork.org). One or more of these local variants contribute to a 2-to 3-fold difference in Gabra2 expression.
In contrast and contrary to our expectations, cis modulation of other subunits is not consistent across platforms, tissues, and mRNA/probe set target regions ( Figure 1, Figure 2, Figure 3, and Table S3) even though these genes and gene clusters are inherited as a single haplotype block. Determining the source of this discordance is inherently complex. For example, the B allele of Gabra1 has high expression in most array data sets, a finding that we cannot validate using either allele-specific expression assays or RNA-seq. In the case of Gabrg2, Gabrb2, and Gabrb3, alleles with higher expression are dependent on tissue, platform, and mRNA target region. A parsimonious explanation is the expression of a diverse repertoire of isoforms-a finding that is supported by the existence of alternative isoforms for many of the cis modulated subunits, including Gabrb2 and Gabrg2 [18,19]. Hundreds of B vs. D polymorphisms in Gabra1, Gabrb2, Gabrb3, and Gabrg2 may modulate isoform usage and expression among strains, highlight-ing the need for a complete molecular dissection to determine their impact.

Distant modulation of GABA(A)R subunits and genetic independence
Trans eQTLs are responsible for a significant level of genetic variation in subunits, even when expression is also under the control of a strong cis eQTL (Figure 1 and Table S2). For example, many of the cis-regulated subunits-Gabra1, Gabra2, Gabrb2, and Gabrb3-are co-controlled by one or more trans eQTLs. In addition, subunits without strong cis-modulation such as Gabrg1 and Gabra5, are strongly modulated by at least one trans eQTL.
Despite substantial coexpression among these subunits-especially for those comprising the major receptor subtype (Gabra1, Gabrb2, and Gabrg2)-individual subunits do not control the expression of other subunits either directly or indirectly (Figure 1).
Due to the complex evolutionary history of vertebrate genomes, most GABA(A)R subunit genes are organized into clusters. In mouse, the major clusters reside on Chrs 5, 7, 11, and X. Each contains a single c, one or two a, and a single b subunit. We expected evidence of receptor-to-receptor genetic co-modulation in the form of trans eQTLs with locations corresponding to the GABA(A) gene clusters. Remarkably, sequence variants and strong cis eQTLS within each cluster do not influence the expression of subunits in other clusters. While we do see strong cis eQTLs in some of these clusters there is no evidence of genetic ''cross-talk'' among clusters. All trans eQTLS were located in regions of the genome that do not contain any GABA(A)R genes ( Figure 2). We also rarely observe a high level of coexpression (r.0.6) between GABA(A)R subunits and genes known to be involved in receptor trafficking, such as gephyrin, GABARAP, Bdnf, neurexin, neuroligin, protein kinase A and C, Trak1, and others (Table S9). Sequence variants in these genes do not contribute to variation in subunit expression. Our findings demonstrate that there is a small set of novel and still undefined genes that modulate GABA(A)R expression.

Novel candidate QTLs and gene variants modulate subunit expression
Rfk, Son, and Kbtbd11 are the three best candidate genes modulating Gabra5, Gabra1, and Gabrb2 expression in the BXD hippocampus, respectively. Rfk is a TNF-receptor-1 and p22 (phox) binding protein that mediates NADPH oxidase activation and reactive oxygen species signaling [20]. Son is thought to be involved in mRNA processing [21,22,23]. Although, highly correlated with the expression of Gabrb2 probe sets in both mouse and humans, the function of Kbtbd11 is unknown. Gabrg1 has three good candidate genes-Atp8a2, Ebpl, and Tnfrsf19. Atp8a2 and Tnfrsf19 play a role in lipid transport across membranes and signaling pathways associated with myelination and axon regeneration, respectively [24,25]. Mutations in ATP8A2 have been associated with mental retardation [26]. The functions of Ebpl are unknown. Finally, Gabra2 is modulated by a trans eQTL on chromosome 10 at ,86 Mb that aligns precisely with a QTL that controls variation in reversal learning among BXD strains [27]. Laughlin and colleagues nominated Hcfc2, Syn3 and Nt5dc3 as strong candidates but Hcfc2, a gene with unknown functions, has the highest correlation with a2 expression in mouse and human.

Master QTLs that modulate expression of multiple subunits
Several loci are trans eQTLs for multiple subunits. In mouse hippocampus (database accession number GN110) three unlinked loci on Chr 5 modulate expression of nine subunits including Gabra1, Gabra2, Gabra5, Gabra3, Gabrb1, Gabrb2, Gabrb3, Gabrg1, and Gabrg2 ( Figure S4). The most proximal QTL modulates Gabra5 and the last coding exon and proximal to middle 39 UTR of Gabrb3. Naaa and Cox18, genes involved in endocannabinoid metabolism [28,29] and cytochrome c oxidase assembly [30] respectively, are excellent candidates. COX18 (227442_at) expression is correlated with expression of the middle 39 UTR of GABRA5 (206456_at, r = 0.41) and the distal 39 UTR of GABRB3 (227690_at, r = 0.42). Although COX18 is thought to be involved in assembly of mitochondrial respiratory chain complexes, it may also be involved in other undiscovered biological processes. In normal human brain (database accession number GN314), NAAA (215178_x_at) expression is highly correlated with both GABRA5 (215531_s_at, 20.62) and GABRB3 (227690_at, 20.62). NAAA, along with FAAH, is involved in the degradation of N-acylethanolaminesendogenous endocannabinoids that include the cannabinoid receptor type 1 (CB1) agonist anandamide [31]. Many GABAergic neurons and synapses express CB1 receptors [32,33]. Endocannabinoids can modulate GABAergic neurotransmission, through both CB1 receptor-dependent [34,35] and -independent mechanisms [36]. FAAH is generally considered to be the major enzyme responsible for endocannabinoid degradation-especially anandamide-in brain [37] but our data suggest a novel role for NAAA at GABAergic synapses.
The middle QTL on Chr 5 modulates expression of the middle to distal parts of the 39 UTRs of four subunits-Gabra1, Gabrb1, Gabrb2, and Gabrg2. Numerous complex traits map to this part of Chr 5, which has been identified previously as a QTL for adult hippocampal neurogenesis and morphology of the infrapyramidal mossy fiber tract [38]. Three strong candidates in this region are involved in cell signaling (Nos1 and Rasal1) [39] and transcriptional regulation (Thrap2 or Med13l) [40,41]. Krebs and colleagues also nominated Nos1 based on previous associations between Nos1 expression and neurogenesis and morphological variation. In human brain (database accession number GN314), multiple NOS1 probe sets targeting different mRNA regions (239132_at, 207309_at, 207310_s_at, and 1560974_s_at) exhibit negative and positive correlations (|r|.0.4) with the expression of multiple GABRA1 (244118_at), GABRB1 (207010_at), GABRB2 (1557122_s_at and 207352_s_at), and GABRG2 (1568612_at) probe sets targeting mainly the middle to distal 39 UTR. Much research has linked nitric oxide synthase activity with GABA levels and GABAergic neuronal subpopulations in hippocampus and other regions [42,43,44,45]. However, the same positive or negative correlation structure (|r|.0.46) is observed for MED13L probe sets (212207_at, 212209_at, 242911_at) with the exception of the GABRB1 probe set which is less well correlated (r = 20.35). In addition, MED13L is a RNA polymerase II transcriptional coactivator and could be involved in the regulation of subunit expression. In contrast, the expression of RASAL1 (219752_at) is strongly and negatively correlated (r.20.5) with GABRB2 (242344_at and 1557122_s_at) only.
he distal QTL on Chr 5 identified by Overall and colleagues is a major trans regulatory locus in the hippocampus (Trans5a) that controls the expression of multiple GABA(A)R subunits and hundreds of other transcripts. Both our study and that of Overall and colleagues highlights Zkscan1 as an excellent candidate gene. Multiple probe sets target different regions of Zkscan1 mRNA. Probe set 1447944_at measures expression from the proximal 39 UTR and is controlled by a strong cis eQTL. This probe set varies significantly among BXD strains (2.45-fold) and is highly correlated (r.0.6) with the expression of multiple GABA(A)R subunits including, Gabra1, Gabrb2, Gabrb3, and Gabrg2. Several insertions and deletions within Zkscan1 and just upstream of the 59 UTR could account for the enormous variation in expression among strains. In human brain (database accession number GN314), ZKSCAN1 (214670_at and 1557953_at) is highly correlated (r.|0.5|) with GABRA1 (244118_at), GABRA2 (207014_at and 1554308_s_at), GABRG1 (1552943_at), and GABRG2 (1568612_at) and modestly correlated (|r|.0.35) with GABRB2 (207352_s_at) and GABRB3 (227690_at and 229724_at). Almost nothing is known about the function of Zkscan1 but our data implicate this gene as a high priority target for functional validation.
In summary, we have identified strong candidate genes that modify the expression of one or more GABA(A)R subunits. We often observe differential genetic control corresponding to different mRNA regions of the same subunit. In general, correlations between candidate genes and GABA(A)R subunits are preserved across species. Identifying the causal gene in trans regulatory regions is likely to be more important for understanding GABA(A)R subunit regulation than identifying specific causal variants underlying a cis eQTL. This is because molecular networks will be much better conserved between species than specific sequence polymorphisms.

Prevalence of 39 UTR isoforms
The Affymetrix M430 array is unusual in that the designer, Dr. David Kulp, intentionally selected multiple probe sets to target both proximal and distal regions of the 39 UTR. This design is ideal for comparative analysis of UTR variants that are associated with major differences in mRNA trafficking [46] and higher order function [8]. CNS transcripts have unusually long 39 UTRs and tend to utilize alternative and more distal polyadenylation sites [47,48]. For example, the 39 UTR of Gabrb2 is nearly 6 Kb long-4 times longer than the protein coding sequence. We also see multiple consensus sites for polyadenylation, and evidence of longer and shorter 39 UTR isoforms for Gabra1, Gabrb2, and Gabrg2 in our whole brain RNA-seq data ( Figure S5, S6, S7).
An intriguing trend emerges from our dissection of GABA(A)R subunit expression in which probe sets targeting similar mRNA regions of different subunits are well correlated (Figure 4 and Table S4). This is especially true for the major receptor subtype-Gabra1, Gabrb2, and Gabrg2. Poor or negative correlations between probe sets representing different mRNA regions of the same subunit and high positive correlations between probe sets from different subunits that target similar regions are replicated across brain tissues from BXD strains as well as across human brain samples (database accession number GN314) (Table S5 and Figure  S2). We asked if these coexpression patterns were part of a more general transcriptional network in the hippocampus of BXD strains (database accession number GN110). Probe sets targeting more distal 39 UTR regions generally have higher correlations with other probe sets targeting distal regions whereas probe sets targeting the last few coding exons and proximal 39 UTR generally have higher correlations with other probe sets targeting these proximal regions (Table S9). This pattern of coexpression-short to short, and long to long-may be related to differences in mRNA processing based on 39 UTR length. Variation at the 39 end of mRNA can expose alternate miRNA and protein binding sites that change localization, stability, or translation efficiency (reviewed in [49]). Another intriguing possibility that could explain the discordant expression between 39 UTR and coding regions is the expression of distinct RNA transcripts from within the 39 UTR [50]. Mercer and colleagues mined RIKEN CAGE data to identify genes producing full-length transcripts from the 39 UTR, including Gabra1, Gabra6, Gabra3, and Gabrb2 (mouse) and GABRA2, GABRB1, and GABRB2 (human). Within this gene category and in 54% of the cases, expression of the 39 UTR and coding exon is discordant. Further molecular analysis of 39 UTR expression and new highthroughput methods for quantifying the 39 end of mRNA will be needed to resolve which mechanisms contribute to the divergent expression patterns we observe among GABA(A)R subunit mRNA regions. How these isoforms contribute to cell and synapse type, circuitry, and disease will be a critical aspect of future studies.

Functional impact of subunit variation and relevance to human disease
There is a lack of association between cis modulated GABA(A)R subunits and classic phenotypes generally linked with GABAergic signaling, including alcohol-related, anxiety, and pharmacological measures. Of course not every trait associated with these receptors has been measured among BXD strains and several physiological phenotypes, including receptor density and neuronal excitability are notable absent. However, the lack of association with available traits is still surprising given that as little as a 6 to 35% reduction in receptor density can produce robust phenotypic differences [51,52]. The sole exception is for b3, a subunit that is genetically linked to mouse behaviors associated with anxiety. Deletion of Gabrb3 leads to an increase in stretch-attend posturing-an phenotype associated with risk assessment [53]. Deletions of the 15q11-13 region in humans-including GABRB3-are associated with Angelman syndrome and severe cognitive deficits, a constellation of phenotypes matched in Gabrb3-deficient mice [54,55].
The variation level of many GABA(A)R subunits is very high among individual strains and yet the functional impact on classic phenotypes-behavior and neuropharmacological traits, in particular-is minimal. In a normal human population we observe even greater differences in expression of these subunits-up to 100-fold in adults. Individual humans are highly polymorphic and a small portion of this variation could result from platform-dependent artifacts caused by hidden variants. However, whole genome and transcriptome sequencing technologies are rapidly advancing and our findings suggest that a high level of expression variation will persist even after these types of technical errors are corrected. This is an important negative finding that demonstrates the impressive ability of higher order network properties to buffer the variation of their constitutes. Much of this buffering may occur during translational processing and a 10-fold difference in transcript may be muted at the protein level. How molecular networks buffer, exploit, or accommodate high levels of variation in transcript expression is a major question that has not yet been addressed.
Why study variation in expression that does not cause prominent differences in phenotype? It is rare to find a natural variant that produces a null or fully penetrant phenotype. The subtle push and pull of modest gene variants on expression is a natural perturbation that exposes the large-scale structure of networks of genes and molecules. Understanding network structure and relations is highly relevant for understanding complex traits and disease risk. We identified members of the GABA(A)R network that are conserved across rodents and humans ( Figure 5). These members are involved in signal transduction (Rgs4, Mapk9, Socs5) and synaptic transmission (Syt1) and are likely to be important for the function of GABAergic synapses. In the hippocampus, this network is under genetic control and expression is strongly modulated by the Trans5a QTL.
We also examined the set of conserved expression correlates of five or more GABA(A)R subunits and identified genes with the highest correlation levels in both human and mouse brain. These genes are likely to be important for subunit expression at GABAergic synapses and include, Lis1, Gucy1a3 and Gucy1b3 (Table S7, Table S8). Lis1 is strongly associated with expression of many subunits in both species and is a member of the conserved GABA(A)R network. During development Lis1-along with GABA and other signaling molecules-is a key mediator of neuronal migration [56,57]. Knockout mice have cortical and hippocampal defects and heterozygotes exhibit alterations in CA1 inhibitory circuits implicating dysfunction of the GABAergic system [58,59]. GUCY1A3 (a1) and GUCY1B3 (b1) are subunits of soluble guanylyl cyclase, the main receptor for the retrograde signaling molecule nitric oxide (NO). NO can enhance GABAergic signaling through activation of presynaptic guanylyl cyclase and production of cGMP [60,61,62,63].

Conclusion
Our comprehensive analysis highlights the way in which expression variation and covariation can be used to identify strong candidates and novel biological processes that play a critical role in GABA(A)R mediated inhibitory transmission. We now have ample evidence of profound alternative 39 UTR expression for many subunits in human and mouse brain. Elucidating the molecular regulation and functional consequences of these expression patterns is of the utmost importance. We also have identified a large and surprising list of regulatory genes and network members that are highly correlated with subunit expression in human and mouse brain. Most of these candidates are novel and suggest new and interesting modes of organization that influence GABA(A)R expression in brain.

Ethics Statement
All animal work was conducted according to an approved animal use protocol (UTHSC680) and in accordance with procedures approved by the Executive Committee of the Institutional Animal Care and Use Committee at The University of Tennessee Health Science Center.
Briefly, human brain tissue data sets include, ''GSE5281 Human Brain Normal Full Liang (Jul09) RMA'' (accession number GN314), ''GSE15222 Human Brain Normal Meyers (Apr09) RankInv'' (accession number GN234), and Harvard Brain Tissue Resource Center ''HBTRC-MLC Human Cerebellum Agilent Normal (Jun11) mlratio'' (accession number GN361). For data set GN314, cases include 16 normal aged subjects [10,66], and gene expression was profiled for six regions on the Affymetrix platform. For data set GN234, expression was profiled for temporal cortex and cortical tissue from 187 normal aged adults on the Illumina Sentrix Bead array (HumanRef-8) using Illumina's rank invariant transform [67]. Data set GN361 includes expression profiles for cerebellum, visual cortex, and prefrontal cortex from 170 normal subjects using a custom-made Agilent 44 K microarray. Only expression in the cerebellum was used for comparison of rodent and human coexpression. This data set was contributed by Merck Pharmaceutical through the Sage Bionetworks Repository.

Removal of probe sets overlapping B6 and D2 sequence variants
All 19 known GABA(A) subunit genes are represented by at least one probe set on the Affymetrix M430 2.0 and Exon 1.0 ST platforms. The former is designed to produce consensus estimates of expression across multiple mRNA isoforms using probes that target the 39 UTR and the last few coding exons. In contrast, the exon array is designed to target every exon. Illumina probes generally target a single region from the last coding exons or the 39 UTR. Probe target sequences corresponding to GABA(A) probe sets for Illumina and exon hippocampus BXD data sets were matched to all B6 and D2 SNPs and small indels. The offending probe set was removed every time an overlap occurred. For the M430 data we filtered out probes overlapping polymorphisms using tools in GeneNetwork. We used the remaining probe values and principal component analysis to construct corrected probe set values from the first principal component (PC1). B6 vs. D2 sequence variants were generated at UTHSC by comparison of the B6 reference genome and the D2 genome (1006 whole genome coverage). This level of coverage was achieved using data generated from both the SOLiD and Illumina sequencing platforms and exceeds that of any other publicly available variant database. Variant data can be visualized by gene at ucscbrowser.genenetwork.org or by probe set using the RNA-seq tool in GeneNetwork. A list of all affected and unaffected probe sets is available in Table S1.

Heritability
We calculated the broad-sense heritability to estimate the influence of genetic factors on variation in GABA(A)R subunit expression for the Hippocampus M430 data set. Heritability is the genetic variance (variance among strain means) divided by the total variance (variance among all measurements).

RNAseq Analysis
Total RNA was isolated from the whole brain of B6, D2, 31 BXD, and two reciprocal F1 hybrid strains using the QIAGEN RNeasy kit. RNA quality was assessed using the Agilent Bioanalyzer 2100. Ribosomal RNA (rRNA) was removed from total RNA using the RiboMinus Eukaryote kit (Invitrogen). The SOLiD Whole Transcriptome Analysis kit was used to generate cDNA libraries for each strain. Each library was subsequently sequenced on the ABI SOLiD platform at the University of Tennessee Health Science Center. Data was analyzed using Applied Biosystems whole transcriptome software tools (http://www.solidsoftwaretools.com/). Reads are mapped to the mouse reference genome (mm9, National Center for Biotechnology Information (NCBI) build 37) with a minimum score of 24. Technical replicates were merged for each strain to produce a single BAM file (data available in our mirror of Galaxy at GeneNetwork, galaxy.genenetwork.org/library). The Cufflinks tool [68] was used within the Galaxy framework [69,70] to generate gene level expression data based on RefSeq gene models (build 37) using a maximum intron length of 300,000, minimum isoform fraction of 0.1, and pre-mRNA fraction of 0.15.

Identification of Candidate Genes in Upstream Regulatory Regions
Genes physically located within trans regulatory regions were filtered based on cis-modulation, coexpression within the tissue of interest, and covariation with the target GABA(A) subunit.
Expression in the hippocampus was confirmed for each candidate using the Allen Brain Atlas resource. Any candidate gene probe sets overlapping SNPs that also showed significantly higher expression associated with the B Allele were removed (see Ciobanu and colleagues 2010 for details). Candidates not polymorphic between B6 and D2 were also removed. Overlapping expression of the control and target gene in the tissue of interest and the presence of polymorphisms in the control gene locus are required conditions of nearly all cis-trans interactions as both genes can act within the same biological process or pathway [71]. If a candidate gene did not demonstrate highly significant cis-regulation (LRS.15) or covariation (r.|0.3|) with the target gene it was excluded from further analysis. Correlational analysis has also proven to be a useful tool to narrow the list of candidate genes [4,72].
In general, multiple candidates were identified within each trans regulatory region using the above approach. Linkage between genes in close proximity can make it difficult to further narrow the list of candidate genes. To facilitate candidate gene identification we used the partial correlation feature available in GeneNetwork to determine the most likely set of cis-modulated genes within upstream regulatory regions that control target GABA(A) subunit expression. A partial correlation is the degree of association between a primary variable and a target variable after removing the effect of one or more control variables [73]. Specific details on partial correlations have been described in detail elsewhere [74]. In this case, the correlation between GABA(A)R subunit probe set expression (primary variable) and the expression of candidate regulatory gene probe sets (target variables) was measured after controlling for genetic variation at the peak trans QTL marker. We reasoned that, in the absence of genetic variation, residual biological variation would retain the correlation between the expression of the GABA(A)R subunit and the candidate regulatory gene. Therefore the candidate regulatory gene with the lowest reduction in the correlation between itself and GABA(A)R subunit expression is a superior candidate. Finally, covariation of the candidate control gene and the target subunit were examined in normal aged human brain samples (database accession number GN314) [10]. Correlations between control and target genes that are preserved across species are more likely to represent real biological networks of high translational relevance and were used to further prioritize candidate genes.

Exon Coexpression
The ''UMUTAffy Hippocampus Exon (Feb09) RMA'' (accession number GN206) data set was used to investigate coexpression between expressed exons and introns for individual GABA(A)R subunits. As described previously, probe sets overlapping SNPs were removed and only those probe sets with an average log2 expression of 10 or greater were included in the analysis. Pair-wise Pearson and Spearman correlation coefficients were calculated among probe sets for each subunit using tools available in GeneNetwork. The correlation structure between probe sets was used as an indirect method to identify possible instances of alternative splicing-observed as a lack of correlation between probe sets targeting different exons of the same subunit. Probe sets from this array do not target exon-exon junctions and do not provide a direct measure of splicing.

Cross-Species Coexpression
Modules (usually identified by a color) containing GABA(A) R subunits were identified from a weighted gene coexpression analysis of human and mouse data sets [16]. Most GABA(A) subunits had the highest correlation with the black and red module eigengenes in human and the tan module eigengene in mouse. Correlation with these three module eigengenes (R.0.3) was used to define coexpressed genes, including the GABA(A) subunits (GABRA1, GABRB2, GABRG2, and GABRB3). Overlap between the human and mouse module membership was determined by gene symbol. The expression of this consensus module of 28 genes was explored using the Hippocampus M430 database. A single probe set was selected for each gene based on correlation within the module. The first principal component captured ,55% of the expression variation of module members and was used for QTL mapping. Genes from the human and mouse modules (r.0.3) were also compared to the results of a limited GeneNetwork cross-species survey to assess the amount of connectivity with other GABA(A) R subunits.
The criterion for consistent cross-species association was defined as: r$|0.3| for three human (GN314, GN234, GN361) and three mouse brain data sets (GN110, GN133, GN268). Because of the number of samples within each database, this correlation is often highly significant (p,0.001). Figure S1 Pairwise comparison of platform expression for GABA(A)R subunits. Scatter plots were constructed for the expression data from each platform shown in Table 1. All probe sets were used for this analysis. Pairwise correlations based on Pearson's r are shown at the top left-hand corner of each scatter plot. M430 (accession number GN110), exon (accession number GN206), and in situ data (Allen Brain Atlas Resource) are from the hippocampus while RNA-seq data is from whole brain. There is good agreement for the expression of most subunits as assayed across platforms. More modest correlations (,0.5) are generally associated with comparisons between different tissue types-whole brain and hippocampus.  Figure S5 Summary of Gabra1 39 UTR expression as assayed using RNA-seq. Data available at ucscbrowser.genenetwork.org. Gabra1 is located on the (2) strand. Tapering of reads from the proximal to more distal 39 UTR indicates expression of longer and shorter 39 UTR isoforms. This is especially prevalent in the striatum, which included only polyadenylated transcripts for RNA-seq. Track 1 shows the chromosomal location at top followed by gene model based on RefSeq as track 2. Larger width indicates coding exons and smaller width indicates 39 UTR. All known variants between the B6 and D2 genome are indicated as part of track 3. Tracks 4 through 7 summarize the normalized (+) and (2) strand reads for the B6 (blue) and D2 (red) parental strains in whole brain (RiboMinus method). Tracks 8 and 9 show normalized reads in the striatum of 10 B6 and 11 D2 mice (PolyA enrichment method, [75]). For tracks 4 through 9 the scale to the left shows the read number, which will be lower for whole brain since it represents only a single RNA-seq run from one animal. Track 10 shows previously characterized mRNA species in mouse. Track 11 shows the degree of sequence conservation in mammals and Track 12 shows the presence of any repetitive DNA which is masked by most RNA-seq alignment algorithms. (TIF) Figure S6 Summary of Gabrb2 39 UTR expression as assayed using RNA-seq. Data available at ucscbrowser.genenet-work.org. Gabrb2 is located on the (+) strand. Unequal distribution of reads in the 39 UTR, especially evident in the striatum, could indicate alternative polyadenylation or splicing. Track 1 shows the chromosomal location at top followed by gene model based on RefSeq as track 2. Larger width indicates coding exons and smaller width indicates 39 UTR. All known variants between the B6 and D2 genome are indicated as part of track 3. Tracks 4 through 7 summarize the normalized (+) and (2) strand reads for the B6 (blue) and D2 (red) parental strains in whole brain (RiboMinus method). Tracks 8 and 9 show normalized reads in the striatum of 10 B6 and 11 D2 mice (PolyA enrichment method, [75]). For tracks 4 through 9 the scale to the left shows the read number, which will be lower for whole brain since it represents only a single RNA-seq run from one animal. Track 10 shows previously characterized mRNA species in mouse. Track 11 shows the degree of sequence conservation in mammals and Track 12 shows the presence of any repetitive DNA which is masked by most RNA-seq alignment algorithms. (TIF) Figure S7 Summary of Gabrg2 39 UTR expression as assayed using RNA-seq. Data available at ucscbrowser.genenetwork.org. Gabrg2 is located on the (2) strand. Tapering of reads from the proximal to more distal 39 UTR indicates expression of longer and shorter 39 UTR isoforms. This is especially prevalent in the striatum, which included only polyadenylated transcripts for RNA-seq. Track 1 shows the chromosomal location at top followed by gene model based on RefSeq as track 2. Larger width indicates coding exons and smaller width indicates 39 UTR. All known variants between the B6 and D2 genome are indicated as part of track 3. Tracks 4 through 7 summarize the normalized (+) and (2) strand reads for the B6 (blue) and D2 (red) parental strains in whole brain (RiboMinus method). Tracks 8 and 9 show normalized reads in the striatum of 10 B6 and 11 D2 mice (PolyA enrichment method, [75]). For tracks 4 through 9 the scale to the left shows the read number, which will be lower for whole brain since it represents only a single RNA-seq run from one animal. Track 10 shows previously characterized mRNA species in mouse. Track 11 shows the degree of sequence conservation in mammals and Track 12 shows the presence of any repetitive DNA, which is masked by most RNA-seq alignment algorithms. (TIF)