3D reconstruction of the mouse cochlea from scRNA-seq data suggests morphogen-based principles in apex-to-base specification

Abstract

In the mammalian auditory system, frequency discrimination depends on numerous morphological and physiological properties of the organ of Corti, which gradually change along the apex-to-base (tonotopic) axis of the organ. For example, the basilar membrane stiffness changes tonotopically, thus affecting the tuning properties of individual hair cells. At the molecular level, those frequency-specific characteristics are mirrored by gene expression gradients; however, the molecular mechanisms controlling tonotopic gene expression in the mouse cochlea remain elusive. Through analyzing single-cell RNA sequencing (scRNA-seq) data from E12.5 and E14.5 time points, we predicted that morphogens, rather than a cell division-associated mechanism, confer spatial identity in the extending cochlea. Subsequently, we reconstructed the developing cochlea in 3D space from scRNA-seq data to investigate the molecular pathways mediating positional information. The retinoic acid (RA) and hedgehog pathways were found to form opposing apex-to-base gradients, and functional interrogation using mouse cochlear explants suggested that both pathways jointly specify the longitudinal axis.

Keywords: 3D reconstruction; organ of corti; retinoic acid; sonic hedgehog; tonotopy.

Figure 1.. Time-space vs morphogen hypothesis testing.

(A) Model of two different genetic timers hypothesized to confer apex-to-base identity in the cochlea (adopted from Negrete & Oates). A genetic timer represents a gene regulatory network that changes gene expression dynamics in a sequential manner. Patterning occurs when the timer is arrested in different cells, which can be triggered by cell cycle exit or morphogen concentration. The arrest results in spatial gene expression signatures, which are represented by different colors for apex (blue) and base (orange). In timers modulated by asymmetrical cell division, the daughter cell’s genetic timer is arrested at the time of division, while the timer of the progenitor cell continues to run. In the morphogen model, the timer is controlled by the concentration of the morphogen, which allows for dynamic adjustment of the positional information as the duct extends over time. Dotted yellow lines between the E12.5 and E14.5 cells indicate highest transcriptional similarity. (B) Experimental paradigm to test the hypotheses. (C) UMAP plots of cells isolated at E12.5 and E14.5 with cell type annotation. Cell types constituting the cochlear duct are color coded: roof (red), cochlear floor segregating into apex (blue) and base (orange). (D) Heatmap showing the E12.5 expression level of the overlapping GEGs between the E12.5 and E14.5 time points. Columns represent individual cochlear floor cells that are rank ordered along the base-to-apex (x-) axis. Rows depict GEGs ordered by hierarchical clustering. Genes highly expressed in the apex are outlined in purple, and genes highly expressed in the base are in yellow. Gene expression levels are color coded from low (blue) to high (red). The heatmap is scaled by rows. (E) Line plots showing the LOWESS regressed expression value of the selected GEGs along the base-to-apex (x-) axis, as represented by PC1. Hmga2 and Hmgb2 are two GEGs with opposing gradients. The lines are color coded by age: E12.5 (green) and E14.5 (brown). (F) Representative E14.5 histological sections of the cochlea stained for Hmga2 (green) and Hmgb2 (red) mRNA with FISH. Upper left: overview. Remaining panels: insets of the cochlear floor at apex, mid, and base locations. A minimum of two samples per probe were analyzed and yielded similar results. (G) Heatmap showing the affinity between E12.5 metacells and E14.5 metacells. Both columns and rows represent metacells ordered along the apex-to-base axis. Color coding indicates similarity from low (blue) to high (red). The heatmap is scaled by rows. Scale bars in (G): 200 µm (overview) and 10 µm (insets). See also Figure S1 and S2.

Figure 2.. 3D spatial reconstruction of the developing cochlear duct.

