Recurrent Loss of Sex Is Associated with Accumulation of Deleterious Mutations in Oenothera Free

Abstract

Sexual reproduction is nearly universal among eukaryotes. Theory predicts that the rarity of asexual eukaryotic species is in part caused by accumulation of deleterious mutations and heightened extinction risk associated with suppressed recombination and segregation in asexual species. We tested this prediction with a large data set of 62 transcriptomes from 29 species in the plant genus Oenothera, spanning ten independent transitions between sexual and a functionally asexual genetic system called permanent translocation heterozygosity. Illumina short-read sequencing and de novo transcript assembly yielded an average of 16.4 Mb of sequence per individual. Here, we show that functionally asexual species accumulate more deleterious mutations than sexual species using both population genomic and phylogenetic analysis. At an individual level, asexual species exhibited 1.8× higher heterozygosity than sexual species. Within species, we detected a higher proportion of nonsynonymous polymorphism relative to synonymous variation within asexual compared with sexual species, indicating reduced efficacy of purifying selection. Asexual species also exhibited a greater proportion of transcripts with premature stop codons. The increased proportion of nonsynonymous mutations was also positively correlated with divergence time between sexual and asexual species, consistent with Muller’s ratchet. Between species, we detected repeated increases in the ratio of nonsynonymous to synonymous divergence in asexual species compared with sexually reproducing sister taxa, indicating increased accumulation of deleterious mutations. These results confirm that an important advantage of sex is that it facilitates selection against deleterious alleles, which might help to explain the dearth of extant asexual species.

Introduction

Sexual reproduction is virtually ubiquitous among higher eukaryotes, despite significant costs in comparison to asexual reproduction (Williams 1975; Maynard Smith 1978). The rarity of asexual eukaryotes has spurred much debate on the advantages of sexual over asexual reproduction (Weismann 1889; Maynard Smith 1978; Bell 1982; Fisher 1930). Theory predicts that asexuality may reduce the efficacy of natural selection, leading to accumulation of deleterious mutations and/or stochastic loss of the fittest genotypic classes (Muller’s ratchet), and reducing genomic rates of adaptive evolution (Muller 1964; Felsenstein 1974; Keightley and Otto 2006). Asexual lineages may therefore have higher extinction risk compared with sexual populations (Lynch et al. 1993), potentially explaining the rarity of asexual species. Numerous theoretical studies support these predictions, but empirical investigations have proved more difficult (Paland and Lynch 2006; Neiman and Taylor 2009; Tucker et al. 2013).

A potentially informative approach for testing these predictions is to compare the abundance of mildly deleterious alleles in closely related sexual and asexual populations or species (Neiman and Taylor 2009). A reduction in the efficacy of selection is predicted to increase the proportion of mildly deleterious alleles, such as nonsynonymous single nucleotide polymorphism (SNPs), segregating within populations or accumulating as fixed differences between species (Gabriel et al. 1993). This approach has been adopted in a handful of systems, including microcrustaceans, mollusks, plants, and insects (Barraclough et al. 2007; Johnson and Howard 2007; Neiman et al. 2010; Henry et al. 2012; Hersch-Green et al. 2012). However, previous studies have examined only a small number of genes and/or taxa, and thus the genome-wide effects of transitions to asexuality, and the timing over which these changes occur, remain uncertain (Johnson and Howard 2007; Neiman et al. 2010).

Ideally, a comparative investigation of this question would make use of a system in which a loss of sex has occurred repeatedly in independent lineages over a range of timescales, allowing for an evaluation of the time-dependent effect of asexuality over multiple comparisons (Muller 1964; Gabriel et al. 1993). Based on these criteria, the plant genus Oenothera (Onagraceae) provides a powerful system for investigating the evolutionary consequences of sex (Ranganath 2008; Johnson et al. 2009). The ancestral state of the genus is sexual reproduction, and sexual species are typically highly outcrossing (Stebbins 1950; Raven et al. 1979). However, a significant proportion (29%) (Johnson et al. 2011) of Oenothera species have evolved to be functionally asexual due to a genetic system called Permanent Translocation Heterozygosity (PTH). PTH evolves through reciprocal chromosome translocations, and is characterized by nearly complete suppression of recombination, and balanced lethal mutations that prevent segregation during meiosis. Combined with self-fertilization, PTH results in functionally asexual individuals (i.e., progeny are genetically identical to parents) that transmit two independently evolving haploid genomes (Cleland 1972; Rauwolf et al. 2008). Importantly, closely related PTH (hereafter asexual/asex) and sexual Oenothera species typically have similar life-history, morphology, and ecological requirements, and in many cases they exhibit overlapping ranges. Moreover, PTH and sexual species are typically diploid, whereas in other systems asexuality is often associated with polyploidy. Therefore, the PTH system allows us to test the consequences of suppressed recombination and segregation specifically, without confounding effects of ecological differences (Ranganath 2008).

