Characterization of microRNA and gene expression in the cochlea of an echolocating bat (Rhinolophus affinis)

Abstract MicroRNAs (miRNAs) are important post‐transcriptional regulators of gene expression and play key roles in many biological processes, such as development and response to multiple stresses. However, little is known about their roles in generating novel phenotypes and phenotypic variation during the course of animal evolution. Here, we, for the first time, characterized the miRNAs of the cochlea in an echolocating bat (Rhinolophus affinis). We sampled eight individuals from two R. affinis subspecies with significant echolocation call frequency differences. We identified 365 miRNAs and 121 of them were novel. By searching sequences of these miRNAs precursors in multiple high‐quality mammal genomes, we found one specific miRNA shared by all echolocating bats but not present in all other nonecholocating mammals. The targeted genes of this miRNA included several known hearing genes (e.g., KCNQ4 and GJB6). Together with the matched mRNA‐seq data, we identified 1766 differentially expressed genes (DEGs) between the two subspecies and 555 of them were negatively regulated by differentially expressed miRNAs (DEMs). We found that almost half of known hearing genes in the list of all DEGs were regulated negatively by DEMs, suggesting an important role of miRNAs in call frequency variation of the two subspecies. These targeted DEGs included several important hearing genes (e.g., Piezo1, Piezo2, and CDH23) that have been shown to be important in ultrasonic hearing of echolocating mammals.

