Behavioral innovation and genomic novelty are associated with the exploitation of a challenging dietary opportunity by an avivorous bat.

Foraging on nocturnally migrating birds is one of the most challenging foraging tasks in the animal kingdom. Only three bat species (e.g., Ia io) known to date can prey on migratory birds. However, how these bats have exploited this challenging dietary niche remains unknown. Here, we demonstrate that I. io hunts at the altitude of migrating birds during the bird migration season. The foraging I. io exhibited high flight altitudes (up to 4945 m above sea level) and high flight speeds (up to 143.7 km h-1). I. io in flight can actively prey on birds in the night sky via echolocation cues. Genes associated with DNA damage repair, hypoxia adaptation, biting and mastication, and digestion and metabolism have evolved to adapt to this species' avivorous habits. Our results suggest that the evolution of behavioral innovation and genomic novelty are associated with the exploitation of challenging dietary opportunities.

Introduction

How ecological niche breadth evolves is central to biological adaptation and has long been a topic in evolutionary ecology (Sexton et al., 2017). Niche expansion can improve the ability of organisms to cope with environmental fluctuations and improve their fitness (McCormack and Smith, 2008; Pagani-Núñez et al., 2019; Weng et al., 2005). However, niche expansion is often assumed to be limited by genetics and phenotypic factors (e.g., morphology, behavior, and physiology) (McCormack and Smith, 2008; Sexton et al., 2017). Thus, studying how organisms break through these limitations is crucial to understanding how niche breadth evolves, but the behavioral mechanisms and genetic architecture underlying the evolution of niche breadth remain unclear.

The migration of birds serves as an ecological dietary opportunity for many animal predators (Canário et al., 2012; Dierschke, 2003). At least 40% of the world’s 10,000 bird species are migratory birds (El-Sayed, 2019). To take advantage of the favorable weather conditions at night and to reduce the risk of daytime predation, billions of normally diurnal songbirds switch to migrating at night and fly across the continents and oceans on their annual migratory journeys (Kerlinger and Moore, 1989; Schmaljohann et al., 2007; Sjöberg et al., 2021). Foraging on nocturnally migrating birds is one of the most challenging foraging tasks in the animal kingdom, as it requires predators to overcome behavioral, sensory, and physiological challenges such as high-altitude flight, low visibility, and low oxygen. Preying on nocturnally migrating birds is impossible for diurnal animals or for animals without the capacity for powered flight.

Only three bat species—Nyctalus lasiopterus, Nyctalus aviator, and Ia io—are known to prey on migratory birds (Gual-Suárez and Medellín, 2021). Dietary analysis based on DNA metabarcoding suggests that N. lasiopterus, N. aviator, and I. io prey on 31, 14, and 22 species of passerine birds, respectively (Gong et al., 2021; Ibáñez et al., 2016, 2020). Avivory in bats is rare in nature. Avivorous bats have experienced a dietary niche expansion from insects to birds (Popa-Lisseanu et al., 2007). However, it remains uncertain whether avivorous bats forage on birds by the aerial-hawking or near-ground gleaning strategies using echolocation and/or other sensory cues (e.g., passive auditory and visual cues). Foraging on birds via the aerial-hawking strategy is more challenging in terms of energy consumption and predation risk than the ground-gleaning strategy. The calls of nocturnal birds and light conditions at night may be used by bats as auditory and visual cues for predation. Moreover, it is unclear how avivorous bats exploit this challenging dietary niche. Our previous study found that a population of I. io in China feeds on insects in summer and on birds during the bird migration season (Gong et al., 2021). Here, we used behavioral, sensory, and genomic analyses to reveal how I. io exploits this challenging dietary ecological opportunity.

Results and discussion

Foraging spatial altitude and flight speed

To investigate whether I. io can prey on migrating birds flying at high altitudes using an aerial-hawking strategy, we fitted I. io with Global Position System (GPS) loggers that recorded their nocturnal foraging spatial locations (three-dimensional locations). We predicted that I. io in the bird migration season would hunt at high altitudes where nocturnally migrating birds normally fly. The GPS recording was performed in Xingyi, Guizhou Province, China, which is located on a bird migration route and is within the range of the East Asian-Australasian flyway (EAAF) (UNEP/CMS, 2012). I. io in this area prey almost exclusively on insects in summer but on a large number of bird species in autumn (Gong et al., 2021). We obtained 1678 spatial location data records of 13 individuals foraging in summer, and 521 spatial location data records of 7 individuals foraging in autumn.

I. io exhibited different foraging spatial activity patterns between the bird-eating season and the insect-eating season. The bats hunt in the same areas in summer and autumn (Figure 1A). However, the bats’ flight altitudes above sea level (asl) and above ground level (agl) were significantly higher in autumn than in summer (Mann–Whitney U test, for asl: Z = −22.343, p < 0.001; for agl, Z = −26.478, p < 0.001; Figures 1A and 1B). The flight altitude of I. io during the autumn night activities was 2068.03 ± 563.90 m asl (mean ± SD for flight altitude) and 493.94 ± 393.51 m agl (Figures 1B and 1E), values that were basically in accordance with the flight altitude of migrating small passerines (about 2000 m asl and 200–800 m agl) (Bruderer et al., 2018; Kerlinger and Moore, 1989). In summer, I. io foraged at a flight altitude of 93.91 ± 85.06 m agl (1489.76 ± 181.25 m asl; Figures 1B and 1D), which was consistent with the flight altitude of most insectivorous bats that prey on insects that are largely concentrated below 100 m agl [e.g., Nyctalus noctula (Roeleke et al., 2016) and Rhinolophus ferrumequinum (Fujioka et al., 2020)]. The high flight altitude of I. io foraging in autumn suggests that they prey on nocturnally migrating birds rather than on near-ground overnighting birds. However, the fact that bats fly high does not necessarily mean that they are hunting birds. Some bats have been shown to fly high in search of migratory insects such as Tadarida brasiliensis (McCracken et al., 2008). Thus, further studies should be conducted during the bird migration season (e.g., autumn) to confirm that avivorous bats fly high to prey on nocturnally migrating birds.

Figure 1

The flight altitude and speed of foraging activity in Ia io

(A) Tracking site and 3D flight site data for bat individuals in summer and autumn. The underlying digital elevation model (DEM) was derived from 12.5 m National Aeronautics and Space Administration imagery; values ranged from 710 m asl to 2173 m asl.

(B) Comparisons of the bats’ flight altitude above sea level (asl) and above ground level (agl) between summer and autumn. In the boxplot, the lower and upper edges of a box represent the 25% (q1) and 75% quartiles (q3), respectively. The horizontal line within a box indicates the median (md). The whiskers extend to the most extreme values within inner borders, md ± 1.5 (q3-q1). Dots represent the raw data points. p values are from Mann–Whitney U tests. ∗∗∗p < 0.001.

