1. Evolutionary Biology
  2. Genetics and Genomics
Download icon

The effect of hybridization on transposable element accumulation in an undomesticated fungal species

  1. Mathieu Hénault  Is a corresponding author
  2. Souhir Marsit
  3. Guillaume Charron
  4. Christian R Landry
  1. Institut de Biologie Intégrative et des Systèmes (IBIS), Université Laval, Canada
  2. Département de biochimie, microbiologie et bioinformatique, Université Laval, Canada
  3. Quebec Network for Research on Protein Function, Engineering, and Applications (PROTEO), Université Laval, Canada
  4. Université Laval Big Data Research Center (BDRC_UL), Canada
  5. Département de biologie, Université Laval, Canada
Research Article
  • Cited 0
  • Views 956
  • Annotations
Cite this article as: eLife 2020;9:e60474 doi: 10.7554/eLife.60474

Abstract

Transposable elements (TEs) are mobile genetic elements that can profoundly impact the evolution of genomes and species. A long-standing hypothesis suggests that hybridization could deregulate TEs and trigger their accumulation, although it received mixed support from studies mostly in plants and animals. Here, we tested this hypothesis in fungi using incipient species of the undomesticated yeast Saccharomyces paradoxus. Population genomic data revealed no signature of higher transposition in natural hybrids. As we could not rule out the elimination of past transposition increase signatures by natural selection, we performed a laboratory evolution experiment on a panel of artificial hybrids to measure TE accumulation in the near absence of selection. Changes in TE copy numbers were not predicted by the level of evolutionary divergence between the parents of a hybrid genotype. Rather, they were highly dependent on the individual hybrid genotypes, showing that strong genotype-specific deterministic factors govern TE accumulation in yeast hybrids.

eLife digest

Hybrids arise when two populations of organisms that are related but genetically different mate and produce offspring. Hybridization has long been regarded as one of the many ways species evolve. Studying the changes in the genome that result from this process can provide insights into evolutionary history and predict the outcome of mixing between genetically different populations. In fact, the inability of two organisms to mate and produce viable and fertile hybrids has been used as a way to define species. It has been speculated that the infertility of many hybrids is due to short sequences of DNA in the genome called transposable elements. These elements are sequences of DNA that, when active, can move to a different position in the genome, causing mutations. It is thought that the process of hybridization may be activating transposable elements leading to the infertility often observed in hybrids.

The activation of transposable elements in hybrids has been studied in animals and plants, and usually, the hybrids studied were either generated in the laboratory or found in the wild. Fungal species, such as the yeast Saccharomyces paradoxus, have hundreds of wild strains, including many hybrids, and can also be crossed in the laboratory to produce new hybrids, allowing a combined approach to studying the activation of transposable elements. Hénault et al. used this yeast to investigate whether hybridization leads to increased activity of transposable elements in fungi.

To test this hypothesis, Hénault et al. analyzed the genomes of hundreds of natural strains of S. paradoxus to count and locate their transposable elements and establish evolutionary relationships between them. Next, they crossed different strains in the laboratory to see how the transposable elements would act upon hybridization.

If transposable elements were activated by hybridization, then hybrids would accumulate more transposable elements. However, the analyses did not show increased numbers of transposable elements in wild hybrids of S. paradoxus. This could be explained by an actual absence of transposable element activation, or by natural selection eliminating individuals that accumulate more transposable elements. To determine which is the case, Hénault et al. next recreated several hybrids in the laboratory and reproduced them for hundreds of generations. Hybrids were grown in the laboratory such that natural selection was almost incapable of favoring some yeasts over others, allowing the hybrids to accumulate transposable elements. These experiments revealed that hybrids accumulated transposable elements at different and largely unpredictable rates. Indeed, closely related hybrids often had highly different numbers of transposable elements in their genomes after being reproduced in the laboratory.

These observations indicate that the accumulation of transposable elements depends on various factors and cannot be easily predicted, and that hybridization may only be a small piece of the puzzle. Additionally, Hénault et al. demonstrated that undomesticated organisms like fungi can provide unique insights into evolutionary hypotheses.

Introduction

