Modeling and simulation of recurrent phenotypic and genomic selections in plant breeding under the presence of epistasis

Mohsin Ali,Luyn Zhng,In DeLy,Vivi Arief,Mrk Dieters,Wolfgng H.Pfeiffer,Jinkng Wng,Huihui Li,*

aNational Key Facility for Crop Gene Resources and Genetic Improvement,Institute of Crop Sciences,Chinese Academy of Agricultural Sciences,Beijing 100081,China.

bSchool of Agriculture and Food Sciences,The University of Queensland,Brisbane,QLD 4072,Australia

cHarvestPlus Challenge Programme,c/o International Food Policy Research Institute(IFPRI),1201 Eye St,NW,Washington,DC 20005,USA

ABSTRACT

1.Introduction

Selection of genetically superior genotypes among the huge amount of recombinant and segregating progenies is an essential but complex procedure in plant breeding[1,2].Development of cultivars involves cyclic crossing and selection procedures over long periods of time.Traditionally,plant breeders relied on phenotypic selection(PS)to determine the genetic potential of individuals or families in the field and chose the best genotypes that simultaneously exhibited multiple desirable traits.PS serves as an efficient strategy to improve complex traits by continuously increasing the frequency of favourable alleles.However,PS requires extensive field experiments,resources,and selection over a number of breeding cycles[3].In some cases,plant breeders have to reduce the number of genotypes that are to be phenotyped in the field due to limited resources.Conventional and modern breeding techniques have pushed the annual genetic gain of wheat grain yield from～0.7%to～1.2%[4,5]and a selection plateau has yet to be reached.It is forecasted that the world population will increase by 50%by the middle of this century and will require a 70% increase in crop productivity.The current annual genetic yield gain in major food crops including wheat is insufficient to meet these predicted future demands[4].New methods and tools must be considered and integrated with the conventional breeding methods in order to speed up the genetic gain.

In the 1980s,development of molecular markers greatly facilitated our understanding of breeding targeted traits and provided great potential to improve selection efficiency.One major application of molecular markers is QTL mapping,which is used to identify genomic regions linked to major genes for targeted traits.The identified linkage information between markers and QTL can be used in selection,which is referred to as marker-assisted selection(MAS).MAS has been successfully used for gene introgression by selecting those individuals that have favourable alleles linked with monogenic or oligogenic traits,for example,for nematode resistance in soybean[6],and Fusarium head blight in wheat[7].However,MAS has some limitations when used for selecting polygenic traits(such as grain yield),which are controlled by many QTL with minor effects[8].To address this limitation,marker-assisted recurrent selection(MARS)has been used for selection of complex traits[9].In MARS,selection is initially based on phenotypic values and marker scores,followed by several cycles of selection based on marker scores alone[10].MARS could accelerate the recurrent selection procedure by saving several seasons of field phenotyping.To acquire a score for selection,MARS relies on ad hoc significance tests for the marker and QTL association.It requires a cut-off criteria for QTL exhibiting major effects and therefore may exclude QTL with minor effects.

Genomic selection(GS)is another marker-based selection method using genome-wide and densely distributed molecular markers to increase the efficiency of improving complex traits[11].There are no significance tests in GS,and all markers contribute to the prediction of phenotypic values.In this way both major and minor QTL for the complex traits are included[12].GS requires the use of a training population(TP)and one or a few breeding populations(BP).To implement GS,the TP provides both phenotypic and genotypic data to train or develop the statistical model and predict genomic estimated breeding values(GEBVs)for selection in the BP,which is only genotyped.High-throughput,cost-effective,and high-density genotyping platforms have made it possible to predict genotypic values based on marker effects[13,14].In the past 10 years,many prediction models have been proposed that differ from each other in their range of assumptions in estimating breeding values and by their computational complexity[15].Due to its robustness and simplicity,GBLUP has been extensively used in animal and plant breeding for prediction and selection.Previous simulation studies on maize have demonstrated that BLUP-based GS produces an 18% to 43% higher response to selection than MARS,across genetic models with different QTL numbers and levels of heritability[12].Different selection methods continue to make progress in improving response to selection.However,gene-to-phenotype architecture(e.g.,epistasis and gene-gene interaction)of complex quantitative trait influences the expected phenotypic performance of new individuals[16].According to the present literature,the presence of epistasis in the underlying genetics of a trait is expected to influence phenotypic performance of progeny.However,the majority of studies have focused on additive genetic variance and paid little attention to the epistasis in breeding procedures.This is because detection of epistasis effects requires heavy computation for pairwise testing of alleles.Large-scale field-testing is also impractical for assessing the effects of epistasis on selection responses due to the time and resource limitation.Thus,computer simulation provides a fast and affordable alternative.

