- Research
- Open access
- Published:
Molecular adaptations underlying high-frequency hearing in the brain of CF bats species
BMC Genomics volume 25, Article number: 279 (2024)
Abstract
Background
The majority of bat species have developed remarkable echolocation ability, especially for the laryngeally echolocating bats along with high-frequency hearing. Adaptive evolution has been widely detected for the cochleae in the laryngeally echolocating bats, however, limited understanding for the brain which is the central to echolocation signal processing in the auditory perception system, the laryngeally echolocating bats brain may also undergo adaptive changes.
Result
In order to uncover the molecular adaptations related with high-frequency hearing in the brain of laryngeally echolocating bats, the genes expressed in the brain of Rhinolophus ferrumequinum (CF bat) and Myotis pilosus (FM bat) were both detected and also compared. A total of 346,891 genes were detected and the signal transduction mechanisms were annotated by the most abundant genes, followed by the transcription. In hence, there were 3,088 DEGs were found between the two bat brains, with 1,426 highly expressed in the brain of R. ferrumequinum, which were significantly enriched in the neuron and neurodevelopmental processes. Moreover, we found a key candidate hearing gene, ADCY1, playing an important role in the R. ferrumequinum brain and undergoing adaptive evolution in CF bats.
Conclusions
Our study provides a new insight to the molecular bases of high-frequency hearing in two laryngeally echolocating bats brain and revealed different nervous system activities during auditory perception in the brain of CF bats.
Background
Echolocation is the ability of an animal to determine its surroundings and target position by emitting sound waves and receiving echoes [1]. The echolocation behavior mainly relies on the emission process of ultrasonic signals and the auditory perception process of high-frequency echo acoustic waves [2, 3]. The bats brain can accurately sense surrounding environmental information by analyzing the signal differences between emitted sound waves and echogenic sound waves, including the acquisition of information, such as the distance from obstacles, shape and size, and texture, and can even accurately determine prey distance and flight speed [4,5,6].
With the developments of molecular phylogeny and high throughput sequencing technology, it is helpful for researchers to gain insights to the molecular bases underlying high-frequency hearing of echolocating bats. Typical echolocating bats relying on laryngeal vibrations for vocalization could be divided into two types, FM bats and CF bats. Compared with FM bats, CF bats have a higher dominant frequency of echolocation calls [7], unique CF components [8], highly sensitive structures to the dominant frequency on the basilar membrane of the cochlea [9], and unique Doppler frequency shift compensation behavior [10]. With the development of sequencing technology, the study on the evolution of high-frequency auditory adaptation in echolocation bats has entered the level of omics. Papers demonstrated that the molecular adaptations can occur not only in the coding sequence but also in the regulation of gene expression, where the level of gene expression is closely related to the demand for the protein and its functional necessity [11, 12].
Previously, significantly different expressed genes, along with different physiological processes, and adaptive evolutionary sites of hearing related genes were detected in the cochlea or inner ear of bats with different echolocating types [13]. Nearly, the cochlear structure and molecular bases underlying high-frequency hearing are well studied, however, lacking knowledge about the molecular adaptations of the bats brain [14].
The brain plays an important role in the processing of echolocating signals, and its auditory cortex contains higher-level auditory neural pathways [15]. At the molecular level, it has been found that the Otof gene functioning crucial in acoustic signal transmission in the brain of echolocating bats, and participating in the release of neurotransmitters at the synapse between inner hair cells and the auditory nerve [16]. Especially for the CF bats, developed different auditory nucleus specializations, including enlarged nuclei and cellular differentiation in the cochlear nucleus (CN) [17], distinctly different in the inferior colliculus (IC) region from other types of echolocating bats [18], medial geniculate nucleus (MG) neurons that are sensitive to CF acoustic signals [19], and delay-tuned neurons located in the auditory cortex are specialized [20]. Therefore, we suppose that significant genes and related physiological processes involving with high-frequency hearing could be present in the bat brain, and those genes could be differentially expressed in the brain of CF bats compared with FM bats.
Herein, we performed comparative brain transcriptome for the CF bat (Rhinolophus ferrumequinum, Rhinolophidae, dominant frequency: 74.70 ± 0.13 kHz) and FM bat (Myotis pilosus, Vesperonidae, dominant frequency: 38.21 ± 1.18 kHz) to identify the differentially expressed genes and associated physiological processes. For further, the key hearing related genes were detected and the following adaptive evolution analyses were conducted. This study provides further support for a comprehensive understanding of the widely molecular adaptations in the bat brain, especially to detect key genes and related physiological processes involving with high-frequency hearing in different types of echolocating bats.
Materials and methods
Sampling, RNA extraction and sequencing
The two bat species were both caught from a wild bat colony during September 2021 on the outskirts of Beijing, China (115°59′ N, 39°43′ E). Generally, transcriptome sequencing requires at least three biological repeats to evaluate the reproducibility among individuals within the same species, especially for wild animals. In order to minimize the influence on the wild bat populations and also meet the requirements of comparative transcriptome and statistical analyses, three individuals were collected as three biological repeats for each bat species, that is a total of six bat individuals were used in this study. To avoid any influence of sex-related differences, only males were selected. The brain tissues from each individual were collected and immediately flash-frozen in liquid nitrogen followed by placement in a − 80 °C freezer until processing for total RNA isolation.
Total RNA was isolated using the Total RNA extraction reagent (Trizol) in accordance with the manufacturer’s protocol. The quantity and quality of total RNA were measured using a Qubit 2.0 Fluorometer and gel electrophoresis. RNA samples of the same volume and concentration were used during the step of converting mRNA into cDNA. Three paired-end cDNA libraries of each species were generated using the mRNA-Seq assay. In total, 6 cDNA libraries were prepared at an equimolar ratio for transcriptome sequencing on the Illumina HiSeq TM platform. All raw reads were deposited in the NCBI Short Read Archive (SRA) Database under SRA accession. The raw sequence data generated were deposited into the NCBI Sequence Read Archive database (SRA run accession numbers: R. ferrumequinum: PRJNA1039161, M. pilosus: PRJNA1039163).
Reference transcriptome assembly and functional annotation
The raw reads were filtered by Trimmomatic software [21] using six criteria: removing reads with adaptors; removing reads with unknown “N” bases; removing low-quality bases (Q value < 20) from Reads 3’ to 5’; removing low-quality bases (Q value < 20) from Reads 5’ to 3’; removing bases with quality values below 20 in the tail of Reads (window size of 5 bp) using the sliding window method; removing Reads with length less than 35nt and their paired Reads. To construct a common and powerful reference transcriptome for the comparative analyses, all high-quality raw reads from 6 individual cDNA libraries were used for de novo assembly by Trinity software [22, 23] with the default parameters. After the transcripts were reduced for sequence redundancy, the longest transcript in each cluster was taken as a unigene for further analysis. All of the remaining contigs are described as unigenes in the following text. The assembly result was evaluated by parameters such as the longest value of the sequence, the shortest value of the sequence, the N50 and the N90 values. To obtain the complete information on unigenes functional annotation, NCBI Blast + software [24] and KAAS [25] were used to compare the annotation results with 7 large databases, including NT, NR, COG, PFAM, CDD, GO and KEGG.
Identification of differentially expressed genes (DEGs)
Considering the differences in library size, we first performed inter-sample normalization and used Bowtie2 software [26] to map clean reads to the assembled reference sequence. The expression levels of genes were calculated using the number of reads that were uniquely aligned with the reference transcriptome. Unique mapped reads were quantified into counts for each unigene and the transcripts per million (TPM) [27] formula method was used to determine the unigene expression level. The raw counts for each gene were converted into TPM value using the Salmon software [27] based on the following formula:
TPM was introduced in an attempt to facilitate comparisons across samples [28]. TPM stands for transcript per million, and the sum of all TPM values is the same in all samples, such that a TPM value represents a relative expression level that in principle should be comparable between samples [29]. Here, Xi denotes reads mapped to a transcript; Li is the transcript length, and the \( {{\sum }_{j}\frac{{X}_{j}}{L}}_{j}\) corresponds to the sum of mapped reads to the transcript normalized by the transcript length.
The correlation coefficient between each pair of replicates for each sample was calculated using the R package (corrplot, version 0.92) to evaluate the reliability of the experimental results as well as the operational stability. To determine the separation of expression patterns across samples, principal component analysis (PCA)was performed on the levels of all unigenes using R package (vegan, version 2.6.4) to ensure more reliable results for subsequent analysis.
In consideration of less accurate genes with low expression levels, only unigenes with TPM ≥ 1 were retained for the following analyses. The DEGs between the R. ferrumequinum and M. pilosus were evaluated by DESeq2 (V1.12.4) software [30, 31]. In consideration of the general divergence among bat species and that the detection of DEGs would be more accurate for those with a greater difference in expression, a threshold of fold change ≥ 2 was in this study. Hence, to create a list of high-confidence DEGs for further analyses, the following stringent criteria were used: fold change ≥ 2, namely |log2 (fold change)| ≥ 1, in detail, the P-value was adjusted using the Benjamini-Hochberg method based on an FDR (False Discovery Rate) cut-off of 0.001 [32]. Visualization of the DEGs in the two bat species was achieved by creating a Volcano plot with the R package (ggplot2, version 2.2.1).
GO category and KEGG pathway enrichment analyses
Downstream functional classification was achieved through the integrated localization of GO [33] and KEGG pathway databases [34]. We conducted the GO and KEGG enrichment analysis on the DEGs in OmicShare (http://www.omicshare.com/). All P-values were computed using the hypergeometric test, and multiple test correction was performed using the Benjamini–Hochberg method [35] based on P-value cut-off of 0.01.
To clarify the interactions between genes and key pathways, a directed acyclic graph [36] was subsequently constructed. Following GO and KEGG analysis of DEGs, the directed acyclic graph was plotted using the R package (topGO, version 2.50.0) [37] to visualize the GO terms generated by GO enrichment. GO enrichment is divided into three categories: cellular component, molecular function, and biological process [33]. For the results of KEGG enrichment analysis, function-function, and function-gene interaction network graphs were constructed by association analysis to identify key genes and key functions.
Species coverage
Based on the above analyses, a candidate hearing related gene, ADCY1 was identified as the core connected gene with highly connectivity with the pathways which were significantly enriched by highly expressed genes in the brain of R. ferrumequinum, implying that the ADCY1 may play important roles in the brain of bats. Therefore, we performed following molecular evolution analyses to detect the adaptive evolutionary levels of the ADCY1 gene in CF bats and also other echolocating bats and whales. In detail, the ADCY1 gene sequence (698 bp) was obtained by RNA-Seq sequencing, after blasting in the NCBI database, the complete coding sequences (3,152 bp) of the ADCY1 gene from 27 representative mammals were obtained. The species comprised of two CF bats (R. ferrumequinum and Rhinolophus sinicus), seven FM bats (Molossus molossus, Pipistrellus kuhlii, Myotis myotis, Desmodus rotundus, Sturnira hondurensis, Phyllostomus hastatus, Phyllostomus discolor), two nonecholocating bats (Pteropus giganteus, Pteropus Alecto), seven echolocating toothed whales (Neophocaena asiaeorientalis asiaeorientalis, Phocoena sinus, Physeter catodon, Tursiops truncates, Orcinus orca, Lipotes vexillifer, Delphinapterus leucas), two nonecholocating baleen whales (Balaenoptera acutorostrata scammony, Balaenoptera musculus) and seven other nonecholocating mammals (e.g., Homo sapiens and Mus musculus) (Supplementary Table S1).
Phylogenetic reconstruction
Nucleotide sequences of the 27 species were aligned and spliced in the software ClustalX [38] and Bioedit [39], respectively. Tree construction used Neighbor-Joining (NJ) [40] as implemented in MEGA 11 [41], maximum likelihood (ML) using IQ-TREE [42], and Bayesian inference (BI) in MrBayes 3.2.2 [43]. For both ML and BI trees, the GTR + G nucleotide substitution model was used as the best model which was selected by the jModeltest 2.1.7 software [44], with bootstrap values of 100,000 replicates. For the BI trees, Markov Chain Monte Carlo (MCMC) data simulation was used to estimate the posterior probabilities and performed for 20 million generations with a sampling frequency of 1,000, including a burn-in step corresponding to the first 25% million generations.
Molecular evolutionary analyses
To explore the heterogeneous selection pressures acting on each species, sliding window analyses were performed for ADCY1 using SWAAP 1.0.2 [45]. We estimated the nonsynonymous substitution (dN) and synonymous substitution (dS) substitution rates (the dN/dS ratio, termed omega ω) according to the Nei and Gojobori method [46]. Window size and step size were set as 30 and 6 codons to determine the variations of selection pressure [47] along the ADCY1 gene sequence for two classes of bat species, namely CF bats and FM bats.
The positive selection analysis needs to be performed based on a true species tree, so we determined the true genealogical occurrence relationships between the species used for this gene and mapped the species tree based on previous studies [3, 48, 49]. By comparing ω among sites and branches, the form and intensity of natural selection can be revealed, with ω < 1, ω = 1, and ω > 1 indicating negative selection, neutral evolution, and positive selection, respectively. Various models were analyzed using the CODEML program [50] of the PAML [51] package, including the Site Model (SM), the Branch Model (BM), and the Branch-Site Model (BSM). In the BM and BSM, five foreground branches were set up: CF bats branch, FM bats branch, non-echolocated bats branch, echolocated toothed whale suborder branch, and non-echolocated baleen whale suborder branch, denoted as Branch a, Branch b, Branch c, Branch d, and Branch e, respectively. In the positive selection analysis, the above five branches were analyzed as foreground branches, respectively, while the corresponding remaining branches were used as background branches.
Each Model includes alternative and null hypotheses, and Likelihood ratio tests (LRTs) are performed on the alternative and null hypotheses of the specific model based on the results of the corresponding model runs. The P-value of the LRTs result is 0.05, and if P-value < 0.05, the null hypothesis is rejected and the alternative hypothesis is accepted. LRTs is the test of significance, which means that the sites are considered positively selected sites when P-value < 0.05 and the Bayesian posterior probability is greater than 0.9.
Localization of important sites
Combining the results of previous studies, we predicted and mapped the protein structure of the ADCY1 gene schematically [52, 53]. The expression pattern of ADCY isoforms is mainly obtained from RNA sequencing analysis (at the mRNA level) [52], the ADCY1 gene in this chapter study is a transmembrane protein gene. The protein domains and transmembrane topology of ADCY1 were predicted and plotted according to InterProScan (https://www.ebi.ac.uk/interpro/) and TMHMM Server v. 2.0 (http://www.cbs.dtu.dk/services/TMHMM/). Subsequently, we mapped all the positively selected sites onto the schematic plot of the ADCY1 protein to illustrate its potential changes in CF bats.
Result
Sequencing, assembly, and functional annotation
As shown by Table 1, there were 134,347,358 and 134,347,358 clean reads were obtained for R. ferrumequinum and M. pilosus, respectively (Table 1). The proportions of clean reads among raw reads in each library were 96.46% and 96.51%, respectively suggesting the high quality of the RNA-Seq data available for further analyses. A total of 346,891 unigenes sequences were obtained for the assembled reference transcriptome, ranging from 201 to 15,577 bp in length, with an N50 of 1,460 bp, N90 of 254 bp (Table 1), and the length of all unigenes were calculated (Supplementary Figure S1). After annotation, there were 17,760, 20,034, 14,299, 22,400, 37,912, 87,465, and 196,500 were respectively annotated by the CDD, PFAM, KEGG, KOG, GO, NR, NT databases, and 3,904 unigenes were successfully annotated in all databases (Supplementary Figure S2). For KOG annotation in particular, the term of signal transduction mechanisms was the most highly represented, followed by the transcription. (Supplementary Figure S3).
Brain molecular differences between R. ferrumequinum and M. pilosus.
Results of PCA analysis, TPM density distribution and the correlation analysis of six bats brain samples (Supplementary Figure S4), consistently suggest that the samples show well repeatability and could be used for the downstream analyses. Accordingly, there were 3,088 DEGs were identified in the brains between the two bat species, including 1,426 highly expressed genes in R. ferrumequinum and 1,662 highly expressed genes in M. pilosus (Fig. 1).
DEGs identified between R. ferrumequinum and M. pilosus were significantly enriched in multiple GO terms and KEGG pathways (Supplementary Figure S5). The 1,426 highly expressed genes in the brain of R. ferrumequinum were significantly enriched in the GO terms of biological processes level, a large number of genes were enriched involving neurodevelopmental or cellular morphogenetic processes, such as nervous system development (GO: 0007399), neurogenesis (GO: 0022008), cellular component morphogenesis (GO:0032989), neuron differentiation (GO:0030182), etc. Refer to the 1,662 highly expressed genes in the brain of M. pilosus, they were significantly enriched in the GO terms of regulation of pharyngeal pumping (GO: 0043051) and cell projection assembly (GO: 0030031), pharyngeal pumping (GO:0043050), single-organism process (GO:0044699) etc., which were different with the corresponding results of R. ferrumequinum (Fig. 2). Furthermore, the directed acyclic graph showed that the nervous system development (GO: 0007399), neurogenesis (GO: 0022008), neuron differentiation (GO: 0030182), neuron projection development (GO: 0031175), neuron projection morphogenesis (GO: 0048812) terms were in the same branch (Fig. 3). In the molecular function level, most of the terms which were significantly enriched by the highly expressed genes detected in R. ferrumequinum were related with enzymatic activity, suggesting that more chemical reactions requiring enzymatic catalysis may exist in the brain of R. ferrumequinum. Whereas, for M. pilous, terms were mostly related with channel activity and receptor activity, etc. (Fig. 2). Whereas, different situations were conducted for the GO terms of cellular component, highly expressed genes detected in R. ferrumequinum and M. pilosus were also significantly enriched in similar terms, including neuron part (GO:0097458), synapse part (GO:0044456), axon (GO:0030424), presynapse (GO:0098793) (Fig. 2). The results of the KEGG enrichment analysis showed that R. ferrumequinum was significantly enriched in pathways related to Thyroid hormone synthesis, etc. In addition, both CF and FM bats were enriched to Insulin secretion (ko04911) (Fig. 4). The network interaction analysis conducted by the significantly enriched KEGG pathways showed that Thyroid hormone synthesis (ko04918) was the key pathway in R. ferrumequinum. Moreover, we found the ADCY1 gene, one candidate hearing related gene, was the most highly connected gene in the brain of R. ferrumequinum revealed by the functional-gene network interaction diagram (Fig. 5).
Adaptive evolution of the ADCY1 gene
In total, coding region sequences of the ADCY1 gene for 27 mammals were successfully collected (Supplementary Table S1). Putative gene tree of the ADCY1 gene based on NJ, ML and BI methods showed similar topological structures and high support was obtained at most of the species branch nodes (Fig. 6). In addition, all methods support that CF bats can be clustered together separately, and then clustered together with FM bats and Click bats, which produced a tree topology slightly different from previously reported mammalian species trees (Fig. 6). The Putative gene tree constructed from the ADCY1 gene set consisting of 27 species clustered CF bats together individually, while FM bats and non-echolocated bats clustered at the periphery, and the above clusters had high statistical support. The true species relationship, however, is that CF bats cluster together with non-echolocating bats before clustering with FM bats, demonstrating that the gene evolved differently in CF bats than in other species.
Putative gene tree for ADCY1 coding region sequences using Neighbor-Joining (NJ), maximum likelihood (ML) and Bayesian inference (BI). Values on the branch indicate support from NJ, ML and BI, respectively.
Well-established species tree [3, 48, 49] based on the basis of previous studies for the 27 species. Red amino acid substitutions with ω > 1 are showed. Letters from a–e indicate branches to be tested under Branch model or Branch-site model.
Molecular evolution analyses for the ADCY1 gene
Sliding window analysis showed that the ω values were greater than 1 in the CF bats clade only, demonstrating that the gene suffered stronger selective pressures than the FM bats (Fig. 7). As the result of PAML analyses, Site Model (SM) found six positively selected sites for ADCY1 gene but only one site (242 K 0.905) were significant (supplementary Table S2). The result of the Branch Model (BM) showed that only the ω values of the CF bats branch (ω1) were greater than those of the other remaining branches, implying different evolutionary rates among the CF bats and other mammals (supplementary Table S3). Two positively selected sites (259 I 0.837 and 842 G 0.518) were found by Branch-Site Model (BSM) when CF bats setted as the foreground branch however the P-value of LRT test was not significant (Supplementary Table S4).
The ADCY1 protein consists of 12 transmembrane segments (S1 – S12), the structural domain Guanylate_cyc contains 161 amino acid sites between the transmembrane regions S6 and S7. Then, we mapped the potential significant positively selected sites onto the protein structure and found that the significant positive selection sites 242 K, 259 I, and 477 S were located close to Guanylate_cyc, while the other positive selection sites were all located at the C-terminus of the intramembrane region of ADCY1 protein (Fig. 8).
Red rectangles cover the twelve transmembrane segments (S1 – S12) and the green hexagon represent one structural domain (Guanylate_cyc) between the S6 and S7. Red points with labeled numbers indicate nominally positively selected sites detected by PAML, and the site of 242 K is a significant positively selected site.
Discussion
As was first proposed several decades ago, alterations (or innovations) in gene expression are regarded as an essential means of generating biological diversity [11, 54]. Therefore, analyses of DEGs can reveal the molecular mechanisms underlying phenotypic diversity and provide a deeper understanding of the relationship between gene expression patterns and the resultant phenotypes [55, 56]. We have performed the comparative brain transcriptome between R. ferrumequinum and M. pilosus, using Illumina sequencing technology. A total of 346,891 reference transcriptomes were obtained. In addition, differences in expressed genes and associated physiological processes were found in the two bat brains with different types of echolocating types and two important insights based on downstream analyses were summarized.
Different brain neural activities between CF and FM bats
The results of this study suggest that there may be different neural activity in brain in R. ferrumequinum (CF bat) compared with M. pilosus (FM bat). On the one hand, the results of GO enrichment analysis biological process categories showed that highly expressed genes in brain of R. ferrumequinum were significantly enriched in GO terms directly related with nervous system development. Especially, some of the GO terms enriched by the highly expressed genes in the brain of R. ferrumequinum, such as nervous system development (GO: 0007399), neurogenesis (GO: 0022008), neuron differentiation (GO: 0030182), neuron development (GO: 0031175), neuron projection morphogenesis (GO: 0048812) terms were in the same branch revealed by the directed acyclic graph. On the other hand, Genes highly expressed in the R. ferrumequinum brain were significantly enriched in pathways related to thyroid hormone synthesis (ko 04918) which play a critical role in the differentiation, development of the central nervous system as well as the formation of various functions [57], indicating that the brains of the CF bats may present more active neural activity rather than the FM bats.
Echolocation principal frequency, as an important acoustic parameter, carries a lot of important individual information, and there are large differences in principal frequencies between different types of echolocation bat species [58,59,60]. Bat species with different echolocation types and frequencies have evolved over a long period to be able to perceive well the echolocation acoustic frequencies emitted by themselves [61, 62]. The dominant frequency of the CF bats is significantly higher than the FM bats. The dominant frequencies of CF bat (R. ferrumequinum) and FM bat (M. pilosus) in this study are 74.70 ± 0.13 kHz and 38.21 ± 1.18 kHz, respectively. The perception of echolocation acoustic signals begins in the cochlear hair cells, is continuously distributed on the auditory nerve and terminates in the auditory cortex of the brain [63]. Therefore, we speculate that differences in echolocation principal frequencies is one of the reasons for the differences in gene expression in the brain of different types of echolocating bat species. In addition, there are neural specializations within the auditory brainstem and cortex of echolocating bats [18], such as delay-tuned neurons [64, 65]. In particular, the tonotopic organization of the auditory regions [19], the dorsal and abdominal hypothalamic regions of CF bats are different from other types of echolocating bats and other mammals [66]. The differences described above in CF bats brain are likely the result of adaptive changes to adapt to their unique high-frequency CF component in the echolocation calls. Compared with FM bats, genes highly expressed in the CF bats brain mostly involving with nervous system development related biological processes, which is likely to be the genetic basis of those specific structures and functions in CF bats brain. It has been shown that the auditory cortex of CF bats has some well-developed physiology and organization to respond to specific neurons compared to the auditory cortex of other mammals [19], since then, we suppose that there may be specific types of neurons in CF bats’ brain, and the above highly expressed genes and associated biological functions may be a direct reflection at the molecular level.
Besides, although the highly expressed genes were different in the CF bats compared with FM bats, the same pathway, Insulin secretion (ko04911), was also detected to be significantly enriched by DEGs detected in CF and FM bats, respectively, indicating there were similar biological processes and molecular mechanisms in the brains for recognizing CF and FM bats. Insulin acts on the CNS to modulate behavior and systemic metabolism [67]. Recently, the WFS1 gene expressed in the brain involving with insulin secretion has attracted extensive attention [68], and the heterozygous mutations in the WFS1 gene may cause deafness [69], hence this gene and its associated pathways, like insulin secretion, may be a possible mechanism for different hearing development.
Adaptive evolution of hearing-related gene ADCY1
Previous studies have demonstrated that the ADCY1 gene is a member of the adenylate cyclase gene family and primarily expressed in the brain, in details, with high expression levels in neurons and moderate in oligodendrocytes, microglia and astrocytes [52, 70]. The ADCY1 gene was also detected to express throughout inner ear development and maturation in mouse [71]. In this study, The ADCY1 gene was identified to be the core connected gene with highly connectivity with the pathways which were significantly enriched by highly expressed genes in the brain of R. ferrumequinum, indicating its important functions in the auditory process in echolocating bat. Therefore, following adaptive evolutionary analyses of the ADCY1 gene were performed to further reveal its potential changes of functional sites along the sequence in echolocating bats. A multi-model positive selection analysis of the ADCY1 gene identified several positive sites with highly significant 242 K 0.905, and the remaining sites were detected with posterior probabilities not reaching 0.9, suggesting that these sites may be evolutionarily important sites for the ADCY1 gene, but have not yet reached significant positive selection at this stage. Not many sites were detected in the positive selection analysis, indicating that the gene is relatively conserved during evolution in CF bats, but functional gene sequences with important roles are usually conserved to ensure the stability of gene function in genetic evolution [72,73,74].
Combined with the results of the sliding window analysis (Fig. 7), the ADCY1 gene was subjected to inconsistent selection pressure in the two species, with ω>1 occurring only in the CF bats category, demonstrating that a locus in this part of the bat may have been subjected to positive selection in CF bats and that most of the positively selected sites were at locations where non-synonymous and synonymous substitutions of ADCY1 amino acids occurred frequently. Previously, hearing-related genes, sk2 and shh were also detected by the sliding window analyses hand with stronger selection pressure in echolocation bats, and their potential positive selection sites were found in the CF bats branch, but both genes were evolutionarily relatively conserved [75].
By the results of structural and molecular evolutionary analysis of ADCY1 protein, we found that all positive selection sites were distributed in parts outside the functional domain. The protein structure-function prediction revealed that most of the sites subject to positive selection are located near the structural domain or at the C-terminus. The main reason may be that since the key functional domains of channel proteins are responsible for exercising important transport functions, the protein sequences of these important functional domains have evolved more conservatively across species to ensure proper ion transport and signaling, while the positive selection sites detected in the non-important functional domain fraction may be the sites that exert an influence on the evolution of high-frequency hearing in CF bats. The structural domains of genes may be more conserved in gene sequence to ensure the stability of gene function during inheritance due to their decisive role in the exercise and maintenance of the overall gene function [74].
All of these results suggest that the auditory-related gene ADCY1 has undergone adaptive evolution in CF bats species. Considering that CF bats have higher vocalization frequencies and auditory abilities compared to other echolocation species, it is hypothesized that potential positive selection sites detected on the ADCY1 gene may play an important role in the acquisition and development of their echolocation abilities or high-frequency auditory abilities. Although ADCY1 is evolutionarily conserved at the sequence level, compared with M. pilosus, it is highly expressed in R. ferrumequinum brain. Therefore, we suspect that ADCY1 may be act by increasing the expression level, rather than through the changes at the sequence level for now.
In addition, we found the auditory gene Otof [16] and the vocal gene FOXP2 [76] were also highly expressed in CF bats, although not significantly so. The above results demonstrate the presence of key genes related to hearing in the brain of CF bats, which may be responsible for the ability of CF bats to receive higher-frequency sound waves. However, the well studied hearing genes, Prestin [12, 77], KCNQ4 [72, 78] and SK2 [75, 79], which were previously reported to have undergo adaptive evolution in the laryngeal echolocation bats cochlea, were not detected in these two bats brains. Together, we speculated different molecular mechanisms may exist in the brain and cochlea in sensing high-frequency sound waves, reflecting by expressed different hearing related genes, different expression modes, and different adaptive sites.
The auditory-related gene ADCY1 in this study is the first to be identified in the comparative transcriptome of the bat brain, complementing previous findings in the auditory system. We have found the ADCY1 gene in the brain, which enriches the previous adaptive evolution research carried out by auditory genes.
Conclusions
Neuron, synapse parts and other related genes were detected to be widely expressed in the brains of two typical laryngeally echolocating bat species. In particular, compared with FM bats, much more highly expressed genes were found to significantly enriched in neurodevelopment and thyroid hormone synthesis suggesting more active nervous activity in the brain of CF bats. Similar results were also demonstrated by the study of the cochleae between CF bats and FM bats. Moreover, a key gene, ADCY1 was identified to be the most highly connected gene, and downstream analyses indicated that the ADCY1 gene has suffered more powerful selective pressures and may act important functions in the brain of CF bats. This study provides further support for a comprehensive understanding of the expressed genes in the brain involved in bats echolocation with a view to uncover the molecular bases of high-frequency hearing in echolocating bats. In future, further analyses and correspondingly functional verification experiments would be performed on other types of bats to explore the accurately molecular mechanisms of high-frequency hearing in laryngeally echolocating bats.
Data availability
Raw data are available at SRA, under the accession number: PRJNA1039161 and PRJNA1039163.
Abbreviations
- DEGs:
-
differentially expressed genes
- CN:
-
cochlear nucleus
- IC:
-
inferior colliculus
- MG:
-
medial geniculate
- LRTs:
-
Likelihood ratio tests
References
Jones G, Echolocation. Curr Biol. 2005;15(13):R484–8.
Moss CF, Surlykke A. Auditory scene analysis by echolocation in bats. J Acoust Soc Am. 2001;110(4):2207–26.
Jones G, Teeling EC. The evolution of echolocation in bats. Trends Ecol Evol. 2006;21(3):149–56.
Kunz TH. Ecology of bats. Springer Science & Business Media; 2013.
Griffin DR. Listening in the dark: the acoustic orientation of bats and men. 1958.
Thomas JA, Moss CF, Vater M. Echolocation in bats and dolphins. University of Chicago Press; 2004.
Suga N. Principles of auditory information-processing derived from neuroethology. J Exp Biol. 1989;146:277–86.
Neuweiler G. Evolutionary aspects of bat echolocation. J Comp Physiol Neuroethol Sens Neural Behav Physiol. 2003;189(4):245–56.
Davies KT, Maryanto I, Rossiter SJ. Evolutionary origins of ultrasonic hearing and laryngeal echolocation in bats inferred from morphological analyses of the inner ear. Front Zool. 2013;10(1):2.
Russell IJ, Kössl M. Micromechanical responses to tones in the auditory fovea of the greater mustached bat’s cochlea. J Neurophysiol. 1999;82(2):676–86.
Carroll SB. Evolution at two levels: on genes and form. PLoS Biol. 2005;3(7):e245.
Li Y, Liu Z, Shi P, Zhang J. The hearing gene prestin unites echolocating bats and whales. Curr Biol. 2010;20(2):R55–56.
Wang H, Zhao H, Huang X, Sun K, Feng J. Comparative cochlear transcriptomics of echolocating bats provides new insights into different nervous activities of CF bat species. Sci Rep. 2018;8(1):15934.
Sulser RB, Patterson BD, Urban DJ, Neander AI, Luo ZX. Evolution of inner ear neuroanatomy of bats and implications for echolocation. Nature. 2022;602(7897):449–54.
Zeng Y. A study of th molecular mechanisms of adaptive evolution in Chiroptera by transcriptomes analyses. East China Normal University; 2016.
Shen YY, Liang L, Li GS, Murphy RW, Zhang YP. Parallel evolution of auditory genes for echolocation in bats and toothed whales. PLoS Genet. 2012;8(6):e1002788.
Zook J, Casseday J. Cytoarchitecture of auditory system in lower brainstem of the mustache bat, Pteronotus parnellii. J Comp Neurol. 1982;207(1):1–13.
Schuller G, Pollak G. Disproportionate frequency representation in the inferior colliculus of Doppler-compensating greater horseshoe bats: evidence for an acoustic fovea. J Comp Physiol. 1979;132:47–54.
Wenstrup JJ. Frequency organization and responses to complex sounds in the medial geniculate body of the mustached bat. J Neurophysiol. 1999;82(5):2528–44.
Butman JA, Suga N. Inhibitory mechanisms shaping delay-tuned combination-sensitivity in the auditory cortex and thalamus of the mustached bat. Hear Res. 2019;373:71–84.
Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30(15):2114–20.
Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, Bowden J, Couger MB, Eccles D, Li B, Lieber M. De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. Nat Protoc. 2013;8(8):1494–512.
Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, Adiconis X, Fan L, Raychowdhury R, Zeng Q. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29(7):644–52.
Altschul SF, Madden TL, Schäffer AA, Zhang J, Zhang Z, Miller W, Lipman DJ. Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. Nucleic Acids Res. 1997;25(17):3389–402.
Moriya Y, Itoh M, Okuda S, Yoshizawa AC, Kanehisa M. KAAS: an automatic genome annotation and pathway reconstruction server. Nucleic Acids Res. 2007;35(suppl2):W182–5.
Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods. 2012;9(4):357–9.
Patro R, Duggal G, Love MI, Irizarry RA, Kingsford C. Salmon provides fast and bias-aware quantification of transcript expression. Nat Methods. 2017;14(4):417–9.
Zhao Y, Li MC, Konaté MM, Chen L, Das B, Karlovich C, Williams PM, Evrard YA, Doroshow JH, McShane LM. TPM, FPKM, or normalized counts? A comparative study of quantification measures for the analysis of RNA-seq data from the NCI patient-derived models repository. J Transl Med. 2021;19(1):269.
Vera Alvarez R, Pongor LS, Mariño-Ramírez L, Landsman D. TPMCalculator: one-step software to quantify mRNA abundance of genomic features. Bioinformatics. 2019;35(11):1960–2.
Love MI, Huber W, Anders S. Moderated estimation of Fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):1–21.
Fröhlich A, Pfaff AL, Bubb VJ, Quinn JP, Koks S. Transcriptomic profiling of cerebrospinal fluid identifies ALS pathway enrichment and RNA biomarkers in MND individuals. Exp Biol Med (Maywood). 2023;248(23):2325–31.
Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J Roy Stat Soc: Ser B (Methodol). 1995;57(1):289–300.
Harris MA, Clark J, Ireland A, Lomax J, Ashburner M, Foulger R, Eilbeck K, Lewis S, Marshall B, Mungall C et al. The gene ontology (GO) database and informatics resource. Nucleic Acids Res. 2004;32:D258–261.
Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30.
Ferreira JA. The Benjamini-Hochberg method in the case of discrete test statistics. Int J Biostat 2007, 3(1):Article 11.
Han SW, Chen G, Cheon M-S, Zhong H. Estimation of directed acyclic graphs through two-stage adaptive lasso for gene network inference. J Am Stat Assoc. 2016;111(515):1004–19.
Alexa A, Rahnenfuhrer J. topGO: enrichment analysis for gene ontology. R Package Version. 2010;2(0):2010.
Larkin MA, Blackshields G, Brown NP, Chenna R, McGettigan PA, McWilliam H, Valentin F, Wallace IM, Wilm A, Lopez R, et al. Clustal W and Clustal X version 2.0. Bioinformatics. 2007;23(21):2947–8.
Hall T, Biosciences I, Carlsbad C. BioEdit: an important software for molecular biology. GERF Bull Biosci. 2011;2(1):60–1.
Saitou N, Nei M. The neighbor-joining method: a new method for reconstructing phylogenetic trees. Mol Biol Evol. 1987;4(4):406–25.
Tamura K, Stecher G, Kumar S. MEGA11: molecular evolutionary genetics analysis version 11. Mol Biol Evol. 2021;38(7):3022–7.
Nguyen L-T, Schmidt HA, Von Haeseler A, Minh BQ. IQ-TREE: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. Mol Biol Evol. 2015;32(1):268–74.
Huelsenbeck JP, Ronquist F. MRBAYES: bayesian inference of phylogenetic trees. Bioinformatics. 2001;17(8):754–5.
Darriba D, Taboada GL, Doallo R, Posada D. jModelTest 2: more models, new heuristics and parallel computing. Nat Methods. 2012;9(8):772–772.
Pride D. SWAPP 1.0. 2: a tool for analyzing substitutions and similarity in multiple alignments. See http://www bacteriamuseum org/SWAAP/SwaapPage htm 2004.
Nei M, Gojobori T. Simple methods for estimating the numbers of synonymous and nonsynonymous nucleotide substitutions. Mol Biol Evol. 1986;3(5):418–26.
Pride D. A tool for analyzing substitutions and similarity in multiple alignments. Distributed by the author 2005.
Lack JB, Roehrs ZP, Stanley CE Jr, Ruedi M, Van Den Bussche RA. Molecular phylogenetics of Myotis indicate familial-level divergence for the genus Cistugo (Chiroptera). J Mammal. 2010;91(4):976–92.
Teeling E, Dool S, Springer M. Phylogenies, fossils and functional genes: the evolution of echolocation in bats. Evolutionary History bats: Fossils Molecules Morphology 2012:1–22.
Gao F, Chen C, Arab DA, Du Z, He Y, Ho SY. EasyCodeML: a visual tool for analysis of selection using CodeML. Ecol Evol. 2019;9(7):3891–8.
Yang Z. PAML 4: phylogenetic analysis by maximum likelihood. Mol Biol Evol. 2007;24(8):1586–91.
Devasani K, Yao Y. Expression and functions of adenylyl cyclases in the CNS. Fluids Barriers CNS. 2022;19(1):23.
Nicol X, Muzerelle A, Bachy I, Ravary A, Gaspar P. Spatiotemporal localization of the calcium-stimulated adenylate cyclases, AC1 and AC8, during mouse brain development. J Comp Neurol. 2005;486(3):281–94.
Dittmar K, Liberles D. Evolution after gene duplication. Wiley; 2011.
Perry GH, Melsted P, Marioni JC, Wang Y, Bainer R, Pickrell JK, Michelini K, Zehr S, Yoder AD, Stephens M, et al. Comparative RNA sequencing reveals substantial genetic variation in endangered primates. Genome Res. 2012;22(4):602–10.
Todd EV, Black MA, Gemmell NJ. The power and promise of RNA-seq in ecology and evolution. Mol Ecol. 2016;25(6):1224–41.
Billon N, Jolicoeur C, Tokumoto Y, Vennström B, Raff M. Normal timing of oligodendrocyte development depends on thyroid hormone receptor alpha 1 (TRα1). EMBO J. 2002;21(23):6452–60.
Siemers BM, Beedholm K, Dietz C, Dietz I, Ivanova T. Is species identity, sex, age or individual quality conveyed by echolocation call frequency in European horseshoe bats? Acta Chiropterologica. 2005;7(2):259–74.
Metzner W, Zhang S, Smotherman M. Doppler-shift compensation behavior in horseshoe bats revisited: auditory feedback controls both a decrease and an increase in call frequency. J Exp Biol. 2002;205(11):1607–16.
Chen S-F, Jones G, Rossiter SJ. Determinants of echolocation call frequency variation in the Formosan lesser horseshoe bat (Rhinolophus monoceros). Proceedings of the Royal Society B: Biological Sciences 2009, 276(1674):3901–3909.
Hiryu S, Mora EC, Riquimaroux H. Behavioral and physiological bases for doppler shift compensation by echolocating bats. Bat Bioacoustics 2016:239–63.
Zhang J-S, Han N-J, Jones G, Lin L-K, Zhang J-P, Zhu G-J, Huang D-W, Zhang S-Y. A new species of Barbastella (Chiroptera: Vespertilionidae) from north China. J Mammal. 2007;88(6):1393–403.
Suga N, Neuweiler G, Möller J. Peripheral auditory tuning for fine frequency analysis by the CF-FM bat, Rhinolophus ferrumequinum: IV. Properties of peripheral auditory neurons. J Comp Physiol. 1976;106(1):111–25.
Simmons JA. The resolution of target range by echolocating bats. J Acoust Soc Am. 1973;54(1):157–73.
Neuweiler G, Bruns V, Schuller G. Ears adapted for the detection of motion, or how echolocating bats have exploited the capacities of the mammalian auditory system. J Acoust Soc Am. 1980;68(3):741–53.
Suga N, Jen PH. Disproportionate tonotopic representation for processing CF-FM sonar signals in the mustache bat auditory cortex. Science. 1976;194(4264):542–4.
Kullmann S, Kleinridders A, Small DM, Fritsche A, Häring HU, Preissl H, Heni M. Central nervous pathways of insulin action in the control of metabolism and food intake. Lancet Diabetes Endocrinol. 2020;8(6):524–34.
Kõks S. Genomics of Wolfram Syndrome 1 (WFS1). Biomolecules 2023, 13(9).
Cryns K, Sivakumaran TA, Van den Ouweland JM, Pennings RJ, Cremers CW, Flothmann K, Young TL, Smith RJ, Lesperance MM, Van Camp G. Mutational spectrum of the WFS1 gene in Wolfram syndrome, nonsyndromic hearing impairment, diabetes mellitus, and psychiatric disease. Hum Mutat. 2003;22(4):275–87.
Brain RNA-S. https://brainrnaseq.org/. Accessed 29 Nov 2021.
Santos-Cortez RLP, Lee K, Giese AP, Ansar M, Amin-Ud-Din M, Rehn K, Wang X, Aziz A, Chiu I, Hussain Ali R. Adenylate cyclase 1 (ADCY1) mutations cause recessive hearing impairment in humans and defects in hair cell function and hearing in zebrafish. Hum Mol Genet. 2014;23(12):3289–98.
Liu Y, Han N, Franchini LF, Xu H, Pisciottano F, Elgoyhen AB, Rajan KE, Zhang S. The voltage-gated potassium channel subfamily KQT member 4 (KCNQ4) displays parallel evolution in echolocating bats. Mol Biol Evol. 2012;29(5):1441–50.
van Ooyen A, Kwee V, Nusse R. The nucleotide sequence of the human int-1 mammary oncogene; evolutionary conservation of coding and non-coding sequences. Embo j. 1985;4(11):2905–9.
Dermitzakis ET, Clark AG. Evolution of transcription factor binding sites in mammalian gene regulatory regions: conservation and turnover. Mol Biol Evol. 2002;19(7):1114–21.
Wang H, Zhao H, Chu Y, Feng J, Sun K. Assessing evidence for adaptive evolution in two hearing-related genes important for high-frequency hearing in echolocating mammals. G3 (Bethesda) 2021, 11(4).
Haesler S, Rochefort C, Georgi B, Licznerski P, Osten P, Scharff C. Incomplete and inaccurate vocal imitation after knockdown of FoxP2 in songbird basal ganglia nucleus area X. PLoS Biol. 2007;5(12):e321.
Li G, Wang J, Rossiter SJ, Jones G, Cotton JA, Zhang S. The hearing gene prestin reunites echolocating bats. Proc Natl Acad Sci U S A. 2008;105(37):13959–64.
Liu Z, Li S, Wang W, Xu D, Murphy RW, Shi P. Parallel evolution of KCNQ4 in echolocating bats. PLoS ONE. 2011;6(10):e26618.
Nie L, Song H, Chen M-F, Chiamvimonvat N, Beisel KW, Yamoah EN, Vázquez AE. Cloning and expression of a small-conductance Ca2+-activated K + channel from the mouse cochlea: coexpression with α9/α10 acetylcholine receptors. J Neurophysiol. 2004;91(4):1536–44.
Acknowledgements
We thank the faculty and students of the Key Laboratory of Jilin Province for the conservation and utilization of animal resources and valuable assistance in the field. We would like to express our gratitude to our colleagues who contributed to this article and to the faculty members who reviewed the article.
Funding
This research was funded by the National Natural Science Foundation of China (grant number 32371561 and 32001089 to HW and grant number 32071492 to JF).
Author information
Authors and Affiliations
Contributions
X.L. and H.W. designed the experiments, analyzed and interpreted the data and organized the manuscript; X.W., M.B., R.S. and W.D. provided great help to the sample collection; H.W. provides supervision during the entire project. K.S. and J.F. assisted with the study design and provided laboratory space and reagents. All authors reviewed the manuscript.
Corresponding authors
Ethics declarations
Ethics approval and consent to participate
Ethics approval and consent to participate According to the regulations of Wildlife Conservation of the People’s Republic of China (Chairman Decree (2016) No. 47), permits are required only for species included in the list of state-protected and region-protected wildlife species. None of the bats used in this study are endangered or region-protected species, so no specific permission is required. All animal experimental procedures were approved by the National Animal Research Authority of Northeast Normal University, China (approval number: Nenu-20080416), and the Forestry Bureau of Jilin Province, China (approval number: [2006] 178). All efforts were made to minimize the suffering of the animals. We confirm that all methods were performed in accordance with the relevant guidelines and regulations mentioned above. We declare that this study is reported in accordance with ARRIVE guidelines.
Competing interests
The authors declare that they have no competing interests.
Additional information
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Electronic supplementary material
Below is the link to the electronic supplementary material.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Li, X., Wang, H., Wang, X. et al. Molecular adaptations underlying high-frequency hearing in the brain of CF bats species. BMC Genomics 25, 279 (2024). https://doi.org/10.1186/s12864-024-10212-6
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s12864-024-10212-6