Here, we present the first genome-wide analysis of polymorphism and divergence at protein-coding genes from 29 species of Oenothera, spanning ten independent sexual/asexual transitions. This data set allows for an unprecedented opportunity to examine the consequences of loss of sex, replicated across thousands of genes and a phylogeny. We test the hypothesis that loss of sex is associated with reduced efficacy of purifying selection, leading to an accumulation of deleterious mutations in asexual lineages.

Results and Discussion

Sampling of Individuals and Short-Read Sequencing

We generated 62 leaf-tissue transcriptomes from 29 Oenothera species using Illumina short-read sequencing, encompassing at least ten independent transitions from sexual to functionally asexual reproduction (fig. 1; supplementary table S1, Supplementary Material online). On average, we sequenced 2.5 billion nucleotides per individual. For all sequenced individuals, we produced de novo assemblies of mRNA transcripts. Individuals had on average 22,098 transcript assemblies of at least 300 nucleotides in length and 16.4 Mb of total assembled sequence, providing a genome-wide view of genetic diversity in protein-coding genes. Our sample included within- and among-species variation, allowing us to investigate both the short-term (population genomic) and long-term (phylogenomic) consequences of asexual reproduction.

Patterns of Individual Heterozygosity

We first examined how transitions between sexual and asexual reproduction affected within-individual heterozygosity and within-species polymorphism. Asexual species are expected to have elevated heterozygosity, both because nucleotide differences inherited from sexual progenitors are fixed in a heterozygous state, and a lack of segregation prevents coalescence between founding alleles. Asexual Oenothera species are self-fertilizing—a mating system that typically leads to near-complete homozygosity because segregation allows sampling of identical alleles from selfing individuals (Nordborg 2000). However, if functional asexuality in Oenothera prevents recombination and segregation across much of the genome, we expect to observe high levels of heterozygosity throughout the transcriptome of asexual species.

To genotype individuals, we selected a single reference transcriptome from each of seven closely related clades in our data set, with each clade containing one or more transitions between sexual and asexual reproduction (fig. 1). We then mapped reads from all individuals within a given clade to the reference, and called SNPs simultaneously for all individuals within each clade. We implemented a maximum-likelihood method to exclude SNPs resulting from spurious mapping of reads from multiple transcripts to single reference loci (See Materials and Methods) (Gayral et al. 2013).

Consistent with the expectation for asexuality, the genome-wide proportion of heterozygous SNPs in asexual individuals was on average nearly twice as high as in sexual individuals (fig. 2a; 22% [n = 34] vs. 12% [n = 28]; t-test P < 0.002). Our Oenothera sample also included within-species sampling of four independent sexual/asexual species pairs (table 1), allowing us to examine the contribution of heterozygosity to within-species polymorphism. For each sexual/asexual pair, we calculated F_IT, a measure of the amount of inbreeding within individuals relative to the total population (table 1). For all comparisons, F_IT was lower among asexual species, reflecting a greater contribution of heterozygosity to total polymorphism compared with sexual sister species.

Within-Species Polymorphism

If the evolution of asexuality is associated with an ongoing reduction in the efficacy of natural selection to purge deleterious mutations, we expect to see an accumulation of segregating nonsynonymous (amino acid changing) polymorphisms (θ_N) compared with synonymous polymorphisms (θ_S) in asexual species. To test this prediction, we annotated protein-coding regions within transcripts, and compared the θ_N/θ_S ratio between four pairs of sexual/asexual sister species for which we had within-species sampling (n = 3–5 per species; table 1). Transcriptome-wide, there was a significantly higher θ_N/θ_S ratio in asexual species for three of four comparisons (fig. 2b; clade 4 nonsignificant; clade 7 P < 10⁻³; clade 2 P < 10⁻¹⁰; clade 3 P < 10⁻⁷). The single exception was clade 4, which exhibits very low sexual/asexual species divergence (0.01%), high inbreeding (F_IT = 0.8), and low effective population size (θ_S = 0.0046 ± 0.0003) in the sexual species (Oenothera grandiflora). These patterns likely reflect the fact that the sexual O. grandiflora is geographically restricted in southeastern United States where it occurs in small fragmented populations. By contrast, its asexual sister taxon O. biennis occurs throughout eastern North America (Dietrich et al. 1997).

Importantly, the observed increase in θ_N/θ_S was positively correlated with the degree of sequence divergence between sexual/asexual pairs (fig. 2b), and with levels of individual heterozygosity as a proportion of segregating sites in the asexual species (supplementary fig. S1, Supplementary Material online, R² = 0.85 and 0.80 for divergence and heterozygosity, respectively). As both measures provide rough estimates of the age of the transition to asexual reproduction, these patterns are consistent with a time-dependent process in which asexuality decreases the efficacy of selection against mildly deleterious alleles (Charlesworth 2012). We acknowledge that this pattern is based on a small number of comparisons (N = 4) and additional sampling within and/or between species (see below) would strengthen this conclusion. Furthermore, higher levels of heterozygosity in some asexual lineages could in part be due to their hybrid origins (see below).

