Transcriptional maturation of the mouse auditory forebrain

Abstract

Background

The maturation of the brain involves the coordinated expression of thousands of genes, proteins and regulatory elements over time. In sensory pathways, gene expression profiles are modified by age and sensory experience in a manner that differs between brain regions and cell types. In the auditory system of altricial animals, neuronal activity increases markedly after the opening of the ear canals, initiating events that culminate in the maturation of auditory circuitry in the brain. This window provides a unique opportunity to study how gene expression patterns are modified by the onset of sensory experience through maturity. As a tool for capturing these features, next-generation sequencing of total RNA (RNAseq) has tremendous utility, because the entire transcriptome can be screened to index expression of any gene. To date, whole transcriptome profiles have not been generated for any central auditory structure in any species at any age. In the present study, RNAseq was used to profile two regions of the mouse auditory forebrain (A1, primary auditory cortex; MG, medial geniculate) at key stages of postnatal development (P7, P14, P21, adult) before and after the onset of hearing (~P12). Hierarchical clustering, differential expression, and functional geneset enrichment analyses (GSEA) were used to profile the expression patterns of all genes. Selected genesets related to neurotransmission, developmental plasticity, critical periods and brain structure were highlighted. An accessible repository of the entire dataset was also constructed that permits extraction and screening of all data from the global through single-gene levels. To our knowledge, this is the first whole transcriptome sequencing study of the forebrain of any mammalian sensory system. Although the data are most relevant for the auditory system, they are generally applicable to forebrain structures in the visual and somatosensory systems, as well.

Results

The main findings were: (1) Global gene expression patterns were tightly clustered by postnatal age and brain region; (2) comparing A1 and MG, the total numbers of differentially expressed genes were comparable from P7 to P21, then dropped to nearly half by adulthood; (3) comparing successive age groups, the greatest numbers of differentially expressed genes were found between P7 and P14 in both regions, followed by a steady decline in numbers with age; (4) maturational trajectories in expression levels varied at the single gene level (increasing, decreasing, static, other); (5) between regions, the profiles of single genes were often asymmetric; (6) GSEA revealed that genesets related to neural activity and plasticity were typically upregulated from P7 to adult, while those related to structure tended to be downregulated; (7) GSEA and pathways analysis of selected functional networks were not predictive of expression patterns in the auditory forebrain for all genes, reflecting regional specificity at the single gene level.

Conclusions

Gene expression in the auditory forebrain during postnatal development is in constant flux and becomes increasingly stable with age. Maturational changes are evident at the global through single gene levels. Transcriptome profiles in A1 and MG are distinct at all ages, and differ from other brain regions. The database generated by this study provides a rich foundation for the identification of novel developmental biomarkers, functional gene pathways, and targeted studies of postnatal maturation in the auditory forebrain.

Background

The development and maturation of the brain is an exceptionally complex biological process that depends on the coordinated expression of many thousands of genes and proteins [1]. In every region of the brain, much remains to be learned about the spatial and temporal properties of their expression patterns, regulation, and functional roles.

An important goal in sensory systems research is to understand the mechanisms that govern maturation, and how these factors affect or are affected by major milestones, such as the emergence of intrinsically generated electrical signals or the onset of activity evoked by extrinsic stimuli from the sensory environment. In the central auditory pathways of rodents, structural and functional development begins during gestation and continues through the first three to four postnatal weeks. Auditory processing capabilities develop rapidly between postnatal days P10–P16, catalyzed by the opening of the ear canals (~P12) and the associated shift from intrinsically to extrinsically generated patterns of electrical activity [2–5]. This window has provided researchers with a unique opportunity to document structural and functional maturation associated with the onset of hearing [6–9], and the formation of critical periods for the plasticity of sound feature encoding and behavior [10–14]. In this context, plasticity refers to the potential for structural and functional change at the level of the synapse or networks of neurons. These changes are mediated by intrinsic mechanisms at the cellular and molecular levels, and shaped by extrinsic factors, such as the onset of sensory experience. Critical periods are windows of time during which the conditions for plasticity are such that the functional properties of a synapse or network can be altered by experience (or lack, thereof) in a manner that has long-lasting or permanent effects [15].

Efforts to characterize the cellular and molecular landscape during maturation, and their relationship to the specific mechanisms that regulate plasticity and critical periods are ongoing. A wide range of factors has been explored. Synaptic inhibition (GABA) and excitation (glutamate) are considered to be central regulators in the maturation of auditory response properties, and continue to be intensively studied [16–20]. Other studies have focused on the influences of neuromodulatory inputs (e.g., cholinergic, dopaminergic, serotonergic) [21–26] and the roles played by ion channels [27]. Myriad structural factors also impact neuronal activity, such as dendritic spine formation [28], gap junctions [29], synaptic morphology [30, 31], and myelin signaling and extracellular matrix formation [32, 14].

These studies, and many more, have contributed much to our understanding of the mechanisms involved. Yet, much remains to be learned, and it may be that important, even essential, mechanisms have not yet been identified. Lacking so far is application of a comprehensive broad-spectrum approach to identify novel mechanisms on a large scale. Among the techniques that could be employed, whole transcriptome sequencing of total RNA is a powerful tool for the generation of gene expression profiles and identification of functional biomarkers. By sequencing samples from different brain regions at several time-points during development, a series of snapshots documenting the influences of age and experience on the entire transcriptome can be acquired. This permits identification of significant changes in the expression of any coding or non-coding gene. So far, targeted profiling of up to about 2000 genes or proteins has been successfully used to identify changes in auditory brainstem nuclei associated with postnatal development [33, 34], hearing loss [35], and auditory cortex lesions [36]. To date, however, whole transcriptome profiles have not been generated for any central auditory structure in any species at any age.

To enhance the foundation for discovery along these lines, we used high-throughput next-generation sequencing of total RNA (RNAseq) to profile RNA expression in the primary auditory cortex (A1) and medial geniculate body (MG) of mice at selected time-points during postnatal development, before and after the onset of hearing (postnatal days P7, P14, P21, and adult). Differential expression analyses were employed to compare the transcriptomes between brain regions and age groups. Functional gene set analyses were performed to create reference libraries of gene families and functional gene ontology categories that have importance for brain function and developmental neurobiology. Several representative genesets were profiled in detail at the single gene level, one of which was explored by pathways analysis. Finally, all of the data were organized into an accessible and searchable database that facilitates the identification of genes that are involved in the maturation of the auditory forebrain.

To our knowledge, this is the first whole transcriptome sequencing study of maturation in the forebrain of any mammalian sensory system. Although the data are most relevant for the auditory system, they are generally applicable to forebrain structures in the visual and somatosensory systems, as well.

Methods

Tissue acquisition

All procedures were approved by the Animal Care and Use Committee at Massachusetts Eye and Ear Infirmary and followed the guidelines established by the National Institutes of Health for the care and use of laboratory animals. The morning that a new litter of pups was first observed was designated as P0. Brains were collected from 24 adult (8–10 weeks) and juvenile (P7, P14, and P21) male and female C57BL/6 J mice (Jackson Labs 000664) (N = 6 per age, equal numbers of males and females, total = 24). Animals were euthanized with a lethal dose of ketamine and xylazine (200/50 mg/kg, respectively) intraperitoneally. Brains were removed immediately, flash frozen on dry ice, and stored at−80⁰ C.

Sample acquisition

Frozen brains from 6 animals in each age group (3 male, 3 female) were sectioned at 40 μm in the coronal plane (rostral to caudal) on a sliding microtome and viewed through a surgical microscope. As established anatomical landmarks [37] became visible in the frozen tissue block, the regions targeted for sampling (A1, primary auditory cortex; MG, medial geniculate body), were extracted using a sterile tissue punch or curette of a size appropriate to the brain region (Additional file 1: Figure S1) (note that Additional Figs and Tables are indicated by inclusion of the letter S before the number). Punches from homologous areas of both hemispheres were combined in sterile tubes containing 400 μl of Trizol, homogenized for 45 s using a mechanized sterile pestle, flash frozen on dry ice, then stored at−80⁰ C. Each of the A1 and MG sample pairs were from the same animals, as indicated by the Sample ID and Animal ID codes in Additional file 2: Table S1.

