Adaptations in wild radish (Raphanus raphanistrum) flowering time, Part 1: Individual-based modeling of a polygenic trait

Abstract

This study investigated replicating six generations of glasshouse-based flowering date selection in wild radish (Raphanus raphanistrum L.) using an adaptation of the population model SOMER (Spatial Orientated Modelling of Evolutionary Resistance). This individual-based model was chosen because it could be altered to contain varying numbers of genes, along with varying levels of environmental influence on the phenotype (namely the heritability). Accurate replication of six generations of genetic change that had occurred in a previous glasshouse-based selection was achieved, without intermediate adjustments. This study found that multiple copies of just two genes were required to reproduce the polygenic flowering time adaptations demonstrated in that previous research. The model included major effect type M1 genes, with linkage and crossing over, and minor effect type M2 genes undergoing independent assortment. Within the model, transmissibility (heritability of each gene type) was parameterized at 0.60 for the M1 genes and 0.45 for the M2 genes. The serviceable parameterization of the genetics of flowering in R. raphanistrum within a population model means that simulated examinations of the effects of external weed control on flowering time adaptations are now more feasible. An accurate and simplified Mendelian-based model replicating the adaptive shifts of flowering time that is controlled by a complex array of genes is useful in predicting life-cycle adaptations to evade weed control measures such as harvest weed seed control, which apply intense adaptive selections on traits that affect seed retention at harvest, including flowering time.

Introduction

Herbicide-resistant (HR) weeds pose the greatest threat to modern food production globally. Despite the adoption of multiple strategies for control of HR weeds, the rate of increase in new HR weed biotypes globally continues to grow (Delye et al. 2013; Gorddard et al. 1996; Heap 2023). To mitigate the selection of multiple HR weed biotypes, agricultural systems employ a range of non-herbicidal weed control tactics targeted at specific timings of the weed’s life cycle.

Evolutionary theory suggests that recurrent management practices will select for traits in weeds that improve survival and reproduction. Similar to herbicide resistance evolving as a consequence of herbicide overreliance, weeds also have a long history of amending their life cycles to evade nonchemical weed control practices (Barrett 1983). Life-cycle adaptations include temporal mimicry of the crop’s germination, growth, and seed set phenology (Barrett 1983), plus the physical mimicry of the crop’s seed color, shape, and size (Baker 1974; Vigueira et al. 2013). These adaptations enable weeds to more successfully grow and reproduce within a crop field.

In response to the ongoing evolution of multiple herbicide resistance in weed populations, harvest weed seed control (HWSC) was developed to capture or destroy weed seeds during crop harvest (Walsh et al. 2022). The various HWSC techniques target weeds that have survived herbicide applications earlier in the season and intercept seeds retained above the cutting height of the harvester before they enter the soil seedbank. Therefore, HWSC reduces the fecundity and seed dispersal of HR weed biotypes (Somerville et al. 2018). However, the long-term success of HWSC is dependent on the continuing opportunity to capture weed seeds at harvest over multiple generations (Norsworthy et al. 2016; Walsh et al. 2022).

To evade capture, weeds can adapt to flower earlier in the growing season, thereby increasing the potential for weed seed shedding before crop harvest. Weeds typically contain significant standing genetic variability in flowering time traits (Cleland et al. 2007; Simpson and Dean 2002). This variability is particularly evident in wild radish (Raphanus raphanistrum L.) populations (Ashworth et al. 2016). In this species, flowering time adaptation could confer an ability to shed seeds earlier and thereby evade seed capture at harvest. However, earlier flowering can come at a cost to the plant. Plants that flower earlier in the season have less time to grow vegetatively, making them likely to be less competitive against crops and produce fewer seeds (Ashworth et al. 2016).

Enhanced knowledge and understanding are needed to predict the potential effects of future polygenic changes in weed populations. Accordingly, the aim of this study was to build a genetic model to replicate induced changes to flowering time in R. raphanistrum. The modeling will simulate individual Mendelian-based genetics (with crossing over and environmental effects on expression), aiming to replicate six generations of glasshouse-based results reported in Ashworth et al. (2016), without intermediate adjustment.

Materials and Methods

Model Design

The SOMER (Spatial Orientated Modelling of Evolutionary Resistance) model adapted here had previously used genetic change to predict herbicide resistance (Somerville and Renton 2015). In this paper, polygenic change in days to first flower (DFF) is investigated. Therefore, the herein newly adapted model has been renamed SOMEF (Spatially Orientated Model of Evolutionary Flowering). This renaming is somewhat presumptive, as there is no spatial component used here, because we replicated a glasshouse study.

