DOI: 10.5061/dryad.rn8pk0ph9
title: Ancient hybridization leads to the repeated evolution of red flowers across a monkeyflower radiation

description: The re-use of old genetic variation can promote rapid diversification in
 evolutionary radiations, but in most cases, the historical events
 underlying this divergence are not known. For example, ancient
 hybridization can generate new combinations of alleles that sort into
 descendant lineages, potentially providing the raw material to initiate
 divergence. In the Mimulus aurantiacus species complex, there is evidence
 for widespread gene flow among members of this radiation. In addition,
 allelic variation in the MaMyb2 gene is responsible for differences in
 flower color between the closely related ecotypes of subspecies puniceus,
 contributing to reproductive isolation by pollinators. Previous work
 suggested that MaMyb2 was introgressed into the red-flowered ecotype of
 puniceus. However, additional taxa within the radiation have independently
 evolved red flowers from their yellow-flowered ancestors, raising the
 possibility that this introgression had a more ancient origin. In this
 study, we used repeated tests of admixture from whole-genome sequence data
 across this diverse radiation to demonstrate that there has been both
 ancient and recurrent hybridization in this group. However, most of the
 signal of this ancient introgression has been removed due to selection,
 suggesting that widespread barriers to gene flow are in place between
 taxa. Yet, a roughly 30 kb region that contains the MaMyb2 gene is
 currently shared only among the red-flowered taxa. Patterns of admixture,
 sequence divergence, and extended haplotype homozygosity across this
 region confirm a history of ancient hybridization, where functional
 variants have been preserved due to positive selection in red-flowered
 taxa but lost in their yellow-flowered counterparts. The results of this
 study reveal that selection against gene flow can reduce genomic
 signatures of ancient hybridization, but that historical introgression can
 provide essential genetic variation that facilitates the repeated
 evolution of phenotypic traits between lineages.