To further infer the extent to which higher θ_N/θ_S ratios in asexual species reflect less efficient purifying selection, we fit a model of the distribution of fitness effects of deleterious nonsynonymous variants for the four sexual/asexual pairs, using polymorphisms at 4-fold degenerate sites as a proxy for neutral evolution. Consistent with the θ_N/θ_S patterns, we inferred an increase in the proportion of effectively neutral (N_eS ≤ 1) nonsynonymous variants in asexual species compared with their sexual sister species for all comparisons except clade 4 (supplementary fig. S2, Supplementary Material online). This indicates that a significant proportion of sites under efficient purifying selection in sexual clades are effectively neutral in asexual lineages, consistent with less efficient purifying selection in the latter.

An additional prediction from the hypothesis of reduced efficacy of selection in asexuals is that deleterious mutations of large effect, such as premature termination codons, should increase in abundance following loss of sex. Although such mutations are under strong purifying selection in sexual populations, a lack of segregation in asexual species ensures that a functional gene copy is consistently inherited, possibly leading to relaxed selection against such mutations. To test this, we quantified the occurrence of premature termination codons as a proportion of polymorphic sites in sexual and asexual species; this proportion was significantly higher in asexuals on average (fig. 2c; P < 0.02). Taken together, these results are consistent with a genome-wide reduction in the efficacy of purifying selection and accumulation of deleterious nonsynonymous mutations in asexual species.

Between-Species Divergence

To investigate whether the observed reduction in the efficacy of purifying selection translates into an accumulation of deleterious mutations across species, we determined orthology for a core set of 4,081 de novo assembled transcripts across all individuals in our data set. We aligned orthologs, inferred locus-specific phylogenetic trees using maximum likelihood, and estimated the ratio of nonsynonymous to synonymous divergence between species (dN/dS = ω) (Yang and Nielsen 1998). To account for possible heterogeneity in ω across the Oenothera phylogeny, for example, from different levels of divergence between asexual species and their sexual progenitors, we fit a “free-ratio” model that estimated ω separately for each locus along each branch. We then compared the distributions of ω along terminal branches leading to sexual/asexual species separately for the seven clades.

In all seven clades, average ω across the 4,081 loci was higher in asexual than sexual lineages, and this difference was statistically significant for four of the seven comparisons (fig. 2d; t-test clade 2: P < 1.0 × 10⁻⁶; clades 1 and 6: P < 0.04; clade 5: P < 0.05). Although the P value of some of these tests is not particularly low, the binomial probability of obtaining four of seven tests with P less than 0.05 by chance is exceedingly unlikely (Binomial expansion test: P = 0.0002; Zar 1984). Furthermore, when we focus on the three clades for which ω was not significantly elevated, and consider only those loci that have accumulated a comparatively larger number of mutations since divergence between asexual and sexual species (i.e., loci with dS > median across loci), we find that ω was also significantly elevated within asexual species in clades 3 and 7 (P < 0.05 and P < 0.03, respectively). Therefore, the hypothesis of greater ω in asexual species is conditionally supported within these clades as well. Importantly, ω was never significantly elevated in a sexually reproducing species, for either small or large divergence levels between species.

Oenothera species are diploid, therefore the evolution of asexuality involves loss of chromosome segregation and recombination during meiosis. The observed increase in ω can therefore be interpreted as fixation of heterozygous nonsynonymous changes captured in de novo transcript assemblies from asexual species. It is important to distinguish this from the familiar definition of fixation of homozygous mutations in sexually reproducing diploids. Homozygous fixation in asexuals would require either gene conversion or recurrent mutation between nonrecombining (nonsegregating) gene copes (Haag and Roze 2007).

Our phylogenetic analysis of ω is consistent with the patterns of within-species polymorphism (see above). This comparative analysis also suggests that the patterns of polymorphism in asexual species did not result simply from demographic effects, such as population bottlenecks associated with speciation. If bottlenecks drove our polymorphism results, we would not expect the observed patterns of ω, because it considers mutations in individual genomes.

It is important to point out that the heterogeneity in ω among sampled Oenothera clades was often greater than the difference between sexual/asexual sister taxa. This highlights that our data do not imply a simple relationship between the ω ratio, reproductive system and extinction risk in the genus, and clearly additional factors are influencing rates of molecular evolution. However, the consistent observation that shifts to asexuality are associated with increased genome-wide ω within a clade provides evidence for the theoretical prediction of accumulation of deleterious mutations in asexuals (Muller 1964; Felsenstein 1974; Keightley and Otto 2006). It also supports the corollary hypothesis that asexuality may increase the risk of extinction relative to closely related sexually reproducing species (Lynch et al. 1993).