Hybridization is increasingly recognized as an important component of species evolution (Mallet, 2005; Mallet et al., 2016; Stukenbrock, 2016). Abundant empirical evidence shows that hybridization frequently leads to transgressive phenotypes in plants, animals, and fungi (Langdon et al., 2019; Rieseberg et al., 1999). Mutation rates are traits for which transgressive segregation can have a major impact on the evolutionary fate of hybrids, with potential implications for the emergence of reproductive barriers and speciation. Insertions caused by transposable elements (TEs) are a class of mutations that could be particularly responsive to hybridization. TEs are diverse and ubiquitous DNA sequences that can duplicate and spread within genomes, forming families of dispersed repeated sequences (Bourque et al., 2018). They can either replicate through ‘copy-and-paste’ (retrotransposons) or ‘cut-and-paste’ (DNA transposons) mechanisms. TEs were first identified in cytological studies of chromosome instability in Zea mays (McClintock, 1950). In her seminal 1984 perspective, McClintock, 1984 summarized evidence suggesting that TEs respond to various types of ‘genome shocks’, one of which she suspected was hybridization. Since McClintock’s early insights, many studies investigated the association between TEs and hybridization across taxonomic groups and documented increased TE activity in hybrids. We refer to this increase as TE reactivation, but our use of this term implies no mechanistic basis and other terms like derepression are used in the literature. Reactivation was evidenced by increased TE abundance, novel TE insertions, TE epigenetic changes or higher TE transcript levels (Dion-Côté et al., 2014; Kidwell et al., 1977; Labrador et al., 1999; Liu and Wendel, 2000; O'Neill et al., 1998; Picard, 1976; Ungerer et al., 2006). However, hybrid reactivation is not universal, as many studies across a similar taxonomic range found no evidence for it (Coyne, 1989; Göbel et al., 2018; Hey, 1988; Kawakami et al., 2011). In fungi, a recent preprint showed that no TE reactivation occurs in experimental interspecific hybrids between the yeasts Saccharomyces cerevisiae and Saccharomyces uvarum (Smukowski Heil et al., 2020). Thus, TE reactivation exhibits no clear taxonomic distribution, and the conditions favorable for its manifestation are still unclear.

Studies of hybrids in controlled crosses revealed genetic bases to TE reactivation, but less focus was put on cases of natural hybrids (but see Kawakami et al., 2011; Ungerer et al., 2009). In addition to the genetic structure of populations, environmental pressures and demographic factors can influence the action of evolutionary forces on TEs, thus impacting the genomic landscapes of TEs in natural populations, including hybrids. Population genetic models have investigated the action of the major evolutionary forces on mean TE copy number (CN) per individual (Charlesworth and Charlesworth, 1983; Charlesworth and Langley, 1989). Three such forces are deterministic (Charlesworth et al., 1994). First, TE CN increases at a certain transposition rate, which incorporates the various forms of transposition regulation (either by the host or by TEs themselves). Second, TEs are removed from the genome at a certain excision rate. Third, most TE insertions are deleterious, such that individual fitness is assumed to be a decreasing function of CN and natural selection acts to purge TE insertions from populations. The fourth important force is genetic drift, a stochastic force that impacts the efficiency of natural selection by increasing the probability of fixation of deleterious mutations (Lynch, 2007). The intensity of drift is inversely proportional to the effective population size Ne. In populations with small Ne, drift renders natural selection inefficient, such that deleterious mutations are more likely to fix by chance. In principle, any of those evolutionary forces could affect hybrids in a way that differs from their parents, leading to hybrid-specific dynamics in the accumulation of TEs.

Population data in Drosophila indicate that deterministic forces have a predominant influence on the allele frequency spectrum of TE insertions (Charlesworth and Langley, 1989), consistent with a minimal contribution from drift. Among the deterministic forces, transposition rate is arguably the most likely to be impacted by hybridization, for instance if genetic incompatibilities impair transposition regulation mechanisms. A classic illustration of this is hybrid dysgenesis in Drosophila, a syndrome caused by mismatches between paternal TEs and egg-loaded pools of Piwi-interacting RNAs (piRNAs) that lead to inefficient TE repression (Brennecke et al., 2008; Bucheton et al., 1984; Erwin et al., 2015; Rubin et al., 1982). Similarly, hybrids between Arabidopsis thaliana and Arabidopsis arenosa show overexpression of a TE family from the A. arenosa subgenome that is dependent on the ploidy level of the parental accessions involved, suggesting that dosage imbalance of epigenetic regulators causes reactivation (Josefsson et al., 2006). It can also be expected that the number of genetic incompatibilities increases with genetic divergence, yielding a positive relationship between evolutionary divergence and TE deregulation (Orr, 1995; Parhad and Theurkauf, 2019). Additionally, evidence from Drosophila suggests that hybridization may impact excision rates for specific TE families (Coyne, 1989). Along with these deterministic factors, reductions in natural selection efficiency should also be considered. Evidence from Arabidopsis (Lockton et al., 2008; Lockton and Gaut, 2010), Caenorhabditis (Dolgin et al., 2008) and Daphnia (Ye et al., 2017) suggests that long-term differences in Ne caused by demography or mating system are associated with variation in TE content (Lynch and Conery, 2003). A reduction in Ne can be expected when a novel hybrid lineage is formed. Indeed, it has been hypothesized that population bottlenecks inherent with hybrid speciation in Helianthus might have played a role in the accumulation of TEs through a decreased efficiency of purifying selection (Ungerer et al., 2009).

An emerging eukaryotic model system in ecological and evolutionary genomics is the yeast Saccharomyces paradoxus, the undomesticated sister species of the budding yeast S. cerevisiae. S. paradoxus has a complex admixture history in North America, making it an ideal model system to test hybridization-related evolutionary hypotheses in natural populations and in the laboratory (Figure 1; Hénault et al., 2017). North American S. paradoxus lineages comprise the Eurasian lineage SpA, the endemic incipient species SpB and SpC (hereafter termed pure lineages), and the admixed lineages SpC*, SpD1 and SpD2 (hereafter termed hybrid lineages; Figure 2a). SpC* is an incipient hybrid species that arose from a post-glacial secondary contact between the endemic lineages SpB and SpC and consists of mostly SpC-like genomes with ~5% introgressions from SpB (Leducq et al., 2016), while SpD1 and SpD2 arose more recently from hybridization between SpC* and SpB and have roughly equal genomic contributions from both parents (Eberlein et al., 2019; Figure 2b).

Population genomic data and laboratory evolution experiments allow to uncover the effects of various factors on transposable element (TE) accumulation in hybrids.

Upon hybridization, the TE content of a hybrid is determined by equal contributions from the TE contents of its parents. Subsequent evolution can either leave the hybrid TE content unchanged, or drive it to increase or decrease. These changes can be driven by biases in transposition or excision rates or by variation in natural selection efficiency. Two complementary approaches can be used to investigate which factors drive TE dynamics in hybrids. First, genomic data from natural populations can be harnessed to understand how population structure, environmental pressures and natural selection efficiency shaped natural variation in TE content. Second, evolution experiments in laboratory controlled conditions (constant environment and relaxed natural selection) allow to test the effect of properties of artificial hybrid genotypes on TE accumulation, namely evolutionary divergence between parents and initial TE content.

Figure 2 with 3 supplements see all
Ty elements showed no reactivation in hybrid lineages and contrasted abundances among pure lineages.

(a) Bayesian phylogenetic tree of the North American Saccharomyces paradoxus lineages based on genome-wide SNPs. (b) Schematic admixture graph showing the evolutionary relationships between pure and hybrid lineages as inferred by Eberlein et al., 2019; Leducq et al., 2016. Dotted lines represent admixture events. (c) Counts of full-length elements (top) and solo LTRs (bottom) in six whole-genome assemblies based on long reads. No full-length Ty4 or Ty5 were found in any assembly, thus these families were considered extinct. (d) Copy number (CN) variation in active Ty families (i.e. families which have full-length elements) measured as log2 normalized read depth (NRD) over Ty reference internal sequences for 207 diploid wild strains. Sample size for each lineage (excluding strains for which there was no coverage) is indicated between parentheses. Whiskers span 1.5 times the interquartile range. The p-values of Kruskal-Wallis tests are shown. Stars on bars show the results of Conover post-hoc tests. *: p<0.05, **: p<0.01, ***: p<0.001.

Yeasts of the Saccharomyces genus are valuable model species to investigate the interplay between hybridization and transposition. Saccharomyces genomes contain many long terminal repeats (LTR) retrotransposons of the Ty1/copia (Ty1, Ty2, Ty4 and Ty5) and Ty3/gypsy (Ty3) families. Ty1/copia and Ty3/gypsy retrotransposon superfamilies (Wicker et al., 2007) respectively correspond to the Pseudoviridae and Metaviridae virus families (Krupovic et al., 2018) and have a broad taxonomic distribution. Ty elements comprise two overlapping open reading frames (ORFs) within a ~5 kb internal sequence, which is flanked by two ~300 bp LTRs in direct orientation. The two LTRs of full-length elements are prone to intra-element recombination, leaving behind solo LTRs in the genome (Parket et al., 1995). Ty1 in the budding yeast S. cerevisiae is among the best characterized eukaryotic TEs, benefiting from decades of research on its molecular and cellular biology and, more recently, on its population genomics (Bleykasten-Grosshans et al., 2013; Curcio et al., 2015). Although Saccharomyces genomes harbor a relatively modest fraction of Ty elements (~1–3% in size), their streamlined genomes (~12 Mbp) and short generation time (~2 hr) provide great power to perform TE comparative genomics and laboratory evolution experiments.

Here, we use population genomic data and laboratory evolution experiments on S. paradoxus and its sibling species S. cerevisiae to investigate the hybrid reactivation hypothesis and, more generally, the factors governing TE accumulation in natural and experimental lineages. Our results show that natural and experimental hybrids exhibit no reactivation of TEs. Rather, we show that deterministic factors like population structure and the properties of individual hybrid genotypes are major determinants of TE content evolution in hybrid genomes.

Results

S. paradoxus lineages exhibit contrasted TE evolutionary dynamics but no reactivation in hybrids

To investigate the TE content and dynamics in North American S. paradoxus pure and hybrid lineages, we annotated six highly contiguous whole-genome assemblies previously generated from Oxford Nanopore long reads (PRJNA514804) (Eberlein et al., 2019). We used a database comprising reference sequences for S. cerevisiae Ty1, Ty2, Ty3, Ty4, and Ty5 families, along with S. paradoxus Ty3-related Ty3p (Carr et al., 2012) and Ty4-related Tsu4, the latter originating from a horizontal transfer from the Saccharomyces eubayanus/Saccharomyces uvarum clade to S. paradoxus (Bergman, 2018). These annotations revealed counts of full-length elements an order of magnitude lower than solo LTRs (Figure 2c, Figure 2—figure supplement 1), consistent with previously reported data on S. cerevisiae and S. paradoxus (Carr et al., 2012; Goffeau et al., 1996; Kim et al., 1998; Yue et al., 2017). To overcome the difficulty of comparing multiple families of short and often degenerated LTR sequences within a phylogenetic framework, we built intra-genome LTR sequence similarity networks. Sequences within the Ty1-Ty2 and Ty3p-Ty3 LTR family pairs were indistinguishable (Figure 2—figure supplement 2), in agreement with Ty2 and Ty3 being specific to S. cerevisiae (Carr et al., 2012; Naumov et al., 1992). Consequently, Ty2 and Ty3 LTR annotations were merged with Ty1 and Ty3p, respectively. Although they were clearly distinct, Ty4 and Tsu4 LTRs showed a close relatedness, which justified their combination into the Ty4/Tsu4 family. No Tsu4 LTR was annotated in the SpA genome, consistent with this family being specific to North American S. paradoxus and absent from Eurasian S. paradoxus (Bergman, 2018).

Based on the assembly annotations, Ty1 was the most abundant family and exhibited the largest CN variation among lineages, with a striking difference between SpB and SpC (Figure 2c). The generality of CN patterns revealed by genome assemblies was confirmed by normalized short-read depth (NRD) over Ty reference sequences for 207 natural diploid isolates (Figure 2d, Figure 2—figure supplement 3, Supplementary file 1a; PRJNA277692, PRJNA324830, PRJNA479851). The hybrid lineages SpC*, SpD1 and SpD2 had intermediate full-length Ty1 and Tsu4 CNs compared to their parental lineages, indicating no apparent persistent reactivation or maintenance of Ty elements in these lineages. Notably, Ty3p appeared extinct in all hybrid lineages, with no signal for internal sequences. The most outstanding profile was that of the pure lineage SpC, which had significantly higher full-length Ty1 and Ty3p CNs compared to the other lineages. For Tsu4, SpB was the lineage with the highest CN. Thus, lineage-wide CNs were consistent with an absence of reactivation in hybrids, the most contrasted profiles being found among pure lineages.

CNs alone may not reflect the actual dynamics of a TE family within a genome. For instance, a recently expanding family may not have had the time to reach the CN of another family (or the same family in a different lineage) that would be more abundant but at equilibrium. We thus explored the evolutionary dynamics of Ty families within each genome. If high CNs were caused by recent family expansions, they should be associated with a signature of low divergence (Bergman and Bensasson, 2007). We used the sequence divergence to the most closely related LTR for each LTR annotation to approximate evolutionary divergence since transposition. The importance of low-divergence peaks (≤1%) for Ty1 was significantly higher in SpC and SpA compared to SpB (Fisher’s exact test, FDR corrected, Figure 3a; see Figure 3—figure supplement 1 for other divergence thresholds). SpA also had a significant underrepresentation of low-divergence LTRs for the composite Ty4/Tsu4 family due to the absence of Tsu4. The height of low-divergence peaks was consistent with CNs of full-length elements (Figure 2d). Hybrid lineages SpC*, SpD1 and SpD2 had Ty1 peaks similar to SpC, although less important. They also had no peak for Ty3p, consistent with the absence of full-length elements. Overall, this result indicated that CN variation reflected variation in family dynamics and confirmed that Ty elements exhibited no reactivation in hybrid lineages.

Figure 3 with 3 supplements see all
Evolutionary dynamics of Ty families in S. paradoxus.

(a) Distributions of minimum nucleotide divergence between long terminal repeat (LTR) sequences. Heatmaps show FDR-corrected p-values for pairwise Fisher’s exact tests between ratios of elements ≤1% and >1% nucleotide divergence. Color maps are centered at the significance threshold of 0.05. Blue cells depict statistically significant tests. See Figure 3—figure supplement 1 for other divergence thresholds. (b) LTR divergence and conservation levels. Horizontal positions depict Ty orthogroups clustered from conserved (left) to private (right). Color map shows minimum nucleotide divergence. Scattered red dots show LTRs that belong to full-length elements. FDR-corrected p-values of Fisher’s exact tests for the ratio of full-length elements in private orthogroups are shown. ns: p≥0.05, *: p<0.05, **: p<0.01, ***: p<0.001. (c) Lineages through time curves for Ty1 ultrametric bayesian phylogenetic trees, using a conserved Ty4 LTR orthogroup as an outgroup. Branching time is scaled such that values of one correspond to the present. From left to right, distributions show all LTRs (with a close-up on recent dynamics), LTRs from private elements and LTRs from full-length elements.

Next, we tested whether orthology between genome-wide Ty insertion loci was consistent with transposition activity. We defined Ty orthogroups from a whole-genome multiple alignment of the six assemblies (Figure 3b). Most orthogroups were either private to one genome or conserved across the six (Figure 3—figure supplement 2), suggesting that a major fraction of Ty insertion polymorphisms segregated at low frequency in North American lineages. This was consistent with studies in S. cerevisiae showing that insertions are mostly fixed or at low frequency (Bleykasten-Grosshans et al., 2013; Carr et al., 2012). Private orthogroups comprised LTRs that had low divergence, consistent with recent transposition activity (Figure 3—figure supplement 2). Private orthogroups were significantly enriched in full-length Ty1 elements (Figure 3b), including in hybrid lineages. This indicated that all lineages had active transposition for at least some Ty families. It also suggested that an absence of transposition in hybrid lineages could not explain the absence of reactivation. We called lineage-wide insertion polymorphisms from discordant mappings of the 207 short read libraries to estimate the distribution of insertion allele frequencies. This calling method is noisy and only allows to call non-reference alleles, but should nevertheless provide rough estimates of the insertion allele frequency spectra. Principal component analysis (PCA) on non-reference insertions yielded the expected grouping of strains into well-defined lineages (Figure 3—figure supplement 3). The allele frequency distributions showed a similar excess of very-low-frequency alleles for the six lineages (Figure 3—figure supplement 3), suggesting that Ty elements were subject to largely uniform evolutionary forces across lineages.

We focused on the most abundant Ty1 to investigate the transposition dynamics in a phylogenetic framework (Le Rouzic et al., 2013). Bayesian phylogenetic trees of Ty1 LTRs were converted into ultrametric phylogenies to approximate the timing of branching events using a conserved Ty4 LTR orthogroup as an outgroup. Branching times revealed that lineages had largely similar Ty1 dynamics (Figure 3c). Due to the shortness of LTR sequences, most branches collapsed with the outgroup, indicating an expected saturation of the molecular substitution signal. A focus on recent dynamics revealed a clear expansion signature in SpC that was shared by SpC* and, to a lesser extent, by SpD1 and SpD2 (Figure 3c), which was reflected by the peaks of young private and full-length elements. This suggests that any recent expansion signature in hybrid lineages can be largely explained by inheritance from their SpC-related parent. The rapid expansion in SpC and SpC* can be contrasted with the more gradual expansion in SpA, which harbored the second highest Ty1 CN. This contrasted with the assumption that TE abundances can be interpreted as equilibrium CNs under the transposition-selection balance hypothesis (Charlesworth and Charlesworth, 1983) and rather suggested a recent transposition rate increase in SpC-related lineages.

Population structure best explains TE content variation, with minor contributions from climatic variation

Our results showed that hybrid lineages exhibited no outstanding Ty CN profiles or sequence dynamics. However, potential hybrid-specific differences in transposition rate might be masked by other factors that can be expected to impact Ty dynamics. To exclude this possibility, we investigated which factors best explain the variation in Ty CN in natural lineages. First, population structure incorporates both host-encoded (Curcio et al., 2015) and Ty-encoded (Czaja et al., 2020; Saha et al., 2015) determinants of Ty transposition rates. Second, environmental factors can influence Ty transposition rate, as demonstrated by temperature (Paquin and Williamson, 1984) and nitrogen starvation (Morillon et al., 2000) for Ty1 in S. cerevisiae. If hybrid lineages experience distinct environmental conditions in nature, this may impact their transposition rates. Finally, the efficiency at which natural selection purges deleterious Ty insertions (which depends on Ne) could vary between lineages, ultimately affecting the transposition-selection balance. To help discriminate global from lineage-specific patterns arising from these factors, analyses with all lineages aggregated were complemented by lineage-specific analyses. However, hybrid lineages were represented by relatively few strains (≤10 in each case), which made it hard to reach meaningful conclusions for some analyses. Thus, some of the lineage-specific analyses were focused on the more abundant SpB and SpC pure lineages.

We first investigated the relative contribution of population structure and environmental factors. We combined genome-wide single nucleotide variants from our population genomic data with climatic (temperature, precipitation, solar radiation, vapor pressure, wind) and water balance data (Abatzoglou et al., 2018). We reduced the dimensionality of both datasets by PCA, using coordinates on the major principal components (PCs) to describe population structure (Figure 4—figure supplement 1) and climatic variation (Figure 4—figure supplement 2, Figure 4—figure supplement 3). We then fitted linear models to compare how genotypes-based PCs (gPCs) and environment-based PCs (ePCs) can explain CN variation in Ty1, Ty3p and Tsu4 (Figure 4a). Globally, gPCs reached higher levels of statistical significance as linear predictors of CN compared to ePCs. gPC1 both explained a high proportion of variance (52%, 9% and 36% for all lineages, SpB and SpC, respectively; Figure 4—figure supplement 1) and was a highly significant linear predictor of Ty1 CN. Linear models that explained Ty CNs with both gPCs and ePCs had generally lower Akaike information criterion (AIC) values than models that included only gPCs (Supplementary file 1b), suggesting that climatic variables have modest although significant contributions in explaining CN variation. Unsupervised clustering of strains along gPCs consistently yielded more homogeneous clusters than ePCs in terms of Ty CN (Figure 4b, Figure 4—figure supplement 4). Overall, these results show that the genetic structure of North American S. paradoxus is generally a much better predictor of Ty CNs across lineages than climatic variables. Since locations where hybrid lineages were sampled have climatic profiles similar to pure lineages (Figure 4—figure supplement 2), exposure to specific environmental conditions appears unlikely to explain the absence of observable Ty reactivation in hybrid lineages.

Figure 4 with 4 supplements see all
Associations between population structure, climatic variables and Ty CNs for all lineages (left), SpB only (middle), or SpC only (right).

(a) Linear models fitted on Ty CNs with gPCs and ePCs as independent variables. For each PC, the p-value in the linear model is plotted against the proportion of variance explained in the respective PCA. (b) Unsupervised clustering performed on either gPC or ePC coordinates. Cumulative distributions show the standard deviation in Ty CN within each gPC-cluster or ePC-cluster. Distributions are also shown as boxplots, with whiskers spanning 1.5 times the interquartile range. FDR-corrected p-values of two-tailed Mann-Whitney U tests between gPC-cluster and ePC-cluster CN distributions are shown. ns: p≥0.05, *: p<0.05, **: p<0.01.

Variation in natural selection efficiency does not explain TE abundance variation in wild lineages

Most TE insertions are expected to have deleterious effects on the individual in which they occur (Wilke and Adams, 1992), although some can be advantageous (Aminetzach et al., 2005; Filteau et al., 2015; Gresham et al., 2008; Van't Hof et al., 2016). Natural selection acts to purge deleterious insertions, and its efficiency depends on the strength of genetic drift opposing it (which is inversely proportional to Ne). Ty reactivation in hybrid lineages might be counterbalanced by higher efficiency of natural selection in these lineages, although there is no a priori reason to expect higher Ne in hybrids. We thus sought to investigate if CN variation across lineages can be explained by variation in Ne. Genome-wide nucleotide diversity at fourfold degenerate codon positions (πS) is related to long-term Ne and was used to approximate this parameter across diverse taxa (Lynch and Conery, 2003). Nucleotide diversity was also used at the intraspecific level in S. paradoxus to estimate Ne (Tsai et al., 2008). While pure lineages did exhibit variation in πS (Supplementary file 1c), there was no obvious negative relationship between Ty CNs and πS (Figure 5—figure supplement 1), as would be expected under a transposition-selection equilibrium (assuming negligible variation in mutation rate). For instance, nucleotide diversity was higher in SpC than in SpA, while SpC exhibited the highest Ty1 and Ty3p CN. While these rough Ne estimates should be interpreted with caution, they suggested that natural selection efficiency is likely not the main determinant in CN variation. However, our data indicated that some Ty families experienced recent expansions and thus may not be at equilibrium. Recent hybrid lineages also cannot be meaningfully characterized by long-term Ne, further limiting the usefulness of πS as a measure of the efficiency of natural selection.

We can expect that most of the fitness effects associated with TE insertions are caused by alterations of host functions encoded in cis. Perhaps the TE insertions with the largest functional alteration potential are those which interrupt ORFs of host genes. Across the six genome assemblies, no Ty insertion was found to interrupt an ORF. Lineage-wide Ty insertion calls revealed four Ty insertions predicted to be within host ORFs, three of which are found in SpC strains and one in SpB. All the ORF-disrupting insertions were private alleles, and in all cases the strain harboring them had a CN higher than the median of its lineage (Figure 5—figure supplement 2). Examination of depth of coverage at these four loci revealed narrow peaks in the corresponding strains that were consistent with target site duplications (TSDs) arising from Ty integration (Figure 5—figure supplements 36). The absence of ORF-disrupting Ty insertions in hybrid lineages and their scarcity in pure lineages suggest that natural selection is globally very efficient at purging these mutations.

The number of ORF-disrupting Ty insertion alleles in SpC (three) compared to SpB (one) suggested a slight overrepresentation in SpC, which could be explained by a bias towards the accumulation of deleterious insertion alleles in this lineage. To test this hypothesis, we examined the functional annotations of the corresponding genes. There were striking patterns of functional relatedness among genes harboring an ORF-disrupting insertion. Two SpC strains had independent Ty1 insertions at different loci within YDL024C (DIA3, which encodes a protein of unknown function [Palecek et al., 2000]), as evidenced by the corresponding TSDs (Figure 5—figure supplement 3, Figure 5—figure supplement 4). The two other genes were involved in resistance to arsenic compounds: 16_037 (SpB) had a Tsu4 insertion in YPR200C (ARR2, which encodes an arsenate reductase [Mukhopadhyay and Rosen, 1998], Figure 5—figure supplement 5), while LL2011_004 (SpC) had a Ty1 insertion disrupting YPR201W (ARR3, which encodes an arsenite/H+ antiporter [Wysocki et al., 1997], Figure 5—figure supplement 6). This functional relatedness could be explained by positive selection; however, this scenario is incompatible with the extremely low frequency at which these insertions segregate. Thus, the most parsimonious explanations for these patterns appears to be either a mutational bias favoring Ty insertions at these loci or fitness effects too small to be visible to natural selection.

We investigated the fitness effect of the ORF-disrupting Ty1 insertion in ARR3 by measuring the growth of a collection of 123 S. paradoxus strains in culture medium supplemented with arsenite (Figure 5a, Figure 5—figure supplement 7, Figure 5—figure supplement 8). Growth in 0.8 mM sodium meta-arsenite (NaAsO2) showed significant variation between lineages (mixed design ANOVA, p=9.54 × 10−39), and more total variance was explained by variation between lineages (lineages: 86.58%, strains: 10.20%, residuals: 3.22%). SpC and SpC* lineages exhibited poor growth in 0.8 mM arsenite compared to SpB and SpA (Figure 5a). This phenotype was generalized to the whole lineages, and while LL2011_004 was in the lower part of the SpC distribution, it was not significantly different from many SpC strains not harboring an insertion (Figure 5b, Figure 5—figure supplement 9). This indicated that the disruption of a major gene involved in arsenite resistance had no further deleterious effect on the poor growth ability of SpC-related strains in presence of arsenite. The absence of growth defect in arsenite for LL2011_004, combined with the extremely low frequencies of all the ORF-disrupting alleles, suggests that these insertions are effectively neutral and are not observed because of an important reduction in purifying selection efficiency.

Figure 5 with 11 supplements see all
Genomic patterns of Ty insertions indicate that natural selection efficiency cannot explain CN variation in natural lineages.

(a) Difference between growth curves in 0.8 mM NaAsO2 (As) and control medium (synthetic complete, SC). Curves depict the difference in culture optical density (OD) through time, averaged across replicates and bins of four timepoints (one hour). 123 strains were assayed. SpC strain LL2011_004, which harbors a Ty1 insertion in ARR3, is represented by a dashed line. (b) Area under the curve (AUC) of the difference between growth curves in As and SC. Replicates of 107 strains that passed the criteria of number of replicates and homoscedasticity are shown. The number of strains per lineage is shown between parentheses. Replicates of the SpC strain LL2011_004 are highlighted with large empty circle symbols. Lowercase letters indicate groups that are not significantly different following Tukey’s HSD test. (c) Counts of essential and non-essential genes in the immediate vicinity of Ty insertions across the six whole-genome assemblies. Essential to non-essential ratios are shown on top of the bars. (d-e) Cumulative distributions of number of protein physical interactions (d) and ω (e) for genes neighboring conserved (dashed line) or private (solid lines) Ty insertions. Genes containing ORF-disrupting insertions are labeled with symbols. Heatmaps show FDR-corrected p-values for pairwise Mann-Whitney U tests between distributions. Color maps are centered at the significance threshold of 0.05. (f) Cumulative distributions of Δk between SpB and SpC. (g) Counts of neighboring genes with positive or negative Δk values.

To gain a more comprehensive view of the strength of selection acting on Ty insertions, we compared the properties of genes in their neighborhood. Ty elements contain many transcription factor binding sites and can affect the expression levels of adjacent genes (Morillon et al., 2000; Coney and Roeder, 1988; Todeschini et al., 2005; Lesage and Todeschini, 2005). Dosage-sensitive and functionally central genes should on average be subject to stronger selection against mutations that affect their regulation. For instance, in S. cerevisiae, genes that are essential or that have many protein physical interactions were shown to be depleted in eQTLs, i.e. variants that affect expression levels (Albert et al., 2018). However, the association of Ty insertions with regulatory effects on adjacent genes is expected to be highly variable, as transcriptional effects are known to depend on the location, orientation and nucleotide sequence of the inserted element (Coney and Roeder, 1988; Lesage and Todeschini, 2005). Nevertheless, if selection efficiency explains a significant part of Ty CN variation, we reason that when CNs are high (i.e. at low Ne), Ty insertions should be more frequently observed near genes which properties reflect stronger dosage and functional constraints, compared to when CNs are low (i.e. at high Ne). Focusing on the Ty insertions private to each lineage, we found no significant difference in the proportion of neighboring genes annotated as essential in S. cerevisiae (p>0.05, pairwise Fisher’s exact tests, FDR correction), although SpC did exhibit the highest essential to non-essential ratio (Figure 5c). Distributions of physical distance to neighboring genes were similar across lineages, as were distributions of relative orientation (Figure 5—figure supplement 10). Additionally, genes neighboring Ty insertions private to SpC show no bias towards more protein physical interactions (Figure 5d, p>0.05, pairwise two-tailed Mann-Whitney U tests, FDR correction). We fitted the ratio of non-synonymous to synonymous substitutions (ω) to each neighboring gene to estimate the strength of purifying selection on protein-coding sequences. The ω distributions for genes neighboring private Ty insertions exhibited no significant bias (p>0.05, pairwise two-tailed Mann-Whitney U tests, FDR correction). In particular, the distribution in SpC exhibited no skew towards lower values, suggesting no association with genes under strong purifying selection (Figure 5e). Notably, genes harboring ORF-disrupting insertions show little protein physical interactions and no signature of purifying selection (Figure 5d–e).

We next investigated whether lineage-specific selection signatures on genes are associated with Ty insertions. We used the RELAX branch-site codon model (Wertheim et al., 2015) to derive estimates of selection intensity (k) along individual gene trees. Within the RELAX framework, k < 1 indicates relaxed selection and k > 1 indicates strengthened (positive or purifying) selection. We define Δk = kSpB−kSpC, corresponding to the difference between the average k values for clades comprising SpB and SpC strains, respectively. Overall, the distribution of Δk values for genes neighboring private SpC insertions was very similar to that for conserved insertions (Figure 5f), while the distribution for private SpB insertions was slightly biased toward genes relaxed in SpB. For both SpB and SpC, there were more private insertions adjacent to genes experiencing relaxed selection in the corresponding lineage, although comparisons with conserved insertions were not statistically significant (p>0.05, Fisher’s exact tests, Figure 5g). This is consistent with efficient purifying selection acting against insertions near genes in both lineages. While genes harboring ORF-disrupting insertions have signatures of stronger selection in SpC (Figure 5e, Figure 5—figure supplement 11), this observation is mitigated by their far-right position in the ω distribution (Figure 5d). Overall, these results show that variation in natural selection efficiency is unlikely to explain a large part of the variation in Ty CNs in natural lineages.

TE content variation in experimental evolution under relaxed selection reveal no systematic effect of hybridization but major effects of individual genotypes

Thus far, our results showed that no reactivation is apparent in natural hybrid lineages, that population structure over environmental variables best explains variation in Ty CNs and that natural selection efficiency is unlikely to explain much variation in Ty genomic landscapes. However, natural selection seems equally efficient at purging deleterious Ty insertions in all lineages, leaving the possibility that strong purifying selection could mask any reactivation in hybrids. Transient reactivation may also be followed by rapid re-establishment of host repression and leave little to no signature in contemporary genomes. Moreover, natural hybrids are somewhat limited in their ability to uncover the interplay between hybridization and transposition, as they are relatively recent introgressed lineages that resulted from admixture between closely related lineages. To systematically explore the effect of hybridization on Ty transposition rates, we performed a large-scale mutation accumulation (MA) experiment on 11 artificial hybrid backgrounds, eight of which were previously published (Charron et al., 2019). A MA experiment with yeast consists in subjecting a population to periodic single-cell bottlenecks (typically during hundreds to thousands of mitotic generations) by repeatedly streaking one randomly chosen isolated colony, thus amplifying genetic drift and making natural selection inefficient (Lynch et al., 2008; Figure 6a). S. paradoxus was inferred to reproduce mitotically 99.9% of the time in the wild (Tsai et al., 2008), making MA experiments a reasonably realistic approximation of its natural life cycle. Pairs of haploid derivatives of natural S. paradoxus and S. cerevisiae strains were chosen to span a range of evolutionary divergence (Figure 6b), yielding four categories of intraspecific crosses (from the least to the most divergent): SpC × SpC (CC), SpB × SpB (BB), SpB × SpC (BC) and SpB × SpA (BA); and one category of interspecific crosses: SpB × S. cerevisiae (BSc; Figure 6c, Supplementary file 1d-f). The crosses used also spanned a range of variation in parental Ty CNs (Figure 6d). For each cross, 48 to 96 MA lines derived from independent matings were propagated for ~770 mitotic divisions.

Figure 6 with 1 supplement see all
Ty copy number (CN) variation in mutation accumulation (MA) lines reveal no effect of parental divergence and major genotype-specific effects.

(a) Principle of an MA experiment on a yeast population. N corresponds to population size. (b) Bayesian phylogenetic tree of the parental haploid strains used to generate artificial hybrids. (c) Design of the crosses for the MA experiment. (d) Ty CNs in the haploid parental strains of each cross in c. (e) CN at the onset (Tini, empty boxes) and at the end (Tend, filled boxes) of the MA experiment. Whiskers span 1.5 times the interquartile range. FDR-corrected p-values of Wilcoxon signed-rank tests and Brown–Forsythe tests for the equality of variances are shown in black and blue, respectively. ns: p≥0.05, *: p<0.05, **: p<0.01, ***: p<0.001. (f) Linear models explaining the mean variation in CN (ΔCN = CNTend−CNTini) for a given cross as a function of the joint CN and genetic divergence of its parental strains. The p-value for each term in the models is shown. Scatter plots (CCPR plots) show the partial residuals for each independent variable, along with the fitted regression.

We selected a random subset of 23 to 40 MA lines per cross and performed whole-genome resequencing on stocks frozen at Tini and Tend, totaling 372 strains sequenced at each timepoint (PRJNA515073). All these strains were confirmed to be diploid at Tini and Tend by quantifying total DNA content using flow cytometry (Charron et al., 2019; Marsit et al., 2020). Focusing on full-length Ty1, Ty3p, and Tsu4, no global trend emerged from the variation of NRD before and after MA (Figure 6e, Figure 6—figure supplement 1, Supplementary file 1f). Mixed effects linear models fitting cross-specific intercepts and slopes showed no statistically significant support for a global effect of time on CNs. Nevertheless, specific combinations of crosses and Ty families showed significant differences in either central tendency or variance. Often, crosses within a category of parental strain divergence showed contrasting patterns. For example, although CC1, CC2, and CC3 are genetically very similar SpC × SpC crosses, only CC1 exhibits a significant decrease in Ty1 and Tsu4 CNs. There were also cases of opposite trends within the same category of cross, for instance Ty1 CN in BC and BA crosses. This showed that individual parental genotypes, even if they are evolutionarily closely related, can have very different dynamics, both in sign and magnitude.

We next sought to understand how these dynamics are determined by properties of the different hybrid genotypes tested. We fitted linear models to the variation in CN between Tini and Tend (ΔCN), using the sum of CNs of both parental strains of a cross and the evolutionary divergence between them as independent variables. Although these models explain little global variance (R2 = 0.40, 0.08, and 0.06 for Ty1, Ty3p, and Tsu4, respectively), we found that neither independent variable was a significant predictor of ΔCN (Figure 6f). Only Tsu4 showed a slightly positive regression slope with respect to evolutionary divergence, and none of them was statistically significant, indicating that divergence between parents of a hybrid has no large influence on its CN dynamics. On the other hand, Ty1 ΔCN had a negative relationship with parental Ty1 CN and was marginally significant (p=0.056), suggesting that the initial CN of a hybrid genotype can dictate the sign and magnitude of its CN change. However, this result should be interpreted with care as some of the crosses shared common parental strains and were thus not entirely independent. This relatedness may in fact explain part of the pattern seen for Ty1. Crosses that shared the SpB parental strain UWOPS-91–202 (BB2, BC2, BA2, BSc2) seemed to have systematically higher Ty1 ΔCN values than the corresponding crosses that shared the SpB parental strain MSH-604 (BB1, BC1, BA1, BSc1). UWOPS-91–202 appeared to be a naturally Ty-less strain, as no full-length element was found for any family in its genome (Figure 6d). Previous reports of such Ty-less strains were shown to be permissive to the accumulation of Ty elements when assayed in the laboratory using a reporter assay (Czaja et al., 2020; Garfinkel et al., 2003). This observation also supports the hypothesis that genotype-specific characteristics can be major determinants of variation in Ty CN.

Discussion

The literature on TE reactivation in hybrids demonstrates its relevance to the field of speciation (Serrato-Capuchina and Matute, 2018). Reactivation can be a source of postzygotic reproductive isolation by compromising genome stability (Dion-Côté and Barbash, 2017), which translates into a spectrum of deleterious hybrid phenotypes spanning sterility, developmental abnormalities and inviability. Many studies investigated the TE reactivation hypothesis at the interspecific level and yielded equivocal evidence. For instance, mixed evidence for TE reactivation was found within the Drosophila genus (Coyne, 1989; Hey, 1988; Labrador et al., 1999; Vela et al., 2014), where reactivation seems to depend on the species pair involved. On the other hand, three well-established hybrid species that emerged from parallel hybridization between the same species pair in the genus Helianthus all exhibit important reactivation of LTR retrotransposons, both at the genomic (Ungerer et al., 2009; Ungerer et al., 2006) and transcriptomic (Renaut et al., 2014) levels. However, a study of hybrid zones between the two parental species revealed that no reactivation is apparent in recent hybrid populations, casting doubt on the idea that hybridization alone triggered reactivation and suggesting the involvement of other factors like abiotic stress (Kawakami et al., 2011). Overall, this data nuances the hypothesis that hybridization could act as a systematic trigger of TE proliferation and reproductive isolation.

In this study, we used population genomic data on North American lineages of the yeast S. paradoxus to test whether TEs were reactivated in natural hybrids. It was previously shown that these lineages experienced two hybridization events that resulted in the establishment of admixed lineages that have partial postzygotic reproductive isolation with their respective parental lineages (Eberlein et al., 2019; Leducq et al., 2016). We found no evidence of TE reactivation in any of these natural hybrid lineages. Our data showed that while some hybrid lineages may have relatively high CNs and exhibit evolutionary dynamics consistent with recent expansions, these signatures were never more pronounced than that of their parental lineages (SpC in the case of SpC* and SpC* in the case of SpD1 and SpD2). Moreover, our analyses suggest that environmental and demographic factors are unlikely to be major factors affecting TE CNs in natural lineages, which therefore likely reflect the underlying transposition rates. Our observations would seem consistent with a role of TE reactivation in speciation, as natural hybrids with high levels of reactivation would likely be unobserved. However, they contrast with long-established hybrid species in Helianthus, showing that reactivation and long-term maintenance can be compatible.

We note that the analysis of contemporary genomes is inherently limited in the ability to uncover transient reactivation that would immediately follow hybridization and rapidly faint, perhaps as a result of a re-establishment of TE repression. However, this scenario may be unlikely in yeast, given the absence of many epigenetic and post-transcriptional repression mechanisms (see discussion below) that can yield such a rapid re-establishment in other species. For instance, the re-establishment of piRNA-mediated P element repression in Drosophila P-M dysgenic hybrids occurs in single female individuals because paternally inherited piRNA clusters become expressed as they age (Khurana et al., 2011). If a hypothetical re-establishment of Ty repression was rather driven by selection for genomic variants, this should result in a gradual improvement of repression with time. Our data provides no support for this hypothesis, as younger hybrid lineages SpD1 and SpD2 do not show more pronounced signatures of recent Ty expansion compared to the older SpC* lineage. Additionally, Ty CNs of artificial hybrids at the initial time point of our MA experiment (i.e. early after hybridization) generally show good agreement with the expectation from their parental haploid strains, which is incompatible with an early reactivation of Ty elements.

Many factors are susceptible to determine whether a hybrid genotype exhibits reactivation or not. According to the Bateson-Dobzhansky-Muller (BDM) model, genetic incompatibilities should accumulate as populations diverge (Orr, 1995). Additionally, the evolutionary conflict between host genomes and TEs results in an arms race that can accelerate the evolution of the genes involved (Blumenstiel et al., 2016; Van Valen, 1973). Thus, it can be expected that incompatibilities affecting transposition rates are more likely to be observed in hybrids between highly diverged genotypes compared to closely related ones. Here, we report data from MA experiments on artificial hybrids that probe the relationship between evolutionary divergence and transposition rates. We show that contrary to the expectation from the BDM model, hybrids between more divergent genotypes exhibit no deregulation of TEs. This conclusion is consistent with recent results from Smukowski Heil and colleagues, who showed that no Ty1 reactivation occurs in experimental interspecific hybrids between S. cerevisiae and S. uvarum (Smukowski Heil et al., 2020). Instead of reactivation, many hybrid genotypes exhibited marked decreases in TE CNs after evolution. This trend is especially important for Ty3p, as eight out of ten genotypes exhibit a negative change in median CN. Ty3p is also the least abundant Ty family in natural lineages, as it is found in rarely more than ~5 copies and extinct in all hybrid strains sampled. Such low abundance families are likely more prone to extinction, as stochastic fluctuation in CN would more frequently lead to the loss of all copies. Stochastic loss of Ty copies may come from intra-element LTR recombination or loss of heterozygosity at a locus where a Ty element is present in one parent only, although the mechanistic causes of CN decrease remain to be investigated. This result also seems consistent with the reported Ty CN decrease in mitotically propagated experimental populations of S. cerevisiae (Bast et al., 2019), suggesting a potential intrinsic bias towards TE loss. However, it should be noted that the latter study used haploid populations in a regime with efficient natural selection (McDonald et al., 2016), thus limiting the relevance of comparisons with our results.

Transposition is a complex cellular process, with multiple steps regulated by both TEs and host genomes. Genetic screens in S. cerevisiae revealed that genes involved in multiple cellular functions (chromatin organization, transcription, translation, RNA and protein metabolism) affect Ty1 transposition rate, either positively or negatively (reviewed in Curcio et al., 2015). Thus, it is possible that divergent hybrids suffer from misregulation at different steps of the transposition process that have antagonistic effects, such that they cancel out at the level of transposition rates. It is also possible that such antagonistic interactions exist between transposition and excision rates, yielding no net change in CN. Answering these questions will require a thorough investigation of novel Ty insertion loci in the MA lines, which is beyond the scope of this study. Transposition regulation being a polygenic trait with a large mutational target is also consistent with the strong genotype specificity we observed in our artificial hybrids. Some of the contrasted dynamics between crosses that share high genetic similarity might originate from cryptic genetic variation in genes involved in Ty regulation, whose effect on transposition rate would be uncovered in hybrids.

TEs are structurally very diverse and populate the genomes of a vast range of species, and the control mechanisms in response to their mobilization may exhibit a correspondingly large diversity. In many eukaryotes, TEs are repressed epigenetically through mechanisms that include DNA methylation, chromatin modification, and mRNA degradation, processes in which small RNAs play central roles (Lisch, 2009; Slotkin and Martienssen, 2007). In many fungi, TE repression acts via multiple mechanisms involving RNA interference (RNAi) or cytosine methylation, including repeat-induced point mutation (RIP) (Gladyshev, 2017). Notably, introgressed species of the fungal genus Coccidioides harbor genomic signatures of efficient RIP (Neafsey et al., 2010), suggesting that this mechanism may remain active in the face of hybridization. Unlike other eukaryotes, S. cerevisiae and its closely related species lost the genes of the RNAi pathway (Drinnenberg et al., 2009) and have no cytosine methylation (Proffitt et al., 1984). Instead, S. cerevisiae Ty1 has a potent self-repression mechanism, with a strength that depends on the actual Ty1 CN in the host genome (thus termed CN control, CNC) (Garfinkel et al., 2003). CNC is a post-transcriptional mechanism that relies on p22, an N-truncated form of the Gag (p49) Ty1-encoded capsid protein which is translated from an alternative transcript of Ty1 (Saha et al., 2015). p22 inhibits transposition by interfering with the assembly of the cytoplasmic virus-like particles that are essential to the replication cycle of Ty1 (Saha et al., 2015). This unique repression mechanism may explain some of the patterns we observed in the artificial hybrids. In particular, our data suggests a negative relationship between the initial Ty1 CN in hybrid genotypes and the change in Ty1 CN between the onset and the end of the MA experiment. This pattern seems consistent with an efficient action of CNC, although our results would then suggest that CNC efficiency is largely independent of the degree of evolutionary divergence between the subgenomes of a hybrid. This contrasts with the recent finding that interspecific sequence divergence in p22 has a dramatic impact on the action of CNC (Czaja et al., 2020), although such impact is likely mitigated by the generally much lower divergence levels in our crosses. Thus, further experimental work is needed to disentangle the mechanistic contributions of Ty1 CN and sequence divergence in the regulation of transposition rates in Saccharomyces hybrids.

In conclusion, our work illustrates how the use of natural model systems that are amenable to experimental investigations can bring novel insights on long-standing evolutionary hypotheses. Our results show that natural hybrid lineages of the wild yeast S. paradoxus exhibit no TE reactivation and suggest that dynamics of TE content evolution in hybrids depend on deterministic factors that evolutionary divergence between parents can hardly predict.

Materials and methods

Key resources table
Reagent type (species)
or resource
DesignationSource or referenceIdentifiersAdditional information
Sequence-based reagentForward primer to construct
the deletion cassette for
the ADE2 gene in SpC
Charron et al., 2019CLOP97-F5acaattaaggaatcaagaaaccgt
gataaaaaattcaagt
CAGCTGAAGCTTCGTACGC
Sequence-based reagentReverse primer to construct
the deletion cassette for
the ADE2 gene in SpC
Charron et al., 2019CLOP97-F6gtaattgttcgctggccaagtata
ttaatacatttatata
GCATAGGCCACTAGTGGATC
Strain, strain background
(Saccharomyces paradoxus)
Haploid parent for MA linesCharron et al., 2014LL2011_001 MATa hoΔ::kanMX4
Strain, strain background
(Saccharomyces paradoxus)
Haploid parent for MA linesThis studyLL2011_001 MATa
hoΔ::kanMX4 ade2Δ::hphNT1
Haploid yeast strain with ho
and ade2 deletions.
See Materials and methods
section Mutation accumulation
Strain, strain background
(Saccharomyces paradoxus)
Haploid parent for MA linesLeducq et al., 2016LL2011_012 MATa hoΔ::kanMX4
Strain, strain background
(Saccharomyces paradoxus)
Haploid parent for MA linesThis studyLL2011_012 MATa
hoΔ::kanMX4 ade2Δ::hphNT1
Haploid yeast strain with ho and
ade2 deletions. See Materials
and methods section
Mutation accumulation
Strain, strain background
(Saccharomyces paradoxus)
Haploid parent for MA linesCharron et al., 2014MSH-587–1 MATa hoΔ::natMX4
Strain, strain background
(Saccharomyces paradoxus)
Haploid parent for MA linesThis studyMSH-587–1 MATa
hoΔ::natMX4 ade2Δ::hphNT1
Haploid yeast strain with ho and
ade2 deletions. See Materials
and methods section
Mutation accumulation
Strain, strain background
(Saccharomyces paradoxus)
Haploid parent for MA linesCharron et al., 2019LL2011_004 MATα
hoΔ::kanMX4 ade2Δ::hphNT1
Strain, strain background
(Saccharomyces paradoxus)
Haploid parent for MA linesCharron et al., 2019LL2011_009 MATα
hoΔ::natMX4 ade2Δ::hphNT1

General tools

Request a detailed protocol

GNU Parallel version 20170922 (Tange, 2011) and SnakeMake version 5.5.4 (Köster and Rahmann, 2012) were used extensively in the analyses. Unless otherwise stated, Python version 3.7.6 (Van Rossum and Drake, 2009) was used for the custom scripts. These scripts were run on JupyterLab version 2.0.1 (Kluyver et al., 2016) and made extensive use of the following packages: BioPython version 1.76 (Cock et al., 2009), DendroPy version 4.4.0 (Sukumaran and Holder, 2010), Matplotlib version 3.2.0 (Hunter, 2007), NetworkX version 2.4 (Hagberg et al., 2008), NumPy version 1.18.1 (van der Walt et al., 2011), pandas version 1.0.2 (McKinney et al., 2010), SciPy version 1.4.1 (Virtanen et al., 2020), seaborn version 0.10.0 (Waskom et al., 2020). Custom scripts can be found in Supplementary file 2

Protein-coding genes annotation and orthogroup definition

Request a detailed protocol

Whole-genome de novo assemblies were previously generated for six North American S. paradoxus strains using a combination of long (Oxford Nanopore) and short (Illumina) reads (PRJNA514804). Protein-coding genes were predicted ab initio using AUGUSTUS version 3.3.1 with parameters species=saccharomyces genemodel=partial singlestrand=false. Predicted ORF sequences for the six assemblies were combined with ORF annotated from the S. cerevisiae strains S288c and YPS128 from Yue et al., 2017 as references. Protein-coding genes orthogroups were inferred from sequence similarity and synteny using the program SynChro version Jan2015 (Drillon et al., 2014) with parameters a=0 DeltaRBH=1. Orthogroups were defined as genes with one-to-one relationships across all pairwise comparisons between the eight genomes. Protein-coding genes systematic names were propagated from S288c to the corresponding orthogroups.

Annotation of Ty LTR retrotransposons

Request a detailed protocol

Ty LTR retrotransposons were annotated from the six whole-genome assemblies using RepeatMasker version open-4.0.7 (Smit et al., 2015) with parameters -s -low -gccalc and the library of reference LTR and internal sequences assembled by Carr et al., 2012 (containing Ty1, Ty2, Ty3, Ty4, and Ty5 from S. cerevisiae and Ty3p from S. paradoxus), additionned with Tsu4 from S. uvarum (Bergman, 2018). RepeatMasker output was defragmented using REannotate version 26.11.2007 (Pereira, 2008) with parameters -s 15000 -d 10000, yielding defragmented elements classified as full-length, solo LTR, truncated or orphan internal sequences. The -f parameter was used to allow defragmentation of Ty1/Ty2 hybrid elements following Carr et al., 2012. The defragmented annotations were manually curated to split defragmented elements spanning entire protein-coding genes or forming concatemers.

Network-based Ty LTR sequence similarity analysis

Request a detailed protocol

We analyzed the sequence diversity of LTR annotations within each of our six S. paradoxus assemblies using a network approach (Levy et al., 2017). The sequences of all LTR annotations within a genome were extracted and all-against-all search was performed between LTR sequences in a genome using BLASTN version 2.7.0 (Camacho et al., 2009) with parameters 'word_size 7 dust no'. Networks were built using the Python package Networkx version 2.4 (Hagberg et al., 2008) by keeping hits with a bit score higher or equal to the 10th percentile of all bit scores for a given LTR. Visualization of the resulting networks was done using Cytoscape version 3.7.1 (Shannon et al., 2003) using the Prefuse Force Directed layout algorithm. Inspection of these networks justified the merging of the following LTR families: Ty1-Ty2, Ty3p-Ty3, and Tsu4-Ty4.

Estimation of evolutionary divergence since transposition

Request a detailed protocol

Evolutionary divergence since transposition was estimated from inter-element LTR sequence identity, assuming that the LTRs of a newly inserted element are identical to those of the source element. From the sequence similarity graphs, the closest relative of each element was retrieved by keeping the hit with the highest bit score. Hits between the two LTRs of full-length elements were excluded, and only one LTR per full-length element was considered. For each element, the evolutionary divergence since transposition was extracted as the percent sequence identity to the closest relative.

Definition of Ty orthogroups from genome assemblies

Request a detailed protocol

Orthology of LTR retrotransposons in the six whole-genome assemblies was assessed as follows. Whole-genome multiple sequence alignment (MSA) with rearrangements was performed using Mugsy version 1r2.3 (Angiuoli and Salzberg, 2011) with parameters -d 20000 -nucmeropts ‘-b 1000 l 15’ -fullsearch -refine mugsy. MSA blocks for which at least one genome had no corresponding sequence were filtered out. The remaining alignment blocks were scanned in 10 kb sliding, non-overlapping windows in both directions and windows were trimmed until at least 80% of sites were informative in a MSA window. Ty annotation coordinates were transformed into an MSA coordinates system and clustered using the DBSCAN algorithm implemented in the Python package scikit-learn version 0.22.2 with parameters min_samples=1 eps=500. This one-dimensional clustering approach allowed to define orthologous clusters containing at least one Ty annotation in at least one genome.

We defined a set of confident orthogroups to derive an empirical distribution of physical distance between start and end coordinates of orthologs. Clusters were selected if they contained at most one annotation per genome and if all annotations were on the same strand and belonged to the same family. Pairwise BLASTN searches within each cluster were used to keep elements that produced at least one hit with ≥80% nucleotide identity and an alignment length between 90% and 110% of both query and subject lengths. From the 38 resulting confident orthogroups, we extracted the distribution of differences between annotation start and end coordinates and used the 95th percentile of that distribution (17 bp) as threshold for the following orthogroup definition procedure. Within each cluster, pairs of annotations meeting the following two criteria were considered orthologous: the pair had to include at least one best BLASTN hit, and either start or end coordinates (or both) had to be within the defined threshold of 17 nt. This method allowed to define 713 Ty orthogroups, 20 of which (3%) had >1 element for at least one genome and were discarded.

Phylogenetic analysis of Ty1 evolution

Request a detailed protocol

For each of our six genome assemblies, we extracted the LTR sequence of all Ty1 annotations, keeping only one LTR at random for full-length elements. As an outgroup, a conserved Ty4 LTR orthogroup was added on the basis of minimal variance in start and end coordinates across genomes and minimal difference between the length of the reference Ty4 sequence and the average length of the annotations. MSAs were computed using MUSCLE version 3.8.31 (Edgar, 2004) and filtered to maximize the proportion of informative sites while keeping the outgroup, first discarding sequences with >25% gaps and then removing sites with ≥25% gaps. Bayesian phylogenetic trees were computed using MrBayes version 3.2.5 (Ronquist et al., 2012) with model parameters Nst=6 Rates=invgamma. Markov chain Monte Carlo parameters were set to 10 million generations for two runs of four chains each with 25% burn-in fraction and temperature parameter of 0.03 to improve the efficiency of Metropolis coupling. Bayesian phylogenetic trees were converted into ultrametric phylogenetic trees using PATHd8 (Britton et al., 2007) and rooted at the Ty4 LTR outgroup. The branching time of each internal node in the ultrametric trees was computed as the distance from the root. Branching events were counted as the number of children nodes for each internal node, as the bayesian phylogenetic reconstructions yielded many polytomies at the base and tips of the trees.

Construction of a library of reference Ty sequences for S. paradoxus

Request a detailed protocol

We combined the Ty LTR sequences of all families and full-length sequences of Ty1, Ty3p and Tsu4 from the defragmented annotations of our six S. paradoxus assemblies and aligned the sequences of each family and type using MUSCLE. The resulting MSAs were filtered to remove sequences and sites that had ≥70% gaps. Consensus sequences from those MSAs were computed using EMBOSS cons version 6.6 (Rice et al., 2000) and the annotation sequence with the lowest BLASTN E value against the resulting consensus was included in the reference library. LTRs of the full-length reference elements were manually trimmed using MEGA version 7.0.26 (Kumar et al., 2016). The S. paradoxus sequences were complemented with S. cerevisiae internal sequences for Ty2, Ty3, Ty4, and Ty5 (Carr et al., 2012).

Ty insertion allele calling and population structure analysis

Request a detailed protocol

Illumina paired-end 150 bp or 100 bp libraries for 207 wild S. paradoxus diploid strains were downloaded from NCBI (PRJNA277692, PRJNA324830, PRJNA479851) and trimmed using Trimmomatic version 0.36 (Bolger et al., 2014) with parameters ILLUMINACLIP:combined_adapters.fa:2:30:10 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36 and a custom library of Illumina adapter sequences. Trimmed reads were mapped against the masked version of our SpB genome assembly (strain MSH-604, from the RepeatMasker output) using bwa mem version 0.7.16a-r1181 (Li, 2013) and sorted using samtools sort version 1.9 (Li, 2011). Duplicates were removed using Picard RemoveDuplicates version 2.18.5-SNAPSHOT (Broad Institute, 2019) with parameter REMOVE_DUPLICATES=true. RetroSeq version 1.5 (Keane et al., 2013) was used to call Ty insertion alleles from discordant pairs of mapped reads. The -discover command was used with parameter -eref and our custom reference Ty sequences library. The -call command was run with parameter -soft. Visual inspection of the resulting calls revealed several instances of multiple, tightly clustered calls of the same family, likely reflecting a single true positive call. To minimize false positives, the calls were transformed into MSA coordinates and clustered using DBSCAN with parameters min_samples=1 eps=500. Resulting clusters were considered true positives if they contained at least one original call of high-quality (flag 8). This strategy provided the best balance between false positive rates and false negative rates in an analysis made with various simulated short-read datasets (data not shown). Insertion allele frequency spectra were computed within each lineage. PCA was performed on Ty alleles using the Python package scikit-learn version 0.22.2.post1.

Population structure based on genome-wide SNPs

Request a detailed protocol

Trimmed reads from the 207 wild strains were mapped against the unmasked version of our SpB genome assembly (strain MSH-604) using bwa mem version 0.7.16a-r1181 and sorted using samtools sort. Duplicates were removed using Picard RemoveDuplicates with parameter REMOVE_DUPLICATES=true. Variants were called using freebayes version 1.3.1–17-gaa2ace8 (Garrison and Marth, 2012). Variants were filtered using bcftools view version 1.9 (Li, 2011) with parameter -e ‘QUAL<=20’ -g ^het, along with the tools vcfallelicprimitives, vcfbreakmulti and vcffilter with parameter -f ‘TYPE = snp’ from vcflib version 1.0.0_rc2 (Garrison, 2018). PCA on genome-wide SNPs was performed using the Python package scikit-allel version 1.2.1 (Miles et al., 2019), keeping only biallelic SNPs and excluding singletons. From the filtered SNPs of the 207 wild strains, 10000 SNPs were randomly selected and converted to NEXUS format using vcf2phylip version 2.0 (Ortiz, 2019). A bayesian phylogenetic tree was computed using MrBayes with model parameters Nst=6 Rates=invgamma. Markov chain Monte Carlo parameters were set to 1 million generations for two runs of four chains each with 25% burn-in fraction and temperature parameter of 0.1.

Depth of coverage of wild strains on Ty reference sequences

Request a detailed protocol

Trimmed reads from the 207 wild strains were mapped against our custom reference Ty sequences library using bwa mem version 0.7.16a-r1181 and sorted using samtools sort version 1.9. Only internal sequences extracted from full-length reference sequences were used. Depth of coverage was extracted using samtools depth. Median depth values of Ty elements were transformed in log2 and normalized by the log2 genome-wide depth computed from the mapped bases on the genome-wide mappings previously described. Kruskal-Wallis tests were performed on normalized log2 depth values for families with full-length elements in natural lineages (Ty1, Ty3p, Tsu4). Conover post-hoc tests were performed to test for pairwise differences between lineages using the Python package scikit-posthocs version 0.6.3 (Terpilowski, 2019).

Climatic data

Request a detailed protocol

We used the TerraClimate climatic dataset (Abatzoglou et al., 2018) (accessed in October 2019), which incorporates high spatial and temporal resolution climatic interpolation data from primary and derived climatic variables. For the sampling locations of our 207 wild strains, we aggregated the 12 climatic variables over the time period from 1958 to 2015 and computed the mean, standard deviation, minimum and maximum values over time for a given variable variable at a given location. The resulting dataset was normalized into z-scores and summarized by PCA using scikit-learn.

Association analysis between Ty CNs, population structure and climatic variables

Request a detailed protocol

SNP-based population structure PCAs and climatic variable PCAs were combined in linear models as predictors of CNs for the three active families Ty1, Ty3p and Tsu4. To avoid spurious effects arising from log transformation of very sparse coverage, CN values lower than −1.25 were replaced by −1.25 as for MA lines (see below). In each case, PCs explaining a cumulative sum of 80% of the total variance in the dataset were kept and named gPCs (from PCAs on SNP data) or ePCs (from PCAs on climatic variables). Ordinary least squares models were fitted using the Python package statsmodels version 0.10.1 (Seabold and Perktold, 2010). Two types of models were fitted: OLS using both gPCs and ePCs as predictors, and OLS using only gPCs as predictors. The quality of the two model types was compared using their AIC values. Strains were clustered according to their coordinates on gPCs using the AffinityPropagation algorithm from scikit-learn with parameters max_iter=10,000 damping=0.8, and subsequently according to their coordinates on ePCs using the KMeans algorithm from scikit-learn with the same number of clusters fitted for gPCs. The inconsistency of each cluster was computed as the standard deviation of CN for each Ty family.

Nucleotide diversity at synonymous sites

Request a detailed protocol

Long-term effective population size for pure lineages was estimated by computing nucleotide diversity at synonymous sites in protein-coding genes. Only SNPs at four-fold degenerate sites in protein-coding genes (third positions of codons CTN, GTN, TCN, CCN, ACN, GCN, CGN, and GGN) were considered. Nucleotide diversity was computed using the mean_pairwise_difference function from scikit-allel, excluding SpD2 strain WX21 since it is known to be heterozygous (Eberlein et al., 2019).

Determination of the optimal NaAsO2 concentration for growth assay

Request a detailed protocol

Three strains were chosen at random in each of the lineages SpA, SpB and SpC* and two SpC strains were chosen at random, totalling 12 strains with the SpC strain LL2011_004. Precultures were grown in synthetic complete medium (SC; 2% glucose, 0.174% yeast nitrogen base [without amino acids without ammonium sulfate], 0.5% ammonium sulfate, 0.134% complete drop-out) until they reached a 595 nm optical density (OD595) of ~0.5–1.0. Precultures were diluted at an OD595 of 0.1 ml−1 in SC medium supplemented with 0, 0.4, 0.6, 0.8 and 1 mM NaAsO2 in a single 96-well microplate, with border columns and rows filled with sterile SC medium. OD595 was monitored during 60 hr at 15 min intervals at 25°C in a Tecan Spark microplate reader (Tecan Group Ltd, Männedorf, Switzerland). OD595 values were normalized by subtracting the initial value for each well.

Growth assay in NaAsO2

Request a detailed protocol

From the 207 wild strains used for whole genome sequence analyses, 123 strains had glycerol stocks available. Five replicates of each strain were randomly arranged into four arrays of 14 × 11 positions, totalling 615 replicates with one empty position, with the constraint of having at most three replicates of any given strain within an array. For each array, glycerol stocks of the 123 strains were spotted on YPD plates (1% yeast extract, 2% tryptone, 2% glucose, 2% agar) and incubated at room temperature for three days. Cells from fresh spots were diluted in 200 μl SC medium. OD595 was monitored and 500 μl precultures were initiated at 0.08 OD595 ml−1 and incubated at room temperature for 10 hr. OD595 was monitored and precultures were diluted at 0.2 OD595 ml−1 in 1 ml SC medium. BRANDplates pureGrade S 384-well microplates (BRAND GMBH + CO KG, Wertheim, Germany) were filled with 40 μl of SC medium (even column numbers) or SC medium supplemented with 1.6 mM NaAsO2 (odd column numbers). 40 μl of cultures at 0.2 OD595 ml−1 were added to the plates using an Eppendorf epMotion 5075 liquid handling platform (Eppendorf, Hamburg, Germany). Replicates were positioned according to the predefined randomized scheme in side-by-side pairs, yielding adjacent replicates diluted at 0.1 OD595 ml−1 in SC and SC with 0.8 mM NaAsO2. Border rows and columns of the 384-well microplates were filled with 80 μl sterile medium. Microplate lids were treated with 0.05% Triton X-100 in 20% ethanol for 1 min and air dried (Brewster, 2003). OD595 was monitored during 70 hr at 15 min intervals at 25°C in a Tecan Spark microplate reader.

OD595 values were normalized by subtracting the initial value for each well, and values in SC were subtracted from values in NaAsO2 for each replicate. The area under the curve (AUC) of the resulting curves was obtained by summing all timepoints. Replicates arising from precultures that failed to grow were excluded, and strains which had less than three replicates left were excluded. The standard deviation of replicates was computed for each strain, and strains which fell outside the [−1, 1] z-score range were filtered out to meet the homoscedasticity requirements of analysis of variance (ANOVA). To assess the significance of lineage effects, a mixed effects linear model (with strains as random variable) was fitted in R version 3.4.3 (R Development Core Team, 2017) using the lmer function from the package lme4 version 1.1–15 (Bates et al., 2015) and an ANOVA table was produced using the package lmerTest version 3.1–2 (Kuznetsova et al., 2017). Variance components were estimated using the R package VCA version 1.3.4 (Schützenmeister and Piepho, 2012). Post-hoc Tukey’s HSD test was performed using scikit-posthocs.

Properties of genes neighboring Ty insertions

Request a detailed protocol

The analysis of properties of neighboring genes considered Ty-centric immediate adjacency relationships with protein-coding genes, meaning that a single gene could be represented more than once if many Ty annotations are immediately adjacent to it. From the whole-genome MSAs, protein-coding genes adjacent to Ty annotations were retrieved on the basis of synteny with the S. cerevisiae S288c genome. Genes annotated as essential in S. cerevisiae were extracted from YeastMine (accessed in August 2019). Physical interaction data for S. cerevisiae S288c were retrieved from the BioGRID Multi-Validated Physical interaction dataset version 3.5.188 (Oughtred et al., 2019) and the interactome was built using Networkx.

The strength of purifying selection acting on individual genes was assessed at two levels: at the species-wide level (i.e. including all North American S. paradoxus lineages) by fitting a single ω ratio (non-synonymous to synonymous substitution rates) per gene; and at the lineage level, by evaluating the relative strength of natural selection in the main pure lineages SpB and SpC. A randomly selected subset of 34 strains from the 207 wild strains were chosen. The VCF matrix from the SNP calling step above was used to derive the DNA sequences of protein-coding genes in these 34 strains. Neighbor-joining gene trees were produced using the nj function from the R package ape version 5.0 (Paradis and Schliep, 2019) and R version 3.4.3 (R Development Core Team, 2017). ω ratios were fitted using codeml from the PAML suite version 4.8 (Yang, 2007) with parameters seqtype=1 model=0. Selection strength at the lineage level was evaluated using the RELAX framework from the HyPhy software suite version develop-0d598ad (Wertheim et al., 2015). The RELAX general descriptive model was fitted using the RELAX-scan script (Pond, 2019). From the model output, the fitted value of the k parameter (describing the strength of selection) was extracted for each branch of the gene trees. For SpB and SpC, an average k value was computed by selecting the clade containing all the strains of the corresponding lineage (from the random subset) and weighting k values by the length of the corresponding branch. The difference in selection intensity between SpB and SpC was defined as Δk = kSpB−kSpC. For the three genes with ORF-disrupting Ty insertions, the RELAX framework was fitted on gene trees comprising all 207 wild strains.

Mutation accumulation

Request a detailed protocol

Crosses BB, BC, BA, and BSc were previously published (Charron et al., 2019). CC crosses were performed following the same methodology. The ADE2 gene was deleted in haploid MATa hoΔ::kanMX4 or MATa hoΔ::natMX4 derivatives of the SpC strains LL2011_001, LL2011_012 and MSH-587–1 (Charron et al., 2014) using a hphNT1 cassette that was amplified from plasmid pFA6a-hphNT1 (Janke et al., 2004) using primers CLOP97-F5 and CLOP97-F6 (Supplementary file 1d,e). The three strains were paired with the previously published SpC strains LL2011_004 MATα hoΔ::kanMX4 and LL2011_009 MATα hoΔ::natMX4 to generate three CC crosses by selection on medium containing 100 μg ml−1 geneticin and 10 μg ml−1 nourseothricin. For each cross, 64 MA lines were initiated by selecting single colonies from independent mating cultures that were spotted on selective medium. Diploid hybrids were propagated mitotically on Petri dishes containing solid YPD medium (1% yeast extract, 2% tryptone, 2% glucose, 2% agar) by streaking for single colonies. Every three days, a single colony was randomly chosen to inoculate a new plate by picking the colony that grew the closest to a predefined mark on the Petri plate. Colonies that had an obvious difference in red pigmentation were excluded to avoid the propagation of mitochondrial petite mutants. MA lines were propagated for 35 passages (approx. 770 mitotic generations), keeping glycerol stocks of all the lines every three passages.

Illumina library preparation for MA lines

Request a detailed protocol

For each cross, 23 to 40 MA lines were randomly chosen among the subset for which diploidy was confirmed at Tini and Tend (Charron et al., 2019; Marsit et al., 2020). The glycerol stocks from passage 1 (Tini) and passage 35 (Tend) were streaked on solid YPD medium. Single colonies were used to inoculate liquid YPD cultures. DNA was extracted from saturated cultures using QIAGEN Blood and Tissue kit with modifications for yeast cell wall digestion. DNA was treated with RNase A and purified on Axygen AxyPrep Mag PCR Clean-up SPRI beads. Genomic DNA libraries were prepared using the Illumina Nextera Tagmentase and Nextera index kits following a modified version of the protocol described in Baym et al., 2015. The quality of randomly selected libraries was controlled using an Agilent BioAnalyzer 2100 electrophoresis system. Sequencing was performed at the CHU Sainte-Justine and Génome Québec Integrated Centre for Pediatric Clinical Genomics (Montreal, Canada) on an Illumina NovaSeq 6000 sequencer with SP flow cells.

Depth of coverage of MA lines on Ty reference sequences

Request a detailed protocol

Illumina paired-end 150 bp libraries for 372 MA lines and their 13 parental haploid strains (deposited at NCBI Sequence Read Archive under PRJNA515073) were trimmed using Trimmomatic version 0.33 with parameters ILLUMINACLIP:nextera.fa:6:20:10 MINLEN:40 and a library of Illumina Nextera adapter sequences. Trimmed reads were mapped against our SpB genome assembly (strain MSH-604) using bwa mem version 0.7.17 and sorted using samtools sort version 1.8. Duplicates were removed using Picard RemoveDuplicates version 2.18.29-SNAPSHOT with parameter REMOVE_DUPLICATES=true. We extracted the list of S. cerevisiae genes with paralogs from YeastMine (accessed in December 2018) (Balakrishnan et al., 2012) and randomly selected 50 genes that did not figure in this list. The protein sequences of these genes were used as TBLASTN queries against our SpB genome assembly and the best hit loci were extracted. Depth of coverage over these loci was extracted using samtools depth version 1.9 and depth values were transformed in log2. Depth values were transformed in log2 and normalized by the log2 genome-wide depth computed from the mapped bases on the genome-wide mappings.

The trimmed reads of the 372 MA lines and their 13 parental haploid strains were mapped against our custom reference Ty sequences library using bwa mem version 0.7.16a-r1181 and sorted using samtools sort version 1.9. Only internal sequences extracted from full-length reference sequences were used. Median depth values of Ty elements were transformed in log2 and normalized by the log2 genome-wide computed from the mapped bases on the genome-wide mappings. Normalized log2 depth values for a given cross and Ty family were discarded if the median was lower than the 10th percentile of the normalized log2 depth values for non-paralogous genes in the corresponding cross, thus considering the family absent from the initial MA genotypes. Similarly, the 10th percentile of normalized log2 depth values for non-paralogous genes in all crosses aggregated (−1.25) was used as a lower threshold for depth values of the parental haploid strains, replacing lower values by −1.25.

Phylogenetic analysis of the MA lines haploid parents

Request a detailed protocol

From the mappings of the MA lines parental haploid strains, variants were called using freebayes. Variants were filtered using bcftools view with parameter -e ‘QUAL <= 20 g ^het’, along with the tools vcfallelicprimitives, vcfbreakmulti and vcffilter with parameter -f ‘TYPE = snp’ from vcflib. From the filtered SNPs, 10000 SNPs were randomly selected and converted to NEXUS format using vcf2phylip. A bayesian phylogenetic tree was computed using MrBayes with model parameters Nst=6 Rates=invgamma Ploidy=haploid. Markov chain Monte Carlo parameters were set to 1 million generations for two runs of four chains each with 25% burn-in fraction and temperature parameter of 0.1. From the resulting consensus tree, the evolutionary divergence between parental strains was measured as the sum of branch lengths connecting each pair.

Association analysis between Ty ΔCN, evolutionary divergence between parents and initial Ty CN

Request a detailed protocol

For each cross, the average ΔCN was analysed in a linear model incorporating parental evolutionary divergence and initial Ty CN of the hybrid genotypes as predictors. For each cross, the Ty CNs of parental haploid strains were summed and used as a measure of the initial Ty CN in all MA lines of the cross. An OLS model was fitted using statsmodels.

References

  1. 1
  2. 2
  3. 3
  4. 4
  5. 5
  6. 6
  7. 7
  8. 8
  9. 9
  10. 10
  11. 11
  12. 12
    What drives positive selection in the Drosophila piRNA machinery? the genomic autoimmunity hypothesis
    1. JP Blumenstiel
    2. AA Erwin
    3. LW Hemmer
    (2016)
    The Yale Journal of Biology and Medicine 89:499–512.
  13. 13
  14. 14
  15. 15
  16. 16
  17. 17
  18. 18
  19. 19
  20. 20
  21. 21
  22. 22
  23. 23
  24. 24
  25. 25
  26. 26
  27. 27
  28. 28
  29. 29
  30. 30
  31. 31
  32. 32
  33. 33
  34. 34
  35. 35
  36. 36
  37. 37
  38. 38
  39. 39
  40. 40
  41. 41
    Post-transcriptional cosuppression of Ty1 retrotransposition
    1. DJ Garfinkel
    2. K Nyswaner
    3. J Wang
    4. JY Cho
    (2003)
    Genetics 165:83–99.
  42. 42
  43. 43
  44. 44
  45. 45
  46. 46
  47. 47
  48. 48
    Exploring network structure, dynamics, and function using NetworkX
    1. A Hagberg
    2. P Swart
    3. S Chult D
    (2008)
    Proceedings of the 7th Python in Science Conference (SciPy2008). Presented at the SciPy. pp. 11–15.
  49. 49
  50. 50
  51. 51
  52. 52
  53. 53
  54. 54
  55. 55
  56. 56
  57. 57
    Hybrid dysgenesis in Drosophila melanogaster: a syndrome of aberrant traits including mutation, sterility and male recombination
    1. MG Kidwell
    2. JF Kidwell
    3. JA Sved
    (1977)
    Genetics 86:813–833.
  58. 58
  59. 59
    Jupyter notebooks – a publishing format for reproducible computational workflows
    1. T Kluyver
    2. B Ragan-Kelley
    3. F Pérez
    4. B Granger
    5. M Bussonnier
    6. J Frederic
    7. K Kelley
    8. J Hamrick
    9. J Grout
    10. S Corlay
    11. P Ivanov
    12. D Avila
    13. S Abdalla
    14. C Willing
    15. JD Team
    (2016)
    Proceedings of the 20th International Conference on Electronic Publishing. Presented at the 20th International Conference on Electronic Publishing. pp. 87–90.
  60. 60
  61. 61
  62. 62
  63. 63
  64. 64
  65. 65
  66. 66
  67. 67
  68. 68
  69. 69
  70. 70
  71. 71
  72. 72
  73. 73
  74. 74
  75. 75
  76. 76
  77. 77
  78. 78
  79. 79
  80. 80
  81. 81
  82. 82
  83. 83
  84. 84
  85. 85
    Data structures for statistical computing in Python
    1. W McKinney
    2. S van der Walt
    3. J Millman
    (2010)
    Proceedings of the 9th Python in Science Conference. pp. 56–61.
  86. 86
  87. 87
  88. 88
  89. 89
  90. 90
  91. 91
  92. 92
    The population genetics of speciation: the evolution of hybrid incompatibilities
    1. HA Orr
    (1995)
    Genetics 139:1805–1813.
  93. 93
  94. 94
  95. 95
    Genetic analysis reveals that FLO11 upregulation and cell polarization independently regulate invasive growth in Saccharomyces cerevisiae
    1. SP Palecek
    2. AS Parikh
    3. SJ Kron
    (2000)
    Genetics 156:1005–1023.
  96. 96
  97. 97
  98. 98
  99. 99
    Recombination of ty elements in yeast can be induced by a double-strand break
    1. A Parket
    2. O Inbar
    3. M Kupiec
    (1995)
    Genetics 140:67–77.
  100. 100
  101. 101
    Non-mendelian female sterility in Drosophila melanogaster hereditary transmission of I factor
    1. G Picard
    (1976)
    Genetics 83:107–123.
  102. 102
  103. 103
  104. 104
    R: A Language and Environment for Statistical Computing
    1. R Development Core Team
    (2017)
    R Foundation for Statistical Computing, Vienna, Austria.
  105. 105
  106. 106
  107. 107
  108. 108
  109. 109
  110. 110
  111. 111
  112. 112
    Statsmodels Econometric and statistical modeling with python
    1. S Seabold
    2. J Perktold
    (2010)
    Proceedings of the 9th Python in Science Conference. Scipy. P.
  113. 113
  114. 114
  115. 115
  116. 116
  117. 117
  118. 118
  119. 119
  120. 120
    GNU Parallel: The Command-Line Power Tool
    1. O Tange
    (2011)
    Login: The USENIX Magazine.
  121. 121
  122. 122
  123. 123
  124. 124
  125. 125
  126. 126
  127. 127
  128. 128
    A new evolutionary law
    1. L Van Valen
    (1973)
    Evolutionary Theory 1:1–30.
  129. 129
  130. 130
  131. 131
  132. 132
  133. 133
  134. 134
  135. 135
    Fitness effects of Ty transposition in Saccharomyces cerevisiae
    1. CM Wilke
    2. J Adams
    (1992)
    Genetics 131:31–42.
  136. 136
  137. 137
  138. 138
  139. 139

Decision letter

  1. Kevin J Verstrepen
    Reviewing Editor; VIB-KU Leuven Center for Microbiology, Belgium
  2. Patricia J Wittkopp
    Senior Editor; University of Michigan, United States
  3. Douda Bensasson
    Reviewer

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

Acceptance summary:

We appreciate that your study furthers our insight into the factors that control transposon activity. It seems that in contrast to what has been suggested previously, hybrids do not always show an elevated transposon activity, and that the genetic background (and complex interactions between alleles/mutations) may be a more important influence on transposon activity.

Decision letter after peer review:

Thank you for submitting your article "The effect of hybridization on transposable element accumulation in an undomesticated fungal species" for consideration by eLife. Your article has been reviewed by three peer reviewers, one of whom is a member of our Board of Reviewing Editors, and the evaluation has been overseen by Patricia Wittkopp as the Senior Editor. The following individual involved in review of your submission has agreed to reveal their identity: Douda Bensasson (Reviewer #3).

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

We would like to draw your attention to changes in our revision policy that we have made in response to COVID-19 (https://elifesciences.org/articles/57162). Specifically, when editors judge that a submitted work as a whole belongs in eLife but that some conclusions require a modest amount of additional new data, as they do with your paper, we are asking that the manuscript be revised to either limit claims to those supported by data in hand, or to explicitly state that the relevant conclusions require additional supporting data.

Our expectation is that the authors will eventually carry out the additional experiments and report on how they affect the relevant conclusions either in a preprint on bioRxiv or medRxiv, or if appropriate, as a Research Advance in eLife, either of which would be linked to the original paper.

All reviewers appreciate the relevance and scientific rigor of your study. Still, all reviewers also agree that, despite the large body of work presented here, it may not be possible to conclude that hybridization never induces Ty activity (or reduces repression of Ty activity), as this may depend on the exact nature (evolutionary distance, specific alleleles…) of the parental genomes. While you are very careful with interpreting the results, we would ask you to more explicitly discuss these possible alternative scenarios.

Two reviewers also suggested using established reporter assays to directly measure Ty activity in hybrids. In principle, these are very easy experiments that would, in our opinion, provide a fantastic way to consolidate the conclusions. However, we want to leave it up to you to decide whether you want to include this in the current manuscript.

Reviewer #1:

In this paper, Hénault and colleagues investigate the hypothesis that transposons (Ty elements) may become more active in hybrids. They first explore the genomes of natural hybrid S. paradoxus lineages and do not find evidence for elevated Ty activity in hybrid lineages, although some strains exhibit signs of higher Ty activity than others. Next, a mutation-accumulation experiment where Ty activity is measured in 6 lineages of hybrids with varying degrees of heterozygosity further confirms these findings.

Overall, this is a well-written and interesting study that helps to gain further insight into the factors that control transposon activity. The results suggest that, in contrast to what has been suggested previously, at least not all hybrids show an elevated transposon activity, and that the genetic background (and complex interactions between alleles/mutations) may be a more important factor. That said, I am a bit puzzled by the way the problem was approached – I appreciate that forensics of natural genomes is perhaps the best way to investigate a natural genetic phenomenon, but it obviously also has shortcomings. The lab experiments are a great addition, but there I cannot help but wondering whether using reporter systems for Ty activity (e.g. the typical systems with fluorescent reporters that have an intron) would have not been a more efficient way to measure Ty activity in thousands of different (artificial and natural) hybrids? Perhaps substance for a follow-up study…

Specific questions and concerns:

I wonder to what extent the genome analysis (i.e. the first part of the paper) depends on how long ago the hybridization events occurred. It seems plausible that if Ty elements become more active, this effect is the biggest immediately after the hybridization, and disappears over time as the genome re-establishes suppression of Ty activity. Would this influence the findings? It seems useful to discuss this in the text.

I suspect that using third-position codon diversity is a rather poor estimate of an effective population size, especially over longer periods of time? If this is the case, perhaps the results using this analysis should be interpreted with even more caution?

Is it necessarily true that the overrepresentation of ORF-disrupting Ty insertions in the SpC lineage points at a bias towards more deleterious insertions? Are the numbers strong enough to exclude noise? Or could selection be lower/more recent here?

Is it necessarily true that we expect a negative correlation between Ty insertions near essential genes and a lower purifying selection? Could noise play a role? And, how strong is the effect that a mutation (or more accurately, a Ty insertion) near an essential gene is on average more deleterious? Do we know how far that effect stretches from the gene, and was this taken into account in the analysis?

Reviewer #2:

I think the topic is very interesting, the work is well done, and the manuscript is well written (maybe can be shorten). Overall, the message is that the hybridisation does not trigger TE expansion in yeast and is supported by population genomics and experimental evolution data. I have two main points:

1) Did the authors were too conservative with their hybrid design and crossed populations that were too similar? I guess their rational is to mimic what they isolated in nature, but maybe if they look at the S. cerevisiae x S. paradoxus crosses that they published before, this might change.

2) I think the natural strains described in the first part of the paper and used for the crossing, would benefit for characterising TE activity, e.g. looking expression or better using existing reporters that enable to calculate rates of retrotransposition and excision. I think this can have major impact in the results obtained from the MALs and is a more direct read out of the hybridisation effect.

Reviewer #3:

This work uses wild yeast as a model and thus integrates many approaches to test competing theories about transposable element (TE) contributions to early speciation. It shows that natural hybrids between S. paradoxus subspecies show TE copy numbers (CNs) that are intermediate between those of their parents. Population genetic structure explains most of the natural variation among strains in TE levels though there could be a small effect of climate on Ty3p CNs. There is no evidence that TEs are at higher levels because of random drift. Mutation accumulation experiments with good replication similarly show no increase in transposition in hybrids compared with non-hybrid controls. Past work showed that the model yeast lineages used here are partially reproductively isolated. This work shows that their hybridization does not lead to increased TE reactivation and is not driving their speciation through Bateson-Dobzhanski-Muller interactions.

The experiments and analyses are convincing, and I have no substantive concerns. Furthermore, this number of integrated analyses is needed and would only be possible in few eukaryotes. The authors present a thorough discussion that interprets findings in the context of the model system and literature for all eukaryotes, from molecular biology to evolutionary genetics.

https://doi.org/10.7554/eLife.60474.sa1

Author response

Reviewer #1:

[…] That said, I am a bit puzzled by the way the problem was approached – I appreciate that forensics of natural genomes is perhaps the best way to investigate a natural genetic phenomenon, but it obviously also has shortcomings. The lab experiments are a great addition, but there I cannot help but wondering whether using reporter systems for Ty activity (e.g. the typical systems with fluorescent reporters that have an intron) would have not been a more efficient way to measure Ty activity in thousands of different (artificial and natural) hybrids? Perhaps substance for a follow-up study…

We thank the reviewer for the positive feedback and for the constructive suggestion of using retrotransposition reporter assays as an orthogonal approach. To avoid repeating our response, we invite the reviewer to see the response to reviewer 2’s second major comment, which expressed basically the same suggestion.

Specific questions and concerns:

I wonder to what extent the genome analysis (i.e. the first part of the paper) depends on how long ago the hybridization events occurred. It seems plausible that if Ty elements become more active, this effect is the biggest immediately after the hybridization, and disappears over time as the genome re-establishes suppression of Ty activity. Would this influence the findings? It seems useful to discuss this in the text.

We agree with the reviewer that the timing of hybridization may have an impact on the ability to observe Ty reactivation. It is entirely possible that a reactivation immediately following hybridization may have been rapidly suppressed by re-establishment of host repressive mechanisms, and that most resulting insertions were further purged by selection or lost by drift. Thus, this approach has limited ability to detect past, transient episodes of reactivation. We are aware that this is an important drawback of the population genomic approach, as highlighted by the reviewer. Thus, we modified the text to emphasize that our ability to detect reactivation from contemporary genomes is limited to persistent reactivation or to the maintenance of the accumulated TEs.

When first presenting the CN profiles from wild lineages, we added:

“The hybrid lineages SpC*, SpD1 and SpD2 had intermediate full-length Ty1 and Tsu4 CNs compared to their parental lineages, indicating no apparent persistent reactivation or maintenance of Ty elements in these lineages.”

When transitioning from the population genomics section to the experimental evolution section, we emphasized the difficulty to detect transient reactivation from contemporary genome analysis by adding this sentence: “Transient reactivation may also be followed by rapid re-establishment of host repression and leave little to no signature in contemporary genomes.”

Finally, we added a paragraph in our Discussion to mention this limitation and discuss the likelihood of transient reactivation scenarios in the light of our results and the known repression mechanisms in yeast:

“We note that the analysis of contemporary genomes is inherently limited in the ability to uncover transient reactivation that would immediately follow hybridization and rapidly faint, perhaps as a result of a re-establishment of TE repression. […] Additionally, Ty CNs of artificial hybrids at the initial time point of our MA experiment (i.e. early after hybridization) generally show good agreement with the expectation from their parental haploid strains, which is incompatible with an early reactivation of Ty elements.”

I suspect that using third-position codon diversity is a rather poor estimate of an effective population size, especially over longer periods of time? If this is the case, perhaps the results using this analysis should be interpreted with even more caution?

We agree with the reviewer that non-synonymous diversity at codon third positions only provides very rough estimates of Ne. While we tried to state this in the text, we understand that readers may easily overestimate the level of confidence we put in this result. We thus changed the main text to emphasize the suggestive nature of this result. We also cited Tsai et al., 2008 as an example of the use of nucleotide diversity in population genomic data to estimate population genetic parameters.

“Genome-wide nucleotide diversity at four-fold degenerate codon positions (πS) is related to long-term Ne and was used to approximate this parameter across diverse taxa (Lynch and Conery, 2003). […] While these rough Ne estimates should be interpreted with caution, this result suggested that natural selection efficiency is likely not the main determinant in CN variation.”

Is it necessarily true that the overrepresentation of ORF-disrupting Ty insertions in the SpC lineage points at a bias towards more deleterious insertions? Are the numbers strong enough to exclude noise? Or could selection be lower/more recent here?

Although many of our results tend to suggest that the SpC lineage could be subject to relaxed selection compared to others, we want to emphasize that these results constitute mixed evidence and we do not believe that the higher count of ORF-disrupting insertions in SpC compared to SpB (3 vs 1) is strong evidence in support for this hypothesis. We cannot exclude that it could be due to chance or sampling bias. The results we obtained in the search for a phenotypic defect associated with the private Ty1 insertion in ARR3 indicate that this polymorphism is probably neutral, since the single strain bearing it has no significant growth defect in arsenite compared to the rest of the SpC lineage. We conclude from this result that relaxed selection cannot explain the occurrence of this rare polymorphism and speculate that this interpretation could apply to other ORF-disrupting insertions. We thank the reviewer for highlighting that this was not easily understood from our previous formulation, and we modified the text to make it clear that we do not claim any selective interpretation about the ORF-disrupting insertions.

We edited the beginning of the first paragraph of this section:

“The number of ORF-disrupting Ty insertion alleles in SpC (three) compared to SpB (one) suggested a slight overrepresentation in SpC, which could be explained by a bias towards more deleterious insertion alleles in this lineage. To test this hypothesis, we examined the functional annotations of the corresponding genes.”

To strengthen our estimation of the phenotypic effect of the private Ty1 insertion in ARR3, we replicated the phenotypic measurement in arsenite at a larger scale. We modified the fourth paragraph of this section as follows:

“We investigated the fitness effect of the ORF-disrupting Ty1 insertion in ARR3 by measuring the growth of a collection of 123 S. paradoxus strains in culture medium supplemented with arsenite (Figure 5A, Figure 5—figure supplement 7, Figure 5—figure supplement 8). […] The absence of growth defect in arsenite for LL2011_004, combined with the extremely low frequencies of all the ORF-disrupting alleles, suggests that these insertions are effectively neutral and are not observed because of an important reduction in purifying selection efficiency.”

Figure 5 was modified. The previous Figure 5A-B was modified and added as a preliminary result as Figure 5—figure supplement 7. Replicates for all phenotypic measurements in NaAsO2 were added as Figure 5—figure supplement 8. Statistical analysis of phenotypic variation within the SpC lineage was added as Figure 5—figure supplement 9:

We added this section to the Materials and methods:

“Growth assay in NaAsO2

From the 208 wild strains used for whole genome sequence analyses, 123 strains had glycerol stocks available. [...] Post-hoc Tukey’s HSD test was performed using scikit-posthocs.”

Is it necessarily true that we expect a negative correlation between Ty insertions near essential genes and a lower purifying selection? Could noise play a role? And, how strong is the effect that a mutation (or more accurately, a Ty insertion) near an essential gene is on average more deleterious? Do we know how far that effect stretches from the gene, and was this taken into account in the analysis?

We agree with the reviewer that our explanation of the rationale we put forward to use adjacency to essential genes as a measure of selection efficiency lacks precision and nuance. Although it is noisy and surely has limited sensibility, this approach is expected to capture large variation in natural selection efficiency acting on the regulation of essential gene expression. Thus, we moderate our claims on the interpretability of these patterns and mention the expected sources of noise, which are indeed important.

Essential genes are known to be dosage-sensitive. Since by definition loss of function is lethal for these genes, important reduction in expression levels is expected to be deleterious for this class of genes. Additionally, they are overrepresented among genes that tolerate little overexpression (Makanae et al., 2013). They also have less genetic variation in mRNA levels (eQTLs) than most genes (Albert et al., 2018), showing that their expression levels are likely constrained by stabilizing selection. If natural selection is less efficient in a given lineage, we expect mutations that affect gene expression levels to be more frequently observed for this subset of genes compared to other lineages.

It is relatively well established that Ty elements can affect neighboring gene expression levels (Morillon, Springer and Lesage, 2000; Coney and Roeder, 1988; Todeschini et al. 2005; Lesage and Todeschini, 2005), either by increasing or decreasing them. Not all Ty insertions affect the expression of neighboring genes, as it depends on characteristics of the individual Ty insertions, notably position and orientation relative to the gene’s regulatory sequences. Although the study of many Ty1 insertions in S. cerevisiae allowed to define rules regarding the effect of position and orientation on neighboring gene expression (Lesage and Todeschini, 2005), to our knowledge, the generality of these rules at the genome scale were not assessed. Moreover, it was shown that single nucleotide substitutions can have dramatic consequences on the repressing or activating ability of a Ty element (Coney and Roeder, 1988), adding complexity to the prediction of gene expression effects caused by Ty insertions.

To address these aspects, we moved the background elements on the expression effects on adjacent genes from the second paragraph to the fifth paragraph of the section, and complemented it as follows:

“To gain a more comprehensive view of the strength of selection acting on Ty insertions, we compared the properties of genes in their neighborhood. […] Nevertheless, if selection efficiency explains a significant part of Ty CN variation, we reason that when CNs are high, Ty insertions should be more frequently observed near genes which properties reflect stronger dosage and functional constraints, compared to when CNs are low.”

In an attempt to add robustness to our analysis, we analyzed the distributions of number of protein-protein interactions (i.e. degree in the interactome) for genes neighboring Ty insertions, which corroborates our results based on gene essentiality, ω and k. We added the following sentences to the fifth paragraph of this section:

“Additionally, genes neighboring Ty insertions private to SpC show no bias towards higher number of physical interactions (Figure 5D, p>0.05, pairwise two-tailed Mann-Whitney U tests, FDR correction).”

We modified Figure 5 by adding a panel showing the distributions of counts of physical interactions, as detailed in the previous comment about ORF-disrupting Ty insertions.

We added the analysis of physical interaction data in the Materials and methods:

“Physical interaction data for S. cerevisiae S288c were retrieved from the BioGRID Multi-Validated Physical interaction dataset version 3.5.188 (Oughtred et al., 2019) and the interactome was built using Networkx.”

We analyzed the distributions of physical distance and relative orientation between Ty insertions and neighboring genes and show that differences among lineages are not significantly different. We added one sentence to the fifth paragraph of this section:

“Distributions of physical distance to neighboring genes were similar across lineages, as were distributions of relative orientation (Figure 5—figure supplement 10).”

We added Figure 5—figure supplement 10 to present the results of this analysis.

Former Figure 5—figure supplement 7 was changed to Figure 5—figure supplement 11.

Reviewer #2:

[…]

1) Did the authors were too conservative with their hybrid design and crossed populations that were too similar? I guess their rational is to mimic what they isolated in nature, but maybe if they look at the S. cerevisiae x S. paradoxus crosses that they published before, this might change.

We agree with the reviewer that restricting our analyses to the intraspecific divergence range excludes a vast range of higher divergence levels in interspecific crosses. The reason why we initially excluded the data from our S. paradoxus x S. cerevisiae crosses is that we judged it was less relevant to our natural system of S. paradoxus lineages (wild S. cerevisiae being largely absent at the latitudes where the analyzed strains were collected). While the choice to focus on intraspecific S. paradoxus crosses may seem conservative, we want to highlight their relevance in uncovering the mechanisms that contribute to early isolation in the process of speciation. Nevertheless, we agree that it makes sense to include our S. paradoxus x S. cerevisiae crosses in our analyses, which we named BSc in this manuscript. We updated the Introduction text to account for this new category of crosses:

“Here, we use population genomic data and laboratory evolution experiments on S. paradoxus and its sibling species S. cerevisiae to investigate the hybrid reactivation hypothesis and, more generally, the factors governing TE accumulation in natural and experimental lineages.”

In the first paragraph of this section, we changed the following:

“To systematically explore the effect of hybridization on Ty transposition rates, we performed a large-scale mutation accumulation (MA) experiment on 11 artificial hybrid backgrounds, eight of which were previously published (Charron et al., 2019).”

”Pairs of haploid derivatives of natural S. paradoxus and S. cerevisiae strains were chosen to span a range of evolutionary divergence (Figure 6B), yielding four categories of intraspecific crosses (from the least to the most divergent): SpC×SpC (CC), SpB×SpB (BB), SpB×SpC (BC) and SpB×SpA (BA); and one category of interspecific crosses: SpB×SpA (BSc; Figure 6C, Supplementary file 1A-F).”

We changed the third paragraph of this section as follows:

“Although these models explain little global variance (R2=0.40, 0.08 and 0.06 for Ty1, Ty3p and Tsu4 respectively), we found that neither independent variable was a significant predictor of ΔCN (Figure 6F). […] Crosses that shared the SpB parental strain UWOPS-91-202 (BB2, BC2, BA2, BSc2) seemed to have systematically higher Ty1 ΔCN values than the corresponding crosses that shared the SpB parental strain MSH-604 (BB1, BC1, BA1, BSc1).”

We modified the fourth paragraph of the Discussion as follows:

“This trend is especially important for Ty3p, as eight out of ten genotypes exhibit a negative change in median CN.”

We modified Figure 6 and Figure 6—figure supplement 1.

2) I think the natural strains described in the first part of the paper and used for the crossing, would benefit for characterising TE activity, e.g. looking expression or better using existing reporters that enable to calculate rates of retrotransposition and excision. I think this can have major impact in the results obtained from the MALs and is a more direct read out of the hybridisation effect.

