A total of 56 BXD strains (one male and one female for most strains) and 54 BXD strains (two male and two female for most strains) plus their parental strains are used for collecting hippocampus gene expression data and hearing screening, respectively. The age of all mice except one BXD101 mouse was between 12 and 32 months. The mice were housed in groups in a temperature- and humidity-controlled vivarium with a constant 12-h light–dark cycle with ad libitum access to food and water. For the hippocampal profiling, the mice were anesthetized by cervical dislocation. For the hearing screening, the mice were anesthetized with an intraperitoneal injection (IP) of ketamine, xylazine, and acepromazine at doses of 40, 5, and 1 mg/kg, respectively. The present study was carried out in accordance with the Guidelines for the Care and Use of Laboratory Animals published by the National Institutes of Health and was approved by the Animal Care and Use Committee at the University of Tennessee Health Science Center (UTHSC, Memphis, TN, United States).

Hearing acuity was assessed using an auditory-evoked brainstem response (ABR) (Intelligent Hearing Systems, Miami, FL, United States) as detailed previously (Zhou et al., 2006). Briefly, the mice were anesthetized with an intraperitoneal injection (IP) of ketamine, xylazine, and acepromazine at doses of 40, 5, and 1 mg/kg, respectively. The body temperature was maintained at 37–38°C. The ABRs were recorded using platinum subdermal needle electrodes inserted at the vertex (active electrode), ventrolateral to the right (reference electrode) and left (ground electrode) ears. The acoustic stimuli were tone-bursts (3-ms duration with a 1.5-ms cosine-gated rise/fall time) that were delivered through a high-frequency transducer. The stimuli were presented in a 5-dB step decrement from 80 dB sound pressure level (SPL) until the lowest intensity that could still evoked a reproducible ABR pattern was detected. If 80 dB could not evoke a reproducible ABR pattern, the stimuli were increased in a 5-dB step until the maximal SPL was reached. All the animals were tested with three frequencies (8, 16, and 32 kHz). For each frequency, the strain ABR is the mean value of individual ABRs, and the ABR thresholds from three frequencies were averaged and used as hearing acuity feature.

Snap frozen hippocampi from BXD mice across 56 strains were used for RNA quantification. RNA was extracted using the RNeasy mini kit (Qiagen, CA, United States) according to the instructions of the manufacturer. The 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, United States) was used to evaluate RNA integrity and quality. Samples with RNA Integrity Numbers (RIN values) > 8.0 were analyzed on Affy MoGene1.0 ST at the UTHSC. Raw microarray data were normalized using the Robust Multichip Array (RMA) method (Chesler et al., 2005; Geisert et al., 2009; King et al., 2015; Lu et al., 2016). The expression data were then re-normalized using a modified Z score described in a previous publication (Chesler et al., 2005). Briefly, RMAs were first transformed into log2-values. Then, the data of each single array was converted to Z-scores, multiplied by 2, and a value of 8 was added. The normalized data is available on GeneNetwork¹ under the “BXD” group and “Hippocampus mRNA” type with the identifier “UTHSC BXD Aged Hippocampus Affy MoGene1.0 ST (May 15) RMA Gene Level.”

The published learning-related traits of the BXD mice used in this study were retrieved from our GeneNetwork. The detailed descriptions can be found in the previous publications (Graybeal et al., 2014; Delprato et al., 2015; Knoll et al., 2016; Neuner et al., 2019)]. The summary statistics and individual values are available on GeneNetwork under the “BXD” group, “trait and cofactors” type, and “BXD published phenotypes” dataset with their corresponding GN accession number listed in Supplementary Table 1.

WEB-based Gene SeT AnaLysis Toolkit (WebGestalt) was used to perform gene set enrichment analysis with the mouse genome reference gene set as the background (Liao et al., 2019)². The over-representation of Gene Ontology (GO) was determined by the hypergeometric test. KEGG Orthology-Based Annotation System (KOBAS) was used to analyze the pathways involved in the ARHL correlated genes (Bu et al., 2021).

