Author response: Joint inference of evolutionary transitions to self-fertilization and demographic history using whole-genome sequences Article Swipe

MAKE IMPACT

Stefan Strütt , Thibaut Sellinger , Sylvain Glémin , Aurélien Tellier , Stefan Laurent · YOU? · · 2023 · Open Access · · DOI: https://doi.org/10.7554/elife.82384.sa2

Full text Figures and data Side by side Abstract Editor's evaluation Introduction Results Discussion Materials and methods Appendix 1 Appendix 2 Data availability References Decision letter Author response Article and author information Metrics Abstract The evolution from outcrossing to selfing occurred recently across the eukaryote tree of life in plants, animals, fungi, and algae. Despite short-term advantages, selfing is hypothetically an evolutionary dead-end reproductive strategy. The tippy distribution on phylogenies suggests that most selfing species are of recent origin. However, dating such transitions is challenging yet central for testing this hypothesis. We build on previous theories to disentangle the differential effect of past changes in selfing rate or from that of population size on recombination probability along the genome. This allowed us to develop two methods using full-genome polymorphisms to (1) test if a transition from outcrossing to selfing occurred and (2) infer its age. The teSMC and tsABC methods use a transition matrix summarizing the distribution of times to the most recent common ancestor along the genome to estimate changes in the ratio of population recombination and mutation rates overtime. First, we demonstrate that our methods distinguish between past changes in selfing rate and demographic history. Second, we assess the accuracy of our methods to infer transitions to selfing approximately up to 2.5Ne generations ago. Third, we demonstrate that our estimates are robust to the presence of purifying selection. Finally, as a proof of principle, we apply both methods to three Arabidopsis thaliana populations, revealing a transition to selfing approximately 600,000 years ago. Our methods pave the way for studying recent transitions to self-fertilization and better accounting for variation in mating systems in demographic inferences. Editor's evaluation This manuscript details an important development of population genetics theory that can be used to infer past changes in the selfing rate of natural populations. The inference procedure is convincing and represents a substantial improvement upon previous methods. The work will be of broad interest to researchers studying mating system evolution and its consequences and will improve demographic inferences drawn from population genetic approaches. https://doi.org/10.7554/eLife.82384.sa0 Decision letter Reviews on Sciety eLife's review process Introduction Hermaphroditism is common in many groups of eukaryotes, especially in plants, and allows uniparental reproduction through selfing. The rate of self-fertilization is known to vary widely between species and populations; from exclusive outcrossing, through mixed mating, to predominant self-fertilization (Stebbins, 1950; Whitehead et al., 2018). Transition from outcrossing to predominant self-fertilization is the most frequent reproductive transition in flowering plants (Barrett, 2010) and is thought to have occurred many hundreds of times. Such a transition has profound ecological and evolutionary consequences affecting the genetic functioning and the demography of a population as well as patterns of dispersal. Taking selfing rate and its possible variations into account is thus key to properly understand the evolution of fate of a species. In flowering plants, cross-fertilization is ensured by diverse molecular and morphological features, some of which being referred to as self-incompatibility (SI) mechanisms (Charlesworth, 2010), which are defined as ‘the inability of a hermaphroditic fertile seed plant to produce zygotes after self-pollination’ (de Nettancourt, 1997). SI mechanisms enable the pistil of a plant to identify and repel self-pollen or pollen of a related genetic type and as a consequence, avoid inbreeding (de Nettancourt, 1997). SI systems are known to be encoded by a small number of genes and can therefore easily be lost. Indeed, genetic disruptions of SI systems through naturally occurring mutations are thought to be a major driver of plant reproductive diversity and have been linked to recent transitions to predominant self-fertilization in several species (Shimizu and Tsuchimatsu, 2015; Mattila et al., 2020). On the short term, whether a new mutation responsible for an SI breakdown will invade the population or be lost depends on the balance between the main advantages of selfing (reproductive assurance and gene transmission advantage) and inbreeding depression, that is the reduced fitness caused by an increased homozygosity under inbreeding (Charlesworth and Charlesworth, 1987; Charlesworth, 2006). On the long term, selfing is predicted to increase extinction rate and to reduce diversification of selfing lineages as it has been observed in a few clades such as Solanaceae (Zenil-Ferguson et al., 2019), Primulaceae (de Vos et al., 2014); but see Landis et al., 2018. A consequence is the supposed recent origin of selfing species due to an excess of transitions on terminal branches. It thus nicely illustrates how the balance between micro-evolutionary and macro-evolutionary processes can generate the observed distribution of mating systems among species (Igic et al., 2008; Goldberg et al., 2010). This peculiar dynamic is also invoked to explain that the genomic effects of selfing are often detected on within-species polymorphism but very rarely on between-species divergence (Glémin et al., 2019). The recent timing of these transitions is thus a central assumption but has not been systematically tested as it remains a challenging task. Despite the increasing number of available methods to reconstruct the evolutionary history of populations using genomic data (Excoffier et al., 2013; Schiffels and Durbin, 2014; Boitard et al., 2016; Terhorst et al., 2017; Speidel et al., 2019), none of them is currently able to explicitly identify and estimate the age of transitions in reproductive strategies. In flowering plants, previous attempts to estimate the age of transitions to selfing were based on limited genetic variation at a single locus directly controlling the reproductive mode in plants: the S-locus (reviewed in Mattila et al., 2020). For example, naturally occurring loss-of-function alleles at the S-locus were shown to be responsible for the loss of SI in Arabidopsis thaliana (Tsuchimatsu et al., 2010). The steady accumulation of non-synonymous alleles following loss of constraint at the S-locus was used to estimate that the age of the transition is at most 1.48 million years old, based on current estimates of the mutation rate (Shimizu and Tsuchimatsu, 2015). Note that the original upper-bound estimate in Bechsgaard et al., 2006, is 413,000 years owing to the use of a different mutation rate (see Figure 4 in Shimizu and Tsuchimatsu, 2015). This approach is limited by the small number of genetic variants upon which the estimation is conducted and can only be used in species for which, as is the case in A. thaliana, the genetic determinism of the loss of SI is known. However, a shift in reproductive system is expected to strongly impact genome-wide polymorphisms patterns, that is not only at the (S-)loci controlling it, thereby leaving a potentially characteristic molecular signature. We use this rationale and build upon previous theoretical work to develop two inference tools allowing to use full-genome polymorphism data of any species in order to (1) reveal the occurrence of past changes in reproduction mode and (2) estimate their age. A classic way to consider selfing assumes a theoretical population of N diploid individuals which produce offspring through selfing or outcrossing with probability σ and 1 – σ, respectively. Under a model of neutral evolution, the distribution of polymorphic sites in a sample of sequenced individuals, that is the frequency of alleles (single nucleotide polymorphism [SNPs]), is determined by the underlying genealogy of this population. A genealogy has for properties its length measured as the time to the most recent common ancestor (TMRCA) and its topology, that is the order and number of branching processes. Along the genome, genealogies change due to the effect of recombination (the so-called ancestral recombination graph [ARG]; Hudson and Kaplan, 1988; Wiuf and Hein, 1999). Two population parameters determine the distribution and characteristics of genealogies observed in a sample of several genomes: the population mutation rate (θ) and the population recombination rate (ρ). In the presence of predominant selfing, the effective population size (Nσ) of a population of N individuals is scaled by the selfing rate (σ) yielding Nσ=N/(1+F) (Fu, 1997; Nordborg and Donnelly, 1997) and the recombination rate (rσ) is scaled as rσ=r(1 – F) (Golding and Strobeck, 1980; Nordborg, 2000; Padhukasahasram et al., 2008) at least when r is not too high (Padhukasahasram et al., 2008), where r is the recombination rate (for example per site) in the genome and F is the inbreeding coefficient. When inbreeding is only due to partial selfing, F=σ/(2 – σ) and takes values between 0 and 1 (0 for outcrossing and 1 for fully selfing). As a consequence, the population recombination rate takes now for value in a selfing population: (1) ρσ=4Nσrσ=4Nr(1−F)/(1+F) With μ being the mutation rate in the genome (e.g. per site), the population mutation parameter accounts for the effect of selfing as follows: (2) θσ=4Nσμ=4Nμ/(1+F) We note that the classic ratio of population recombination by population mutation rate ρ/θ=r/μ in outcrossing species becomes now with selfing ρσ/θσ=r(1 – F)/μ (Nordborg and Donnelly, 1997; Möhle, 1998; Nordborg, 2000). Taken together, expressions (Equations 1; 2) suggest that selfing does not amount to a simple change in effective population size (from N to Nσ), and that the reduction of the effective recombination rate is more severe than the reduction in effective population size. Following these insights, two key predictions can be derived. First, a characteristic and specific signal of selfing, in contrast to outcrossing, is expected to be present in SNP data due to the joint action of recombination along the genome (rate ρσ) and of the genealogical (coalescence) process (rate θσ). The first prediction underlies the previous development of a sequentially Markov coalescent method (eSMC) to estimate a fixed selfing rate using estimations of the ratio ρσ/θσ from pairs of genomes (Sellinger et al., 2020). Second, evolutionary changes in reproductive systems (transition from selfing to outcrossing or vice versa) are reflected in variations of the ratio ρσ/θσ in time and are identifiable and distinguishable from changes in population sizes alone. The latter suggests that a characteristic signature of the change in selfing rate in time should be observed in polymorphism data along the genome, if genetic variation can be summarized in a way that is informative about the joint effect of genetic drift and recombination. Indeed, it is desirable to (1) estimate changes in population size which occur independently of a transition, for example when species colonize new habitats as facilitated by selfing, and (2) disentangle the possible confounding effect of population-size changes on the estimation and detection of a transition. To our knowledge no statistical inference method exists that takes advantage of these theoretical predictions to jointly estimate temporal changes in selfing rates and in population sizes. Two types of model-based methods are used to draw inference of past demographic events using full-genome polymorphism data. First, the distributions of the time to the most recent common ancestor (TMRCA) along a pair of chromosomes sampled from a sexually reproducing population can be modeled and approximated assuming the sequentially Markov coalescent (SMC) (McVean and Cardin, 2005). SMC-based methods infer model parameters, such as demographic histories, from fully sequenced genomes and have been implemented in several statistical software used to estimate changes in effective population size through time under the assumption of a Wright-Fisher (WF) model. Li and Durbin, 2011, developed a pairwise SMC estimating the distribution of TMRCA over two sequences, and Schiffels and Durbin, 2014, extended the framework to consider multiple (more than two) sequences at a time, albeit only estimating the most recent pairwise TMRCA. These SMC-based methods are based on a hidden Markov model (HMM) and the forward-backward algorithm for estimating the ARG (i.e. the genealogies along the sequence) and a Baum-Welch algorithm to estimate the parameters of the model. In our previously developed SMC method, eSMC, we included the effect of constant-in-time seed banks and self-fertilization in order to estimate these parameters (i.e. dormancy and selfing rates) jointly with past population sizes (Sellinger et al., 2020). eSMC uses also pairs of sequences and the estimated transition matrix which summarizes the transition between two consecutive genealogies along the genomes (due to a recombination event), in contrast to the computationally intensive approach based on the actual series of coalescence times (Gattepaille et al., 2016). The SMC methods are characterized by their easy applicability to datasets as these do not require the user to specify a given demographic model. However, the SMC methods rely on an analytical expression for the transition probability between genealogies and do not provide a measure of uncertainty of the parameters or allow for hypothesis testing between different models/scenarios (though in principle the transition matrix can be used to do so, Palacios et al., 2019). Second, approximate Bayesian computation (ABC) is a computational approach to estimate posterior probabilities for models and parameters that is well suited for demographic modeling with many parameters and without any analytically derived likelihood function (Beaumont et al., 2002; Csilléry et al., 2010). Two advantages of the ABC method are that it allows to compare competing demographic hypotheses on the basis of Bayes factors and it does not require bootstrapping the data to generate measures of uncertainty for the inferred parameters. A critical aspect of ABC is that it requires a careful summarization of the genomic data into a set of summary statistics that carry information about the parameters of interest (Beaumont et al., 2002). This step is specially challenging when using genome-wide data (Boitard et al., 2016). With the aim to use the information on genetic drift and recombination present in full-genome polymorphism data, we develop therefore both an SMC and an ABC approach to infer past demographic events and the time of transition to selfing. We first provide analytical and simulation results explicating the consequences of a transition to selfing on genomic variation, thereby confirming our second prediction that temporal changes in selfing rate leave observable specific patterns in polymorphism data across the genome. Then, we analyze the statistical accuracy of the two newly developed methods to identify and estimate the age of transitions to selfing using a small number of sequenced genomes. Third, as a proof of principle, we apply these methods to estimate the age of the transition to selfing in A. thaliana, in which it has been documented and for which full-genome polymorphism data exist. These new methods are useful toolkits for dating and understanding the evolution of the selfing syndrome and shifts in breeding systems and reproduction modes. Results The consequences of a transition to selfing on patterns of genomic variation We consider a theoretical population of N diploid individuals (equivalent to 2N chromosomes), which produce offspring through selfing or outcrossing with probability σ and 1 – σ, respectively (following the notations in Nordborg, 2000). At some time (tσ) in the past, the previously outcrossing population (with σ=0) undergoes a transition to selfing and remains selfing until present (with a selfing rate σ>0). Independently of the change in selfing rate, the population size can change once from NANC (ancestral) to NPRES (present) at time tN (measured from the present). We implement the selfing model both in the forward WF framework, in which selfing can be simulated explicitly, and using the coalescent-with-selfing, in which selfing is modeled through a rescaling of the effective population size and of the recombination rate at tσ (Nordborg and Donnelly, 1997; Möhle, 1998; Nordborg, 2000) (see Materials and methods). In other words, selfing rate is not constant in time in our model. In the spirit of the SMC approaches, we obtain for our model a first analytical result extending the previous theoretical work to consider the transition between two consecutive genealogies along the genome due to a recombination event (here reduced to the coalescence time of a sample of size two). In contrast to the classic assumption of SMC methods (Li and Durbin, 2011; Schiffels and Durbin, 2014; Palacios et al., 2019) the ages of recombination events do not follow here a uniform distribution along the genealogy, but rather are functions of the selfing and recombination rates at each time point. As recombination and selfing rates are assumed non-constant through time, the probability for a recombination event to occur follows an inhomogeneous Poisson process along the coalescent tree. We thus compute the probability for an effective recombination event to occur in a sample of size two (p(rec|s)), conditioned to the current coalescence time s, and with the selfing and recombination rates at time k being, respectively, σk and rk. We find (using Equation 1 for recombination above): (3) p(rec|s)=(1−e−∫0s 2(1−σk)(2−σk)2rkdk) Note that this value depends in an inhomogeneous way on the current coalescent time, which requires to integrate along the branch length of the coalescent tree. In other words, the variation of population size through time does not affect the probability of a recombination to occur conditioned on the coalescence time, that is the variation of population size affects the coalescence rate only. In Equation 3 it is visible that the probability for a recombination event to occur and to modify the genealogy is affected by the selfing rate (σk). Thus, this result opens the possibility to jointly infer piecewise functions of selfing or recombination rates and population sizes through time simultaneously. In practice, a change in selfing rate can be detected when the span of the genealogy along the genome does not decrease homogeneously with increasing coalescence time as demonstrated in Figure 1—figure supplement 1 (based on Equation 3). A further formal description of the coalescence process with variable recombination and/or selfing rates through time in the context of the SMC is provided in Appendix 1, especially providing the required expressions for the transition and emission probabilities needed to compute the transition matrix. Following the spirit of SMC methods, we now turn to simulations to explicit the importance of Equation 3 for inference. We thereafter verify the second prediction - temporal changes in selfing rate leave observable specific patterns in polymorphism data across the genome - by quantifying the expected distributions of TMRCA-segments (hereafter TL-distribution) under our model. These segments are defined as successive and adjacent sets of nucleotidic positions sharing the same time to the most recent common ancestor (TMRCA) and are separated by the breakpoints of ancestral recombination events (McVean and Cardin, 2005). Segments are summarized by their lengths and TMRCA and When the selfing rate is constant the has a constant on segments are and by a number of recombination events to segments In the case of a transition from outcrossing to predominant selfing, the rate at which the segments with age is increased in the outcrossing to the selfing to a characteristic change in the of the joint distribution at This of the model is using our coalescent model rescaling and or by explicitly using forward WF simulations 1—figure supplement this specific genomic signature of a transition to selfing is not observed when the selfing rate is constant and only the population-size changes from to where be the effective population size of a population with selfing rate σ The simulations thus the results from Equation 3 and Figure 1—figure supplement Figure 1 with 3 see of a transition to selfing on the genealogies of simulated and distributions of ages in generations on a and lengths of TMRCA-segments on a in a selfing population with a change from NANC to NPRES population size. The population sizes were to to the rescaling of the effective population size by the selfing rates used in and distributions of ages (TMRCA) and lengths of TMRCA-segments in a population with a constant population size and a shift from outcrossing σ=0) to predominant selfing distribution along the genome of a of the TMRCA-segments simulated in The transition matrix of ages (TMRCA) between adjacent segments along the genome for the data simulated in A. This matrix summarizes the probabilities that the with a given age is by the of age The the transition probabilities The transition matrix of ages (TMRCA) between adjacent segments along the genome for the data simulated in The recombination rate for the simulations was set to The data was by (see Materials and methods). the same increase in in principle be for by a very ancestral population size (i.e. NANC – such a model also have a ancestral to the selfing and thus produce the same as a transition to selfing thus is not a confounding for a transition to predominant The is determined by the probability of recombination which, as in Equation is affected by the change in selfing (σ) when on a given under a of values for tσ the of the change in the between the age and length of TMRCA-segments on the age of the transition 1—figure supplement 3). Thus, this suggests that genome-wide polymorphism data information about shifts to selfing when the age of the transition well the distribution of TMRCA. We also the important that the segments that in the outcrossing than are along simulated chromosomes This effect can also be by the of transitions between different TMRCA for a number of successive and adjacent and In the case of a transition to selfing, the TMRCA transition matrix that segments with TMRCA than tσ are more to be by segments that are also than this a between successive TMRCA also exists when selfing is the of the effect is more in the case of a shift to selfing. The specific genomic of a transition to selfing and its to two temporal changes in effective recombination rate and population size us to two new statistical methods to estimate the age of transitions to selfing. methods to estimate the age of a transition to teSMC on Equation we eSMC into allowing the estimation of selfing or recombination rates through time, jointly with population size (see Materials and methods and Appendix In order to account for two are implemented for parameter (1) the in which each hidden has its rate, and (2) the mode in which teSMC estimates only three parameters (the current and ancestral and the transition time between both a constraint the number of inferred parameters and well suited for the of recent and shifts from outcrossing to predominant First, to demonstrate the theoretical accuracy of our model and inference method, we analyze its when sequences of TMRCA are given as (the age of coalescent for two This is the of teSMC (Sellinger et al., We data from a population a and a transition to selfing or a change in recombination In both the population size and the past selfing or recombination values are with high accuracy supplement Second, to understand the properties of we analyze its under a simple assuming a constant population size and a constant selfing value of given different amount of data. We compare the eSMC method, which estimate a constant rate of selfing in time with which estimates selfing through time supplement When selfing is known to be constant the value of this parameter is with high accuracy and with the amount of given data supplement When each hidden can have its selfing value constant selfing rate is but a amount of data is required to reduce the in the estimation supplement We now the statistical accuracy of teSMC on polymorphism data from genome of simulated under a model with constant population size with mutation and recombination rates of and with an change from outcrossing to predominant selfing at time tσ (see Materials and methods). the mode and the mode estimation well over the of tσ values Figure 2 with see of teSMC on simulated polymorphism data. of times of transition from outcrossing to predominant selfing using neutral The the value of tσ in of and the the values of tσ estimated by was using the mode and the mode of teSMC and per time point. Under constant population size. were with an change in population the the age of the change in population size. NANC to NPRES at generations or generations in the NANC to NPRES at generations or generations in the The inference process was times for each independently simulated data sizes estimated under the assumption of a constant selfing rate are than the value in the outcrossing and in the selfing which be for past population-size supplement 3). However, when teSMC is allowed to account for the change in selfing rate, population-size estimates to the We note that the increased in N in the recent selfing are caused by a number of available Finally, we the of teSMC to jointly estimate the age of a transition to predominant selfing and the time of a change in population size. To do this we use simulated data as with the of a single population-size reduction and or and In both our results that especially the mode inference method, is able to estimate the age of the shift to selfing, of the timing of the population-size change and the transition to selfing. in most the population sizes inferred by teSMC are to the simulated values supplement However, when the transition is recent and the present population size is the of the population-size estimates supplement We note that teSMC in the latter to the population a of data (Sellinger et al., These results demonstrate that transitions to predominant self-fertilization and more changes in recombination rate through time can be by teSMC and the estimations can

