Post-hatch day 1 (P1) chickens were euthanized by CO₂ inhalation followed by decapitation. The head was hemi-sectioned and the cochlea duct was isolated in cold DEPC-treated PBS. The basilar papilla was dissected from the cochlear duct and divided into three pieces: base, middle, and apex. Each piece was collected in a separate Eppendorf tube containing Trizol reagent (Invitrogen, Carlsbad, CA, United States) and quick-frozen to minimize RNA degradation. Total RNA from 12 chickens (24 cochleas) were used for 3 biological replicates for each region (i.e., 4 chickens or 8 cochleas per replicate). RNA was isolated using Trizol reagent according to the manufacturer’s instructions.

We performed an RNA-seq experiment for three sample groups on the base, middle, and apex regions of chick cochlea. Following cochlea dissection and mRNA extraction, the quality and quantity of RNA were examined using Tapestation2100–RNA Screen Tape (Agilent) and a Victor3 Multilavel plate reader (PerkinElmer), respectively. 1μg of total RNA from each sample was reverse transcribed into cDNA and amplified using the Illumina TruSeq RNA Sample Prep Kit v2. For library preparation, mRNA was enriched by oligo-dT to minimize DNA contamination. The PCR products for the cDNA library were purified and selected using the AMPure XP bead (Beckman Coulter). The DNA was quantified using Tapestation2100–D1000 Screen Tape. The cDNA libraries were multiplexed and sequenced on the Illumina Hi-Seq 4,000 platform, and 101-bp paired-end reads were generated. The files from the multiplexed RNA-seq samples were demultiplexed, and fastq files for each sample were generated. Approximately 57–82 million read pairs per replicate were generated. To quantify the transcript abundance, sequencing reads were aligned to the chicken genome (galGal5) and the transcriptome (Ensembl GTF release 91), using the STAR aligner (v2.5.4b) (Dobin et al., 2013), allowing up to 3 mismatches per read.

We identified DEGs on base, middle, and apex regions using Cuffdiff (v2.2.1) (Trapnell et al., 2012). We calculated gene expression level as fragments per kilobase of transcript per million mapped reads (FPKM) of each gene based on the length of the gene and the number of reads mapped to the gene. The values of FPKM for genes are included in the Supplementary Table S1. Significant DEGs between two sample groups were identified at FDR <5%, >1.5-fold difference between mean FPKMs, and mean FPKM ≥1.0 in at least one of all three sample groups. These were further categorized depending on the changes of the values of FPKM as Up-trend and Down-trend as below (see Figure 1B for a diagram).

• Up-trend includes genes that show at least a 1.5-fold increase in mean FPKM values with directionality from base to middle to apex, allowing up to 15% variability against the directionality.

We searched differential AS events from the base, middle, and apex regions from RNA-seq data. We identified AS events between two sample groups using rMATS (v3.2.5) (Shen et al., 2012; Park et al., 2013; Shen et al., 2014), which can detect five major types of AS events. Briefly, rMATS uses a modified version of the generalized linear mixed model to detect differential AS from RNA-seq data with replicates. It accounts for exon-specific sequencing coverage in individual samples as well as variation in exon splicing levels among replicates. For each AS event, we used both the reads mapped to the splice junctions and the reads mapped to the exon body as the input for rMATS. In each rMATS run, the first sample group was compared to the second sample group to identify differentially spliced events with an associated change in the percent spliced in (ΔPSI or Δψ) of these events. rMATS was run using the -c 0.0001 parameter to compute p values and FDRs of splicing events with a cutoff |Δψ| of >0.01%, then significant splicing events were identified at FDR <5% and |Δψ| ≥ 5% as described previously (Yang et al., 2016).

To validate DEGs and AS events detected from RNA-seq data, we performed reverse transcription (RT)-PCR using extracted total RNA as described previously (Warzecha et al., 2010). Total RNA from chicken P1 cochlear base, middle, apex were isolated using Trizol reagent (Invitrogen, Carlsbad, CA, United States). Reverse transcription (RT) was performed with 1 μg of RNA using ImProm-ll™ Reverse transcriptase (Promega, Madison, WI, United States). A duplicate set of PCR was performed using IP-Pro DNA Taq polymerase (Cosmogenetech, Seoul, Korea). The primers used for DEG and AS event validation, and the PCR product sizes are summarized in the Supplementary Table S2 and Supplementary Table S3 respectively.