Gene co-expression network analysis has been widely used to explore the key genes in the gene sets. Therefore, we deployed this approach to identify the key genes from the gene set that correlated with hearing loss (Zhang and Horvath, 2005). Briefly, a network was constructed with Pearson’s correlation coefficient matrix (Supplementary Figure 1). In the network, each node stands for a gene, and the correlation coefficient was set as the edge. Binomial correlation higher than 0.3 or lower than −0.3 was defined as connected. The connection weight was calculated for each node, which is the sum of the binominal correlation coefficient connected to each node. The gene with the most connectivity and connection weight in the network can be considered as the central hub gene.

The gene–phenotype and gene–gene correlations were performed on the GeneNetwork (see text footnote 1) online platform by using Pearson’s correlation. A p-value lower than 0.05 was considered as statistically significant. A correlation coefficient higher than 0.3 or lower than −0.3 was considered as moderate correlation. To further validate the false discovery rate (FDR) of the genes-of-interest, we performed a further permutation test based on the Westfall and Young’s multiple testing procedure (Chaubey, 1993), Briefly, we randomly permuted the ABR measurements 1,000 times. For each permutation, we computed the p-value of Pearson’s correlation between the randomized ABR value and those transcripts that are associated with ABR. The adjusted p-value was determined by ranking the correlation coefficient.

A total of 26 strains that overlapped between hearing and hippocampus transcriptomic data were used for association analysis. Those strains are mainly in the similar age (from 19 to 25 months), except for BXD101 (11 months), BXD55 (27 months), and BXD45 (30 months) (Supplementary Table 2, Supplementary Figure 2, and Figure 1A). The tested animals represent a diverse degree of hearing loss determined with ABR thresholds between 33 (BXD74) and 100 dB (BXD101). The mean ABR thresholds was 71 dB SPL (Figure 1B), and the median ABR thresholds was 76 dB SPL. These data were used for the association analysis with hippocampus transcriptomic profile. In total, 1,435 genes presented significant correlation with hearing acuity (p < 0.05, Pearson’s correlation). With the FDR threshold of 0.05 determined by the Westfall and Young’s multiple testing procedure, equivalent to the r-value of 0.388, all the 1,435 transcripts were achieved significance. GO analysis of this gene set showed that four of the top 10 categories were associated with synapse organization (Figure 1C), and these genes were abundantly enriched in cellular compartment of synapse (ratio = 0.07, p = 6.17E−13) (Figure 1D). Further pathway analysis indicated that the glutamatergic synapse pathway is significantly enriched (ratio = 0.08, p = 0.0038) (Figure 1E). Nine genes that correlated with ABR thresholds were involved in this pathway, including four genes that showed positive correlation (Adcy4, Mapk3, Shank3, and Dlg4; n = 26, r > 0.4, p < 0.05, Pearson’s correlation), and five genes that showed negative correlation with ABR thresholds (Adrbk2, Slc38a1, Slc38a2, Gria3, and Gls, n = 26, r < −0.4, p < 0.05, Pearson’s correlation) (Figure 1F). To better clarify the gene expression correlation to each frequency, we supplemented the gene expression of synapse pathway to each of three frequencies (Supplementary Figure 3), which showed that all of those genes were significantly associated with the ABR thresholds at all three frequencies.

The mouse strains in our study are a senile population. Even though most of the strains are in similar age, there is still age variance among the strains. To exclude potential influence of such variance on the gene expression, we performed linear regression analysis between age and various target transcripts expression. The results indicated that the transcript expression variance is dominated by strain background but not age variance. Similarly, the age variation in our senile population has little effect on ABR threshold (Supplementary Tables 3, 4).

