Tissue source and age are associated with differences in transcription
Sensory epithelia were dissected from the cochlea and vestibule of mice at two stages of development, embryonic day 16.5 (E16.5) and postnatal day 0 (P0). This was followed by RNA-seq, as previously described [24, 25]. Our analysis identified 39,178 Ensembl genes (including non-coding genes and pseudogenes), 15,206 of which have at least one read per million in three or more of the samples. A principal component analysis (PCA) plot demonstrated four well separated groups (Fig. 1). The first principal component (PC1) explained almost half the variance and is associated with the age of the sample, whereas PC2 explained about a quarter of the variance and is associated with the originating tissue (F test on associations, p values = 1.99 × 10−5, 1.31 × 10−5, respectively). Additional PCs were not associated with either tissue or age (p value ≥ 0.05). The E16.5 genes displayed a lower intra-group variability than at P0. This might reflect differences in the rate of development of the different organs between mice from the same population in the period between E16.5 and P0.
We used linear mixed models to estimate the percentage of variance that can be attributed to age, tissue, or the interaction of age and tissue (Additional file 1: Supplementary Methods). According to our estimates, the majority of variance can be attributed to either age (44.0 ± 6.5) or tissue (39.6 ± 5.4) (mean percentage ± standard deviation). The remaining non-negligible percentage can be attributed to the interaction term (8.0 ± 1.5), and a model with this interaction term describes the data better according to a restricted likelihood ratio test (p value ≤ 2.2 × 10−16). Less than 10% of the variance was left unexplained (8.4 ± 1.08).
We selected genes that were differentially expressed between tissues or between ages, and genes for which the interaction of tissue and age was significant in determining expression. Our results identified 3306 upregulated genes and 6890 downregulated genes at P0 compared to E16.5. Four thousand one hundred fifty-nine genes were found to be upregulated and 2382 were downregulated in the vestibule compared to the cochlea. For 745 genes, the cochlea to vestibule expression ratio increased over development, and it decreased for 1211 genes. We performed gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses on genes from the six identified sets (Additional file 2: Table S1). The enrichment results are summarized below.
Expression changes with age
Genes that were upregulated at E16.5 are enriched for terms related to cell cycle, DNA replication, cytoskeleton organization, and other terms that are in accordance with a highly proliferative state. In contrast, genes that were upregulated at P0 are enriched for ribosomes, indicating high protein synthesis, mainly of plasma membrane and extracellular matrix proteins. The lipid and oxphos-related metabolic activities are also high in this group. The cells at this stage of development are more adhesive, communicate more with one another, and are more responsive to external cues. They are also responsive to a variety of signaling receptors, including calcium signaling, and have high ion transport activity. The upregulated terms are typical of a less proliferative environment, where the highly expressed genes promote homeostatic processes and inhibit peptidase activity. Some terms show signs of cell specialization, in terms of sensory perception, cartilage-related metabolism, and the regulation of ossification; the last might indicate a cross-talk between sensory epithelium cells and endochondral cells. Another marker for the more differentiated state is an up-regulation of the MHC protein complex. In summary, the enrichment suggests that the inner ear is in a more proliferative state at E16.5 than at P0, whereas at P0 the tissues are more differentiated and exhibit specialization for sensory perception.
KEGG enrichment generally confirmed the aforementioned differences and provided more details regarding specific metabolic processes activated at P0. For example, we could attribute the enriched lipid metabolism to sphingolipids, arachidonic acid, and retinol, the enriched aminoglycan metabolism to glycan degradation, and the biosynthesis of chondroitin and keratan sulfate. Pathways enriched at P0 suggest that the activity of the immune system increases during development, with leukocytes migrating into the tissue and intercellular communication using cytokines. As the complement and coagulation cascades and the renin-angiotensin system are also enriched at P0, we can hypothesize that the inner ear is more exposed to blood circulation at this age.
Expression change between tissues
According to the enrichment analysis (Additional file 2: Table S1), a number of the differentially expressed (DE) genes in both the cochlea and vestibule are involved in signal transduction. In the cochlea, the majority of the signaling is mediated by voltage- and ligand-gated ion channels and can be attributed to neuron-neuron synaptic transmission. In agreement with this finding, other upregulated activities are neurogenesis and neuron projection. In contrast, the signaling in the vestibule is probably required for the coordination of both innate and acquired immune responses, an observation that relates to the main function enriched in this tissue. The signaling, some of which involves purinergic receptors, plays a role in the response to external stimulus and stress, and also in taxis. Another function enriched in the vestibule is locomotion, with the cilium and the axoneme being two enriched cellular components related to the movement of the hair cells’ stereocilia. The vestibule is richer in blood vessel formation and hematopoiesis, and the extracellular matrix is more evolved than in the cochlea. Together with the high immune-related activity, these factors may explain why the vestibular cells are more adhesive. We also detected enrichment for replacement ossification, suggesting the development of bone. As a generalization, upregulated genes were associated with neurological terms in the cochlear, but to vascular, structural, and immunological terms in the vestibule. This partitioning was not perfect as we could detect enrichment for mesenchymal cell differentiation in the cochlea, and 3.1% of the upregulated genes in the vestibule were annotated for a role in sensory perception.
The KEGG enrichment data also agreed with the characterization of the cochlea as more neurological versus a more vascular vestibule. In addition, the data provided more information about the typical signaling in each apparatus. Neuroactive ligand signaling was identified in both, although the cochlea was associated with the TGF-beta, MAPK, and ErbB signaling pathways, while cytokine-mediated, calcium, and Toll-like receptor signaling were more important in the vestibule. Three pathways shown to be unique to the cochlea affect cell proliferation, survival, differentiation, and migration [26,27,28], suggesting that these developmental processes are more activated in the cochlea. Other unique metabolic pathways enriched in the cochlea were O-glycan and chondroitin sulfate biosynthesis. The vestibule, on the other hand, was enriched for glycan degradation and metabolic pathways concerning arachidonic acid, retinol, and glutathione.
Tissue expression ratio change with age
Genes for which the cochlea to vestibule expression ratio increased with age \( \left(\frac{Cochlea}{Vestibule}\uparrow \right) \) were enriched for processes related to sensory perception and central nervous system development, as well as signaling through G-coupled receptors, ligand-gated ion channels, or calcium. Accordingly, a significant number of genes were annotated to be in the apical part of the cell. Other genes annotated to the extracellular region might mediate the biological adhesion, which increases during development. Another enriched component was identified as the sarcomere, which most closely resembles the stereocilia in the inner ear.
We can envision two possible scenarios for each of these enrichments. The first option is that genes annotated for enrichment are upregulated in the cochlea at E16.5 and the gap between the cochlea and the vestibule increases during development. The second option is that these genes are upregulated in the vestibule at E16.5 and the gap between the cochlea and the vestibule decreases during development. To distinguish between the two, we compared the expression of all genes that are annotated for each GO term. The median expression log-ratio between the cochlea and the vestibule at P0 was plotted against the value of the same parameter at E16.5 (Fig. 2, circles). The plot only contains the terms for which the gap between the cochlea and the vestibule significantly increases with age. More precisely, only terms for which the log-ratios at P0 were larger than their paired values at E16.5 were included (Wilcoxon signed rank test, q values ≤ 0.05).
Interestingly, the vestibule is appeared to be more specialized for sensory perception at E16.5 than the cochlea, as manifested by a negative median log-ratio for terms sensory perception, mechanoreceptor differentiation, and detection of stimulus involved in sensory perception. However, by P0, the cochlea surpassed the vestibule in all of these fields. In contrast, ligand-gated ion channel activity was already higher in the cochlea at E16.5, and the gap only increased with development.
Genes for which the vestibule to cochlea ratio increased with age \( \left(\frac{Vestibule}{Cochlea}\uparrow \right) \) were enriched for signaling, neuron projection, neurotransmitter transport, and secretion. These are all functions that are higher in the cochlea at E16.5, and for which the difference between the vestibule and the cochlea decreases with time (Fig. 2, triangles).
Deafness genes can be predicted using expression patterns
A list of 140 genes associated with human deafness was compiled from a public dataset (http://hereditaryhearingloss.org/; Additional file 3: Table S2). Expression data for 130 orthologous mouse genes are available in our dataset. Of these genes, mutations in 25 orthologs are associated with syndromic deafness in human, 96 with non-syndromic deafness, and nine with both types of deafness and are treated as syndromic deafness genes (DGs) in subsequent analyses. It should be noted that we found no ortholog for any of the five mitochondrial DGs.
We observed general patterns of expression for these syndromic and non-syndromic DGs. First, when comparing vestibular and cochlear expression, the absolute values of the fold change (FC) of the DGs were higher than for the background FCs (p value = 1.98 × 10−5, one-sided Wilcoxon rank sum test; Fig. 3, upper subfigure). In addition, the absolute FCs of non-syndromic DGs were slightly higher than the FCs of syndromic DGs (p value = 7.00 × 10−2, same test). That is, DGs tend to be tissue-specific, with the non-syndromic genes possibly being even more specific. Interestingly, the majority of the DE DGs were higher in the vestibule than in the cochlea, despite the acknowledged role of the cochlea in hearing (57 out of 76, p value = 2.19 × 10−5, two-sided proportion test).
Second, when comparing P0 and E16.5 expression, DGs tended to have higher FCs compared to background FCs (p values = 6.32 × 10−6, one-sided Wilcoxon rank sum test; Fig. 3, middle subfigure). This indicates that their expression tends to increase with development. Third, their cochlea to vestibule expression ratio tended to increase with development compared to background (p values = 5.15 × 10−6, same test; Fig. 3, lower subfigure). Moreover, the increase in the ratio of non-syndromic DGs was higher than that for syndromic genes (p value = 3.48 × 10−3, same test).
Deafness gene prediction
We used the three types of FC and the averaged expression (see the “Methods” section), to build a classifier that can predict whether a gene is a DG. The classifier achieved a ROC score of 0.66 ± 0.04 across repeated training/test splits. A ROC score of greater than 0.5 indicates that the expression data have some predictive value for the relation of a gene to deafness. This classifier performs better than a similar version that used the averaged RPKM values in each condition (ROC score 0.60 ± 0.05). Removing one or more of the four feature types from the original classifier resulted in a lower score.
It must be appreciated that genes not marked as DGs might still represent undiscovered DGs. For this reason, it was important to train our classifier to distinguish between known DGs and genes with an unknown role in deafness. The first group of genes is termed positive, and those in the second group are classified as unlabeled. We wished to adapt our positive unlabeled (PU) classifier to output the probability that an unlabeled gene is a positive gene. This type of classification is referred to as transductive PU learning [29]. Supposing that the known DGs are a random subset of all DGs, i.e., the features we explore impose no bias over which of the positive genes are labeled, then, the probability that the PU classifier assigns to the positivity of new genes both (i) correctly ranks the genes and (ii) the probabilities are only off by a constant factor (see [30] for details). We used a bagging-like algorithm similar to the one presented in [29] in order to calculate the probabilities for the set of unlabeled genes. Some modifications in our system are described in the “Methods” section. One main difference between our approach and the previously reported version in [29] was that we kept the same proportion of positive (labeled) samples in the training set as in the test set, whereas in [29], all positive samples were included in training. This property allowed us to address the issue of biases in the probabilities, albeit at the price of losing some predictive power. One source of bias was due to undersampling in the learning process [31]. A second source of bias was the one described above for a PU classifier. We addressed the latter using methods presented in [30].
To gain some insight about the accuracy of our estimator, in spite of the lack of a definitive classification of the unlabeled set, we downloaded lists of genes associated with hearing loss according to the text mining tools DigSeE [32], DisGeNET [33], and DISEASES [34]. We refer to these genes as deafness-associated genes (DAGs). By these means, we obtained 1313 genes that were associated with deafness according to at least one tool. These included 115 known DGs, accounting for 82% of all reported DGs. The respective numbers of mouse orthologs were 1021, 106, and 82%. See Additional file 1: Supplementary Results, Figure S1 for a comparison of the lists of genes provided by the tools.
Applying our bagging-like algorithm resulted in a PU classifier with a ROC score of 0.694 where the probabilities from this native classifier were probably biased upward due to undersampling. Correcting for this bias resulted in a better calibration of the probabilities, as demonstrated by a calibration plot (Additional file 1: Figure S2, left), and the lowering of the Brier score (BS) from 2.07 × 10−1 to 8.47 × 10−3. We employed three different methods (e1, e2, e3; see [30]) to correct the bias in the probabilities caused by the PU scenario. In order to perform the calibration, all three methods first estimate the probability that a known DG is labeled p(s = 1| y = 1). The estimates for this probability, according to e1, e2, and e3, were 0.032 ± 0.014, 0.022 ± 0.007, and 0.518 ± 0.248, respectively. The estimates made by e1 and e2 support the existence of a few thousand DGs, compared to the few hundred predicted according to e3 (4.1 × 103, 5.9 × 103, and 2.5 × 102, respectively). We believe that given the status of deafness research, the last estimate is the most reasonable. To investigate the issue further, we re-evaluated the calibration of the probabilities produced by each method. For this purpose, we assumed that all the DAGs are in fact deafness genes. With this assumption, the e3 method resulted in the best calibration, as demonstrated by a calibration plot (Additional file 1: Figure S2, right), and the lowest BS (scores 6.64 × 10−2, 1.20 × 10−1, 2.84 × 10−1, 6.45 × 10−2 for no fix, e1, e2, and e3, respectively). Hence, we decided to use e3 probabilities in all subsequent analyses and let p
g
be the probability that gene g is positive according to e3.
We then reran our bagging-like algorithm, but this time, we chose to treat a gene g as positive with probability p
g
, and as negative with probability 1 − p
g
. This reassignment was performed before each iteration. Finally, we recalculated the ROC score of our classifier. In this case, we ignored known DGs in order to make a proper separation between training and test stages. The rerun achieved a slightly better ROC score (0.602 vs 0.600, p < 0.05, DeLong’s test for two correlated ROC curves [35]). We chose to continue with the rerun classifier and added a correction for undersampling to the resultant probabilities. The predictions for both human genes and mouse orthologs are available in Additional file 4: Table S3. The 20 mouse genes with the highest predicted probabilities include the known non-syndromic DGs Smpx, and Ptprq, seven DAGs (Gfi1, Lhx3, Erbb4, Ephx1, Il33, Slc52a3, and Ttr), and nine genes not associated with deafness (Mlf1, Nell1, Espnl, Rbm24, Lrrc10b, Agr3, Tgm2, Id4Cd164l2, and Faim2).
For the purpose of selecting a discrimination threshold for our binary classifier, we can consider two plots, which demonstrate how well our classifier predicts DAGs (again while ignoring known DGs). The first is a ROC curve, which visualizes the balance between specificity and sensitivity (Fig. 4, top). The threshold maximizing the sum of these two parameters is suggested as a candidate threshold. A disadvantage of a ROC curve, in our context, is that it ignores the association scores provided by the text mining tools. In order to account for these scores, we can consider a range of values of the threshold and use a non-parametric test (one-tailed Wilcoxon rank sum test) to compare the association scores of the genes with probabilities higher than the threshold, with all the others. We hypothesized that genes above the “right” threshold would tend to have higher association scores. We analyzed the association scores from each tool separately and together (see the “Methods” section) and plotted −log2P − value against the threshold (Fig. 4, bottom). The value giving the lowest p value for the combined scoring was proposed as a candidate threshold. Four thousand six hundred seventy-four and 1934 genes passed the thresholds suggested by the ROC curve (0.027) and the Wilcoxon test (0.043), respectively. Other thresholds may also be considered, depending on the required number of candidates, specificity, and sensitivity. We recommend choosing thresholds that give local maxima on either curve (available in Additional file 4: Table S3).
Transcription factors affecting expression
When we screened for enrichment of transcription factor (TF) binding sites in three sets of DE genes (Additional file 5: Table S4) we could identify six motifs that were associated with changes in expression during development, 43, between tissues, and 10, across an age-tissue interaction (i.e., the change in the cochlea to vestibule expression ratio throughout development). This 7-fold increase in tissue-specific motifs over those associated with a developmental stage was very surprising, in view of the fact that the absolute number of tissue-specific DE genes identified was about 35% less than the number that changed during development. In total, we identified 50 unique motifs across all comparisons and manually connected them to 64 mouse TFs (i.e., a few motifs were associated with multiple TFs).
For each TF, we tested whether the TF gene itself was DE under the same conditions as the gene it regulates (Additional file 5: Table S4). This property interests us for three reasons: (i) It indicates whether the regulation of the TF activity is (at least partially) transcriptional. Knowing how a TF is regulated makes it a better candidate for experimental interventions. (ii) The direction (upregulation or downregulation) in which a TF is DE implies whether it functions as a repressor or an activator. (iii) It strengthens our faith that the associated motif is important for regulation, and not a false-positive. In our analysis, 30 of the 64 TFs identified were DE (in at least one comparison).
In order to investigate how the levels of the TFs affect their targets, we plotted the median FC of all targets of a specific motif, against the median FC of the TFs associated with that motif (Additional file 1: Figure S3). In all cases, we observed a positive, although insignificant, correlation between the two values (Pearson’s r = 0.51, 0.05, or 0.52 for the comparisons between tissues, across age, and for age-tissue interaction, respectively; combined p value (22) = 0.15) [36]. Among the factors that contribute to the incomplete correlation is the post-transcriptional regulation of TFs, which reduces the correlation between the transcript levels and TF activity. In addition, while most TFs activate the transcription of the targets, certain factors can repress the transcription of some or all of their targets. Moreover, taking the median FC of the TFs associated with a motif ignores the possibility of complex interrelationships, such as the ability of a subset of the TFs to activate transcription alone (an example for the motif AHRHIF is discussed later).
The TFs identified as being associated with development were compared with TF genes shown to change their expression during avian regeneration of inner ear sensory epithelia in one of two experiments conducted in chick. The first experiment measured expression of TFs after either laser “wounding” cultured sensory epithelia or treating inner ear organs with the ototoxic antibiotic neomycin [22]. The sampling time points after the laser lesion (30 min, 1 h, 2 h and 3 h) were chosen in order to provide insights into the very early signaling events that occur after epithelial injury, while the sampling time points after a 24-h incubation of inner ear organ cultures with neomycin (immediately after incubation, and after 24 and 48 h in neomycin-free medium) were chosen to cover the period of the S-phase entry by supporting cells, which peaks at about 48 h after ototoxic injury in vitro. Unlike the first experiment, which measured responses to lesions in both cochlear and vestibular cultures, the second regeneration experiment [23] from the same group only measured expression in chick utricles. Also, the second experiment only explored the effect of an ototoxic antibiotic on gene expression. Still, this study had several advantages over [22]. First, it followed the expression changes across a 7-day time course, with more frequent measures in the 48–72 h window of regeneration, a period characterized by highly dynamic patterns of expression for many genes. Second, it employed RNA-seq instead of microarrays for the expression measuring, allowing the authors to obtain a comprehensive transcriptome, instead of specifically focusing on transcription factors expression as was done in the first experiment.
The comparison of the TFs was designed to detect pathways that are common to inner ear development and regeneration and specifically to reveal genes essential for either proliferation of supporting cells, or for transdifferentiation to hair cells. Out of 712 DE TFs in the first regeneration experiment [22], we mapped 596 to orthologous mouse genes. Intersection with our list of 64 TFs yielded 33 TFs that are involved in both development and regeneration (Additional file 1: Figure S4). Significantly, eight of these are also DE. The overlap with the later avian transcriptome experiment [23] was far more limited. Out of 212 DE TFs found in the experiment, we mapped 208 to orthologous mouse genes, of which only six appear also in our list of TFs (Additional file 1: Figure S5), and five of them are also DE. The TFs SMAD9 and SPI1 were DE in both avian experiments and were associated with enriched motifs in our developmental study.
Finally, we performed a comprehensive literature search for the motifs found in the context of inner ear development. For a small subset, the results are detailed in the following sections, with a more complete list available in Additional file 1: Supplementary Results.
Transcription factors affecting expression changes with age
The set of genes that were upregulated at E16.5 was enriched for binding sites for the motifs: Elk-1, Nrf-1, E2F-1, E2F, NF-Y, and AHRHIF of which, the subset Elk-1, Nrf-1, NF-Y, E2F-1 and some TFs associated with the motifs AHRHIF (Arnt and AhR) were upregulated together with their regulated genes. Upregulation of Hif1a, another motif associated with AHRHIF, could be detected at P0, suggesting that the upregulation of Arnt controlled genes is achieved by an increase in the formation of the heterodimer Arnt:AhR and not Arnt:Hif1a [37]. There was no enrichment of binding sites detected in the genes upregulated at P0.
ELK1 and TFs associated with AHRHIF are known to change their expression during regeneration. The expression of ELK1 was reported to increase 30 min after wounding cochlear hair cells with a laser, marking an early signaling event that occurs after epithelial damage [22]. Another TF that was increased after cellular insult was the AHRHIF TF, ARNT whose expression increased 24 h after exposing cochlear hair cells to neomycin, only to decrease again by 48 h, together with HIF1A and AHR. These time points reflect a change of expression in the supporting cells [22]. The transient increase of ARNT during regeneration resembles its transient expression pattern during normal inner ear development, between E13 and E17 in mouse cochlear epithelial cells [38]. Interestingly, the three TFs (ARNT, HIF1A, and AHR) were also reported to respond to tissue damage caused by a different toxic compound (TCDD [39]). E2F1 is an important pro-apoptotic TF [40], and under some mitochondrial stress, it engages apoptotic signals to cause deafness [41]. Regulation of transcription during the cell cycle is under the control of E2 factors (E2Fs), often in cooperation with nuclear factor Y (NF-Y) [42], another TF highlighted in this comparison. In utricle hair cell regeneration [23], E2F1 is changing its expression in a pattern that is associated with cell cycle genes. As mentioned, this TF and its targets are upregulated in E16.5, an age when we see enrichment for cell cycle activity.
Transcription factors affecting expression change between tissues
In contrast to the genes differentially expressed during the development, the set of genes upregulated in the cochlea was enriched for binding sites for the motifs: HIC1, E2F, ZNF219, ZF5, UF1H3BETA, MOVO-B, MAZ, VDR, MAZR, MTF-1, c-Myc:Max, AP-2, CAC-binding protein, ETF, E47, Lmo2 complex, RREB-1, LBP-1, CP2/LBP-1c/LSF, and Spz1. TFs associated with E2F, ZF5, and MAZ were significantly upregulated in the cochlea, while TFs associated with MOVO-B, VDR, and Lmo2 complex were upregulated in the vestibule. TFs associated with 11 of these 20 motifs (LBP-1, Lmo2 complex, E47, E2F, ZNF219, ZF5, VDR, MTF-1, c-Myc:Max, AP-2, CP2/LBP-1c/LSF) have been previously reported to change their expression during the regeneration of inner ear sensory epithelia [22, 23].
Focusing on genes that are altered both during development and regeneration, the enrichment of E2F noted indicates the presence of proliferation in the cochlea at the relevant period of development. Given the role of this TF family in inducing proliferation, their involvement in hair cell regeneration is not surprising and is currently the focus of active research [43]. In contrast, ZF5 is known primarily as a repressor of transcription and specifically as a regulator of cell cycle progression (through c-myc [44]), and cognitive development (through FMR1 [45]). Thus, the upregulation of expression in the cochlea, where its targets are also upregulated, was unexpected. This might indicate the existence of an additional activating role for ZF5, or that another TF is activating the transcription of these targets, and ZF5 is upregulated as part of a negative feedback loop. In avian hair cell regeneration, the expression of ZF5 in the cochlea increases late in the recovery from neomycin damage, suggesting a role in cochlear hair cell differentiation. Another TF with the same pattern of expression during hair cell regeneration is LMO2. Enrichments for LMO2 binding sites were found in the list of upregulated genes from both the cochlea and the vestibule. While the results of the regeneration experiment support a function for LMO2 in the cochlea, the expression of the TF in our experiment was higher in the vestibule. A possible explanation for this duality could be that LMO2 interacts with different partners in the two tissues, and thus, a different subset of genes is increased in each case. The Lmo2 complex typically contains a single GATA factor and a single TAL1/E47 heterodimer, but the GATA factor can be replaced by an additional TAL1/E47 heterodimer, resulting in a change in the genes regulated [46]. As Gata2 and Gata3 are upregulated in the vestibule and Tal1 is upregulated in the cochlea (DE q-values = 7.32 × 10−18, 1.67 × 10−175, and 5.29 × 10− 7, respectively), the complexes formed in each tissue might differ in composition. VDR is a transcription factor regulated by vitamin D levels [47]. Hypo- and hypervitaminosis D can cause sensorineural hearing loss [48]. The downregulation of the gene in the cochlea, where its targets are upregulated, suggests a repressor role for this TF, which is supported by existing literature [49]. In utricle hair cell generation [23], VDR’s expression peaks in the 54–72 h window after the aminoglycoside damage. This pattern makes it a candidate for playing a role in the phenotypic conversion process from supporting cells to hair cells in the vestibule. According to our experiment, it might fulfill a similar role in the vestibular development.
In the set of genes upregulated in the vestibule, we could detect enrichment for binding sites for the 21 motifs: HNF4, SREBP-1, NF-1, PEA3, TEF-1, AP-2rep, NF-kappaB (p65), LBP-1, LUN-1, E2A, PU.1, MyoD, Nrf2, Lmo2 complex, COUPTF, ISRE, HEB, E47, SMAD, AML-1a, and c-Ets-1. TFs associated with five motifs (TEF-1, PU.1, Nrf2, Lmo2 complex, and ISRE) were significantly upregulated in the vestibule, while TFs associated with four other motifs (PEA3, COUPTF, most SMADs, and AML-1a) were upregulated in the cochlea. TFs associated with 15 of the 21 enriched motifs (HNF4, SREBP-1, PEA3, NF-kappaB p65, LBP-1, E2A, PU.1, MyoD, Nrf2, Lmo2 complex, COUPTF, ISRE, HEB, E47, SMAD, c-Ets-1) displayed a change in expression during the regeneration experiments [22, 23].
The upregulation of TFs associated with the motifs Spi1 [PU.1] and Nfe2l2 [Nrf2] in the vestibule supports their role as inducers of transcription. However, the decrease seen in the cochlear expression in late (48 h) recovery from neomycin [22] suggests that their repression is required for proper differentiation of supporting cells to cochlear cells. SPI1 is known to be involved in hematopoietic development and induces proliferation of immune cells [50] and therefore might upregulate the immune functions that are enriched in the vestibule. Similarly, NFE2L2 can upregulate functions related to stress response and specifically to antioxidant defense [51]. The expression pattern of Nr2f1 and Nr2f2 associated with the COUPTF motif is in agreement with their suggested role as repressors of transcription, as they are downregulated in the vestibule, although the motif as a whole is enriched in the genes upregulated in the vestibule. Following laser damage, the expression of NR2F2 increases in the cochlea for 3 h and an increase in cochlear expression is also evident in late (48 h) recovery from neomycin [22]. Nr2f2 is known to work as a repressor of myogenesis, inhibiting MyoD [52], another TF whose targets are upregulated in the vestibule. Our data suggest that their repressive effect might have a role in cochlea development.
SMADs are intracellular proteins that transduce extracellular signals from transforming growth factor beta (TGF-β) ligands to the nucleus, where they activate downstream gene transcription [53]. Although TGF-β signaling is thought to be active in the cochlea, our results show rather that the downstream targets of this pathway are enriched in the vestibule. In order to address this issue, we examined the expression levels of individual SMADs. Most receptor-regulated SMADs (R-SMADs) were upregulated in the cochlea (Smad1, Smad2, Smad5, Smad9), in agreement with the hypothesis of higher TGF-β activity in the cochlea. However, inhibitors of this signaling pathway (Smad6 and Smad7) were also upregulated in the cochlea, and with relatively high FCs (1.9 and 1.6, respectively), and may be responsible for decreasing the transcription of the downstream genes in the cochlea compared to the vestibule. The story becomes more complex with the two intracellular pathways involving SMADs. The R-SMADS Smad2 and Smad3 mediate the response to TGF-β ligands, which participate in the regulation of inner ear development by retinoic acid [54]. Smad2 was upregulated in the cochlea, while Smad3 was upregulated in the vestibule. In the regeneration experiment [22], SMAD2 expression in the vestibule increased in a late response to neomycin damage in the utricle, emphasizing the importance of TGF-β signaling for vestibular differentiation. In a different pathway, the R-SMADS Smad1, Smad5, and Smad9 mediate the response to bone morphogenetic proteins (BPMs), which are involved in generation of inner ear sensory epithelia [55], as well as chondrogenesis [56]. All three were upregulated in the cochlea, with Smad9 showing a very impressive FC of 3.4. SMAD9 was also increased in response to late neomycin damage in the cochlea [22]. This, together with its high cochlear levels, implies that it plays a role in cochlear differentiation.
SPI1 and SMAD9 also change in expression during utricle hair cell regeneration [23]. The patterns of the expression are complex. Notably, their maximal deviations from the control are at time point 66 h, where SPI1 is upregulated and SMAD9 is downregulated. These changes are of opposite directions to those observed in cochlear regeneration in [22], agreeing with the tissue-specific roles of the two TFs.
Transcription factors affecting expression ratio change with age
In the set of genes for which the cochlea to vestibule expression ratio increases with age \( \left(\frac{Cochlea}{Vestibule}\uparrow \right) \), we could detect enrichment for the binding sites for the motifs HNF4, E47, and a group of nuclear receptors (LXR, PXR, CAR, COUP, RAR), AP-4, and SMAD. The expression ratio of Nr2f1, a COUP TF, increased significantly in the same direction as its targets, which might have a positive downstream effect on retinoic acid receptor (RAR) signaling [57]. Interestingly, TFs associated with all the motifs changed their expression during the regeneration of inner ear sensory epithelia [22].
Retinoid signaling is critical during inner ear embryonic development, as well as in the postnatal maintenance of its function [58]. Both vitamin A deficiency and intake of excess retinoic acid (RA) during pregnancy have been shown to cause malformations in ear development. In rodents, in utero exposure of fetuses to RA negatively affected the semicircular canals and the cochlea. Key components in retinoid signaling show spatiotemporal expression patterns, and the interactions that excess RA interferes with are dependent on the developmental stage. KEGG enrichment of our DE genes showed that metabolism of RA was higher in the vestibule and at P0. Taken together with the motif enrichment, we deduce that retinoid signaling is important to both cochlear and vestibular development, with its role in the cochlea becoming more prominent in the period between E16.5 and P0. In the hair cell regeneration experiment, the cochlear expression of the retinoid receptor RARA decreased 24 h after neomycin damage, but by 48 h, NR2F1 expression increased [22]. This later increase might mimic the increase in retinoid signaling seen in normal development.
Interestingly, we could detect enrichment of binding sites for AML-1a, LEF1, LBP-1, HEB, and POU6F1 in the set of genes for which the vestibule to cochlea expression ratio increases with age \( \left(\frac{Vestibule}{Cochlea}\uparrow \right) \). The expression ratio of Runx1 [AML-1a] increased significantly in the same direction as its targets. TFs associated with LBP-1 and HEB also changed their expression during the regeneration of inner ear sensory epithelia.
Comparison with transcription factors known to enhance mammalian hair cell regeneration
Previous studies that induced hair cell regeneration by coordinated manipulation of multiple factors, showed a better efficacy for the ectopic expression of ETV4, TCF3, GATA3, MYCN, or ETS2 in combination with ATOH1 over the overexpression of ATOH1 alone [18, 19]. In retrospect, our method identifies some of the TFs that were mentioned earlier for their ability to induce differentiation. Specifically, it singles out Etv4 [PEA3] and Tcf3 [E47]; Gata3 is a partner of the highlighted Lmo2; and Mycn and Ets2 have similar targets as c-myc [59] and Ets1 [60], respectively. Moreover, four out of five of the TF genes are DE between the cochlea and the vestibule (Etv4, Gata3, Mycn, and Ets4; q value ≤ 0.05), which makes them good candidates for experimental interventions, as explained above.
Change in the proportion of hair cells in sensory epithelia
Because of the difficulty of dissecting out the sensory epithelia and separating the hair cells from the adjacent supporting cells, all tissue samples of this type contain varying amounts of both hair cells and supporting cells. This complicates conclusions as to whether differential expression can be attributed to differences in the expression profiles or to variability in the cell mixture composition. In order to address this issue, we produced expression signatures of hair cells and supporting cells from a previous experiment [21] and used them to compute the proportion of each type in each preparation. We also evaluated the heterogeneity of cell types assuming that the cochlear sample is contaminated by cells from the vestibule (or utricle) and vice versa.
A different subset of genes was used to create the signatures for E16.5 and P0. For each age, we ranked the genes in decreasing order of expression variance across the four reference samples (cochlear and vestibular GFP+ and GFP− samples). We then took the expression of the first k genes in the list, where k equals 453 for E16.5, and 193 for P0. The value of k was chosen so that it minimized the estimated percentage of contamination in our mixed data, i.e., the estimated percentage of cochlear cells in vestibular samples plus the estimated percentage of vestibular cells in cochlear samples. We predicted that this heuristic approach would improve the overall prediction accuracy, although it did not directly optimize the precision of estimation of the percentage hair cells, which was our main goal. We used DeconRNASeq to estimate the mixing proportions [61].
The estimated proportions of hair cells were similar in both scenarios (Fig. 5), which allows us to ignore the issue of possible tissue contamination. The estimated percentages (±SD) are 32.6 (± 1.6) and 23.8 (± 1.0) in the cochlea and the vestibule at E16.5 and 44.0 (± 1.1) and 40.1 (± 0.2) in the cochlea and the vestibule at P0, respectively. These results indicate that the percentage of hair cells is higher in the cochlea at both ages and increases with development in both tissues, with the increase in the vestibule being more prominent (1.9-fold increase compared to 1.4-fold in the cochlea). Strikingly, in all estimations, the percentage of supporting cells was higher than 50%, suggesting that these cells have a dominant influence on the expression profiles.
Even with the value of k selected to minimize the contamination, our calculation gives 51.8% contamination in the cochlea at P0. We are unsure how to interpret this high number. Possible causes for an overestimate of contamination could be (1) experimental noise, either in our data or the data used to generate the expression signatures at P0, or (2) inaccuracy of the deconvolution method when the signatures are similar. The similarity of the signatures of the same cell type in the cochlea and the vestibule can be seen by the high correlation values (r = 0.65, or 0.83 for signatures of hair cells and supporting cells, respectively).