The Age of Asexual Lineages

Asexual lineages are hypothesized to be evolutionarily short-lived, which lends support to the hypothesis that fitness costs of asexuality accrue over time. However, for the costs of asexuality to be apparent in Oenothera, asexual species must have persisted long enough for such costs to manifest themselves at the molecular level. For the four species pairs that had within-species sampling (table 1), we obtained an approximate age of asexual species by estimating the frequency of fixed differences in synonymous sites between sister sexual/asexual species and published mutation rates. Assuming a neutral mutation rate of 7.09 × 10⁻⁹ from Arabidopsis (Ossowski et al. 2010), and a generation time scaled to species’ life-history, point estimates for divergence time between sexual and asexual lineages are 506,477, 8,906, 5,534, and 12,014 years for clades 2 (O. speciosa/O. rosea; perennial/annual), 3 (O. filiformis/O. gaura; biennial), 4 (O. grandiflora/O. biennis; biennial), and 7 (O. grandis/O. laciniata; annual), respectively (figs. 1 and 2).

It must be emphasized that there is considerable uncertainty in estimates based on fixed differences. On the one hand, these are estimates of divergence between sister species, which is likely greater than the time since the evolution of asexuality. On the other hand, in the absence of recurrent mutations or gene conversion events between alleles, new mutations can no longer fix in asexual lineages due to permanent heterozygosity, which will lead to a downward bias in estimates of divergence (though mutations can still fix in sexual species). These estimates also assume a constant generation time for species pairs, an assumption known to be violated in O. speciosa (perennial)/O. rosea (annual).

Despite these caveats, our estimates do indicate that asexual lineages have likely evolved independently for thousands to hundreds of thousands of generations, enough time for the evolutionary consequences of asexual reproduction to become apparent (Lynch et al. 1993). It is important, however, to separate the consequences of asexuality in Oenothera from the evolutionary processes that give rise to asexual species, and we explore this issue below.

The Impact of Hybrid Origins on Sequence Variation in Asexual Oenothera

The evolution of asexuality by hybridization of divergent genomes is widespread in both plants and animals (Rieseberg 1997; Simon et al. 2003). A hybrid origin introduces initial allelic heterozygosity in asexual species. Additional heterozygosity accumulates in asexual lineages over evolutionary time, leading to high levels of divergence between homologous loci in ancient asexual lineages (Simon et al. 2003). Asexual Oenothera species are heterozygous for two haploid genome complexes that have accumulated reciprocal translocations within or between species (Cleland 1972). Because recombination between translocated chromosomes is suppressed, they diverge from one another over time. Among the functionally asexual Oenothera species we sampled, levels of heterozygosity were nearly twice as high as in sexual individuals. Some proportion of heterozygous sites may have arisen from allelic divergence after the evolution of asexuality, but some may have been inherited from historical divergence between parental genotypes (Levy and Levin 1975). Therefore, it is important to consider whether hybrid origins, rather than ongoing reductions in selection efficacy, could explain the elevated rate of nonsynonymous substitution and polymorphism in asexual Oenothera species.

To explore the effect of hybridization, we took advantage of the extensive prior knowledge of genome structure and hybrid origins in subsection Oenothera (fig. 1; clade 4). Species in this subsection are differentiated by three well-defined genomic complexes, designated A, B, and C (Stubbe 1964; Dietrich et al. 1997). Sexual species are homozygous for AA (O. elata, O. jamesii, O. longissima), BB (O. grandiflora), or CC (O. argillicola). Asexual species are either hybrids of two divergent genomic complexes (intercomplex hybrids; e.g., O. biennis, AB; O. parviflora, BC; O. oakesiana, AC) or are derived from a single ancestral complex that has accumulated structural heterozygosity through chromosomal translocations (intracomplex hybrids; e.g., O. wolfii and O. villosa, AA; O. nutans, BB) (Dietrich et al. 1997). Thus, we were able to explore the impact of hybrid origin, and divergence between founding genomes, on the proportion of nonsynonymous (relative to synonymous) mutations in the sample of transcriptomes in clade 4.

Among clade 4 asexual individuals, intercomplex hybrids harbored twice as much genomic heterozygosity (mean 0.006/site/individual) as intracomplex hybrids (mean 0.003/site/individual) at silent sites (P < 0.002). In contrast, intercomplex hybrids had a lower ratio of nonsynonymous to synonymous heterozygosity than intracomplex hybrids (mean 0.40 and 0.48, respectively, P < 0.04). Thus, asexual species derived from more divergent parental genomes have lower proportions of nonsynonymous variation compared with asexuals derived from more closely related parental genomes.

