The Legacy of Recurrent Introgression during the Radiation of Hares

Abstract

Hybridization may often be an important source of adaptive variation, but the extent and long-term impacts of introgression have seldom been evaluated in the phylogenetic context of a radiation. Hares (Lepus) represent a widespread mammalian radiation of 32 extant species characterized by striking ecological adaptations and recurrent admixture. To understand the relevance of introgressive hybridization during the diversification of Lepus, we analyzed whole exome sequences (61.7 Mb) from 15 species of hares (1-4 individuals per species), spanning the global distribution of the genus, and two outgroups. We used a coalescent framework to infer species relationships and divergence times, despite extensive genealogical discordance. We found high levels of allele sharing among species and show that this reflects extensive incomplete lineage sorting and temporally layered hybridization. Our results revealed recurrent introgression at all stages along the Lepus radiation, including recent gene flow between extant species since the last glacial maximum but also pervasive ancient introgression occurring since near the origin of the hare lineages. We show that ancient hybridization between northern hemisphere species has resulted in shared variation of potential adaptive relevance to highly seasonal environments, including genes involved in circadian rhythm regulation, pigmentation, and thermoregulation. Our results illustrate how the genetic legacy of ancestral hybridization may persist across a radiation, leaving a long-lasting signature of shared genetic variation that may contribute to adaptation. [Adaptation; ancient introgression; hybridization; Lepus; phylogenomics.].

Figure 1

Hare (Lepus spp.) evolutionary history and biogeography. a) The distribution of the 15 hare species studied in this work obtained from the IUCN database (http://www.iucnredlist.org); b) Coalescent species tree estimated with ASTRAL assigning individuals to species groups ( 1 except for the white-sided jackrabbit); c) Ancestral range reconstruction under the DIVALIKE+J model implemented in BioGeoBEARS on the divergence time tree estimated with MCMCtree and calibrated with deep fossil calibration points. Pie charts represent the probability of each potential range, and squares represent the current range of extant taxa. Node labels represent estimated divergence times in millions of years. Confidence intervals, an alternative biogeographic reconstruction based on Lepus fossil calibrations and ancestral ranges for the shoulders of the tree are provided in Supplementary Material available on Dryad.

Figure 2

The hare (Lepus spp.) species tree is underlined by extensive gene tree incongruence. Split network constructed from 8,889 gene trees (5% threshold) shows discordance among the gene tree topologies (cuboid structures represent alternative topologies) in deeper nodes of the species tree. Species are marked in accordance with the continents where they are distributed.

Figure 3

Admixture events are distributed across the hare (Lepus spp.) species tree. The “f-branch” statistic represents excess allele sharing between branches (y-axis) and (x-axis) of the species tree in Figure 1b. The gradient represents the score, dark gray represents tests not consistent with the species tree (for each branch , having itself or a sister taxon as donor ) and asterisks denote block jackknifing significance at 0.05 (after Bonferroni correction). Tips of the tree are marked according to their current distribution, and ancestral tips (dashed lines and labeled with letters) are marked according to the ancestral range reconstruction in Figure 1c [two colors represent the inference of two ranges with equivalent (0.5) probability]. Ancestral tips are labeled from A to L corresponding to labels in Supplementary Table S11 available on Dryad.

Figure 4

Admixture proportions decrease with genetic divergence between allopatric species. We plot the “f-branch” ( values against exome-wide divergence (d for extant species pairs, differentiating species with overlapping (sympatric/parapatric) or nonoverlapping (allopatric) distributions. The tendency line represents a linear regression relating d and all (significant and non-significant) values calculated with function lm() in R.

Figure 5

The impact of ancestral introgression on extant northern latitude species. a) Events of past and recent admixture inferred in this study involving snowshoe, Alaskan hares, mountain hares and white-tailed jackrabbits. The arrow indicates the direction of ancestral introgression inferred by PhyloNet and “f-branch” and dashed lines indicate recent introgression inferred with “f-branch.” Values above the lines represent admixture proportion (f) estimated with “f-branch” or inheritance probability (IP) estimated with PhyloNet; b) Genomic distribution of ancestral blocks of introgression (circles) inferred as shared outlier windows among fraction of admixture ( analysis testing for admixture among snowshoe hares and each one of the three other northern latitude species. Tests were run in 50 kb genomic sliding windows, and outlier windows are in the top 0.5% of the distribution.