• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    Genomic prediction of Fusarium head blight resistance in early stages using advanced breeding lines in hard winter wheat

    2022-12-02 01:01:00JinfengZhngHrsimreepGillNvreetBrrJyotirmoyHlerShuktAliXiotinLiuAmyBernroPulStAmnGuihuBiUpinerGillBrentTurnipseeSunishSehgl
    The Crop Journal 2022年6期

    Jinfeng Zhng,Hrsimreep S.Gill,Nvreet K.Brr,Jyotirmoy Hler,Shukt Ali,Xiotin Liu,Amy Bernro,Pul St.Amn,Guihu Bi,Upiner S.Gill,Brent Turnipsee,Sunish K.Sehgl,*

    a Department of Agronomy,Horticulture & Plant Science,South Dakota State University,Brookings,SD 57007,USA

    b School of Computing,Queen’s University,Kingston,ON K7L 3N6,Canada

    c USDA-ARS,Hard Winter Wheat Genetics Research Unit,Manhattan,KS 66506,USA

    d Department of Plant Pathology,North Dakota State University,Fargo,ND 58108,USA

    Keywords:FHB GBS Genomic selection Multi-trait models Winter wheat Wheat scab

    ABSTRACT Fusarium head blight(FHB),also known as scab,is a devastating fungal disease of wheat that causes significant losses in grain yield and quality.Quantitative inheritance and cumbersome phenotyping make FHB resistance a challenging trait for direct selection in wheat breeding.Genomic selection to predict FHB resistance traits has shown promise in several studies.Here,we used univariate and multivariate genomic prediction models to evaluate the prediction accuracy(PA)for different FHB traits using 476 elite and advanced breeding lines developed by South Dakota State University hard winter wheat breeding program.These breeding lines were assessed for FHB disease index(DIS),and percentage of Fusarium damaged kernels(FDK)in three FHB nurseries in 2018,2019,and 2020(TP18,TP19,and TP20)and were evaluated as training populations(TP)for genomic prediction(GP)of FHB traits.We observed a moderate PA using univariate models for DIS(0.39 and 0.35)and FDK(0.35 and 0.37)using TP19 and TP20,respectively,while slightly higher PA was observed(0.41 for DIS and 0.38 for FDK)when TP19 and TP20(TP19+20)were combined to leverage the advantage of a large training population.Although GP with multivariate approach including plant height and days to heading as covariates did not significantly improve PA for DIS and FDK over univariate models,PA for DON increased by 20% using DIS,FDK,DTH as covariates using multi-trait model in 2020.Finally,we used TP19,TP20,and TP19+20 in forward prediction to calculate genomic-estimated breeding values(GEBVs)for DIS and FDK in preliminary breeding lines at an early stage of the breeding program.We observed moderate PA of up to 0.59 for DIS and 0.54 for FDK,demonstrating the promise in genomic prediction for FHB resistance in earlier stages using advanced lines.Our results suggest GP for expensive FHB traits like DON and FDK can facilitate the rejection of highly susceptible materials at an early stage in a breeding program.

    1.Introduction

    Several fungal pathogens continuously constrain global wheat production and food security.Fusarium head blight(FHB),also known as scab,is a devastating fungal disease of wheat that causes significant losses in grain yield and quality[1].FHB is expanding its horizons throughout major wheat-producing areas due to climate change,an increased wheat acreage under no-till cultivation,and adoption of maize-wheat rotations[2,3].Though several Fusarium species can cause FHB,Fusarium graminearum is the prominent pathogen for FHB in the United States,Canada,China,and some European countries[4].In addition to grain yield losses,FHB results in reduced quality due to contamination by mycotoxins,primarily deoxynivalenol(DON),posing serious health consequences to humans and animals if ingested in certain quantities[4,5].In 2014,the revenue losses for hard wheat due to FHB were estimated to be around $600 million in the Great Plains region of the US[6].

    Even though fungicides are frequently used to reduce FHB damage,the utilization of resistant varieties is considered the most effective and economical approach to combat diseases like FHB[7,8].Being quantitative in nature,resistance to FHB is governed by multiple quantitative trait loci(QTL)and highly influenced by changing environments.Resistance to FHB has been categorized as type I(resistance to initial infection),type II(resistance to spread of FHB symptoms within a spike),or type III(low mycotoxin accumulation)[9].However,type II resistance is more stable and utilized in breeding programs as compared to type I and III FHB resistance[7,9].Conventional QTL mapping and genome-wide association study(GWAS)have been used to dissect the genetic basis of FHB resistance and a large number of QTL have been identified across all the wheat chromosomes,including seven named QTL,Fhb1 to Fhb7[10-13].However,only a handful of QTL have a major effect on type II resistance(resistance to FHB symptom spread in a spike)and have been effectively utilized in wheat breeding globally,in particular Fhb1[9,14].In the US hard winter wheat region,most of the variation in FHB resistance are from native sources including cultivars‘Everest’,‘Overland’,‘Lyman’,and‘Expedition’.Furthermore,phenotypic selection for FHB resistance in both field and greenhouse is very complicated and some of the measurements such as DON content can only be obtained after harvest and are costly.Thus,genomic selection(GS)could be a promising approach to improve FHB resistance in wheat with reduced phenotyping efforts and costs.

    GS is an approach that employs linkage disequilibrium(LD)and estimates the genetic worth of an individual using genomewide markers[15,16].GS addresses the primary limitation of QTL mapping and marker-assisted selection by using a joint estimate of all marker effects[17].Thus,GS is useful for predicting and selecting complex traits controlled by several minor QTL that are difficult to map using QTL mapping[18].The genomic prediction(GP)models are developed by using genotypic and phenotypic data in a training population(TP)to predict the genomicestimated breeding value(GEBV)of individuals in the breeding population(BP)[16].GS has shown immense potential in plant breeding,and several studies have reported successful implementation of these strategies in different crops in recent years[19-21].GS is particularly useful for these traits where phenotyping is cumbersome or costly[22,23].

    Predictive ability of the GS model refers to the correlation between estimated GEBVs and the actual phenotypic values of the individuals in the validation set.The PA mainly depends on heritability of the traits,TP nature and size,and choice and optimization of the statistical models[23,24].Several studies have evaluated the GP models for predicting different FHB traits including disease index(DIS),fusarium damaged kernels(FDK),and DON in wheat.Most of these studies employed a cross-validation approach to evaluate the PA of GP in spring and soft winter wheat using various prediction models and strategies[25-30].However,the focus of these studies was limited to improving the PA within the TP being evaluated,rather than validating the improved models in forward prediction.Unlike other traits,such as yield,only a few studies reported the implementation of GP to select preliminary breeding lines for FHB resistance[31,32].So far,the potential of GS in improving the FHB resistance in early generations of a hard winter wheat breeding program has not been examined.

    Most of the previous studies compared univariate GP approaches,including ridge-regression best linear unbiased prediction(rrBLUP),and genomic best linear unbiased prediction(GBLUP),LASSO,Random Forest(RF),and several Bayesian approaches[25,26,29].In most cases,the performance of these GP models varied with FHB traits and cross-validation schemes used in the analyses.On the other hand,multi-trait(MT)approaches are used to improve the prediction ability for a primary trait,when secondary traits genetically correlated to the primary trait are available[33].MT models are of particular importance when the primary trait is difficult or expensive to phenotype and has low heritability.Plant height(PH)and days to heading(DTH)are often associated with FHB resistance in bread wheat and durum wheat[30,32].Thus,several studies have used PH and DTH as covariates to predict FHB resistance[14,30,32,34],and only two studies suggested improvement in PA using PH or DTH as covariates[30,32].Furthermore,different FHB traits,including DIS,FDK,and DON are known to have moderate positive genetic correlations[35].Recent studies have evaluated the use of different combinations of FHB traits along with DTH and PH as secondary traits to utilize this genetic correlation for improving prediction ability of MT models[30].For instance,FHB traits like DIS and FDK can be used as secondary traits for improving the prediction of DON,which is cumbersome and costly to phenotype.Thus,there is a need to evaluate the usefulness of these covariates or to figure out the best combinations of secondary traits in MT models to predict FHB traits,especially in hard winter wheat.

    Another aspect of the successful application of GS in a breeding program is establishing a training population.Previous studies have evaluated strategies to optimize training populations using advanced breeding lines to predict GEBVs of preliminary breeding lines for various traits in wheat,including FHB[28,31,36,37].The application of GS in this scenario can be handy as it is challenging to phenotype a large set of preliminary lines in expensive FHB nurseries.

    The primary objective of this study was to use different sets of advanced breeding lines as training populations to predict FHB traits in preliminary breeding lines of our hard winter wheat breeding program.Secondly,we wanted to evaluate the performance of MT models to predict FHB traits,including DIS,FDK,and DON,using different combinations of secondary traits.For this,we evaluated FHB nurseries comprising advanced breeding lines from three years for their usability as training populations to predict these traits.Thus,our specific objectives for this study were to(a)evaluate the usability of FHB nurseries comprising advanced breeding lines for predicting FHB traits and assessing the improvement in predictive ability if the best-performing TPs from individual years were combined based on the lines shared among the FHB nurseries over years,(b)compare predictive abilities of MT models when different combinations of secondary traits are used to predict DIS,FDK,and DON,and(c)to validate the selected models and TPs in a forward prediction scheme to calculate the GEBVs for FHB traits in an independent breeding population comprising the preliminary breeding lines.

    2.Materials and methods

    2.1.Plant materials

    A panel of 476 wheat breeding lines from the winter wheat breeding program at South Dakota State University was used in this study.A set of breeding lines that is tested every year in the advanced yield trial(AYT)and the elite yield trial(EYT)was evaluated for FHB resistance in a mist irrigated field nursery.Among the 476 breeding lines,153 were evaluated in the 2018 nursery,169 in the 2019 nursery,and 154 in the 2020 nursery to optimize a training population for predicting preliminary breeding lines.Further,65 lines were shared between the 2018 and 2019 nurseries;and 58 lines were shared between the 2019 and 2020 nurseries.The majority of the breeding lines were either F4:7or F4:8filial generation.Lines without genotypic data and lacking consistency between replications were excluded and the final analysis was conducted on 152,161,and 153 lines from 2018,2019,and 2020 nurseries,respectively.In addition,a set of 200 breeding lines from the preliminary yield trial(PYT)was used as a breeding population for the prediction of FHB indices(independent validation)using the optimized training population.A random set of 60 lines was selected from the BP and evaluated in the 2020 nursery to validate the prediction accuracy(PA)in forward selection.

    2.2.Experimental design and studied characters

    Plant materials were planted in the FHB nurseries at Brookings,South Dakota(44.3114°N,96.7984°W)during the 2018,2019,and 2020 growing seasons(Table S1)using a randomized complete block design with 2 or 3 replicates for different sets of lines with corresponding checks.We used cultivars‘Lyman’and‘Emerson’as resistant checks while‘Overley’and‘Flourish’were used as susceptible checks.The experimental unit was a single-row plot(40 plants per meter row)for each line.FHB resistance and related traits were evaluated in 2018,2019,and 2020.DTH was recorded using Julian date when 50% of the main tillers in the row had completely emerged heads.PH was measured from the soil surface to the top of main tiller spikes excluding awns when the plant materials matured.

    All FHB nurseries were artificially inoculated using both cornspawn and spraying a spore suspension of F.graminearum isolates(SD-FG1)as described in Halder et al.[13].Briefly,Fusariuminfested corn kernels were scattered on soil surface once at the boot stage(Feekes 10)and again at heading stage(Feekes 10.1)to enhance the chances of infection in the field.At anthesis,a conidial suspension containing 100,000 spores/ml was sprayed on the heads of each line at 50%anthesis to avoid any escape.The nursery was misted using sprinklers every night(19:00-7:00)to maintain the humidity for disease development.Disease incidence and severity was recorded by scoring FHB symptoms on 20 heads/replication/line 21 days post-flowering using a visual scale described by Stack and McMullen[38].The DIS was calculated as(Incidence(INC)×Severity(SEV))/100.Percentage of Fusarium damaged kernel was evaluated using grain samples harvested by a low airspeed harvester.Sampled kernels were compared against a set of known FDK standards(https://agcrops.osu.edu/newsletter/corn-newsletter/2015-21/rating-fusarium-damaged-kernels-fdkscabby-wheat)to estimate the FDK values in two replications per sample.We used DIS data from three seasons,whereas FDK from 2019 and 2020 seasons for evaluating GP.For DON estimation,the samples were analysed at Department of Plant Science at North Dakota State University using gas chromatography-mass spectrometry method.As limited number of unreplicated samples were analysed for DON during 2018 and 2019,we used DON data only from the 2020 season.

    2.3.Genotyping-by-sequencing

    Fresh leaf tissues were taken from each line for DNA isolation using the hexadecyltrimethylammonium bromide(CTAB)method[39].The seedlings for each line in advanced or elite trials in their respective year of evaluation were grown in small pots/cones for tissue sampling and DNA extraction.A genotyping-by-sequencing(GBS)library was constructed using the PstI and MspI restriction enzymes[20]and sequenced in an IonTorrent Proton sequencer(Thermo Fisher Scientific,Waltham,MA,USA)at the USDA Central Small Grain Genotyping Lab,Manhattan,KS,USA.DNA sequence data was used to call single-nucleotide polymorphisms(SNP)using the previously described approach employing TASSEL v5.0(Trait Analysis by aSSociation,Evolution and Linkage)[40].After removal of SNPs with more than 30% missing genotypes,minor allele frequency(MAF)of less than 0.05 and unmapped on any chromosome,9321 high-quality SNPs were imputed using BEAGLE v4.1[41]for further analysis.

    2.4.Phenotypic data analysis

    The phenotypic data for FHB DIS and FDK was analyz ed as best linear unbiased estimates(BLUEs)for individual nurseries.The following model was used to estimate the BLUEs:

    where yijis the trait of interest,μis the overall mean,Riis the effect of the ithreplicate,Gjis the effect of the jthgenotype,and eijis the residual error effect associated with the ithreplication and jthgenotype.The lines shared between 2019 and 2020 seasons were used to combine the DIS and FDK data,BLUEs was estimated across environments using the following statistical model:

    where yijkis the trait of interest,μis the overall mean,Eiis the effect of the ithenvironment,Rj(i)is the effect of the jthreplicate nested within the ithenvironment,Gkis the effect of the kthgenotype,GEikis the effect of the genotype×environment(G×E)interaction,and eijkis the residual error associated with the replication and genotype effects.

    The broad-sense heritability(H2)for DIS and FDK was estimated for independent nurseries as follow:

    where σ2gandσ2e,are the genotype and error variance components,respectively.We used META-R[42]based on the LME4 R-package[43]for the linear mixed model analysis and heritability estimation.The Pearson correlations among traits and environments were estimated based on the BLUEs for each trait using R environment[44].

    The principal component analysis(PCA)was conducted using the genotypic data from 457 lines(257 lines from 2019 and 2020 nurseries and 200 lines from the breeding population)to study the relationship between the training and breeding populations.The‘prcomp’function in R was selected to perform PCA using 9321 SNP markers,and the first two principal components were used for the scatterplot.

    2.5.Genomic prediction models

    The univariate genomic prediction for DIS and FDK was performed using four different algorithms.The rrBLUP model[45]is the widely used GS model in plant breeding.In rrBLUP,we assume a normal distribution of marker effects with equal variance.The GEBVs for DIS and FDK were estimated for each training population using the trait BLUEs.A linear mixed model was implemented using the following model:

    where y is the vector(n×1)of adjusted means(BLUEs)from n genotypes for a given trait;μis the overall mean;Z is the design matrix(n×p)with known values of p markers for n genotypes;u is a genotypic predictor with u~N(0,),where G is positive semidefinite matrix,obtained from markers using‘A.mat’which is an additive relation matrix function andis the additive genetic variance;εis the residual error with e~N(0,).The model was implemented in the‘rrBLUP’R package[45]for one trait at a time.

    The rrBLUP model assumes common variance across the marker effects,which causes an underestimation of the large-effect QTL.However,the Bayesian method assumes unequal variances across marker effects and uses different priors to estimate these variances to overcome the limitations of rrBLUP[16].We used two Bayesian algorithms,BayesA(BA)and BayesB(BB),to estimate GEBVs for given traits.BA assumes that all markers have a non-zero effect by treating the proportion of markers with no effect(π)as zero.The markers are included in the model after shrinking their estimates to a normal distribution.BB is an extension of BA,which employs an inverse chi-square distribution for marker effects and assumes that some markers have no effect,which are excluded from the model.Thus,BB considers the presence of some large effect QTL controlling the given trait[16,46].A detailed description of Bayesian models can be found in Pérez and de los Campos[47].The Bayesian models were implemented in the‘BGLR’package using a Gibbs sampler with 5000 burn-in and 15,000 iterations for each run[47].Random Forest(RF)is an ensemble method that uses a collection of classification or regression trees to conduct prediction.The idea is that by combining a large number of smaller decision trees,RF can reduce the variance of prediction.The bias of RF models converges to a limiting value in the limits.RF methods have been applied for both genetic association studies and phenotype predictions[48].In this study,we extended the method used in Grinberg et al.[49]with 1000 iterations and 64 random states.The algorithm was implemented in Python using Sklearn library.A multivariate model was used to predict DIS and FDK by including DTH and PH as secondary traits in the model.Furthermore,we evaluated the performance of multivariate model to predict DON content using different combinations of secondary traits including DIS,FDK,DTH,and PH.A Bayesian Multivariate Gaussian model with an unstructured variance-covariance matrix was used for MT model[50].The MT model predicts FDK/DIS/DON using the secondary traits as described in the following:

    where y is the vector with a length of n×t(n genotypes and t traits);μis the means vector;Zrepresents the incidence matrix of order[(n×t)p],;u[(n×t)p]is genotypic predictor for all individuals and traits with u~N(0,∑?G).The matrix G represents the positive semidefinite matrix obtained from markers.The residuals of the MT model are represented by the vectorε,withε~N(0,R?I).The matrices∑and R are the variance-covariance matrices for depicting the genetic and residual effects for each individual in all traits,respectively.∑was estimated as an unstructured matrix and R as a diagonal matrix following Lado et al.[50].The variancecovariance matrices were estimated using a Gibbs sampler with 15,000 iterations,where the first 5000 iterations were used for burn-in.The MT model was implemented in R package‘MTM’[51].

    2.6.Training population design and cross-validation

    We evaluated the breeding lines from 2018,2019,and 2020 nurseries as training populations to predict the performance of preliminary breeding lines.Furthermore,we combined and evaluated the lines from 2019 and 2020 as one training population.Though the data for these two nurseries come from different years,a set of 58 lines that were common between the two nurseries was used to obtain the trait BLUEs.The resulting TPs were referred to as TP18,TP19,TP20,and TP19+20 for further genomic prediction.

    We first used each of these TPs for predicting single-year DIS and FDK.Briefly,each TP was randomly divided into five sets of equal size.Four of the five sets(80%)were used as a training set(phenotyped and genotyped)to train the model,and the remaining set(20%)was used as a testing set(genotyped only)for prediction(Fig.S1).Predictive ability was estimated as Pearson’s correlation between the GEBVs and observed phenotypes for the testing set.The predictions were assessed using five different models,as discussed earlier.The cross-validation process was repeated 500 times for rrBLUP and 100 times for Bayesian models,where each iteration included different lines in the training and testing sets.Cross-validation was used to evaluate the ability of four TPs using different models to estimate GEBVs for DIS and FDK.

    For the MT model,the lines were randomly split into a training set(80%)and a testing set(20%).To train the model,we used phenotypic data of secondary traits(PH and DTH)from both the training and testing sets,but the phenotypic data of the target trait(DIS or FDK)from the training set only(Fig.S1).

    As mentioned earlier,we used phenotypic data for DON from 2020 season to evaluate the performance of MT models and compare it with standard ST model(rrBLUP).The lines were randomly split into a training set(80%)and a testing set(20%)as in earlier case(Fig.S1).In this case,we evaluated the inclusion of different secondary traits(DIS,FDK,DTH,and PH)in different combinations to evaluate the predictive ability of MT model for DON.For instance,we used a MT model which included all four secondary traits to predict DON,then only three traits as secondary traits,and overall evaluating eight such combinations.

    2.7.Forward prediction of breeding lines for DIS and FDK

    Based on cross-validation analysis,we selected TP19,TP20,and TP19+20 for independent prediction of DIS and FDK in the preliminary lines from the breeding population.TP18 was not used for forward predictions as we did not observe a good PA for DIS and FDK in cross-validation.As all the prediction models yielded comparable results for DIS and FDK using cross-validation,we selected rrBLUP over other models for independent predictions based on its easy and less-intensive implementation.The model was trained using genotypic and phenotypic data from TP19,TP20,and TP19+20 in the‘rrBLUP’package to predict the GEBVs of 200 individuals in the breeding population.To assess the predictive ability,we randomly selected 60 lines of the BP,phenotyped these lines for DIS and FDK in the 2020 nursery,and used the observed phenotypic values to compare to the GEBVs from the TPs.

    3.Results

    3.1.Phenotypic analysis

    A significant genotypic variation in DIS and FDK(P less than 0.001)was observed in FHB nurseries from three seasons(Table S2).The largest variation in DIS was observed in 2019 ranging from 16.0 to 91.2 with a mean DIS of 49.1(Table 1),while the smallest variation was observed in 2018,ranging from 21.0 to 59.7(Table 1;Fig.1).The mean percent FDK were also different between 2019 and 2020,with 80% in 2019 and 58.6% in 2020.Despite high mean FDK values,sufficient phenotypic variation was observed among evaluated lines in both years(Table 1).Overall,2019 had the highest disease occurrence and 2018 the lowest.

    Significant phenotypic correlations(0.45 and 0.46)were observed between DIS and FDK in both 2019 and 2020.DIS and FDK exhibited negative correlations with plant height in all three nurseries(Figs.S1,S2 and S3).A negative correlation was observed between DIS and days to heading in 2018.However,both DIS and FDK showed weak but positive correlations with days to heading in 2019 and 2020(Fig.S2).As DON was also estimated in 2020,we observed significant positive phenotypic correlations for DON with DIS(0.34),FDK(0.33),and DTH(0.33)(Fig.S4).Broad-sense heritability for both FHB traits were moderate(Table 1)with the highest heritability for DIS(0.77)in 2020,and for FDK(0.75)in 2019.

    3.2.Relatedness between TP and BP

    The principal component analysis was conducted using 9321 SNP markers and 457 lines.The first two principal components(PCs)explained 6.9%(PC1)and 5.0%(PC2)of the genetic variance(Fig.1).The PCA revealed two primary clusters for the 457 lines,including 257 unique lines from the training populations(TP19 and TP20)and 200 lines of the breeding population(BP).The mixed distribution of the lines from both TPs and BP in the two clusters suggested a close relationship between the two populations,therefore it will be useful to use the TP for forward prediction in a breeding program.

    Table 1Descriptive statistics of two Fusarium head blight(FHB)resistance traits,disease index(DIS),and Fusarium-damaged kernels(FDK)for advanced winter wheat breeding lines evaluated in three independent FHB nurseries from 2018,2019,and 2020.

    Fig.1.The scatterplot for principal component analysis(PCA)of 457 lines based on 9321 SNP markers.The 457 breeding lines include 257 lines from training populations(TP19 and TP20)and 200 lines of the breeding population(BP)used in the forward prediction.The blue triangles represent the lines from TPs,and the green circles represents lines from BP.The first two PCs explained 6.9%and 5.0%of the total variation,respectively.

    3.3.Cross-validation within TPs for DIS and FDK

    Various GP models were used to evaluate the prediction accuracy within each TP using a cross-validation approach.We masked the observed phenotype in 20%of the lines in each TP,and treated them as untested new lines.The PA of DIS and FDK was measured as the Pearsons correlations between the predicted and the observed phenotypes of the masked lines.A general comparison of DIS and FDK prediction using different models is presented using boxplots in Figs.2 and 3.TP18 provided a low mean predictive ability of 0.20 for DIS,ranging from 0.15 to 0.23 using different GP models(Table 2).The TP20 provided a moderate mean predictive accuracy of 0.35 for DIS,and the TP19 yielded the highest mean prediction accuracy(0.39)for DIS,ranging from 0.37 to 0.41 using various models.Furthermore,for FDK,the mean predictive ability was 0.35 ranging from 0.32 to 0.37 in TP19,whereas,in TP20 mean predictive ability was slightly higher ranging from 0.34 to 0.40 with a mean of 0.37(Table 2).Among four different ST models,rrBLUP outperformed other models in all the TPs for predicting DIS and FDK(Figs.2 and 3).

    3.4.Combining two TPs into a large TP

    To improve prediction accuracies for different traits,we combined two individual populations(TP19 and TP20)into one large population(TP19+20).The BLUE values based on the shared set of lines evaluated in both years were used to assess PA for DIS and FDK using a cross-validation procedure.TP18 had poor performance in the cross-validation,therefore,it was excluded from the combined TP.The resulting population produced a slightly higher average PA of 0.41 for DIS and 0.38 for FDK using different GP models(Table 2;Figs.3,4)than for either TP19 or TP20.Besides a slight increase in the PA using larger TP,we observed that PA of all the models was quite similar,indicating a stable PA irrespective of the models.

    3.5.Univariate v/s multivariate prediction models for DIS and FDK

    PA also varied with types of GP models used for predicting DIS and FDK.Among four univariate models used for both DIS and FDK predictions,the rrBLUP model performed slightly better than BA,BB,and RF models in all TP scenarios(Table 2).When the univariate GP models were compared with a MT model with plant height and days to heading as covariates,the MT model did not show much improvement in DIS predictions than univariate models using TP18,TP19,and TP19+20.However,the prediction accuracy of MT was improved by 20%(0.35 to 0.42)for DIS in TP20 when two secondary traits(DTH and PH)were included in the model(Table 2).Improvement using the MT model in TP20 can be attributed to the moderate correlation between DIS and secondary traits in the model(Fig.S4).For FDK,the MT model did not improve prediction accuracy in any of the TPs(Table 2)as FDK seemed to be less correlated with DTH or PH in any of the growing seasons.

    3.6.Evaluating MT model to predict DON using different combinations of secondary traits

    The lines from 2020 nursery were also evaluated for DON along with other FHB traits.As DON is a costly trait to phenotype and a smaller number of samples are analyzed for DON each year,we were interested to see if we can use other FHB traits in MT models to predict DON.For this,we used the MT model with different combinations of DIS,FDK,DTH,and PH as covariates to predict DON.We used ST rrBLUP model as a benchmark to compare the performance of MT model(Fig.4).Using TP20,the rrBLUP model yielded a predictive accuracy of 0.49 for DON,which was higher as compared to DIS and FDK using the same model.The MT model showed an improvement of up to 20%,yielding PA ranging from 0.54 to 0.59 with different combinations of secondary traits(Table S3).The MT model having DIS,FDK,DTH,and PH as covariates had PA of 0.56,whereas the MT model with DIS,FDK,and DTH had the highest PA(0.59)among all the combinations(Table S3).We also evaluated the MT model if only single trait(DIS,FDK,or DTH)is used as covariate to predict DON.Interestingly,the MT model with only FDK as covariate had PA of 0.56,which is comparable to the combination having all traits as covariates(Table S3).

    Fig.2.The predictive ability(PA)for Fusarium head blight(FHB)resistance disease index(DIS)in different sets of training populations(TPs)used in the study.Boxplots compare the PA using five genomic prediction models:rrBLUP,ridge-regression best linear unbiased prediction;BA,BayesA;BB,BayesB;RF,Random Forest;and MT,multitrait model.TP18,training population based on 2018 FHB nursery;TP19,training population based on 2019 FHB nursery;TP20,training population based on 2020 FHB nursery;TP19+20,training population combining 2019 and 2020 FHB nurseries.

    Fig.3.The predictive ability(PA)for Fusarium-damaged kernels(FDK)in different sets of training populations(TPs)used in the study.Boxplots compare the PA using five genomic prediction models:rrBLUP,ridge-regression best linear unbiased prediction;BA,BayesA;BB,BayesB;RF,Random Forest;and MT,multi-trait model.TP19,training population based on 2019 FHB nursery;TP20,training population based on 2020 FHB nursery;TP19+20,training population combining 2019 and 2020 FHB nurseries.

    Table 2Mean prediction accuracy and standard error for two Fusarium head blight(FHB)resistance traits,disease index(DIS)and Fusarium-damaged kernels(FDK)with cross-validation in different Training Populations using different genomic prediction models.

    Fig.4.Boxplots comparing the predictive ability of ST model and the MT model with different combinations of secondary traits.DIS,FHB disease Index;FDK,fusarium damaged kernel percentage;DTH,days to heading;PH,plant height;DON,deoxynivalenol content.

    3.7.Accuracy of independent predictions for DIS and FDK in forward breeding

    To validate usefulness of different FHB nurseries as possible TPs,GEBVs for DIS and FDK of a random set of preliminary breeding lines from PYT were predicted using TP19,TP20,and TP19+20 with the rrBLUP model.Moderate prediction accuracies(Table 3)were observed for DIS and FDK using the three TPs.The TP19 provided the highest prediction accuracy(0.59)for DIS,following the trend observed in our cross-validation.The TP19+20 produced the highest prediction accuracy of 0.54 for FDK,followed by TP19 and TP20(0.50 and 0.49,Table 3).Overall,independent predictions provided better PA than the cross-validation.Furthermore,we used a scatterplot to compare the breeding lines which were rejected based on estimated GEBVs but retained based on observed data(Figs.5,S5).Interestingly,we observed there was low probability of rejecting lines with lower DIS or FDK ratings as most of lines rejected based on GEBVs were having a moderate DIS or FDK observed value in both 2019 and 2020(refer to top-left quadrant of scatterplots in Figs.5,S5).These results demonstrate that the genomic prediction can be implemented to improve FHB resistance in wheat breeding programs.

    4.Discussion

    Several studies have evaluated the inclusion of GS in wheat breeding programs to predict FHB resistance in recent years[25,26,29,31];however,most of these studies were done in soft winter wheat.The current study used a cross-validation strategy to evaluate the potential of hard winter wheat breeding lines as training populations for GS of FHB resistance.The multi-trait GP model was also evaluated for predicting different FHB traits.Finally,the current study demonstrates the use of different sets of advanced breeding lines as training sets to predict preliminary breeding lines in a forward prediction scheme.

    Table 3Prediction accuracy for two Fusarium head blight(FHB)resistance traits,disease index(DIS)and Fusarium-damaged kernels(FDK)in forward prediction scheme.

    4.1.Phenotypic response to FHB

    We used an FHB disease index estimated based on incidence and severity,and FDK percentage to evaluate FHB resistance in advanced breeding lines.We observed a wide variation for FHB resistance in the advanced breeding lines from the South Dakota State University winter wheat breeding program in three disease nurseries from 2018 to 2020.For example,variations for DIS were from 16.0%to 91.2%and FDK from 38.1%to 99.0%in the 2019 nursery(Table 1).A similar trend was observed in two other nurseries(Table 1).The majority of the advanced and elite breeding material from our program does not carry Fhb1 likely due to yield drag and negative agronomic potential,suggesting minor genes from native sources predominantly govern FHB resistance in our breeding program.This is similar to other hard winter wheat breeding programs in the region as only one variety(TAM 205)carrying Fhb1 has been released till date.Therefore,GS seems to be a suitable approach to breed for the FHB resistance.

    Although the FHB nursery had a controlled mist system,the prevailing environment in respective years is believed to play a significant role in varying disease pressure.The 2018 season was very dry in South Dakota,but it was reasonably wet in 2019.These environmental factors lead to fluctuations in the disease pressure affecting the spread of the data over the years,consistent with previous studies[28,30,31,36].However,these fluctuations did not alter the ranking of different check genotypes over different nurseries.For instance,the DIS indices were 27.2,29.8,and 32.6 in 2018,2019,and 2020,respectively,for resistant check‘Lyman’,and 84.6 and 72.4 in 2019 and 2020,respectively,for susceptible check‘Overley’.The consistent performance of checks shows uniformity and reliability of phenotyping across nurseries,which is supported by higher heritability for DIS(0.54-0.77)and FDK(0.66-0.75)(Table 1).The moderate H2estimates were in similar range as that of related studies using different types of population in spring-or winter wheat[30,36,52-54].

    Fig.5.Scatterplots showing observed v/s predicted values for FHB disease index(DIS).The observed DIS estimates were based on phenotypic evaluation of 60 independent lines,and the predicted values are GEBVs using(A)TP19,(B)TP20,and(C)TP19+20.The red dashed line represents the cutoff(50%)to discard genotypes based on the observed data and estimated GEBVs.Genotype to the right side of the red dashed line would be discarded based on observed phenotype,and genotype above the red dashed line would be discarded based on GEBVs.

    4.2.Within-population cross-validation for DIS and FDK

    Several studies successfully used advanced or preliminary breeding lines from a breeding program[28,36,37]and from unrelated regional nurseries[31]as training populations to predict various traits in wheat.We evaluated advanced breeding lines in three FHB nurseries as three possible training sets for genomic prediction using cross-validation(Table S1),and obtained moderate prediction accuracies when TP19 and TP20 were used to predict DIS(up to 0.41 using TP19 and 0.42 in TP20)and FDK(up to 0.37 in TP19 and 0.40 in TP20)(Table 2).The prediction accuracy for both traits was consistent with several previous studies in soft winter or spring wheat[29,31,36,54].Similarly,prediction accuracies for FDK in our study were also similar to Rutkoski et al.[25]and Adeyemo et al.[36],but lower than those reported by Arruda et al.[26]and Verges et al.[31].The poor performance of TP18 for DIS prediction could result from low disease pressure and phenotypic variation observed in the 2018 nursery,hence,such a nursery is not recommended for a forward genomic prediction.Overall,based on mutliple years of data our study demonstrates the usefulness of advanced breeding lines as TPs for genomic prediction,given that the quality of phenotyping is robust.

    Training population size is another crucial factor that affects the PA of GP models.Previous studies have reported an increase in PA with an increased TP size[55-57].Lorenz et al.[57]and Arruda et al.[26]obtained higher PA for FHB traits in barley and wheat when the TP size contained 250 to 300 lines.In the current study,the three populations,each with around 150 lines,were used as independent TPs for FHB trait prediction.Also,the performance of a larger TP(TP19+20)developed by combining two TPs(TP19 and TP20)were evaluated.Although TP19 and TP20 were phenotyped in two independent nurseries,they shared a subset of 58 breeding lines.Thus,a large TP(TP19+20)was formed with 265 unique lines and a higher PA was observed for DIS(0.41)and FDK(0.38)when TP19+20 was used as the TP.Apart from improvement in PA using TP19+20,we observed that all prediction models yielded similar results which was not the case when using smaller TPs(Table 2).Furthermore,a lower standard deviation for PA based on several repeats of cross-validation suggests consistency of prediction when using a large TP(Figs.3 and 4).The results of this study suggest that selection of the right TP is crucial to improve PA.TP19+20 showed the best PA due to the large TP size and correction of BLUEs for DIS and FDK across two different nurseries/years,thus it can be used as the TP in forward genomic prediction of preliminary breeding lines.

    4.3.Model comparison for DIS and FDK

    We compared the PA for DIS and FHB using four univariate and one multivariate GP model.Among the univariate models,rrBLUP outstripped the other three models for predicting DIS and FDK in all individual TPs(Table 2).Further,rrBLUP is preferred model for predicting FHB traits owing to its better performance and computational advantage.Previous studies have also reported that rrBLUP has better PA than Bayesian models[25,26]is one of the most often used methods in GS for FHB resistance.

    Using a MT GP model is another approach to increase the PA for FHB traits.This model includes correlated secondary traits such as PH and DTH as covariates to predict DIS and FDK.Schulthess et al.[32]and Larkin et al.[30]reported an increase in PA by including PH or DTH in the GP model.In the current study,PH and DTH were included as secondary traits in a MT-GBLUP model to predict DIS and FDK,however,improvement of PA were not considerable except for DIS in TP20 where the PA was improved by at least 20%over the other univariate models(Table 2).The better PA using MT model in the TP20 resulted from the moderate correlation of DIS with PH and DTH(Fig.S1).In other TPs,a lower correlation was observed among evaluated traits.The MT models have been used to improve the PA of low-heritability traits by using the information from correlated traits with high heritability[23,33,58].However,moderate to high heritability estimates for DIS and FDK in this study could be another reason for MT models showing no advantages over univariate models.Thus,our results suggest that the MT models could be useful and employed in forward prediction in early generation nurseries if the observed correlations between FHB traits and the covariates are moderate to high.

    4.4.Optimization of MT model with different combinations of secondary traits to predict DON

    DON is an important and primary FHB trait;however,most of the winter wheat breeders across the US are unable to make decisions based on DON estimates.For instance,the South Dakota State University winter wheat program harvests FHB nursery in August and then plants the next cycle in September,making the turnaround cycle very short.Contrary to this,phenotyping for DON is mostly outsourced and it takes several months before the breeder gets the data.Thus,it will be of great importance if we can predict DON and utilize the predictions to advance resistant lines to the next generation.Furthermore,previous studies have suggested the use of different pre-harvest traits,such as DIS,FDK,DTH or PH,in MT models to predict DON with better accuracy[30].We used DON estimates available for TP20 to evaluate the predictive ability for DON using ST and different versions of MT model based on cross-validation.The ST(rrBLUP)model predicted DON with PA of 0.49,which was better than DIS and FDK using any of the models(Table S3).Further,we compared MT model having different combinations of traits as covariates and observed that trait combinations of DIS,FDK,DTH(0.59)and FDK,DTH(0.58)had the highest PA(Fig.4),similar to the results reported by Larkin et al.[30].Interestingly,it was found that using FDK as‘only’covariate in the MT model had high PA(0.56),suggesting that using only FDK can improve the performance of MT model over the ST model.Overall,the results suggest that we can use GP for predicting DON in earlier stages and the PA can be further improved by using secondary traits such as FDK.

    4.5.Genomic predictions in forward breeding

    Relatedness among the individuals in TP and BP is considered crucial for getting higher PA in genomic prediction.In this study,we aimed to evaluate the usefulness of advanced breeding lines as TPs to predict earlier generations lines from our breeding program.The PCA showed a good association between the lines from TP and BP(Fig.2).Hence,we obtained moderate prediction accuracy for DIS and FDK in an independent BP(Table 3),comparable to the other reports[31,59,60].Though a moderate PA for DIS and FDK were achieved in the current study,the phenotypic data of advanced lines evaluated from FHB nurseries could be used to predict GEBVs for the early generation breeding lines in wheat breeding programs(Table 3).The predicted GEBVs can be used to discard susceptible lines at the earlier stages of the breeding cycle(Figs.3,5).For instance,we achieved a PA of 0.59 for DIS using TP19 and discarded 50% of the most susceptible lines based on the TP19 based GEBVs and selected the remaining 50% lines for further selection(Table 3).Among these discarded lines,87%were highly susceptible based on their observed disease index in the mist irrigated inoculated FHB nursery,and only 13%susceptible genotypes were carried forward based on GEBVs(Fig.5).Similar results were observed for FDK(Fig.3),which suggests that the estimated GEBVs from the TP can be used to discard the most susceptible lines at an early stage to reduce phenotyping costs.Our results also suggest that these lines discarded based on GEBVs had an extremely low chance to be highly resistant in terms of DIS and FDK.

    In summary,our study demonstrated that advanced breeding lines evaluated in the FHB nurseries can serve as TPs and predict GEBVs for untested earlier generation lines in the breeding programs.Further advanced lines from several years can be combined to increase the PA in the forward breeding.However,we recommend evaluating the performance of individual advanced breeding nurseries through cross-validation before pooling multiple nurseries to develop larger TPs for forward prediction.Furthermore,the MT models using secondary traits as covariates could be useful in predicting cumbersome FHB traits(DON and FDK).Finally,our results suggest genomic prediction can be successfully applied in a wheat breeding program in discarding the most FHB susceptible lines at an early stage as opposed to laborious phenotypic selection in most years and especially in abnormal years when phenotypic evaluation is unreliable or unavailable due to environmental conditions.

    CRediT authorship contribution statement

    Jinfeng Zhang:Conceptualization;Investigation;Formal analysis;Visualization;Software implementation;Writing-original draft.Harsimardeep S.Gill:Conceptualization;Investigation;Formal analysis;Data curation;Visualization;Software implementation;Writing-original draft.Navreet K Brar:Investigation.Jyotirmoy Halder:Formal analysis;Investigation.Shaukat Ali:Investigation;Writing-review & editing.Xiaotian Liu:Formal analysis.Amy Bernardo:Genotyping and SNP discovery analysis.Paul St Amand:Genotyping and SNP discovery analysis.Guihua Bai:Genotyping and SNP discovery analysis;Interpretation of results;Writing-review&editing.Upinder S.Gill:Interpretation of results;Writing-review&editing.Brent Turnipseed:Interpretation of results;Writing-review & editing.Sunish K.Sehgal:Conceptualization;Data curation;Visualization;Methodology;Supervision;Resources;Writing-original draft,review & editing.

    Declaration of competing interest

    The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

    Acknowledgments

    This project was collectively funded by the USDA hatch projects SD00H695-20,USDA-ARS agreement 59-0206-0-177(USDAUSWBSI),the Agriculture and Food Research Initiative Competitive Grant 2022-68013-36439(WheatCAP)from the USDA National Institute of Food and Agriculture,and South Dakota Wheat Commission Grant 3X1340.The funders had no role in the study design,data collection and analysis,decision to publish,or manuscript preparation.The authors would like to thank the South Dakota Agriculture Experimental Station(Brookings,SD,USA)for providing the resources to conduct the experiments.

    Appendix A.Supplementary data

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

    亚洲片人在线观看| 18禁裸乳无遮挡动漫免费视频| 欧美精品一区二区免费开放| 国产熟女午夜一区二区三区| 亚洲一区高清亚洲精品| 国产精品久久电影中文字幕 | 国产av精品麻豆| 国产有黄有色有爽视频| 嫁个100分男人电影在线观看| 操出白浆在线播放| 自拍欧美九色日韩亚洲蝌蚪91| 成人国语在线视频| 亚洲情色 制服丝袜| 亚洲色图综合在线观看| 无限看片的www在线观看| 一级a爱视频在线免费观看| 国产精品99久久99久久久不卡| 成人av一区二区三区在线看| av在线播放免费不卡| 人妻一区二区av| 亚洲精品粉嫩美女一区| 99riav亚洲国产免费| 中文欧美无线码| 久久这里只有精品19| 欧美一级毛片孕妇| 精品国产乱子伦一区二区三区| 啦啦啦视频在线资源免费观看| 好看av亚洲va欧美ⅴa在| 亚洲色图av天堂| 在线十欧美十亚洲十日本专区| 亚洲中文字幕日韩| 涩涩av久久男人的天堂| 亚洲成人免费电影在线观看| 一级毛片高清免费大全| 一级作爱视频免费观看| 一边摸一边抽搐一进一小说 | 国产午夜精品久久久久久| 亚洲欧美色中文字幕在线| 韩国精品一区二区三区| 在线永久观看黄色视频| 男女之事视频高清在线观看| 亚洲一卡2卡3卡4卡5卡精品中文| 高清av免费在线| 免费在线观看日本一区| 高清欧美精品videossex| 建设人人有责人人尽责人人享有的| 午夜视频精品福利| 国产主播在线观看一区二区| 99精品欧美一区二区三区四区| 午夜福利免费观看在线| 欧美亚洲 丝袜 人妻 在线| 欧美日韩视频精品一区| 国产熟女午夜一区二区三区| svipshipincom国产片| 高潮久久久久久久久久久不卡| 丝袜美足系列| 亚洲一区中文字幕在线| 嫩草影视91久久| 国产又色又爽无遮挡免费看| 欧美激情高清一区二区三区| 18禁裸乳无遮挡免费网站照片 | 免费女性裸体啪啪无遮挡网站| 午夜亚洲福利在线播放| 亚洲一卡2卡3卡4卡5卡精品中文| 十八禁人妻一区二区| 在线观看日韩欧美| 精品久久久久久久久久免费视频 | 国产精品 国内视频| 久久香蕉激情| 亚洲欧美激情综合另类| 一区在线观看完整版| 国产亚洲欧美98| 久久ye,这里只有精品| 香蕉久久夜色| 身体一侧抽搐| 日韩欧美一区视频在线观看| 国产精品.久久久| cao死你这个sao货| 国产精品.久久久| 亚洲av日韩精品久久久久久密| 日本欧美视频一区| 身体一侧抽搐| 中文亚洲av片在线观看爽 | 久久久久久久精品吃奶| 99香蕉大伊视频| 久久久久久亚洲精品国产蜜桃av| 国产成人啪精品午夜网站| 欧美另类亚洲清纯唯美| bbb黄色大片| 亚洲五月婷婷丁香| cao死你这个sao货| 免费看十八禁软件| 女人精品久久久久毛片| 精品一区二区三区av网在线观看| xxx96com| 9热在线视频观看99| 热re99久久精品国产66热6| 亚洲美女黄片视频| 国产高清激情床上av| 十八禁网站免费在线| 麻豆国产av国片精品| 亚洲精品久久午夜乱码| 久热爱精品视频在线9| 老熟女久久久| 国产一区二区激情短视频| 黄频高清免费视频| 亚洲综合色网址| 亚洲av片天天在线观看| 大型黄色视频在线免费观看| 成年人免费黄色播放视频| bbb黄色大片| 亚洲在线自拍视频| 国产人伦9x9x在线观看| 国内久久婷婷六月综合欲色啪| 在线av久久热| 亚洲人成电影免费在线| 免费观看人在逋| 91av网站免费观看| 亚洲黑人精品在线| 99久久精品国产亚洲精品| 久久青草综合色| 国产精品乱码一区二三区的特点 | 99精品久久久久人妻精品| 亚洲国产中文字幕在线视频| 欧美不卡视频在线免费观看 | 免费看a级黄色片| 国产精品九九99| 国产精品成人在线| 国产精品久久电影中文字幕 | 侵犯人妻中文字幕一二三四区| av线在线观看网站| 99热只有精品国产| 天天操日日干夜夜撸| 亚洲人成电影免费在线| 一级毛片女人18水好多| 三上悠亚av全集在线观看| 80岁老熟妇乱子伦牲交| 久久天躁狠狠躁夜夜2o2o| 操美女的视频在线观看| 亚洲国产欧美网| 亚洲av片天天在线观看| 亚洲成人免费av在线播放| 一二三四社区在线视频社区8| 国产成人精品无人区| 亚洲成国产人片在线观看| 人人妻人人添人人爽欧美一区卜| 一二三四在线观看免费中文在| 国产亚洲精品第一综合不卡| 亚洲中文av在线| 国产精品成人在线| 伦理电影免费视频| 国产亚洲精品久久久久久毛片 | 日韩人妻精品一区2区三区| 操出白浆在线播放| 久久久国产一区二区| 王馨瑶露胸无遮挡在线观看| 啦啦啦视频在线资源免费观看| 亚洲自偷自拍图片 自拍| 好男人电影高清在线观看| 免费在线观看完整版高清| 麻豆av在线久日| 99精品久久久久人妻精品| 精品国内亚洲2022精品成人 | 99精品欧美一区二区三区四区| 男女免费视频国产| 国产精品一区二区在线不卡| 国产亚洲一区二区精品| 少妇被粗大的猛进出69影院| 国产亚洲精品久久久久5区| 首页视频小说图片口味搜索| 免费在线观看视频国产中文字幕亚洲| 国产精品综合久久久久久久免费 | 少妇裸体淫交视频免费看高清 | 国产精品免费大片| 超碰成人久久| 午夜福利欧美成人| 十八禁人妻一区二区| 搡老乐熟女国产| 亚洲精品一卡2卡三卡4卡5卡| av在线播放免费不卡| 亚洲av熟女| 久久草成人影院| 搡老熟女国产l中国老女人| 大型黄色视频在线免费观看| 国产有黄有色有爽视频| 国产男女内射视频| 国产精品免费视频内射| 99re在线观看精品视频| 亚洲aⅴ乱码一区二区在线播放 | 亚洲国产精品合色在线| 99re在线观看精品视频| 最近最新免费中文字幕在线| 亚洲国产中文字幕在线视频| 欧美日韩av久久| 啦啦啦免费观看视频1| 久久人妻福利社区极品人妻图片| 777米奇影视久久| 女性被躁到高潮视频| av福利片在线| 亚洲 欧美一区二区三区| 国产精品99久久99久久久不卡| 久热这里只有精品99| 在线看a的网站| 女人被狂操c到高潮| 99久久99久久久精品蜜桃| 午夜精品国产一区二区电影| www.999成人在线观看| 两性夫妻黄色片| 久久午夜综合久久蜜桃| 国产成人欧美| 日本黄色日本黄色录像| 成年人黄色毛片网站| 欧美久久黑人一区二区| 亚洲成人免费电影在线观看| 美女高潮到喷水免费观看| 在线观看日韩欧美| 一本大道久久a久久精品| 亚洲成人免费av在线播放| 俄罗斯特黄特色一大片| 一区在线观看完整版| 脱女人内裤的视频| 自线自在国产av| 免费少妇av软件| 亚洲aⅴ乱码一区二区在线播放 | 欧美日韩福利视频一区二区| 国产激情欧美一区二区| 大型黄色视频在线免费观看| 亚洲中文av在线| 老司机午夜十八禁免费视频| 韩国av一区二区三区四区| 男女午夜视频在线观看| 免费女性裸体啪啪无遮挡网站| 久久午夜亚洲精品久久| 人妻久久中文字幕网| 一边摸一边抽搐一进一出视频| 国产精品香港三级国产av潘金莲| 老鸭窝网址在线观看| 免费久久久久久久精品成人欧美视频| 老司机在亚洲福利影院| 国产av一区二区精品久久| 久久精品91无色码中文字幕| 天天操日日干夜夜撸| 十八禁网站免费在线| 在线观看午夜福利视频| 99久久国产精品久久久| 日日摸夜夜添夜夜添小说| 国产精品1区2区在线观看. | 91成人精品电影| 超色免费av| 国产在视频线精品| www.999成人在线观看| 国产成人欧美在线观看 | 国产伦人伦偷精品视频| 窝窝影院91人妻| 美女视频免费永久观看网站| 男女之事视频高清在线观看| 韩国av一区二区三区四区| 亚洲成av片中文字幕在线观看| 成熟少妇高潮喷水视频| 动漫黄色视频在线观看| 又大又爽又粗| 欧美国产精品va在线观看不卡| 国产又色又爽无遮挡免费看| 亚洲一区高清亚洲精品| 精品人妻1区二区| 久久性视频一级片| 免费在线观看视频国产中文字幕亚洲| 一区福利在线观看| 国产日韩欧美亚洲二区| 一级,二级,三级黄色视频| 日韩免费高清中文字幕av| 亚洲成av片中文字幕在线观看| 国产精品免费大片| 亚洲国产看品久久| 亚洲精品乱久久久久久| 十八禁高潮呻吟视频| 麻豆国产av国片精品| 国产成+人综合+亚洲专区| 亚洲男人天堂网一区| 成人国语在线视频| 国产亚洲精品第一综合不卡| 日本vs欧美在线观看视频| 国产99白浆流出| 亚洲黑人精品在线| www日本在线高清视频| 欧美日韩中文字幕国产精品一区二区三区 | 欧美最黄视频在线播放免费 | 99re在线观看精品视频| 国产97色在线日韩免费| 伊人久久大香线蕉亚洲五| 亚洲人成电影观看| 老司机在亚洲福利影院| 精品福利永久在线观看| 亚洲午夜精品一区,二区,三区| 极品教师在线免费播放| 午夜91福利影院| 欧美黑人欧美精品刺激| 91老司机精品| 黄片小视频在线播放| 欧洲精品卡2卡3卡4卡5卡区| 国产xxxxx性猛交| 亚洲成a人片在线一区二区| 老熟妇仑乱视频hdxx| 国产成人av教育| 亚洲欧美日韩高清在线视频| 精品亚洲成a人片在线观看| 久久人妻熟女aⅴ| 久久久久视频综合| 日本五十路高清| 久久久久久久精品吃奶| 中文字幕高清在线视频| 丰满的人妻完整版| 少妇被粗大的猛进出69影院| 日韩一卡2卡3卡4卡2021年| 老汉色∧v一级毛片| 天天添夜夜摸| 欧美日韩视频精品一区| 国产精品亚洲一级av第二区| 高清视频免费观看一区二区| 精品一区二区三区四区五区乱码| 免费日韩欧美在线观看| 99热只有精品国产| 很黄的视频免费| 曰老女人黄片| 欧美黄色片欧美黄色片| 校园春色视频在线观看| 国产成+人综合+亚洲专区| 午夜福利在线免费观看网站| 国产aⅴ精品一区二区三区波| 久久精品熟女亚洲av麻豆精品| 黄片大片在线免费观看| 午夜福利影视在线免费观看| 乱人伦中国视频| 亚洲五月天丁香| 啦啦啦视频在线资源免费观看| 日韩欧美三级三区| 黄色丝袜av网址大全| 国产深夜福利视频在线观看| 99精品久久久久人妻精品| 亚洲自偷自拍图片 自拍| av视频免费观看在线观看| 欧美成人免费av一区二区三区 | 亚洲成人手机| 久热这里只有精品99| 后天国语完整版免费观看| 激情视频va一区二区三区| 一个人免费在线观看的高清视频| 亚洲欧美激情综合另类| 日韩中文字幕欧美一区二区| 精品人妻熟女毛片av久久网站| 欧美黄色淫秽网站| 一级a爱视频在线免费观看| 国产一区二区三区综合在线观看| 成人黄色视频免费在线看| 国内久久婷婷六月综合欲色啪| 三上悠亚av全集在线观看| 狠狠狠狠99中文字幕| 国产高清videossex| 中文字幕制服av| 国产亚洲精品久久久久5区| 国产精品久久电影中文字幕 | 精品少妇一区二区三区视频日本电影| 成人亚洲精品一区在线观看| 欧美大码av| 99国产极品粉嫩在线观看| 在线观看午夜福利视频| 91在线观看av| 久久天堂一区二区三区四区| 女人精品久久久久毛片| 久久国产精品影院| 欧美激情极品国产一区二区三区| 日本wwww免费看| 欧美丝袜亚洲另类 | 超碰成人久久| 一本一本久久a久久精品综合妖精| 亚洲成a人片在线一区二区| 欧美国产精品va在线观看不卡| 波多野结衣av一区二区av| 高清视频免费观看一区二区| 岛国毛片在线播放| 狠狠狠狠99中文字幕| 18禁美女被吸乳视频| 99riav亚洲国产免费| 成熟少妇高潮喷水视频| 无限看片的www在线观看| 欧美av亚洲av综合av国产av| 国产不卡av网站在线观看| 麻豆av在线久日| 在线播放国产精品三级| 亚洲专区字幕在线| 女人高潮潮喷娇喘18禁视频| 91成年电影在线观看| 久久性视频一级片| 欧美中文综合在线视频| 女同久久另类99精品国产91| 午夜福利免费观看在线| 国产人伦9x9x在线观看| 久久国产精品男人的天堂亚洲| 后天国语完整版免费观看| 欧美不卡视频在线免费观看 | 50天的宝宝边吃奶边哭怎么回事| 国产精品香港三级国产av潘金莲| 飞空精品影院首页| 精品国产乱码久久久久久男人| 国产精品欧美亚洲77777| 久久久久久久久免费视频了| 国产一区有黄有色的免费视频| 国产精华一区二区三区| 欧美日韩国产mv在线观看视频| 亚洲视频免费观看视频| 国产精品亚洲av一区麻豆| 婷婷成人精品国产| 男女床上黄色一级片免费看| 老司机午夜十八禁免费视频| 国产精品美女特级片免费视频播放器 | 日日爽夜夜爽网站| 亚洲av美国av| 亚洲色图 男人天堂 中文字幕| 亚洲一区中文字幕在线| 亚洲av日韩精品久久久久久密| 亚洲欧洲精品一区二区精品久久久| 欧美+亚洲+日韩+国产| 日韩大码丰满熟妇| 一本综合久久免费| 久久久水蜜桃国产精品网| 午夜免费鲁丝| 高清视频免费观看一区二区| 看片在线看免费视频| 黄片小视频在线播放| 熟女少妇亚洲综合色aaa.| 新久久久久国产一级毛片| 亚洲免费av在线视频| 老熟妇仑乱视频hdxx| 丝袜人妻中文字幕| 国产成人av激情在线播放| 欧美日韩中文字幕国产精品一区二区三区 | 久久热在线av| 亚洲片人在线观看| 两个人看的免费小视频| av欧美777| 99久久精品国产亚洲精品| 久热爱精品视频在线9| 午夜精品在线福利| 一区二区三区激情视频| 精品福利永久在线观看| 一区二区三区激情视频| 母亲3免费完整高清在线观看| 久久精品国产综合久久久| 一个人免费在线观看的高清视频| 正在播放国产对白刺激| 19禁男女啪啪无遮挡网站| 午夜日韩欧美国产| 婷婷成人精品国产| 国产成人免费观看mmmm| 精品一区二区三卡| 精品一区二区三区av网在线观看| 日韩视频一区二区在线观看| 亚洲色图av天堂| 成人手机av| 国产深夜福利视频在线观看| 亚洲精品自拍成人| 99久久99久久久精品蜜桃| 国产午夜精品久久久久久| 久久人人爽av亚洲精品天堂| 一区二区日韩欧美中文字幕| 欧美日韩亚洲国产一区二区在线观看 | 久久久久国内视频| 老司机亚洲免费影院| 1024香蕉在线观看| √禁漫天堂资源中文www| 国产真人三级小视频在线观看| 亚洲人成电影免费在线| 12—13女人毛片做爰片一| 他把我摸到了高潮在线观看| 多毛熟女@视频| 又黄又爽又免费观看的视频| 少妇裸体淫交视频免费看高清 | 又大又爽又粗| 黄色成人免费大全| 亚洲中文av在线| 久久久久久久久免费视频了| 日韩欧美国产一区二区入口| 国产一区二区三区在线臀色熟女 | 狠狠狠狠99中文字幕| 91大片在线观看| 欧美在线黄色| 国产97色在线日韩免费| 女人高潮潮喷娇喘18禁视频| 热99re8久久精品国产| 国产亚洲av高清不卡| 亚洲一码二码三码区别大吗| 国产精品乱码一区二三区的特点 | 91大片在线观看| 精品国产亚洲在线| 亚洲五月婷婷丁香| 黄色成人免费大全| 99国产精品一区二区蜜桃av | 久久影院123| 韩国精品一区二区三区| 男人的好看免费观看在线视频 | 夜夜爽天天搞| 中文字幕最新亚洲高清| 悠悠久久av| 国产精品免费一区二区三区在线 | 欧美日韩国产mv在线观看视频| 热re99久久精品国产66热6| 久久久精品免费免费高清| 欧美色视频一区免费| 又黄又爽又免费观看的视频| 国内毛片毛片毛片毛片毛片| 91九色精品人成在线观看| 久久久久精品人妻al黑| 后天国语完整版免费观看| 午夜两性在线视频| 在线观看午夜福利视频| 亚洲专区国产一区二区| 视频在线观看一区二区三区| svipshipincom国产片| 美女午夜性视频免费| 国产麻豆69| 超碰97精品在线观看| 人人妻人人澡人人爽人人夜夜| 黄色a级毛片大全视频| 国产精品香港三级国产av潘金莲| 欧美一级毛片孕妇| 久久久国产欧美日韩av| 国产成人免费观看mmmm| www.自偷自拍.com| 午夜精品久久久久久毛片777| 老熟女久久久| 首页视频小说图片口味搜索| 国产在线观看jvid| 大片电影免费在线观看免费| 999久久久国产精品视频| 在线观看免费日韩欧美大片| 一区二区三区激情视频| 一边摸一边抽搐一进一小说 | 在线十欧美十亚洲十日本专区| 欧美激情高清一区二区三区| 99国产精品一区二区三区| 精品国产一区二区三区久久久樱花| 午夜免费鲁丝| 大香蕉久久网| 免费观看精品视频网站| 91精品三级在线观看| 欧美日韩瑟瑟在线播放| 99精国产麻豆久久婷婷| 亚洲一区二区三区不卡视频| 久久国产精品大桥未久av| 麻豆乱淫一区二区| 久久天躁狠狠躁夜夜2o2o| 久久香蕉精品热| 免费人成视频x8x8入口观看| 欧美+亚洲+日韩+国产| 亚洲欧美激情综合另类| 国产高清视频在线播放一区| 亚洲综合色网址| 另类亚洲欧美激情| 搡老乐熟女国产| 国产不卡av网站在线观看| 一本综合久久免费| 最近最新中文字幕大全免费视频| 午夜日韩欧美国产| 日韩三级视频一区二区三区| 窝窝影院91人妻| 亚洲av电影在线进入| 亚洲av片天天在线观看| 欧美日本中文国产一区发布| 国产成人欧美在线观看 | 久99久视频精品免费| 欧美日韩亚洲高清精品| 一进一出抽搐动态| 嫩草影视91久久| 精品无人区乱码1区二区| 精品国产一区二区三区四区第35| svipshipincom国产片| 另类亚洲欧美激情| 丰满人妻熟妇乱又伦精品不卡| 激情在线观看视频在线高清 | 免费一级毛片在线播放高清视频 | 亚洲性夜色夜夜综合| av电影中文网址| 日韩人妻精品一区2区三区| 激情在线观看视频在线高清 | 热re99久久精品国产66热6| 18禁裸乳无遮挡免费网站照片 | 波多野结衣av一区二区av| 亚洲五月天丁香| 9热在线视频观看99| 久久人妻熟女aⅴ| 欧美 亚洲 国产 日韩一| 久久99一区二区三区| 一边摸一边抽搐一进一出视频| 亚洲黑人精品在线| 久热爱精品视频在线9| 操美女的视频在线观看| 午夜免费成人在线视频| 女人被狂操c到高潮| 日本黄色视频三级网站网址 | 黄色片一级片一级黄色片| 精品亚洲成国产av| 欧美丝袜亚洲另类 | 亚洲欧美一区二区三区黑人| 91精品国产国语对白视频| 一区二区日韩欧美中文字幕| 人人妻人人澡人人看| 久久香蕉精品热| 大型av网站在线播放| 日韩免费高清中文字幕av| av视频免费观看在线观看| 国产成人欧美在线观看 | 中文字幕制服av| 天天操日日干夜夜撸| 国产99白浆流出| 两人在一起打扑克的视频|