description: Genome resequencing and variant
 calling  Leaf tissue was
 collected from the wild (Table S1), and DNA from the 10 new samples was
 extracted, prepared, and sequenced as described in Stankowski et al.
 (2019). Variant calling, filtering, and haplotype phasing were performed
 as described in Stankowski et al. (2019) after mapping reads from all 47
 samples to the <em>M. aurantiacus</em> reference genome. We
 then further filtered the VCF file for biallelic SNPs and removed sites
 with missing genotype data. The final data set contained 12,730,568 SNPs
 across all 47 samples.   PCA, network, and phylogenetic
 analyses Due to evidence
 of previous gene flow among taxa in this radiation (Stankowski et al
 2019), we generated a neighbor-joining splits network
 using <em>SplitsTree4
 v</em>.<em>4</em>.<em>17</em>.1<em>
 </em>(Huson, Kloepper, & Bryant 2008) from all 47 samples. A
 splits network provides a way to visualize historical events, such as
 hybridization, in a tree-like framework without requiring a fully
 bifurcating tree. We filtered SNPs for linkage disequilibrium in
 <em>Plink </em>version 1.90 (Chang et al. 2015) to remove
 variant sites with an <em>r<sup>2</sup></em>
 greater than 0.2 in 50 kb windows, sliding by 10 kb, resulting in a data
 set containing 1,640,850 SNPs. For comparison, we also used IQ-TREE
 v1.6.12 (Nguyen et al. 2015) to generate a maximum likelihood consensus
 tree from a concatenated dataset of all 12,730,568 biallelic SNPs across
 the 47 samples. Support for each node was generated from 1000 bootstrap
 replicates, visualized using FigTree v1.4.4 (Rambaut, 2009), and rooted on
 the branch leading to the <em>M. clevelandii</em> samples. To
 assess clustering patterns among more closely related samples, we ran a
 principal components analysis (PCA) in <em>Plink</em> using
 only the 27 samples from clade D<em>. </em> Tests for genome-wide admixture 
 We calculated Patterson’s <em>D
 </em>for all possible groups of three in-group taxa using the Dtrios
 command in the <em>Dsuite</em> software package (Malinsky,
 Matschiner, & Svardal 2021) based on the relationships inferred
 from genome-wide data in Stankowski et al. (2019), with
 <em>M</em>. <em>clevelandii </em>used as the
 outgroup for all tests. The four samples from subspecies
 <em>grandiflorus</em> were removed prior to this analysis, as
 the calculation of <em>D</em> using these samples does not
 provide information about the history of introgression between clade D and
 <em>aridus. </em>We corrected for multiple tests using a
 Bonferroni-corrected alpha of 0.0009.  <em>f<sub>d</sub></em> was calculated in 50 kb non-overlapping genomic windows using the ABBABABAwindows.py Python script (https://github.com/simonhmartin/genomics_general) and smoothed across the genome using the <em>loess</em> function from the <em>scales</em> package (Wickham and Seidel 2022) in R with a span of 0.02. To determine the average level of admixed genetic variation between <em>aridus</em> and the red or yellow ecotypes using P1 taxa at different levels of taxonomic divergence, we estimated the mean <em>f<sub>d</sub></em> from each test across 50 kb windows. To estimate the average level of sequence divergence between the red or yellow ecotypes and the taxa from clades C and D, we calculated average <em>d<sub>a</sub></em> in 50 kb windows for each pair of taxa (Nei & Li 1979). <em>d<sub>a</sub></em> accounts for genetic divergence that predated species divergence by subtracting the average intraspecific pairwise differences (π in both species) from the observed interspecific value (<em>d<sub>xy</sub></em>). π and <em>d<sub>xy</sub></em> were calculated using PIXY version 1.2.5 (Korunes & Samuk 2021) with both variant and invariant sites included. To test for significant differences among the mean <em>f<sub>d</sub></em> values, we fit a linear model with <em>f<sub>d</sub></em> as the dependent variable and the P1 taxon as the independent variable using the <em>lm</em> function in R and then used the <em>emmeans</em> package (Lenth, 2019) to perform pairwise comparisons of the estimated marginal mean <em>f<sub>d</sub></em> values for the different P1s from the linear model. We estimated the Spearman’s rank correlation coefficient between <em>f<sub>d</sub></em>, <em>F<sub>ST</sub></em>, and π using the <em>cor.test</em> function in R. <em>F<sub>ST </sub></em>was calculated in each 50 kb window between the red or yellow ecotype and one of the five remaining taxa from clades C and D using the <em>popgenwindows.py </em>Python script (https://github.com/simonhmartin/genomics_general), with only variant sites included. Finally, we sorted the <em>f<sub>d</sub></em> values calculated in 50 kb windows into quantile bins based on recombination rates estimated in 500 kb windows in Stankowski et al. (2019). We then calculated the mean and 95% confidence intervals of the <em>f<sub>d</sub></em> values within each recombination rate quantile bin. To test for differences in <em>f<sub>d</sub></em> among recombination bins, we fit a linear model with <em>f<sub>d</sub></em> as the dependent variable and recombination bin as the independent variable and then perform pairwise comparisons of the estimated marginal mean <em>f<sub>d</sub></em> values for the different recombination bins from the linear model.  Admixture and genetic divergence across the <em>MaMyb2</em> region To determine the history of the red allele at <em>MaMyb2</em> in the different lineages of clade D, we calculated <em>f<sub>d</sub></em> and <em>d<sub>xy</sub></em> as described above, but in overlapping 10 kb windows with 100 bp steps. <em>f<sub>d</sub></em> was calculated three times, each time with <em>aurantiacus</em> as P1 and <em>aridus</em> as P3, but varying the different red samples (red ecotype, red-OC, and red <em>longiflorus</em>) as P2. <em>d<sub>xy</sub></em> was calculated between each of the red taxa and<em> aridus</em>, each of their closest yellow-flowered partners and <em>aridus</em>, between each of the red-flowered taxa, and between each of the red taxa and their most closely related yellow partners. Haplotype Network We used VCFx version: 2.0.6b to generate a FASTA file of the entire <em>MaMyb2</em> gene (positions 12,317,113 to 12,318,500 on chromosome 4) from the VCF file from all 47 samples. This FASTA file was then used to construct a haplotype network from all recovered haplotypes based on an infinite site model and uncorrected distances using <em>Pegas</em> version 1.1 (Paradis 2010) in R version 3.6.3.   Extended Haplotype Homozygosity We calculated the site-specific extended haplotype homozygosity statistic for a SNP in the middle of the <em>MaMyb2</em> gene (position 12,317,808) using <em>rehh</em> version 3.2.2 (Gautier & Vitalis 2012). Separate VCFs were created that contained samples from the red ecotype, the yellow ecotype, red <em>OC</em>, or red<em> longiflorus</em>. Haplotype and marker information was extracted from each VCF file using the <em>data2haplohh</em> function of the <em>rehh</em> package. EHH was calculated from the marker and haplotype information for each of the taxa using the <em>calc_ehh</em> function in <em>rehh</em>.