One limitation of the above analysis is that asexual species exhibit patterns of variation that reflect their genetic origins as well as subsequent sequence evolution. To infer the unique impact of historical divergence between parental genomes on the proportion of nonsynonymous variants in a hybrid (asexual) genome, we utilized polymorphism data from sexual Oenothera species in clade 4. We simulated “hybrids” between O. elata ssp. hookeri cv. Johansen Standard (JS), which has the putative ancestral genomic complex (Cleland 1972), and other sexual Oenothera species in clade 4, each of which varies in genomic divergence from JS (fig. 1; see Materials and Methods). In these simulated hybrid genomes, the ratio of nonsynonymous to synonymous heterozygosity was strongly negatively correlated with the homozygous divergence between “parental” genomes (r = −0.91, P < 0.001; supplementary fig. S3, Supplementary Material online). Taken together, these analyses indicate that hybridization per se does not increase the ratio of nonsynonymous to synonymous polymorphism, but rather decreases it, particularly in hybrid crosses of more divergent genotypes.

There is a straightforward explanation for why this should be the case. Most heterozygous sites in hybrid individuals are derived from high frequency polymorphisms, or fixed differences, that have accumulated between parental populations or species. Because purifying selection generally prevents fixation of nonsynonymous variants and constrains them to low frequencies, hybrid individuals generated from divergent parents are expected to initially have a reduced ratio of nonsynonymous to synonymous polymorphism at heterozygous sites. Hybrids derived from more closely related parents, conversely, combine more recent mutations (i.e., low-frequency polymorphisms), and therefore have higher proportions of nonsynonymous mutations (Mugal et al. 2014). Thus, to the extent that hybridization drives the relationship between polymorphism and θ_N/θ_S ratios among clades, we would expect the opposite pattern to what we actually observed (fig. 2b and supplementary fig. S1, Supplementary Material online).

In contrast, a time-dependent increase in the proportion of deleterious alleles within asexual species is predicted under Muller’s ratchet or a “mutational meltdown” model (Muller 1964; Lynch et al. 1993). We see such a time-dependent correlation among asexual Oenothera species (fig. 2 and supplementary fig. S1, Supplementary Material online), which suggests that the increased proportion of nonsynonymous variation is driven by suppressed recombination and segregation associated with asexuality, rather than by hybridization or other demographic influences. Genome-wide linkage in asexuals likely decreases the efficacy of purifying selection on deleterious mutations due to interference among mutations. However, a reduction in the strength of selection may also contribute, because asexual species are permanently heterozygous, which could mask the deleterious effects of partially recessive mutations.

Asexuality and Heterozygous Mutation Accumulation

We compared pN/pS ratios spanning four sexual-asexual transitions. These comparisons indicate increased accumulation of deleterious mutation in asexuals, compared with sister sexual species. Although we did not sample multiple individuals for all Oenothera species in our data set, we were able to compare the ratio of heterozygosity at nonsynonymous compared with synonymous sites (hN/hS) for all individuals and species. Because of a loss of recombination and segregation, diploid asexual species are predicted to accumulate deleterious heterozygous mutations over a large range of evolutionary parameters. To what extent this leads to decline in mean fitness in asexuals depends in important ways on the distribution of dominance coefficients, population size, and population structure (Haag and Roze 2007). In addition, as shown above, the impact of phylogenetic relationships and diversity at silent sites could also have strong effects on levels of nonsynonymous diversity (heterozygosity).

To evaluate evidence for accumulation of deleterious heterozygous mutations among all sampled Oenothera individuals, we fit a linear model testing the effects of clade, mode of reproduction (asexual vs. sexual), and silent site heterozygosity (hS) on nonsynonymous heterozygosity (hN). This model also took into account the interaction between clade and hS. Overall, the model explained greater than 99% of the variance in hN (adjusted R² = 0.9927; P << 10⁻¹⁰) across individuals. Taking into account the effects and interaction of hS and clade, sexual reproduction was significantly negatively correlated with hN (P < 0.00005). Thus, patterns of heterozygosity across the full data set further support the theoretical prediction that asexuality is associated with accumulation of deleterious mutations.

Conclusions

Theoretical work predicts that interference among linked mutations may select for increased meiotic recombination and segregation associated with sexual reproduction (Felsenstein 1974; Barton and Otto 2005; Keightley and Otto 2006; Agrawal 2009). These studies further indicate that genome-wide linkage in asexual populations can cause reduced effective population size, allowing deleterious mutations to rise to high frequencies and fix within populations (Muller 1964; Charlesworth 2012). Our analyses confirm the prediction of accumulation of deleterious mutations following repeated loss of sex in the genus Oenothera. These results were robust and consistent at both microevolutionary and macroevolutionary time-scales. Given these patterns, repeated loss of sex in this genus likely results from short-term benefits of reproductive assurance (Baker 1955; Barrett 2002), whereas the ephemerality of functionally asexual reproduction (Johnson et al. 2011) may be associated with the long-term costs of accumulating deleterious mutations. This work provides strong genome-wide support for the link between asexuality and mutation accumulation (Muller 1964; Lynch et al. 1993), providing a plausible mechanism underlying the rarity and ephemerality of asexual eukaryotic lineages due to increased extinction risk.