Type

peer-review

Language

en

Landing Page

https://doi.org/10.7554/elife.82384.sa2

OA Status

gold

Related Works

10

OpenAlex ID

https://openalex.org/W4382394009

Raw OpenAlex JSON

OpenAlex ID

https://openalex.org/W4382394009

Canonical identifier for this work in OpenAlex

DOI

https://doi.org/10.7554/elife.82384.sa2

Digital Object Identifier

Title

Author response: Joint inference of evolutionary transitions to self-fertilization and demographic history using whole-genome sequences

Work title

Type

peer-review

OpenAlex work type

Language

en

Primary language

Publication year

2023

Year of publication

Publication date

2023-04-27

Full publication date if available

Authors

Stefan Strütt, Thibaut Sellinger, Sylvain Glémin, Aurélien Tellier, Stefan Laurent

List of authors in order

Landing page

https://doi.org/10.7554/elife.82384.sa2

Publisher landing page

Open access

Yes

Whether a free full text is available

OA status

gold

Open access status per OpenAlex

OA URL

https://doi.org/10.7554/elife.82384.sa2

Direct OA link when available

Concepts

Inference, Joint (building), Evolutionary biology, Human fertilization, Biology, Computational biology, Computer science, Artificial intelligence, Genetics, Engineering, Architectural engineering

Top concepts (fields/topics) attached by OpenAlex

Cited by

0

Total citation count in OpenAlex

Related works (count)

10

Other works algorithmically related by OpenAlex

Full payload