(A-D) Visualizations of the cochlear duct with spatial domains color coded for E12.5 (1^st row) and E14.5 (2^nd row). (A) 3D schematic representation of spatial domains asymmetrically distributed across the cochlear duct. Shown are 3D representations of the cochlear duct (top) and 3D approximations of the cochlear floor without the roof (bottom). (B) Spatial domains color coded on 2D scaled PCA plots. PC1, corresponding to the longitudinal axis, was projected onto the y-axis for intuitive visualization. PC2, scaled by PC3, was projected onto the x-axis to distinguish the spatial domains of the cochlear duct. The inset shows the same 2D scaled PCA plot with the original cluster IDs and color code as in Figures 1B, C. (C) Spatial domains visualized with 3D cylinder projection. The height of the cylinder corresponds to the apex-to-base axis of the cochlea with the apex facing up. The inset shows the original cluster IDs and color code projected onto the same cylinder. (D) Spatial domains visualized on 2D circular projections with a flattened apex-to-base axis. The lower half corresponds to the cochlear floor, whereas roof cells form the upper half of the reconstruction. Basal cells are located at the outer perimeter of the projection, and apical cells are located at the center. Medial cells contribute to the left side, whereas lateral cells form the right side. The inset highlights original cluster IDs with original color code. (E-G) Three major body axes of the cochlear duct, visualized according to canonical marker gene expression. (E) Oc90 and Lum segregate along the roof-floor axis. (F) Fgf10 and Bmp4 expression indicate the medial-lateral axis in the cochlear floor. (G) The apex-base axis is represented by opposing expression of Fst and Mki67. Gene expression levels are color coded from absent (gray) to low (purple) to high (yellow). See also Figure S3 – S5.

Figure 3.. RA signaling components in the developing cochlea.

(A) Heatmaps showing enrichment scores of the top 100 differentially enriched GO terms for the E12.5 and E14.5 time points. Columns represent single cells ordered along the base-to-apex (x-) axis. Rows represent differentially enriched GO terms with ordering based on hierarchical clustering. Significance: adjusted P < 0.05 (two-sided Wilcoxon rank-sum test with BH correction). Heatmaps are scaled by rows. (B) Heatmap showing statistical significance for ten RA-associated GO terms. Statistical comparisons include apex versus base, and floor versus roof compartments for the E12.5 and E14.5 time points. High significance, green; low significance, yellow; not significant (n.s.), gray. (C) Retinoic acid biosynthetic process (GO:0002138) GO term enrichment scores projected onto E12.5 (1^st row) and E14.5 (2^nd row) circular plots. The inset is color coded by floor (pink) and roof (red). The boxplots show significantly different activity scores between the floor and roof compartments for both time points. E12.5: ***P=1.62E-18. E14.5: ***P=9.92E-33 (two-sided Wilcoxon rank-sum test with BH correction). (D-E) mRNA expression levels projected onto E14.5 circular plots. (D) Rdh10 is differentially expressed in roof cells. ***P=2.58E-12 (two-sided Wilcoxon rank-sum test with Bonferroni correction). (E) Aldh1a3 is differentially expressed in roof cells. ***P=2.88E-46 (two-sided Wilcoxon rank-sum test with Bonferroni correction). (F) Crabp2 is differentially expressed in cochlear floor cells. ***P= 7.65E-15 (two-sided Wilcoxon rank-sum test with Bonferroni correction). (G) Cellular response to retinoic acid (GO:0071300) GO term projected. The inset is color coded by apex (blue), base (orange), and roof (red). Boxplots show significantly different activity scores between the apical and basal compartments for both time points E12.5: ***P=2.42E-05. E14.5: ***P=7.28E-06 (two-sided Wilcoxon rank-sum test with BH correction). (H) Cyp26b1 is gradually expressed in a narrow band of cells centered in the prosensory domain extending along the tonotopic axis at E14.5. The difference was not significant when all cells annotated as cochlear floor were considered in the statistical comparison after Bonferroni correction. (I) Significantly higher Dhrs3 expression in apical cells than basal cells. **P=4.87E-03 (two-sided Wilcoxon rank-sum test with Bonferroni correction). A minimum of three sections was analyzed and yielded similar results. The box plots show the interquartile range (box limits), median (centre line), minimum to maximum values (whiskers). (J) β-galactosidase histochemical staining of E12.5, E14.5, and P0 cochleae. RARE-lacZ signal was localized to the basal end of the cochlea. RARE-lacZ staining was also found in the spiral ganglion along the entire length of the cochlea, though at lower levels compared to the basal cochlea. Three samples were analyzed and yielded similar results. Scale bar in (J) 200 μm. See also Figure S6.

Figure 4.. RA gradient in vivo and functional RA signal transduction in E14.5 cochlear explants.

