Shared Selection and Genetic Architecture Drive Strikingly Repeatable Evolution in Long-Term Experimental Hybrid Populations

Abstract

The degree to which evolution repeats itself has implications regarding the major forces driving evolution and the potential for evolutionary biology to be a predictive (vs. solely historical) science. To understand the factors that control evolutionary repeatability, we experimentally evolved four replicate hybrid populations of sunflowers at natural sites for up to 14 years and tracked ancestry across the genome. We found that there was very strong negative selection against introgressed ancestry in several chromosomes, but positive selection for introgressed ancestry in one chromosome. Further, the strength of selection was influenced by recombination rate. High recombination regions had lower selection against introgressed ancestry due to more frequent recombination away from incompatible backgrounds. Strikingly, evolution was highly parallel across replicates, with shared selection driving 88% of variance in introgressed allele frequency change. Parallel evolution was driven by both high levels of sustained linkage in introgressed alleles and strong selection on large-effect quantitative trait loci. This work highlights the repeatability of evolution through hybridization and confirms the central roles that natural selection, genomic architecture, and recombination play in the process.

Introduction

A central question in evolutionary biology is whether the path of evolution is primarily driven by deterministic processes such as natural selection and genetic constraint or stochastic processes such as historical contingency and genetic drift (Stern and Orgogozo 2008; de Visser and Krug 2014; Orgogozo 2015; Blount et al. 2018). If natural selection and constraint dominate, then short-term evolution can be a predictable, repeatable process (Bolnick et al. 2018). Understanding the forces that affect evolutionary repeatability is thus a critical problem in biology.

Evolutionary repeatability has been explored in the context of admixture between species during hybridization and introgression. Genomic analyses have shown that interspecies ancestry is seen in many lineages including fungi, plants, mammals, birds, and lizards (Schumer et al. 2013; Brandvain et al. 2014; Figueiró et al. 2017; Taylor and Larson 2019; Langdon et al. 2020; Yang et al. 2021; Owens et al. 2023). While all hybridization originates with a 1:1 mix in the F₁ generation, backcrossing or mating between hybrids produces offspring with segregating ancestry. Unbalanced population sizes, natural selection, and/or drift lead most admixed lineages to greater ancestry from one parental species (the major parent) and less from the other (the minor parent) (Green et al. 2010; Elgvin et al. 2017; Jones et al. 2018; Schumer et al. 2018). Where in the genome minor parent ancestry remains in an admixed lineage (i.e. introgression) has been studied in several lineages; understanding the evolutionary forces governing its distribution is key for predicting the outcome of future hybridization events (Sankararaman et al. 2014; Corbett-Detig and Nielsen 2017; Schumer et al. 2018; Martin et al. 2019; Duranton and Pool 2022; Nouhaud et al. 2022; Vilgalys et al. 2022; Langdon et al. 2024).

Studies so far have suggested that minor parent ancestry is typically deleterious either due to incompatibilities with the major parent genome or greater deleterious load (Harris and Nielsen 2016; Juric et al. 2016; Moran et al. 2021). Consequently, minor parent ancestry is more likely to be lost in regions of low recombination (Schumer et al. 2018; Martin et al. 2019; Calfee et al. 2021), because long minor parent haplotypes are more likely to contain multiple deleterious alleles and be purged more efficiently than shorter haplotypes that occur in regions of high recombination (Veller et al. 2023). Interestingly, in some scenarios where the minor parent contributes beneficial alleles, the reverse relationship between ancestry and recombination rate is predicted, but this has not been empirically shown (Duranton and Pool 2022). Supporting the general pattern of minor parent ancestry being removed due to selection, minor parent ancestry is typically reduced in regions with more functional or conserved basepairs (Brandvain et al. 2014; Juric et al. 2016; Schumer et al. 2018; Calfee et al. 2021).

These genomic patterns suggest that the outcome of hybridization may be a repeatable process. This has been tested in several systems. In Xiphophorus fishes, recent (<300 generations) hybrid populations show repeatable minor parent ancestry (Schumer et al. 2018; Langdon et al. 2024). Interestingly, the repeatability is higher in crosses between more diverged species, suggesting that the repeatability is driven by the selection on incompatibilities, which increase with evolutionary distance. For more diverged species, incompatibilities are likely to be stronger and more common. Similarly, in experimental or natural hybrid populations of Drosophila and ants, minor parent ancestry patterns were found to be highly repeatable (Brennan et al. 2019; Matute et al. 2020).

Here, we test the repeatability of hybrid genome evolution using four experimental hybrid populations of Texas sunflowers (Helianthus) planted in nature. Previous studies on phenotypes of these populations have shown the hybrid populations evolved increased fitness, while nonhybrid controls did not, and that hybrid populations had higher repeatability of phenotypic evolution than controls (Mitchell et al. 2019, 2022). In new genomic analyses, we specifically track parental haplotype frequency across the genome for up to 14 generations and find that repeatability is higher for minor parent haplotypes than for major parent haplotypes. We use quantitative trait locus (QTL) mapping of hybrid populations to show that the minor parent haplotype has stronger effect QTLs, including major incompatibilities, and lower recombination, both of which increase repeatability in simulations. Lastly, we show that the relationship between recombination and the amount of minor parent ancestry can be flipped depending on whether minor parent ancestry is selected for or against.