Materials and Methods

Plant Materials

All plant material used in the experiment originated from seeds taken from either natural populations or previous experiments. We sampled a total of 62 individual plants from 32 species and subspecies (supplementary table S1, Supplementary Material online). These species were sampled to maximize phylogenetic diversity across the Oenothera genus, where the functionally asexual PTH genetic system is most prevalent (Raven 1979; Wagner et al. 2007; Johnson et al. 2009). Our sampling included 32 individuals known to have a PTH genetic system (referred to as “asexual species” hereafter) and 30 individuals that are non-PTH (referred to as “sexual species” hereafter). This sampling represented a minimum of ten independent transitions between sexual and functionally asexual reproduction (see fig. 1). To gain insight into population genomic patterns of evolution in sexual and asexual species, we sampled one individual from each of 3–5 populations in each of eight species. These eight species represented four pairs of closely related sexual and asexual species. The species’ names, genetic system, and the origin of samples are included in supplementary table S1, Supplementary Material online. Seed material and vouchers for all samples studied here are curated by M.T.J.J.

Germination, Plant Growth, and Tissue Harvesting

Seeds were initially soaked at 4 °C for 24 h in 1 ml of 0.1% agarose solution in filtered water. All species except those with indehiscent fruits (i.e., O. gaura, O. filiformis, O. suffulta) were then germinated for 1–3 weeks in 55-mm Petri dishes on wet filter paper left in a windowsill exposed to full sun for 3–4 h per day. Plant species with indehiscent fruits were planted directly 1 cm below the soil surface and placed under fluorescent lights at room temperature until seedlings emerged. Upon germination, seedlings were transplanted individually into plastic plug tray wells containing approximately 150 ml of steam-sterilized potting soil (Ready Earth Seedling and Plug Soil Mixture, Sungro Horticulture, Bellevue, WA). We grew all plants in a single reach-in growth chamber set to a 16 h (25 °C):8 h (20 °C) light:dark cycle. Lighting was provided by fluorescent and incandescent bulbs at an intensity of 450 μmol m⁻² s⁻¹. Plants were watered ad libitum with deionized water containing a one-third strength standard Hoagland’s solution that provided macro- and micronutrients. We harvested the fourth true leaf from each plant, which we excised when it reached one-quarter of the size of the first true leaf. This criterion insured that all plants were harvested at the same developmental stage, and at a time when the leaf was undergoing active growth and development. Leaves were flash-frozen and stored at −80 °C until extraction of total RNA. For details on Oenothera cultivation, see Greiner and Köhl (2014).

RNA Extraction

Total RNA was extracted using a hybrid CTAB/acid phenol/silica membrane method developed specifically for this project. The full details of our extraction protocol are available in Johnson et al. (2012) (see protocol 9), and we provide a brief summary here. We extracted total RNA from 50 to 100 mg of tissue per sample. Tissue was first ground to a powder with a prechilled mortar and pestle sterilized with RNase Zap (Ambion, Austin, TX). CTAB extraction buffer (2–5 ml) was added and mixed with the sample in the mortar. We then transferred the sample to a 15-ml tube, vortexed the sample (15 s), and incubated it at 65 °C (15 min). The sample was centrifuged at 14,000 × g (3 min) and the supernatant was transferred to a new 2-ml tube. An equal volume of 24:1 chloroform:isoamyl alcohol was added to each sample, vortexed, and then centrifuged at 14,000 × g (3 min). The upper aqueous phase was transferred to a new 2-ml tube and the solvent extraction was repeated. After transferring the upper aqueous phase to a new tube we added an equal volume of 5:1 phenol:chloroform (pH 4.5), followed by vortexing and centrifugation at 14,000 × g (3 min). The upper aqueous phase was moved to a new 2-ml tube and the phenol:chloroform solvent extraction was repeated multiple times until we saw a clear boundary layer between the upper and lower aqueous solutions (usually 2–3 extractions). After transferring each sample to a new tube we added an equal volume of 24:1 chloroform:isoamyl alcohol, vortexed the sample, and then centrifuged samples at 14,000 × g (3 min). The upper aqueous phase was transferred to a new tube and the extraction was repeated once. After estimating the volume of each sample, 0.5 volumes of RLT solution from the Qiagen RNeasy Plant minikit (Qiagen, Valencia, CA) was added to each sample and thoroughly mixed by inversion. We precipitated RNA by adding 0.5 volumes of 95% EtOH followed by inverting tubes. Total RNA was bound to a silica membrane by transferring the entire sample to a pink Qiagen miniRNA spin column followed by centrifugation at 8,000 × g (15 s). We added 350 μl of RW1 buffer wash to the spin column and centrifuged the column at 8,000 × g (15 s). DNase solution (80 μl) was applied to the membrane and incubated (15 min), followed by repeating the RW1 wash step. The column was then washed with 500 μl of RPE and centrifuged at 6,300 × g (15 s); this step was repeated once. After additional centrifugation to remove all remaining liquid from the column we added 44 μl of RNase-free water onto each column. Samples were incubated for 1 min to allow RNA to fully dissolve before centrifuging samples at 6,300 × g into a RNase free tube. The purity of total RNA was assessed on a NanoDrop 8000 spectrophotometer (Thermo Scientific, Wilmington, DE). The concentration, integrity, and yield of total RNA were determined by assaying 1 μl of diluted total RNA on an Agilent 2100 Bioanalyzer using the Plant RNA Nano or Pico chip following the manufacturer’s instructions (Agilent Technologies, Sanata Clara, CA).