DFF

The glasshouse study of Ashworth et al. (2016) is used here as a template to parameterize DFF in R. raphanistrum. That detailed study was chosen because it contained records of data spanning a wide range of genetic change. In addition, procedures were followed in that study to minimize environmental influence on DFF, meaning that the observed phenotypic results may better reflect the underlying genotypes. Different combinations of genetics were trialled within SOMEF, with the aim of building a model capable of replicating (without intermediate adjustment) eight generations of phenotypic changes reported in Ashworth et al. (2016).

Flowering time–related studies often express results in relation to growing degree-days (GDD), which allows comparisons across different environments. The effect of naturally occurring changes in day length in the sequential, out of growing season studies in the glasshouse experiments (Ashworth et al. 2016) were accounted for in that study by using GDD to assess flowering date. However, the current study required the use of temporal reproductive data in order to fit the temporal design within SOMER (Somerville and Renton 2015), meaning that DFF was needed within SOMEF as the unit of measure. Original records from the glasshouse study were available, which included the data for days required for the initiation of flowering. Adjustments were needed to account for the different growing seasons used in that study, with parameterization based on the same climatic conditions as the glasshouse study. Fortunately, two datasets from southern Western Australia were enabled (see Supplementary Material) conversion of the glasshouse data to a DFF unit of measure that would replicate the weed R. raphanistrum germinating within a southern Australia autumn-sown crop. This relationship enabled the use of the GDD in Ashworth et al. (2016) to estimate DFF and to thereby predict fecundity of individuals throughout the temporal growing season (used within SOMEF).

The Genes of Days to Flowering

The aim here was to build a model with a manageable number of genes influencing flowering date. These modeled genes are arbitrarily labeled M1 and M2. The modeled genes were assumed to act in a semidominant way within SOMEF. Semidominance was chosen because it produces a more diverse pool of phenotypes, meaning that SOMEF would be better at replicating a polygenic characteristic that is, in reality, under the control of hundreds of genes (Nie et al. 2016).

M1 and M2 Genes

Due to the complex nature of the genetics behind DFF (Nie et al. 2016), it was impractical to attempt to replicate these genetics exactly within an individual-based model. Instead, the preexisting Mendelian-based genetic inheritance within the SOMER model (Somerville and Renton 2015) was used, which was expanded herein to contain a manageable number of duplicated genes. To encompass the large range of values in the original DFF data (Ashworth et al. 2016) while maintaining model simplicity, it was necessary to use three copies each of just two types of genes.

The simulation genes created were (1) linked type M1 genes that each had a major effect on the flowering date and (2) unlinked type M2 genes that had a minor effect on the flowering date. The large effect of the M1 genes was included to allow coverage of the range in DFF (28–135 d) seen in Ashworth et al. (2016), whereas the M2 genes were designed to fill the gaps in the phenotypic range generated by using the major effect M1 genes. Although the M2 genes behaved differently from the M1 genes during reproduction, the R code (R Core Team 2008) was written to ensure that the linkages between the M1 genes (those that did not cross over) were maintained during meiosis. In addition, haploid M2 genetic combinations created by each plant at meiosis are retained together in the same ovum (and in the same pollen grain) as haploid M1 genetic combinations from the same plant. This was incorporated into the modeling and meant that within SOMEF no single ovum could be pollinated with M1 genes from one plant and M2 genes from another plant. For reproductive purposes, the genotypes were kept distinct (e.g., AaAaAa ≠ AAAaaa during meiosis). Notwithstanding this distinction during breeding and assuming equal effects for the same value of M1 alleles, the values of the diploid genotypes on flowering dates could be expressed more simply by counting the number of semidominant M1 alleles (with similar treatment of the M2 alleles). For both types of genes, minimal DFF were assumed to be due to the absence of any semidominant (uppercase type) alleles.

Crossing Over

Genetic recombination within this modeling was assumed to be “crossing over,” in which genetic material moves from one chromosome to another during meiosis. Crossing over was investigated as a method of creating new genetic combinations that would be phenotypically new and somewhat stable in the new populations, while allowing for further genetic change under more vigorous selection. Four extra semidominant M1 alleles were assumed capable of being sequentially added via crossing over (two per chromosome) to generate plants with up to 10 (5 × 2) semidominant M1 alleles. Interactions between crossing over events and the number of recessive M1 alleles were ignored; any number or position of recessive M1 alleles was assumed to result in no change in the crossing over frequency or on the weed’s DFF.