We thank the reviewer for the suggestion of complementing our MA experiments with reporter assays that measure retrotransposition rate. We think this would be an entirely valuable approach to the question, as shown by Smukowski Heil and collaborators (Smukowski Heil et al., 2020), who recently demonstrated that no Ty reactivation occurs in interspecific hybrids between S. cerevisiae and S. uvarum.

As we understand the reviewer's comment, the occurrence of genomic changes during MA following a hypothetical initial and transient reactivation of retrotransposition might be the major factor making reporter assays a “more direct read out” of retrotransposition rates compared to MA. However, we would like to argue that in the time scale of our MA experiment (less than 800 mitotic generations), and given that we sequenced the strains at the initial timepoint of the MA experiment, it appears unlikely that signatures of an early and transient reactivation would have been purged from the genomes we sequenced. While we agree that many types of genomic changes can yield variation in Ty CNs that is unrelated to retrotransposition per se, such a “masking” of reactivation would require a strong mutational bias towards elimination of genomic loci containing Ty elements, or systematically very strong deleterious effects associated with novel Ty insertions that would drive their counter-selection in a MA regime characterized by minimal Ne. We cannot find any good reason to expect any of these scenarios. As we mention in the discussion of the manuscript, we believe that confident assessment of the Ty insertion loci in the MA lines would be essential in gaining further insight into the genomic processes involved in Ty CN change, but we believe it would require new data collection and thorough analyses that are beyond the scope of this manuscript.