Evolution of new miRNAs is sometimes coincided with the origin of novel phenotypes, such as the cortical development in primates (Kosik & Nowakowski, 2018). More studies in other wild groups are needed to test whether miRNAs play important roles in generating novel traits.
Gene expression regulation has been considered as a key mechanism underlying phenotypic variations (Carroll, 2008;Jacob & Monod, 1961;Mank, 2017). Thus, miRNAs, as important posttranscriptional regulator of gene expression, may also play important roles in phenotypic evolution of wild organisms. Up to now, few such studies have been conducted, such as in rapid phenotypic diversification of the Midas cichlids (Franchini et al., 2016(Franchini et al., , 2019. Bats have multiple novel traits, such as flight, echolocation, and extreme longevity (Huang et al., 2016). Recently, a number of novel miRNAs in bats has been identified to have potential contributions to their extreme longevity (Huang et al., 2016. However, little is known about the role of miRNAs in echolocation. Bats use echolocation to explore environments and detect prey (Schnitzler et al., 2003).
Horseshoe bats (Rhinolophidae) in particular have a specialized echolocation system. They emit long constant frequency components in their echolocation calls with the most energy contained in the second harmonic of the pulses (Siemers et al., 2005). These bats have auditory foveae responsible for response to frequencies of echolocation calls (Schuller & Pollak, 1979). When stationary, they omit calls at frequencies (resting frequency) that match to the frequencies of their hearing system (i.e., acoustic foveae) (Jones & Siemers, 2011).
Thus, the resting frequency of horseshoe bats' echolocation calls can be considered as a morphological parameter associated with their acoustic foveae (Siemers et al., 2005). Considerable variations in peak frequencies of echolocation calls have been found among different horseshoe bats (Zhang et al., 2009) and also among different subspecies or geographic populations of several horseshoe bats (Jacobs et al., 2017;Sun et al., 2013). Using RNA-seq data of cochlea tissue, the genetic basis of intraspecific variation in echolocation call frequency has been investigated among three R. ferrumequinum geographic populations (Zhao et al., 2019) and among three recently diverged subspecies of R. affinis (Sun et al., 2020). These two studies supported important roles of gene expression changes in intraspecific call frequency variations. However, little is known about the role of miRNAs in echolocation call frequency variation.
In this study, we focused on our previous study system, R. affinis.
This species has three subspecies in China which diverged recently, <1 million years ago (Mao & Rossiter, 2020 there are no significant differences in body size (e.g., forearm length) among the three subspecies (Mao et al., 2014). Specifically, R. a. hainanus and R. a. macrurus have similar echolocation call frequencies (71-74 kHz), consistent with the close phylogenetic relationships between them (Mao et al., 2010(Mao et al., , 2013, whereas R. a. himalayanus shows over 15 kHz higher call frequency than each of them (Mao et al., 2010(Mao et al., , 2013(Mao et al., , 2014. Our previous phylogenetic and phylogenomic studies have detected extensive introgression of mtDNA between R. a. himalayanus and R. a. macrurus in their secondary contact region (Mao et al., 2013;Mao & Rossiter, 2020). Thus, we chose two allopatric subspecies (R. a. himalayanus and R. a. hainanus) as the system to assess the role of miRNAs during the divergence of echolocation call frequency to reduce the effect of introgression on patterns of gene expression (Dannemann et al., 2017;McCoy et al., 2017).
Here by characterizing miRNAs and mRNAs in the cochlea tissue from two subspecies of an echolocating bat (R. a. himalayanus and R. a. hainanus), we test the following two specific hypotheses. First, we hypothesize that there might be specific miRNAs which are essential in the evolution of echolocation. If the answer is yes, we would expect to identify shared miRNAs in all echolocating bats or echolocating mammals. Second, we hypothesize that miRNAs might be involved in regulating echolocation call frequency variation of the two subspecies.
Following this hypothesis, we would expect to see multiple known hearing genes in the list of differentially expressed genes that are regulated by differentially expressed miRNAs identified between the two subspecies.

| MiRNA discovery
We used the following pipeline to identify miRNAs ( Figure 2). Raw small RNA sequencing data were processed using Cutadapt v2.6 (Martin, 2011) by trimming adapter sequences and removing reads with unknown bases. Low-quality reads were further filtered out using TRIMMOMATIC version 0.38 (Bolger et al., 2014)  to the known category. The final novel miRNAs in this study mean those that have not been reported in any organisms until now.

| Identification of specific miRNAs in echolocating mammals
This is the first study to characterize miRNAs in cochlea from an echolocating mammal. To identify candidate miRNAs specific to echolocating mammals, we searched all miRNAs identified here in genomes of other echolocating and nonecholocating mammals ( Figure 3). First, all miRNA precursors were mapped to genomes of one nonecholocating bat (Rousettus aegyptiacus) and five other nonecholocating mammals (human, mouse, dog, goat, and horse).
Second, unmapped miRNA precursors were aligned to recently published high-quality genomes of five echolocating bats (Rhinolophus

| miRNA differential expression analysis
To perform differential expression analysis (Figure 2), trimmed small RNA reads from each sample were mapped to the mature and precursor miRNAs identified in R. affinis above using the script quanti-  Figure 4b). Prior to differential expression analysis, we filtered out the lowly expressed miRNAs with the average CPM < 1 among all samples. Then, differentially expressed miRNAs between the two subspecies were determined using both DESeq2 and edgeR (Robinson et al., 2010) with significant results (Benjamini & Hochberg, 1995, padj < 0.05 in DESeq2 and FDR < 0.05 in edgeR) and |log 2 (fold change)| > 1. To test whether different amounts of sequencing data generated for each sample affected the current results, we repeated the differential expression analysis based on similar amounts of data for each sample. For this, we reduced sequencing data of each sample to 4.6 million reads (the lowest amount of data, Dryad file T3). Echolocating mammal (dolphin):

Platanista minor
Zero miRNA in all echolocating mammals

| mRNA differential expression analysis
We used the following pipeline to analysis mRNA-seq data ( Figure 2). Specifically, raw reads from mRNA-seq were processed using TRIMMOMATIC v0.38 (Bolger et al., 2014) (Figure 4d). Then, we filtered out those lowly expressed transcripts with the average CPM < 1 among all samples. Differentially expressed genes (DEGs) between the two subspecies were also identified using both DESeq2 and edgeR with the same criteria as in the miRNA above.

| MiRNA target prediction
To assemble transcripts expressed in cochlea, we applied a reference-guided assembly method using StringTie v2.1.5 (Pertea et al., 2015) with R. affinis genome as the reference. Filtered RNAseq data from all eight samples (see below mRNA differential  (Trapnell et al., 2012). Then, the transcripts were annotated by searching against the coding sequences (CDSs) of 20,180 genes retrieved from the R. affinis genome using BLASTN with e-value 10 −6 , sequence identify >90%, and alignment length > 100 nt. Sequences downstream of the end of the alignment were extracted from the transcript of each gene as 3′UTRs. We only kept those 3′UTRs with the length of 25-2500 nts, resulting in a total of 10,212 3′UTRs used for the prediction of miRNA target sites.
We predicted target sites of differentially expressed miRNAs at 3′UTRs using miRanda v3.3a (Enright et al., 2003) and TargetScan v 7.0.2 (Agarwal et al., 2015). In order to obtain more reliable results, we only retained those miRNA-mRNA interactions predicted by both types of software. For miRanda, the entire miRNAs sequences were used as the input and the following parameters were applied: with energy threshold −10 and strict 5′seed pairing. For TargetScan, the seed region (2-8) of miRNAs was extracted and used as the input to predict the miRNA targets with default parameters. We used the conserved seeds for target prediction and retained prediction type of 7mer-1a, 7mer-m8, and 8mer.

| Identification of differentially expressed miRNA targets
With the lists of differentially expressed miRNA and their differentially expressed target genes, we manually extracted miRNA-gene pairs with positively and negatively correlated expression level. As an alternative method, we used Pearson function in MIRLAB (Le et al., 2015) package to test for Pearson correlations of expression level between differentially expressed miRNAs and their differentially expressed target genes. This analysis returns a matrix with the correlation coefficients/scores from −1 to 1, indicating negative (from −1 to 0) and positive (from 0 to 1) miRNA-gene pairs. To reduce the false positives, we only retained the miRNA-gene pairs identified by both methods.

| Functional enrichment analysis
We used Metascape (http://metas cape.org) to perform Gene Ontology (GO) enrichment analysis with the Custom Analysis module (Zhou et al., 2019). A total of 13,409 expressed genes in cochlea were included as the background list. Significantly enriched GO terms were identified using q-value <0.05 adjusted for multiple tests in the Benjamini-Hochberg procedure (Hochberg & Benjamini, 1990). We conducted GO analyses on target genes in negatively correlated miRNA-gene pairs. Redundancy of the GO terms were reduced using the REVIGO clustering algorithm (http:// revigo.irb.hr/). We then used the scatterplots to visualize the clustered GO terms based on their semantic similarities.
In addition, we used a candidate gene approach to identify loci associated with echolocation call frequency variation following Sun et al. (2020). Specifically, we compared the target genes to the list of genes that are known to cause hearing loss and/or deafness in human or mouse (called hearing genes). We collected

| Identification of miRNAs in R. affinis
To annotate miRNAs genes in the R. affinis genome, we used small RNA-seq reads of the cochlea tissue from the eight samples (two subspecies, four samples per subspecies). After removing adaptors and low-quality reads, we obtained an average of 12,160,468 reads per sample with the length of 18-35 nts (Dryad file T3 and

F I G U R E 5
Length distribution of sequenced miRNAs for each sample. Sample information has been described in Table 1 3.2 | One specific miRNA was identified in all echolocating bats Following Jebb et al. (2020), we used Bowtie with the -n 1 parameter to test for the miRNA genes existing in genomes of all echolocating mammals but absent in genomes of nonecholocating mammals. For all 390 miRNA precursors, 10 were present in at least two echolocating bats but absent in nonecholocating mammals (Table 2). More importantly, we found that one miRNA (raf-miR-281) was shared by all echolocating bats (Table 2). A total of 435 target genes were predicted by this miRNA and 16 of them were in the list of hearing genes (Dryad file T5). However, no miRNA was shared by all echolocating mammals (bats and dolphin).

| Identification of differentially expressed miRNAs
In this study, we conducted differential expression analysis on cochlea miRNA-seq data from two subspecies of R. affinis (himalayanus and hainanus), which show over 15 kHz difference of echolocation call frequency (Table 1). Using both DESeq2 and edgeR, we identified 48 differentially expressed miRNAs (DEMs) between the two subspecies (Dryad file T6). To test whether the relative low amount of sequencing reads of Hai-05 and Hai-19 comparing with other samples (Dryad file T3) affected the current results, we repeated the differential expression analysis by reducing the sequencing data of each sample to the one in Hai-05, which resulted in 50 DEMs.
Almost all of these DEMs (45 of 50) were overlapped with the ones based on the whole data above (Dryad file T7), indicating that the retained sequencing reads for each sample were enough to produce reliable results.
We predicted the target genes of all 48 DEMs using two commonly used programs (miRanda and TargetScan). A total of 7741 target genes and 27,067 miRNA-gene relationship pairs were predicted by both programs. Similar to the results in Franchini et al. (2019), almost all retained miRNA-gene pairs were from the predictions of miRanda (98.4%, Figure S5).

| Identification of differentially expressed genes
In this study, we generated the matched mRNA-seq data for each sample with an average of 20,976,159 read pairs 150 bp long per sample, 96.88% of which were retained after quality control (Dryad file T3). A total of 1766 DEGs were identified between the two subspecies using both DESeq2 and edgeR (Dryad file T8). To obtain the Ensembl Gene ID of each DEG, we performed BLAST searches against nt (Nucleotide Sequence Database) and nr (Non-Redundant Protein Sequence Database) with the coding sequences of these DEGs. A total of 1418 DEGs have Ensembl Gene IDs and 1015 of them have 3′UTRs. Among them, 44 were found to be in the list of hearing genes (Dryad file T8). It is to be noted that only 135 DEGs identified in this study were also shown in the list of 799 DEGs identified between the same two subspecies in Sun et al. (2020).
Functional enrichment analysis on these 135 overlapped genes revealed zero significant GO terms.

| Identification of differentially expressed miRNA targets
It is to be noted that as stated above, not all DEGs were annotated with 3′UTRs and those DEGs without 3′UTRs would be missed in the ex- A further GO enrichment analysis on these target DEGs revealed 74 significant GO terms and a majority of GO terms were found to be associated with immunity (Dryad file T11 and Figure 6). TA B L E 2 Ten miRNAs shared in at least two echolocating bats with reference quality genomes (Jebb et al., 2020). Yes and no means presence and absence in the genome of species, respectively. MiRNA shared in all echolocating bats was shown in bold

| DISCUSS ION
As far as we know, this is the first study to characterize miRNAs of the cochlea in an echolocating mammal. The number of mature miRNAs (365) identified here was either more or less than that in other bats (e.g., 468 in Myotis myotis, Huang et al., 2016 and 404 in Myotis lucifugus, Biggar & Storey, 2014;196 in Myotis ricketti, Yuan et al., 2015 and 217 in Eptesicus fuscus, Platt et al., 2014).
This contrast of the number of miRNAs across studies might result from different tissues used because miRNAs have been shown to be expressed with tissue specificity (Lagos-Quintana et al., 2002;Xu et al., 2007). In addition, expression of miRNAs also depends on ecological or environmental context, for example, of torpid and active states in hibernating mammals (Biggar & Storey, 2014Yuan et al., 2015).
The current miRNA-seq data from the cochlea tissue of an echolocating bat (R. affinis) allow us to investigate whether miRNAs play an important role in the evolution of echolocation. We found 10 miRNAs present in at least two echolocating bats but not in other nonecholocating mammals. Interestingly, one of them (raf-miR-281) was shared by all echolocating bats and targeted genes of this miRNA included several hearing genes, which have been proved to be important in echolocation of bats (e.g., KCNQ4, Liu et al., 2012).
Therefore, raf-miR-281 might be a candidate miRNA that is essential in the evolution of echolocation in bats, which can be tested by functional tests of this miRNA and its targeted genes in the future.
However, this candidate miRNA was not present in other echolocating mammals (e.g., dolphin), suggesting independent origins of echolocation in bats and dolphin.
Our current results also support an important role of miRNAs in echolocation call frequency variation between the two subspecies. First, we found that almost half of known hearing genes in the list of all DEGs were regulated negatively by differentially expressed miRNAs (DEMs). Second, target DEGs regulated by DEMs included multiple important hearing genes (Piezo1, Piezo2, and CDH23). PIEZO (PIEZO1/2) has been identified as vertebrate mechanically gated ion channels (Coste et al., 2012) and may play an essential role in mechnotransduction in auditory hair cells . A recent study has shown that PIEZO2 mediates ultrasonic hearing during social behaviors in mice (Li et al., 2021).
In the future, functional assays can be conducted to test for the interactions of these two miRNAs and Piezo1/2 and more importantly to assess their roles in ultrasonic hearing of echolocating mammals. Another one, CDH23, has been shown to be important in the ultrasonic hearing of echolocating mammals (e.g., Cdh23, Shen et al., 2012).
Although it is not one of our main aims, we also compared our results of mRNA differential expression analysis to a previous study that used the same two subspecies as the system (Sun et al., 2020).
To our surprise, we found a low number of shared DEGs between the current study and Sun et al. (2020). This low level of consistency might result from: (1) different individuals sampled; (2) different sampling times; (3) different references used in reads mapping (a high-quality reference genome in this study while a de novo assembled transcriptome in Sun et al., 2020); (4) the high degree of heterogeneity of cell types in the cochlea. Nevertheless, the shared DEGs between these two studies could be considered as the most candidates related to echolocation call frequency variation between the two subspecies.

| CON CLUS ION
In this study, we characterized miRNAs of the cochlea in an echolocating bat. By searching these miRNAs in high-quality genomes of other echolocating and nonecholocating mammals, we found one specific miRNA in all echolocating bats that might be essential in the evolution of echolocation in bats. In addition, our current study adds to a growing number of works that support an important role of miRNAs in phenotypic diversifications of wild animals (e.g., Franchini et al., 2016Franchini et al., , 2019.

AUTH O R CO NTR I B UTI O N S
Qianqian Li: Formal analysis (lead); writing -original draft (equal).

ACK N OWLED G M ENTS
We thank Sun Haijian, Wang Jiaying, and Ding Yuting for assistances with the field data collection. We also thank the associated editor and three anonymous reviewers whose comments significantly improved the manuscript.

CO N FLI C T O F I NTE R E S T
The authors declared that they have no competing interests.