Heritability of Genetic Traits

This study assumed that the mean phenotype equaled the mean genotype for each population, and that environmental effects were random and normally distributed. Heritability between each pair of temporally adjacent generations was calculated from the original glasshouse data, where narrow sense heritability (h ²) was calculated from average phenotypic measurements of the parent (P ₀) and offspring (P ₁) generations (Equation 1):

(1) $${h^2} = {{{\rm{mean}}({P_0}) - {\rm{mean}}({P_1})} \over {{\rm{mean}}({P_0}) - {\rm{mean}}({\rm{selection}}\;{\rm{from}}\;{P_0})}}$$

An incorporation of heritability was needed within SOMEF, as individual plant phenotypes (P _i) were used for selection. A new term, “transmissibility” (T ²), is introduced here. Transmissibility describes the degree of influence of change in each type of gene on the modeled plant’s phenotype. A reworking here of Equation 1 (using the abovementioned assumptions) yielded Equation 2. Transmissibility (T ²) was derived for the M1 and M2 genes by using Equation 2 and iteration.

(2) $${P_{\rm{i}}} = {\rm{mean}}({P_0})\; + {{[{G_{\rm{i}}} - {\rm{mean}}({P_0})]} \times {{T^2}\left( {1 + e} \right)\;}}$$

Equation 2 was needed so that each plant’s genotype (G _i) could be used to stochastically assign it a phenotype (P _i). These phenotypes were then used as the basis for the plant’s selection to reproduce. The best construction of Equation 2 was chosen as the one giving the best fit to the data set from Ashworth et al. (2016). The random environmental factor (e) in Equation 2 was generated from a standard normal distribution. The constant T ² was calculated via iteration for each gene (M1 and M2) to give the best fit to the data. This particular construction of Equation 2 incorporates an influence whereby larger intergenerational change increases the phenotypic variability.

Model Parameterization

This investigation looked at the number and strength of the major effect, linked M1 genes, and their crossing over frequency. The minor unlinked M2 genes were also investigated, as was the transmissibility (T ²) of both the M1 and M2 genes.

Biologically based values within the ranges identified for each variable were simultaneously trialed within concentric program loops. Each simulation began with a randomly selected G ₀ population (as defined in Ashworth et al. [2016]), and used R. raphanistrum life cycle, breeding, and selection parameters replicating conditions in the same original glasshouse study. Two hundred replications were averaged to reduce environmental variability in the results. Best fit to the original data was determined simultaneously and equally for early flowering date–selected EF1, EF3, EF4, and EF5 populations and long flowering time–selected LF1 and LF3 populations (as defined in Ashworth et al. [2016]), using least-square regression analysis in R (R Core Team 2008). The modeling implemented allowed no intergenerational adjustments to any variables.

Results and Discussion

The successful replication of six generations of glasshouse-based flowering date selection was achieved in an individual-based model using crossing over and Mendelian-based inheritance of only two types of genes. This simulation was achieved even though the true genetics controlling the initiation of flowering in R. raphanistrum is complex, controlled by more than 100 genes (Nie et al. 2016) and influenced by environmental cues (Ehrlén 2015). The use of just two types of semidominant genes in this model in no way replicates the complex genetics involved in flowering initiation. However, this study lends support to the concept that the exact reproduction of genetics in models is unnecessary to successfully model individual-based inheritance of polygenic characteristics.

The simulation of eight generations of DFF selection (six of which were matched with the environmental conditions and data from Ashworth et al. [2016]) was achieved, using two types of genes: linked M1 genes with high transmissibility and unlinked smaller-effect M2 genes with lower transmissibility (Figure 1). The specific parameterizations required included the number and effect size of the M1 and M2 genes, their transmissibility, and the crossing over probabilities of the various M1 gene combinations. Genetic combinations were chosen that gave the best fit to the data from the glasshouse study of all the options that were trialed. Two glasshouse-grown generations (EF2 and LF2) described in Ashworth et al. (2016) reversed the intergenerational shifts in DFF (Figure 1B) for reasons unknown, and therefore were not matched with modeled data in this study.

Figure 1. Goodness of fit of the modeled data (A) when compared with the recorded glasshouse data (B). Graphs show cumulative days to first flowering (DFF) adaptations in Raphanus raphanistrum populations as a result of repeated early and late days to flowering selection. In both graphs, the basal population is the black solid line in the center, the darker lines are the matched generations of early flowering (EF1, EF3, FE4, and EF5; colored orange), and late flowering (LF2 and LF3; colored purple), with and late flowering (LF2 and LF3), with the two unmatched generations (EF2 and LF2) shown in lighter tones. The far-left population (early flowering EF5) displays very little phenotypic variability, whereas the far-right population (late flowering 3) is very diverse.

