https://orcid.org/0000-0001-5956-1695
Abstract
Bats possess a range of distinctive characteristics, including flight, echolocation, impressive longevity, and the ability to harbor various zoonotic pathogens. Additionally, they account for the second-highest species diversity among mammalian orders, yet their phylogenetic relationships and demographic history remain underexplored. Here, we generated de novo assembled genomes for 17 bat species and 2 of their mammalian relatives (the Amur hedgehog and Chinese mole shrew), with 12 genomes reaching chromosome-level assembly. Comparative genomics and ChIP-seq assays identified newly gained genomic regions in bats potentially linked to the regulation of gene activity and expression. Notably, some antiviral infection-related gene under positive selection exhibited the activity of suppressing cancer, evidencing the linkage between virus tolerance and cancer resistance in bats. By integrating published bat genome assemblies, phylogenetic reconstruction established the proximity of noctilionoid bats to vesper bats. Interestingly, we found 2 distinct patterns of ancient population dynamics in bats and population changes since the last glacial maximum does not reflect species phylogenetic relationships. These findings enriched our understanding of adaptive mechanisms and demographic history of bats.
Introduction
Bats are a diverse group of mammals found worldwide, occupying nearly every ecosystem except for polar regions. Their diversity peaks in tropical areas, and they exhibit remarkable variations in morphology, behavior, and life history characteristics (Vaughan et al. 2013). In addition, bats can host numerous zoonotic viruses, including viruses with high sequence similarity to pandemic-causing viruses, such as Middle East respiratory syndrome, severe acute respiratory syndrome coronavirus (SARS-CoV), and SARS-CoV-2 (Li et al. 2005; Letko et al. 2020; Zhou et al. 2020; Irving et al. 2021). Consequently, bats represent essential model systems for investigating adaptive radiations and the interplay between genetic and phenotypic traits, with a special focus on understanding the traits underlying antiviral and disease resistance. With the expansion of genomic resources for bats, pioneering research has uncovered key facets of bat biology, including echolocation abilities (Parker et al. 2013; Dong et al. 2017), flight mechanics (Zhang et al. 2013; Eckalbar et al. 2016), diet variations (Zepeda Mendoza et al. 2018; Wang et al. 2020), antiviral defenses (Pavlovich et al. 2018; Jebb et al. 2020), and extreme longevity (Seim et al. 2013). Alongside these findings, there has been substantial progress in bat genomics recently, with numerous genomes assembled and some even reaching chromosome-level resolution (Jebb et al. 2020; Blumer et al. 2022; Scheben et al. 2023; Tian et al. 2023).
However, there is ample scope for further discovery by incorporating a broader array of bat species because they also represent the second-largest mammalian order, with >1,400 species. This diverse group includes 21 families, distributed across 2 suborders: Yinpterochiroptera and Yangochiroptera (Wilson and Mittermeier 2019). Key questions which would benefit from the inclusion of more species in analysis would include understanding what constitutes the expanded components when genome size has contracted overall. Additionally, the exploration into bats' potential resistance to tumors is an emerging field of interest. Comparative genomic analysis has revealed that certain tumor suppressor genes show signs of positive selection in bats (Scheben et al. 2023). Experimental work, such as the activation of oncogenes in the primary cells of 7 bat species, has presented evidence of cancer resistance, particularly in Rickett's big-footed bat (Myotis pilosus) (Hua et al. 2024). Aside from these adaptive features, the demographics and phylogenetics of bats remain relatively underexplored, highlighting the necessity for more comprehensive phylogenetic investigations and historical reconstructions.
To address these knowledge gaps, we performed de novo genome assembly for 17 bat species, representing 15 genera and 6 families. By integrating these new genomes with 18 previously published bat genomes, we constructed a time-calibrated phylogenetic tree, analyzed genomic components and gene family dynamics, and explored positively selected genes and species population histories. Overall, our findings offer a deeper understanding of bat evolution and adaptation, as well as providing new insights into genetic conservatism and selection in different lineages.
Results and Discussion
Sequencing and Annotation of 17 Novel Bat Genomes
We performed de novo assembly of 17 bat genomes using single-tube long fragment read (stLFR) technology (Wang et al. 2019), effectively doubling the available genomic resources alongside 18 previously published bat genomes (up to March 2022). Additionally, we assembled genomes from 2 outgroups: the Amur hedgehog (Erinaceus amurensis) and the Chinese mole shrew (Anourosorex squamipes) (Fig. 1a; supplementary tables S1 to S3, Supplementary Material online), which also fall in Laurasiatheria. The initial assembly size of the bat genomes was 1.86 to 2.39 Gb (supplementary table S4, Supplementary Material online), consistent with k-mer-based estimations (supplementary table S5, Supplementary Material online) and the contig N50 ranged from 23.45 to 170.84 kb in the 17 bat species. Using Hi-C data and the Juicer and 3D-DNA pipeline, we anchored 10 bat and 2 outgroup genomes to pseudochromosomes (Durand et al. 2016; Dudchenko et al. 2017) (Fig. 1b), leading to an improvement of at least 4.59-fold in the assembly N50 compared with the initial draft scaffold level for each bat species (supplementary table S4, Supplementary Material online). The average N50 of the chromosome-level assembly reached 121.60 Mb, surpassing the short-read assembly of our self-sequenced bats (9.41 Mb) (supplementary fig. S1, Supplementary Material online). Totally, approximately 83.60% of the genome across 10 bat species achieved chromosome-level assembly. Evaluation of these genomes using Benchmarking Universal Single-Copy Orthologs (BUSCO) highlighted that, on average, 86.68% of genes were complete, as assessed by conserved mammalian genes (supplementary fig. S2 and table S4, Supplementary Material online). Although the sequence completeness did not meet the high benchmarks set by long-read sequencing platforms, most of our annotations were comparable to those reported in recent bat studies (Moreno et al. 2021). Additionally, all 17 newly generated bat genomes displayed greater than 94.92% intactness for 197 nonexonic ultraconserved elements (Bejerano et al. 2004) (supplementary fig. S3, Supplementary Material online).
The divergence time, assembly strategy, genome size evolution, and newly gain DNA in bats. a) The maximum likelihood phylogenetic topology from 4,213 orthologous genes of 35 bat species and 4 outgroups. Newly sequenced 17 bat species in this study are indicated by bold font, with “*” denoting sequencing using stLFR and Hi-C reads (n = 10) and “#” denoting sequencing based on stLFR reads (n = 7). All nodes in the tree have 100% support under the 200 bootstraps. Divergence times were estimated using the 4d sites and 4 fossil calibrations. Light blue bars along the branches represent 95% highest posterior density intervals of divergence time for internodes. The bat superfamily is depicted by 5 different colors on the tree branches. Nine families are distinguished by the vertical gray bars with a head image of a representative bat. b) Genome assembly pipeline methods used in 2 hierarchical assembly processes. The Hi-C maps (right part) of T. melanopogon before (left below) and after (left above) manual calibration, visualized using Juicebox software to show contact density. c) Comparison of the average size of introns, exons, and intergenic regions across avian, reptilian, chiropteran, and flightless mammalian genomes. d) DNA loss and gain in Chiroptera, Yinpterochiroptera, and Yangochiroptera ancestor. The size of TEs and intergenic and gene regions show on the left and the corresponding contents of different TE types show on the right, e.g. DNA, LINE, LTR, RC, and SINE. e) A 800 bp fragment of bat DNA gain compared to human genome (top), located in PI4K2A Intron1 of reference genome R. pusillus (center). The vertical orange box labels the region of this DNA gain (bottom). The H3K27ac of ChIP-seq tracks are represented in 10 bp intervals (light blue peak, bottom), and enhancer marks are indicated by light blue boxes below the x axis (homologous coordinate in R. pusillus genome) for liver and kidney of the 3 bats. f) Above: evaluation of promoter activity for 3 fragments using dual-luciferase reporter assay. The 2 bat fragments contained special DNA gain and were inserted into the 1.2 kb control fragment of human genome. Data are the means ± SD of 3 independent experiments. Significance values for each pairwise comparison were calculated using the Mann–Whitney text. *P < 0.05; ***P < 0.001. Bottom: gene expression of PI4K2A by RNA-seq in the liver and kidney of 36 bats (orange) and 70 flightless mammals (FLM, blue). **P represents P < 0.01 calculated by the Wilcoxon test. The bat picture in b) and the liver and kidney picture in f) were from BioRender (https://www.biorender.com/).
After repeat masking, 3 methods (transcriptome sequencing data, ab-initio, and homology annotation) were used to predict protein-coding genes in bat genomes (supplementary tables S6 and S7, Supplementary Material online). The number of protein-coding genes in the 17 bat species was 17,764 to 20,978 (supplementary table S7, Supplementary Material online), which appeared to be irrespective of assembly metrics (supplementary fig. S3, Supplementary Material online). To investigate gene compaction in bats, we compared gene structure among reptiles (n = 3), birds (n = 3), flightless mammals (n = 9), and the newly sequenced bats (n = 17) (supplementary tables S8 and S9, Supplementary Material online). On average, bats exhibited gene lengths that were 19.45% shorter than those of flightless mammals and 18.87% shorter than those of reptiles (Fig. 1c). Specifically, the intron length in bats (4.19 kb) was reduced compared to flightless mammals (4.96 kb) and reptiles (4.50 kb), supporting the notion that volant species tend to have shorter introns than their non-flying counterparts (Zhang and Edwards 2012). Notably, no significant differences in intron size were found between the Yangochiroptera and Yinpterochiroptera suborders.
In examining the repeat elements within bat genomes, the Yangochiroptera suborder showed a heightened accumulation of hAT elements within the last 30 million years. On the other hand, the Yinpterochiroptera suborder generally exhibited an earlier surge in accumulation, occurring around 50 to 75 million years ago (Mya), with the exception of great roundleaf bat (Hipposideros armiger) (supplementary figs. S4 to S6 and table S10, Supplementary Material online). Intriguingly, when compared to Yangochiroptera, 3 species within the Rhinolophoidea group demonstrated a substantial uptick in TcMar-Mariner transposons over the last 50 million years (supplementary fig. S5, Supplementary Material online). Thus, the timing of insertion for certain DNA transposons shows marked variability between these 2 bat suborders.
Newly Gained Genome Regions Have Been Recruited for Gene Regulation
To investigate the regions of genomic variation characterized by DNA gain and loss in bats, genomes of 9 representative bats from distinct families were selected and underwent multispecies alignments with 6 other mammalian species (supplementary fig. S7 and table S11, Supplementary Material online). We observed a 13.19-fold DNA depletion (deletion/insertion rate) with 145.56 and 11.03 Mb of DNA loss and gain in the bats' ancestral lineage, respectively, indicating that the bat ancestor underwent substantial genomic shrinkage rather than expansion (Fig. 1d). Additionally, individual species experienced diverse DNA loss, ranging from 744.28 Mb in the great roundleaf bat to 1039.94 Mb in the brown long-eared bat (Plecotus auritus) (supplementary figs. S8 to S10, Supplementary Material online). In addition to DNA loss, transposable elements played a significant role in the acquisition of new DNA in bats. The expansion of transposable element lengths within the Chiroptera, Yangochiroptera, and Yinpterochiroptera lineages was measured at 7.42, 43.50, and 24.03 Mb, respectively (Fig. 1d). Notably, LINEs were the predominant type of insertion during this expansion (Fig. 1d).
To explore the functional implications of these expanded regions, we performed ChIP-seq assays for histone modifications (H3K4me3 and H3K27ac) in representative bat species (supplementary table S12, Supplementary Material online). Intersection analysis of gained DNA and epigenetic modification regions revealed that ∼4.63% of bat-specific DNA gains represented potential active promoter or enhancer regions. Correspondingly, approximately 371 and 878 genes may be regulated by these potential active promoters and enhancers, respectively, enriching essential metabolic processes, including RNA processing, nucleic acid binding, Rho protein signal transduction, cation transport, and hydrolase activity (supplementary tables S13 and S14, Supplementary Material online). These newly gained potential promoters and enhancers may promote the expression of corresponding genes. For example, an 800 bp insertion exclusively occurred in intron1 of PI4K2A in bats (Fig. 1e), consisting of 2 endogenous retrovirus (ERV) elements. Dual-luciferase reporter assays showed that plasmids constructed using 2 bat insertion sequences exhibited increased luciferase activity relative to the control group and plasmids containing human sequences (Fig. 1f; supplementary table S15, Supplementary Material online). Accompanying these, both H3K27ac and H3K4me3 histone modifications assay suggested a promoter activity for this newly acquired region (supplementary table S16, Supplementary Material online). Consistently, the gene PI4K2A was found to be more highly expressed in bats than in other mammals (Fig. 1f; supplementary tables S17 and S18, Supplementary Material online). Notably, high levels of PI4K2A expression are associated with tumorigenesis and play a role in various signaling pathways (Huang et al. 2023). Considering that bats are known for their resistance to tumors, further research is necessary to elucidate the role of this newly acquired genomic region.
Three other bat-specific gained regions are identified with enhancer signals and coincided with the presence of ERV elements. Consistently, their nearby genes, namely, SLC25A3, SDHB, and DHTKD1, exhibited higher expression levels in the liver or kidneys of bats compared with other mammals (supplementary fig. S11, Supplementary Material online) and were related to mitochondrial energy metabolism (Mayr et al. 2007; Xu et al. 2013; Goncalves et al. 2021). For example, SLC25A3 encodes a mitochondrial copper transporter that is essential for the biogenesis of cytochrome c oxidase (Boulet et al. 2018). Dysfunction of SLC25A3 has been linked to impaired mitochondrial ATP synthesis (Peoples et al. 2021). Given that high-speed flight requires efficient energy production, these elevated gene expressions suggest that they represent adaptations to flight in bats.
Improved Resolution of Bat Family Phylogeny
The use of molecular approaches has led to fundamental changes within Chiroptera classification, dividing them into 2 groups: Yangochiroptera and Yinpterochiroptera, in contrast to previous analysis which assumed that Pteropodidae (known as megabats) were distinctly clustered from all other bats (formerly microbats) (Amador et al. 2018). However, the family level relationships within Yangochiroptera, specifically among Emballonuridae (sheath-tailed bats), Noctilionoidea (noctilionoid bats), and Vespertilionoidea (vesper bats), remain uncertain. Currently, 3 debated hypotheses exist concerning these groups: (i) the sheath-tailed bats hold the most basal position, while the noctilionoid and vesper are sister taxa (topology T1). This topology garners support from multidisciplinary research, which includes the morphological characteristics (Springer et al. 2001) and nuclear gene sequences (Van Den Bussche and Hoofer 2004; Meredith et al. 2011). (ii) Sheath-tailed bats served the sister clade to noctilionoid bats (topology T2), with this relationship being primarily inferred from DNA sequence analysis (Teeling et al. 2000, 2002; Eick et al. 2005). (iii) Topology T3 contends that sheath-tailed and vesper bats are sister groups, a hypothesis that is upheld by gene tree analyses alone (Upham et al. 2019). In this study, strong support (bootstrap value [BS] = 100%) for the close relationship between noctilionoid and vesper bats was observed in maximum likelihood trees for most data sets, including whole-genome alignment (∼446 Mb) (supplementary fig. S12, Supplementary Material online), orthologous genes (4,213 genes of 39 species and 4,223 genes of 40 species; supplementary figs. S13 to S15, Supplementary Material online), and CNEs (∼14.68 Mb) (supplementary figs. S16 and S17, Supplementary Material online). This topology remained unaffected by variation in substitution rates among codon positions (supplementary figs. S13 to S15, Supplementary Material online), substitution saturation (supplementary fig. S18 and tables S19 to S23, Supplementary Material online), or specific genes with the highest log-likelihood (supplementary figs. S19 and S20, Supplementary Material online). Furthermore, the relative frequency of gene trees in these data sets and concordance factors consistently supported noctilionoid bats as a sister group to vesper bats (Fig. 2a and b; supplementary figs. S21 and S22, Supplementary Material online).
Improved resolution and quantification of introgression among bat phylogeny. a) The quartet gene trees of conflict superfamily relationship within Yangochiroptera from ASTRAL and ILS detected by QuIBL software. The percent on the quartet gene tree is the relative frequency of trees supporting this topology, and the percent above was ILS (yellow) and introgression (gray, but the percent not shown) contributed to this topology. Noc is an abbreviation of Noctilioidea; Ves is Vespertilionoidea; Emb is Emballonuroidea; Yin is Yinpterochiroptera. b and c) The relative frequency of gene trees in b) supports the topology shown in c). The bar color in b) corresponds to the pie color in c). The outside cycle above tree c) is the site concordance factor (sCF), and inner cycle is gene concordance factor (gCF). Pteropodidae ER represents Eonycteris spelaea and Rousettus aegyptiacus. Pteropodidae P represents Pteropus vampyrus and Pteropus alecto. Pteropodidae C represents Cynopteru sphix. d) The introgression detected on the superfamily (left) and family (right) by Dsuite software. The tree along the y axis is in an “expanded” form, including the internal branches (the dotted line), as the ancestral branch of the descendant branch. The lineages adapt on the x axis. The gray box indicates no information on f-branch values. e) The introgression schematic diagram detected by HyDe on the superfamily (left) and family (right) levels within Yangochiroptera. Noc is an abbreviation of Noctilioidea; Ves is Vespertilionoidea; Emb is Emballonuroidea. Vesp is an abbreviation of Vespertilionidae; Mini is Miniopteridae; Molo is Molossodea; Phyl is Phyllostomidea; Emba is Emballonuridea.
A coalescent-based “species tree” analysis and the relative frequency of quartet gene trees revealed inconsistent topologies across different data sets (Fig. 2b and c; supplementary figs. S21 to S27, Supplementary Material online), indicating a role for incomplete lineage sorting (ILS) in these conflicts. However, CoalHMM analyses (Dutheil et al. 2009) revealed that ILS alone could not fully explain these conflicts (100% supported T3 topology). A further examination of the whole-genome alignments showed that ILS accounted for 88.36% of the conflicting relationships (Fig. 2a; supplementary tables S24 and S25, Supplementary Material online), with gene introgression accounting for the remaining portion, supported by evidence of introgression events between sheath-tailed and vesper bats (Fig. 2d and e; supplementary fig. S28 and tables S26 and S27, Supplementary Material online). Thus, ILS primarily contributed to conflicting relationships in the superfamily, complemented by introgression.
Another contentious topology in previous studies was the position of the greater short-nosed fruit bat (Cynopterus sphinx) within the Pteropodidae family. In our study, the maximum likelihood tree of coding alignment strongly supported Cynopterinae as the sister group of Rousettinae (BScd12-39 = BScd123-39 = BScd12-40 = BScd123-40 =100%; supplementary figs. S14 to S17, Supplementary Material online). Furthermore, approximately unbiased tests of coding alignments demonstrated significant divergence from alternative topologies (supplementary table S28, Supplementary Material online). Conversely, analyses of CNEs, mitochondrial genomes, and 4 degenerated (4d) sites supported the divergence of Cynopterinae from other Pteropodidae bats (BSCNE-39 = BSCNE-40 = BSmito = BS4d-39 = BS4d-40 = 100%; supplementary figs. S16, S17, and S29 to S31, Supplementary Material online). This relationship aligned with coalescent-based “species tree” analyses (local posterior probability [LPP], LPPcd123-39 = LPPcd123-40 = LPPCNE-39 = LPPCNE-40 = 1.0; BScd123-39 = BScd123-40 = BSCNE-39 = BSCNE-40= 100%; supplementary figs. S24 to S27, Supplementary Material online], the relative frequency of quartet gene trees (Fig. 3c; supplementary figs. S21 and S22, Supplementary Material online), and CFs (Fig. 2b and c; supplementary fig. S32, Supplementary Material online). This outcome contradicted a recent genomic study on Pteropodidae family relationships (Nesi et al. 2021), possibly due to the limited sample size of Pteropodidae in our study.
Adaptive genetic changes in bats. a) The blue shades indicated amino acid sites in MAVS, RNase L, and SLC39A4 under positive selection in the Chiroptera ancestor. b) Diagram showed positively selected genes MAVS and RNase L act in response to virus infection and cancer resistance. c) Migration abilities of RNase L in A549 cells. Histogram showed the migrated cell number per field with example pictures below. The first column is pEGFP control. The following 3 columns show cells with Mus musculus type (M-RnaseL) and 1-point mutated types from M. musculus type to Taphozous melanopogon type (M-RnaseL-D556N and M-RnaseL-D561N), respectively. The last 3 columns are shown T. melanopogon types (T-RnaseL) and 1-point mutated types from T. melanopogon type to M. musculus type (T-RnaseL-N556D and T-RnaseL-N561D), respectively.
Selection Signatures on Antiviral- and Antitumor-Related Genes
Under the robust phylogeny reconstructed above, we identified 111 genes undergoing positive selection in the ancestral branch of bats in PAML (supplementary tables S29 and S30, Supplementary Material online). Within this set, 51 PSGs also showed signatures of positive selection occurring at 1 site or more on the branch using the BUESED method (Murrell et al. 2015). The 111 PSGs were significantly enriched in 167 Gene Ontology (GO) terms (supplementary table S31, Supplementary Material online), including response to the virus (GO:0009615; P = 1.47 × 10−3; NCBP3, OASL, AGBL4, RNASE6, IFNGR1, RNase L, and MAVS), as well as KEGG pathways (supplementary table S32, Supplementary Material online) related to viruses or viral infection, including influenza A (P = 2.13 × 10−4), Epstein–Barr virus infection (P = 3.52 × 10−3), hepatitis C (P = 8.26 × 10−3), and human immunodeficiency virus 1 infection (P = 2.27 × 10−2) (supplementary fig. S33, Supplementary Material online). Evidence of positive selection was also observed in immunity-related genes in bats, including TLR4, NR1I2, FAS, CXCL16 (Jebb et al. 2020), and IFNGR1 (supplementary tables S31 and S32, Supplementary Material online). Importantly, 1 of the PSGs, Rnase L, an endoribonuclease, is induced by interferon (IFN) activation in response to RNA virus infections (Fig. 3a and b). This gene is known to be under positive selection in the Egyptian fruit bat (Egyptian rousette) (Pavlovich et al. 2018). Activation of Rnase L can trigger autophagy (Chakrabarti et al. 2012) and apoptosis (Zhou et al. 1997), both of which contribute to its antiviral effects. The direct antiviral effects of Rnase L involve the degradation of viral and cellular RNA, generating RNA cleavage products that stimulate the formation of the NLRP3 complex and the mitochondrial adaptor protein MAVS (Chakrabarti et al. 2015), which also underwent positive selection in bats (Fig. 3b). Notably, the NLRP3 inflammasome, which plays a role in inflammatory response, is known to be significantly dampened in bat primary immune cells compared with human or mouse counterparts, which has implications for the longevity and unique viral reservoir status in bats (Ahn et al. 2019). Furthermore, the viral double-stranded RNA recognition receptors RIG-I and MDA5 undergo a conformational change to promote binding to MAVS (Oshiumi et al. 2010; Hou et al. 2011; Peisley et al. 2014). This interaction triggers the nuclear translocation of NF-κB and interferon regulatory factor 3 pathways (Fig. 3b). Meanwhile, the antiviral Rnase L exhibited antitumor activities through its cell cycle arrest and apoptosis functions, which are induced by IFNs (Fig. 3b) (Sangfelt et al. 2000; Liang et al. 2006; Bisbal and Silverman 2007; Andersen et al. 2009). Individuals with heterozygous and homozygous RnaseL mutations (R462Q) have been shown to exhibit a 150% and >2-fold increased risk of prostate cancer, respectively (Carpten et al. 2002; Casey et al. 2002; Rökman et al. 2002; Wang et al. 2002). Thus, we used a transwell assay to explore the potential cancer-suppressing ability of bat-specific residues at Rnase L. Specifically, we used Rnase L gene sequences from the black-bearded tomb bat (Bat-RnaseL) and mouse (Mouse-RnaseL) as representatives and express them in A549 cells, respectively. The number of migrated cells expressing Bat-RnaseL was significantly lower than that of Mouse-RnaseL and control cells (A549 cells that express PEGFP) (P < 0.001; Fig. 3c). If we mutated the positively selected sites in Bat-RnaseL to mouse phenotype (i.e. Bat-RnaseL-N556D and Bat-RnaseL-N561D), cells exhibited increased migration relative to unmutated Bat-RnaseL cells (P = 4.99 × 10−2; P = 8.00 × 10−3; Fig. 3c; supplementary fig. S34, Supplementary Material online). Conversely, if we mutated the positively selected sites in Mouse-RnaseL (i.e. Mouse-RnaseL-D556N and Mouse-RnaseL-D561N), cells displayed decreased migration compared with unmutated Mouse-RnaseL cells (P = 4.60 × 10−2; P = 3.00 × 10−3; Fig. 3c; supplementary fig. S34, Supplementary Material online). Therefore, residues at positively selected sites enhance the ability of the RNase L gene in bats to suppress cancer cell migration, potentially contributing to cancer resistance.
Additionally, 4 PSGs are KEGG annotated into PD-L1 expression and PD-1 checkpoint pathway (P = 1.06 × 10−3) (IFNGR1, PDCD1, CD3G, and TLR4) (supplementary fig. S33, Supplementary Material online). ANP32A, SLC39A4, FAS, COPS2, EMP1, and CDH17 were related to tumor suppressor in the Tumor Suppressor Gene database. Therefore, The PSGs TLR4, NR1I2, FAS, IFNGR1, and Rnase L not only participated in antiviral response but also related to cancer resistance. Collectively, this evidence indicates a linkage between virus tolerance and cancer resistance in bats. This finding is particularly intriguing given that recent comparative genomic analyses across multiple bat species have identified signatures of positive selection in numerous cancer-related genes and receptors for coronaviruses in the ancestral line of bats (Scheben et al. 2023; MacDonald et al. 2024). In parallel, there appears to be species-specific or lineage-specific genetic responses to cancer resistance in bats. For instance, the TP53 gene has undergone selection in the David’s myotis (Myotis davidii), while the Rickett's big-footed bat (M. pilosus) shows downregulation of HIF1A, COPS5, and RPS3 (Zhang et al. 2013). These findings suggest that the functional pleiotropy of positively selected genes along the ancestral branch of bats, along with distinct evolutionary strategies, has both contributed to their cancer resistance.
Contrast Patterns of Ancient Population Dynamics in Bats
The history of effective population size (Ne) over time was reconstructed using pairwise sequential Markovian coalescent (PSMC) analysis (Li and Durbin 2011) using the substitution rate (supplementary table S33, Supplementary Material online) from the divergence time estimation with fossil calibration (supplementary table S34, Supplementary Material online). These demographic profiles exhibited substantial variation in terms of population dynamic trajectories (supplementary fig. S35, Supplementary Material online). However, this variability did not align with the phylogenetic relationships among species (Fig. 4a), suggesting that closely related taxa did not share similar demographic histories. We also conducted generalized additive models (GAMs) to examine potential correlations between effective population size and surface air temperature or global sea level (supplementary figs. S36 and S37, Supplementary Material online). The historical demographies and their associations with climate fluctuations display several contrasting patterns. For example, 6 bat species experienced only 1 population growth event (Pattern 1; Fig. 4b; supplementary fig. S35, Supplementary Material online) including Barbastella leucomelas, Rhinolophus pusillus, Chaerephon plicatus, Scotomanes ornatus, Harpiocephalus harpia, and S. kuhlii. In contrast, another 6 bat species exhibited 2 periods of population expansion, including C. sphinx, H. armiger, Hipposideros larvatus, Tylonycteris fulvida, P. javanicus, and M. armiger (Pattern 2; Fig. 4b; supplementary fig. S35, Supplementary Material online). Notably, these species display very different life histories, especially with regard to colony size, roosting type, and the degree of cosmopolitan associations within roosts. In general, bats exhibited increased Ne since 0.2 Mya and reached their Ne peaks at approximately 0.05 to 0.30 Mya during the middle-to-late Pleistocene (Cohen et al. 2013), i.e. the beginning of the Eemian interglacial period (also known as the last interglacial, around 0.12 Mya). The Eemian climate was warmer than the current Holocene, with Arctic temperatures about 2°C to 4°C higher than current temperatures (Nikolova et al. 2013; Schüpbach et al. 2018). The warmer climate would have provided more food and new habitats, and the flying ability of bats would have facilitated colonization, contributing to population growth. Fossil records from the northern Caucasus support the notion that the Eemian interglacial period was a time of rich bat species diversity (Rossina et al. 2006).
The population histories of bats. a) The normalized Ne of bats shows an inconsistent pattern of the demographic history with the phylogenetic relationships. Each species' Ne was reconstructed using PSMC (Li and Durbin 2011). b) The line depicts Ne fluctuations for Barbastella leucomelas, which has 1 expansion (top) and C. sphinx which has 2 expansions (bottom). The dotted line shows the surface air temperature, and the black dashed line shows the global sea level (Bintanja et al. 2005). The left y axis corresponds to Ne, and the right y axis corresponds to surface air temperature or global sea level. The shade represents the last glacial maximum.
Four other species displayed a fluctuating demographic pattern, including P auritus, Kerivoula hardwickii, Murina hilgendorfi, and Taphozous melanopogon (supplementary fig. S35, Supplementary Material online). Notably, 3 of these species have smaller body masses (4.55 to 10.87 g) compared with the median body mass of the other bat species for which body mass data were available (n = 10, median = 19.88, IQR = 11.38 to 37.99), indicating that bat populations with smaller body sizes may be more vulnerable (Carstens et al. 2018). Interestingly, the large myotis (Myotis chinensis) exhibited a Ne peak around 9 to 10 Mya, followed by minor fluctuations between 0.03 and 3 Mya, and later experienced a sharp decline (supplementary fig. S35, Supplementary Material online). This pattern contrasted with the patterns of other Myotis bats (e.g. Myotis brandtii and Myotis rufoniger) (Bhak et al. 2017). Therefore, M. chinensis requires special attention and conservation efforts owing to its relatively low population size from 3 to 0.01 Mya, followed by a dramatic decline at 0.03 Mya.
Our study uncovered that the majority of the species examined (15 out of 17) followed comparable demographic patterns, marked by a pronounced reduction in Ne during the last glacial maximum. This contraction reached its nadir approximately 0.01 to 0.03 Mya, corroborating findings from a prior investigation of 11 bat species that also reported minimal Ne during the same period (Chattopadhyay et al. 2019). This decrease in effective population size could be attributed to the cold and arid climate during the last glacial maximum, which severely impacted bat habitat and food sources. Notably, 2 exceptions to this pattern were Hilgendorf's tube-nosed bat (M. hilgendorfi) and brown long-eared bat (P. auritus), the populations of which reached their historical peaks during this period. These 2 species are distributed in regions approximately 10° farther north (>35° N) compared with other bats, which are typically found around the Tropic of Cancer and the equator. This suggests that glaciation had little influence on the population dynamics of bats adapted to high-latitude and low-temperature environments. Additionally, the recorded elevation of Hilgendorf's tube-nosed bat reaches up to 4,000 m, further supporting its adaptation to cold habitats (Wilson and Mittermeier 2019). Overall, it becomes evident that the ancestral dynamics of bats display a significant degree of diversity, influenced by both climatic conditions and their physiological capabilities.
Conclusions
This research advances our knowledge of chiropteran evolution through the introduction of 17 newly assembled bat genomes, marking the most considerable expansion in available bat genomic data to date. We uncover that bats have uniquely gained DNA sequences with apparent enhancer capabilities. Intriguingly, our work suggests an intriguing interplay between bats' innate antiviral defense and their resistance to cancer. Furthermore, we have clarified the phylogenetic placement of the sheath-tailed bats within the Yangochiroptera, resolving several previous uncertainties. Population dynamic analyses conducted herein demonstrate diverse evolutionary patterns among species, underlining the extensive variation in bat evolutionary lineage. However, our study acknowledges certain limitations, including incomplete taxonomic sampling, notably from groups like the sucker-footed bats, which might affect the precision of our phylogenetic inferences and the detection of gene flow signals. Additionally, only 10 of the 17 genomes achieved chromosome-level resolution. Future studies should thus focus on generating complete genomic assemblies for a wider selection of bat species to deepen our understanding of their diversification, evolution, and ecological adaptations.
Materials and Methods
The following materials of methods were briefly discribed, and the detailed information is in the Supplementary material.
Samples and DNA Sequencing
Genomic DNA was extracted from 17 bat species and 2 outgroups (the Amur hedgehog and Chinese mole shrew) to construct the stLFR library (Wang et al. 2019). Notably, 10 of these bats were used to construct Hi-C libraries. Subsequently, all samples were sequenced on BGISEQ 500 and DNBSEQ T1 platforms. All animal care and research protocols used in this study were approved by the Institute of Zoology, Chinese Academy of Sciences.
Genome Assembly
All samples were assembled using the stlfr2supernova pipeline (https://github.com/BGI-Qingdao/stlfr2supernova_pipeline). Notably, draft scaffolds of 10 bats were further anchored at the chromosome level using the 3D-DNA v170123 pipeline (Dudchenko et al. 2017). The karyotype of species was integrated as an input parameter for clustering (supplementary table S35, Supplementary Material online). Genome size and heterozygosity were estimated using GenomeScope v1.0 (Vurture et al. 2017) (supplementary fig. S38, Supplementary Material online), and the assembly’s completeness was evaluated via analyses of GC depth distribution (supplementary fig. S39, Supplementary Material online), read coverage (supplementary table S36, Supplementary Material online), ultraconserved elements, and BUSCO v3 (Simão et al. 2015). Additionally, the mitochondrial genome for all samples was assembled and annotated using NOVOPlasty v4.2 (Dierckxsens et al. 2017) and MitoZ v2.4 (Meng et al. 2019) (supplementary table S37, Supplementary Material online).
Genome Annotation
Three methods were used for protein-coding gene prediction: (i) ab initio gene prediction via Augustus v3.1.0 (Stanke et al. 2008) and Genscan v1.0 (Burge and Karlin 1997); (ii) homology annotation through GeneWise v2.4.1 (Birney et al. 2004); and (iii) transcript annotation using HISAT2 v2.1.0 and StringTie v1.3.4d (Pertea et al. 2016). Glean (Elsik et al. 2007) was used to integrate these 3 pieces of evidence to generate the final protein-coding gene sets. Noncoding RNA was annotated in bat assemblies using tRNAscan-SE v1.3.1 (Chan and Lowe 2019) and INFERNAL v1.1.1 (Nawrocki and Eddy 2013) (supplementary table S38, Supplementary Material online). RepeatMasker v4.0.9 (Smit et al. 2015) was used for repeat annotation, whereas RepeatClassifier v2.0.1 and TEclass v2.1.3 (Abrusán et al. 2009) were employed to predict TE types.
Genome Size Variation and Identification of HCEs
We explored genome size variation from 3 perspectives: bat-specific regions of DNA loss, DNA gain, and HCEs. We created whole-genome alignments for 9 representative bats and 7 FLMs using Lastz v1.04.03 (Schwartz et al. 2003), with humans used as the reference. First, we identified deletion segments in bat genomes (supplementary fig. S40, Supplementary Material online), lost syntenic gene blocks (supplementary table S39, Supplementary Material online), and gene losses to determine bat-specific DNA losses (supplementary table S40, Supplementary Material online). Second, ChIP-seq data from bats were used to predict the functional domains of enhancers and promoters within bat-specific regions of DNA gain (supplementary table S41, Supplementary Material online). A dual-luciferase reporter assay was employed to validate the enhancer activity of specific DNA gain regions. Finally, the HCEs were parsed using phastCons (in the v1.5 PHAST package) (Siepel et al. 2005) and genomic evolutionary rate profiling software (Cooper et al. 2005; Davydov et al. 2010) based on a multispecies genome alignment of 35 bats. We annotated regulatory elements in noncoding conserved regions, including transcription factor binding sites (supplementary table S42, Supplementary Material online) and long noncoding RNAs (supplementary table S43, Supplementary Material online).
Phylogenetic Relationship Reconstruction
To reconstruct the phylogenetic tree, we generated a whole-genome alignment, orthologous genes, and 4-fold degenerate sites, as well as CNEs and a mitochondrial genome tree (supplementary tables S44 and S45, Supplementary Material online). We employed RAxML-ng v1.1.0 (Kozlov et al. 2019) to reconstruct the ML tree under the best-fit substitution model determined via IQ-TREE v2.1.3 (Minh et al. 2020). To mitigate bias, we examined genes with extreme gene-wise log-likelihood scores via RAxML v8.2.9 (Stamatakis 2014) and assessed substitution saturation using PhyloMAd v1.2 (Duchêne et al. 2018). The AU test (Shimodaira 2002) was conducted in IQ-TREE v2.1.3 (Minh et al. 2020) (supplementary table S46, Supplementary Material online), which was also used to implement CFs (Minh et al. 2020). For species tree reconstruction, we used ASTRAL v5.7.3 (Zhang et al. 2018) and MP-EST v2.0 (Liu et al. 2010).
Gene Introgression and ILS
We depicted gene trees using DensiTree v2.01 (Bouckaert 2010) (supplementary fig. S41, Supplementary Material online) and identified gene introgressions via Dsuite v0.5 (Malinsky et al. 2021) and HyDe v1.0 (Blischak et al. 2018). ILS fragment was identified using CoalHMM v1.0.3 (Dutheil et al. 2009). We used QuIBL (Edelman et al. 2019) to explore which model (ILS or ILS 0+ introgression) is more fit to explain the conflict relationship. To infer the species tree under a gene introgression and ILS scenario, we employed PhyloNet v3.0.1 (Than et al. 2008) and PhyloNetworks v0.12.0 (Solís-Lemus et al. 2017) (supplementary fig. S42, Supplementary Material online). To distinguish the effects of ILS and introgression on controversial relationships, we used quantifying introgression via branch length analysis (Edelman et al. 2019).
Gene Selection Analysis and Cellular Assay
The alignment results of the longest transcripts in all species were clustered using OrthoMCL v1.4 (Li et al. 2003). First, we employed a genome-alignment pipeline based on a gene-alignment method (Zhou et al. 2008) to identify orthologous genes and detect new lineage-specific genes at different interior nodes (supplementary fig. S43, Supplementary Material online). Second, to identify positive selection genes, we used the branch-site model in PAML v4.8 (Yang 2007) and BUSTED v2.5.17 (Murrell et al. 2015). We conducted cell migration experiments to explore the roles of RnaseL residues at positively selected sites in relation to cancer resistance. The cell transwell migration assay is an essential in vitro method commonly used to investigate the migratory capacity and mechanisms of tumor cells under specific conditions, as well as how these mechanisms affect tumor growth, invasion, and metastasis. By regulating the expression or knockout of the specific genes, we can investigate their impact on the migration of tumor cells. Finally, we scrutinized research on 736 cancer genes in the COSMIC database to investigate the evolutionary rate and copy number of cancer genes in bats (supplementary figs. S44 and S45, Supplementary Material online).
Divergence Time Estimation and Demographic History
SNPs for each sample were detected using SAMtools v1.9 (Li and Durbin 2009). We estimated divergence times using r8s (Sanderson 2003) and MCMCTree in PAML v4.8 (Yang 2007) (supplementary figs. S46 and S47, Supplementary Material online) with fossil calibration. The mutation rate derived from r8s was used as an input parameter to reconstruct the demographic history of 19 species using PSMC v0.6.5-r67 (Li and Durbin 2011). We also introduced the data of surface air temperature and global sea level in Fig. 4b (Bintanja et al. 2005). Considering the effects of a wide range of possible tao, u, and g on the estimation of Ne, we have used 4 different values for tao (11, 13, 17, and 19) and 3 different values for generation time (1/2g, g, and 2g) and for mutation rates (1/2u, u, and 2u). The identical Ne patterns were shown, though the values changed for each parameter. With the increasing generation time or mutation rate, the pattern line shifts to the right or left along the “time” x axis, respectively (supplementary fig. S48, Supplementary Material online).
Supplementary Material
Supplementary material is available at Molecular Biology and Evolution online.
Acknowledgments
We thank Pengcheng Wang, Xing Chen, Alice Catherine Hughes, Yanhua Chen, Zhiqian Yu, Xuanjing Li, and Yun Huang for supporting sample collection and preparation. We thank Weiwei Zhai and Alice Catherine Hughes for the comment on the early draft. We thank Xiaocheng Wang for drawing the animal picture. This work was supported by China National GeneBank (CNGB).
Author Contributions
X. Zhou. conceived the study and designed the project. G.L., Q.P., and P.Z. performed data curation, data analysis, and visualization and wrote the manuscript; X.G., Y.Z., J.W., and X.L. performed data curation, data analysis, and visualization; Z.L. prepared the cell cultures and implemented the analysis, data analysis, and visualization; Z.Z. performed data curation and visualization; X. Zhou., G.L., Q.P., P.Z., X.W., Y.Y., X. Zhang., C.H., W.C., and W.L. discussed the results; I.S. and G.F. revised this manuscript. All authors contributed to the data interpretation.
Funding
This project was granted by the National Key Research and Development Projects of the Ministry of Science and Technology of China (2021YFC2301300), the National Natural Science Foundation of China (32070528 and 82050002), the key program of Chinese Academy of Sciences (KJZD-SW-L11), the Beijing Natural Sciences Foundation (JQ19022), and the Institute of Zoology, Chinese Academy of Sciences (2023IOZ0104).
Conflict of Interest
None declared.
Data availability
The genome sequencing data, assemblies, and annotations of species have been deposited in Sciencedb with doi number https://doi.org/10.57760/sciencedb.11285. The data have also been deposited into CNGB Sequence Archive (CNSA) of China National GeneBank DataBase (CNGBdb) with accession number CNP0004892. All custom code has been made available on GitHub at https://github.com/panqi113/bats_project. Any additional information required to reanalyze the data reported in this paper is available upon request.
References
Author notes
Gaoming Liu, Qi Pan, Pingfen Zhu and Xinyu Guo contributed equally to this work.