A1 samples were obtained using a 0.5 mm diameter punch, with the ventral edge beginning approximately 1 mm dorsal to the rhinal fissure. Samples were centered on A1, but potentially also included some tissue in the adjacent auditory field dorsal to A1. MG samples were harvested with a curette after using a micro-dissecting scalpel to circumscribe its perimeter (Additional file 1: Figure S1b). For the MG, the microdissection procedure was intended to exclude the lateral geniculate nucleus (LGN), which was achieved by identification of the septum between the MG and LGN dorsolaterally, in the rostral third of the MG. Because there are no remnants of the LGN caudal to the MG at this point, the LGN was easily excluded by the dissection. Additional evidence that the LGN was successfully excluded is supported by comparison of our results with a prior study comparing the expression of a subset of genes in LGN and MG [38]. In that study, 10 genes had moderately-high to high levels of expression in the LGN (Zic4, Zic5, Ecel1, Isl1, Npy, Arx, Pvalb, Pmch, Pax6, Zfp503). All of these genes had very low to nominal expression at all age in the MG of our samples. Therefore, we conclude that there was no significant contamination by LGN. The dissection was also intended to exclude adjoining nuclei located medial, and ventral to the MG, but some tissue from these nuclei may have been included (e.g., suprageniculate, peripeduncular). The extreme rostral and caudal poles of the MG were largely excluded from these samples.

RNA extraction and sequencing

For each Trizol lysate, 100 μl of Reagent Grade Chloroform (Fisher Scientific, S25248) was added. The samples were centrifuged for 3 min on a desktop centrifuge to fractionate the aqueous and organic layers. After centrifugation, the resulting aqueous layer was carefully removed and transferred to 2.0 ml Sarstedt tubes (Sarstedt, 72.694) which were run on the QIAsymphony SP (Qiagen Corporation, Germany) using the QIAsymphony RNA Kit (Qiagen, 931636) and protocol RNA_CT_400_V7 which incorporates DNAse treatment. Prior to each run, the desk was uv-irradiated using the programmed cycle. The resulting RNA was eluted to 100 μl of RNase free water and stored at−80 °C in 2.0 ml Sarstedt tubes until use. Samples were initially quantitated using a Qubit fluorometric RNA assay (Life Technologies, Grand Island, NY). Additional analyses of purity and the quantitation of total RNA were performed using a NanoDrop spectrophotometer (Thermo Scientific) and Agilent RNA 6000 Pico chip (Agilent) according to the manufacturer’s protocol using the reagents, chips, and ladder provided in the kit. Quality control data for the 48 sequenced samples are contained in Additional file 2: Table S1.

RNAseq was performed by the Vanderbilt Technologies for Advanced Genomics core (VANTAGE). First, ribosomal reduction was performed on 1 μg total RNA using the Ribo-Zero Magnetic Gold Kit (Human/Mouse/Rat) (Epicentre), following the manufacturer’s protocol. After ribosomal RNA (rRNA) depletion, samples were purified using the Agencourt RNAClean XP Kit (Beckman Coulter) according to the Epicentre protocol specifications. After purification, samples were eluted in 11 μl RNase-free water. Next, 1ul ribosomal depleted samples were run on the Agilent RNA 6000 Pico Chip to confirm rRNA removal. After confirmation of rRNA removal, 8.5 μl rRNA-depleted sample was input into the Illumina TruSeq Stranded RNA Sample Preparation kit (Illumina) for library preparation. Libraries were multiplexed six per lane and sequenced on the HiSeq 2500 to obtain at least 30 million paired end (2x50 bp) reads per sample.

RNAseq data processing

The RNAseq data went through multiple stages of thorough quality control as recommended by Guo et al. [39]. Raw data and alignment quality control were performed using QC3 [40], and gene quantification quality control was conducted using MultiRankSeq [41]. Raw data were aligned with TopHat2 [42] against mouse mm10 reference genome, and read counts per gene were obtained using HTSeq [43]. Default settings were used for MultiRankSeq, TopHat2, and HTSeq. Normalized counts (used in all plots) were obtained by normalizing each gene’s count against the sample’s total read count, then multiplying by a constant (1 X 10⁶). Hierarchical clustering analysis and heatmaps were produced using the Heatmap3 [44] package from R (Fig. 1). For all samples, quality control data are contained in Additional file 2: Tables S2 – S3. The raw counts are contained in Additional file 2: Table S4. Differential expression analyses between all postnatal ages and brain regions were performed using MultiRankSeq [41], which combines three independent methods for RNAseq analysis: DESeq [45]; EdgeR [46]; BaySeq [47]. These three methods were chosen based on results of several previous studies in which multiple RNAseq differential analysis methods were compared for accuracy and sensitivity of read count-based data [48–52]. In analyses of the same dataset, the methods typically differ in numbers of differentially expressed genes identified in a comparison of any two samples, and also in direction of expression (up- or down-regulation). The false discovery rate (FDR < 0.05) was used to correct for multiple testing, and a given comparison was considered to be significant if all three methods identified it as significant. The differential expression data associated with each pairwise comparison (4 ages X 2 brain areas) are summarized in the Results section, with complete data for all genes for all comparisons contained in Additional files 3, 4, 5, 6, 7: Tables S5 – S20. These Additional files are Excel workbooks, organized by tabs corresponding to each supplementary Table. Within each of these files, the listing of single genes is ordered from the smallest to highest numerical ranking (i.e., highest to lowest degree of differential expression), based on p-values from DESeq, EdgeR, and BaySeq. The order can be changed with sorting and filtering functions in Excel.

Fig. 1

Grand summary of global gene expression in MG and A1 from P7 to adult. (Top) Unsupervised hierarchical clustering of samples by sex, brain region, and age. (Bottom) Heatmap summarizing total gene expression for each sample, arranged in columns by cluster. Each bar represents one gene. Color code denotes expression level

Validation of sequencing

Validation of sequencing was accomplished by in situ hybridization (ISH) of 4 genes (Gapdh, Slc32a1, Slc17a6, Slc17a7) in A1 and MG, also profiled in a related study of their maturational trajectories, regional patterns of expression, and co-expression within single neurons. Full methodological descriptions of the tissue processing, primer sequences, in situ hybridization, and quantification are available in Hackett et al. (2015, in press) [53]. Briefly, 3 animals in each age group from the same breeding colony were euthanized and perfused transcardially with 4 % phosphate buffered paraformalin. The extracted brains were sectioned at 50 μm in the coronal plane. Single colorimetric in situ hybridization was performed on sequential tissue sections processed for each gene (Additional file 8: Figure S2, top). Expression levels were quantified by densitometric measurements in regions of interest confined to A1 and MG. Raw grayscale intensity of the three target genes (Slc32a1, Slc17a6, Slc17a7) was background corrected and normalized by Gapdh grayscale intensity, which did not change significantly during development in either region by RNAseq or ISH.

These 4 genes are particularly useful for validation as they have distinct patterns of expression in A1 and MG and a documented maturational time-course. The housekeeping gene (Gapdh) is widely expressed in all neurons. It had a flat maturational trajectory and was used for normalization of ISH for the other genes. Slc32a1 (aka VGAT) is expressed at moderate levels in A1 and low levels in the MG. Slc17a7 (aka VGluT1) is expressed at high levels in cortex and low in the MG, whereas the expression of Slc17a6 (aka VGluT2) is complementary in these structures. Additional file 8: Figure S2 (bottom) contains plots of expression levels comparing quantification of the sequencing and in situ hybridization. These data indicate good agreement between methods with respect to both regional differences in expression and the maturational trajectories.

Functional gene set analyses (GSEA)