We sought to identify binding sites of splicing factors and other RBPs that were significantly enriched in differential exon skipping events. To assess the enrichment of RBPs near alternatively spliced exons, the differentially spliced exons from rMATS runs were compiled. These differentially spliced exons were scanned for RNA binding protein motifs using rMAPS2 (Hwang et al., 2020) with default parameters.

We conducted two sets of GO analyses: the one with DEG results and the other with AS results. 1,479 Up-trend genes and 797 Down-trend genes from DEG analysis were compiled for the first set of GO analysis using PANTHER (Mi et al., 2021) with GO Ontology database release 2021–07–02 (The Gene Ontology, 2019). The host genes of significant skipped exon events from rMATS output were collected for the second set of GO analysis. For both sets of GO analysis, we examined BP (biological processes), MF (molecular functions), and CC (cellular components).

To evaluate phenotypes associated with different regions along the longitudinal axis of the cochlea duct, we extracted mRNA from three different regions of the chick cochlea, base, middle and apex, and performed an RNA-seq experiment. This yielded 57 to 82 million read pairs per replicate, which were then aligned to the chicken genome. The transcript abundance was quantified, and FPKM values for each sample were obtained. Of 24,880 Ensembl-annotated genes, there were 13,045 (52.4%) in which at least one of the 9 samples had an FPKM ≥1. To assess position-specific gene expression that might be responsible for configuration and maintenance of the tonotopic axis of the cochlea duct, we identified differentially expressed genes on the three regions. We compared gene expression profiles from two different regions of cochlea at a time for all three regions. Since gene expression profiles originate from the same cochlea tissue despite being extracted from different regions, they should show intrinsically much similar gene expression patterns compared to those from different organs. Thus, we searched for DEGs with more relaxed criteria. We considered a 1.5 fold-change in expression level as a significant difference along three different regions.

To identify genes that are responsible for tonotopic configuration along the cochlea duct, we searched for genes that gradually increase their expression by a factor of at least 1.5 from base to middle to apex regions and designated them as Up-trend. Likewise, we searched genes that gradually decrease their expression at least 1.5-fold from base to middle to apex regions and designated them as Down-trend. For these gene expression gradients, we identified 1,479 Up-trend genes and 797 Down-trend genes. The gene expression patterns along the regions of the cochlea are represented in Figure 2A which clearly shows up-trend (blue to red) and down-trend (red to blue).

Among those genes, we selected several genes that have been associated with cochlear development and hearing function for validation by RT-PCR. For Up-trend genes, we chose BMP7, DNER, KCNJ2, TECTB, and TMC2 (Figure 2B). BMP7 and DNER were shown to be expressed gradually higher towards the apex and regulate hair bundle morphology along the tonotopic axis in chicken basilar papilla (Kowalik and Hudspeth, 2011; Mann et al., 2014). KCNJ2, which encodes an inward rectifier potassium ion channel, is known to be a bona fide marker of apical hair cells in chicken cochlea (Son et al., 2015). TECTB, which is categorized as a supporting cell marker in chicken cochlea (Janesick et al., 2021), was shown to be expressed in a base-to-apex increasing gradient in mouse cochlea (Son et al., 2012; Yoshimura et al., 2014). TMC2, which encodes a component of the stereocilia mechanotransduction channel complex, was shown to be expressed higher at the apical regions of neonatal chicken and mouse cochlea (Kurima et al., 2015; Beurg et al., 2018; Janesick et al., 2021). For Down-trend genes, BMPR1B, CHRDL1, GJA1, and IRX2 were selected. BMPR1B, encoding one of the receptors transducing BMP signaling, was shown to be expressed in the chicken hair cells (Lewis et al., 2018). CHRDL1 was shown to be expressed higher at the base and promote basal hair bundle phenotypes in chicken basilar papilla (Frucht et al., 2011; Mann et al., 2014). Mutation in GJA1 is shown to be associated with non-syndromic autosomal recessive deafness in human (Liu et al., 2001). IRX2 was shown to be expressed higher in the base of chicken basilar papilla (Kowalik and Hudspeth, 2011; Janesick et al., 2021). We confirmed that their gene expression gradient pattern along the different regions of cochlea observed in RNA-Seq data agrees with the RT-PCR result as shown in Figures 2B,C. These results suggest that their differential expression along the cochlea duct is closely associated with the configuration and maintenance of the characteristic functions in tonotopy.