Glutamate receptors are the primary mediators of excitatory transmission in the central nervous system and play a pivotal role in learning and memory. To further investigate the glutamate receptor profiling associated with ARHL, we performed a correlation screening of ABR thresholds to 32 glutamate receptors (Figure 2A and Supplementary Table 5). Besides Gria3, which had the most significant correlation (n = 26, r = −0.41, p = 0.0401) (Figure 2B), six other receptors also showed moderate correlation to ABR thresholds (Figures 2C–H, r < −0.3), including Grm7, Grik1, Gria4, Grm5, Grm8, and Gria2. Notably, all these were negative correlations to ABR thresholds. These results indicated that ARHL is associated with altered glutamate receptor expression profiling in the hippocampal synapse.

Two genes involved in glutamate synthesis are significantly correlated with ABR thresholds, which are Gls (n = 26, r = −0.46, p = 0.0193) and Slc38a2 (n = 26, r = −0.49, p = 0.0109) (Figures 3A,B). Both genes presented negative correlation to ABR thresholds. In the glutamatergic synapse pathway, the glutamine is primarily released from astroglia cells and further imported into presynaptic neurons with glutamine receptor Slc38a2. The imported glutamine is further converted to glutamate with Gls in the presynaptic neurons. The negative correlation of these two genes with ABR thresholds indicated that glutamate synthesis is impaired in ARHL.

Two glutamate receptor-interacting protein chaperones—Dlg4 and Shank3—were positively correlated with hearing loss (n = 26, r = 0.39, p = 0.0493); (n = 26, r = 0.44, p = 0.0247) (Figures 3C,D). Both Dlg4 and Shank3 act as chaperones to assist glutamate-receptor reorganization. Thus, we performed a further correlation screening between hearing loss and 10 known glutamate receptor chaperone proteins (Figure 3E). Besides Dlg4 and Shank3, Nos1 (n = 26, r = 0.45, p = 0.0198) (Figure 3F) also presented significant positive correlations with ABR thresholds. These data showed that ARHL is associated with enhanced glutamate receptor chaperone expression.

The cyclic adenosine monophosphate (cAMP) signaling is a critical second messenger signaling in the glutamatergic synapse pathway. Our correlation analysis indicated that the gene expression of Adcy4 is significantly correlated with ABR thresholds (n = 26, r = 0.52, p = 0.0064) (Figure 4A). cAMP is synthesized by adenylate cyclase (Figure 4B), which is a protein family including 10 different members. Thus, a correlation screen was performed between ABR thresholds and AC family members (Figure 4C). Among the AC family, only Adcy4 presented a significant correlation with ABR thresholds, suggesting that ARHL is associated with cAMP signaling pathway through Adcy4.

A gene co-expression network was constructed based on nine hearing loss-associated genes in the synaptic signaling pathway. Among these genes, Gls showed the highest connectivity of eight and an average connection weight of 0.49 as a central node (Figure 5A and Supplementary Tables 6, 7). Notably, Gls presented a significant correlation with five other hub genes; Gria3 (n = 26, r = 0.76, p < 0.0001), Adcy4 (n = 26, r = −0.55, p < 0.0001), Dlg4 (n = 26, r = −0.35, p = 0.008), and Shank3 (n = 26, r = −0.56, p < 0.0001) (Figures 5B–E). Taken together, these results suggest that Gls is a central hub gene in the glutamate signaling network.

To reveal the association between Gls expression and learning behavior, we performed a phenotype-wide association analysis between Gls and published phenotypes in BXD strains. The expression of Gls in the hippocampus is significantly correlated with 31 learning-related phenotypes (Supplementary Table 1). Briefly, the expression of Gls was significantly correlated to the learning latency in a touch screen test (n = 19, r = 0.66, p = 0.0020), the percentage of time spent freezing in contextual fear conditioning test (n = 16, r = 0.52, p = 0.0418), the percentage of successful alternations in the Y-maze test (n = 22, r = 0.62, p = 0.0021), and the number of new entries during the first 8-arm choices in an 8-arm radial maze test (n = 41, r = 0.38, p = 0.0151) (Figures 6A–D). These data collectively proved that the hippocampal Gls expression is associated with learning and memory behavior.