Results

We explored the repeatability of evolution using synthetic hybrid populations of sunflower grown in their natural habitat. A single F₁ hybrid between Helianthus annuus and Helianthus debilis was backcrossed to a second individual of H. annuus, and the resulting BC₁ offspring were planted in four locations across Texas and allowed to evolve naturally (Mitchell et al. 2019) (Fig. 1a and c). Previous work in these populations found that the experimental hybrids had repeatable phenotypic evolution (Mitchell et al. 2022). Since all individuals shared the same two parents (i.e. the F₁ and the H. annuus backcross parent), we were able to track the four individual parental haplotypes across generations. We collected leaf tissue for 7 to 14 generations per location and genotyped a total of 2,504 samples using genotyping-by-sequencing (supplementary table S1, Supplementary Material online). We used linkage patterns between markers in the BC₁ generation to identify parental haplotype-specific alleles and then called the diploid copy number of each parental haplotype across the genome for each sample using a hidden Markov model (HMM, ancestry_HMM v1.0.2) (Corbett-Detig and Nielsen 2017). We used this approach because single-nucleotide polymorphism (SNP)-based analysis suggested that outside gene flow brought in nonparental alleles (see Materials and Methods). Our HMM approach allowed us to quantify parental haplotype number across the genome and estimate the relative frequency of the four parental haplotypes while ignoring nonparental alleles, which would bias estimates of natural selection or repeatability. See Materials and Methods for details on haplotype identification. In addition, for the Katy Prairie Conservancy (KPC) and Houston Coastal Center (HCC) populations in generation 2 a significant portion of individuals did not contain outside alleles and were genetically pure. For a subset of tests, we repeated our analysis using SNPs from only these pure individuals, comparing generation 1 to 2.

To explore whether the replicate populations evolved in parallel, we first generated a principal component analysis of inferred parental markers every 1 Mbp using SNPRelate (Zheng et al. 2012) (Fig. 1b). We found highly similar population trajectories in all four replicate locations and that the largest changes occurred in the first few generations (supplementary fig. S1, Supplementary Material online). The most similar trajectories occurred in the Brackenridge Field Laboratory (BFL) and Lady Bird Johnson Wildflower Center (LBJ) locations, at 14.5 km apart the two most geographically and environmentally close lineages, although since we used four populations, we were unable to test the link between environmental similarity and evolutionary repeatability (supplementary fig. S2, Supplementary Material online). To better understand the contribution of natural selection to the parallelism observed, we used a method by Buffalo and Coop to quantify the percent variation in allele frequency change due to convergent natural selection (Buffalo and Coop 2020). Two populations evolving independently without shared selection will change allele frequencies independently, but shared selection will drive covariance in allele frequency changes. This method compares the average covariance in temporal allele frequency shifts across replicates with the total variance in allele frequency, while controlling for variance due to sampling. We selected the first and sixth generations as comparison points for consistency and because the largest changes in allele frequency occurred within the first six generations. We analyzed each parental haplotype separately and thus have a measure of contribution for both major and minor parent haplotypes. We found that the percent variation due to shared natural selection was 27% to 59% for the H. annuus haplotypes, but 88% for the introgressed H. debilis haplotype (Fig. 1d). While the parallelism for the H. annuus haplotypes is in line with previous estimates from wild Drosophila simulans populations (37%) or artificial selection on mice (32%) (Marchini et al. 2014; Castro et al. 2019; Kelly and Hughes 2019; Buffalo and Coop 2020), the introgressed H. debilis haplotype is likely experiencing stronger and more consistent natural selection and has a more repeatable evolutionary trajectory. When using diagnostic SNPs, we found higher variance explained by selection for all haplotypes (supplementary fig. S3, Supplementary Material online), but these values are inflated by outside gene flow reducing parental allele frequency genome wide. This issue is controlled by our HMM-based genotypes that reflect the frequency of the parental haplotypes excluding outside gene flow. The highly parallel genetic evolution in Fig. 1b recalls the parallel phenotypic evolution seen in these populations, in which hybrids showed stronger patterns of repeatability than nonhybrid controls (Mitchell et al. 2022).

What factors are causing shared selection to have greater effects on frequency change of alleles derived from the minor parent haplotype, relative to those from the major parent haplotypes? We hypothesize that there are two main reasons: limited recombination and larger QTL effect sizes, including reproductive incompatibilities. Natural selection acts on organisms and not individual alleles; therefore, the selection coefficient of an allele depends not only on its effect, but also on the genetic background in which it is found. This background will shift each generation due to segregation and recombination during meiosis. Genetic linkage is a measure of how often alleles co-occur and an allele with tight linkage to surrounding sites (i.e. limited recombination) is more consistently found in the same genetic background. We expect that higher consistency of genetic background will lead to more consistent selection and repeatable evolution. To test this hypothesis, we simulated two species that differed by a pair of Dobzhansky–Muller incompatibility (DMI) loci that are crossed to produce two BC₁ populations and evolved for six generations, mimicking our experimental design. We varied the recombination rate and measured the proportion of change in allele frequency due to natural selection as above. We found that lower recombination rates led to more variation being driven by natural selection (Fig. 2a). Consistent with our predictions, we found that linkage in the introgressed haplotype decayed more slowly than in the other haplotypes in our experimental data, which means reduced recombination, likely due to a combination of small- and large-scale chromosomal rearrangements, including known inversions and translocations (Chandler et al. 1986; Ostevik et al. 2020) (Fig. 2b).

