Inferring phylogeny and introgression using RADseq data: an example from flowering plants (Pedicularis: Orobanchaceae)

Abstract

Phylogenetic relationships among recently diverged species are often difficult to resolve due to insufficient phylogenetic signal in available markers and/or conflict among gene trees. Here we explore the use of reduced-representation genome sequencing, specifically in the form of restriction-site associated DNA (RAD), for phylogenetic inference and the detection of ancestral hybridization in non-model organisms. As a case study, we investigate Pedicularis section Cyathophora, a systematically recalcitrant clade of flowering plants in the broomrape family (Orobanchaceae). Two methods of phylogenetic inference, maximum likelihood and Bayesian concordance, were applied to data sets that included as many as 40,000 RAD loci. Both methods yielded similar topologies that included two major clades: a "rex-thamnophila" clade, composed of two species and several subspecies with relatively low floral diversity, and geographically widespread distributions at lower elevations, and a "superba" clade, composed of three species characterized by relatively high floral diversity and isolated geographic distributions at higher elevations. Levels of molecular divergence between subspecies in the rex-thamnophila clade are similar to those between species in the superba clade. Using Patterson's D-statistic test, including a novel extension of the method that enables finer-grained resolution of introgression among multiple candidate taxa by removing the effect of their shared ancestry, we detect significant introgression among nearly all taxa in the rex-thamnophila clade, but not between clades or among taxa within the superba clade. These results suggest an important role for geographic isolation in the emergence of species barriers, by facilitating local adaptation and differentiation in the absence of homogenizing gene flow.

Figure 1

Sampling localities and photographs of taxa in Pedicularis sect. Cyathophora. Codes for taxon names are listed in Table 1.

Figure 2

Coalescent model used to simulate RADseq-like data. Unidirectional migration (gene flow) from taxon C_i into B_i (solid gray arrow) occurred over 10,000 generations, beginning 50,000 generations in the past. Three data sets were created with different rates of gene flow, and one additional data set with gene flow from both C and C_i into B_i (dashed gray arrow).

FIGURE 3

A schematic description of Patterson’s four-taxon D-statistic test for hybridization and the five-taxon partitioned D-statistic test. a) The four-taxon test measures asymmetry in the occurrence of two incongruent allele patterns (shown in the box), which are expected to arise with equal frequencies in the absence of introgression, but to deviate from symmetry if P3 exchanged genes with either P1 or P2 to the exclusion of the other. b) The five-taxon test samples an additional population within the P3 lineage, the two samples now denoted P3₁ and P3₂, and measures asymmetry in three sets of allele patterns (the three boxes shown), thus measuring three separate D-statistics. The first (D₁) measures asymmetry in incongruent alleles where the derived allele is present in P3₁ but not P3₂, the second (D₂) where it is present in P3₂ but not P3₁, and the third where the derived allele is shared by both P3₁ and P3₂, having arisen in their ancestor (branch indicated by **).

FIGURE 4

Phylogeny of Pedicularis sect. Cyathophora inferred from RADseq data. ML trees were estimated from the sparse (minimum-taxa) supermatrix a) and the densely sampled (full-taxa) supermatrix b), yielding high bootstrap support (100 except where indicated). Primary concordance trees c, d) and population trees e, f) were inferred on 945 variable loci from the full-taxa supermatrix, at either α = 0.01 (c, e) or 100 (d, f). The 95% CI for CFs are shown on primary concordance trees, those in bold did not overlap with any conflicting CF. Quartet CFs are shown on population trees.

Figure 5

Censored data sets, in which one or more taxa have been removed to minimize the effect of their introgressed alleles, are used to compare alternative topologies within the Shimodeira–Hasegawa test. The ML topology with P. rex subsp. rex and P. rex subsp. lipskyana removed a) is compared with an alternative topology in which P. thamnophila is monophyletic b) Similarly, the ML topology with both subspecies of P. thamnophila removed c) is compared with two alternative resolutions for the three subspecies of P. rex d) and e).

Figure 6

Temporal and geographic reconstructions of population divergence and gene flow, as inferred from the partitioned D-statistic. Codes for taxon names are listed in Table 1. a) Phylogeny of rex-thamnophila based on censored subtree comparisons which account for hybrid introgression. Gray bars indicate distinct time periods discussed in the text, and arrows are the direction of introgressive gene flow. b) Schematic reconstruction of historical biogeography in the rex-thamnophila clade, inferred from present-day distributions and patterns of historical introgression. Overlapping edges indicate geographic overlap, and arrows represent introgression.