DFF and Fecundity

Goodness of conversion from GDD data in the glasshouse study (Ashworth et al. 2016) to southern Western Australia DFF data was assumed, as genetically similar populations grown at different times of the year had similar DFF after independent mathematical conversions to autumn germination (results not shown). Two original Western Australian data sets from Cheam (1986) and Taghizadeh et al. (2012) were scaled so that they could be merged (Supplementary Figure S1). A quasi-Poisson function gave the best fit to the adjusted the best fit to the adjusted curve, yielding Equation 3 (Akaike information criteria, AIC=2.98) using nonlinear least-squares multiple regression:

(3) $${{{\rm{Seeds}}} \over {{\rm{Plant}}}} = {\rm{exp}}\left( { - 5.117 + 0.299{\rm{DFF}} - 0.002{\rm{DF}}{{\rm{F}}^2}} \right)$$

Equation 3 was used to calculate fitness (seeds per plant), based on DFF for each R. raphanistrum plant in the simulation. Random selection was then applied in the model to seeds from the phenotypically selected plants to reflect the seed-selection criteria used in the glasshouse study (Ashworth et al. 2016).

The Genes of Days to Flowering

A satisfactory fit to the four generations of early-flowering data (populations EF1, EF3, EF4, and EF5) in the glasshouse study (Ashworth et al. 2016) was possible using three, four, five, or six pairs of linked M1 genes (along with the M2 genes described later; results not shown). However, the wide temporal spread of the late-selected DFF generations in the long flowering day glasshouse study (populations LF1 and LF3) could not be duplicated using the same linked Mendelian-based genetics for any number and combination of these genes when starting from the same initial G ₀ population (results for the extensive range of non-fitting modeling are not shown). For greatest simplicity, three pairs of diploid M1 genes (initially adding 6 d to DFF per semidominant allele) were utilized in the next part of this study. The second set of genes, the M2 genes, were also fit using three pairs, with each semidominant M2 allele causing a 3-d delay in DFF initiation.

Further expansion to the genetic modeling within SOMEF was needed to successfully replicate the late-selected DFF data (populations LF1 and LF3) in the glasshouse study (Ashworth et al. 2016). Two changes to the behavior of the M1 genes were instigated. First, there was a need to increase the impact of each additional M1 allele on a sliding scale from 6 to 15 d per allele (Figure 2). Second, the M1 genes, which were initially linked (all on the same chromosome), were given the ability to undergo crossing over, resulting in an increase in the maximum number of M1 alleles on any chromosome. Crossing over events between paired chromosomes with similar M1 allelic strength were given an increased probability (Table 1).

Figure 2. The relationship between the number of semidominant larger M1 alleles and days to first flower (DFF), needed to model the long-day selections. Note that six smaller-effect M2 genes will add another 0 to 18 days to DFF. If each M1 allele added the same amount in the long-day selections, the relationship would be linear.

Table 1. This table can be read similar to a Punnett square and shows the results of chromosomal crossing over of the linked larger-effect M1 alleles. Rates were determined using concentric program loops to achieve best fit^a.

^a The diploid parental genes are expressed as counts of M1 alleles on their two haploid chromosomes, one chromosome presented vertically and the other horizontally. The gamete alleles within the table are counts of M1 alleles in their haploid form after chromosomal crossing over. For example, in the shaded squares, the parental genotype “1,3” (e.g., Aaa and AAA) can generate the two reproductive haploid chromosomes “2,2” (e.g., AaA and AAa) at a comparative rate of 0.259. These comparative rates were reduced by a factor of 10 for use in the simulations.

In addition, a small degree of “central tendency” was added to the model so that the crossing over events that increased the percentage of plants with the median DFF were slightly more likely to occur, compared with crossing over events that reduced the percentage of plants with median DFF (analogous to stabilizing selection [Donohue 2002]). The crossing over probabilities introduced here, for creating new linkages between M1 genes, were initially based on the number of possible different forms of each genotype. However, both these initial recombination probabilities, and their multiplicative factor were adjusted to enable results to better match the data in the glasshouse study, thereby yielding the crossing-over probabilities in Table 1. This expanded genetic variability, incorporating a sliding scale and crossing over, was able to replicate all six generations in the glasshouse study without intermediate adjustment.