Functional gene set enrichment analysis (GSEA) is widely used to characterize enrichment of functionally-related sets of genes in a sample [54]. In this study, GSEA was used to rank genesets by enrichment magnitude and indicate whether the geneset was up- or down-regulated. GSEA was conducted on geneset listings drawn from two sources as of September 2014: (1) the Gene Families database maintained by the HUGO Gene Nomenclature Committee at the European Bioinformatics Institute (http://www.genenames.org) (HUGO) [55]; and (2) Mouse Genome Informatics (MGI) Gene Ontology Browser, maintained by Jackson Laboratories (www.informatics.jax.org) [56]. The HUGO database organizes the genome by gene family (e.g., ion channels, receptors, zinc finger proteins, etc.). The MGI database organizes genes into functional categories, where each geneset may contain genes from multiple gene families. Both of these databases are constructed and updated by consortium contributors based on review of the primary literature.

GSEA was applied to 111 gene families from the HUGO database and 51 Gene Ontology (GO) categories from the MGI database. Categories were selected for relevance to brain development and structure, synaptic transmission, and synaptic plasticity. Additional file 9: Table S21 contains the normalized counts of all samples for the 19,826 genes currently listed in the entire HUGO database, organized alphabetically by gene name. From this listing, a subset of 1557 genes within 111 gene families related to brain maturation and function were used to generate Additional file 9: Table S22, which contains the complete GSEA results for these gene families, ordered by FDR value. From the MGI database, Additional file 9: Table S23 contains normalized counts for 1402 genes distributed within the 51 GO categories selected (note that some genes are members of more than one GO category). The complete GSEA results are contained in Additional file 9: Table S24.

Construction of a gene families database

A major goal of this study was to develop a repository of the entire dataset formatted in a manner that simplified the screening and extraction of data at the global and single gene levels. The intent was to enable users to identify regional and age-related expression in genes of interest without extraction and analysis of the raw data (although the raw data are also available for such purposes). One of the most extensive and accessible resources provided is organized by gene family, and is contained within a single file (Additional file 10: Table S25). Approximately 4700 genes within 237 gene families listed in the HUGO database are organized into 20 functional groups, segregated by tabs. For each gene family, the normalized read counts of each member gene are tabulated and plotted as a function of postnatal age and brain region. There are 145 graphs, each plotting normalized read counts in A1 and MG as a function of postnatal age. The unique format of Additional file 10: Table S25 permits rapid inspection of the maturation trends for individual genes and gene families for both brain regions. As a guide to the use of this resource, an index with instructions is included under the first tab entitled, “Read Me + Index”. Two of the gene families contained in Additional file 10: Table S25 are highlighted in the Results section below.

Look-Up tool for generating maturational profiles at the single gene level

To facilitate screening and extraction of profiles from the database, a Look-Up tool was developed (Additional file 11). The tool automatically plots the maturational profiles and correlation matrices for any single gene or list of genes (up to 25 at a time). It also generates a listing of the normalized counts for all samples for extraction for other purposes.

Results

Data quality

RNAseq data were obtained from 48 samples and quality controlled. Samples information (sample ID, brain region, age, sex, and quality assessments) are contained in Additional file 2: Table S1. On average, sequencing produced 33.8 million reads per sample (range: 27.6–45.1 million). Sample 10 failed sequencing with less than half million reads produced, and sample 41 had relatively low read counts. Both were removed from subsequent analyses. No other quality issue was observed. The raw data statistics are contained in Additional file 2: Table S2. Alignment quality control was conducted, revealing that an average of 77.19 % of all reads (range: 51.86–83.01 %) were aligned to coding RNA regions (Additional file 2: Table S3). The complete raw read count information can be found in Additional file 2: Table S4.

Hierarchical clustering and differential expression analyses

Comparative transcriptomic analyses over all samples indicated that global gene expression patterns varied by region and postnatal age. Unsupervised hierarchical clustering analysis (Fig. 1) revealed several global trends. First, gene expression patterns were dominated by age at P7, but by brain region from P14 to adult. At P7, samples were clustered by age, then by brain region. From P14 to adult, samples were almost perfectly separated into two large clusters by region, and then by age within each regional cluster. Only a single sample (P21, A1) clustered with another age group (Adult, A1). By comparison, samples did not cluster by sex within any brain region or postnatal age. Overall, then, clustering was dominated by brain region and postnatal age.

To further explore these observations, differential expression analyses were systematically conducted comparing brain region and postnatal age (Figs. 2 and 3). Complete results for all comparisons (including fold change and p-values for all genes) are contained in Additional files 3, 4, 5, 6 and 7: Tables S5-S20. These analyses revealed several trends, described in the next two sections.

Fig. 2

Differential expression in A1 and MG. a The total numbers of differentially expressed genes between A1 and MG are plotted for each postnatal age. b Overlapping differential expression in A1 and MG. The Venn diagram depicts the total numbers of genes that were differentially expressed (MG vs A1) at only one postnatal age, and the numbers that were commonly expressed in all age group combinations. See text for proportions

Fig. 3

Differential expression between age groups. The total numbers of differentially-expressed genes are plotted for each of the six possible comparisons. Comparisons with P7 yielded the largest numbers of differentially expressed genes, and totals declined with increasing age

Differential expression between brain regions

First, comparing A1 and MG, regional differences in expression were substantial at all ages (Fig. 2a). The total numbers of differentially expressed genes were comparable from P7 to P21 (P7: 6773; P14: 7056; P21: 6629), then dropped to nearly half by adulthood (Adult: 3613). This indicates that regional differences in gene expression are greatest during postnatal development, but remain significant in adulthood. The Venn diagram in Fig. 2b depicts the total numbers of differentially expressed genes that were unique to each age, and those that were also differentially expressed in at least one other. The totals indicate that about 20 % of differentially expressed genes were unique to one age group, while a majority of those identified was common to more than one age group (P7, 67 %; P14, 83 %; P21, 82 %; Adult, 91 %; overall, 80 %). Fewer genes were differentially expressed in more than 2 age groups, however (3 ages, 10 %; all ages, 6 %). Overall, these data reflect robust regional differences in gene expression at all ages, and account for much of the regional clustering observed in Fig. 1.

In Table 1, the top 50 differentially expressed genes between MG and A1 are listed by postnatal age. Ranking was based on the p-values returned by DESeq2, EdgeR, and BaySeq analyses. Genes in the first four columns (left) were more highly expressed in the MG, whereas genes in the next four columns (right) were more highly expressed in A1. Several notable trends were observed. First, the genes in each listing (column) represent multiple gene families. Rarely was a single gene family represented more than once. This indicates that regional differences in gene expression were broadly distributed across multiple gene families. Second, about one-third of the genes listed in each column (i.e., age group) were also listed in at least one other (MG > A1, 33 %; A1 > MG, 35 %). Comparing P7 with adult in both regions yielded the least number of duplicated genes (N = 8, 16 %), suggesting greater diversity in the most highly expressed genes for that age interval compared with the others. Only two genes were more highly expressed in the MG than A1 at all ages (Slitrk6, Vav3). Four genes had higher expression in A1 than MG at all ages (Met, Efcab6, Hs3st2, Scube1). These genes are unique in that they ranked among the top differentially expressed genes between regions in the entire transcriptome across the entire age range. Of these, Slitrk6 and Met have been intensively studied and found to be critical for normal development in the forebrain [57–63].

Table 1 Top 50 differentially expressed genes in MG and A1

Differential expression between age groups

To further elucidate the age-related changes in A1 and MG, differential expression analyses were conducted between age groups within each region. In Fig. 3 and Table 2 (top panel), the numbers of differentially expressed genes are given for comparisons of each age with all other ages. Several trends were observed. First, comparisons of all other ages with P7 yielded the largest numbers of differentially expressed genes. The fewest differentially expressed genes were found in comparisons between older animals. Similarly, for the comparisons between successive age groups (i.e., P7-P14, P14-P21, P21-Adult), the number of differentially expressed genes declined steadily with increasing age in both regions. These trends reflect the stabilization of gene expression levels with increasing age. Second, with one exception (P7-Adult), the numbers of differentially expressed gene in MG and A1 was comparable for each age comparison. This suggests that gene expression matures at a comparable rate in both regions.

Table 2 Differential expression by age group and brain region

In the middle and bottom panels of Table 2, the totals reflect the number of genes that were differentially expressed at more than one age interval. That is, how many of the differentially expressed genes identified from comparisons between each age group (e.g., Fig. 3) were also differentially expressed at the other age intervals? Trends in these data are less obvious, but two observations are worth noting. First, substantial numbers of the same genes are differentially expressed in more than one age interval (A1 range: 275–8473; MG range: 292–6203). Secondly, these totals were lowest for comparisons involving the older age groups (e.g., P21, adult), and declined with age in a manner that was proportional to the numbers of differentially expressed genes compared. This is most clearly visualized from interactions with P21-Adult (last column), where totals in both brain regions reach minimum values. Overall, these findings support the general conclusion that changes in gene expression decline with age in both regions.

In Tables 3 and 4, the top 50 up-regulated and down-regulated genes in A1 and MG are listed for each of 4 maturation intervals (P7-Adult, P7-P14, P14-P21, P21-Adult). Rankings were based on the p-values from the 3 differential expression analysis methods (see methods). Although many hundreds of genes had increasing or decreasing trajectories (see Fig. 4), these truncated listings are instructive in that they reflect patterns in the regional and maturational changes observed.

Table 3 Top 50 up-regulated genes between age groups in A1 and MG

Table 4 Top 50 down-regulated genes between age groups in A1 and MG

Fig. 4

Expression trend analysis. The number of genes with increasing, decreasing, static or other maturational trajectories is plotted for A1 and MG. Also plotted are the numbers of genes with these profiles in both A1 and MG

First, note that the genes in each listing represent multiple gene families. Rarely was any single gene family represented more than once for a single comparison. Second, in all maturation intervals, a minority of genes was up- or down-regulated in both A1 and MG (range: 5 to 19). This matches global trends in Fig. 4 and is a further indication of regional specificity in expression trends. Third, relatively few genes were up- or down-regulated in more than one maturational interval (range: 1 to 12). This indicates specificity between age intervals in genes that are differentially expressed. As an example, only 1 gene (Plekhb1) was in the top 50 up-regulated genes of all 3 consecutive age intervals (MG, but not A1). Only 2 genes were down-regulated in A1 in all 3 consecutive age intervals (Cd24a, Nrep).

Expression trend analyses

Inspection of expression levels by postnatal age at the single gene level revealed that their maturational trajectories from P7 to adulthood had different profiles. To capture the main patterns, expression trend analyses were carried out to identify and tally genes with four different profiles types: monotonically increasing or decreasing, static, and other (Fig. 4, Table 5). A profile was monotonically “increasing” if expression increased successively at each time point and the change between P7 and adult was statistically significant. Monotonically “decreasing” genes were defined in the same fashion, but with decreased expression at each time point. Genes with flat trajectories across all ages were defined as “static”, and those with other patterns of expression (e.g., increasing, then decreasing or decreasing, then increasing) were categorized as “other”. The total numbers of genes with monotonically increasing or decreasing profiles was comparable in A1 (15.3 %) and MG (20.4 %). Of these, nearly equal numbers of genes had increasing and decreasing trajectories in A1, whereas 85.2 % of genes in MG had decreasing trajectories. In comparison to the monotonically changing profile types, the numbers of genes with “static” or “other” profiles were much greater, and similar in both regions. As will be noted in Figs. 5, 6, 7 and 8, a frequently observed profile in the “other” category was characterized by upregulation between P7 and P14 or P14 and P21, followed by downregulation at a subsequent age. Finally, in the third data series (A1 | MG), the number of genes that were differentially expressed in both A1 and MG (i.e., common to both regions) was given for each profile type. These numbers were a variable fraction (between 15 % and 64 %) of the total numbers in either region, depending on the profile. A possible interpretation is that expression of genes with the same maturational trajectory in both regions may be governed by similar factors.

Table 5 Expression trend analyses

Fig. 5

Gene expression profiles of the synaptic vesicle exocytosis gene ontology category. Gene expression profiles are plotted for a subset of genes from one gene ontology category, selected from the GSEA analysis in Table 7 (GO: 0016079, synaptic vesicle exocytosis). For each gene, mean normalized counts and % of maximum counts are plotted by postnatal age (P7, P14, P21, Adult) and brain region (A1, MG). Expression trajectory is indicated by arrows (up, down, none). Arrows were included only when differential expression from P7-Adult was significant (p < 0.05) by all three methods (DEseq2, EdgeR, and Bayseq)

Fig. 6

Gene expression profiles of the extracellular matrix proteoglycan family. Gene expression profiles are plotted for a subset of genes from the proteoglycan gene family, which contribute to formation of the extracelluar matrix. For each gene, mean normalized counts and % of maximum counts are plotted by postnatal age (P7, P14, P21, Adult) and brain region (A1, MG). Expression trajectory is indicated by arrows (up, down, none). Arrows were included only when differential expression from P7-Adult was significant (p < 0.05) by all three methods (DEseq2, EdgeR, and Bayseq)

Fig. 7

Gene expression profiles of 4 neurotransmitter receptor families. Gene expression profiles are plotted for selected genes from 4 receptor families with roles in neurotransmission and neuromodulation (glutamate, GABA, acetylcholine, serotonin). For each gene, mean normalized counts and % of maximum counts are plotted by postnatal age (P7, P14, P21, Adult) and brain region (A1, MG). Expression trajectory is indicated by arrows (up, down, none). Arrows were included only when differential expression from P7-Adult was significant (p < 0.05) by all three methods (DEseq2, EdgeR, and Bayseq)

Fig. 8

Gene expression profiles of a custom gene ontology category. A subset of genes with established roles in critical period formation in the visual cortex [64, 15] is profiled for A1 and MG. The listing spans multiple gene families. For each gene, mean normalized counts and % of maximum counts are plotted by postnatal age (P7, P14, P21, Adult) and brain region (A1, MG). Expression trajectory is indicated by arrows (up, down, none). Arrows were included only when differential expression from P7-Adult was significant (p < 0.05) by all three methods (DEseq2, EdgeR, and Bayseq)

Functional gene set analyses (GSEA)

Functional Gene Set Analysis (GSEA) were performed on 111 HUGO gene families (1557 genes) and 51 MGI gene ontology categories. The data from both resources are described separately below.

GSEA of MGI gene ontology categories

The complete MGI database contains scores of gene ontology (GO) categories arranged a priori by function, rather than by gene family. We performed GSEA on 51 GO categories selected for relevance to brain development, neurotransmission and plasticity (Additional file 9: Table S24). Of these, 20 had a false discovery rate (FDR) of q < 0.25 (Table 6). The majority of these had upward maturational trajectories, and mainly included genes related to synaptic plasticity and transmission. Among the categories with downward trajectories were those that included genes involved in cell migration, layer formation, axon extension, regionalization, and cell proliferation. The complete listing of genes (including counts) within the MGI GO categories chosen for this study is located in Additional file 9: Table S23.

Table 6 GSEA of selected MGI gene ontology categories

GSEA of HUGO gene families

From the 111 HUGO gene families (1557 genes) analyzed with GSEA (Additional file 9: Table S22), 27 families had a FDR of q < 0.25 (Table 7). Of these, 5 families had upward maturational trajectories, and contain genes related to neurotransmission or neuronal activity. Of the 22 families with downward trajectories, most are related to intracellular and extracellular structure. Note, however, that Additional file 9: Table S22 contains many gene families that did not reach the 25 % FDR criterion, even though they may contain several genes that were highly expressed, (e.g., all neurotransmitter receptor families other than adrenergic). This typically occurred when multiple members of that family were expressed at low or nominal levels. Similarly, the developmental trajectory of a gene family may not represent all members of that family. Thus, categorization of a family by trajectory or ranking by GSEA may not reflect the profiles of all genes in that family. The complete listing of genes (including counts) within the MGI GO categories chosen for this study is located in Additional file 9: Table S21.

Table 7 GSEA of selected HUGO gene family categories

Geneset profiling at the single gene level

To fully characterize the data contained within all relevant GO categories and gene families at the single gene level exceeds the scope of any single paper, as many thousands of genes are involved. As an alternative, we created searchable database tables and a look-up tool to permit viewing of the profiles of any single gene or geneset. In addition, the gene families database contains maturational profile plots for about 4700 genes within 237 selected gene families, organized into a searchable database for rapid screening of any gene family (Additional file 10: Table S25).

Examples of these profiles were generated from selected genes within one GO category (GO:0016079, synaptic vesicle exocytosis) and two gene families (extracellular matrix proteoglycans; neurotransmitter receptors) that were found to be enriched by GSEA (Figs. 5, 6 and 7). These three genesets were chosen because they are involved in different aspects of brain maturation (structure and function), and also exemplify the type of information contained within the database. Further, since none of these genesets have previously been profiled in the developing auditory forebrain, the data are both novel and informative. Finally, the profiles selected for illustration are typical of the regional and age-related diversity observed among members of the genesets analyzed in this study, and highlight the importance of evaluating expression patterns individually.

In each of these figures (Figs. 5, 6 and 7), normalized counts (mean and % of maximum) are plotted for each gene by postnatal age and brain region. The mean counts (white bars) convey information about expression magnitude of each gene by age. The % of maximum values (colored bars) facilitate visualization of the maturational trajectories, which are difficult to resolve when genes with high and low expression levels are plotted together. Arrows denote whether expression from P7 to adult was significantly up- or down-regulated.

Synaptic vesicle exocytosis (GO: 0016079)

Synaptic vesicle exocytosis is the process by which membrane-bound vesicles containing neurotransmitters are directed to their contents at a neuronal synapse. From this GO category (37 genes), the profiles of 15 were plotted in Fig. 5. Overall, this group was up-regulated (Table 6) from P7 to Adult, as might be expected following the onset of auditory experience. At the single gene level, however, profiles were mixed. In A1, 12 of 15 had upward profiles, 2 were static, and only one was significantly down-regulated (Unc13b). In MG, only 1 gene was significantly up-regulated (Cplx1), 4 were down-regulated, and the remaining 10 were static. The diversity of maturational profiles within this GO category highlights the importance of evaluating trends in expression at the single gene level, as well as by group or family.

One additional comment should be made here before considering other genesets. Note that genes with overall static profiles may exhibit significant changes in expression between one or more age comparisons. For example, the asterisk comparing P7 with P14 for Syt1 in the MG denotes a significant change for that age interval, whereas comparisons between other age groups were not significant. We did not calculate all 6 of the possible age comparisons for all of the genes profiled in Figs. 5, 6, 7 and 8, but included the Syt1 example here to draw attention to the potential for such changes at the single gene level of analysis. Also, recall that we used a very strict criterion (the difference must be significant by DESeq2, EdgeR and BaySeq methods), which results in the identification of fewer significant changes.

Extracellular matrix: proteoglycan family

The proteoglycans are a large class of glycoproteins that contribute to formation of the extracellular tissue matrix surrounding neurons and glia in the brain. Several other gene families are also involved in the formation and maintenance of the extracellular matrix (e.g., collagens, contactins, cadherins, laminins, neural cell adhesion molecules) (see Additional file 10: Table S25). In Fig. 6, 10 proteoglycan genes are profiled. As a family, the proteoglycans were downregulated from P7 to adult, but the profiles of individual genes were diverse. Comparing regions, expression levels and trajectories were fairly symmetric for this set of genes. 5 genes were significantly downregulated in one or both regions (Ncan, Vcan, Agrn, Sdc3, Gpc1). Two genes were significantly upregulated in both regions (Sdc4, Spock2), and two others had disparate regional profiles (Gpc1, Spock3). Acan had nominal levels of expression, which was unexpected based on prior studies of visual and somatosensory cortex (see Discussion).

Receptor families with roles in neurotransmission and neuromodulation

Several major classes of neurotransmitters and neuromodulators are involved in signaling between networks of neurons in the brain (e.g., glutamate, GABA, glycine, acetylcholine, dopamine, serotonin, noradrenaline). Multiple receptor types are associated with each class, forming a large and diverse collection of proteins. None of these families reached the 25 % FDR criterion by GSEA (see Additional file 9: Table S22). Their relatively low ranking in this listing was due to the mixed expression profiles of the genes in these families. Profiled in Fig. 7 are 5 representative genes from each of 4 receptor families with upward (GABA, glutamate) and downward (serotonin, acetylcholine) trajectories by GSEA. The genes were selected to represent the regional and maturational diversity within each family, and for the neurotransmitter receptors, in general. As for the gene families described above, expression levels and maturational trajectories typically varied between genes and brain region. For the GABA and glutamate receptors, expression levels and trajectories were comparable between regions. That is, most of the genes had symmetric expression levels in A1 and MG, and the maturational trajectories of about half were in the same direction (Gabra1, Gabra4, Gabrb2, Grm3, Grm5). The remainder had profiles that were opposing (Gria2, Grin2b) or mixed (Gabrb3, Gabbr1, Gria1). In contrast, there was greater diversity among the serotonin and acetylcholine receptors. For example, expression levels between A1 and MG were asymmetric for most of the genes in these two families. Maturational trajectories were also in opposite directions for three of these genes (Htr2a, Chrm3, Chrm4) and non-matching for four others (Htr1a, Htr5a, Chrm1, Chrnb2). Overall, these examples illustrate the diversity of expression patterns within each of the receptor families, and reveal that multiple receptor types from several receptor families are expressed in the same brain region. The functional roles of relatively few have been studied in detail, and many are contained within functional GO categories from which their function might be inferred. See Additional file 10: Table S25 to view the profiles of all of the genes within 7 neurotransmitter and 12 neuropeptide receptor families (“Receptors” tabs).

Custom gene ontology categories

Profiling by gene family is a convenient means of determining which genes within a possibly large family are expressed in the sampled region of interest and how their expression changes with age or some other manipulation. One limitation of this approach is that functionally-related sets of genes belong to multiple families and must be pieced together by additional analysis. By profiling established functional GO categories, genes that are functionally-related (based on prior experimental work) can be explored individually and as a group. Upon closer inspection of the GO categories profiled in Tables 6 and Additional file 9: Table S24, we noticed several trends that could be considered as caveats. First, there is typically some overlap across categories. That is, multiple categories index the same or similar functions, and the same gene or set of genes may be listed in several categories. Second, some GO categories contain a small number of genes (e.g., less than 10), or appear to be constructed from a limited or selective sampling of the literature. Third, the listings are almost always compiled from studies of other brain regions (i.e., non-auditory), some of which have very different patterns of organization (e.g., hippocampus). Thus, for some purposes, it may be advantageous to generate a custom GO category from an established model or a focused literature review.

Profiled in Fig. 8 is a custom GO listing of 16 genes associated with the opening and closing of critical periods plasticity in the visual cortex, based on the model advanced by Takesian and Hensch [64]. Because a comparable model does not exist for auditory cortex, we used this established model to probe our dataset. In this context, plasticity refers to the capacity for structural and functional change in some part of the brain (e.g., synapse, circuit, network), as regulated by intrinsic mechanisms or extrinsic factors, such as the onset of sensory experience. Critical periods denote periods of time during which the capacity for plasticity is high, and the functional properties of a brain region can be strongly shaped by experience in a manner that has long-lasting or permanent effects [15]. The model described by Takesian and Hensch elucidates the molecular mechanisms associated with the opening and closing of critical periods, which are themselves subject to modification. Because regions of sensory cortex share many features of organization, application of this model to the auditory forebrain may have higher relevance than a GO category based on studies of other brain regions.

As observed for nearly all of the GO categories related to synaptic plasticity (see Table 6), most of the genes had an upward maturational trajectory in one or both regions. 12 of the 16 genes had upward trajectories in A1, but only 5 had upward profiles in both regions (Nptx1, Lynx1, Camk2a, Mag, Mobp) (Fig. 8). This appears to reflect regional differences in the genes that govern critical periods. The genes that were up-regulated in A1 are involved in the opening and closing of critical periods in the visual cortex. Nine of these genes were significantly up-regulated between P7 and P14 (Bdnf, Nptx1, Nptx2, Nptxr, Lynx1, Camk2a, Mapk1, Mag, Mobp), which correspond to ages before and after hearing onset, and the opening of a critical period for auditory plasticity [13]. One of the genes that was down-regulated in both regions (Ncam1, neural cell adhesion molecule) is involved in preventing precocious plasticity prior to the opening of the critical period in visual cortex. The sharp down-regulation of Ncam1 between P7 and P14 may signal a reduction in its role as an attenuator of plasticity in the auditory forebrain. Conversely, the significant up-regulation between P7 and P14 of two genes related to myelination (Mag, Mobp) appears to be related to the ultimate closure of the critical period, as myelin formation dampens plasticity.

To further probe this custom geneset, we used the 16 critical periods genes from Fig. 8 as seeds to generate a functional association network of known and predicted interactions using the GeneMANIA tool (http://genemania.org) [65, 66], based on gene ontology (GO) biological function annotations. The network illustrated in Fig. 9a includes the 16 critical periods genes (nodes with black circles) and 50 interacting genes (nodes with gray circles). The number of interacting genes included is a user-selected option. Connection type is denoted by line color (see legend), and strength (line thickness) is weighted by linear regression-based computations of the functional association data in the databases indexed. In addition to a dense plexus of connections, the 66 genes in this network were cross-listed in 91 GO categories. Figure 9b depicts the connections of 13 genes that were contained within one of these GO categories: regulation of synaptic transmission. This includes 3 critical periods genes from Fig. 8 (Bdnf, Ntrk2, Camk2a), and 10 interacting genes from the pathways analysis. The normalized counts of these genes are plotted below. Note the expected close association between Bdnf (brain-derived neurotrophic factor) and its tyrosine kinase receptor, TrkB (Ntrk2).

Fig. 9

Pathways analysis of critical periods genes. The 16 critical periods genes from Fig. 8 were used as seeds to generate a functional association network of known and predicted interactions using the GeneMANIA tool (http://genemania.org). Analysis based on gene ontology (GO) biological function annotations. a Network generated from the 16 critical periods genes (nodes with black circles) and 50 interacting genes (nodes with gray circles). Connection type is denoted by line color (see legend), and strength (line thickness) is weighted by linear regression-based computations of the functional association data in the databases indexed. b Connections of the 13 genes (out of 66) from the entire network that were listed in the GO category: regulation of synaptic transmission. The normalized counts of these genes are plotted below. c – d Detailed connections of 2 genes from panel b: JPH3 (junctophilin-3), CSPG5 (chondroitin sulfate proteoglycan 5). See text for details

In Fig. 9c-d, the connections of 2 genes from Fig. 9b are illustrated in detail: Jph3 (junctophilin-3), Cspg5 (chondroitin sulfate proteoglycan 5). Both have significant expression in A1 and MG, and are broadly connected with other nodes in the network. Jph3 expression increased significantly from P7 to adult in A1 and trended upward in MG. Figure 9c reveals interactions between Jph3 and 12 of the 13 genes (not Ntf3) within the GO category, as well as 56 of 66 genes across the entire network. Cspg5 expression decreased significantly in both auditory regions. Figure 9d reveals interactions between Cspg5 and all 13 genes in the GO category, as well as 54 of 66 genes in the network. To date, neither gene has been explicitly linked to critical periods plasticity and functional roles in the auditory forebrain are unknown. Studies of other brain regions, however, indicate that both genes are essential for normal neuronal function. Jph3 mutations are linked to neuropathological conditions in the brain, such as Huntington’s disease [67, 68]. Cspg5 plays a role in normal neuronal growth and differentiation [69]. Given these data and their position in the network, further investigation may reveal heretofore unknown roles for one or more of these genes in critical periods plasticity or a related function.

Overall, the analysis illustrated in Fig. 9 provides an example of how transcriptome profiles may be used in conjunction with pathways analysis to guide discovery and generate testable hypotheses. In addition to identification of novel interactions, single gene expression profiles permit identification of genes that are not expressed at significant levels in the brain region of interest. Regionally-specific profiling at the single gene level is essential given the significant differences in expression patterns between regions.

Discussion

In the present study, we set out to achieve two main goals. The first was to generate complete transcriptome profiles of A1 and MG during postnatal development using next generation sequencing of total RNA. The second was to construct an accessible database in a format that permits extraction and screening of the data for any purpose (Additional files 3, 4, 5, 6, 7, 9, 10 and 11: Tables S5 – S26). Overall, the analyses of global gene expression revealed significant differences between brain regions at all ages, and changes within both regions with postnatal age. Geneset enrichment analyses revealed how those changes manifested within functional categories and gene families. The differential expression and gene families databases permit screening and extraction of the entire dataset down to the single gene level, aided by application of a look-up and plotting tool. To further illustrate the functional relevance and potential applications of the dataset, some of the results are discussed in more detail below.

Regional differences in gene expression

Regional differences in gene expression are well known in the forebrain [70–72], including numerous genes that have been intensively studied for their roles in brain development [63, 73, 74]. Those roles include both structural (e.g., regionalization, axon guidance, cell differentiation, synapse formation), and functional (e.g., neurotransmission, synaptic plasticity) features, which vary significantly between brain regions. Much progress has been made, but characterization of these features is far from complete for most brain areas.

The present study is the first to profile the transcriptome of the auditory forebrain in any species at any age. In sequencing two anatomically interconnected regions at the same time from the same subjects, we were able to directly compare the magnitude and trajectory of expression in juvenile and mature animals at key time points before and after the onset of hearing. Both hierarchical clustering and differential analyses revealed clear differences in global gene expression between A1 and MG at all ages. While global differences between regions located in the telencephalon and diencephalon are not surprising or especially informative on their own, they reflect important differences in the underlying patterns at the group (gene ontology categories, gene families) and single gene levels. In the present study, regional asymmetries in expression levels and maturational trajectories were commonplace, even for genes of the same family or functional group. Therefore, accurate characterization of the functional roles played by single genes and groups of genes must also account for brain location. This is important, because the regional differences in gene expression are likely to subserve important differences in function. For example, the asymmetrical expression of neurotransmitter receptors in A1 and MG implies that excitatory, inhibitory, and modulatory inputs to each region are mediated by a unique blend of receptor types that variably influence activity [26, 18, 22]. In addition to unique receptor profiles, we also observed regional differences in other genes that directly impact activity (e.g., ion channels, calcium binding proteins). Thus, in addition to differences that are conferred by unique input connectivity, regional differences in function may also be influenced by a rather large set of other factors. Expression profiling provides a way to screen for these factors, and narrow the range of targets for further study.

In addition to differences between the major divisions of the auditory forebrain, we also observed that regional differences in gene expression exist between closely-related structures in the brain. As an example, one of the most widely expressed of the proteoglycans is aggrecan (Acan), which has multiple isoforms that are differentially distributed in cortex [75–78]. In studies of visual or somatosensory cortex, Acan is upregulated during development and involved in the regulation of critical periods [79, 80]. Matthews et al. (2002) [76] found that Acan mRNA levels peak in somato-motor cortex at about P21, where a small subset of cells expressed the gene. Higher expression was noted in subsets of neurons in subcortical motor nuclei and the cochlear nucleus. Ye and Maio (2013) used immunohistochemistry to study development of several perineuronal net components, including Acan, in mouse visual cortex from P10 to adulthood. They also found a steady increase in expression that reached a plateau around P28. Based on these studies, we expected to see substantial and increasing levels of Acan with age, but this gene was expressed at very low levels in A1 and MG at all ages (Fig. 6). By comparison, the other lecticans (Bcan, Ncan, Vcan) were expressed at relatively high levels that changed substantially with age. Turning to the Allen Brain Atlas for additional validation, we observed that Acan is expressed in sizable subpopulations of neurons in somatosensory, motor and retrosplenial cortex, moderately in visual cortex, but rarely in auditory cortex or MG. Otherwise, the atlas revealed that Acan is highly expressed in the thalamic reticular nucleus (TRN) and subpopulations of neurons in major auditory brainstem nuclei (IC, SOC, CN). Thus, Acan is expressed at relatively low levels in the auditory forebrain, compared to other central auditory nuclei and sensory cortical areas. More generally, such findings may signal the existence of numerous important differences in structure and function between brain regions – including different sensory systems. While such differences are reminders to use caution when using the findings of one brain region (e.g., visual cortex, hippocampus) to draw conclusions about another (e.g., auditory cortex), the differences offer exciting opportunities for discovery.

Maturational changes in genes related to neuronal activity and brain structure

Differential analyses of the entire dataset revealed broad-based maturational changes in gene expression in A1 and MG. Globally, samples in both regions were almost perfectly clustered by age group, reflecting distinct patterns of gene expression at each age. The differences between age groups remained robust, despite a steady decline in the total numbers of differentially expressed genes with age.

The greatest changes in global gene expression occurred between P7 and P14, spanning the period before and after eye and ear canal opening in mice. Presumably, some of the changes observed are linked to the onset of sensory experience, especially in gene families related to neurotransmission or synaptic plasticity [81–84]. Also anticipated were changes in genes involved in brain structure, which may be altered as the architecture of each region becomes established. Indeed, the GSEA analyses revealed genes involved in synaptic transmission and plasticity were up-regulated on the whole, while groups of genes involved in establishing brain structure were down-regulated. For example, 14 of the 27 gene families listed in Table 7 contain genes primarily related to structure. All of these families had downward maturational trajectories.

As noted above, however, the expression patterns of individual genes are not always in line with the group as a whole. That is, most groups contained genes with a blend of increasing, decreasing, static and other trajectories. In addition, the distinction between genes with structural and functional roles is not absolute, as multiple functions may be attributed to the products of the same gene. As an example, the proteoglycan gene family (highlighted in the results) is a large and diverse class of glycoproteins known to be major structural components of perineuronal nets (PNNs) and the extracellular matrix in the brain and other tissues [85–88, 76]. The most widely studied are those that carry chondroitin sulfate (CSPG) or heparin sulfate (HSPG) side chains (e.g., hyalectans or lecticans; glypicans; syndecans). In addition to their structural support roles, the proteoglycans have also been intensively studied for their roles in brain development and plasticity [89–93, 80, 94–97]. Note that three of these genes (Bcan, Vcan, Ncan) were cross-listed in Fig. 6 (proteoglycans) and Fig. 8 (critical periods). In cortex, PNNs surround specific subsets of neurons, most notably parvalbumin expressing interneurons [77, 98, 99], where they are involved in the transfer of Otx2, which is apparently essential for the opening and closing of critical periods in visual cortex [91, 32]. In the present study, the majority had downward developmental trajectories, consistent with prior evidence that these genes reach peak expression levels at around the time of birth [97, 85, 100, 101]. Exceptions were Sdc4 and the testicans (e.g., Spock2, Spock3), especially Spock2, which was expressed at high levels and increased significantly in both regions. Thus, down-regulation of gene expression may signify a decreasing demand for the production of structural proteins as the architecture is laid down, while their persistent expression at relatively modest levels may reflect alternative roles, such as structural maintenance or synaptic plasticity [15, 64].

Comparison with prior studies of the central auditory pathway

To date, transcriptomic and/or proteomic profiling of the central auditory system remains rather limited. Although no previous studies have focused specifically on the forebrain, it is fortunate that we can at least compare some of our results with studies of auditory brainstem nuclei, where activity-dependent mechanisms were a theme (i.e., postnatal development, hearing loss). These comparisons reveal interesting regional differences in expression within the central auditory pathway, overall, that can improve understanding of the underlying circuits and guide future research.

Kaltwasser et al. (2013) used a proteomic approach to profile over 1200 proteins in the superior olivary complex (SOC) and IC during development (P4 and P60) [34]. The number of differentially regulated proteins between the regions was high (>75 %), whereas less than 20 % had common regulation patterns. Among the upregulated proteins were synaptophysin (SYPH or SYP) and two synaptic vesicle proteins (SV2A, SV2B), which are involved in the regulation of neurotransmission. Our data indicated that all three were greatly upregulated between P7 and P14 in both A1 and MG, and especially in A1 (Additional file 4: Tables S9-S10). The timecourse corresponds to the period of increased synaptogenesis during the early postnatal period [102–104]. In contrast, the cytoskeletal protein Stathmin-1 (STMN1) and calpain-6 (CAN6, CAPN6), involved in cytoskeletal remodeling, were greatly downregulated in the SOC. By comparison, we also observed rapid downregulation of STMN1 between P7 and P14, followed by slower reductions thereafter. Curiously, CAPN6 was only detectable at nominal levels across the entire age range in A1 or MG in the present study, whereas several other calpain genes were strongly expressed in A1 and MG, and exhibited age-related changes in expression (CAPN1, CAPN2, CAPN3, CAPN5, CAPN7, CAPN15, CAPNS1) (Additional file 10: Table S25: Calcium Binding & EF Hand tab). As noted earlier, such trends are likely to reflect a fraction of the regional similarities and differences in gene expression that are present.

Ehmann et al. (2013) used microarray analysis to profile gene expression in the rat SOC at P0, P14, P16, and P25. The extensive dataset profiled some 2000 genes, which limits detailed comparisons with the present study within this manuscript [33]. Among the top upregulated genes between P4 and P25 in the SOC were Mog and Mobp, which are related to myelination. These were also significantly upregulated in A1 and MG in the present study. Several potassium channel genes were significantly upregulated in both studies (Kcna1, Kcna2, Kcnab3), whereas others were upregulated or flat in A1 and MG (Kcna4, Kcnb1, Kcnk2, Kcnh1, Kcnt2, Kcnv1). In addition, several of the upregulated SOC genes were expressed at only nominal levels in A1 or MG (Kcn15, Kcns3, Kcnk5, Kcnj8). These differences signify regional differences in the distribution of potassium channels between the SOC and auditory forebrain, and are also likely to reflect many other differences that we did not take time to compare here.

Finally, two studies profiled changes in gene expression after manipulations of activity in the central auditory pathways. Holt et al. (2005) used targeted microarrays to track transcript expression in the rat IC at 3, 21 and 90 days after bilateral deafening. Variable trajectories in expression were observed in GABA, glycine, glutamate, and serotonin receptor genes, among the many others profiled [35]. The expression of several genes increased after deafening (their nomenclature: GluR2, GABA-A A1, GABA-A B2, GABA-A B3, 5HT2C). By way of comparison with the present study (refer to Fig. 7), we found that GluR2 (Gria2) was steadily upregulated after hearing onset in A1, but downregulated in the MG; GABA-A A1 and B2 (Gabra1, Gabrb2) had upward trajectories from P7 to adult in A1 and MG; GABA-A B3 (Gabrb3) decreased in A1 and MG; and 5HT2C (Htr2c) increased between P7 and P14, then declined in A1 and MG. Although auditory activity, per se, could be considered a common variable between studies, comparison of the just a portion of the findings highlights the complex and rather unpredictable regional differences that could also be associated with activity-dependent regulation of the same genes. Clarkson et al. (2012) used microarrays and qPCR to track gene expression in the IC after unilateral lesions of the auditory cortex (affecting the cortico-tectal projections to both hemispheres) [36]. Their dataset cannot be succinctly summarized, but a general trend relevant to this discussion was that genes related to neurotransmission and synaptic growth (among many others) were downregulated in the ipsilateral IC and upregulated contralaterally, presumably reflecting changes in the balance of excitation and inhibition. Thus, in both studies, transcript profiling was sensitive to changes induced by the manipulation in each of the brain regions studied. Although not yet applied in this manner, whole transcriptome sequencing would be expected to be at least as sensitive as microarrays to such changes, with the added advantage that the entire transcriptome can be profiled.

Applications and future directions

RNAseq is a powerful tool for mRNA profiling and transcriptome analyses, with broad potential applicability in neurobiology [105]. Relatively small amounts of starting material (<10 ng) are sufficient to conduct whole transcriptome sequencing of discrete brain areas or cell populations. The reduction in sequencing costs, development of bioinformatics tools, and availability of genomic libraries add further to the attractiveness of this approach [106, 107]. As mentioned above, an important advantage of targeted profiling of selected genes using qPCR or microarrays is that data from the entire transcriptome is obtained. This removes limits on the identification of novel or unexpected changes, and broadens the scope of pathway analyses.

The dataset generated by this study comprises an extensive reference library that indexes the expression of any gene or gene family in the A1 and MG from P7 to adult. In addition to information about these structures during postnatal development, the dataset is also a rich source of information about mature animals. We envision several potential uses of this dataset by those interested in the structure and function of the auditory forebrain.

One application would be as a screening tool to guide hypothesis formation and streamline the design of neuroanatomical and neurophysiological studies [108]. The dataset provides a priori knowledge of the expression levels of genes and gene pathways that may be targeted for studies of expression and co-expression patterns in intact tissue sections. That is, what are the subpopulations of cells in A1 or MG that express the genes of interest, and to what extent are they colocalized, or not? We used this approach in our recent study of the vesicular transporters of glutamate and GABA in developing mouse auditory forebrain [53]. In addition to augmentation of anatomical libraries, that information could also be used to guide selection of transgenic models for neurophysiological and behavioral assays, and improve the specificity of optogenetic or pharmacological studies where a particular cell population or receptor type is targeted.

A second, and related, application is to provide a baseline for experimental studies (e.g., altered sound exposure during development, hearing loss, aging, other pathology) [109, 110, 35, 36]. Transcriptomic analyses of global or targeted gene expression are powerful means to identify genes that are changing the most (or the least). Relationships between members of different gene families and functional pathways can easily be extracted and incorporated into analyses of functional pathways and construction of models.

Finally, we also envision that transcriptome profiles could be used to conduct comparative genomic studies in other animal models, including humans. We have long been aware of species differences in gene and protein expression in the auditory pathways [111–113] – an observation that is in line with a growing number of studies in the brain [114–119, 72, 120]. Documentation of species differences is absolutely essential to make informed conclusions and predictions about the roles of particular genes, and we must be vigilant to consider such differences in the interpretation of and application of profiling data.

Caveats and limitations

GSEA analysis

GSEA fostered the identification of GO categories and gene families that were enriched in each region of interest. A general advantage of this approach is that well-developed analytical tools are available to identify robust patterns in large datasets, such as those produced by RNAseq. A potential weakness is that important differences in the expression of some genes or entire gene families may be overlooked because they were ranked low by GSEA.

Application of GSEA to GO categories permits profiling of the expression patterns among genes that were grouped into functional categories based on published data. In that sense, they are a useful guide to the identification and exploration of pathways involved in a particular process. One disadvantage is that the underlying data may not apply to all regions of interest (i.e., auditory forebrain). That is, the data used to develop a GO category may be derived from a brain region or cell population for which transcriptome profiles differ substantially from the current region of interest. The differences observed between A1 and MG are examples of this consideration. A second limitation is that GO categories are often incomplete and subject to change as new data are obtained. In that sense, they are not invariant descriptors of the genes and regulatory factors that contribute to a given function. Finally, a limitation of all such analyses is that genes are often multifunctional and may appear in multiple GO categories.

Application of GSEA to gene families facilitates identification of novel patterns and relationships between genes in the same family or related families. One limitation is that co-expression or co-variation of genes within a family does not imply that a functionally-significant relationship exists between those genes. Most functional pathways involve interactions between genes across multiple families. Secondly, since the members of gene families often have different expression profiles, some genes with functional relevance may be missed in the categorical GSEA analysis.

Thus, application of GSEA is useful for identification of major trends, but is best used in combination with tools that enable construction of custom GO categories and pathway analyses.

Single-gene profiling

An important goal of the analyses that were employed in this study was to distill some of most important expression patterns identified at the global transcriptome level into manageable gene-level assemblages with functional relevance. Unfortunately, an exhaustive treatment of the interesting trends and relationships contained within this dataset are rather far beyond the scope of any single paper. Instead, we highlighted subsets of the data through analyses of selected genesets, along with single gene analyses to fill in the details. In so doing, many hundreds of genes with significant functionality were not profiled in the main body of the manuscript. Fortunately, information about any gene in the library can be easily extracted and plotted from the supplementary Tables or by using the Look-Up tool.

One proviso concerns the interpretation of gene expression magnitude (read counts) in such data. For example, if a gene has low or even undetectable expression in a sample, or if changes in expression are not significant, several explanations may apply [121]: (1) the gene is expressed at very low levels across all cells in the sample; (2) the gene is expressed at moderate to high levels, but in a small subpopulation of cells; and (3) the somata containing the transcripts are located outside the region of harvest (e.g., endogenous ligands, vesicular transporters). Conversely, if a gene is expressed at high levels, comparable questions about specificity also apply. For example, it may be that a gene is expressed at high levels only within certain cell subpopulations or anatomical subdivisions within the sample. Therefore, interpretation of all results must be scrutinized within a valid anatomical framework. That would include, but not be limited to, knowledge of the brain region, layer or nucleus, and specific cell types that express each gene [122, 71, 123], as well as their patterns of co-expression with other related genes.

A second limitation is that, in their present form, our data do not directly address the identity of the elements that are regulating or driving changes in gene expression within a region. Are genetic factors driving the formation of functional networks? Are the developing networks or external stimuli driving de novo gene expression adaptively? Does a statistically-significant change in expression have functional significance? The biological importance of a regional or age-related change in gene expression must be determined by other means. Therefore, the data contained in this manuscript may best be used to drive hypotheses and inform the design of experiments that can provide answers to these questions.

Conclusions

Transcriptome profiling at the global, group, and single gene levels revealed that gene expression is in constant flux from P7 to adult in both A1 and MG. Globally, expression profiles are strongly clustered by brain region and age. The greatest changes in global gene expression occur between P7 and P14 – before and after the onset of hearing. Thereafter, differences between age groups declined with age, consistent with an increase in the stability of gene expression patterns with postnatal age. At the group level, functional GO categories and gene families were up- or down-regulated from P7 to adult in a manner consistent with their functional roles. Overall, genesets related to the establishment of brain structure (architecture) were down-regulated as the brain matured. In contrast, genesets involved with neuronal activity and plasticity were up-regulated, especially after the onset of hearing. At the single gene level, maturational profiles varied by brain region and age, and in a manner that was not always predicted by analyses at the group level. Although functional studies are lacking at present, the database generated by this study provides a foundation for the identification of pathways that operate specifically within the auditory forebrain.

Availability of Supporting Data

The complete set of raw sequencing files is available from the National Center for Biotechnology Information (NCBI) database under accession number SRP053237 (http://www.ncbi.nlm.nih.gov/projects/geo/). All other supporting data are included in the Additional files section.