We agree that reporter assays are an attractive and certainly powerful way of characterizing transposition rates. However, we believe that many limitations make them hard to apply to the current study in a straightforward and meaningful way. First, to our knowledge, all the retrotransposition marker constructs available are based on a Ty1 copy from S. cerevisiae. While this choice is certainly defendable in the case of interspecific hybrids involving S. cerevisiae as a parent (as in the case of Smukowski Heil and collaborators and our BSc crosses), most of our results are on intraspecific hybrids of S. paradoxus. The large sequence divergence level that separates these species is likely to have functional consequences on the ability of the marker to retrotranspose and impact retrotransposition rates in an unpredictable way, thus prompting the adaptation of these tools to S. paradoxus Ty elements. Moreover, MA experiments offer the advantage of exploring retrotransposition activity of all families at the time. Thus, exploring the dynamics of all families would ideally require generating constructs that account for both inter-lineage and intra-strain variation in Ty content. Second, while its interpretation might seem straightforward at first sight, retrotransposition marker assays require the function of many cellular processes that can exhibit phenotypic diversity at the intraspecific level and thus need to be carefully controlled for. In order to yield sufficiently high numbers of mutants to reach appreciable sensitivity in a reasonable experimental setup, most assays require transcriptional induction using the GAL1 promoter. In S. cerevisiae, it was shown that very low levels of intraspecific genetic variation results in wide phenotypic variation in the regulation of galactose genes (Lee et al., 2017). Additionally, splicing efficiency of the artificial introns inserted into HIS3 or GFP can also exhibit intraspecific variation. Ruling out potential effects from these confounding processes requires many controls in the form of assays and constructs that are not readily available to us. Finally, our choice of experimental condition for the MA experiment (YPD) is incompatible with most reporter assays, which require induction in galactose and selection of a plasmid, either by auxotrophy complementation in synthetic medium or antibiotic resistance in complex medium. All these requirements are important deviations from the basic YPD medium used in the MA experiment. Thus, the results from reporter assays would be hard to interpret, as it would be difficult to disentangle the orthogonal experimental approach from environmental effects.

