DOI: 10.5061/dryad.bzkh18967
title: The legacy of recurrent introgression during the radiation of hares

description: 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 within and among
 species.
description: We generated new genome-wide resequencing data targeting 207,691
 exonic and non-coding regions [totaling 61.7 Megabases (Mb) total] from 14
 hare species (30 individuals)and the outgroup pygmy rabbit
 (<em>Brachylagus idahoensis</em>; 2 individuals). We combined
 these data with published whole exomes from four snowshoes hares
 (<em>Lepus americanus</em>; NCBI Sequence Read Archive
 BioProject PRJNA420081) and extracted data from the reference genome of
 European rabbit (<em>Oryctolagus cuniculus</em>; OryCun2.0) to
 use as a second outgroup. Our total sample of 15 hare species (34
 individuals, 1 to 4 individuals per species) and 2 outgroup species (3
 individuals) included species from all major regions of the
 <em>Lepus</em> native distribution. A smaller subset of 32
 hare individuals (15 hare species), 1 pygmy rabbit individual and the
 European rabbit reference were used for subsequent analyses (check
 Supplementary Table S1 and Materials and Methods). We
 reconstructed the evolutionary history of <em>Lepus </em>using
 maximum likelihood with a concatenated alignment and multispecies
 coalescent approaches, with ASTRAL and SVDquartets. We also studied
 phylogenetic incongruence among exome regions underlying the species tree.
 We infered hybridization among extant and ancestral branches of the
 species tree using network based inferences (PhyloNet and TreeMix) and
 summary statistics (minimum D-statistic, 'f-branch' statistic,
 admixture proportion and fraction of admixture).
