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).