Another potential reason for increased predictability in introgressed haplotypes is stronger selection. Indeed, QTLs with strong effects are more likely to evolve in parallel than QTLs with weak effects (MacPherson and Nuismer 2017). To test this prediction, we used the simulation detailed above and varied the effect size, in this case the selection coefficient, for the DMI loci. We found that larger effect size increased the proportion of alleles controlled by natural selection, and the predictability of evolution (Fig. 2c). Both simulations were repeated using underdominant incompatibilities in the place of DMIs and found similar results (supplementary fig. S4, Supplementary Material online). Based on these results, we predicted that the alleles of the introgressed H. debilis haplotype would have a larger phenotypic effect than alleles from the H. annuus backcross parent. To test this prediction, we QTL-mapped 25 phenotypes in 1,000 BC₁ plants grown in a common garden and genotyped using genotyping-by-sequencing (supplementary tables S2 and S3, Supplementary Material online). We mapped QTL using the H. annuus genetic map separately for each BC₁ parent using haplotype-specific markers identified earlier. We then compared the maximum effect size for each trait between the F₁ parent, containing the H. debilis haplotype, and the recurrent H. annuus parent, which contained only H. annuus haplotypes. We found that the F₁ parent had an average maximum QTL effect size 79% larger than the backcross parent (t(24) = 4.51, P = 0.00014), supporting the idea that QTL alleles originating from the introgressed H. debilis haplotype had much stronger QTL effects than the QTL alleles segregating within H. annuus haplotypes (Fig. 2d).

Although our simulations considered the effects of recombination and QTL effect size separately, they can have synergistic effects, especially if QTL effects are mostly in the same direction. That is, reduced recombination will combine the effects of many alleles, potentially leading to larger QTL effect sizes, and selection on such larger “combined” QTLs may contribute to the elimination of recombinant haplotypes.

While we expect that larger-effect QTLs are under stronger selection, we can also directly measure natural selection across the genome using our time-series data. We estimated the selection coefficient for every 1 Mbp across the genome for each parental haplotype and location using a linear least square regression (Taus et al. 2017). This revealed strong negative selection across most of the H. debilis haplotype for many chromosomes and consistent positive selection for one (Fig. 3a). The negative selection for the H. debilis haplotype on chromosomes 6, 12, 15, 16, and 17 is likely due to QTL for viable seed production found on the corresponding chromosomes (supplementary fig. S5, Supplementary Material online). In each case, the H. debilis haplotype reduces viable seed production in the heterozygote and is likely caused by chromosomal translocations between chromosomes 6 and 15, as well as chromosomes 12, 16, and 17. The translocations were identified using linkage information in the BC₁ generation and are consistent with rearrangements seen in related species (Ostevik et al. 2020) (supplementary fig. S6a, Supplementary Material online). The causes of reduced seed production could trace to a form of hybrid sterility caused by issues in chromosomal segregation during meiosis, as has been seen previously in sunflowers (Chandler et al. 1986; Quillet et al. 1995; Lai et al. 2005). We additionally see a weaker QTL for the same trait, as well as consistent negative selection on the H. debilis haplotype, on chromosome 4 which has a small translocation with chromosome 7 (supplementary fig. S6b, Supplementary Material online). Alternately or in addition, other forms of hybrid dysfunction could contribute, e.g. reduced acquisition of resources leading to lower seed production. The chromosome that experienced consistent positive selection (chromosome 2) is discussed below.

Helianthus annuus contains large polymorphic inversions on several chromosomes that have been linked to environmental adaptation and shown to differ between Texas H. annuus populations and the rest of the range (Todesco et al. 2020; Owens et al. 2021). To test whether inversion haplotypes are under selection in our experiment, we inferred inversion genotypes in the BC₁ generation using our sequence data. We then linked inversion and parental haplotype by testing cosegregation in inversion regions, identifying which inversion haplotype each parental haplotype contains. Lastly, we identified the inversion haplotype more common in Texas and limited our analysis to inversions that are both larger than 10 Mbp (allowing for accurate genotyping with genotype-by-sequencing [GBS]) and are in the top 2% most differentiated loci between Texas and the rest of range (Owens et al. 2021) (supplementary table S4, Supplementary Material online). With these data, we compared the average selection coefficient for all windows in each inversion for the three H. annuus parental haplotypes. We found that for all four inversions, a parental haplotype that matches the local Texas inversion haplotype had the most positive average selection coefficient (supplementary fig. S7, Supplementary Material online). This is consistent with a locally adaptive inversion haplotype, although it is also possible that contamination from local populations is artificially increasing the frequency of the matched parental haplotype. Our parental diagnostic markers are likely enriched for inversion haplotype-specific SNPs, which may make this error stronger in inversion regions.