(C) The maximum flight altitude above ground level for seven bat species [Rousettus aegyptiacus (Tsoar et al., 2011), Tadarida teniotis (O'Mara et al., 2021), Taphozous theobaldi (Roeleke et al., 2018), Rhinophylla microphyllum (Cvikel et al., 2015), Nyctalus noctula (Roeleke et al., 2016), Rhinolophus ferrumequinum (Fujioka et al., 2020), and Nyctalus lasiopterus (Naďo et al., 2019)] recorded by GPS loggers as well as the maximum flight altitude above ground level of I. io in summer and autumn.

(D–G) The kernel density distribution of flight altitudes above sea level and flight speeds for all bat individuals. Each individual bat’s kernel density was normalized to sum to one. The kernel density distribution of individual flight altitude above sea level of I. io in summer (D) and in autumn (E). The kernel density distribution of individual flight speeds of I. io in summer (F) and in autumn (G).

I. io exhibited flight behavioral innovation to adapt to preying on nocturnally migrating birds. First, I. io exhibited high flight altitude associated with foraging on migrating birds. The flight altitude of nocturnally migrating songbirds can exceed 2000 m asl, reaching as high as 5000 m asl (Sjöberg et al., 2021). The maximum foraging flight altitudes of insectivorous bats and an avivorous bat (N. lasiopterus) recorded by GPS loggers were below 1700 m agl (Figure 1C). We recorded three bat individuals that flew at altitudes over 4500 m asl in autumn, and the maximum flight altitude was 4945 m asl (3301 m agl; Figures 1B, 1C, and 1E). Thus, the high-altitude flight behavior of I. io in autumn enables them to meet the spatial accessibility requirements of preying on migratory birds flying at high altitudes. Second, I. io exhibited a high flight speed adapted to hunting migrating birds. The migratory flight speed of most songbirds was about 32–48 km h−1 (El-Sayed, 2019), which is faster than the flight speed of most known insectivorous bats (14–28 km h−1) (Fujioka et al., 2020; Roeleke et al., 2016; Voigt et al., 2019). We recorded the flight speed of I. io and found that its speeds ranged from 1.7 to 65.8 km h−1 (mean ± SD: 21.02 ± 9.06 km h−1) in summer and from 1.7 to 143.7 km h−1 (21.06 ± 11.17 km h−1) in autumn (Figures 1F and 1G); speeds that can meet the requirements for preying on migrating birds. Although no significant difference was observed in the average flight speed between summer and autumn (Mann–Whitney U test, Z = 0.616, p = 0.538), the standard deviation and maximum flight speed of I. io were larger in autumn than in summer, suggesting that flight speed functioned to meet the requirements for chasing fast-flying birds. Taken together, the results suggest that I. io exhibits flight behavioral innovations that may allow them to overcome the challenges of flight altitude and speed to prey on nocturnally migrating birds.

Sensory cues for preying on birds

In addition to the constraint of high-altitude flight, predators preying on nocturnally migrating birds are also subject to the sensory constraint of low visibility. Many nocturnally migrating songbirds produce flight calls in the frequency range of 1–10 kHz; these calls have the function of maintaining flock contact and stimulating migratory restlessness during migration (Farnsworth, 2005; Hamilton, 1962). The calls of birds and light conditions at night may be used by predators as auditory and visual cues for predation. As nocturnal predators, bats can locate their target prey through an active echolocation system, passive hearing, or vision (Fenton et al., 2016; Jacobs and Bastian, 2016). We performed behavioral experiments to confirm the sensory cues that I. io uses to prey on flying birds at night.

We investigated whether I. io can actively prey on birds in flight under the cover of darkness. Our previous study showed that Phylloscopus inornatus and Zosterops japonicus had the highest and medium frequencies of occurrence in the diet of I. io, respectively (Gong et al., 2021). Thus, we hung one P. inornatus and two Z. japonicus bird specimens in a completely dark flight room (see STAR Methods and Figure S1). We placed bats in the flight room and recorded their foraging behavior and echolocation calls. All 11 tested bats actively attack the hanging bird specimens (for example, see Video S1). The number of attacks ranged from 1 to 6 times (Table S1). The bats emitted echolocation calls when they approached and attacked the bird specimens. The echolocation calls of the bats that directly attacked the specimens included a feeding buzz (defined as the pulse repetition rate being increased) (Figure 2A). Echolocating bats use feeding buzzes only for prey capture (Simmons et al., 1979; Surlykke et al., 2003). Our results suggested that I. io can actively locate and prey on birds in flight by using echolocation cues.

Figure 2

Active aerial flight and sensory cues of Ia io preying on birds

(A) The echolocation call sequences of bats attacking bird specimens and the corresponding reconstructed three-dimensional flight trajectory. Waveform and spectrogram showing approach phase and final feeding buzz. Top views of the experimental setup for auditory cues (B) and visual cues (C). The dual-choice experimental flight cage was divided into four selection chambers (Chamber 1, 2, 3, and 4) at both ends by a sliding baffle. Playback files of silence and bird calls, bird specimens, and LED table lamps were alternately presented at each end of the flight cage.

(D) For auditory cues, no significant differences were found in the percentage of bats’ choices in the experimental groups of silence versus bird calls and specimen + silence versus specimen + bird calls.

(E) The low-frequency hearing of one of the bats (range 4–7 kHz). Bats did not demonstrate the typical ABR waveform curve when the frequency of sound playback was lower than 6 kHz. All three individual bats’ ABR waveform curves are shown in Figure S4.

(F) For visual cues, no significant differences were found in the percentages of bats’ choices in the three experimental groups of different light intensities (darkness, moonlight, and dim light). p values are from two-tailed binomial tests.

We conducted dual-choice acoustic playback (Figures 2B and S2A) and electrophysiological experiments to examine whether I. io may use passive auditory cues (eavesdropping on bird calls) to prey on birds at night. First, we studied the behavioral responses of 28 individuals of I. io to the playback of sound files of P. inornatus and Z. japonicus (Figure S3) and of silence files in a dual-choice experimental flight cage (Figures 2B and S2A; see Video S2 as an example). We found no significant difference in the number of flights to the selection chambers in bats that responded to the playback files between silence and bird calls (two-tailed binomial test, all p > 0.05; Figure 2D). There was no significant difference in the number of bats choosing the chambers that presented bird specimen + silence versus bird specimen + bird calls (two-tailed binomial test, all p > 0.05; Figure 2D). In addition, model analysis showed that none of the variables (i.e., flight latency, flight time, the number of echolocation pulses, playback position, and playback side) had an effect on bats’ choice preference (GLMM, all p > 0.05; Table S2) in the silence versus bird calls experimental group. In the specimen + silence versus specimen + bird calls experimental groups, the flight latency of bats that chose the specimen + bird call chamber was significantly shorter than that of bats that chose the specimen + silence chamber when the calls of P. inornatus were broadcast (GLMM: estimate = −0.019, z = −2.514, p = 0.012; Table S3). The number of echolocation pulses emitted by bats that chose the specimen + bird call chamber was significantly greater than that of bats that chose the specimen + silence chamber when the calls of Z. japonicus were broadcast (GLMM: estimate = 0.154, z = 2.117, p = 0.030), whereas no differences were observed in other variables (GLMM: all p > 0.05; Table S3). These results showed that although bats flying to the chamber that presented specimen + bird calls would make some behavioral changes (such as changes in flight latency and the number of emitted echolocation pulses), these behaviors had no significant impact on the bats’ choice preference. The results suggested that the behavioral choice preference of I. io to bird calls was not stronger than that of I. io to silence. In addition, we measured the hearing of three I. io individuals based on the auditory brainstem response (ABR) at a frequency of 1–8 kHz. This frequency range covers the call frequency of five representative birds (including P. inornatus and Z. japonicus) in the diet of I. io, with a frequency range of 3.6–7.3 kHz and a peak frequency range of 3.9–6.1 kHz (Table S4). The ABR results showed that when the frequency of sound playback was 6–8 kHz and the sound intensity was 70–90 dB peSPL, the three bats showed typical ABR waveforms, of which the first (P1), second (P2), and fourth peaks (P4) were relatively clear; however, when the sound intensity was lower than 70 dB peSPL, the typical ABR waveform curve did not appear, and when the sound frequency was lower than 6 kHz, the bats did not show the typical ABR waveform curves (Figures 2E and S4 and S5). The call intensity of Z. japonicus was about 75 dB peSPL at a distance of 1 m (temperature 20°C, humidity 90%, atmospheric pressure 852 hPa). When the propagation distance was greater than 1.8 m, the intensity of the bird calls was less than 70 dB peSPL. Taken together, the behavioral and electrophysiological results suggest that I. io may not use passive auditory cues to prey on nocturnal birds, because bats were insensitive to bird calls below 6 kHz, and it was difficult for the bats to perceive bird calls when they were far away from the birds.

Video S2. The behavioral choice preferences of I. io in the experimental group of silence versus bird call and specimen + silence versus specimen + bird call in the auditory cues experiment, related to Figure 2

The bird calls cannot be heard on the video because of the high-speed camera recording.

Furthermore, we conducted a dual-choice experiment to examine whether I. io may prey on nocturnal birds by using visual cues (Figures 2C and S2B). I. io has relatively large eyes compared with many other aerial-hawking insectivorous bats [e.g., Eptesicus nilssonii (Eklöf et al., 2002)], implying that they may have relatively good vision for finding prey (Land and Nilsson, 2002). We simulated the behavioral response preferences of I. io to bird specimens under different light intensities in the field through experiments at three light intensities: darkness (0–0.02 lux), moonlight (0.05–0.40 lux), and dim light (1.50–4.00 lux; see STAR Methods and Video S3 as an example). We found no significant difference in the choice preference of bats on bird specimens between dim light and darkness (two-tailed binomial test, all p > 0.05; Figure 2F). Similarly, in the choice experiments of darkness versus moonlight and moonlight versus dim light, no significant differences were observed in the bats’ choices (two-tailed binomial test, all p > 0.05; Figure 2F). In addition, model analysis showed that none of the examined variables (i.e., flight latency, flight time, the number of echolocation pulses, playback position, and playback side) had an effect on bats’ choice preference (GLMM, all p > 0.05; Table S5) in the three experimental groups of different light intensities. These results suggested that the increased lighting conditions did not increase the probability of I. io finding prey. Thus, I. io may not use vision to locate and prey on nocturnal birds in the field.

Taken together, our results showed that I. io uses active echolocation cues that are able to detect and locate birds in flight in the night sky, regardless of the magnitude of the light intensity and whether the migratory birds give calls. However, this study was an indoor behavioral simulation experiment that may have several limitations. For example, the bats could not show realistic flight foraging ability because of the limitation of the experimental space; the bird specimens could not replace real live birds and therefore could not represent the natural foraging scenarios. Thus, future research should use more advanced techniques to investigate the foraging sensory cues utilized by these bats in the wild. Nevertheless, the results provide the first experimental evidence to support the hypothesis that avivorous bats use acoustic methods similar to those they use to locate large moths to detect small flying passerine birds (Ibáñez et al., 2003), and also opposed the accidental consumption hypothesis that bats inadvertently ingest feathers shed by nocturnally migrating birds (Bontadina and Arlettaz, 2003).

Genomic adaptation

I. io preying on nocturnally migrating birds face physiological challenges, such as low oxygen levels, in addition to requiring specializations in feeding, digestion, and metabolism. Flight is a form of locomotion that requires large amounts of oxygen and energy, and also produces reactive oxygen species (ROS) that may lead to DNA damage. Bats in sustained flight consume more than 20 times the oxygen they do during resting (Maina, 2000). Here, our GPS results showed that the flight altitude of I. io during preying on nocturnally migrating birds reached 2068.03 ± 563.90 m asl, and the maximum flight altitude was 4945 m asl. The atmospheric pressure decreases with altitude, and the partial pressure of oxygen (PO2) dropped to about 50% at 5000 m asl (Hill et al., 2016). Thus, I. io faces significant physiological stress when flying at high altitudes to prey on birds owing to the additive effects of energy consumption during flight and low PO2 at high altitudes. Furthermore, I. io when biting and chewing birds may require greater bite force than that used for insects because of the denser and tighter muscles of birds. Additionally, iron ions and fats are higher in the blood and tissues of birds than in insects (Orkusz, 2021). An excessive intake of iron may result in iron poisoning (Zepeda Mendoza et al., 2018) that can affect the normal functioning of organs. Considering these factors, we hypothesized that I. io may have genomic adaptations to meet the above challenges. To test this hypothesis, we obtained a high-quality genome of I. io and compared it with the genomes of 16 other bats with different diets (Table S6).

We performed high-throughput whole-genome sequencing of a wild-caught I. io by using PacBio continuous long reads (20–40 kb insert size) and Illumina pair-end 150 bp reads (350 bp insert size), and further extension of the genome to the scaffold-level using Hi-C data (STAR Methods). Finally, we obtained the genome of I. io at ∼2.1 Gb (2,099,682,539 bp), and approximately 98.20% of the 25 sequences were mounted to the chromosome-level (Figure S6 and Tables S7 and S8). The scaffold N50 of the genome was 105.83 Mb (105,833,403 bp; Figure S7). BUSCO (Waterhouse et al., 2018) indicated a high integrity of the assembled genome (92.5%; Figure S8). In addition, 99.30% of the Illumina pair-end reads library sequences covered 98.71% of the genome (Table S9), indicating a high sequence identity of the assembly. We successfully annotated 19,551 protein-coding sequences after integrating three annotation strategies (Table S10): Homology prediction of six bat species closely related phylogenetically, de novo prediction, and transcriptome sequences from 10 tissue samples. The annotated protein-coding sequences achieved 90% completeness of 4104 BUSCO mammalian genes set (Figure S8). These results showed that our assemblies and annotations were complete for protein-coding sequences. The functions of 18,732 (95.81%) sequences were successfully annotated by searching against several databases (Table S11).

We screened the genomes of 20 mammalian species to identify single orthologous genes (Table S6). The species included insectivorous bats, frugivorous or nectarivorous bats, sanguivorous bats, and avivorous bats as well as humans (Homo sapiens), cats (Felis catus), and dogs (Canis lupus familiaris) (Figure 3A). A total of 4027 orthologous genes were identified for downstream analysis. The protein sequences, cDNA sequences, and codons (1st+2nd, 3rd) were used to construct a maximum likelihood phylogenetic tree. In the present analysis using various types of orthologous genes, we obtained the same phylogenetic topology with 100% support for all nodes (Figures S9–S11). Using this phylogenetic tree with I. io as the foreground branch, we obtained 187 rapidly evolving genes (REGs) (Table S12), and 100 positively selected genes (PSGs) (Tables S13, S14, and S15; see STAR Methods). Analyzing gene family expansion and contraction in I. io revealed that there were 207 expanded and 31 contracted gene families (Figure 3B and Table S16). Given the similarity in the specific diets of I. io and Desmodus rotundus (containing blood), we examined the convergence between the two species. Our results showed that eight genes had significant convergent evolution in I. io and D. rotundus (Table S17). Functional enrichment analysis of the PSGs, REGs, and expanded gene families in the I. io branch showed enrichment in DNA damage and repair, hypoxic response, biting and mastication, and digestion and metabolism (Tables S18, S19, S20, and S21).

Figure 3

Phylogenetic and gene family analyses

(A) The phylogenetic tree of I. io and 19 other species. The phylogenetic tree was constructed by RAxML using 4027 orthologous genes among 20 species. Phylogenetic trees with the same topology and 100% node support were constructed using protein, cDNA, and codon data from orthologous genes. The blue insects represent insectivorous bats; the yellow bananas represent frugivorous or nectarivorous bats; the red blood represents a sanguivorous bat, and the orange bird represents an avivorous bat. In this phylogeny, the branches of I. io and D. rotundus are shown in bold. Similar divergence times were estimated by r8s and MCMCtree in this study (Figures S14–S16). This plot illustrates the results of MCMCtree. The unit of divergence time is millions of years ago (MYA). Images of bats were downloaded from the iNaturelist website. All images are CC-BY-NC 4.0 licensed.

(B) The bidirectional histogram represents the numbers of expanded and contracted gene families. Red and blue represent the expanded and contracted gene families, respectively. The X-axis represents the number of gene families. The Y-axis of the histogram represents different species arranged in the same order as the 20 species in the phylogenetic tree.

Bats consume a substantial amount of oxygen and energy during flight, thereby producing ROS (Zhang et al., 2013) that have detrimental effects on DNA. We found two important genes, DDIT3 and INO80B, to be positively selected in I. io (Figure 4A). Moreover, we found that 10 REGs (EMSY, RUVBL2, TIMELESS, FEN1, DDB1, FANCF, PHF1, PDRX1, RAD23B, and ERCC6) (Table S12) and two expanded genes (UBE2V2 and UBE2N) were involved in the detection and repair of DNA damage. Further analysis of the functions of these genes revealed that they also play a role in hypoxic responses in addition to DNA damage repair. For example, DDB1 (Katiyar et al., 2009), PHF1 (Millonig et al., 2009), and ERCC6 (Filippi et al., 2008) are involved in the adaptation to hypoxia mediated by HIF1. Then, we focused on genes related to hypoxic response and observed positive selection on both the WDR83 and SFPTB genes (Figure 4A). The amino acid changes at the WDR83 positive selection site (Ser70Gln) may negatively affect its function (PROVEAN score = −2.948 < −2.5) and lead to changes in the three-dimensional structure of the WDR83 protein (RSMD = 0.067; Figure 4B). WDR83 is known to promote the degradation of the protein encoded by HIF1A at normal levels (Su et al., 2012). At hypoxic levels, decreased WDR83 expression reduces the inhibition of HIF1A activity, reducing the effect of low oxygen levels in tissues (Boulahbel et al., 2009; Hopfer et al., 2006). In addition, SFTPB encodes the pulmonary-associated surfactant protein B (SPB) (Whitsett and Weaver, 2002). The surface tension and alveolar volume of the alveoli are maintained by SPB (Ikegami et al., 2009). Of interest, compared to non-flying mammals of similar size, the lung volume of bats is about 72% larger (Canals et al., 2005). In addition to the above findings, we also obtained evidence of convergent evolution to support hypoxic adaptation in I. io. In I. io and Pantholops hodgsonii, we identified SIRT6 and UHRF1 (Table S17). SIRT6 encodes a protein belonging to the sirtuin family, and its over-expression can activate the pAMPK alpha pathway, increasing the protein level of BCL2 and decreasing the protein levels of pAkt and ROS to protect cardiomyocytes of mice during hypoxia (Maksin-Matveev et al., 2015). UHRF1 encodes a ubiquitin-like protein with PDH and ring finger domains and is expressed at high levels in high-altitude Tibetan pigs (Ban et al., 2015). In addition, we found that ALDH5A1 has undergone convergent evolution in I. io, Bos mutus, and Rhinopithecus bieti (Table S17). ALDH5A1 encodes a protein that belongs to the aldehyde dehydrogenase family. ALDH5A1 plays a role in myocardial hypoxia damage in TOF patients (Liu et al., 2021). Thus, these results suggest that adaptive evolution of the above genes enables I. io to better tolerate DNA damage and hypoxic response during flight and hunting at high altitudes.

Figure 4

Positively selected genes (PSGs) and expanded gene families in connection with preying on nocturnally migrating birds by Ia io

(A) The MSA plot of genes that have undergone positive selection and convergent evolution. Here, “.” represents the same residue as in the human protein. The red background represents residues that are positively selected sites. The blue background indicates amino acid sites that have undergone convergent evolution between I. io and D. rotundus. The ancestral amino acids in parentheses were reconstructed using CodeML.

(B) WDR83 protein structure. The inferred 3D structure for WDR83 in I. io was determined using SWISS-MODEL with template-based modeling. To explore the effect of the substitution on the protein stability, the positively selected sites of I. io (red) were replaced with the amino acids of humans (blue). The root mean square deviation of atomic positions (RSMD) was calculated using PyMol.

(C) The functions of PSGs and expanded genes in the mitochondria of cardiomyocytes. The genes with a blue background are expanded genes, whereas those with a red background are PSGs.

(D) The functions of expanded genes in iron metabolism. The genes with a blue background are expanded genes.

Bats have an extremely variable heart rate to ensure adequate oxygen delivery to tissues via blood flow (Canals et al., 2011). In addition, bats also have the highest level of ATP in the heart muscle (Neuweiler and others, 2000). Here, expanded gene families were enriched in items pertaining to myocardial contraction and mitochondria. Genes enriched in the mitochondria are primarily related to the mitochondrial inner membrane and electron transport chains. Genes enriched in myocardial contraction, including CACNA2D3, VDAC2, and RYR2, are primarily involved in the influx of Ca2+ into the mitochondria of cardiomyocytes (Figure 4C). Interestingly, NCLX, a gene related to calcium efflux, has undergone positive selection (Leu134Ser) (Figure 4A). Ca2+ flux and homeostasis play important roles in myocardial energy production and regulation of cardiomyocyte contractility (Cao et al., 2019). CACNA2D3 encodes a subunit of the protein in the L-type calcium channel (Bodi et al., 2005). After sarcolemma depolarization, Ca2+ enters the cell via the opened L-type calcium channels (Cao et al., 2019). Then, ryanodine receptor 2, encoded by RYR2, is stimulated by Ca2+, resulting in a greater release of Ca2+ from the sarcoplasmic reticulum into mitochondrial-associated membranes (Min et al., 2012). Ca2+ enters the mitochondria via voltage-dependent anion channel 2 encoded by VDAC2 (Min et al., 2012). Consequently, Ca2+ binds to troponin C (through ATP consumption), resulting in cardiomyocyte contraction (Modesti et al., 2021). However, the dynamic equilibrium between Ca2+ influx and efflux is essential for maintaining Ca2+ homeostasis in cells. Here, the IMM-localized protein encoded by NCLX is responsible for excreting Ca2+ in exchange for either Li+ or Na+ (Luongo et al., 2017).

Because birds have bones and muscles, I. io may require more developed masticatory muscles and a greater bite force to consume birds, as opposed to insects. Our previous study demonstrated that I. io has a higher bite force compared to the same-sized insectivorous bats (Shi et al., 2020). However, bite force is determined by the structure of the skull and the muscles, especially the masticatory muscles (Santana, 2018). A selective pressure analysis of the I. io branch revealed that EGR2 underwent positive selection (Thr259Gly, Pro260Gln) (Figure 4A), and EGR1 underwent rapid evolution. EGR2 has been shown to affect the development of masticatory muscles by regulating the number of trigeminal nerves (Turman, 2007), and EGR2 deficiency leads to abnormal bite force and feeding problems in mice (Turman et al., 2001). EGR1 and EGR2 work together to control type I collagen and thereby maintain normal muscle function (Lejard et al., 2011). Thus, the adaptive evolution of the two genes in I. io may facilitate the development of stronger masticatory muscles to produce the high bite force that may enable I. io to feed on birds, which are larger than insects.

Expansion of the dietary niche of I. io from insects to birds may have required the adaptive evolution of metabolism-related genes. It is unknown whether I. io evolved adaptations in genes related to iron metabolism to cope with iron toxicity caused by excess iron ions in predated birds. Our results showed that the genes FTH1, PCBP1, and HEPHL1 involved in iron metabolism were expanded (Figure 4D). HEPHL1encodes a copper-binding glycoprotein with ferroxidase activity that can oxidize Fe2+ to Fe3+, and this may be involved in the regulation of intracellular iron content (Sharma et al., 2019). The heavy chain of ferritin is encoded by FTH1 (Hentze et al., 1986), and ferritin stores iron ions transported by the iron chaperone protein encoded by PCBP1 (Pozzi et al., 2015; Shi et al., 2008), thereby regulating the dynamic balance of iron ions in the tissues (Shi et al., 2008). A previous study has also shown that the FTH1 gene of D. rotundus was also expanded (Zepeda Mendoza et al., 2018).

Furthermore, our results indicated that TM4SF5 and TMEM74 in I. io and D. rotundus exhibited convergent evolution (Figure 4A). TMEM74 encodes Transmembrane Protein 74, and TM4SF5 encodes the Transmembrane 4 L Six Family Member 5. Both genes are involved in cell proliferation in the liver and the development of liver cancer (Ryu et al., 2021; Sun et al., 2019). TM4SF5 also plays a role in fat accumulation and the metabolism of fatty acids in the liver (Park et al., 2021). This suggests the adaptive evolution of I. io genes to cope with the challenges of lipid metabolism for digestion of bird prey. A sanguivorous diet facilitates exposure to blood-borne viruses, and an avivorous diet facilitates exposure to avian-borne viruses. In the convergence analysis of I. io and D. rotundus, we also found four convergent genes involved in the immune response, GBF1, LRCH4, CD5, and NCF1 (Table S17). GBF1 plays an important role in the replication of different RNA viruses (Martinez and Arias, 2020). LRCH4 can regulate the innate immune response (Aloor et al., 2019), and CD5 is a key component of adaptive immunity (Pers et al., 1999). Moreover, NCF1 is involved in innate immunity, adaptive immunity, and autoimmune diseases (Efimova et al., 2011; Imai et al., 2008; Kraaij et al., 2010). These results indicated that adaptation of immune-related genes may have helped I. io to cope with the immune challenges of an avivorous diet. Thus, our results suggest that the adaptive evolution of genes associated with metabolism and immune response may facilitate I. io to better digest and metabolize bird prey.

Taken together, the results of our study indicate that I. io has experienced adaptive evolution in genes associated with DNA damage repair, hypoxia adaptation, bite force, and digestion and metabolism. These adaptions demonstrate the evolutionary novelty in I. io genome that may allow them to adapt to the low-oxygen environment of flying at high altitudes to prey on birds, as well as the physiological challenges of feeding, digesting, and metabolizing bird prey.

In conclusion, niche expansion may be limited by the genetics and phenotypes of organisms (McCormack and Smith, 2008; Sexton et al., 2017). Evolutionary novelty and innovation (e.g., behavioral innovation) are believed to represent the ability of species to exploit a wider range of resources, successfully colonize new habitat types or niches, and adapt to novel environments (Erwin, 2015; Erwin et al., 2011; Wang and Liu, 2021). However, evolutionary novelty and innovation are rare (Erwin, 2021). Ecological opportunity is a primary driving force in the origin and success of evolutionary novelties and innovations (Erwin, 2021). In this study, we demonstrated how I. io adapts to preying on nocturnally migrating birds (ecological opportunity), one of the most challenging foraging tasks in the animal kingdom regarding behavior, sensory modalities, and physiology. Our results suggested that I. io exhibits flight behavioral innovations that may allow them to overcome the challenges of flight altitude and speed for preying on nocturnally migrating birds. I. io uses the active echolocation system that can to detect and locate birds, breaking through the visual sensory limitations of low visibility at night. In addition, I. io exhibits genomic novelties that may enable them to adapt to the low-oxygen environment of flying at high altitudes to prey on birds and the physiological challenges of feeding, digesting, and metabolizing birds. These results suggest that I. io is able to prey on nocturnally migrating birds based on existing flight and echolocation behavior, as well as the evolution of flight behavioral innovations and genomic novelties. This study suggests that when animals are faced with challenging ecological opportunities, they can evolve behavioral innovations and genomic novelties to enhance or innovate their foraging abilities, thereby exploiting novel dietary niches.

Limitations of the study

The high flight altitude of I. io foraging in autumn may not necessarily indicate hunting for birds, and thus, future studies should be performed to clarify the fact using on-board sound recordings or advanced animal-borne sensors. Moreover, the behavioral experiments in this study were indoor experiments, and thus future studies should be performed in natural scenarios in the wild, although this is difficult and challenging. Finally, our genomic study only focused on one of the avivorous bats, I. io, and thus, further studies should include the other two species of avivorous bats (N. lasiopterus and N. aviator) to verify the universality of genomic adaptation to bird prey.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Tinglei Jiang (jiangtl730@nenu.edu.cn).

Materials availability

This study did not generate new unique reagents.

Data and code availability

The genomic data were deposited in the BioProject database of the NCBI under accession no. PRJNA782023. All other data are available in the manuscript or the supplemental information.

This paper does not report original code.

Any additional information required to reanalyze the data reported in this paper is available from the lead contact on request.

Experimental model and subject details

Animals

A total of 69 adult I. io males were captured using mist nets at the entrance to the Feilong Cave (Xingyi City, Guizhou Province, China) for this study. Bats were brought to a field temporary laboratory and held in a small flight cage (1.5 × 1.5 × 1.5 m). Each bat was weighed and measured for the forearm length and then marked with a numbered split aluminum alloy bat ring (5.2 mm, Porzana Ltd., Icklesham, UK) for individual identification. The temperature in the field temporary laboratory ranged from 15 to 25°C; the relative humidity was between 70 and 90% and was recorded by a meteorograph (Kestrel 5000, Nieisen-Kellerman, Boothwyn, Pennsylvania). Following the local climate conditions, a naturally fluctuating day/night light cycle (natural day-night cycle) was presented to the bats. All bats were given ad libitum access to fresh water and Zophobas morio larvae enriched with vitamins and minerals.

One wild adult bird species, the Mountain White-eye (Zosterops japonicus, n = 5), was purchased from several breeders at a local Flower & Bird Market. These birds were also brought back to the field temporary laboratory and housed in a birdcage (0.5 × 0.3 × 0.4 m). The birds were given ad libitum mixed seeds, Tenebrio molitor larvae, fruit (apples and peaches), and water, and baths were regularly provided. Additionally, the adult specimens of two representative bird species in the diet of I. io, the Yellow-browed Warbler (Phylloscopus inornatus, n = 2) and Z. japonicus (n = 2), were purchased from Guangzhou Lingchuang Specimen Stationery Co., Ltd., China.

After the experiment, all bats were released back into the cave, and the birds were released back into the wild. Animal use and care were performed under the Guidelines of the Use of Animals in Research (ASAB/ABS) and the Northeast Normal University Guidelines for Animal Research and conformed to all applicable governmental regulations concerning the ethical use of animals and biodiversity rights. This research was approved by the Laboratory Animal Welfare and Ethics Committee of Jilin Agricultural University and the Wildlife Conservation Office of the Jilin Forestry Department, China.

Method details

Field GPS experiments

GPS attachments

The experiment was performed in November (autumn) 2020 and June-July (summer) 2021. Sixteen adult male I. io individuals were captured in summer and ten in autumn. The body mass of I. io was 50.23 ± 3.30 g in summer and 58.70 ± 4.97 g in autumn. We gently removed part of the fur between the bats’ shoulder blades and used skin glue (Surgical Cement, China) to attach Global Positioning System loggers (GPS, HQBG0603, Hunan Global Messenger Technology Co., Ltd., Hunan, China) with the Global System for Mobile Communication (GSM, MC20, Quectel Wireless Solutions Co., Ltd., Shanghai, China) to the backs of the bats. The GPS data could be transmitted back via the GSM module. We kept bats in the hand for 60 s to prevent the logger from shifting. After each bat was fully tagged, they were placed separately in a cloth bag and adapted for 30 min to observe whether the GPS logger was completely attached to the bat’s back. Subsequently, we released the bats to the roost. To reduce additional interference to the bats, we completed the attachment of GPS loggers and release of all bats within one hour. Bats had about 12 h in the roost to adapt to the additional weight of the GPS loggers before they conducted night foraging activities.

All 26 bats were tagged with GPS loggers, and the average weight of the coating and the total GPS logger was 4.7 ± 0.1 g, approximately 7.36% to 9.94% of the body mass of the bats. The best practical rule suggests that bats’ load weighing should not exceed 5% of their body mass (Aldridge and Brigham, 1988), while some previous studies have suggested that a ratio less than 10% may not significantly affect the normal activity of bats (Cvikel et al., 2015; O'Mara et al., 2014). Moreover, our preliminary field experiments also showed that the GPS loggers did not significantly affect the activities of I. io. The sampling time in each season 1 h before sunset to 1 h before sunrise, i.e., 19:00–7:00 in autumn and 20:00–7:00 in summer (UTC+8h). The logging intervals were set to every 10 min in autumn and one min or 10 min in summer (Table S22).

Data collection and analyses

GPS data from 20 bat individuals were obtained, although the data of six bats are missing. This may be because vespertilionid bats (i.e., I. io) prefer to roost together, which may lead to tags falling into cracks in rocks during crawling or fighting. For GPS data, we first manually deleted any incorrect sites that were due to the difficulty in obtaining satellite signals. Finally, we recorded GPS data for 28 nights in summer and 22 nights in autumn, with an average of 2.5 ± 1.2 nights per bat. We obtained 1700 location data records in summer and 528 in autumn.

We tested the errors in GPS locations in two-dimensional latitude and longitude and vertical altitude by the Dilution of Precision (DOP), a unit used to measure the geometric effects of satellites that can be employed to evaluate the quality of locations. The test considered the position of each satellite relative to other satellites to predict the position accuracy that can be obtained. A lower DOP value indicated that the obtained sites have better repair quality and often have small positioning errors and high positioning accuracy. To use high-quality sites as much as possible, we used the 2D plane coordinate precision factor (HDOP) and the vertical coordinate precision factor (VDOP) to define the quality of location data. The maximum of the two values can be used to determine the grade of the site. Finally, we rated the grade of the data as A–E (see Table S23), and the data with D and E grades were discarded. Here, about 1.2% of the sites were discarded. Finally, we obtained 1678 location data points in summer and 521 in autumn for final analysis.

To show the flight altitude of bats at night more intuitively, we used sea level and ground level as the respective base levels to obtain the bats’ flight altitude above sea level and above ground level. When calibrating the flight altitude of I. io above ground level in the two seasons by using the high-precision 12.5 m Digital Elevation Model (DEM) elevation data, which were derived from NASA (National Aeronautics and Space Administration; https://search.asf.alaska.edu), there were still 12.1% of the sites that produced a false negative altitude of 4.3 m below the ground (median) due to the error of the altitude recorded by GPS and the deviation between the whole grid elevation and the actual height. In autumn, only six sites (1.1%) of the bats’ flight altitude above ground level were affected, while in summer 199 sites (11.7%) were affected. This had no effect on the flight altitude of bats above ground level in autumn, but the actual flight altitude above ground level in summer should be lower than the results of our analysis suggest.

We analyzed and plotted the flight altitude above sea level and above ground level of I. io in two seasons based on sea level and ground level as the respective base levels. The flight speeds of I. io individuals were directly obtained from the instantaneous speed of bats flying at each site through the speed sensor carried by the GPS logger. The Mann–Whitney U test was used to analyze the differences in flight altitude and flight speed between summer and autumn. To compare with other known recorded flight altitude data for bats, we obtained the maximum flight altitude above ground level of the bats that were recorded by using GPS technology and plotted the data as a histogram. The bats were one fruit bat, Rousettus aegyptiacus (Tsoar et al., 2011), five insectivorous bats, Tadarida teniotis (O'Mara et al., 2021), Taphozous theobaldi (Roeleke et al., 2018), Rhinophylla microphyllum (Cvikel et al., 2015), Nyctalus noctula (Roeleke et al., 2016), and Rhinolophus ferrumequinum (Fujioka et al., 2020), and one avivorous bat, Nyctalus lasiopterus (Naďo et al., 2019). All statistical analyses were performed in R 3.5.0 (R Core Team, 2018); maps were drawn using ArcGIS 10.2 (Esri), and the three-dimensional map is presented with 2× vertical exaggeration.

Behavioral experiments with sensory cues

Echolocation cues (active acoustic detection) experiment

To verify whether the bats used an aerial-hawking strategy combined with echolocation cues to prey on birds, we carried out mimicking behavioral experiments under darkness in November 2019. A total of eleven adult male I. io individuals were brought to the field temporary laboratory and used for this study. In the training phase, we used cotton threads to hang bird specimens (two P. inornatus and two Z. japonicus) at different heights in a flight room (9.0 × 3.7 × 3.9 m). All bats were released into the flight room in turn, and once they learned to successfully attack tethered bird specimens, the formal experiment started the next night. A successful attacking behavior was defined as bats using the wings to flap or directly attack bird specimens and/or emitting echolocation calls with a feeding buzz. Feeding buzzes of bats were defined as echolocation calls emitted just before prey capture, where the bats increase the pulse repetition rate and progressively shorten individual calls; these calls can indicate whether a bat has successfully captured a prey item (Simmons et al., 1979; Surlykke et al., 2003). In the experiment, a P. inornatus specimen (Bird 1) was suspended from cotton threads at 2 m in front of the horizontal direction of the center microphone in the microphone array (see below) and 2.1 m above the ground in the flight room. Meanwhile, we hung two Z. japonicus specimens (Bird 2 and Bird 3) at 2 m in front of the horizontal direction of Bird 1 and 1.5 m apart from each other at the same height (Figure S1). We randomly selected one individual bat at a time per night (20:00–24:00) and released the bat into the flight room. Each bat performed successive multiple trials until it stopped attacking and started hanging at the wall and grooming itself.

A planar star-shaped four-microphone array (condenser ultrasound microphone CM16/CMPA, Avisoft Bioacoustics, Berlin, Germany) was placed at the middle of one end of the flight room at a distance of 2.1 m from the ground to record echolocation calls during the attack (Figure S1). Three peripheral microphones in the microphone array were set 60 cm apart from the center microphone, and the included angle was 120 degrees (Figure S1). The microphone array was inclined at an angle of 20 degrees to the ground. The recorded echolocation calls of attacking bats were subsequently used to calculate their three-dimensional positions. Echolocation calls were recorded using an UltraSoundGate 416H (Avisoft Bioacoustics, Berlin, Germany) connected to a desktop computer running the Avisoft Recorder software with a 500 kHz sampling rate and 16-bit resolution. The microphones were calibrated by a GRAS 42 AB sound calibrator. We recorded the bats’ attacking behavior from two high-speed, infrared-sensitive cameras (MV1-D1312IE-240-CL-8, Photonfocus AG, Lachen, Switzerland) with a 120 frames/s sampling rate. A high-speed camera was placed 1.5 m above the ground and parallel to the microphone array plane, and another camera was placed at the same height and position on the other end of the wall (Figure S1). Simultaneously, we placed four customized infrared lights (KTJ-GY-300W-42V, RockTech Technology Corp., Ltd., Hunan, China) at 1.5 m above the ground in the four corners of the flight room for sufficient illumination of the infrared high-speed cameras. Audio and video recordings were synchronized using a synchronous controller (AcuteSync-20K, RockTech Technology Corp., Ltd., Hunan, China) connected to the same desktop computer. Additionally, we also used an infrared camera (HDR-CX 760E, Sony Corp., Tokyo, Japan, 50 frames/s sampling rate) to simultaneously film bats’ attacking behavior. A weather instrument (Kestrel 5000, Nieisen-Kellerman, Boothwyn, Pennsylvania) was used to record the air temperature and relative humidity in the flight room at an interval of 1 min.

Auditory cues (passive acoustic detection) experiment

To investigate whether the bats used auditory cues (eavesdropping) to prey on birds, we captured a total of 28 adult male I. io individuals used for experiments conducted from May to June (n = 14) and October to November (n = 14) 2019. The calls (flight calls) of two representative bird species in the diet of I. io, P. inornatus (highest consumption frequency), and Z. japonicus (medium consumption frequency) (Gong et al., 2021) were selected for playback. Flight calls are usually the primary vocalizations given by many species of birds during sustained flights, especially during migration; calls are generally in the 1–10 kHz range, and they have the function of maintaining flock contact during migration and stimulating migratory restlessness (Farnsworth, 2005; Hamilton, 1962). The flight calls of P. inornatus (six individuals) and Z. japonicus (five individuals) were collected from public sources of flight call recordings (Xeno-Canto: http://www.xeno-canto.org/). We used Format Factory (4.5.5, Shanghai Geshi Network Technology Co., Ltd., China; http://www.pcgeshi.com/) to convert the MP3 format of flight calls into WAV format. In the absence of acquired wild P. inornatus, the downloaded flight calls were directly used for subsequent editing and broadcasting. For Z. japonicus, according to the induced vocalization method in darkness from Lanzone et al. (2009) and Landsborough et al. (2018), we exposed wild birds to acoustic stimuli of downloaded flight calls to induce calling by using an ultrasonic loudspeaker (Ultrasonic Dynamic Speaker Vifa, Avisoft Bioacoustics, Berlin, Germany). We then used an UltrasoundGate 116H (Avisoft Bioacoustics, Berlin, Germany) to record the calls of five wild adult Z. japonicus, with a sampling frequency of 50 kHz and a resolution of 16 bits. A sound level meter (Casella CEL-633C1/K1, Kempston, Bedford, UK) was used to measure the intensity of the calls; this was about 75 dB SPL re. 20 μPa at 1 m. The call intervals of the flight calls from two bird species (460–1420 ms in P. inornatus, 120–1340 ms in Z. japonicus) were measured using Avisoft-SASLab Pro (version 5.2.07, Avisoft Bioacoustics, Berlin, Germany). In addition, we randomly selected one or two high-quality calls from each individual and obtained a total of 10 calls from six individuals of P. inornatus (Figure S3A) and three high-quality calls from each individual, and a total of 15 calls from five individuals of Z. japonicus (Figure S3B). The calls were then edited to generate two 10 s call sequences within their call intervals using Avisoft-SASLab Pro 5.2.07. The 10 s call sequences from the two bird species were repeated to generate two 3-min playback files. In addition, a 3 min silence stimulus was generated as a control.

The bats were tested in a dual-choice experimental flight cage (3.0 × 1.8 × 2.1 m). The dual-choice experimental flight cage was divided by a sliding baffle into four selection chambers (Chamber 1, 2, 3, and 4) at both ends, each 0.9 × 0.9 × 2.1 m (Figures 2B and S2A). We used two loudspeakers (Ultrasonic Dynamic Speaker Vifa, Avisoft Bioacoustics, Berlin, Germany) at each end of the flight cage to broadcast bird calls and silence playback files. Each loudspeaker was placed 1.65 m above the ground in the center of the outer wall of each selection chamber. This experiment consisted of two groups: 1) silence vs. bird calls; 2) specimen + silence vs. specimen + bird calls (Figure S2A). Bird specimens (Z. japonicus in May–June and P. inornatus in October–November) were hung in the center of each selection chamber 45 cm below the top of the cage (Figure S2A). Similarly, the bird specimens were alternated at both ends of the flight cage during the experiment. The bats were placed at a take-off spot (15 × 15 cm) connected to the top of the flight cage 20 cm from the wall of the cage, and the take-off spot was alternated between the ends of the cage. Before trials began, the individual bats were allowed to acclimatize to the dual-choice experimental flight cage for one night. The calls of two bird species were broadcast at 75 dB SPL re. 20 μPa at 1 m and measured with a sound level meter (Casella CEL-633C1/K1, Kempston, Bedford, UK) from the front of the loudspeaker. In each experiment, the playback files were broadcast after the tested bat was hung on the take-off spot and stabilized for 30 s. The call playbacks were continued for three min and then stopped. Trials for each bat were randomly repeated five times at both ends of the cage (Chambers 1 and 2, Chambers 3 and 4) with an interval of 2 to 3 d between two repeated experiments, resulting in 10 total trials per bat for the two treatment groups. Bats’ choices and behavioral and acoustic responses were recorded simultaneously using two ultrasonic microphones (CM16/CMPA, Avisoft Bioacoustics, Berlin, Germany) connected to a four-channel Avisoft UltraSoundGate 416H (sampling at 250 kHz) and two infrared high-speed cameras (Photonfocus MV1-D1312IE-240-CL-8; 120 frames/s sampling rate) (Figures 2B and S2A).

Measurement of low-frequency hearing

To test whether the bats could eavesdrop on the calls of birds, we measured the auditory brainstem response (ABR) of three adult male I. io individuals. Each bat was anesthetized by intraperitoneal injection of 0.75% pentobarbital sodium solution (0.07 mL/10 g) before surgery. Additional small doses of pentobarbital sodium solution (0.03–0.05 mL) were superadded during the experiments as appropriate. The standard for bats after anesthesia was steady breathing rhythm and no restlessness. During the experiments, the bats were placed in a shielded room and wrapped in a cotton gauze heating pad to maintain their body temperature. The bats’ scalp fur was partially shaved to expose the skin and muscle over the cranium. The inferior colliculus at the top of the skull was used as the recording electrode; the ear mastoid was used as the reference electrode, and the ground electrode was connected to the tip of the nose. Each location was pierced with a silver-plated special electrode for about 5 mm. The inter-electrode resistance between the needle electrodes was approximately 1 kΩ and not more than 3 kΩ, and the resistance was displayed by the resistance detector of the RA4LI module of the TDT system (Tucker-Davis Technologies, Alachua, FL, USA). Bats were fixed on a recording table in the shielded room during the experiments. Each bat was continuously tested for 3–4 h. After the experiments, bats were given drinking water and glucose and in time returned to normal.

Sound stimuli were generated by the TDT BioSigRZ software program (sample rate 200 kHz). Sound stimuli consisted of clicks and 2 ms tone pips (cosine window, 1 ms rise/fall times) ranging in frequency from 1 to 8 kHz. The stimuli were delivered through an A/D converter (analog-digital converter, 16-PCM, Tucker-Davis Technology, Alachua, FL, USA) and an amplifier (RA16PA, Tucker-Davis Technology, Alachua, FL, USA), and then presented to the bat using an electrostatic speaker (ES1, Tucker-Davis Technology, Alachua, FL, USA) positioned 3 cm in front of the tested ear of each bat. Sound stimuli were automatically attenuated by a programmable attenuator (PA5, Tucker-Davis Technology, Alachua, FL, USA) from 90 dB peSPL to 0 dB peSPL in 5 dB peSPL gradients/steps. ABRs were recorded using stimuli presented at a rate of 21.01 per second, and 512 repetitions were collected over a time window of 10 ms and averaged. The final signal was filtered with a high-pass 300 Hz and a low-pass 3000 Hz by a TDT RP2.1 Real-Time Processor and stored in a computer.

For the off-line analysis, a digital band-pass filter (0.3–3 kHz) was applied to reduce background noise. The averaged waveform responses over the 512 repetitions of each frequency at each sound level were calculated using the TDT BioSigRZ Windows application. We presented the waveform that occurred within the first 7 ms after the sound reached the bat’s ears, and four to five positive peaks were visible within this time window. The typical ABR waveform curve follows the principle that the ABR wave appears the most stable and has the highest appearance rate as well as the ABR wave peaks being able to appear at lower sound stimulus intensity.

Visual cues experiment

To investigate whether bats used visual cues to prey on birds, we obtained data from 14 adult male I. io individuals (these were the individuals used in the auditory cues experiment from May to June) to perform an indoor simulation experiment in May and June 2019. We simulated the illuminance of bats when they encountered birds at night. In the corresponding periods of different illuminances, we used an illuminometer (accuracy = 0.01, Pro’sKit MT-4617LED-C, Pro’sKit Industries Co. Ltd., Taiwan, China) to measure the light intensity in the local area. The measurements were made on multiple days and during multiple periods to finally obtain range values for the simulation experiment. The actual measurements were darkness (0–0.02 lux), moonlight (0.05–0.40 lux), and dim light (1.50–4.00 lux; i.e. the illuminance during the period of the bats flying out after dusk and returning to the cave before dawn). We used 6-W LED table lamps (TL8012, Glonew Lighting, Shenzhen Glonew Lighting Co., Ltd., Guangdong, China) with intelligent remote control and placed them in the central position below each selection chamber of the dual-choice experimental flight cage at 40 cm above the ground to simulate different light intensities. An illuminometer (Pro’sKit MT-4617LED-C, Pro’sKit Industries Co. Ltd., Taiwan, China) was used for synchronous measurement to record the actual ranges of the three lighting conditions (i.e., darkness, moonlight, and dim light). This experiment comprised three groups: 1) specimen + darkness vs. specimen + moonlight; 2) specimen + darkness vs. specimen + dim light; 3) specimen + moonlight vs. specimen + dim light (Figure S2B). For each experiment, two Z. japonicus specimens were hung in the center of two selection chambers on the same side of the dual-choice experimental flight cage 45 cm below the top of the cage (Figure S2B). The bat individuals were allowed to acclimatize to the dual-choice experimental flight cage for one night and to the experimental set-up before trials began. The LED table lamps were turned on to provide the corresponding illumination for 5 min after the tested bat rested for 30 s at the take-off spot. Each bat was tested using randomly repeated trials five times at both ends of the cage (Chambers 1 and 2, Chambers 3 and 4) with an interval of 2 to 3 d between two experiments, resulting in 15 total trials per bat for three experimental groups. Similarly, we used two ultrasonic microphones (CM16/CMPA, Avisoft Bioacoustics, Berlin, Germany) connected to a four-channel Avisoft UltraSoundGate 416H (sampling at 250 kHz) and two infrared high-speed cameras (Photonfocus MV1-D1312IE-240-CL-8; 120 frames/s sampling rate) to synchronously record bats’ choice preferences and behavioral and acoustic responses (Figures 2C and S2B).

Behavioral and statistical analyses

Audio recordings were examined using Avisoft SASLab Pro (version 5.2.07, Avisoft Bioacoustics). All videos were viewed and analyzed using high-speed image acquisition and playback system (RockTech Technology Corp., Ltd., China) and QvodPlayer (Shenzhen Qvod Technology Co., Ltd., China). We used the microphone array recordings to reconstruct the flight paths. We estimated the bat’s position for each echolocation call by determining the time-of-arrival differences of the call between the four microphones by cross-correlation. We defined the choice decision as the bat individuals successfully flying to the playback selection chamber, attacking, or detecting targets (bird calls and specimens), and then landing or leaving the selection chamber. The numbers of bat individuals that selected playback files or control files or did not make a choice were counted via observing the recorded high-speed videos of bats making choice decisions after experiencing bird calls and being presented with bird specimens under different lighting conditions. By analyzing synchronously recorded high-speed video and audio files, we calculated the bats’ response time (flight latency, the time from when starting to playback the sound file and turning on the LED table lamp to leaving the take-off spot), flight time (time from when bats left the take-off spot to landing or leaving the selection chamber after making a choice) and counted the number of echolocation pulses of the bats that made a choice during the flight.

To exclude the familiarity effect (i.e., side biases) of bats preferring to make choices at either end of the flight cage, we first analyzed the number of choices made by all bats in the first round, and then analyzed the number of choices of all five rounds. A two-tailed binomial test was used to test for differences in the preference of bats’ choices in response to broadcasting different sound files and the presentation of bird specimens under different lighting conditions. We found that the results of analyzing only the first round were consistent with those of analyzing all five rounds (i.e., the two analyses showed no differences in the preference of bats in all experimental groups; Figures 2D, 2F, S12, and S13), In the present study, we therefore chose to analyze all five rounds. Additionally, generalized linear mixed models (GLMM) with a binomial distribution were employed to assess the effects of other variables on the proportion of bats’ choice preference using the R package lme4 (Bates et al., 2015). Models always included the bats’ choice preference (0/1) as the dependent variable; the flight latency, flight time, number of echolocation pulses, playback position, and playback side were assigned as fixed variables, and the number of experiments (rounds) and bat ID were assigned as random effects. All statistical analyses were conducted in R 3.5.0 (R Core Team, 2018).

Genomics sequencing, assembly, and annotation

Genomic DNA isolation and library preparation (PacBio, Illumina, Hi-C)

All the genomic sequences were generated by Novogene Inc, Beijing, China.

Phenol-chloroform extraction of genomic DNA

Several fresh samples of one adult male I. io individual muscle tissue were stored in liquid nitrogen. Snap-frozen muscle tissue was pulverized into a fine powder and then lysed overnight in high-salt tissue lysis buffer at 55°C. RNA was removed via RNase A for 1 hour at 37°C. DNA was extracted from muscle tissue via two extractions using Phenol-Chloroform-Isoamyl alcohol (25:24:1; molecular biology grade; RNase- and Protease-Free; PH = 8.0) solution into an aqueous solution. The aqueous solution was then extracted twice with chloroform-isoamyl alcohol (24:1; molecular biology grade). Finally, DNA was precipitated in ice-cold 100% ethanol. The precipitated DNA was washed three times with 80% ethanol. The precipitate was centrifuged using a tabletop centrifuge to remove the ethanol solution, and the filamentous DNA was collected, dried at room temperature for 20 min, and re-suspended in elution buffer (Tris-EDTA).

PacBio long insert library and Illumina short insert library preparation

We constructed both PacBio CLR long (20–40 kb) insert sizes and Illumina short (300 bp) insert size libraries. The libraries were contracted following manuals of SMRTbell Template Prep Kit 1.0 and Illumina TruSeq Nano DNA Library Prep Kits, respectively. Short insert libraries were sequenced using Illumina HiSeq X Ten platforms via the 150 bp paired-end method, and long insert libraries were sequenced using PacBio Sequel II platforms.

Hi-C conformation capture

DNA was extracted in a manner similar to that described above. After labeling the DNA ends with biotin and incubating at 37°C for 45 min, the enzyme was inactivated with 20% SDS. The T4 DNA ligase (NEB) was added to the DNA followed by a 4–6 h incubation at 16°C. Proteinase K was added for reverse cross-linking during incubation at 65°C overnight. Purified DNA fragments were dissolved in 86 μl of water. Unligated ends were then removed. After the purified DNA was fragmented to a size of 300–500 base pairs, the ends of the DNA were repaired. DNA fragments labeled by biotin were finally separated on Dynabeads® M−280 Streptavidin (Life Technologies). Hi-C libraries were controlled for quality and sequenced on an Illumina Hiseq X Ten platform.

Total RNA extraction

Fresh samples of this bat individual tissues (muscle, brain, ear, heart, kidney, stomach, liver, lung, spleen, and tongue) were placed in liquid nitrogen and pulverized into powder. All powdered tissues were lysed in TRIzol reagent for 5minat room temperature. Then, all RNA was extracted and precipitated using chloroform-isopropanol. The total RNA sample was dried for 8 min after the top solution had been removed. The quantity and quality of total RNA were measured using a Bioanalyzer 2100 Tapestation (Agilent Technologies, Santa Clara, CA).

RNA library preparation

All libraries with an insert size of 300 bp of RNA from individual tissues were constructed using an Illumina NEBNext Ultra RNA Library Prep Kit. The resulting 300 bp libraries were subsequently sequenced using Illumina noveseq 6000 instruments.

Genome assembly

PacBio subreads and Illumina paired-end reads de novo assembly

To assemble the high-quality genome of the bat, the reads from PacBio were assembled using wtdbg2 (Ruan and Li, 2020). The wtdbg2 assembles the subreads and generates a “ctg.lay.gz” file with the contig layout and edge sequence. Next, the information from the layout and edge sequences of the contigs was utilized by “wtpoa-cns” to form contig-level genomes. Racon (Vaser et al., 2017) was used to polish the contig-level genome by using raw SMRT reads. To further correct the errors, Illumina paired-end reads were utilized to enhance the quality of genome primary assembly. The paired-end clean reads from the Illumina platform were aligned to the contig-level genome using BWA (Li, 2013). The assembly was corrected for base-pair errors using Pilon (Walker et al., 2014) with the results of the BWA alignment to correct the errors of single insertions and deletions in homopolymer-enriched regions. We excluded contigs shorter than 10 kb to avoid obtaining incorrect assembly results.

Assembly genome assessment

BUSCO was used to assess the completeness of the assembled genome with 4104 single-copy conserved orthologs sets (mammalia_odb9) (Waterhouse et al., 2018). We also assessed the accuracy of the assembly genome by aligning it to 248 highly conserved protein sequences using CEGMA (Parra et al., 2007). To assess the accuracy of the assembled genome, high-quality Illumina reads were mapped to the assembled genome using the BWA method (Li, 2013). The mapping rate of reads and coverage of the genome were measured.

Hi-C scaffolding

A further step was taken to extend the primary assembly genome length by aligning the Hi-C Illumina reads with the primary assembly genome using BWA (Li, 2013). Alignment results were filtered with SAMtools following the Arima filtering pipeline (https://github.com/ArimaGenomics/mapping_pipeline) (Li et al., 2009). AllHiC was used to scaffold the high-quality alignments into an assembly genome at the chromosome level (Zhang et al., 2019). The chromosome-level assembly genome was regarded as the final genome of I. io.

Genome annotation

Repeat annotation

In the current analysis, a combined strategy using homology alignment and de novo approaches was used in the repeat annotation. In the homology alignment method, the assembly genome was analyzed using RepeatMasker (Tarailo-Graovac and Chen, 2009) and RepeatProteinMask (the in-house scripts of RepeatMasker). The repeats sequence library used in homology alignment was Repbase (Jurka, 2000). TRF (Benson, 1999) was used to extract tandem repeats from the genome by ab initio prediction. LTR-FINDER (Xu and Wang, 2007) and RepeatScout (Price et al., 2005) were used to create de novo repetitive elements databases. We removed repeat sequences that were shorter than 100 bp and had more than 5% N bases; the remaining sequences made up the de novo repetitive element library. We combined the de novo repetitive elements library and Repbase library to create a custom repeat sequences library. This library was used in RepeatMasker to identify DNA-level repeats and to create a masked genome.

Structure annotation

Three strategies were used to predict assembly genome structure: homology-based prediction, ab initio-based prediction, and RNA-seq based prediction.

Predictions based on homology

In homology-based prediction, the genomes used in the prediction followed two rules: (i) the species providing the genome should have a closer distance in a phylogenetic tree; (ii) the genome of a species should have higher quality. In homology-based prediction, therefore, six Chiropteran species were used, the big brown bat (Eptesicus fuscus), the little brown bat (Myotis lucifugus), Brandt’s bat (Myotis brandtii), David’s myotis (Myotis davidii), the Natal long-fingered bat (Miniopterus natalensis), and the great roundleaf bat (Hipposideros armiger). The homology protein sequences of six bats were downloaded from NCBI. Protein sequences were aligned to the masked assembly genome using tBLASTn (Camacho et al., 2009). The best blast hit was re-aligned to the masked assembly genome using GeneWise (Birney et al., 2004) to identify precise alignments.

Predictions based on RNA-seq

In this study, the transcriptome sequences of ten tissues (muscle, heart, liver, spleen, lung, kidney, stomach, tongue, brain, and ear) were assembled using Trinity (Grabherr et al., 2011) with the de novo mode. Furthermore, genome-guided-based assembly was also applied to optimize the results of the genome annotation. Two different methods were used in genome-guided assembly. We first assembled transcriptome sequences using Hisat2 (Kim et al., 2019) and Trinity in Trinity’s genome-guided mode. Furthermore, the transcriptome sequences were assembled following Hisat2 and Stringtie (Pertea et al., 2015) protocols. De novo and genome-guided assemblies generated by Trinity and Stringtie were used as inputs for PASA (Haas et al., 2008) for generating a comprehensive transcriptome database. TransDecoder (http://transdecoder.github.io/) was used to predict protein-coding regions based on the results of PASA.

Predictions based on ab initio

PASA was used to align the transcriptome assembly based on Trinity’s genome-guided assembled genome. The gff3 format result of PASA was converted to GeneBank format, and protein sequences with high identity (≥80%) were excluded. After filtration, the protein sequences were divided into two sets (a training set and a test set) for performing multiple rounds of training in Augustus (Keller et al., 2011). The alignments of RNA-seq in the previous analysis were extracted hints using “bam2hints” and “wig2hints.pl” scripts from Augustus. Augustus predicts the structure of the assembled genome based on the results of hints and training. SNAP can also be used to predict the structure of the assembled genome (Korf, 2004). Input sequences were obtained from the genome using the PASA results. Input sequences with high identity or errors were discarded. After filtration, sequences were broken into fragments with one gene per sequence; these could be up to 1000 bp on either side of a gene using the “fathom” of SNAP (Korf, 2004). The HMM file estimated by “fathom” was used for downstream analysis. The “trainGlimmerHMM” of GlimmerHMM (Majoros et al., 2004) was used to train the model using homologous protein sequences. Then, “glimmerhmm” was used to predict the structure of the genome using the trained model. Genescan (http://hollywood.mit.edu/GENSCAN.html) and Geneid (Blanco et al., 2003) were also used in prediction based on ab initio results. EvidenceModeler (Haas et al., 2008) was utilized to integrate the different results from different prediction strategies and to produce consensus results of assembled genome structure prediction. The weight values were set as:

“PROTEIN GeneWise 20;

TRANSCRIPT PASA 20;

ABINITIO_PREDICTION Augustus 5; ABINITIO_PREDICTION Genescan 1; ABINITION_PREDICTION Geneid 1; ABINITIO_PREDICTION SNAP 1; ABINITIO_PREDICTION GlimmerHMM 1;

OTHER_PREDICTION TransDecoder 100”.

Functional Annotation

The program blastp (Camacho et al., 2009) was used to analyze predicted sequences from the structure annotation using SwissProt and KEGG databases. The gene functions were assigned according to the best match. InterProScan (Zdobnov and Apweiler, 2001) was used to annotate the motifs and domains. Six publicly available databases were searched by InterProScan, including PANTHER, SMRT, Pfam, GO, and PROSITE. The GO term ID was assigned to the corresponding gene based on the annotation of InterPro. We also mapped the protein sequences set to the NR database with the threshold “E-value ≤ 1e-5”.

Non-coding RNA annotation

We used tRNAscan-SE (http://lowelab.ucsc.edu/tRNAscan-SE/) (Chan and Lowe, 2019) to predict tRNAs. The related species’ rRNA sequences were used for annotation of rRNA using BLAST. Other ncRNAs, including miRNAs and snRNAs, were identified by searching against the Rfam 4.1 databases with default parameters using the program Infernal (http://infernal.janelia.org/) (Nawrocki and Eddy, 2013).

Genome phylogenetic analysis

Identification of one-to-one orthologous proteins

To understand the adaptation of I. io to seasonal aerial-hawking of birds, we undertook a comparative analysis between I. io to other chiropteran species. To study the adaptations to dietary expansion in I. io, we selected species of Chiroptera with different diets and an available high-quality reference genome. A total of 16 public genomes of bats were used in the comparative analysis: the big brown bat (E. fuscus), Kuhl’s pipistrelle (Pipistrellus kuhlii), Brandt’s bat (M. brandtii), David’s myotis (M. davidii), the Natal long-fingered bat (M. natalensis), Pallas’s mastiff bat (Molossus molossus), the yellow-shouldered bat (Sturnira hondurensis), the Jamaican fruit bat (Artibeus jamaicensis), the pale spear-nosed bat (Phyllostomus discolor), the Indian flying fox (Pteropus giganteus), the large flying fox (Pteropus vampyrus), the Egyptian fruit bat (Rousettus aegyptiacus), the greater horseshoe bat (Rhinolophus ferrumequinum), the great roundleaf bat (H. armiger), the black flying fox (Pteropus alecto), and the common vampire bat (Desmodus rotundus). As outgroups, we used the domestic cat (Felis catus), domestic dog (Canis lupus familiaris), and human (Homo sapiens) (Table S6).

We used three different methods to identify one-to-one orthologous proteins. These consensus sets of the three methods were considered as final orthologous proteins. The protein sequences of 19 species were downloaded from NCBI, and the protein sequences of I. io were extracted from the results of the genome annotation. All protein sequences were filtered using in-house scripts, and only the longest protein sequences were retained. Then, the longest protein sequences were used as input for OrthoFinder (Emms and Kelly, 2019), OrthoMCL (Li et al., 2003), and Inparanoid (O'Brien et al., 2005) to identify one-to-one orthologous proteins. We used the R package ape to extract the topology of phylogenetic trees for 20 species used in OrthoFinder (Paradis et al., 2004).

Phylogenetic tree construction

The protein and corresponding gene sequences of orthologous proteins among the 20 species were used to perform multiple sequence alignment via PRANK (Löytynoja, 2014). The poorly aligned regions of protein sequences were removed using trimAl (Capella-Gutiérrez et al., 2009). The conserved aligned regions of gene sequences were extracted using Gblocks (Talavera and Castresana, 2007). The orthologous proteins and genes were concatenated in the order of first and last used for the construction of the phylogenetic tree, respectively. The program modelFinder (Kalyaanamoorthy et al., 2017) was used to determine the best models for proteins and genes. PartitionFinder (Lanfear et al., 2017) was used to determine the best model in partition mode. The phylogenetic tree was constructed using RAxML with 200 bootstrap replicates (Stamatakis, 2014). The RAxML commands for different sequences type are displayed as follows:

Protein sequences in the normal model:

“raxmlHPC-PTHREADS-AVX2-s [protein alignment file]-n [output]-m PROTGAMMAIHIVBF-f a-N 200-p 12345-x 12345”

Protein sequences in the partitioned model:

“raxmlHPC-PTHREADS-AVX2-s [protein alignment file]-n [output]-q [partition file]—m PROTGAMMAWAG-f a-N 200-p 12345-x 12345”

Gene sequences:

“raxmlHPC-PTHREADS-AVX2-s [protein alignment file]-n [output]-m GTRGAMMAI-f a-N 200-p 12345-x 12345”

1st+2nd codon:

“raxmlHPC-PTHREADS-AVX2-s [protein alignment file]-n [output]-m GTRGAMMAI-f a-N 200-p 12345-x 12345”

3rd codon:

“raxmlHPC-PTHREADS-AVX2-s [protein alignment file]-n [output]-m GTRGAMMA-f a-N 200-p 12345-x 12345”

Divergence time estimation

The programs r8s (Sanderson, 2003) and PAML (Yang, 2007) were used to estimate the divergence time of I. io. A total of 6,168,375 bases and six pieces of calibration-time evidence from TimeTree were used in r8s. Then, the fourfold degenerate synonymous (4D) sites were extracted from the orthologous genes using an in-house script. A total of 2,056,125 4D sites were used in MCMCtree analysis with the same calibration times evidence used in r8s. The MCMCtree was run for 1,000,000 steps, with samples drawn every 2,000 steps after a burn-in of 500,000 steps. The ESS score of MCMCtree was analyzed using Tracer in BEAST (Bouckaert et al., 2019), and a value greater than 200 was considered as a convergence of the stationary distribution.

The expansion and contraction of gene families

The gene families identified by OrthoFinder (Emms and Kelly, 2019) were used in an expansion/contraction analysis. In CAFÉ (Han et al., 2013), gene families with copy numbers less than 100 and at least one copy were analyzed, and a random birth-and-death model was used to determine the expansion and contraction of a gene family. In the present analysis, we used divergence times estimated by r8s and MCMCtree in CAFÉ (Han et al., 2013). The consensus results of gene families were extracted and assigned PANTHER IDs using hmmer (Johnson et al., 2010). The GO and KEGG enrichment analyses of expanded and contracted gene families were applied using KOBAS (Xie et al., 2011).

Identification of positively selected genes (PSGs)

The CodeML of PAML (Yang, 2007) was used to calculate the ratio of non-synonymous to synonymous substitutions (dN/dS). The PSGs were identified by CodeML with the branch-site model (model A vs a null model). I. io was set to be the foreground branch. Model A was set with parameters “model = 2; NSsites = 2; fix_kappa = 0; fix_omega = 0”; and the null model was set with parameters “model = 2; NSsites = 2; fix_kappa = 0; fix_omega = 1; omega = 1”. The likelihood ratio test (LRT) was used to compare the dN/dS of the branch of I. io between model A and the null model using chi2 in PAML (Yang, 2007) with parameter “chi2 1 |lnLmodeA - lnLmodelnull|∗2”. Genes with a P value less than 0.05 were treated as candidates for positive selection. The aBSREL model of HyPhy (Kosakovsky Pond et al., 2020) was also used to detect whether a gene had undergone positive selection. Genes with significant positive selection (Holm-Bonferroni corrected p value < 0.05) were treated as candidates for positive selection. Finally, PSGs were identified in the consensus results of CodeML and aBSREL. The positively selected sites (PSSs) were identified with the result of Bayes Empirical Bayes (BEB) analyses in PAML (Yang, 2007). The GO and KEGG enrichment analysis of PSGs were performed via KOBAS (Xie et al., 2011).

Identification of rapidly evolving genes (REGs)

CodeML was also used to analyze rapidly evolving genes based on the orthologous genes and phylogenetic trees used in the positive selection analysis. In this analysis, a null model and an alternative model were used. The former assumes that all branches evolve at the same rate, while the alternative model assumes that the foreground branch evolves at a different rate than background branches (Yang, 2007). I. io was set to be the foreground branch. The null model used parameters “model = 0; NSsites = 0; ncatG = 4; fix_kappa = 0; fix_omega = 0”; the alternative model used parameters “model = 2; NSsites = 0; ncatG = 4; fix_kappa = 0; fix_omega = 0”. The comparisons of both models of each gene were calculated using LRT with “df = 1”. Significant (p < 0.05) genes with higher foreground ω than background d branches were considered as rapidly evolving genes. The GO and KEGG enrichment analyses of REGs were performed with KOBAS (Xie et al., 2011).

Genomic convergence analysis

The same challenge faced by I. io and D. rotundus is the nutritional composition of their specialized diet, i.e., blood. Therefore, we assumed that some metabolism-related genes of D. rotundus and I. io underwent convergent evolution. Hypoxia is another challenge faced by bats and other high-altitude mammals. Thus, we analyzed convergent genes among I. io and high-altitude mammals. Ten species were selected for the convergence analysis, including some Chiroptera species previously analyzed as well as mammals found at different altitudes. Mammals from different altitudes were the wild yak (Bos mutus), the Tibetan antelope (Pantholops hodgsonii), the black-and-white snub-nosed monkey (Rhinopithecus bieti), the giant panda (Ailuropoda melanoleuca), the dog (Canis lupus familiaris) and humans (Homo sapiens). In this study, convergent amino acid sites were identified by the following criteria: (i) the convergent amino acid sites of D. rotundus and I. io should be identical; (ii) the ancestral amino acid sites of D. rotundus and I. io should not be identical; (iii) changes should have occurred between D. rotundus and inferred ancestor amino acid sites, I. io the same as D. rotundus. The ancestor sequence was reconstructed using CodeML of PAML with “NSsites = 3; ndata = 10; model = 3; aaRatefile = jones.dat; ncatG = 8” for each orthologous gene (Yang, 2007). The unreliable (posterior probability < 0.9) reconstructed ancestor amino acid sites were discarded (Wang et al., 2015). The methods developed by Zou and Zhang (Zou and Zhang (2015) were used to identify the genes with convergent amino acid sites. The Poisson test was used to identify the observed convergent sites whose number was significantly higher than expected under random substitution (Zou and Zhang, 2015). We used JTT-fgene and JTT-fsite models and selected the results under the JTT-fgene model because the JTT-fsite model was affected by the number of species and sequences analyzed (Zou and Zhang, 2015). The sites with p value less than 0.05 were retained.

Quantification and statistical analysis

We have described the data analysis of this study in detail in the “METHOD DETAILS” section. Statistical tests are presented as asterisks in Figure 1 and as p values in Figure 2. The asterisks are defined in the relevant figure legend, together with the name of the statistical test.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Deposited data
Raw and analyzed data	This paper	N/A
Experimental models: Organisms/strains
Great evenning bats, Ia io	Wild-caught	N/A
Mountain White-eye, Zosterops japonicus	Wild-caught	N/A
Software and algorithms
ArcGIS	Environmental Systems Research Institute, Inc.	Version 10.2
R	R Core Team	Version 3.5.0
Avisoft SASLab Pro	Avisoft Bioacoustics, Berlin, Germany	Version 5.2.07
PAML	https://doi.org/10.1093/molbev/msm088	http://abacus.gene.ucl.ac.uk/software/paml.html
Wtdbg2	https://doi.org/10.1038/s41592-019-0669-3	Version 2.5
Racon	https://doi.org/10.1101/gr.214270.116	Version 1.4.3
Pilon	https://doi.org/10.1371/journal.pone.0112963	Version 0.7
BWA	http://github.com/lh3/bwa	Version 1.22
BUSCO	https://doi.org/10.1093/molbev/msx319	Version 3.0.2
CEGMA	https://doi.org/10.1093/bioinformatics/btm071	Version 2.5
SAMtools	https://doi.org/10.1093/bioinformatics/btp352	Version 1.11
AllHiC	https://github.com/tangerzhang/ALLHiC	Version 0.9.8
RepeatMasker	https://doi.org/10.1002/0471250953.bi0410s25	Version 4.1.0
Repbase	https://www.girinst.org/repbase/	Version 20181026
TRF	https://doi.org/10.1093/nar/27.2.573	Version 4.09
LTR-FINDER	https://doi.org/10.1093/nar/gkm286	Version 1.0.6
RepeatScout	https://doi.org/10.1093/bioinformatics/bti1018	Version 1.0.5
tBLASTn	https://doi.org/10.1186/1471-2105-10-421	Version 2.2.26
GeneWise	https://doi.org/10.1101/gr.1865504	Version 2.4.1
Trinity	https://github.com/trinityrnaseq/trinityrnaseq	Version 2.1.1
Hisat2	https://doi.org/10.1038/nmeth.3317	Version 2.0.4
Stringtie	https://doi.org/10.1038/nbt.3122	Version 1.3.3
PASA	https://doi.org/10.1093/nar/gkg770	Version 2.2.0
TransDecoder	https://github.com/TransDecoder/TransDecoder/releases	Version 5.5.0
Augustus	https://doi.org/10.1002/cpbi.57	Version 3.2.3
SNAP	https://github.com/KorfLab/SNAP	Version 2013.11.29
GlimmerHMM	https://doi.org/10.1093/bioinformatics/bth315	Version 3.04
Genescan	http://hollywood.mit.edu/GENSCAN.html	Version 1.0
Geneid	https://doi.org/10.1002/0471250953.bi0403s00	Version 1.4
EvidenceModeler	https://doi.org/10.1186/gb-2008-9-1-r7	Version 1.1.1
blastp	https://doi.org/10.1186/1471-2105-10-421	Version 2.2.26
InterProScan	https://www.ebi.ac.uk/interpro/	Version 5.31
tRNAscan-SE	https://doi.org/10.1007/978-1-4939-9173-0_1	http://lowelab.ucsc.edu/tRNAscan-SE/
Rfam	http://infernal.janelia.org/	Version 4.1
OrthoFinder	https://doi.org/10.1186/s13059-019-1832-y	Version 2.5.2
OrthoMCL	https://doi.org/10.1101/gr.1224503	Version 2.0.9
Inparanoid	https://doi.org/10.1093/nar/gkm1020	Version 4.2
PRANK	https://doi.org/10.1007/978-1-62703-646-7_10	Version 140110
trimAl	https://doi.org/10.1093/bioinformatics/btp348	Version 1.2rev59
Gblocks	https://doi.org/10.1093/oxfordjournals.molbev.a026334	Version 0.91b
modelFinder	https://doi.org/10.1038/nmeth.4285	Version 1.6.12
PartitionFinder	https://doi.org/10.1093/molbev/msw260	Version 2.1.1
RAxML	https://doi.org/10.1093/bioinformatics/btu033	Version 8.2.12
R8s	https://doi.org/10.1093/bioinformatics/19.2.301	Version 1.81
BEAST	https://doi.org/10.1371/journal.pcbi.1006650	Version 2.6.4
CAFE	https://doi.org/10.1093/molbev/mst100	Version 4.2.1
hmmer	https://doi.org/10.1371/journal.pcbi.1002195	Version 3.3.2
KOBAS	https://doi.org/10.1093/nar/gkr483	Version 3.0
HyPhy	https://doi.org/10.1093/molbev/msz197	Version 2.5.31