We next examined the expression patterns of AS factors including PTBPs, ELAVL, and ESPRs. RNA-seq data show that PTBP3, ESRP1 and ESRP2 are included in Up-trend genes and RT-PCR results confirmed their increasing expression gradient from base to apex (Figures 2B,C). Other AS factors including PTBP1, PTBP2, and ELAVL4 were not included either Up- or Down-trend genes and did not show obvious expression gradient along the tonotopic axis by RT-PCR (Figures 2B,C). We reasoned that these AS factors take important roles in regulating gene expression in those regions in a sophisticated manner. In addition, considering subtle changes in controlling tonotopic differentiation along the longitudinal axis of the cochlea, we hypothesized that AS contributes to varying gene expression.

We carried out differential AS analysis on paired-end RNA-Seq data from three regions of the chick cochlea. We used rMATS to identify five common types of AS events by comparing samples from two different regions, base vs middle, middle vs apex and base vs apex. Then, significant splicing events were identified at FDR <5% and |Δψ| ≥ 5% as described previously. The most common type of AS event was skipped exon events (SE), in which a cassette exon around two closely located exons is either included or skipped. From the three comparisons, we identified 275, 114, and 325 significant SE events, respectively (Figure 3A and Supplementary Table S4). Of these, 122 SE events were identified as significant on both base vs middle and base vs apex comparisons. Likewise, 32 SE events were significant on both base vs apex and middle vs apex comparisons and 28 SE events were significant on both base vs middle and middle vs apex comparisons.

We hypothesized the genes with AS events along the longitudinal axis of chick cochlea are related to tonotopic configuration and maintenance. Similar to Up-trend and Down-trend DEGs, some SE events took place along the axis of the cochlea in either an increasing or decreasing manner, from base to middle to apex, as summarized in Supplementary Table S4. We further verified these cochlear position-dependent SE events by RT-PCR for genes known to be crucial for hair cell development and function including CDH23, EPB41L3, KCNMA1, and LMO7 (Figures 3B,C). CDH23 is a component of the tip-link that is directly connected to the stereocilia mechanotransduction channel complex, mutations of which causes nonsyndromic (DFNB12) and syndromic (USH1D) hearing loss in humans (Di Palma et al., 2001). Interestingly, CDH23 shows a base-to-apex increasing inclusion level and decreasing skipping level of the exon 32, suggesting that the protein region encoded by the exon 32 may provide the tonotopic properties of the tip-link crucial for frequency discrimination. EPB41L3, which was shown to be expressed in the tips of tallest stereocilia of mouse hair cells (Okumura et al., 2010), displayed higher inclusion levels in the middle of the cochlea compared to the base and apex (Figure 3B). KCNMA1 encodes for the alpha subunit of the large conductance calcium-activated potassium channel (BK channel; also known as Slo), and its alternative splicing has been suggested to be important for frequency tuning along the tonotopic axis by influencing the channel gating properties (Navaratnam et al., 1997; Rosenblatt et al., 1997; Frucht et al., 2011). We also observed that KCNMA1 showed differential exon usage along the tonotopic axis (Figure 3B). Our AS analysis also revealed that several different AS events can occur in the same gene. For example, LMO7, which was shown to be crucial for normal cuticular plate of hair cells and hearing function (Du et al., 2019), exhibited base-to-apex decreasing inclusion levels for exon 24 and the opposite base-to-apex increasing inclusion levels for exon 16 (Figures 3B,C).