Interestingly, although the introgressed chromosomal regions containing incompatibilities are almost entirely purged early on, the strength of negative selection is often reduced or reversed on the tips of the chromosomes, where recombination rates are highest (Fig. 3a and b). This occurs because in high recombination regions, neutral or beneficial H. debilis alleles are more likely to recombine onto a different background and become unlinked to the reproductive incompatibility QTL (Fig. 4b). There is a weak but positive relationship between recombination rate and the selection coefficient (P < e⁻¹⁶, percent variance explained = 0.5%), and this relationship is stronger when limiting to chromosomes with strong identified reproductive incompatibilities (P < e⁻¹⁶, pve = 11.2%) (Fig. 4d). Higher introgression in high recombination regions has been predicted through modeling and observed in hybrid populations (Duranton and Pool 2022), but here we observe the process in action. Strikingly, on chromosome 2 we observe consistent positive selection for the introgressed haplotype, suggesting adaptive introgression, and a negative relationship between selection and recombination rate (P < e⁻¹⁶, pve = 17.3%) (Fig. 4c). In this case, higher recombination regions have less positive selection because loci within them are more likely to become unlinked to the adaptive allele or alleles (Fig. 4a).

To corroborate our selection estimates and exclude the effect of outside gene flow, we measured the change in parental allele frequency (either diagnostic SNPs or parental windows) from generation 1 to 2 including only uncontaminated samples from the HCC and KPC lineages. We found a significant correlation between selection estimates and allele frequency change in parental windows (KPC: Pearson’s r = 0.72, P < e⁻¹⁶, n = 3,056, 95% CI [0.70, 0.73]; HCC: Pearson’s r = 0.19, P < e⁻¹⁶, n = 3,056, 95% CI [0.15, 0.22]). Using H. debilis diagnostic SNPs, we found a weak positive relationship between allele frequency change and recombination rate in all chromosomes (P = 4.7e⁻³, percent variance explained = 0.02%), a stronger relationship in chromosomes with strong reproductive incompatibilities (P < e⁻¹⁶, pve = 3.2%), and no relationship on chromosome 2 (P = 0.23). The lack of a positive relationship for chromosome 2 may be due to insufficient time for selection or recombination after only a single generation.

Although we measured several fitness-related traits during QTL mapping, including both reproductive and defense traits, the only QTL to map to chromosome 2 is for inflorescence disk diameter (supplementary fig. S8, Supplementary Material online). While the H. debilis value for this trait may be adaptive, it is one of several QTLs across the genome for the same trait so we think it is implausible that disk diameter is driving this pattern. Recent theory has suggested that selection against introgression may be mediated by stabilizing selection on phenotypic traits (Ragsdale 2024; Veller and Simons 2024). This is consistent with our findings that the minor parent haplotypes have larger QTLs and tend to be under negative selection, as would be predicted if there was stabilizing selection. Under this model, the fact that chromosome 2 is bereft of strong QTL may reduce negative selection on the H. debilis version of it, but on its own cannot explain the evidence of positive selection seen.

Discussion

Our works shows that evolution can be surprisingly repeatable not only at the phenotypic level, but also in terms of genotype (Mitchell et al. 2022). However, there are several reasons why the levels of repeatability observed in our experiment may be unusual. First, each population was initiated with the same alleles and at the same allele frequency. In natural experiments, the starting material is unlikely to be so uniform across lineage origins. Second, each population was founded from the progeny of just two individuals. Thus, a maximum of four alleles were segregating at a locus, which is expected to limit the range of adaptive responses, and reduces effective recombination. Third, because these populations were segregating for major chromosomal and ecological hybrid incompatibilities, selection was exceptionally strong. Lastly, strong selection, combined with structural differences in chromosomes between the parental species, reduced the extent of recombination between haplotypes, further enhancing repeatability.

Although the conditions in our experiment favored parallel evolution, we gained several important insights. Our work suggests that, as hypothesized, hybrid genotypic evolution is more repeatable than nonhybrid evolution due to inherent genetic constraints. Our work also supports the hypothesis that repeatability is positively correlated with evolutionary distance between genomes. That is to say, crosses between populations will have more repeatable outcome than crosses within a population, while between-species crosses will be more repeatable still (Langdon et al. 2024). We also confirmed previous predictions that limited recombination and large-effect QTL enhance repeatability (MacPherson and Nuismer 2017). Lastly, our work highlights that the early stages of evolution in a hybrid population are more predictable than later stages due to the combined effect of large QTL, high linkage, and the resolution of reproductive incompatibilities.

A positive correlation between recombination rate and minor parent ancestry is well established in numerous systems (Moran et al. 2021), but we believe this is the first experimental test showing a negative correlation under adaptive introgression. For chromosome 2, minor parent alleles in high recombination regions on the chromosome ends were more likely to recombine off the adaptive haplotype. Although we measured fitness-related traits including growth, defense against herbivory, and seed production, we found no plausible QTL explaining positive selection for the minor parent haplotype on chromosome 2. Positive selection may arise from QTLs for an unmeasured fitness trait (e.g. pollen production), or it could be due to a transmission distorter like meiotic drive (Lindholm et al. 2016; Finseth 2023). Due to the lack of recombination in the minor parent, the region containing the adaptive QTL or QTLs includes most of the chromosome, so we were unable to identify candidate genes.