Transcriptome Library Preparation and Sequencing

Library preparation and sequencing were performed by BGI-Shenzhen (Shenzhen, China) as part of the One Thousand Plants Consortium project (1KP, www.onekp.com), as well as at North Carolina State University’s (NCSU) Genomic Sciences Lab (see supplementary table S1, Supplementary Material online). The methods for library preparation and sequencing at BGI were previously published and are briefly described below.

For the BGI samples, we isolated polyA mRNA from 20 mg of total RNA using the Dynabeads mRNA Purification Kit (Life Technologies, Grand Island, NY). mRNA was fragmented to a size of 200–300 nt using the RNA Fragmentation Reagents kit (Life Technologies). The first and second strands of cDNA were synthesized using reverse transcriptase and DNA polymerase, respectively. Double-stranded DNA strands were then isolated using a QIAquick PCR Purification Kit (Qiagen), followed by end-repair and the addition of a 3′ adenosine (A base). The Illumina PE Adapters were ligated onto the A base and polymerase chain reaction (PCR) amplified for 15 cycles using “indexed” paired-end PCR primers. Successful library preparation was confirmed by gel electrophoresis and the transcriptome library was further size-selected by gel-purification to ensure that insert size was 200 bp and total fragment size was 322 bp. The NCSU samples were prepared in a similar manner except that we used the Illumina TruSeq RNA protocol to isolate mRNA and we only attached single-read adapters.

All paired-end sequencing was performed on the Illumina HiSeq platform. We sequenced samples as paired-ends at a read length of 90 bp with up to 11 samples multiplexed into each lane of the flow cell. Single-read sequences at NCSU were performed on the Illumina GAIIx platform at a read length of 72 bp with one sample per lane. Prior to assembly, sequences were filtered to minimize contamination and maximize sequence quality. All sequences met the following four criteria: 1) Contained a unique indexed barcode, 2) less than 15 bp of continuous adapter sequence, 3) less than 0.5% uncalled bases on either paired read, and 4) ≥80% of bases with a Q-score of greater than 10 in both paired reads.

De Novo Assembly of Transcripts

Illumina reads from each individual were assembled using Velvet Oases (Schulz et al. 2012). VelvetOptimiser version 2.2.4 (Zerbino and Birney 2008) was invoked to choose the best kmer size (between 27 and 61 nt) for each individual transcript. To avoid missing low coverage transcripts, the final total number of bases in each assembly was used to evaluate the best kmer size. Oases version 0.2.08 was run under default parameters. For each set of transcript isoforms, the longest was chosen as the final transcript.

Read-Mapping and Genotyping