In conclusion, the difficulty to identify a clear drawback to our MA approach given our sequencing strategy, and the complexity of adapting reporter assays and their necessary controls to S. paradoxus, prompt us to argue against the inclusion of reporter assays as a revision to our manuscript. While this is certainly a valid and important orthogonal approach to our question, we believe that its complexity deserves to be addressed in a separate future study.

https://doi.org/10.7554/eLife.60474.sa2

Article and author information

Author details

  1. Mathieu Hénault

    1. Institut de Biologie Intégrative et des Systèmes (IBIS), Université Laval, Québec, Canada
    2. Département de biochimie, microbiologie et bioinformatique, Université Laval, Québec, Canada
    3. Quebec Network for Research on Protein Function, Engineering, and Applications (PROTEO), Université Laval, Québec, Canada
    4. Université Laval Big Data Research Center (BDRC_UL), Québec, Canada
    Contribution
    Conceptualization, Data curation, Software, Formal analysis, Funding acquisition, Investigation, Visualization, Methodology, Writing - original draft, Writing - review and editing
    For correspondence
    mathieu.henault.1@ulaval.ca
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-0760-7545
  2. Souhir Marsit

    1. Institut de Biologie Intégrative et des Systèmes (IBIS), Université Laval, Québec, Canada
    2. Département de biochimie, microbiologie et bioinformatique, Université Laval, Québec, Canada
    3. Quebec Network for Research on Protein Function, Engineering, and Applications (PROTEO), Université Laval, Québec, Canada
    4. Université Laval Big Data Research Center (BDRC_UL), Québec, Canada
    5. Département de biologie, Université Laval, Québec, Canada
    Contribution
    Funding acquisition, Investigation, Methodology, Writing - review and editing
    Competing interests
    No competing interests declared
  3. Guillaume Charron

    1. Institut de Biologie Intégrative et des Systèmes (IBIS), Université Laval, Québec, Canada
    2. Quebec Network for Research on Protein Function, Engineering, and Applications (PROTEO), Université Laval, Québec, Canada
    3. Université Laval Big Data Research Center (BDRC_UL), Québec, Canada
    4. Département de biologie, Université Laval, Québec, Canada
    Present address
    Centre de foresterie des Laurentides, Ressources naturelles Canada, Québec, Canada
    Contribution
    Funding acquisition, Investigation, Methodology, Writing - review and editing
    Competing interests
    No competing interests declared
  4. Christian R Landry

    1. Institut de Biologie Intégrative et des Systèmes (IBIS), Université Laval, Québec, Canada
    2. Département de biochimie, microbiologie et bioinformatique, Université Laval, Québec, Canada
    3. Quebec Network for Research on Protein Function, Engineering, and Applications (PROTEO), Université Laval, Québec, Canada
    4. Université Laval Big Data Research Center (BDRC_UL), Québec, Canada
    5. Département de biologie, Université Laval, Québec, Canada
    Contribution
    Conceptualization, Supervision, Funding acquisition, Methodology, Writing - review and editing
    Competing interests
    Reviewing editor, eLife
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-3028-6866