In this work, we tracked the evolution of individual haplotypes, while previous repeatability studies tracked ancestry, which included many individual haplotypes per ancestry class. This likely increased the repeatability we measured here, since each copy of the haplotype was identical, barring de novo mutations, while ancestry estimates include segregating variation within the ancestry class. Although our work focused on evolutionary repeatability during hybridization, the principles tested here are relevant to short-term evolution in any sexual population. For example, based on our results, we expect more predictable evolution of alleles in low recombination regions like centromeres. Future studies using new methods to fully phase genomes (Cheng et al. 2021; Li and Durbin 2024) combined with ancestral recombination graph methods (Rasmussen et al. 2014; Brandt et al. 2024) may be able to track the repeatability of haplotypes in more complicated experimental designs or even natural populations.

Materials and Methods

Establishment of Experimental Hybrid Populations

Our experimental populations were started with a cross between a single H. debilis from Texas and a wild H. annuus from Oklahoma to create the F₁ generation (see details in Whitney et al. 2006). A single F₁ plant was vegetatively propagated to produce 14 F₁ clones, which were mated to a single H. annuus pollen donor from north Texas (Fig. 1a). Seeds from this cross were germinated on filter paper, transplanted to peat pots, and grown in the greenhouse for 1 month before transplantation to the field. We planted 500 seedlings in each of four different locations in Texas: LBJ (30.184°N, −97.877°W), BFL (, 30.282°N, −97.780°W), HCC (29.39°N, −95.04°W), and KPC (29.9232°N, −95.9236°W) (Fig. 1c). The first two plantings occurred in 2003, while the latter two plantings occurred in 2008.

Populations were allowed to evolve naturally except for rototilling the soil each winter, as annual sunflowers are early-successional species and require disturbance to ensure sufficient germination. Additionally, each year wild sunflowers were removed within 250 m to reduce outside gene flow. For most generations, we collected seed for common garden phenotype measurements (Mitchell et al. 2019, 2022) and leaf tissue for genetic analysis. We attempted to track each population for as many years as possible, until changes in land-use permission (LBJ) and funding constraints (all other sites) terminated data collection. While our earlier exploration of repeated phenotypic evolution (Mitchell et al. 2022) focused on the BFL, LBJ, and HCC populations, here we add consideration of a fourth population (KPC) for which we had genotypic but not phenotypic data.

To quantify environmental distance between population locations, we extracted variables #1 to #19 from WorldClim v2.1 (Fick and Hijmans 2017). We then calculated pairwise distance including variables between each population using the dist() function in R.

Sample Preparation and Sequencing

We sequenced a total of 1,000 individual BC₁ samples grown BFL and LBJ using single-enzyme GBS as well as 1,313 BC₁ and advance generation experimental hybrids from all four locations using two-enzyme GBS (Elshire et al. 2011; Poland et al. 2012) (supplementary table S1 and supplementary fig. S9, Supplementary Material online). A subset of advance generation samples were sequenced using both GBS protocols, and the data were pooled before genotype calling. For all samples, DNA was extracted from 40 mg samples of fresh leaf tissue using Qiagen DNeasy 96 plant kit. Sequencing libraries of 96 barcoded samples were created using the following protocol. For samples with <100 ng of genomic DNA, whole-genome amplification for 6 h was performed using the Qiagen Repli-g kit first. For single-enzyme GBS, we digested 100 ng of DNA with the endonuclease PstI-HF (NEB Inc., MA, USA). Fragments were ligated to adaptors including the Illumina sequencing adaptor and individual barcodes with T4 ligase (NEB Inc.). Phusion high-fidelity polymerase (NEB Inc.) was used for PCR amplification. Size selection retained amplified fragments of length between 300 and 500 bp. Further DNA cleaning was done using Sera-Mag SpeedBeads Carboxyl Magnetic Beads (Cytiva, MA, USA). The libraries were run on an Illumina HiSeq 2000 pair-ended 100-bp platform at Genome Quebec sequencing facility pooling 192 samples per lane of sequencing.

For two-enzyme GBS libraries, we extracted DNA as above, but the digestion protocol used both PstI-HF and MspI (NEB Inc.). Kapa HIFI Hotstart was used for PCR amplification. To reduce the representation of repetitive sequences in the libraries, we performed a depletion step by treating the enriched libraries with duplex-specific nuclease (Evrogen, Moscow, Russia). As before, we pooled 192 samples per lane but instead sequenced the libraries on an Illumina HiSeq 2500 pair-ended 125-bp platform at Genome Quebec sequencing facility.

Parentage Analysis

We demultiplexed and trimmed reads using a custom perl script (Owens et al. 2016). Reads were then aligned to the HA412HOv2 H. annuus genome using NextGenMap (v0.5.5) (Sedlazeck et al. 2013) and processed using samtools (v1.10) (Danecek et al. 2021). We then called variants together using FreeBayes (v1.3.4) for all samples (Garrison and Marth 2012). This sample set included primarily three types of samples: (i) multiple generations of experimentally evolved hybrid populations prepared using two-enzyme GBS, (ii) first-generation backcross samples phenotyped for traits prepared using single-enzyme GBS, and (iii) a small number of reference local H. annuus, and H. debilis samples prepared using single-enzyme GBS. All backcross and experimental hybrid samples were derived from a single BC₁ family and share four parental haplotypes. Tissues from the original BC₁ parents were lost during a lab move, so we instead used the reference genome and inheritance patterns in the BC₁ population to identify diagnostic alleles in each parental haplotype.