To genotype individuals, a single reference transcriptome was chosen from among the individuals in each clade (1–7) on the basis of N50 contig length. Reads from all individuals in each clade were then mapped using Stampy version 1.0.21 (Lunter and Goodson 2011). Postmapping alignment processing was performed using Picard Tools version 1.4.8 (http://picard.sourceforge.net). Regions surrounding insertion/deletion polymorphisms were realigned, and genotyping was performed using the Genome Analysis Tool-Kit version 2.39 (McKenna et al. 2010). For all read-mapping and genotyping steps, default program settings were used.

Base-Call Filtering

After genotyping, sites were filtered for genotype quality (genotype likelihood ratios ≥60 for all individuals) and read coverage (≥10× for all individuals). We implemented an additional filter following Gayral et al. (2013). This filter is based on a likelihood ratio test of a single-locus versus two-locus model, where the former represents the likelihood that reads from the same orthologous locus in all individuals map to a single reference locus, and the latter represents the likelihood that reads from multiple loci map to a single reference locus. The two-locus scenario might take place when, for example, reads from two closely related duplicate genes are assembled as a single locus in the reference transcriptome. This would result in excess heterozygosity (compared with Hardy–Weinberg expectations) for a given allele frequency, and uneven representation of segregating bases among reads mapping to the site. We filtered sites that were a significantly better fit to the two-locus model, and filtered entire contigs for which ≥50% of SNPs were a significantly better fit. Models were fit using data from the outcrossing samples only, because asexual species are not expected to fit Hardy–Weinberg expectations even in the absence of paralogous mapping. The majority of sites filtered using this method were fixed heterozygotes in both sexual and asexual samples. The results presented above were qualitatively equivalent considering data filtered only by sequencing depth and genotype quality scores.

Genomic Analysis

Open reading frames (ORFs) were identified using the getorf program from the EMBOSS suite version 6.3.1. For all loci in all individuals, a custom PERL script was used to identify synonymous and nonsynonymous variants from the read-mapped data (available upon request). Per-site estimates of the number of synonymous and nonsynonymous polymorphisms were calculated following Li (1993). For comparison of within-species polymorphism, we restricted analysis to subsets of within-species samples that were monophyletic with respect to sister taxa (see fig. 1 and table 1). Statistical analysis was performed using R version 3.0. Watterson’s θ was estimated following Tajima (1989) using a custom R function. We estimated the distribution of fitness effects of nonsynonymous variants using the method of Eyre-Walker and Keightley (2009), implemented in a PERL script kindly provided by Dan Halligan. For this analysis, we used 4-fold degenerate positions as a proxy for neutrally evolving sites.

Ortholog Identification

We identified orthologous genes within each clade separately by all-against-all BLASTn searches among individuals within a clade. We subsequently identified orthologs across clades by all-against-all BLASTn searches using reference assemblies from each clade. We filtered loci by BLAST e value ≤ 10⁻¹⁰ and reciprocal best hits over ≥80% of sequence length. In total, we assigned orthology for 4,081 genes.

Sequence Alignment

Orthologous sequences were compiled into single fasta-formatted files. All alignments were done using the program MUSCLE version 3.8.31 (Edgar 2004). For dN/dS analysis, ORFs were aligned for each ortholog set. ORF alignments were then used to guide in-frame alignments of nucleotide sequences.

Phylogenetic Analysis

Maximum-likelihood phylogenetic trees were estimated using RaXML version 7.0.4 (Stamatakis 2006), using default settings (including substitution model = GTRGAMMA). For the phylogeny presented in figure 1, we identified orthologous transcripts by all-against-all BLASTp of the Oenothera ortholog set against a de novo assembly of Chamerion angustifolium (produced using the same methods as above). We identified 939 orthologous outgroup loci. Orthologs and outgroup sequences were aligned as above. Alignments were concatenated into a single large alignment 1,976,729 nucleotides in length. The phylogeny for the sample was estimated based on the concatenated sequence. For each ORF-guided ortholog alignment, a separate maximum-likelihood tree was also estimated.

ω (dN/dS) Analysis

Individual gene trees were used as input for ω analysis in PAML version 3.4 (http://abacus.gene.ucl.ac.uk/software/paml.html; last accessed December 26, 2014). For each of the 4,081 orthologs, we fit a “free-ratio” model, in which a separate ω value was estimated for every branch on the phylogeny. For the free-ratio model, branch-specific ω values from the PAML output file were extracted and reformatted to conform to Newick format and read into R using the “phytools” package (Revell 2012). Terminal branch values were then extracted in R, and values across loci recorded for each individual. We restricted the analysis to loci with ω ratios ≤2, to examine patterns among genes under selective constraint.

Simulation of Hybrid Genomes

To simulate the initial effect of hybridization on the ratio of nonsynonymous to synonymous heterozygous mutations, we simulated hybridization between O. elata (JS) and all other sexual clade 4 Oenothera species (fig. 1). The proportion of genotypic differences at nonsynonymous and synonymous sites was calculated, simulating the proportion of nonsynonymous and synonymous heterozygous sites (hN/hS) that would initially be present in a hybrid between O. elata (JS) and other (clade 4) sexual species with varying levels of divergence. The hN/hS ratio was evaluated as a function of hS across all simulated hybrids (supplementary fig. S3, Supplementary Material online).

Acknowledgments

The authors thank H. Myburg for performing RNA extractions, the NC State Phytotron staff for providing technical and infrastructural support, L. Yaneva-Roder and E. Carpenter for providing technical and logistical support, and M. Deyholos for outgroup data. A.A. Agrawal, K. Krakos, R. Raguso, and the Ornamental Plant Germplasm Centre provided seed material. This work was supported by an NSERC Discovery Grant, NSF (DEB-0919869), the Early Researcher Award (Ontario Government), the Ontario Research Fund, the Canadian Foundation for Innovation to M.T.J.J.; Alberta Enterprise and Advanced Education, Genome Alberta, Alberta Innovates Technology Futures iCORE, Musea Ventures, and BGI-Shenzhen to G.K.S.W.; the Max Planck Society and the Deutsche Forschungsgemeinschaft to S.G.; and an NSERC Discovery Grant and an NSERC Accelerator Award to S.W.

References

Author notes