(A) Representative confocal micrographs of E14.5 cochlear explants cultured for 3 h in control medium and medium supplemented with 500 nM ectopic RA. Specimens were stained for Cyp26b1 mRNA, counterstained with DAPI, imaged in whole mount format, and analyzed in Fiji software. Insets represent higher magnification images that were collected at representative apex and base locations. (B) Quantification of Cyp26b1 puncta normalized to 1000 µm² indicated significantly higher Cyp26b1 transcript counts for the control apex (n=6) than the control base (n=6). **P=2.8E-03 (two-sided paired t-test). In addition, Cyp26b1 expression significantly increased in the basal compartment after exposure to 500 nM RA for 3 h (n=6) compared with the control treatment (n=6). **P=4.1E-03 (two-sided unpaired t-test). The box plots show the interquartile range (box limits), median (centre line), minimum to maximum values (whiskers), and biologically independent samples (circles). Scale bars in (A), first panels: 200 μm. Scale bar in remaining panels, same for all insets, 10 μm.

Figure 5.. SAG enhances Cyp26b1 expression in E14.5 cochlear explants.

(A) PCA-based rank ordering of individual cells of the E14.5 cochlear floor along the apex-to-base (x-) axis. Data points are randomly spread along the y-axis for better visualization (top panel). Color code: cluster identity apex (blue) and base (orange). Activity scores for the smoothened signaling pathway (GO:0060831) and retinoic acid biosynthetic process (GO:0002138), visualized with LOWESS regression lines, show opposing gradients along the tonotopic (x-) axis. In addition, Cyp26b1 mRNA expression levels show a positive correlation with SHH related GO term (GO:0060831). (B) Representative micrographs of E14.5 cochlear explants cultured for 72 h in control medium or medium supplemented with the smoothened agonist SAG (1 µM). Explants were stained for Cyp26b1 transcripts and counterstained with anti-SOX2 antibody and DAPI to enable imaging of the prosensory domain. High resolution images were taken at representative apex, mid, and base locations, and Cyp26b1 transcripts were quantified in Fiji software. (C) Cyp26b1 expression levels averaged across the total length of the explant were significantly higher in the SAG group (n=27) than the DMSO control group (n=18). Data are normalized to 1000 µm². ***P=1.7E-04 (two-sided unpaired t-test). The box plots show the interquartile range (box limits), median (centre line), minimum to maximum values (whiskers). Scale bars in (B): upper panels, 200 μm; lower panels, same for all insets, 10 μm.

Figure 6.. Manipulation of RA and SHH signaling results in mirrored phenotypes.

(A-E) Representative confocal micrographs of E14.5 cochlear explants cultured for 72 h in (A) control medium or medium supplemented with (B) ectopic RA, (C) the RA receptor blocker AGN 193109, (D) smoothened agonist SAG and (E) smoothened antagonist SANT-1. Explants were stained for Hmga2 mRNA and counterstained with anti-SOX2 antibody and DAPI. Sox2⁺ cells of the prosensory domain were imaged at representative locations in apical and basal compartments of the explant, and Hmga2 puncta were analyzed in Fiji software. (F) Quantification of Hmga2 puncta normalized to 1000 µm² indicated significantly higher Hmga2 transcript counts for the control apex (n=11) than the control base (n=11). *P=0.0247 (two-sided paired t-test). With the DMSO control as a baseline, differential Hmga2 expression was determined for the following conditions: exposure to 0.5 µM RA for 72 h resulted in significantly lower Hmga2 expression in the apical compartment (n=11) than the apical control compartment (n=11). *P=0.0195 (one-way ANOVA with post-hoc Tukey’s test). Similarly, blocking SHH signaling with SANT-1 (n=10) resulted in significantly lower Hmga2 in the apical compartment than in the apical control compartment (n=11). *P=0.018 (one-way ANOVA with post-hoc Tukey’s test). In contrast, application of the RA-receptor blocker AGN 193109 (n=10) resulted in significantly higher Hmga2 in the basal compartment than in the control basal compartment (n=11). *P=0.0221 (one-way ANOVA with post-hoc Tukey’s test). The SHH agonist SAG (n=9) also resulted in greater Hmga2 expression in the basal compartment than in the control basal compartment (n=11). *P=0.0326 (one-way ANOVA with post-hoc Tukey’s test). The box plots show the interquartile range (box limits), median (centre line), minimum to maximum values (whiskers), and biologically independent samples (circles). (G) Summary diagram of the inverse relation between RA and SHH signaling gradients, which is a key component in controlling the tonotopic expression of Hmga2. Scale bars (A-E): 10 µm.