To do this, we filtered our initial variant table to only retain BC₁ samples sequenced using the two-enzyme GBS and required each called genotype have read depth ≥ 5 and each site have ≤50% missing data. Markers unique to a single parental haplotype should have 25% minor allele frequency, so we then filtered the dataset to retain only sites with minor allele frequency ≥ 15% and ≤ 35%. This variant dataset was divided by chromosome, reformatted using perl scripts, and then loaded into rqtl (v1.6) (Broman et al. 2003; Arends et al. 2010). We then did a second round of filtering to remove samples with ≥50% missing data and then markers with ≥50% missing data. We then used the est.rf and formLinkageGroups (max.rf = 0.1, min.lod = 16) functions to build linkage groups. Since all markers were from a single chromosome (barring rearrangements in the H. debilis genome), the linkage groups represent the four parental haplotypes. In most cases, we observed four major linkage groups of which two pairs were in repulsion phase, consistent with a BC₁ (supplementary fig. S10a, Supplementary Material online). In a minority of cases, additional linkage groups were observed with unusually small recombination fraction to one of the four core linkage groups. In this case, we manually merged linkage groups to achieve four groups representing the four parental chromosomes.

At this point, we had identified four haplotypes but did not know which haplotype was from the H. debilis parent. To identify this haplotype, we used a set of previously published whole-genome sequence data for H. annuus and H. debilis (Owens et al. 2023). Sunflowers have high amounts of shared variation between species and relatively few fixed differences; therefore, we identified debilis-like markers where the alternate allele was common (allele frequency ≥ 50%) in H. debilis and not common (allele frequency < 50%) in H. annuus. This set of debilis-like alleles was matched against markers in each linkage group, and in all cases, one of the four major linkage groups had a distinctly higher proportion of debilis-like markers and was dubbed the deb haplotype (supplementary fig. S10b, Supplementary Material online). The haplotype in repulsion phase with H. debilis was set as ann1, and the two remaining haplotypes were set as ann2 and ann3. At the end, we kept 10,080, 4,596, 4,114, and 3,090 haplotype-informative markers for the deb, ann1, ann2, and ann3 haplotypes, respectively. We recognize that ann2 and ann3 haplotypes were not consistently named across chromosomes, due to the lack of linkage between chromosomes.

QTL Mapping

Phenotypes were measured for 500 BC₁ plants in each of two common gardens, one at BFL and the other at LBJ, and were combined for analysis. Further details on the gardens are given in Whitney et al. (2015), while details on the traits measured are given in Whitney et al. (2006, 2010). To directly link phenotypes with parental haplotypes, we used diagnostic markers for parental haplotypes. For each marker, we treated all homozygous alternate (1/1) genotypes as heterozygous (0/1) because 1/1 genotypes are not possible in a BC₁ for our markers and are likely to be due to allele dropout. We combined marker sets for haplotypes in repulsion phase (deb + ann1, ann2 + ann3) and flipped genotypes (0/0 to 0/1, 0/1 to 0/0) for one haplotype of each pair. Genetic map positions were set using the H. annuus genetic map and marker position on the H. annuus reference genome (Todesco et al. 2020; Huang et al. 2023). Genotypes were filtered to require at least three reads, and we additionally filtered markers to require they have between 20% and 80% heterozygous genotypes and be genotyped in > 40 samples. The relatively loose filtering is because of the poor overlap of sequence coverage between one- and two-enzyme GBS protocols. Markers were then imputed using the sim.geno (step = 0, n.draws = 100, err = 0.05) in rqtl (Broman et al. 2003). Logarithm of the odds (LOD) scores were generated using the scanone function for each phenotype. We generated LOD thresholds using 200 permutations and recorded the 5% and 10% thresholds. QTL effect sizes were estimated using the effectscan function. We used a 95% Bayesian threshold to determine the size interval for significant QTL.

For the deb/ann1 parents, we found at least one significant QTL for all traits except leaf chewing damage (supplementary fig. S8, Supplementary Material online). In general, we found that QTL effect sizes were in the direction expected by the differences between species. For example, H. debilis alleles led to smaller disk diameter, earlier flowering time, lighter seeds, and reduced glandular trichome density. Interestingly, we also observed strong QTL for seed production on chromosomes 6, 12, 15, 16, and 17. In each case, the H. debilis allele led to lower seed production. For the ann2/ann3 cross, we found fewer and weaker QTLs.

Hybrid Population Evolution

In the absence of outside gene flow, the experimental hybrid populations should only contain the alleles present in the initial BC₁ population as well as a very small number of new mutations. To test this, we filtered to genomic positions that were sequenced in at least 25 BC₁ samples and asked if advance generation samples had at least 1,000 novel alleles. This threshold value was determined based on examining the distribution of novel alleles in all samples. We found that the proportion of samples with novel alleles increased over time and by generation 3 most samples had at least some novel alleles and likely outside haplotypes (supplementary fig. S11, Supplementary Material online). This suggests that although populations were grown in areas that were at least 250 m away from local sunflower populations, this isolation was not absolute for the length of the experiment, and outside pollen entered populations early on.