Advances in computer modeling and simulation provide advantages compared to conventional plant breeding[17]and can help breeders make critical decisions in the design of their breeding programs[12,18].Several computer simulation software packages(e.g.,QU-GENE,AlphaSim,DeltaGen,and BreedingSchemeLanguage)are currently available to plant breeders to support decision-making for cultivar development programs,especially with the integration of MAS methods[19–21].Sun et al.[22]have discussed the importance and application of several computer simulation software programs to assist plant breeders in their critical decisions in developing new cultivars.QU-GENE is one genetics and breeding simulation platform that can evaluate different selection and breeding strategies using complex genetic modeling scenarios[23].The QU-GENE simulation tool has a two-stage architecture(Fig.1).In the first stage,the QU-GENE engine is used to define the genotype-by-environment(GE)system(i.e.,all necessary genetic and environmental information for the simulation experiment)and to generate an initial parental population or base germplasm.

Fig. 1 – Workflow of the QuMARS breeding simulation tool. Three input files are needed to run QuMARS. In the first stage, the *.QUG input file is defined. It contains information regarding breeding targeted traits, genetic architecture (i.e. linkage map, QTL effects, and epistatic networks) of traits. Two output files are generated after running the QU-GENE engine as follows:genotype-by-environment information (*.GES) and a starting population (*.POP) for next step of simulation experiments. In the second stage, one other input called breeding strategy (*.QMP) needs to be defined. It contains information regarding crosses,planting, harvest, and selection details. In this stage, the information in the GES file will be used to simulate phenotypes of breeding traits. Parents for crosses come from the POP file. The selection procedure executes as defined in the QMP file.

The second stage includes application modules to investigate,analyze or manipulate the parental population within the GE system defined by the engine.Application modules are used to evaluate the efficiency of breeding strategies and identify ways to optimize the specific breeding procedure[17,18,24,26].Currently,QU-GENE tools include QuLine for self-pollinating crops[17],QuHybrid for inbred line development and hybrid performance prediction[24],and QuLinePlus,which is an expansion of QuLine to breeding cross-pollinating crop species[18].QU-GENE tools are freely available and can be downloaded from https://sites.google.com/view/qu-gene/.QU-GENE tools have been used in major crops such as wheat and maize to make strategic breeding decisions,like comparison of selection strategies and design breeding[17,25,26],and tactical decisions,such as optimization of crossing and selection methods[2,27,28].

However,to our knowledge,existing simulation software platforms lack the ability and flexibility to simultaneously simulate plant breeding programs with recurrent selection procedures with respect to PS,MAS,and GS.Recently,another QU-GENE tool called QuMARS has been developed to simulate recurrent PS,MARS and GS[29].From a breeding perspective,QuMARS allows breeders to simulate various strategies in phenotypic and genomic recurrent selection,and therefore to evaluate the genetic gain,change in genetic variance,and change in gene frequency etc.,over one or several cycles of selection.To the best of our knowledge,this is the first such simulation study conducted using the QuMARS software.Our objectives for the present study were:(1)to introduce the QuMARS application module;and(2)to model and simulate phenotypic,marker-assisted and genomic selection in the context of recurrent selection breeding programs.