Over the short- and long-day selections, the number of M1 alleles per plant varied over the available range, from 0 up to 10, which resulted in a large adaptative change in flowering times. In contrast, it was desirable for the M2 alleles to be retained as long as possible within the short-day selections, as they smoothed gaps, rather than contributing to significant adaptive change. Lower heritability in the M2 alleles was determined (using regression) and resulted in less selection for the M2 alleles. The lower elimination rate of M2 alleles meant that only three pairs were needed in the initial population, therefore simplifying the individual-based modeling of these Mendelian-type inheritances.

Within the glasshouse study (Ashworth et al. 2016), there was an exceptionally high level of phenotypic change following just five generations of selection, resulting in flowering time life-cycle changes that were not present in the initial population of 1,000 plants (Figure 1B). Although the shorter-day selections were easily duplicated, modeling a long delay in flowering time using just two types of genes was only successful with the implementation of genetic crossing over, plus a sliding scale of effect for the predominant M1 genes. Physical limitations exist for very early flowering; however, large diversity in flowering can be found in R. raphanistrum, including cool-season biennial growth (Singh et al. 2020). This biennial tendency (MBA, personal observation, 2014) may be a factor in the requirement for the predominant M1, genes to exhibit crossing over and a sliding scale of effect. The effect of the sliding scale was to give a multiplicative rather than an additive effect on flower initiation date. This sliding scale and multiplicative effect may be symptomatic of the long-day selections within the newly designed SOMEF model and cannot be taken as indicative of any true genetic changes in the glasshouse-based R. raphanistrum populations.

Heritability

Equation 1 was used to calculate heritability values for each selection within the glasshouse study (Ashworth et al. 2016). The four early-flowering (EF1, EF3, EF4, and EF5) generations replicated here had h ² values of 0.41, 1.10, 0.76, and 0.92, respectively, with the two late-flowering generations (LF1 and LF3) having h ² values of 0.92 and 1.45, respectively (in the glasshouse study). The local success of Equation 2 was confirmed using sequential applications of Equations 1 and 2 (with these h ² values replacing T ²) to intergenerational data (and assumed genotypes) based on data from the glasshouse study (results not shown).

Within the SOMEF simulations, the closest fit was achieved with consistent (stochastically applied) transmissibility values of 0.60 for the more strongly linked M1 genes and 0.45 for the independently assorting weaker M2 genes. Semidominant M2 alleles were still present in simulated EF5 populations, even after five generations of selection for early flowering.

The gene type transmissibility values giving best fit to the data (0.6 and 0.45) were within the range of heritability (h ²) calculated for the 2014 to 2016 glasshouse study data (Ashworth et al. 2016) and also those published in other flowering time studies (Burgess et al. 2007; Conner and Via 1993; Mazer and Schick 1991) In this study, the modeled data fit well using consistent transmissibility (T ²) values across the eight generations of selection, despite the typical real-world intergenerational variability in the h ² values in Ashworth et al. (2016).

The random environmental effect on phenotypic diversity that best fit the glasshouse data (Ashworth et al. 2016) was simulated using a new equation (Equation 2), which, as well as introducing a new variable (T ²), also had an unusual configuration. Within Equation 2, a link was created so that larger intergenerational genetic changes resulted in greater phenotypic diversity. This resulted in increased phenotypic similarity (or a central tendency) in populations exhibiting reduced intergenerational genetic change and increased phenotypic diversity when intergenerational genetic change was high. Increased phenotypic diversity (i.e., increased differences from the genotypes) in highly selected populations (exposed to stressful selection criteria) may have an advantageous biological basis. Further investigations into a link between larger environmental effects (on phenotype) and stress-induced epigenetic changes (Henderson and Jacobsen 2007) may be warranted.

This study clearly demonstrates that a complex, genetically controlled phenotypic trait such as flowering time can be significantly simplified and modeled to replicate measured data (from Ashworth et al. [2016]). This modeling included a new derived constant: “transmissibility,” useful in modeling polygenic change. Transmissibility is designed to moderate the effect of predicted genotypes on each plant’s phenotype, so that phenotypic characteristics can be used for selection. Weed population responses can be modeled by using the modeling methodology outlined in this study, weed population responses can be modeled to predict the rate and scale of phenotypic change in the field. For example, genetic change in response to the repetitive use of novel weed control tactics such as HWSC (Walsh et al. 2022). Predictions such as those derived from using the SOMEF model can help optimize the application of integrated weed management programs, to ensure important weed control techniques remain effective.