Outside gene flow complicated our analyses of selection in two ways. Outside gene flow acts as migration into the population and reduces the allele frequency of parental alleles in the absence of selection. As gene flow was not expected, we have no way of knowing the amount of gene flow per generation and controlling for it. Additionally, outside haplotypes may contain parental specific alleles, and so, migration can cause a parental allele to decrease in frequency, even if the actual parental haplotype is not under negative selection. This is seen in measurements of allele frequency changes over time for parental specific alleles. We found that there were consistent changes in allele frequency for long stretches of the chromosome, but a small minority of markers had very different allele frequencies (supplementary fig. S12, Supplementary Material online). While this may be due to multiple recombination events unlinking diagnostic alleles from their original haplotypes, the fact that these are often single markers and not blocks suggests that it is more likely to be due to outside haplotypes sharing diagnostic markers.

To overcome the uncertainty of individual markers, we used the program ancestry_HMM to identify parental ancestry blocks within the genome (Corbett-Detig and Nielsen 2017). Ancestry_HMM was designed to identify ancestry blocks in recently admixed populations. For each parental haplotype, we subset our variant table to only diagnostic markers for that parent. We then converted the VCF to ancestry_HMM read depth format using a custom perl script. We coded markers as allele 0 for the diagnostic allele and 1 for all other alleles. Genetic map positions for each marker were imputed from the H. annuus genetic map (Todesco et al. 2020). Ancestry_HMM requires allele frequencies for source populations. For this, we assumed that population 1, representing the parental haplotype, had an allele frequency of 100%, while population 2, representing all other possible sources, had an allele frequency of 5%. This recognizes that although the alleles are diagnostic within the four parental haplotypes, they may occur from outside gene flow at a low probability. Each population and generation was run separately, and initial ancestry proportions were estimated based on diagnostic allele counts. We modeled two ancestry pulses, one 10,000 generations ago with 75% contribution representing the other three nontargeted parents and one N generations ago contributing 25%, where N is the hybrid population generation + 2. We used forward–backward decoding and called ancestry when the likelihood was ≥50% for one state and unknown if otherwise. This was translated into genomic blocks based on the likelihood at markers. Where there was transition between states, we chose the midpoint between markers as the boundary of genomic blocks. For each block, the sample has 0, 1, or 2 copies of the parental haplotype.

Since ancestry_HMM inferred parentage for different sets of markers for each parental haplotype, the genomic blocks had different boundaries for each parent. To facilitate easier comparisons of parental composition, we picked positions every 1 Mbp and called parental ancestry based on genomic blocks. We named these parental markers. In this way, we had parental haplotype counts for all possible haplotypes at consistent positions across the genome. In some cases, parentage informative markers did not cover the entire chromosome for all haplotypes (i.e. some were missing diagnostic markers at the distal chromosome tips) so we limited analysis to regions covered for all parents.

If parentage estimations are accurate, we expect to see evidence of recombination between paired haplotypes in ancestry composition in BC₁ samples. We checked this by plotting parentage across the genome (supplementary fig. S13, Supplementary Material online). We found that BC₁ samples often showed clear switches between parental haplotypes as expected by recombination. This is notable because we ran each parental haplotype separately, so this represents agreement between different sets of parental SNPs. We expect that most of the genome should have two haplotypes when all parental copies are summed at any one position. Outside gene flow should not be interpreted as any parental haplotype, so it should produce regions with < 2 haplotype copies. Errors in ancestry inference, including allele dropout, could cause increases or decreases in haplotypes. For each parental marker, we counted the total number of haplotypes and compared this across locations and generations.

We found that for the BC₁ generation nearly all of the genome had two haplotype copies. In further generations, we found higher proportions of regions with fewer than two haplotypes (supplementary fig. S14, Supplementary Material online). We summed the total number of haplotype copies across the genome and divided by two times the number of regions to estimate the proportion of the genome contributed by outside haplotypes and found that outside ancestry was lower in KPC and parental haplotypes averaged ∼60% of the genome in advanced generations (supplementary fig. S15, Supplementary Material online).

Parallel Evolution

We explored how haplotype frequencies changed during experimental evolution by calculating haplotype frequencies based on the 1-Mbp spaced haplotype markers derived from ancestry_HMM. We were primarily interested in the relative change in frequency of the four parental haplotypes, so we calculated frequency only considering those haplotypes. This means that the total number of alleles at any one position may be less than two times the number of samples since some outside genotypes were not counted. For each marker and location, we used a linear least squares method in PoolSeq to estimate the selection coefficient (Taus et al. 2017). For the selection estimates, we set the population size as the harmonic mean of the recorded census sizes for each population, excluding 1 year for KPC which had a census size of zero. The 95% CI of selection estimates was calculating using 1,000 simulations. To explore the relationship between recombination rate and selection on the introgressed haplotype, we used a three-factor ANOVA with the formula s = location + chromosome + recombination rate in R (R Development Core Team 2013; Wickham et al. 2019). We then repeated this using only chromosomes with strong identified reproductive incompatibilities (chromosomes 6, 12, 15, 16, and 17), and again only using chromosome 2 which was under consistent positive selection.

To validate that selection estimates were not driven by outside gene flow effects, we selected a set of “pure” individuals that had<1,000 nonparental alleles. Only generation 1 and generation 2 in HCC and KPC had enough pure individuals to estimate allele frequency. From these samples, we estimated the change in allele frequency between generations 1 and 2 for both populations, both for diagnostic SNPs and for the HMM-estimated parental markers. We then calculated Pearson's correlation between our measure of allele frequency change and our selection estimate for parental markers. Following this, we repeated our previous ANOVA replacing selection for H. debilis SNP allele frequency change from generation 1 to 2.

