Skip to main content

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 [25]. This window has provided researchers with a unique opportunity to document structural and functional maturation associated with the onset of hearing [69], and the formation of critical periods for the plasticity of sound feature encoding and behavior [1014]. 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 [1620]. Other studies have focused on the influences of neuromodulatory inputs (e.g., cholinergic, dopaminergic, serotonergic) [2126] 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−800 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−800 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 106). 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 [4852]. 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
figure1

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
figure2

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
figure3

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 [5763].

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
figure4

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
figure5

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
figure6

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
figure7

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
figure8

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
figure9

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 [7072], 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 [7578]. 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 [8184]. 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 [8588, 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 [8993, 80, 9497]. 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 [102104]. 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 [111113] – an observation that is in line with a growing number of studies in the brain [114119, 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.

References

  1. 1.

    Kang HJ, Kawasawa YI, Cheng F, Zhu Y, Xu X, Li M, et al. Spatio-temporal transcriptome of the human brain. Nature. 2011;478(7370):483–9. doi:10.1038/nature10523.

    PubMed Central  CAS  PubMed  Google Scholar 

  2. 2.

    Willott JF, Shnerson A. Rapid development of tuning characteristics of inferior colliculus neurons of mouse pups. Brain Res. 1978;148(1):230–3.

    CAS  PubMed  Google Scholar 

  3. 3.

    Ehret G. Development of absolute auditory thresholds in the house mouse (Mus musculus). J Am Audiol Soc. 1976;1(5):179–84.

    CAS  PubMed  Google Scholar 

  4. 4.

    Ehret G. Postnatal development in the acoustic system of the house mouse in the light of developing masked thresholds. J Acoust Soc Am. 1977;62(1):143–8.

    CAS  PubMed  Google Scholar 

  5. 5.

    Ehret G, Romand R. Development of tone response thresholds, latencies and tuning in the mouse inferior colliculus. Brain Res Dev Brain Res. 1992;67(2):317–26.

    CAS  PubMed  Google Scholar 

  6. 6.

    Froemke RC, Jones BJ. Development of auditory cortical synaptic receptive fields. Neurosci Biobehav Rev. 2011;35(10):2105–13. doi:10.1016/j.neubiorev.2011.02.006.

    PubMed Central  PubMed  Google Scholar 

  7. 7.

    Sanes DH, Bao S. Tuning up the developing auditory CNS. Curr Opin Neurobiol. 2009;19(2):188–99. doi:10.1016/j.conb.2009.05.014.

    PubMed Central  CAS  PubMed  Google Scholar 

  8. 8.

    Sanes DH, Woolley SM. A behavioral framework to guide research on central auditory development and plasticity. Neuron. 2011;72(6):912–29. doi:10.1016/j.neuron.2011.12.005.

    PubMed Central  CAS  PubMed  Google Scholar 

  9. 9.

    Bao S. Perceptual learning in the developing auditory cortex. Eur J Neurosci. 2015;41(5):718–24. doi:10.1111/ejn.12826.

    PubMed  Google Scholar 

  10. 10.

    Kral A. Auditory critical periods: a review from system's perspective. Neuroscience. 2013;247:117–33. doi:10.1016/j.neuroscience.2013.05.021.

    CAS  PubMed  Google Scholar 

  11. 11.

    de Villers-Sidani E, Merzenich MM. Lifelong plasticity in the rat auditory cortex: basic mechanisms and role of sensory experience. Prog Brain Res. 2011;191:119–31. doi:10.1016/B978-0-444-53752-2.00009-6.

    PubMed  Google Scholar 

  12. 12.

    Chun S, Bayazitov IT, Blundon JA, Zakharenko SS. Thalamocortical long-term potentiation becomes gated after the early critical period in the auditory cortex. J Neurosci. 2013;33(17):7345–57. doi:10.1523/JNEUROSCI.4500-12.2013.

    PubMed Central  CAS  PubMed  Google Scholar 

  13. 13.

    Barkat TR, Polley DB, Hensch TK. A critical period for auditory thalamocortical connectivity. Nat Neurosci. 2011;14(9):1189–94. doi:10.1038/nn.2882.

    CAS  PubMed  Google Scholar 

  14. 14.

    Yang EJ, Lin EW, Hensch TK. Critical period for acoustic preference in mice. Proc Natl Acad Sci U S A. 2012;109 Suppl 2:17213–20. doi:10.1073/pnas.1200705109.

    PubMed Central  CAS  PubMed  Google Scholar 

  15. 15.

    Hensch TK. Critical period plasticity in local cortical circuits. Nat Rev Neurosci. 2005;6(11):877–88. doi:10.1038/nrn1787.

    CAS  PubMed  Google Scholar 

  16. 16.

    Chang EH, Kotak VC, Sanes DH. Long-term depression of synaptic inhibition is expressed postsynaptically in the developing auditory system. J Neurophysiol. 2003;90(3):1479–88. doi:10.1152/jn.00386.2003.

    CAS  PubMed  Google Scholar 

  17. 17.

    Oswald AM, Reyes AD. Development of inhibitory timescales in auditory cortex. Cereb Cortex. 2011;21(6):1351–61. doi:10.1093/cercor/bhq214.

    PubMed Central  PubMed  Google Scholar 

  18. 18.

    Venkataraman Y, Bartlett EL. Postnatal development of synaptic properties of the GABAergic projection from the inferior colliculus to the auditory thalamus. J Neurophysiol. 2013;109(12):2866–82. doi:10.1152/jn.00021.2013.

    PubMed Central  CAS  PubMed  Google Scholar 

  19. 19.

    Sun YJ, Wu GK, Liu BH, Li P, Zhou M, Xiao Z, et al. Fine-tuning of pre-balanced excitation and inhibition during auditory cortical development. Nature. 2010;465(7300):927–31. doi:10.1038/nature09079.

    PubMed Central  CAS  PubMed  Google Scholar 

  20. 20.

    Dorrn AL, Yuan K, Barker AJ, Schreiner CE, Froemke RC. Developmental sensory experience balances cortical excitation and inhibition. Nature. 2010;465(7300):932–6. doi:10.1038/nature09119.

    PubMed Central  CAS  PubMed  Google Scholar 

  21. 21.

    Kilgard MP, Merzenich MM. Plasticity of temporal information processing in the primary auditory cortex. Nat Neurosci. 1998;1(8):727–31.

    PubMed Central  CAS  PubMed  Google Scholar 

  22. 22.

    Metherate R, Hsieh CY. Regulation of glutamate synapses by nicotinic acetylcholine receptors in auditory cortex. Neurobiol Learn Mem. 2003;80(3):285–90.

    CAS  PubMed  Google Scholar 

  23. 23.

    Froemke RC, Carcea I, Barker AJ, Yuan K, Seybold BA, Martins AR, et al. Long-term modification of cortical synapses improves sensory perception. Nat Neurosci. 2013;16(1):79–88. doi:10.1038/nn.3274.

    PubMed Central  CAS  PubMed  Google Scholar 

  24. 24.

    Hurley LM, Sullivan MR. From behavioral context to receptors: serotonergic modulatory pathways in the IC. Front Neural Circ. 2012;6:58. doi:10.3389/fncir.2012.00058.

    Google Scholar 

  25. 25.

    Bao S, Chan VT, Merzenich MM. Cortical remodelling induced by activity of ventral tegmental dopamine neurons. Nature. 2001;412(6842):79–83.

    CAS  PubMed  Google Scholar 

  26. 26.

    Edeline JM, Manunta Y, Hennevin E. Induction of selective plasticity in the frequency tuning of auditory cortex and auditory thalamus neurons by locus coeruleus stimulation. Hear Res. 2011;274(1–2):75–84. doi:10.1016/j.heares.2010.08.005.

    PubMed  Google Scholar 

  27. 27.

    Brown MR, Kaczmarek LK. Potassium channel modulation and auditory processing. Hear Res. 2011;279(1–2):32–42. doi:10.1016/j.heares.2011.03.004.

    PubMed Central  CAS  PubMed  Google Scholar 

  28. 28.

    Schachtele SJ, Losh J, Dailey ME, Green SH. Spine formation and maturation in the developing rat auditory cortex. J Comp Neurol. 2011;519(16):3327–45. doi:10.1002/cne.22728.

    PubMed Central  PubMed  Google Scholar 

  29. 29.

    Sutor B, Hagerty T. Involvement of gap junctions in the development of the neocortex. Biochimica et biophysica acta. 2005;1719(1–2):59–68. doi:10.1016/j.bbamem.2005.09.005.

    CAS  PubMed  Google Scholar 

  30. 30.

    O'Neil JN, Connelly CJ, Limb CJ, Ryugo DK. Synaptic morphology and the influence of auditory experience. Hear Res. 2011;279(1–2):118–30. doi:10.1016/j.heares.2011.01.019.

    PubMed Central  PubMed  Google Scholar 

  31. 31.

    Yu WM, Goodrich LV. Morphological and physiological development of auditory synapses. Hear Res. 2014;311:3–16. doi:10.1016/j.heares.2014.01.007.

    PubMed  Google Scholar 

  32. 32.

    Sugiyama S, Di Nardo AA, Aizawa S, Matsuo I, Volovitch M, Prochiantz A, et al. Experience-dependent transfer of Otx2 homeoprotein into the visual cortex activates postnatal plasticity. Cell. 2008;134(3):508–20. doi:10.1016/j.cell.2008.05.054.

    CAS  PubMed  Google Scholar 

  33. 33.

    Ehmann H, Hartwich H, Salzig C, Hartmann N, Clement-Ziza M, Ushakov K, et al. Time-dependent gene expression analysis of the developing superior olivary complex. J Biol Chem. 2013;288(36):25865–79. doi:10.1074/jbc.M113.490508.

    PubMed Central  CAS  PubMed  Google Scholar 

  34. 34.

    Kaltwasser B, Schulenborg T, Beck F, Klotz M, Schafer KH, Schmitt M, et al. Developmental changes of the protein repertoire in the rat auditory brainstem: a comparative proteomics approach in the superior olivary complex and the inferior colliculus with DIGE and iTRAQ. J Proteomics. 2013;79:43–59. doi:10.1016/j.jprot.2012.11.018.

    CAS  PubMed  Google Scholar 

  35. 35.

    Holt AG, Asako M, Lomax CA, MacDonald JW, Tong L, Lomax MI, et al. Deafness-related plasticity in the inferior colliculus: gene expression profiling following removal of peripheral activity. J Neurochem. 2005;93(5):1069–86. doi:10.1111/j.1471-4159.2005.03090.x.

    CAS  PubMed  Google Scholar 

  36. 36.

    Clarkson C, Herrero-Turrion MJ, Merchan MA. Cortical auditory deafferentation induces long-term plasticity in the inferior colliculus of adult rats: microarray and qPCR analysis. Front Neural Circuits. 2012;6:86. doi:10.3389/fncir.2012.00086.

    PubMed Central  CAS  PubMed  Google Scholar 

  37. 37.

    Hackett TA, Barkat TR, O'Brien BM, Hensch TK, Polley DB. Linking topography to tonotopy in the mouse auditory thalamocortical circuit. J Neurosci. 2011;31(8):2983–95. doi:10.1523/JNEUROSCI.5333-10.2011.

    PubMed Central  CAS  PubMed  Google Scholar 

  38. 38.

    Horng S, Kreiman G, Ellsworth C, Page D, Blank M, Millen K, et al. Differential gene expression in the developing lateral geniculate nucleus and medial geniculate nucleus reveals novel roles for Zic4 and Foxp2 in visual and auditory pathway development. J Neurosci. 2009;29(43):13672–83. doi:10.1523/JNEUROSCI.2127-09.2009.

    PubMed Central  CAS  PubMed  Google Scholar 

  39. 39.

    Guo Y, Ye F, Sheng Q, Clark T, Samuels DC. Three-stage quality control strategies for DNA re-sequencing data. Brief Bioinform. 2013. doi:10.1093/bib/bbt069.

    Google Scholar 

  40. 40.

    Guo Y, Zhao S, Sheng Q, Ye F, Li J, Lehmann B, et al. Multi-perspective quality control of Illumina exome sequencing data using QC3. Genomics. 2014;103(5–6):323–8. doi:10.1016/j.ygeno.2014.03.006.

    CAS  PubMed  Google Scholar 

  41. 41.

    Guo Y, Zhao S, Ye F, Sheng Q, Shyr Y. MultiRankSeq: multiperspective approach for RNAseq differential expression analysis and quality control. Biomed Res Int. 2014;2014:248090. doi:10.1155/2014/248090.

    PubMed Central  PubMed  Google Scholar 

  42. 42.

    Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg SL. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 2013;14(4):R36. doi:10.1186/gb-2013-14-4-r36.

    PubMed Central  PubMed  Google Scholar 

  43. 43.

    Anders S, Pyl PT, Huber W. HTSeq-a Python framework to work with high-throughput sequencing data. Bioinformatics. 2014. doi:10.1093/bioinformatics/btu638.

    PubMed Central  PubMed  Google Scholar 

  44. 44.

    Zhao S, Guo Y, Sheng Q, Shyr Y. Advanced heat map and clustering analysis using heatmap3. Biomed Res Int. 2014;2014:986048. doi:10.1155/2014/986048.

    PubMed Central  PubMed  Google Scholar 

  45. 45.

    Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010;11(10):R106. doi:10.1186/gb-2010-11-10-r106.

    PubMed Central  CAS  PubMed  Google Scholar 

  46. 46.

    Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26(1):139–40. doi:10.1093/bioinformatics/btp616.

    PubMed Central  CAS  PubMed  Google Scholar 

  47. 47.

    Hardcastle TJ, Kelly KA. baySeq: empirical Bayesian methods for identifying differential expression in sequence count data. BMC Bioinformatics. 2010;11:422. doi:10.1186/1471-2105-11-422.

    PubMed Central  PubMed  Google Scholar 

  48. 48.

    Dillies MA, Rau A, Aubert J, Hennequet-Antier C, Jeanmougin M, Servant N, et al. A comprehensive evaluation of normalization methods for Illumina high-throughput RNA sequencing data analysis. Brief Bioinform. 2013;14(6):671–83. doi:10.1093/bib/bbs046.

    CAS  PubMed  Google Scholar 

  49. 49.

    Kvam VM, Liu P, Si Y. A comparison of statistical methods for detecting differentially expressed genes from RNA-seq data. Am J Botany. 2012;99(2):248–56. doi:10.3732/ajb.1100340.

    Google Scholar 

  50. 50.

    Robles JA, Qureshi SE, Stephen SJ, Wilson SR, Burden CJ, Taylor JM. Efficient experimental design and analysis strategies for the detection of differential expression using RNA-Sequencing. BMC Genomics. 2012;13:484. doi:10.1186/1471-2164-13-484.

    PubMed Central  CAS  PubMed  Google Scholar 

  51. 51.

    Soneson C, Delorenzi M. A comparison of methods for differential expression analysis of RNA-seq data. BMC Bioinformatics. 2013;14:91. doi:10.1186/1471-2105-14-91.

    PubMed Central  PubMed  Google Scholar 

  52. 52.

    Guo Y, Li CI, Ye F, Shyr Y. Evaluation of read count based RNAseq analysis methods. BMC Genomics. 2013;14 Suppl 8:S2. doi:10.1186/1471-2164-14-S8-S2.

    PubMed  Google Scholar 

  53. 53.

    Hackett TA, Takahata T, Clause AR, Hackett NJ, Polley DB. Differential maturation of vesicular glutamate and GABA transporter expression in the mouse auditory forebrain during the first weeks of hearing. Brain Structure & Function. 2015:in press.

  54. 54.

    Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005;102(43):15545–50. doi:10.1073/pnas.0506580102.

    PubMed Central  CAS  PubMed  Google Scholar 

  55. 55.

    Gray KA, Yates B, Seal RL, Wright MW, Bruford EA. Genenames.org: the HGNC resources in 2015. Nucleic Acids Res. 2015;43(Database issue):D1079–85. doi:10.1093/nar/gku1071.

    PubMed Central  PubMed  Google Scholar 

  56. 56.

    Blake JA, Bult CJ, Eppig JT, Kadin JA, Richardson JE. Mouse genome database G. The mouse genome database: integration of and access to knowledge about the laboratory mouse. Nucleic Acids Res. 2014;42(Database issue):D810–7. doi:10.1093/nar/gkt1225.

    PubMed Central  CAS  PubMed  Google Scholar 

  57. 57.

    Matsumoto Y, Katayama K, Okamoto T, Yamada K, Takashima N, Nagao S, et al. Impaired auditory-vestibular functions and behavioral abnormalities of Slitrk6-deficient mice. PLoS One. 2011;6(1):e16497. doi:10.1371/journal.pone.0016497.

    PubMed Central  CAS  PubMed  Google Scholar 

  58. 58.

    Tekin M, Chioza BA, Matsumoto Y, Diaz-Horta O, Cross HE, Duman D, et al. SLITRK6 mutations cause myopia and deafness in humans and mice. J Clin Invest. 2013;123(5):2094–102. doi:10.1172/JCI65853.

    PubMed Central  CAS  PubMed  Google Scholar 

  59. 59.

    Aruga J, Mikoshiba K. Identification and characterization of Slitrk, a novel neuronal transmembrane protein family controlling neurite outgrowth. Mol Cell Neurosci. 2003;24(1):117–29.

    CAS  PubMed  Google Scholar 

  60. 60.

    Beaubien F, Cloutier JF. Differential expression of Slitrk family members in the mouse nervous system. Dev Dynamics. 2009;238(12):3285–96. doi:10.1002/dvdy.22160.

    CAS  Google Scholar 

  61. 61.

    Eagleson KL, Milner TA, Xie Z, Levitt P. Synaptic and extrasynaptic location of the receptor tyrosine kinase met during postnatal development in the mouse neocortex and hippocampus. J Comp Neurol. 2013;521(14):3241–59. doi:10.1002/cne.23343.

    PubMed Central  CAS  PubMed  Google Scholar 

  62. 62.

    Qiu S, Lu Z, Levitt P. MET receptor tyrosine kinase controls dendritic complexity, spine morphogenesis, and glutamatergic synapse maturation in the hippocampus. J Neurosci. 2014;34(49):16166–79. doi:10.1523/JNEUROSCI.2580-14.2014.

    PubMed Central  PubMed  Google Scholar 

  63. 63.

    O'Leary DD, Sahara S. Genetic regulation of arealization of the neocortex. Curr Opin Neurobiol. 2008;18(1):90–100. doi:10.1016/j.conb.2008.05.011.

    PubMed Central  PubMed  Google Scholar 

  64. 64.

    Takesian AE, Hensch TK. Balancing plasticity/stability across brain development. Prog Brain Res. 2013;207:3–34. doi:10.1016/B978-0-444-63327-9.00001-1.

    PubMed  Google Scholar 

  65. 65.

    Warde-Farley D, Donaldson SL, Comes O, Zuberi K, Badrawi R, Chao P, et al. The GeneMANIA prediction server: biological network integration for gene prioritization and predicting gene function. Nucleic Acids Res. 2010;38(Web Server issue):W214–20. doi:10.1093/nar/gkq537.

    PubMed Central  CAS  PubMed  Google Scholar 

  66. 66.

    Mostafavi S, Ray D, Warde-Farley D, Grouios C, Morris Q. GeneMANIA: a real-time multiple association network integration algorithm for predicting gene function. Genome Biol. 2008;9 Suppl 1:S4. doi:10.1186/gb-2008-9-s1-s4.

    PubMed  Google Scholar 

  67. 67.

    Holmes SE, O'Hearn E, Rosenblatt A, Callahan C, Hwang HS, Ingersoll-Ashworth RG, et al. A repeat expansion in the gene encoding junctophilin-3 is associated with Huntington disease-like 2. Nat Genet. 2001;29(4):377–8. doi:10.1038/ng760.

    CAS  PubMed  Google Scholar 

  68. 68.

    Wilburn B, Rudnicki DD, Zhao J, Weitz TM, Cheng Y, Gu X, et al. An antisense CAG repeat transcript at JPH3 locus mediates expanded polyglutamine protein toxicity in Huntington's disease-like 2 mice. Neuron. 2011;70(3):427–40. doi:10.1016/j.neuron.2011.03.021.

    PubMed Central  CAS  PubMed  Google Scholar 

  69. 69.

    Nakanishi K, Aono S, Hirano K, Kuroda Y, Ida M, Tokita Y, et al. Identification of neurite outgrowth-promoting domains of neuroglycan C, a brain-specific chondroitin sulfate proteoglycan, and involvement of phosphatidylinositol 3-kinase and protein kinase C signaling pathways in neuritogenesis. J Biol Chem. 2006;281(34):24970–8. doi:10.1074/jbc.M601498200.

    CAS  PubMed  Google Scholar 

  70. 70.

    Lein ES, Hawrylycz MJ, Ao N, Ayres M, Bensinger A, Bernard A, et al. Genome-wide atlas of gene expression in the adult mouse brain. Nature. 2007;445(7124):168–76. doi:10.1038/nature05453.

    CAS  PubMed  Google Scholar 

  71. 71.

    Ko Y, Ament SA, Eddy JA, Caballero J, Earls JC, Hood L, et al. Cell type-specific genes show striking and distinct patterns of spatial expression in the mouse brain. Proc Natl Acad Sci U S A. 2013;110(8):3095–100. doi:10.1073/pnas.1222897110.

    PubMed Central  CAS  PubMed  Google Scholar 

  72. 72.

    Bernard A, Lubbers LS, Tanis KQ, Luo R, Podtelezhnikov AA, Finney EM, et al. Transcriptional architecture of the primate neocortex. Neuron. 2012;73(6):1083–99. doi:10.1016/j.neuron.2012.03.002.

    PubMed Central  CAS  PubMed  Google Scholar 

  73. 73.

    Borello U, Pierani A. Patterning the cerebral cortex: traveling with morphogens. Curr Opin Genet Dev. 2010;20(4):408–15. doi:10.1016/j.gde.2010.05.003.

    CAS  PubMed  Google Scholar 

  74. 74.

    Rash BG, Grove EA. Area and layer patterning in the developing cerebral cortex. Curr Opin Neurobiol. 2006;16(1):25–34. doi:10.1016/j.conb.2006.01.004.

    CAS  PubMed  Google Scholar 

  75. 75.

    Virgintino D, Perissinotto D, Girolamo F, Mucignat MT, Montanini L, Errede M, et al. Differential distribution of aggrecan isoforms in perineuronal nets of the human cerebral cortex. J Cell Mol Med. 2009;13(9B):3151–73. doi:10.1111/j.1582-4934.2009.00694.x.

    PubMed Central  PubMed  Google Scholar 

  76. 76.

    Matthews RT, Kelly GM, Zerillo CA, Gray G, Tiemeyer M, Hockfield S. Aggrecan glycoforms contribute to the molecular heterogeneity of perineuronal nets. J Neurosci. 2002;22(17):7536–47.

    CAS  PubMed  Google Scholar 

  77. 77.

    Ye Q, Miao QL. Experience-dependent development of perineuronal nets and chondroitin sulfate proteoglycan receptors in mouse visual cortex. Matrix Biol. 2013;32(6):352–63. doi:10.1016/j.matbio.2013.04.001.

    CAS  PubMed  Google Scholar 

  78. 78.

    Morawski M, Bruckner G, Arendt T, Matthews RT. Aggrecan: Beyond cartilage and into the brain. Int J Biochem Cell Biol. 2012;44(5):690–3. doi:10.1016/j.biocel.2012.01.010.

    CAS  PubMed  Google Scholar 

  79. 79.

    Hockfield S, Kalb RG, Zaremba S, Fryer H. Expression of neural proteoglycans correlates with the acquisition of mature neuronal properties in the mammalian brain. Cold Spring Harb Symp Quant Biol. 1990;55:505–14.

    CAS  PubMed  Google Scholar 

  80. 80.

    Sur M, Frost DO, Hockfield S. Expression of a surface-associated antigen on Y-cells in the cat lateral geniculate nucleus is regulated by visual experience. J Neurosci. 1988;8(3):874–82.

    CAS  PubMed  Google Scholar 

  81. 81.

    Vazdarjanova A, McNaughton BL, Barnes CA, Worley PF, Guzowski JF. Experience-dependent coincident expression of the effector immediate-early genes arc and Homer 1a in hippocampal and neocortical neuronal networks. J Neurosci. 2002;22(23):10067–71.

    CAS  PubMed  Google Scholar 

  82. 82.

    Yuge K, Kataoka A, Yoshida AC, Itoh D, Aggarwal M, Mori S, et al. Region-specific gene expression in early postnatal mouse thalamus. J Comp Neurol. 2011;519(3):544–61. doi:10.1002/cne.22532.

    CAS  PubMed  Google Scholar 

  83. 83.

    Sato H, Fukutani Y, Yamamoto Y, Tatara E, Takemoto M, Shimamura K, et al. Thalamus-derived molecules promote survival and dendritic growth of developing cortical neurons. J Neurosci. 2012;32(44):15388–402. doi:10.1523/JNEUROSCI.0293-12.2012.

    CAS  PubMed  Google Scholar 

  84. 84.

    Lyckman AW, Horng S, Leamey CA, Tropea D, Watakabe A, Van Wart A, et al. Gene expression patterns in visual cortex during the critical period: synaptic stabilization and reversal by visual deprivation. Proc Natl Acad Sci U S A. 2008;105(27):9409–14. doi:10.1073/pnas.0710172105.

    PubMed Central  CAS  PubMed  Google Scholar 

  85. 85.

    Bandtlow CE, Zimmermann DR. Proteoglycans in the developing brain: new conceptual insights for old proteins. Physiol Rev. 2000;80(4):1267–90.

    CAS  PubMed  Google Scholar 

  86. 86.

    Yamaguchi Y. Lecticans: organizers of the brain extracellular matrix. Cell Mol Life Sci. 2000;57(2):276–89.

    CAS  PubMed  Google Scholar 

  87. 87.

    Kwok JC, Warren P, Fawcett JW. Chondroitin sulfate: a key molecule in the brain matrix. Int J Biochem Cell Biol. 2012;44(4):582–6. doi:10.1016/j.biocel.2012.01.004.

    CAS  PubMed  Google Scholar 

  88. 88.

    Giamanco KA, Matthews RT. Deconstructing the perineuronal net: cellular contributions and molecular composition of the neuronal extracellular matrix. Neuroscience. 2012;218:367–84. doi:10.1016/j.neuroscience.2012.05.055.

    PubMed Central  CAS  PubMed  Google Scholar 

  89. 89.

    Wang D, Fawcett J. The perineuronal net and the control of CNS plasticity. Cell Tissue Res. 2012;349(1):147–60. doi:10.1007/s00441-012-1375-y.

    PubMed  Google Scholar 

  90. 90.

    Bartus K, James ND, Bosch KD, Bradbury EJ. Chondroitin sulphate proteoglycans: key modulators of spinal cord and brain plasticity. Exp Neurol. 2012;235(1):5–17. doi:10.1016/j.expneurol.2011.08.008.

    CAS  PubMed  Google Scholar 

  91. 91.

    Beurdeley M, Spatazza J, Lee HH, Sugiyama S, Bernard C, Di Nardo AA, et al. Otx2 binding to perineuronal nets persistently regulates plasticity in the mature visual cortex. J Neurosci. 2012;32(27):9429–37. doi:10.1523/JNEUROSCI.0394-12.2012.

    PubMed Central  CAS  PubMed  Google Scholar 

  92. 92.

    McRae PA, Rocco MM, Kelly G, Brumberg JC, Matthews RT. Sensory deprivation alters aggrecan and perineuronal net expression in the mouse barrel cortex. J Neurosci. 2007;27(20):5405–13. doi:10.1523/JNEUROSCI.5425-06.2007.

    CAS  PubMed  Google Scholar 

  93. 93.

    Carulli D, Pizzorusso T, Kwok JC, Putignano E, Poli A, Forostyak S, et al. Animals lacking link protein have attenuated perineuronal nets and persistent plasticity. Brain. 2010;133(Pt 8):2331–47. doi:10.1093/brain/awq145.

    PubMed  Google Scholar 

  94. 94.

    Lee JS, Chien CB. When sugars guide axons: insights from heparan sulphate proteoglycan mutants. Nat Rev Genet. 2004;5(12):923–35. doi:10.1038/nrg1490.

    CAS  PubMed  Google Scholar 

  95. 95.

    Fox AN, Zinn K. The heparan sulfate proteoglycan syndecan is an in vivo ligand for the Drosophila LAR receptor tyrosine phosphatase. Curr Biol. 2005;15(19):1701–11. doi:10.1016/j.cub.2005.08.035.

    CAS  PubMed  Google Scholar 

  96. 96.

    Aricescu AR, McKinnell IW, Halfter W, Stoker AW. Heparan sulfate proteoglycans are ligands for receptor protein tyrosine phosphatase sigma. Mol Cell Biol. 2002;22(6):1881–92.

    PubMed Central  CAS  PubMed  Google Scholar 

  97. 97.

    Avram S, Shaposhnikov S, Buiu C, Mernea M. Chondroitin sulfate proteoglycans: structure-function relationship with implication in neural development and brain disorders. Biomed Res Int. 2014;2014:642798. doi:10.1155/2014/642798.

    PubMed Central  PubMed  Google Scholar 

  98. 98.

    Bruckner G, Seeger G, Brauer K, Hartig W, Kacza J, Bigl V. Cortical areas are revealed by distribution patterns of proteoglycan components and parvalbumin in the Mongolian gerbil and rat. Brain Res. 1994;658(1–2):67–86.

    CAS  PubMed  Google Scholar 

  99. 99.

    Pizzorusso T, Medini P, Berardi N, Chierzi S, Fawcett JW, Maffei L. Reactivation of ocular dominance plasticity in the adult visual cortex. Science. 2002;298(5596):1248–51. doi:10.1126/science.1072699.

    CAS  PubMed  Google Scholar 

  100. 100.

    Zhou XH, Brakebusch C, Matthies H, Oohashi T, Hirsch E, Moser M, et al. Neurocan is dispensable for brain development. Mol Cell Biol. 2001;21(17):5970–8.

    PubMed Central  CAS  PubMed  Google Scholar 

  101. 101.

    Meyer-Puttlitz B, Milev P, Junker E, Zimmer I, Margolis RU, Margolis RK. Chondroitin sulfate and chondroitin/keratan sulfate proteoglycans of nervous tissue: developmental changes of neurocan and phosphacan. J Neurochem. 1995;65(5):2327–37.

    CAS  PubMed  Google Scholar 

  102. 102.

    De Felipe J, Marco P, Fairen A, Jones EG. Inhibitory synaptogenesis in mouse somatosensory cortex. Cereb Cortex. 1997;7(7):619–34.

    PubMed  Google Scholar 

  103. 103.

    Micheva KD, Beaulieu C. Quantitative aspects of synaptogenesis in the rat barrel field cortex with special reference to GABA circuitry. J Comp Neurol. 1996;373(3):340–54. doi:10.1002/(SICI)1096-9861(19960923)373:3<340::AID-CNE3>3.0.CO;2[-‐]2.

    CAS  PubMed  Google Scholar 

  104. 104.

    Goldman JS, Ashour MA, Magdesian MH, Tritsch NX, Harris SN, Christofi N, et al. Netrin-1 promotes excitatory synaptogenesis between cortical neurons by initiating synapse assembly. J Neurosci. 2013;33(44):17278–89. doi:10.1523/JNEUROSCI.1085-13.2013.

    CAS  PubMed  Google Scholar 

  105. 105.

    Han L, Vickers KC, Samuels DC, Guo Y. Alternative applications for distinct RNA sequencing strategies. Brief Bioinform. 2014. doi:10.1093/bib/bbu032.

  106. 106.

    McGettigan PA. Transcriptomics in the RNA-seq era. Curr Opin Chem Biol. 2013;17(1):4–11. doi:10.1016/j.cbpa.2012.12.008.

    CAS  PubMed  Google Scholar 

  107. 107.

    Sengupta S, Bolin JM, Ruotti V, Nguyen BK, Thomson JA, Elwell AL, et al. Single read and paired end mRNA-Seq Illumina libraries from 10 nanograms total RNA. J Vis Exp. 2011;56:e3340. doi:10.3791/3340.

    PubMed  Google Scholar 

  108. 108.

    Pollock JD, Wu DY, Satterlee JS. Molecular neuroanatomy: a generation of progress. Trends Neurosci. 2014;37(2):106–23. doi:10.1016/j.tins.2013.11.001.

    PubMed Central  CAS  PubMed  Google Scholar 

  109. 109.

    Sun W, Zhang L, Lu J, Yang G, Laundrie E, Salvi R. Noise exposure-induced enhancement of auditory cortex response and changes in gene expression. Neuroscience. 2008;156(2):374–80. doi:10.1016/j.neuroscience.2008.07.040.

    PubMed Central  CAS  PubMed  Google Scholar 

  110. 110.

    Sharma V, Nag TC, Wadhwa S, Roy TS. Temporal distribution of mRNA expression levels of various genes in the developing human inferior colliculus. Neurosci Lett. 2009;461(3):229–34. doi:10.1016/j.neulet.2009.06.049.

    CAS  PubMed  Google Scholar 

  111. 111.

    Bush AL, Hyson RL. Effects of lithium and deafferentation on expression of glycogen synthase kinase-3beta, NFkappaB, beta-catenin and pCreb in the chick cochlear nucleus. Brain Res. 2008;1203:18–25. doi:10.1016/j.brainres.2008.01.076.

    PubMed Central  CAS  PubMed  Google Scholar 

  112. 112.

    Wang H, Brozoski TJ, Turner JG, Ling L, Parrish JL, Hughes LF, et al. Plasticity at glycinergic synapses in dorsal cochlear nucleus of rats with behavioral evidence of tinnitus. Neuroscience. 2009;164(2):747–59. doi:10.1016/j.neuroscience.2009.08.026.

    PubMed Central  CAS  PubMed  Google Scholar 

  113. 113.

    Wang H, Turner JG, Ling L, Parrish JL, Hughes LF, Caspary DM. Age-related changes in glycine receptor subunit composition and binding in dorsal cochlear nucleus. Neuroscience. 2009;160(1):227–39. doi:10.1016/j.neuroscience.2009.01.079.

    PubMed Central  CAS  PubMed  Google Scholar 

  114. 114.

    Mashiko H, Yoshida AC, Kikuchi SS, Niimi K, Takahashi E, Aruga J, et al. Comparative anatomy of marmoset and mouse cortex from genomic expression. J Neurosci. 2012;32(15):5039–53. doi:10.1523/JNEUROSCI.4788-11.2012.

    CAS  PubMed  Google Scholar 

  115. 115.

    Nehme B, Henry M, Mouginot D, Drolet G. The Expression Pattern of the Na(+) Sensor, Na(X) in the Hydromineral Homeostatic Network: A Comparative Study between the Rat and Mouse. Front Neuroanat. 2012;6:26. doi:10.3389/fnana.2012.00026.

    PubMed Central  CAS  PubMed  Google Scholar 

  116. 116.

    Van der Zee EA, Keijser JN. Localization of pre- and postsynaptic cholinergic markers in rodent forebrain: a brief history and comparison of rat and mouse. Behav Brain Res. 2011;221(2):356–66. doi:10.1016/j.bbr.2010.11.051.

    PubMed  Google Scholar 

  117. 117.

    Watakabe A, Komatsu Y, Sadakane O, Shimegi S, Takahata T, Higo N, et al. Enriched expression of serotonin 1B and 2A receptor genes in macaque visual cortex and their bidirectional modulatory effects on neuronal responses. Cereb Cortex. 2009;19(8):1915–28. doi:10.1093/cercor/bhn219.

    PubMed Central  PubMed  Google Scholar 

  118. 118.

    Zeng H, Shen EH, Hohmann JG, Oh SW, Bernard A, Royall JJ, et al. Large-scale cellular-resolution gene profiling in human neocortex reveals species-specific molecular signatures. Cell. 2012;149(2):483–96. doi:10.1016/j.cell.2012.02.052.

    PubMed Central  CAS  PubMed  Google Scholar 

  119. 119.

    Shukla R, Watakabe A, Yamamori T. mRNA expression profile of serotonin receptor subtypes and distribution of serotonergic terminations in marmoset brain. Front Neural Circuits. 2014;8:52. doi:10.3389/fncir.2014.00052.

    PubMed Central  PubMed  Google Scholar 

  120. 120.

    Lin S, Lin Y, Nery JR, Urich MA, Breschi A, Davis CA, et al. Comparison of the transcriptional landscapes between human and mouse tissues. Proc Natl Acad Sci U S A. 2014;111(48):17224–9. doi:10.1073/pnas.1413624111.

    PubMed Central  CAS  PubMed  Google Scholar 

  121. 121.

    Pongrac J, Middleton FA, Lewis DA, Levitt P, Mirnics K. Gene expression profiling with DNA microarrays: advancing our understanding of psychiatric disorders. Neurochem Res. 2002;27(10):1049–63.

    CAS  PubMed  Google Scholar 

  122. 122.

    Grange P, Bohland JW, Okaty BW, Sugino K, Bokil H, Nelson SB, et al. Cell-type-based model explaining coexpression patterns of genes in the brain. Proc Natl Acad Sci U S A. 2014;111(14):5397–402. doi:10.1073/pnas.1312098111.

    PubMed Central  CAS  PubMed  Google Scholar 

  123. 123.

    Hirokawa J, Watakabe A, Ohsawa S, Yamamori T. Analysis of area-specific expression patterns of RORbeta, ER81 and Nurr1 mRNAs in rat neocortex by double in situ hybridization and cortical box method. PLoS One. 2008;3(9):e3266. doi:10.1371/journal.pone.0003266.

    PubMed Central  PubMed  Google Scholar 

Download references

Acknowledgements

The authors gratefully acknowledge the support of NIH/NIDCD grants K18 DC012527 to T.A.H. and R01 DC009836 to D.P. Special thanks to Tia Hughes and Cara Sutcliffe in the VANTAGE core at Vanderbilt University for expert assistance with RNA isolation and sample quality assessment; Dr. Holli Hutcheson-Dilks in the VANTAGE core for design and supervision of RNA sequencing; the Vanderbilt Kennedy Center infrastructure grant P30-HD015052-33.

Author information

Affiliations

Authors

Corresponding author

Correspondence to Troy A. Hackett.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

TH, DP and KM conceived of and designed the study, and coordinated the workflow. AC, KG and TH harvested tissues and prepared samples for sequencing. YG and PZ conducted sequence alignment, quality assessment, differential analyses, and GSEA. TH and NH designed the database structure and analyzed the gene families datasets. All authors helped to draft, critique, edit, and approve the final manuscript.

Additional files

Additional file 1: Figure S1.

Sample harvesting. Low magnification coronal images at the level of A1 and MG. (A) Gapdh in situ hybridization; (B) photograph of a frozen brain during harvesting of samples from A1 and MG for sequencing. The location of A1 within the auditory cortex (AC) is shown, along with a sketch of the 0.5 mm punch used to obtain samples. Note that the size and shape of the punch compresses tissue outside of the punched volume. The left MG has been circumscribed prior to extraction. Scale bars, 1 mm all panels.

Additional file 2: Table S1.

Population and quality control data for all A1 and MG samples obtained for total RNA sequencing. For each sample, the following information is provided: region of interest (ROI), age, sex, RNA integrity number (RIN), 28 s:18 s ratio, 260/280 ratio. Samples are arranged sequentially by an identifier that is used in all other tables. Table S2. Raw data quality control matrix. Contains information about the sequencing platform, total reads, and sample quality measures. Table S3. Quality control of alignments. For each sample, information about the sequencing platform and mapping details are listed. Table S4. Raw counts of all genes for the 48 samples sequenced.

Additional file 3:

Differential expression comparing A1 and MG at at P7, P14, P21, and Adult. Table S5. Differential expression analyses comparing MG with A1 at P7 for all genes using DESeq2, EdgeR, and BaySeq. For each gene, the log2 fold-change, p-values, and rank are listed. The final column (All Rank) is the sum of rankings obtained by all three methods. Gene listing is given in. Table S6. Differential expression analyses comparing MG with A1 at P14 for all genes using DESeq2, EdgeR, and BaySeq. For each gene, the log2 fold-change, p-values, and rank are listed. The final column (All Rank) is the sum of rankings obtained by all three methods (lowest number = highest rank). Gene listing is ordered from highest to lowest rank. Table S7. Differential expression analyses comparing MG with A1 at P21 for all genes using DESeq2, EdgeR, and BaySeq. For each gene, the log2 fold-change, p-values, and rank are listed. The final column (All Rank) is the sum of rankings obtained by all three methods (lowest number = highest rank). Gene listing is ordered from highest to lowest rank. Table S8. Differential expression analyses comparing MG with A1 in adult animals for all genes using DESeq2, EdgeR, and BaySeq. For each gene, the log2 fold-change, p-values, and rank are listed. The final column (All Rank) is the sum of rankings obtained by all three methods (lowest number = highest rank). Gene listing is ordered from highest to lowest rank.

Additional file 4:

Differential expression in the MG for age group comparisons with P7 (P7-Adult, P7-P14, P7-P21). Table S9. Differential expression analyses in the MG comparing P7 with Adult using DESeq2, EdgeR, and BaySeq. For each gene, the log2 fold-change, p-values, and rank are listed. The final column (All Rank) is the sum of rankings obtained by all three methods (lowest number = highest rank). Gene listing is ordered from highest to lowest rank. Table S10. Differential expression analyses in the MG comparing P7 with P14 using DESeq2, EdgeR, and BaySeq. For each gene, the log2 fold-change, p-values, and rank are listed. The final column (All Rank) is the sum of rankings obtained by all three methods (lowest number = highest rank). Gene listing is ordered from highest to lowest rank. Table S11. Differential expression analyses in the MG comparing P7 with P21 using DESeq2, EdgeR, and BaySeq. For each gene, the log2 fold-change, p-values, and rank are listed. The final column (All Rank) is the sum of rankings obtained by all three methods (lowest number = highest rank). Gene listing is ordered from highest to lowest rank.

Additional file 5:

Differential expression in the MG for age group comparisons between P14, P21 and Adult (P14-P21, P14-Adult, P21-Adult). Table S12. Differential expression analyses in the MG comparing P14 with P21 using DESeq2, EdgeR, and BaySeq. For each gene, the log2 fold-change, p-values, and rank are listed. The final column (All Rank) is the sum of rankings obtained by all three methods (lowest number = highest rank). Gene listing is ordered from highest to lowest rank. Table S13. Differential expression analyses in the MG comparing P14 with Adult using DESeq2, EdgeR, and BaySeq. For each gene, the log2 fold-change, p-values, and rank are listed. The final column (All Rank) is the sum of rankings obtained by all three methods (lowest number = highest rank). Gene listing is ordered from highest to lowest rank. Table S14. Differential expression analyses in the MG comparing P21 with Adult using DESeq2, EdgeR, and BaySeq. For each gene, the log2 fold-change, p-values, and rank are listed. The final column (All Rank) is the sum of rankings obtained by all three methods (lowest number = highest rank). Gene listing is ordered from highest to lowest rank.

Additional file 6:

Differential expression in A1 for age group comparisons with P7 (P7-Adult, P7-P14, P7-P21). Table S15. Differential expression analyses in A1 comparing P7 with Adult using DESeq2, EdgeR, and BaySeq. For each gene, the log2 fold-change, p-values, and rank are listed. The final column (All Rank) is the sum of rankings obtained by all three methods (lowest number = highest rank). Gene listing is ordered from highest to lowest rank. Table S16. Differential expression analyses in A1 comparing P7 with P14 using DESeq2, EdgeR, and BaySeq. For each gene, the log2 fold-change, p-values, and rank are listed. The final column (All Rank) is the sum of rankings obtained by all three methods (lowest number = highest rank). Gene listing is ordered from highest to lowest rank. Table S17. Differential expression analyses in A1 comparing P7 with P21 using DESeq2, EdgeR, and BaySeq. For each gene, the log2 fold-change, p-values, and rank are listed. The final column (All Rank) is the sum of rankings obtained by all three methods (lowest number = highest rank). Gene listing is ordered from highest to lowest rank.

Additional file 7:

Differential expression in A1 for age group comparisons between P14, P21 and Adult (P14-P21, P14-Adult, P21-Adult). Table S18. Differential expression analyses in A1 comparing P14 with P21 using DESeq2, EdgeR, and BaySeq. For each gene, the log2 fold-change, p-values, and rank are listed. The final column (All Rank) is the sum of rankings obtained by all three methods (lowest number = highest rank). Gene listing is ordered from highest to lowest rank. Table S19. Differential expression analyses in A1 comparing P14 with Adult using DESeq2, EdgeR, and BaySeq. For each gene, the log2 fold-change, p-values, and rank are listed. The final column (All Rank) is the sum of rankings obtained by all three methods (lowest number = highest rank). Gene listing is ordered from highest to lowest rank. Table S20. Differential expression analyses in A1 comparing P21 with Adult using DESeq2, EdgeR, and BaySeq. For each gene, the log2 fold-change, p-values, and rank are listed. The final column (All Rank) is the sum of rankings obtained by all three methods (lowest number = highest rank). Gene listing is ordered from highest to lowest rank.

Additional file 8: Figure S2.

Comparison of sequencing with in situ hybridization. (Top) In situ hybridization of Slc32a1, Slc17a6, Slc17a7 from P7 to adult in A1 and MG. Coronal sections. Roman numerals indicate cortical layers. Abbrevations: d, dorsal division of MG; v, ventral division of MG; sp, subplate layer. Scale bars, 250 μm all panels. (Bottom) Normalized counts from RNAseq are compared to expression levels derived from quantitative densitometry of colorimetric in situ hybridization (ISH) assays performed at each maturational age (P7, P14, P21, Adult). Results are plotted separately for MG and A1. The housekeeping gene, Gapdh, had a flat maturational trajectory for both RNAseq and ISH, and was used for normalization of the ISH grayscale intensity measurements.

Additional file 9: Table S21.

Read counts for genes in all HUGO gene families. Normalized read counts for 19,826 genes in all HUGO gene families currently listed in that resource, listed alphabetically by gene family. Source: HUGO Gene Families database, HUGO Gene Nomenclature Committee at the European Bioinformatics Institute (http://www.genenames.org). Table S22. Gene set enrichment analysis of 111 HUGO gene families. Genesets are listed by functional group and family (from Table S21). For each geneset, the number of genes (Size), false discovery rate (FDR q-value), and direction of expression (up, down) are listed. Shading denotes gene families enriched above the 25 % FDR cutoff. Table S23. Normalized read counts for genes in selected gene ontology (GO) categories from the MGI database. Genes are arranged alphabetically by GO category, along with gene symbols, names, and annotations. Source: Mouse Genome Informatics (MGI) Gene Ontology Browser, Jackson Laboratories (www.informatics.jax.org). Table S24. Gene set enrichment analysis of the 51 gene ontology categories. Genesets are listed by functional group and GO designation (from Table S23). For each geneset, the number of genes (Size), false discovery rate (FDR q-value), and direction of expression (up, down) are listed. Shading denotes categories enriched above 25 % FDR cutoff.

Additional file 10: Table S25.

Gene families database. Normalized read counts and plots of ~4700 genes in 237 gene families from the HUGO database. Gene families are grouped by 20 functional categories arranged in tabs. For each family, normalized read counts (mean and % of maximum) are provided for A1 and MG as a function of postnatal age to permit rapid screening of the relationships between genes by brain region and age.

Additional file 11:

Look-up and plotting tool. This tool supports automated extraction of normalized counts for A1 and MG, and generates a plot of the counts and correlation matrix comparing maturational trajectories between each age group. A list of up to 25 genes can be entered into Column A. Use standard gene names only. Not case sensitive, but omit empty spaces. The normalized counts will populate in columns B – AU, and the data plotted below (profile plots, correlation matrices). Erase or overwrite gene names to generate a new listing. Counts are extracted from the Normalized Reads tab. This tab is locked and should not be edited.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Hackett, T.A., Guo, Y., Clause, A. et al. Transcriptional maturation of the mouse auditory forebrain. BMC Genomics 16, 606 (2015). https://doi.org/10.1186/s12864-015-1709-8

Download citation

Keywords

  • Synapse
  • Plasticity
  • Development
  • Critical period
  • Cortex
  • Thalamus
  • Neurotransmission
  • Neuromodulation
  • Extracellular matrix
  • Myelination
  • RNAseq
  • Pathway analysis
  • Sequencing
  • RNA