Funding

Natural Sciences and Engineering Research Council of Canada (RGPIN-2015-03755)

  • Christian R Landry

Natural Sciences and Engineering Research Council of Canada (NSERC Alexander Graham-Bell doctoral scholarship)

  • Mathieu Hénault
  • Guillaume Charron

Fonds de recherche du Québec – Nature et technologies (FRQNT doctoral scholarship)

  • Mathieu Hénault
  • Guillaume Charron

Fonds de Recherche du Québec - Santé (FRQS postdoctoral scholarship)

  • Souhir Marsit

Canada Research Chairs (Canada Research Chair in Evolutionary Cell and Systems Biology)

  • Christian R Landry

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank Hélène Martin and Anna Fijarczyk for contributing to sequencing data analyses and Philippe Després, Anne-Marie Dion-Côté, Marika Drouin, Alexandre Dubé, Clara Dubernard and Anna Fijarczyk for reading and providing useful comments on draft versions of this manuscript.

Senior Editor

  1. Patricia J Wittkopp, University of Michigan, United States

Reviewing Editor

  1. Kevin J Verstrepen, VIB-KU Leuven Center for Microbiology, Belgium

Reviewer

  1. Douda Bensasson

Publication history

  1. Received: June 27, 2020
  2. Accepted: September 21, 2020
  3. Accepted Manuscript published: September 21, 2020 (version 1)
  4. Version of Record published: October 23, 2020 (version 2)