2.Materials and methods

2.1.Quantitative genetics and breeding simulation platform of QU-GENE

In QU-GENE,the genotype-to-phenotype simulation is based on an E(N:K)model[23],where E stands for number of different environment types;N stands for number of QTL for breeding targeted traits;and K stands for level of epistasis,such as digenic or trigenic interactions.QU-GENE incorporates GE interaction effects and epistasis networks into the basic genotype-to-phenotype model[30].The phenotypic value of a trait is modelled in QU-GENE by:

where gijrepresents the genotypic value of the ith individual in the jth mega-environment,andεijis the micro-environmental random effect.During the simulation,the allele constitution of any individual was known and therefore its genotypic value could be calculated from the E(N:K)model defined by QU-GENE engine.The error effect in phenotypic values was randomly assigned from a normal distribution with a mean of zero and a variance equal to the user-specified error variance or calculated from the user-specified heritability level.

2.2.Development of the QuMARS application module

QuMARS is one of the application modules based on the QUGENE software[29]and is written in Fortran 90/95(freely available from https://sites.google.com/view/qu-gene/Download-page).It was originally developed to simulate marker-assisted recurrent selection,but can also simulate PS and GS.In the era of molecular breeding,QuMARS can act as a decision-making platform by simulating and optimizing the integration of GS and MARS in an ongoing conventional breeding program.Three input files are needed to run QuMARS.Two are outcomes from the QU-GENE engine(i.e.,*.GES and*.POP;Fig.1).The*.GES file contains the required information that can be used to predict genotypic and phenotypic values of any individuals during simulation,and the*.POP file defines a genetic population which can be used as parents for making crosses.The third input file*.QMP defines the breeding strategy that will be simulated(Fig.1)[29].The number of breeding cycles and selection criteria are defined by users.In simulation,single crosses are first made between two parental lines to derive training populations(e.g.DH,F2,and so on).Phenotyping and genotyping can be conducted for later or earlier generations.For example,when the objective is to select pure lines as cultivars(e.g.,in wheat and beans),dominance and heterosis are less important.One advanced selfing generation can be used as the training population for both phenotyping and genotyping.For hybrid breeding(e.g.,in maize and sorghum),one early segregating generation should be used as the training population where the dominance and heterosis can be fitted in the prediction model.To improve general combining ability,in some cases phenotype may need to be based on a testcross rather than the performance of a single genetic line.It should be noted that phenotyping and genotyping can be conducted within the same generation or in two separate generations.For example,if F2is used as the training population,genotyping can be conducted for F2individuals,but phenotyping should be conducted using F2-derived F3or even F4families.These scenarios can be simulated using QuMARS.

2.3.Genetic models used in the simulation

We assumed a genome consisting of five chromosomes,each with 100 evenly distributed markers at 2-cM intervals(Table 1).We simulated traits that are controlled by simple to complex genetic models,represented by three levels of QTL number(i.e.,1,2,and 5 QTL per chromosome)and three levels of broad-sense heritability(i.e.,0.1,0.4,and 0.8).One additive and two epistasis models were considered.For each model,QTL effects were set as random in the*.QUG input file.During simulation,the QTL effects were randomly generated from uniform distribution.Linked QTL may have effects in the same direction,representing linkage in the coupling phase;or in opposite directions,representing linkage in the repulsion phase.

Epistasis models were only considered for those situations when there were either 2 or 5 QTL per chromosome.For EP1,the interaction was between QTL located on different chromosomes.For the case of 2 QTL per chromosome,the five epistasis networks were Q1:1×Q2:2,Q2:1×Q3:2,Q3:1×Q4:2,Q4:1×Q5:2,and Q5:1×Q1:2(the subscripts represent chromosome number and number of QTL per chromosome,respectively).For the case of 5 QTL per chromosome,the five epistasis networks were Q1:1×Q2:2×Q3:3×Q4:4×Q5:5,Q-2:1×Q3:2×Q4:3×Q5:4×Q1:5,Q3:1×Q4:2×Q5:3×Q1:4×Q2:5,Q4:1×Q5:2×Q1:3×Q2:4×Q3:5,and Q5:1×Q1:2×Q2:3×Q3:4×Q4:5.For EP2,the interaction was between QTL located on the same chromosome.For the case of 2 QTL per chromosome,the five epistasis networks were Q1:1×Q2:1,Q1:2×Q2:2,Q1:3×Q2:3,Q1:4×Q2:4,and Q1:5×Q2:5.For the case of 5 QTL per chromosome,the five epistasis networks were Q1:1×Q2:1×Q3:1-×Q4:1,Q1:2×Q2:2×Q3:2×Q4:2,Q1:3×Q2:3×Q3:3×Q4:3,Q1:4×Q2:4×Q3:4×Q4:4,and Q1:5×Q2:5×Q3:5×Q4:5.

2.4.Simulation of training and base populations

Two inbred parents were generated by the QU-GENE engine.In each genetic model,we assumed a two-locusmodel in which 1 and 2 were used to represent the two alleles at each locus.At each locus,one parent has allele 1 and other parent has allele 2.The breeding strategy to generate the training population was defined in the*.QMP file(Fig.2,Table S1).In the crossing block(CB),a single cross was made between the two inbred parents to generate 10 F1individuals.F1individuals and their selfed seeds were harvested in bulk.In the next generation,500 F2individuals were planted and genotyped;and the selfed seed was harvested individually to produce the next generation.The 500 F2:3families were planted and phenotyped.The top 10%F2:3families were selected based on their phenotypic performance for the trait of interest as defined in the*.QUG file and harvested in bulk.The 50 selected F2:3families were grown and randomly mated to generate the base population for recurrent selection(Fig.2,Table S1).During recurrent selection,the top 10%of individuals were selected by different methods(i.e.,PS,MARS,and GS).The selected 50 individuals were randomly mated to generate seed for the next cycle of recurrent selection(Fig.2,Table S1).The recurrent selection continued for a total of 15 cycles.

Table 1–Summary of marker,QTL,and selection information used to define the genetic model.

Fig.2–Flowchart of the selection methods to be simulated and compared:phenotypic,marker-assisted and genomic selections.A single cross was made between two inbred parents to generate 10 F1 individuals.From the F1 generation,500 F2 individuals were acquired,planted,and genotyped.In the next generation,500 F2:3 families were planted and phenotyped.At each cycle(C)of recurrent selection,the top 10%(50 individuals)families were selected by different methods(Table 1),randomly mated(RM)and grown to serve as the next generation.This procedure continues for 15 cycles.GP represents genotype-to-phenotype;RM represents random mating;P represents parent;and C represents recurrent selection cycle.

2.5.Genotype-to-phenotype prediction models implemented in QuMARS

A total of six linear regression or genotype-to-phenotype(GP)prediction models were implemented in QuMARS.Three GP models were implemented for MARS including stepwise regression(Rstep),regression by forward selection(Forward),and regression by backward selection.These models relied on ad hoc tests for significant markers.Regression by backward selection was excluded from the present simulation due to poor prediction ability in our simulation experiments.Three GP models were implemented for GS including genomic best linear unbiased predictor(GBLUP),ridge regression(Ridge),and regression by the Moore-Penrose general inverse(InverseMP),which considered all markers in the prediction model.For Ridge,the shrinkage penalty of lambda for the tuning factor was set at 0.001 to avoid the exclusion of QTL with small effects.For Rstep,the probabilities of entering(PIN)and removing markers were set at a common level of 0.05 and 0.10,respectively.For InverseMP,a score of 0.99 was used,so that the majority of genetic variation was included in the prediction model.

甘薯淀粉因其分子结构及直、支链淀粉含量与绿豆淀粉相比存在巨大差异，导致以甘薯淀粉为原料生产的粉条品质与绿豆粉条相差甚远。因此，在传统甘薯粉条生产中通过加入明矾［KAl(SO4)2·12H2O］来改善粉丝品质以提高消费者可接受性。但铝的过量摄入对人体健康存在潜在危害，如老年痴呆等疾病。为减少明矾对人体的危害，同时达到改善甘薯粉条食用品质的目的，寻找明矾替代物意义重大。

2.6.Design and outcomes of the simulation experiment

The three genetic models previously described were built into three input files(*.QUG;Fig.1).The QU-GENE engine ran on each input file and produced two output files,one with the required information to define the GE system(*.GES)and the other with the required information to define the parental population(*.POP).The GE systems had a single environment type,three levels of QTL per chromosome,three levels of heritability,and two levels of epistasis.One reference population comprised of 100 homozygous individuals with a gene frequency of 0.5 for all loci was used to convert the specified level of heritability into error variance,which was subsequently used to assign the random effect associated with the phenotypic value of a given individual during simulation.Five GP methods were considered;Rstep,Forward,GBLUP,Ridge,and InverseMP.For comparison,PS was also included as a control method.A total of six selection methods were defined in one single QuMARS input file(*.QMP).Each selection method was replicated for 15 cycles and 50 replications.Each replication differed in QTL effects,genotypes sampled,and phenotypic values.

Three of the QuMARS outcomes were used in this study:adjusted genetic value(*.FIT),genetic variance component(*.VAR),and Hamming distance(*.HAM).Adjusted genetic value(also called fitness or population mean;Fad)was defined as:

3.Results

3.1.Selection responses from PS,MARS and GS under the additive model

Under the ADD model, GBLUP-based GS consistently resulted in higher adjusted genetic values (abbreviated as genetic value or population mean hereafter) than Rstep, Forward, and PS, especially in early recurrent selection cycles (i.e., cycles 1–5) across different numbers of QTL per chromosome and heritability levels (Fig. 3). As expected, high levels of heritability, i.e. H2= 0.8, resulted in a higher population mean regardless of number of QTL per chromosome and selection method (Fig. 3 row-wise). For example, when there was one QTL per chromosome, the population mean from GBLUP increased from 2.3% to 8.5% (P＜0.01) when the heritability level increased from 0.1 to 0.8 across 1–15 recurrent cycles(Fig. 3a, c). For each heritability level, changes in number of QTL per chromosome had a big impact on the population mean (Fig. 3 column-wise) and on the number of cycles to reach the selection plateau. For example, for a heritability level of 0.8 and GBLUP, the population mean decreased by 9.4%–21.9% (P ＜ 0.001) when the number of QTL per chromosome increased from 1 to 5, across 1–15 recurrent cycles (Fig.3). In addition, three, four, and more than 15 cycles were needed to reach the selection plateau for 1, 2, and 5 QTL per chromosome, respectively (Fig. 3c, f, and i).

In contrast to the population mean,Hamming distance was more affected by the number QTL per chromosome rather than the level of heritability for each selection method(Fig.S1).As for population mean,higher levels of heritability resulted in faster achievement of the target genotype,i.e.,lower value of Hamming distance,regardless of the number of QTL per chromosome(Fig.S1 row-wise).Under low to moderate levels of heritability,most of the selection methods(except Ridge and InverseMP)were effective in achieving the target genotype in early cycles(Fig.S1 column-wise).As mentioned previously,greater numbers of QTL per chromosome needed more breeding cycles to achieve the target genotype,which was also confirmed by the results based on Hamming distance.PS took more breeding cycles to reach the minimum Hamming distance as compared with GBLUP and MARS(Fig.S1),indicating that PS may have some advantages in achieving long-term breeding objectives.

Selection of superior genotypes at each recurrent cycle resulted in decreases in additive variance(Fig.4)and total genetic variance(Fig.S2)for both short and long terms.For example,the high heritability level(0.8)resulted in the fastest decrease in total genetic variance for all cycles and selection methods regardless of number of QTL per chromosome(Fig.S2 row-wise).Similarly,the additive variance was also decreased for all selection methods except for PS,InverseMP,and Rstep at cycle 1(Fig.4a,d).

Fig. 3 – Average population means over 50 replications for 15 cycles of recurrent selection and six methods for the additive(ADD) model for three levels of heritability (0.1, 0.4 and 0.8, column-wise) and three values of QTL per chromosome (1, 2 and 5,row-wise).

3.2.Selection responses from PS,MARS and GS under epistasis models

Under the EP1 model,PS consistently resulted in a higher population mean than that from GBLUP and MARS regardless of number of QTL per chromosome,heritability level,and selection cycles,except for cycle 1(Fig.5).For example,when there were 5 QTL per chromosome and the level of heritability was 0.8,the population mean from PS significantly increased(P＜0.001)from 7.8 to 50%in 1–15 recurrent cycles.Furthermore,at cycle 1 when there were 2 QTL per chromosome and the level of heritability was 0.1,the population mean increased by 0.4%(P＜0.001)for GBLUP,whereas population mean significantly decreased(P＜0.001)by 0.4,1.4,2.1 and 4.9%for Ridge,Forward,Rstep and InverseMP,respectively,as compared with PS(Fig.5a).In addition,for a heritability level of 0.1,when the number of QTL per chromosome increased from 2 to 5,the population mean significantly decreased(P＜0.05)from 4.9%to 21.6%for PS,8.7% to 20.0% for GBLUP,and 6.9% to 11.7% for Rstep in 1–15 cycles(Fig.5,row-wise).However,a minor difference in population mean was observed between GS and MARS,except for cycle 1(Fig.5).

Hamming distance under the EP1 model was more affected by the number of QTL per chromosome and heritability level than when using the ADD model(Fig.S3),indicating that the target genotype was more difficult to achieve even after 15 cycles of selection,especially when using GS and MARS(Fig.S3).For 2 QTL per chromosome,the Hamming distance was more favorable when PS was used versus GS and MARS(Fig.S3a,b,c),but no advantage was observed for 5 QTL per chromosome(Fig.S3d,e,f).In contrast,the total genetic variance was higher for MARS(Ridge,Rstep,and Forward)for early cycles,except when the heritability level was 0.4 and 5 QTL were located on each chromosome(Fig.S4).For a heritability level of 0.8,a rapid decrease in total genetic variance was observed for PS in early cycles regardless of the number of QTL per chromosome(Fig.S4c,d,e,f).Low to moderate levels of heritability resulted in a decrease in total genetic variance for PS and InverseMP regardless of the number of QTL per chromosome.For 2 QTL per chromosome,additive variance increased in the early cycles and started decreasing in later cycles except for Ridge,for the three different heritability levels(Fig.6a,b,c).For 5 QTL per chromosome,additive variance increased for all selection methods in the early cycles(Fig.6d,e,f).

Fig.4–Average additive variances over 50 replications for 15 cycles of recurrent selection and six methods for the additive model(ADD)for three levels of heritability(0.1,0.4,and 0.8,column-wise)and three values of QTL per chromosome(1,2,and 5,row-wise).

2Major results from the EP2 model(Figs.S5,S6,S7,S8)were similar to those observed from EP1(Figs.5,6;Figs.S3,S4),but differences were also observed.When the heritability level was moderate to high,and 2 QTL per chromosome were assumed,the EP2 model resulted in a higher population mean as compared with the EP1 model for all selection methods,except for cycle 0 and high heritability levels(Figs.5,S5).Hamming distance results confirmed that the EP2 model resulted in higher genetic values under low to moderate levels of heritability and reached the target genotype faster than the EP1 model when 2 QTL per chromosomes were considered(Figs.S3,S6).PS resulted in higher a Hamming distance for EP2 as compared with EP1,except for cycle 1(Figs.S3,S6).Total genetic variance at cycle 0 showed differences between the two epistasis models(Figs.S4,S7),for all the numbers of QTL per chromosome and levels of heritability.Similarly,results from additive variance also indicated that EP1 and EP2 models differed greatly from each other at cycle 0.

4.Discussion

4.1.Factors affecting genetic gains in simulation

A major task for plant breeding is to improve the population mean and increase genetic gain.These parameters are influenced by population type and size,mating system,genetic architecture,heritability of the trait of interest,and＞selection methods,such as PS,MARS and GS[12,13,31,32].Previous reports indicated that complex genetic architecture(including number of QTL per chromosome,QTL effects,interand intra-locus gene interactions,pleiotropy,coupling,and repulsion linkage phases)of breeding targeted traits influenced phenotypic variation and selection response[16].Present simulations revealed that the population mean decreased as the number of QTL per chromosome increased for all selection methods used in the simulation.In addition,improvement of the population mean for the polygenic trait required more recurrent cycles to reach the selection plateau than for a simple trait.Heritability is another important factor that influences the population mean.Population means from all selection methods increased as the level of heritability increased,regardless of the number of QTL per chromosome,in accordance with previous simulation studies[12,33].Under the ADD model,GS resulted in a higher population mean compared with MARS and PS in early cycles(Fig.3a,d,g),suggesting that GS is more efficient,particularly for complex traits with low levels of heritability[12].In another simulation study,Muleta et al.[34]used the AlphaSimR package to compare genetic gains from genomic-assisted recurrent selection and phenotypic recurrent selection.Simulation results indicated that genomic-assisted recurrent selection caused an increase in the genetic gain when breeding polygenic traits with low levels of heritability,regardless of population size.Further,Muleta et al.[34]discussed that higher genetic gain in early cycles may be due to concomitant changes in the genetic architecture of the trait under selection because short-and long-term recurrent cycles are expected to segregate large and small effect loci,respectively.In contrast,the poor performance of PS for low levels of heritability was due to large random errors(environmental noise)associated with phenotype(Fig.3a,d,g).This problem could be addressed by using GS and MARS(Fig.3g).Hamming distance results provided further evidence that GS and MARS could reach the highest genotypic value faster than PS when the level of heritability is low(Fig.S1).

Fig.5–Average population means over 50 replications for 15 cycles of recurrent selection and six methods for the epistasis network of QTL located on different chromosomes(EP1,Table 1)for three levels of heritability(0.1,0.4,and 0.8,column-wise)and three values of QTL per chromosome(1,2,and 5,row-wise).

Messina et al.[35]have demonstrated that drought resistance in crops could be improved if epistasis was captured by the GP models.In this study,as expected,the presence of QTL interactions(EP1 and EP2)underlying the genetics of traits significantly affected the population mean,particularly for marker-based selection methods.This is because epistatic effects were ignored in prediction models for both GS and MARS,which led to a considerable loss in response to selection[36].Exclusion of epistasis effects in prediction models is a limitation of the current version of QuMARS software.We are working on considering epistatic effects and including more advanced prediction models based on the Bayesian algorithm,machine learning,and deep learning in the next version of QuMARS.

Fig. 6 – Average additive variances over 50 replications for 15 cycles of recurrent selection and six methods for the epistasis network of QTL located on different chromosomes (EP1, Table 1) for three levels of heritability (0.1, 0.4, and 0.8, column-wise)and three values of QTL per chromosome (1, 2, and 5, row-wise).

4.2.Enhancing total genetic and additive variances for improving traits through modeling

The efficiency of a selection method can be quantified by genetic variance and additive variance components in the breeding population.The results from the ADD model confirmed that selection of superior individuals in each recurrent cycle increased the population mean at the loss of genetic and additive variance,across all variations in number of QTL per chromosome and levels of heritability(Figs.4,S2).The loss of variance was due to directional selection of individuals[37],which ultimately fixed favourable alleles[37,38].More recently,Muleta et al.[34]conducted a simulation experiment using AlphaSimR and found that the decline in total genetic variance for phenotypic recurrent selection and genomic-assisted recurrent selection was faster in early selection years than that from later selection years,irrespective of genetic architecture,level of heritability of the trait,and population size,which is also in agreement with the results of this study(Fig.S2).Under the ADD model,increased additive variance was observed in cycle 1 for PS,Ridge and InverseMP(Fig.4a,c,d).This trend also was observed for Hamming distance.In cycle 0,only two parents existed.The population mean and additive variance from cycle 0 shown in Figs.3 and 4,etc.,were only used as the starting point of recurrent selection and may not be comparable with corresponding parameter values after selection(i.e.cycles 1–15).

Presence of epistasis QTL in the genetic architecture largely influenced total genetic and additive variances regardless of selection method[39],especially for low to moderate levels of heritability[26].Our results also confirmed that the presence of epistasis underlying the genetics of polygenic traits could increase or maintain additive genetic variance in the breeding population,regardless of selection method(Figs.6,S8)[40–43].It was also hypothesized that GS has the capability to deliver accurate predictions because the presence of epistatic QTL action could be converted into additive variance,especially when the additive by additive interaction is a major part of epistasis[40,44].Results from the present simulation provide some evidence to support such a hypothesis(Figs.6,S8).

4.3. Comparison of PS with other selection methods

Under the ADD model, GBLUP, Forward and Rstep had a higher population mean and slower decrease in genetic variance,which indicated that GS and MARS favoured the polygenic traits(QTL number = 5) with low levels of heritability for short-term selection (Fig. 3). In long-term selection, population means from PS were lower or equal to GS and MARS, but additive and genetic variance from PS was not the lowest.

GBLUP is useful for monogenic traits where a short-term response can be maximized (Fig. 3). In comparison, PS was more responsive to selection under the epistasis models (Figs. 5, S5).The lower response of marker-based selection under epistasis was due to inconsistent marker effects and allele frequencies.Furthermore, as mentioned previously, prediction models have excluded epistatic effects for simplicity and computational efficiency. However, in the future, if both additive and epistatic effects can be fitted in prediction models, selection efficiency from GS could be further improved.

4.4. Practical applications in breeding

Increasing response to selection (i.e., genetic gain) is the first priority of any breeding program. Breeders rely on various selection methods to enhance response to selection to improve target traits. Simulation results presented in this study suggested that understanding the genetic architecture of the breeding traits of interest is extremely helpful to determine an efficient breeding strategy. If the complex traits are controlled by many additive QTL and with levels of low heritability (e.g.,grain yield), GS would be more useful compared with PS and MARS; if they are controlled not only by additive QTL but also by inter- and -inter locus interactions, PS would be more useful regardless of cost. A few of selection cycles (4–6) are enough to achieve maximum response with modest population size (i.e.,500 individuals reflect small breeding programs in developing countries) for traits with low to moderate levels of heritability.Genetic variance components of the traits of interest should be taken into account to ensure enough genetic diversity for subsequent recurrent cycles. This step is usually neglected in practice. It is notable that there are various factors, such as biotic and abiotic factors, available resources, population size,selection intensity, etc., which may affect the overall response to selection, regardless of the selection methods.

Supplementary data for this article can be found online at https://doi.org/10.1016/j.cj.2020.04.002.

Declaration of competing interest

Authors declare that there are no conflicts of interest.

Acknowledgments

This work was financially supported by the National Key Research and Development Program of China(2015BAD02B01-2-2),and the HarvestPlus Challenge Program(www.harvestplus.org).

Author contributions

Huihui Li and Jiankang Wang developed the QuMARS software.Mohsin Ali conducted the simulation experiments and wrote the draft.Jiankang Wang and Huihui Li revised the manuscript.Luyan Zhang,Ian DeLacy,Vivi Arief,Mark Dieters,and Wolfgang H.Pfeiffer edited the manuscript.