In addition, we examined whether AS events are related to nonsense-mediated mRNA decay (NMD) pathway by introducing stop codons unexpectedly out of the original frame due to non-canonical AS events. As transcripts that have retained intron (RI) often contain premature stop codons, which can trigger the NMD (Yap et al., 2012), we searched differentially spliced RI events, then investigated the steady-state mRNA levels on different regions of the cochlea using our RNA-seq data. At FDR<=5%, we compiled 6 differentially spliced RI events and examined expression level changes of the host genes along the tonotopic axis. Genes that contain introns in the base region (5 genes with low expression level in the base region) showed an increase in their expression levels in the mid and apex regions while the gene that does not contain intron in the base region showed a decrease in its expression levels in the mid and apex region (Supplementary Figure S1). This pattern suggests that splicing change can regulate expression levels of genes by inducing the NMD pathway.

These results suggest that differential AS events along the cochlear duct act as a fundamental mechanism for establishing the tonotopy by altering morphological and functional properties of the proteins essential for frequency tuning. Although these genes represent gene expression gradients in AS events, they are not necessarily directly related to gene expression gradients at the transcription level.

We further evaluated the characteristics of the genes detected from DEG analysis with Gene Ontology (GO) enrichment analysis using PANTHER (Mi et al., 2021). For the 1,479 Up-trend DEGs, enriched GO biological process terms are related to inner ear development, sound perception and microtubule (cilium)-related processes. Consistent with these terms, enriched molecular function terms include ion channel activity and transmembrane transporter activity and enriched cellular component terms are related to stereocilium and cilium. Meanwhile, for the 797 Down-trend DEGs, enriched GO biological process terms are related to the overall system developmental process, signal transduction, and extracellular structural organization. Similarly, enriched molecular function and cellular component terms include signaling receptor binding and extracellular component binding (Figures 4A,B). These results, which show that different GO terms are enriched between Up-trend and Down-trend DEGs, suggest that distinct groups of genes play more prominent roles in frequency tuning depending on the frequency ranges along the tonotopic axis. Genes with these enriched GO terms are listed in Supplementary Table S5.

We next performed GO enrichment analysis on genes with differential skipped exon events using PANTHER (Mi et al., 2021) (Figure 4C). Enriched GO biological process terms are mainly related to actin cytoskeleton organization, ion transport, and synapse organization. Consistent with these, enriched molecular function terms are related to cytoskeletal protein binding, and enriched cellular component terms are related to cell projection, the transport vesicle, and the synapse. These results suggest that differential AS events are essential for establishing the tonotopy by regulating the structural and functional characteristics of actin-based cytoskeletal structures such as stereocilia, ion channels essential for mechanoelectrical transduction and synaptic architectures transmitting the electrical signals to the brain. Genes with these enriched GO terms are listed in Supplementary Table S6.

To examine the relationship between expression gradients of AS events and regulatory RNA binding proteins (RBPs), we performed motif enrichment analysis using rMAPS2 (Hwang et al., 2020) with differential AS data generated from rMATS. We uploaded skipped exon events data from each comparison along the tonotopic axis to the rMAPS2 web server and generated RNA maps for over 200 RBPs. These RNA maps represent enrichment of motif sequence near the target exons of the SE events. Interestingly, PTB like (CT rich) binding motifs were enriched near alternatively spliced exons (Figure 5). PTB like binding motifs were enriched within the target exon and the downstream flanking intron regions when the inclusion of the target exon is promoted (solid red line shows peaks with small p-value) while the PTB like motifs were enriched in the upstream flanking intron regions when the inclusion of the target exon is suppressed (solid blue line shows peaks with small p-values). We also examined the enrichment of ESRP 1/2 motifs (Supplementary Figure S2). There was a strong enrichment for these ESRP motifs (solid red line) in the upstream intron of upregulated exons. The differential expression of the ESRP 1/2 and PTBP3 together with the enrichment of their binding motifs and their positional enrichment relative to the regulated exons further suggest that the tonotopic related AS may be mediated by these splicing factors.