Parallelism was calculated using the convergence correlation method of Buffalo and Coop (2020). Specifically, we calculated the proportion of the total variance in allele frequency change from convergent selection pressure. This method uses the average covariance of allele frequency changes divided by the total variance. Since each population had different numbers of generations sampled, we chose generation 1 and generation 6 as comparison points. For each population, parent and marker, we calculated the change in haplotype frequency (as described above) and then calculated variance and covariance in allele frequency change between locations. CIs were calculated using 250 bootstrap replicates.

To confirm that our method of inferring parental haplotypes was not inflating repeatability estimates, we estimated the percent variance explained by selection on parentage informative SNPs (supplementary fig. S3, Supplementary Material online).

Chromosomal Inversions

We extracted genotypes that overlapped with diagnostic markers for the chromosomal inversions ann13.01, ann13.02, ann14.01, and ann15.01 (Todesco et al. 2020). Based on these markers, we assigned inversion genotypes to each BC1 individual sequenced using the two-enzyme GBS protocol following the same protocol as Todesco et al. We then chose a window in the middle of the inversion and plotted the counts for inversion genotype and parental genotype. For example, for ann13.01, when the ann1 parental haplotype was present, inversion genotypes were always 0/1 or 1/1 and never 0/0. This suggests that ann1 contained the “1' inversion haplotype. We repeated this with each inversion and parental genotype to fully assign inversion haplotypes to each parental haplotype. To identify which inversion allele was the local Texas version, we compared population allele frequencies from Todesco et al. To characterize selection on inversion haplotypes, we calculated the average and SE selection coefficient for all windows within the inversion for each H. annuus parental haplotype.

Linkage Decay

To detect how linkage within the parental haplotypes decays, we subset out parental diagnostic markers and calculated genotype correlation coefficients between all markers within a chromosome for each population and generation using vcftools (v0.1.16) (Danecek et al. 2011). We then grouped values by the distance between markers in 5-Mbp bins and calculated the average and SE for correlation coefficients in each bin for each parent, location, and generation. We found that linkage disequilibrium (LD) decayed at increased physical distance and over generations, although notably less so for the H. debilis haplotype.

Interchromosomal LD

To confirm chromosomal translocations, we measured interchromosomal LD between all H. debilis-specific SNPs in the BC₁ samples. From that, we took the highest LD value between each pair of chromosomes and found strong LD between chromosomes 6 and 15 as well as between chromosomes 12, 16, and 17 (supplementary fig. S6, Supplementary Material online). Additionally, we visualized pairwise SNP LD between chromosomes 4 and 7 for H. debilis diagnostic sites, which showed a localized pattern of interchromosomal LD.

Simulations

We used non-Wright Fisher, forward-time genetic simulations in SLiM (v3.7.1) to explore how recombination rate and QTL effect sizes affect the repeatability of genotype evolution (Haller and Messer 2019). In the simulation, we modeled a simple two-locus DMI. Each species was fixed for one derived DMI allele, and derived DMIs negatively interacted in a partially dominant manner (supplementary table S5, Supplementary Material online).

The genome was set to 200,000 bp, and the two DMI loci were evenly spaced at 50,000 and 150,000 bp. Recombination rate was set such that the entire chromosome was 100 cM long and carrying capacity was set to 275, based on harmonic mean population size in the experimental populations. In the simulation, we created a BC₁ population between the two species and then allowed it to evolve for six generations. We then extracted the allele frequency for 100 species-specific neutral markers evenly spaced across the genome. This was repeated, and between the two independent simulations, we estimated the proportion of variation in neutral markers due to selection. This was repeated 100 times per parameter value. To test for the effect of effect size, we varied the BDM selection coefficient from 0 to 0.9, at 0.1 intervals. To test for recombination rate, we set the BDM selection coefficient to 0.5 and varied the recombination rate such that the genome was 0.5 to 1,000 cM long.

We repeated the simulations replacing BDM loci with underdominant loci to better mimic a chromosomal translocation. For each locus, heterozygotes received a fitness penalty (s), while homozygotes of either type had equal fitness. All other aspects of the simulation were retained between models.

Supplementary Material

Supplementary material is available at Molecular Biology and Evolution online.

Acknowledgments

Special thanks to Rebecca Randall, Steve Hovick, and Loren Albert for their many contributions. Thanks to many dozens of undergraduate assistants and research technicians for assistance in the field and laboratory over the years. The authors thank the HCC, Ladybird Johnson Wildflower Center, the KPC, and BFL for access to plot space and maintenance.

Funding

This work was funded by the National Science Foundation DEB-0716868 and DEB-1257965 to K.D.W. and L.H.R. and by the NSERC Discovery to G.L.O. This research was enabled in part by support provided by the BC DRI Group and the Digital Research Alliance of Canada (alliancecan.ca).

Data Availability

All sequence data are available on the Sequence Read Archive under project #PRJNA445853. All codes used in this project are available on GitHub (https://github.com/owensgl/texanus_repeatability).

References