description: The pipeline and custom scripts used to recapitulate the analyses
 in this work can be found
 at https://github.com/evochange/hare-phylogenomics. <strong>Supplementary Files:</strong> <strong>Supplementary_Materials_and_Methods.pdf</strong> Supplementary Materials and Methods. <strong>Supplementary_Figures.pdf</strong> Supplementary Figures S1 to S13. <strong>Supplementary_Tables.xlsx</strong> Supplementary Table S1 to S15. Data Files: Pseudo-references <b>BID_CG461_pseudoref.gatk.iteration4.consensus.FINAL.fa.gz</b> <b>LCL_FXG64H.gatk.iteration4.consensus.FINAL.fa.gz</b> <b>LCP_TAN3023.gatk.iteration4.consensus.FINAL.fa.gz</b> <b>LCR_ITA1958.gatk.iteration4.consensus.FINAL.fa.gz</b> <b>LCS_1891.gatk.iteration4.consensus.FINAL.fa.gz</b> <b>LER_ALM1580.gatk.iteration4.consensus.FINAL.fa.gz</b> <b>LFG_FAG114_pseudoref.gatk.iteration4.consensus.FINAL.fa.gz</b> <b>LGR_CRE2553.gatk.iteration4.consensus.FINAL.fa.gz</b> <b>LHB_HAB68_pseudoref.gatk.iteration4.consensus.FINAL.fa.gz</b> <b>LMS_PRI2460.gatk.iteration4.consensus.FINAL.fa.gz</b> <b>LOT_2109.gatk.iteration4.consensus.FINAL.fa.gz</b> <b>LST_STA65_pseudoref.gatk.iteration4.consensus.FINAL.fa.gz</b> <b>LTM_FIN2191.gatk.iteration4.consensus.FINAL.fa.gz</b> <b>LTW_JRRK3.gatk.iteration4.consensus.FINAL.fa.gz</b> Pseudo-reference genomes for each species generated by iteratively mapping (four iterations and allowing ambiguity codes) cleaned single and paired-end reads from one individual per species (Supplementary Table S1 available on Dryad) to the European rabbit reference genome (OryCun2.0). Names correspond to sample IDS in Supplementary Table S1. Consensus sequences: <b>consensus_fasta_sequences.zip</b> Consensus fasta sequences for hares and pygmy rabbit used as input for all analyses in this study. Inside the zip file there is a file per individual (identified by Sample ID in Supplementary Table S1) containing one fasta sequence per chromosome.  ASTRAL: <b>astral5.6.3_chrX_TCs5_gene_trees.tre</b> ASTRAL species tree obtained from 181 chromosome X gene trees, displayed in Supplementary Figure S8. <b>astral5.6.3_lineage_TCs5_gene_trees.tre</b> <b>astral5.6.3_species_TCs5_gene_trees.tre</b> ASTRAL species tree obtained with 8889 gene trees without (“lineage”) or with (“species”) assigning individuals to species. These trees are displayed in Supplementary Figure S2. <b>chr_X_bestTrees_TCs5.txt.gz</b> Maximum likelihood gene trees from 181 X-linked regions used as input to generate the ASTRAL species tree astral5.6.3_chrX_TCs5_gene_trees.tre displayed in Supplementary Figure S8. <b>gene_trees_TCs5.tree.gz</b> Maximum likelihood gene trees from 8889 regions used as input to generate the ASTRAL species trees astral5.6.3_lineage_TCs5_gene_trees.tre and astral5.6.3_species_TCs5_gene_trees.tre displayed in Supplementary Figure S2. <b>partitions_fasta_alignments.zip</b> <b>gene_trees_raxml_output.zip</b> Fasta alignments generated from partitioning the genome in 50 kilobase (kb) windows with no overlap. Only windows inside target regions in the capture and 200 bp flanking regions were retained. Positions with more than 30% missing data were excluded from the alignments. These fasta alignments were used as input for RAxML to generate maximum likelihood gene trees (gene_trees_raxml_output.zip) that were later used for the multispecies coalescent analysis implemented in ASTRAL SVDquartets: <b>snp_alignment_33indvs_chr_X.nex</b> <b>snp_alignment_33indvs_WithX.nex</b> SNP alignments in nexus format containing only 1473 X-linked independent SNPs (snp_alignment_33indvs_chr_X.nex) or containing 45,779 exome-wide independent SNPs (snp_alignment_33indvs_WithX.nex) used as input for SVDquartets.  <b>lineage_33indvs_chr_X_qall_b1000</b> <b>lineage_33indvs_chr_X_qall_b1000.tre</b> <b>lineage_33indvs_WithX_qall_b1000</b> <b>lineage_33indvs_WithX_qall_b1000.tre</b> <b>species_33indvs_chr_X_qall_b1000</b> <b>species_33indvs_chr_X_qall_b1000.tre</b> <b>species_33indvs_WithX_qall_b1000</b> <b>species_33indvs_WithX_qall_b1000.tre</b> SVDquartets output log files (no extension) and tree files (.tre). The analyses were performed with (“species”) or without (“lineage”) assigning individuals to species, and with X-linked (“chr_X”) or whole-exome (“WithX”) SNPs. Exome-wide trees are displayed in Supplementary Figure S3 and Chromosome X trees are displayed in Supplementary Figure S8. Concatenated whole-exome maximum likelihood tree: <b>exome_alignment_33indvs_f100.fas.gz</b> Fasta alignment containing 11,949,529 exome-wide concatenated positions with no missing data used to generate the Maximum Likelihood tree in Supplementary Figure S1.  <b>RAxML_bipartitions.exome_alignment_33indvs_f100.tree</b> <b>RAxML_bootstrap.exome_alignment_33indvs_f100.tree</b> <b>RAxML_info.exome_alignment_33indvs_f100.tree</b> RAxML output files used to estimate a single bifurcating phylogeny using the fasta alignment exome_alignment_33indvs_f100.fas.gz and using a maximum likelihood (ML) search and rapid bootstrapping run under the GTR+GAMMA model of<b> </b>sequence evolution. Bayesian estimation of divergence times: <b>9k_CDS_123_alignments.phy</b> PHYLIP alignment with data for 15 hare species (1 individual per species) and 1 outgroup species containing 3 partitions corresponding to the 1<sup>st</sup>, 2<sup>nd</sup> and 3<sup>rd</sup> codon positions of 9015 protein coding genes used as input for Bayesian estimation of divergence times of the <i>Lepus </i>species tree with MCMCtree. <b>FigTree_mcmc_9k_123_CDS_tree3_run1.tre</b> <b>FigTree_mcmc_9k_123_CDS_tree3_run2.tre</b> <b>FigTree_mcmc_9k_123_CDS_tree5_run1.tre</b> <b>FigTree_mcmc_9k_123_CDS_tree5_run2.tre</b> Output tree files obtained after Bayesian estimation of divergence times with MCMCtree. We run MCMCtree twice (run 1 and run2) per calibration strategy (a) molecular-based dates extrapolated from deep fossil calibrations in the lagomorphs (tree 3) or (b) fossil calibrations for the divergence of the genus Lepus and the Lepus species (tree5). See Materials and Methods for details. Run1 trees are displayed in Supplementary Figure S5. <b>mcmc_9k_123_CDS_tree3.zip</b> <b>mcmc_9k_123_CDS_tree5.zip</b> Folders containing all input files (alignment and species with age priors), control files (mcmc.ctl) and output files (trees and log files for hessian and approximate maximum likelihood estimations) obtained from two MCMCtree runs with different calibration strategies (a) molecular-based dates extrapolated from deep fossil calibrations in the lagomorphs (tree 3) and (b) fossil calibrations for the divergence of the genus Lepus and the Lepus species (tree5). <b>Majority Rule Consensus Tree</b> <b>RAxML_MajorityRuleExtendedConsensusTree_IC.gene_trees_TCs5.tree</b> Majority Rule Consensus Tree obtained with RAxML<i> </i>(-L MRE option) with the set of gene trees (gene_trees_TCs5.tree.gz) used to generate the ASTRAL species tree. This tree is displayed in Supplementary Figure S7. PhyloNet: <b>gene_trees_TCs5_m0_plopt.nex.gz</b> <b>gene_trees_TCs5_m1_plopt.nex.gz</b> <b>gene_trees_TCs5_m2_plopt.nex.gz</b> <b>gene_trees_TCs5_m3_plopt.nex.gz</b> <b>gene_trees_TCs5_m4_plopt.nex.gz</b> Nexus file containing 8889 gene trees and PhyloNet command used to produce pseudo-maximum likelihood networks with 0 to 4 reticulation events displayed in Supplementary Figure S10.  <b>gene_trees_TCs5_m0_plopt.out.gz</b> <b>gene_trees_TCs5_m1_plopt.out.gz</b> <b>gene_trees_TCs5_m2_plopt.out.gz</b> <b>gene_trees_TCs5_m3_plopt.out.gz</b> <b>gene_trees_TCs5_m4_plopt.out.gz</b> PhyloNet raw outputs obtained from pseudo-maximum likelihood network analyses with 0 to 4 reticulation events. Outputs contain the results of 10 interactions and the five best networks. The best networks from each run are displayed in Supplementary Figure S10. TreeMix: <b>snp_alignment_LepusOnly_hap_treemix_inputtreemix_file_out.txt.gz</b> Allele frequencies for 30,709 biallelic SNPs from 15 hare species used as input for ancestral population graph reconstruction with TreeMix.  <b>treemix_outputs_m0-9_reticulations.zip</b> Output files for 10 TreeMix runs with 0 to 9 reticulations events using allele frequency data of 30,709 biallelic SNPs for 15 hare species snp_alignment_LepusOnly_hap_treemix_inputtreemix_file_out.txt.gz. The networks in these files and corresponding residuals are displayed in Supplementary Figures S11 and S12. Genetic divergence and nucleotide diversity: <b>wg_filtered_wHead.geno.gz</b> Geno file containing the whole-exome alignment used to calculate genetic divergence between and nucleotide diversity and heterozygosity within hare species with <i>Genomic Generals </i>(https://github.com/simonhmartin/genomics_general) scripts. This was also the input file used to calculate f<sub>hom</sub> among black-tailed jackrabbits and snowshoe hare populations. <b>wg_filtered_wHead.pi</b> <b>wg_filtered_wHead.dxy</b> Nucleotide diversity (“pi”) and genetic divergence (“d<sub>xy</sub>”) in 1 Megabase overlapping windows calculated for 15 hare species with the script <i>popgenWindows.py </i>from <i>Genomics General </i>(https://github.com/simonhmartin/genomics_general). Results average across windows can be found in Supplementary Table S8 and S9.  Minimum D-statistics (<em>D<sub>min</sub></em>): <b>by_species_wg_filtered_wgHead.freq.tsv.gz</b> Derived allele frequencies for 15 hare species calculated with <i>freq.py</i> from <i>Genomics General </i>(https://github.com/simonhmartin/genomics_general) and the input file wg_filtered_wHead.geno.gz. <b>D_stat_by_species_results_FINAL.txt</b> <b>Dmin_by_species_sig_FINAL.txt</b> D-statistic values (D_stat_by_species_results_FINAL.txt)<b> </b>for all possible trios<b> </b>of species in the dataset of 15 hare species analyzed in this study and the minimum D-statistic per trio (Dmin_by_species_sig_FINAL.txt). These results are represented in Supplementary Figure S9. ‘f-branch’ statistics (<em>f<sub>b</sub>(C)</em>): <b>all_frequencies_filtered.freq.tsv.zip</b> <b>fb_C_comparisons_f_G.txt</b> Derived allele frequencies for several taxa configurations (per species, per individual, per population or per haploid chromosome of the white-sided jackrabbit) used to calculate admixture proportion for all f(A,B,C,O) conformations necessary for the ‘f-branch’ statistic. In the table fb_C_comparisons_f_G.txt we list the ABBA-BABA configuration used to calculate admixture proportion for each branch<b> </b><i>b</i> and <i>C</i> pair in the species tree (Figure 1b). The order of the columns is 1) branch b<i>; </i>2) branch a (sister to b); 3) all B descendants of b; 4) all A descendants of a; 5) donor branch C; 6) taxa serving as Ca (or P3a); 7) taxa serving as Cb (or P3b). <b>fG_results_species.txt</b> Admixture proportion results for all f(A,B,C,O) configurations in fb_C_comparisons_f_G.txt used as input to calculate ‘f-branch’ statistics (<em>f<sub>b</sub>(C)</em>) among all pairs of <i>b </i>and <i>C </i>branches<b> </b>in the species tree (Figure 1b). Scans of fraction of admixture (<em>f<sub>d</sub></em>): <b>L_othus_vs_L_americanus_fdscan_50k5k.out</b> <b>L_othus_vs_L_americanus_fdscan_99.5_sig_annotation.txt</b> <b>L_townsendii_vs_Lamericanus_fd_scan_99.5_sig_annotation.txt</b> <b>L_townsendii_vs_Lamericanus_fdscan_50k5k.out</b> <b>L_timidus_vs_Lamericanus_fdscan_50k5k.out</b> <b>L_timidus_vs_Lamericanus_fdscan_99.5_sig_annotation.txt</b> Fraction of admixture (<i>f<sub>d</sub></i>) between snowshoe hares and Alaskan hares, mountain hares or white-tailed jackrabbits for 50 kb overlapping (5 kb) sliding windows (.out files) and European rabbit annotations for the significant outlier windows in each scan (.txt files). Fraction admixture files were obtained with <i>ABBABABAwindows.py </i>from <i>Genomics General</i> (https://github.com/simonhmartin/genomics_general). <b>collective_sig_fd.txt</b> European rabbit annotations of the windows that are significant in all three fraction of admixture (<em>f<sub>d</sub></em>) scans performed between snowshoe hares and Alaskan hares, mountain hares or white-tailed jackrabbits. This annotation was used for Gene Ontology enrichment analysis with g:Profiler. <b>dxy_windows50k_sscc.out</b> Genetic divergence (d<sub>xy</sub>) among snowshoe hares, Alaskan hares, mountain hares and white-tailed jackrabbits in 50 kb overlapping (5 kb) windows obtained with <i>popgenWindows.py </i>from <i>Genomics General </i>(https://github.com/simonhmartin/genomics_general), and used to determine d<sub>xy</sub> of the fraction of admixture (<em>f<sub>d</sub></em>) outlier windows. Results of this scan are displayed in Supplementary Figure S13.