Copyright

© 2020, Hénault et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 956
    Page views
  • 133
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Evolutionary Biology
    2. Genetics and Genomics
    Richard Benton et al.
    Short Report

    The insect chemosensory repertoires of Odorant Receptors (ORs) and Gustatory Receptors (GRs) together represent one of the largest families of ligand-gated ion channels. Previous analyses have identified homologous 'Gustatory Receptor-Like (GRL)' proteins across Animalia, but the evolutionary origin of this novel class of ion channels is unknown. We describe a survey of unicellular eukaryotic genomes for GRLs, identifying several candidates in fungi, protists and algae that contain many structural features characteristic of animal GRLs. The existence of these proteins in unicellular eukaryotes, together with ab initio protein structure predictions, provide evidence for homology between GRLs and a family of uncharacterized plant proteins containing the DUF3537 domain. Together, our analyses suggest an origin of this protein superfamily in the last common eukaryotic ancestor.

    1. Evolutionary Biology
    2. Genetics and Genomics
    Fan Han et al.
    Research Article

    Atlantic herring is widespread in North Atlantic and adjacent waters and is one of the most abundant vertebrates on earth. This species is well suited to explore genetic adaptation due to minute genetic differentiation at selectively neutral loci. Here we report hundreds of loci underlying ecological adaptation to different geographic areas and spawning conditions. Four of these represent megabase inversions confirmed by long read sequencing. The genetic architecture underlying ecological adaptation in herring deviates from expectation under a classical infinitesimal model for complex traits because of large shifts in allele frequencies at hundreds of loci under selection.