Homology inference
BUSCO analysis results indicated that the bat genome assemblies contained between 68.5 and 96.5% of the single–copy orthologs present among mammals (Figure 1). Orthologs were grouped into 42,441 groups, of which 1,193 were single copy. In total, 5,528 orthogroups had at least one representative in each of the entire set of 43 species that were analyzed. In contrast, 1,055 orthogroups were represented in at least 50% of bat species but missing from the six outgroup taxa (Supplementary table 3). To annotate diets, we used the semi–quantitative database compiled by Rojas, Ramos, Fonseca and Dávalos (2018), which focuses on neotropical noctilionoids (Yangochiroptera), supplemented with summaries from Animal Diversity Web (https://animaldiversity.org/).

Enrichment in chiropteran gene losses

We inferred the first densely sampled chiropteran phylogeny based on hundreds of loci (Figure 1). Our results confirmed the monophyly of the suborders Yinpterochiroptera and Yangochiroptera but the phylogeny of the neotropical leaf-nosed bats (family Phyllostomidae) differed from previous phylogenies (Davalos, Velazco, & Rojas, 2020), in the paraphyly of plant-eating lineages. As the obtained phylogeny is the best supported by all genome-scale analyses available thus far (S.J. Rossiter and M. Hiller pers. obs.), we used this phylogeny for gene family evolution analyses.
A total of 1,115 genes (Supplementary Table 4) were identified as missing in bats, even after filtering BLAST searches against the genomes and transcriptomes. Based on this list, we identified eight over-represented pathways in BioPlanet (Supplementary Table 5) and 63 GO terms in GOnet (Supplementary Table 6). While the former included 104 genes, of which 49 were unique, the latter included 339 unique missing genes. As expected, over-represented categories included chemosensory gene losses in the categories of olfactory transduction, G-protein–coupled receptors (GPCR), and signal transduction. BioPlanet pathways were also enriched for less common categories including immune system pathways that include alpha and beta defensins, antigen process and presentation, and graft-versus-host disease (Supplementary Table 5). GOnet analyses also identified the expected enrichments in chemosensory gene losses and general response to stimuli categories, but also included many more immune categories. Of the latter, the categories comprising the most genes were defense response (58 genes), defense response to other organism (54), response to bacterium (53), innate immune response (46), defense response to bacterium (44), humoral immune response (34), adaptive immune response based on somatic recombination of immune receptors built from immunoglobulin superfamily domains (23), lymphocyte mediated immunity (23), and leukocyte mediated immunity (23). Although these categories share many genes across them, a preponderance of immune system losses is evident in Supplementary Table 6. We used BioRender to summarize the immune gene ontology categories and connections, highlighted in Figure 2.

Gene family evolution

To determine branches and gene families with significant gene family expansions and contractions, we analyzed 14,171 orthogroups under two models: a global rate of gene family evolution, and a three multi–λ model. The three–rate model best fit the data (p < 0.01), this analysis estimated a higher rate of gene family turnover (λYangochiroptera = 0.0048) in the ancestral Yangochiroptera lineage than in the Yinpterochiroptera ancestral lineage (λYinpterochiroptera = 0.0024), with the lowest turnover rate for outgroup lineages (λOutgroups = 0.0017).
With an estimated error distribution of 0.049 (i.e., 4.9% of gene families showed an error in gene size), we identified 2,555 orthogroups with significant expansions or contractions along at least one of the branches in the species tree. Given our focus on immune system and metabolic evolution, we extracted PANTHER annotations for the most frequent (900 orthogroups) biological process categories: immune response, metabolic process, and cellular process. All GOnet annotations were used and binned into immune, metabolic, and two additional processes: response to stress (271 orthogroups) and autophagy (19). PANTHER and GOnet annotations were mostly complementary; orthogroups were often annotated in one database but not the other (1,268 orthogroups). When annotations were available from both databases, these tended to agree on both immune and metabolic categories (594 orthogroups), or to agree on one or the other (404), with only 48 orthogroups disagreeing completely in immune and metabolic annotations between the databases. The remaining 241 were not annotated in either database. Categories, locations, and size of significant gene family changes were summarized using tools in the R package ggtree (Yu, Smith, Zhu, Guan, & Lam, 2017) and are shown in Figure 3. Although several pairs of sister species showed apparently large differences along corresponding tips (e.g., Rhinolophus , Miniopterus ), such variation is common in analyses that include genome assemblies of varying quality (Denton et al., 2014; Tsagkogeorga et al., 2017). Therefore, we focus our discussion on the more robust inference of gene family expansions and contractions for non-sister lineages. Taking advantage of relatively dense taxon sampling within bats (Figure 1), we focused on parallel adaptation to plant–rich diets across suborders Yinpterochiroptera and Yangochiroptera.

Selection tests

Branch–site selection tests identified 37 of 268 single–copy genes with evidence for positive selection, of which 27 remained after false discovery rate correction (Table 1). This subset included genes involved in interferon-gamma (IFNG) signaling, inflammatory response, as well as cytokines, chemokines, and interleukins. A total of 16, 979 branches across 268 genes were analysed using the aBSREL model in HyPhy. After FDR correction, 683 branches from 191 gene trees were found to be significant, 25 of which were consistent with CODEML results (Supplementary table 8).