Recombination-Aware Phylogenomics Reveals the Structured Genomic Landscape of Hybridizing Cat Species

Abstract

Current phylogenomic approaches implicitly assume that the predominant phylogenetic signal within a genome reflects the true evolutionary history of organisms, without assessing the confounding effects of postspeciation gene flow that can produce a mosaic of phylogenetic signals that interact with recombinational variation. Here, we tested the validity of this assumption with a phylogenomic analysis of 27 species of the cat family, assessing local effects of recombination rate on species tree inference and divergence time estimation across their genomes. We found that the prevailing phylogenetic signal within the autosomes is not always representative of the most probable speciation history, due to ancient hybridization throughout felid evolution. Instead, phylogenetic signal was concentrated within regions of low recombination, and notably enriched within large X chromosome recombination cold spots that exhibited recurrent patterns of strong genetic differentiation and selective sweeps across mammalian orders. By contrast, regions of high recombination were enriched for signatures of ancient gene flow, and these sequences inflated crown-lineage divergence times by ∼40%. We conclude that existing phylogenomic approaches to infer the Tree of Life may be highly misleading without considering the genomic architecture of phylogenetic signal relative to recombination rate and its interplay with historical hybridization.

Keywords: Felidae; X chromosome; hybridization; phylogenomics; recombination.

Fig. 1.

Previous hypotheses of felid phylogenetic relationships based on 1) maternally inherited mitogenome sequences (Li, Davis, et al. 2016), 2) biparentally inherited SNP markers (Li, Davis, et al. 2016), and 3) paternally inherited Y chromosome sequences (Johnson et al. 2006; Luo et al. 2014).

Fig. 2.

(A) Left: the top three most frequently recovered phylogenies for six felid lineages, inferred from window-based ML supermatrix analysis of the 1.5-Gb whole-genome alignment. Right: the frequency of each topology within autosomes and chrX, based on partitioning into low (<0.5 cM/Mb) and high (>2 cM/Mb) recombination rates (Li, Hillier, et al. 2016). Red arrows indicate three lineages in which the most frequent topology in the LRchrX was not the same as the most frequent genome-wide topology. (B) Distribution of different topologies on chrX, and their relationship to recombination rate (bottom).

Fig. 3.

Left: frequency distribution (y axis) of genomic window-based estimates for divergence time (t2) (x axis) for the three felid lineages in which the most common autosomal and chrX phylogenies differed in figure 2. Right: comparison between divergence time estimates for the sister-species pair found in the two most frequent topologies. In the Asian leopard cat (A) and lynx (B) lineages, the topology enriched in the low-recombination cold spot of chrX (gray and orange trees, respectively) had a significantly older divergence time (Mann–Whitney U test, ***P < 0.001, **P < 0.01, *P < 0.05, ns = not significant, na = not analyzed due to too few windows) than topology 1 (blue tree) that was most commonly observed across the genome, and enriched in regions of higher recombination (see fig. 4). Within the Panthera genus (C), trees supporting jaguar + lion are enriched in low-recombining regions of the genome (see fig. 4). In the low recombination regions of chrX, divergence time estimates for the lion + jaguar node from LRchrX are much lower than elsewhere, consistent with published evidence for selective sweeps across this region (Figueiró et al. 2017).

Fig. 4.

Local recombination rate correlates with topology within three felid lineages. Individual plots show the frequency of a particular topology (y axis) within windows of increasing recombination rate (x axis). The tree inferred here to reflect speciation history in each lineage (topology 3 for the Leopard cat clade, topology 2 for the Lynx and Panthera clades) is enriched in regions of low recombination on the autosomes and chrX, whereas the frequency of the most prevalent tree (topology 1 in all three lineages), inferred here to derive from postspeciation admixture, increases with recombination rate.

Fig. 5.

The influence of local chrX recombination rate on topology and divergence times. (A) ML trees are shown for the autosomes (concatenation of sites selected every 10 kb across the chromosome alignments) and from all sites on chrX within windows of low (<0.5 cM/Mb) and high (>2.5 cM/Mb) recombination rates. Clade coloring follows figure 1. (B) Top: regional variation in phylogenetic signal density across chrX. The distribution of windows supporting the four most common topologies (color coded according to the key underneath) is shown in relation to recombination rate (in cM/Mb, y axis), plotted underneath the chromosome ideogram. Bottom: timetrees inferred from the four most common trees on chrX, showing the impact of recombination rate on topology and divergence time estimates. Node ages shown for Trees 1–4 are the average of point estimates derived from all 100-kb windows that supported a particular topology: Tree 1 (n = 457), Tree 2 (n = 100), Tree 3 (n = 70), and Tree 4 (n = 64). The timescale below each of the four timetrees indicates millions of years before present (MYBP). Shown above each timetree are two parameters that are distorted by hybridization and recombination, the ratio of external/internal branch length (ebl/ibl), or tree compression, and the shape parameter (α) of the gamma distribution of rate variation among sites. Red arrows indicate the crown node of each felid clade with more than one representative and highlight the differences in the estimated clade ages between regions of low- and high-recombination rate.

Fig. 6.

(A) Remarkable conservation of gene order and recombination rate across eutherian mammal orders. Comparison of colinear recombination maps of chrX in cat (Li, Hillier, et al. 2016) and human (Nagaraja et al. 1997). The human recombination rate image is reproduced (with permission) from the original manuscript figure. In each image, the recombination map (in cM) is plotted along the y axis relative to the physical sequence coordinates on the x axis (in Mb). The red-bounded boxes highlight conserved, syntenic regions of elevated recombination rate, whereas yellow boxes highlight conserved syntenic regions of reduced recombination rate. Lettering in the human chrX image refers to shifts in recombination rate noted in the article. (B) Corresponding boundaries of the conserved chrX low and high-recombination windows in human and pig (data from Fernández et al. [2014]) mapped relative to the domestic cat X chromosome sequence.

Fig. 7.

Heatmap showing statistical support for introgression between pairs of species, based on Patterson’s D-statistic. Each cell in the heatmap indicates the average pairwise Z-score inferred from all possible combinations of species trios including that pair of taxa (see supplementary table S4, Supplementary Material online), assuming the species tree inferred from chrX that is shown above and to the left of the matrix. The tree on the left includes dashed red lines indicating postspeciation admixture identified in the heatmap. Intralineage Z-scores are displayed inside black boxes along the diagonal. Z-scores for the autosomes (upper, red shades) and the X chromosome (lower, blue shades) are shown separately. Using a significance threshold (α = 0.05) that was Bonferroni-corrected for 20 comparisons (the maximum number of trios involving a given pair), Z-scores > 2.81 are considered significant. We conservatively display only chrX Z-scores >5, and autosomal Z-scores >10.

Fig. 8.

Quantification and chromosomal distribution of introgressed blocks in the genomes of two felid species (Asian golden cat and marbled cat), resulting from hybridization of their common ancestor with the progenitors of three different lineages of the Felidae. (A) Red and blue arrows show the estimated number of windows which support gene flow between two lineages. (B) Chromosome distribution of windows with signatures of gene flow for different species combinations. (C) Spatial distribution of introgressed blocks derived from the hybridization events shown are shown for chromosome C1 (see supplementary figs. S2–S4, Supplementary Material online, for full results).