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

    The superiority of multi-trait models with genotype-by-environment interactions in a limited number of environments for genomic prediction in pigs

    2021-04-11 12:18:22HailiangSongQinZhangandXiangdongDing

    Hailiang Song,Qin Zhang and Xiangdong Ding*

    Abstract Background:Different production systems and climates could lead to genotype-by-environment(G×E)interactions between populations,and the inclusion of G×E interactions is becoming essential in breeding decisions.The objective of this study was to investigate the performance of multi-trait models in genomic prediction in a limited number of environments with G×E interactions.Results:In total,2,688 and 1,384 individuals with growth and reproduction phenotypes,respectively,from two Yorkshire pig populations with similar genetic backgrounds were genotyped with the PorcineSNP80 panel.Single-and multi-trait models with genomic best linear unbiased prediction(GBLUP)and BayesCπwere implemented to investigate their genomic prediction abilities with 20 replicates of five-fold cross-validation.Our results regarding between-environment genetic correlations of growth and reproductive traits(ranging from 0.618 to 0.723)indicated the existence of G×E interactions between these two Yorkshire pig populations.For single-trait models,genomic prediction with GBLUP was only 1.1%more accurate on average in the combined population than in single populations,and no significant improvements were obtained by BayesCπfor most traits.In addition,single-trait models with either GBLUP or BayesCπproduced greater bias for the combined population than for single populations.However,multi-trait models with GBLUP and BayesCπbetter accommodated G×E interactions,yielding 2.2%-3.8%and 1.0%-2.5%higher prediction accuracies for growth and reproductive traits,respectively,compared to those for single-trait models of single populations and the combined population.The multi-trait models also yielded lower bias and larger gains in the case of a small reference population.The smaller improvement in prediction accuracy and larger bias obtained by the single-trait models in the combined population was mainly due to the low consistency of linkage disequilibrium between the two populations,which also caused the BayesCπmethod to always produce the largest standard error in marker effect estimation for the combined population.Conclusions:In conclusion,our findings confirmed that directly combining populations to enlarge the reference population is not efficient in improving the accuracy of genomic prediction in the presence of G×E interactions,while multi-trait models perform better in a limited number of environments with G×E interactions.

    Keywords:Combined population,Genotype-by-environment interaction,Linkage disequilibrium,Multi-trait model,Pig

    Background

    Genomic selection[1],which relies on linkage disequilibrium(LD)between single nucleotide polymorphisms(SNPs)and causative variants,has become a useful tool in animal breeding[2]and plant breeding[3].Generally,the accuracy of genomic selection increases as the number of animals in the reference population increases[4].For small reference populations,combining populations of the same breed or related breeds reportedly increases the accuracy of genomic selection in cattle and pig[2,5-8].However,genomic prediction using combined populations has not shown a clear advantage over genomic prediction using single populations,possibly because the LD between SNPs and causative variants is not sufficiently consistent between populations[5-7].The failure to consider genotype-by-environment(G×E)interactions between populations is another important reason.Thus,exploiting G×E interactions in combined populations could be an attractive and meaningful approach for increasing the accuracy of genomic prediction.

    A G×E interaction is defined as different genotypes reacting unequally to environmental changes[9],and ignoring possible G×E interactions could lead to a reduction in genetic gains.Two models are widely used to detect G×E interactions.One is a multi-trait model,which assumes that the phenotypic expressions of a trait in various environments are different traits[10].In this case,the genetic correlation between traits in different environments is used as the indicator of a G×E interaction[9,11].The other model is the reaction norm model[12],which models the trajectory of animal performance as a function of an environmental gradient.However,in a small number of environments,e.g.,two or three environments,the reaction norm model captures only part of the G×E interaction due to the limited amount of environmental variation.Furthermore,it is difficult to define a suitable environmental covariate in the reaction norm model[12,13].

    With the development of genomic selection in pig breeding,more farms could be included to enlarge the reference population size in China,and joint genomic evaluation across countries or breeding organizations is expected.In dairy cattle,G×E interactions are usually ignored in combined-population genomic prediction;for example,Zhou et al.[14]reported that the consistency of LD was very high(0.97)between Chinese and Nordic Holsteins,indicating a high level of genetic similarity between the two populations.Therefore,it is necessary to consider the influence of G×E interactions,as the consistency of LD between pig populations is relatively low compared to that in dairy cattle.The objectives of this study were to evaluate G×E interactions between two Yorkshire pig populations with similar genetic backgrounds and to investigate the performance of multitrait models in genomic prediction in a limited number of environments with G×E interactions.

    Methods

    Ethics statement

    The whole procedure for collecting pig blood samples was carried out in strict accordance with the protocol approved by the Animal Welfare Committee of China Agricultural University(Permit Number:DK996).

    Population and phenotypes

    Yorkshire pig populations were sampled from two breeding farms in the north(Beijing)and south(Fujian Province)of China,designated Beijing and Fujian,respectively,for convenience.The pigs in the Beijing and Fujian populations are American Yorkshire progeny with the same genetic background but no pedigree links.Phenotypic records of reproductive traits,namely,the number of piglets born alive(NBA)and the total number of piglets born(TNB),and growth traits,namely,days to 100 kg(AGE)and backfat thickness at 100 kg(BFT),were examined.The measurements of AGE and BFT were reported in Song et al.[7].The phenotypic information is presented in Table 1.For the growth traits,phenotypic information from 28,827 and 13,860 pigs born in 2008-2017 and 2007-2018,respectively,in the Beijing and Fujian populations was available.For the reproductive traits,the farrow records of 5,968 and 3,115 pigs born in 2007-2017 and 2007-2018,respectively,in Beijing and Fujian were obtained.

    Construction of corrected phenotypes

    In pig genomic prediction analyses,corrected phenotypes derived from pedigree-based estimated breeding values(EBVs)are usually used as response variables.Conventional EBVs and genetic parameters were estimated based on a repeatability model of reproductive traits implemented separately in each population.In themodel,the fixed effects included herd,year and season,and the random effects included the additive genetic,random residual and permanent environmental effects.For the growth traits,a bivariate animal model was implemented,and the fixed effects included herd,year,season and sex.In addition to the additive genetic effect of each individual and the random residual effect,litter effect was also taken into account as a random effect.A total of 31,529 and 32,175 animals were traced back to construct a pedigree relationship matrix for the Beijing and Fujian populations,respectively.Afterwards,corrected phenotypic values(yc)for reproductive traits were calculated as the EBV plus the average of estimated residuals over parities for a sow,andycvalues for growth traits were computed as the EBV plus the estimated residual for each individual in each population.EBVs and genetic parameters were computed using the DMUAI procedure in DMU software[15].

    Table 1 Summary statistics for the two Yorkshire pig populations,including the numbers of genotyped animals and estimates of heritability(h2)

    Genotype data and quality control

    Genomic DNA was extracted from blood samples using a TIANamp Blood DNA Kit DP348(Tiangen,Beijing,China).Genotyping was performed using the PorcineSNP80 BeadChip(Illumina,San Diego,CA,USA),which includes 68,528 SNPs distributed across the entire pig genome.The number of genotyped animals for each trait is presented in Table 1.

    Missing genotypes for SNPs with known chromosomal positions were imputed by BEAGLE with default parameter settings[16],and those for SNPs with unknown positions were discarded.Based on the imputed dataset,SNPs were excluded from the analysis if the minor allele frequency was less than 0.05,the call frequency score was less than 0.90,or the genotype frequencies deviated from Hardy-Weinberg equilibrium with aP-value lower than 10-7.The animals with a call rate less than 0.90 or with an EBV reliability less than 0.3 were removed.After quality control,2,688 and 1,384 genotyped individuals remained for the growth and reproductive traits,respectively,and 56,445 SNPs were ultimately used.

    Principal component analysis

    Principal component analysis(PCA)was performed on the G matrix using GCTA software[17].This resulted in a matrix of eigenvectors in descending order that represented principal components(PCs),where PC1 had the largest eigenvalue.The overall structuring of genetic variation was visualized in a scatterplot of the top few PCs.

    Measure of linkage disequilibrium

    LD between a pair of SNPs was measured as[18],andwas calculated as follows:

    wheref(AB),f(A),f(a),f(B)andf(b)are the observed frequencies of haplotype AB and alleles A,a,B and b,respectively.The consistency of LD in the two populations was measured by the correlation ofrLDvalues of adjacent marker pairs on each autosome between the two populations.

    Model

    Four methods,namely,single-trait genomic BLUP(GBLUP),multi-trait GBLUP,single-trait BayesCπand multi-trait BayesCπ,were implemented to predict the genomic EBV(GEBV)for each genotyped individual.

    Single-trait GBLUP(ST-GBLUP)model

    The ST-GBLUP[19]model was used to predict the GEBVs of all genotyped individuals.

    where y is the vector of corrected phenotypic values;u is the overall mean;1 is a vector of 1;and a is the vector of genomic breeding values,following a normal distribution N(0,Gσa2),in whichσa2is the variance of the addictive genetic effect and G is the marker-based genomic relationship matrix[19].e is the vector of random errors,following a normal distribution N(0,Iσe2),in which σe2is the residual variance.

    Multi-trait GBLUP(MT-GBLUP)model

    The MT-GBLUP model was defined as.

    Single-trait BayesCπ(ST-BayesCπ)model

    In ST-BayesCπ[20],marker effects on phenotypic traits were sampled from a mixture of null and normal distributions,

    where yiis a vector of phenotypes,μis a vector of overall means,mijis the genotype covariate at locusjfor individuali(coded as 0,1 and 2),pis the number of genotyped loci,αjis a vector of allele substitution effects,and eiis a vector of random residuals for individualiwith a distributionThe markers in the model shared a common variance.The prior for the allele substitution effects of each molecular markerαjdepends on the varianceand the probabilityπthat markers do not have a genetic effect.

    Multi-trait BayesCπ(MT-BayesCπ)model

    In MT-BayesCπ,where each locus can have an effect on any combination of traits[21],the prior forαjk,the allele substitution or marker effect on traitkfor locusj,is a mixture with a point mass at zero and univariate normal distribution conditional on

    wheremij,jandpare the same as in ST-BayesCπ.yiis a vector of phenotypes ofktraits for individuali(corrected phenotypic values of the same trait in different populations),μis a vector of overall means forktraits,andαjis a vector of marker allele substitution effects onktraits for locusj.The residuals,ei,are a priori assumed to be independently and identically distributed multivariate normal vectors with a null mean and covariance matrix R(R is the same as described for the multi-trait GBLUP model),which,in turn,is assumed to have an inverse Wishart prior distribution,The covariance between effects for traitkandk′at the same locus,i.e.,αjkandαjk′is

    Imputation of missing phenotypic data was implemented in each Markov chain Monte Carlo(MCMC)iteration in MT-BayesCπ.

    For the ST-BayesCπand MT-BayesCπmethods,the MCMC chain was run for 20,000 cycles,and the first 10,000 cycles were discarded as burn-in.Every 10thsample of the remaining 10,000 iterations was saved to estimate SNP effects and the variance components.The Julia package JWAS was used for BayesCπanalyses[22],and the DMUAI procedure,implemented in DMU software,was used for GBLUP analyses.

    Cross-validation and prediction accuracy

    In this study,the accuracy and unbiasedness of prediction were obtained through 5-fold cross-validation(CV)for the Beijing and Fujian populations,where the Beijing(Fujian)population was randomly split into 5 folds,predicting one fold based on the other 4 folds,the combined population was divided into 4 folds and another population,or a multi-trait model was used in which the values of the same trait in 4 folds and another population were considered different traits.In each round of CV,phenotypes from one fold(validation population)were removed from the dataset,and the remaining folds(reference population)were used to predict the future performance of animals in the validation population.This 5-fold CV was replicated 20 times,and the results are presented as the mean and standard deviation for the 20 replicates.Here,the single-trait models(ST-GBLUP and ST-BayesCπ)based on single populations and the combined population were termed ST-GBLUP(BayesC π)_single and ST-GBLUP(BayesCπ)_combined,respectively.For ST-GBLUP(BayesCπ)_single,GBLUP(-BayesCπ)_combined and MT-GBLUP(BayesCπ),the validation population was the same for the four approaches in each replicate of 5-fold CV.The accuracy of genomic prediction was evaluated asr(yc,GEBV),the correlation between GEBVs and corrected phenotypic valuesycin the validation population.In addition,b(yc,GEBV),the regression ofycon GEBVs,was also calculated to assess the possible inflation or unbiasedness of predictions.A two-tailed paired t-test was used to compare prediction accuracy between a pair of scenarios.

    Results

    Population structure and genetic parameters

    To assess the population structure of the two Yorkshire pig populations,PCA was performed.Figure 1 illustrates that the genetic backgrounds of the Beijing and Fujian populations were similar,and both populations were composed of American Yorkshire progeny.In addition,the Yorkshire population(874 pigs)with a British origin was selected for PCA,which also showed that Beijing and Fujian had similar genetic backgrounds(Additional file:Fig.S1).However,no pedigree connection between these populations was detected according to their pedigree records due to a lack of genetic exchange.

    LD between a pair of SNPs was measured asr2LDandrLDin the Beijing and Fujian populations.The LD between adjacent markers in the Beijing and Fujian populations was also investigated,as shown in Table 2.The meanr2LDof adjacent SNP pairs within a chromosome ranged from 0.52 to 0.63 in the Beijing population and from 0.51 to 0.63 in the Fujian population.The meanr2LDacross all chromosomes was 0.56 in both the Beijing and Fujian populations,also indicating the similar genetic backgrounds of the populations.However,the consistency of LD between the two populations was not high.The correlation ofrLDfor adjacent SNPs between the two populations ranged from 0.47 to 0.60 across all chromosomes,with a mean of 0.55 at an average marker distance of 41 Kb,indicating the insufficient consistency in LD between these two pig populations.

    Estimates of heritability for the growth and reproductive traits in the two Yorkshire pig populations are shown in Table 1.In the Beijing population,the heritabilities of AGE and BFT were 0.33 and 0.34,respectively,with a standard error of 0.01.The heritabilities of AGE and BFT were 0.42 and 0.44,respectively,with a standard error of 0.02 in the Fujian population.For NBA and TNB,the estimated heritabilities were 0.07 and 0.08 with a standard error of 0.01 in the Beijing population and 0.09 and 0.11 with a standard error of 0.02 in the Fujian population.

    Fig.1 Principal component analysis(PCA)of two Yorkshire pig populations.Beijing and Fujian represent two Yorkshire pig populations from two Chinese pig breeding farms;PC1(3.5%)=first principal component(variance explained by PC1=3.5%);PC2(2.3%)=second principal component(variance explained by PC2=2.3%);PC3(1.3%)=third principal component(variance explained by PC3=1.3%)

    Table 2 Distance and LDof adjacent SNPs for each autosome(CHR)

    Table 2 Distance and LDof adjacent SNPs for each autosome(CHR)

    a Cor:the correlation of rLD values of adjacent SNP pairs between two populations;b Sum over 18 autosomesc Yorkshire pig populations from Beijing and Fujian with similar genetic backgrounds

    ?

    G×E interactions

    Figure 2 shows the estimated genetic correlations between the Beijing and Fujian populations based on the MT-GBLUP model using all genotyped animals.For AGE and BFT,the estimates of genetic correlations were 0.618 and 0.623,respectively,with standard errors of 0.145 and 0.134.The genetic correlations of NBA and TNB were 0.714 and 0.723,respectively,with standard errors of 0.153 and 0.159.These genetic correlations indicated that G×E interactions most likely exist between the Beijing and Fujian populations,as Robertson[23]suggested that 0.80 was the threshold of biological importance for G×E interactions.

    Tables 3 and 4 and Fig.3 demonstrate the accuracies of genomic prediction for growth and reproductive traits achieved by applying single-trait and multi-trait models with GBLUP and BayesCπ.Generally,in the same scenario,the prediction accuracies of each method were the largest for BFT,as it had the highest estimated heritability among the four traits.Lower prediction accuracies were obtained for two reproductive traits,NBA and TNB,due to their low heritabilities.For the two pig populations,higher accuracies for growth traits were acquired for the Beijing population,as it had a much larger reference population than the Fujian population(Table 1),while their small reference populations for reproductive traits resulted in comparably low accuracies of genomic prediction.Tables 3 and 4 and Fig.3 further show the superiority of the multi-trait model for genomic prediction when dealing with G×E interactions,as discussed below.

    Comparison of the combined population with the single populations

    Table 3 presents the accuracy of genomic prediction for growth and reproductive traits achieved by applying the GBLUP method,as measured by 20 replicates of fivefold CV.For the Beijing population,the prediction accuracies obtained with ST-GBLUP_single were 0.315 and 0.331 for AGE and BFT and 0.142 and 0.172 for NBA and TNB,respectively.When ST-GBLUP_combined was used,the genomic prediction accuracies were improved by 1.1% and 1.7% for NBA and TNB,respectively,compared with those obtained with ST-GBLUP_single,while there were no significant improvements for the growth traits AGE and BFT.In contrast,for the Fujian population,the genomic prediction accuracies for both growth and reproductive traits increased when the single populations were combined to enlarge the reference population.On average,ST_GBLUP_combined yielded 1.1% and 2.0% higher accuracies than STGBLUP_single for the growth and reproductive traits,respectively.Obviously,the gains in accuracy obtained for the Fujian population were larger.

    However,the single-trait model with the BayesCπ method did not yield higher prediction accuracies for the combined population compared with the single populations for most traits,even though the reference population was larger.As shown in Table 4,when predicting the traits of pigs from the Beijing population,the prediction accuracies obtained with STBayesCπ_single were 0.306,0.328,0.156 and 0.179 for AGE,BFT,NBA and TNB,respectively.STBayesCπ_combined produced almost the same accuracies as ST-BayesCπ_single in all scenarios,and the accuracy obtained by ST-BayesCπ_combined was even slightly lower in some cases;e.g.,the prediction accuracies of ST-BayesCπ_single and STBayesCπ_combined for NBA were 0.156 and 0.154,respectively.A trend of similar accuracies for STBayesCπbetween the combined and single populations was also found in the Fujian population,as demonstrated in Table 4.The average prediction accuracies obtained with ST-BayesCπ_single and ST-BayesCπ_combined were almost the same for the growth traits and reproductive traits.

    Fig.2 Genetic correlations between the two Yorkshire pig populations obtained using the multi-trait GBLUP method with all genotyped animals.AGE:days to 100 kg;BFT:backfat thickness at 100 kg;NBA:number of piglets born alive;TNB:total number of piglets born.The red line represents the threshold of 0.8.The error bar represents the standard error

    Table 3 Accuracy and unbiasedness of genomic prediction of growth and reproductive traits performed with the GBLUP method as assessed with 20 replicates of five-fold CV

    Table 4 Accuracy and unbiasedness of genomic prediction of growth and reproductive traits performed with the BayesCπmethod as assessed with 20 replicates of five-fold CV

    Fig.3 Comparison of genomic prediction accuracies of GBLUP and BayesCπmethods when implementing a single-trait model in single populations and a multi-trait model.a Predicting pigs from the Beijing population using a single-trait model;b predicting pigs from the Fujian population using a single-trait model;c predicting pigs from the Beijing population using a multi-trait model;d predicting pigs from the Fujian population using a multi-trait model.AGE:days to 100 kg;BFT:backfat thickness at 100 kg;NBA:number of piglets born alive;TNB:total number of piglets born;BayesCpi:BayesCπ

    The unbiasedness of genomic prediction of growth and reproductive traits as assessed with 20 replicates of five-fold CV is shown in Tables 3 and 4.When predicting the traits of pigs from the Beijing population based on the GBLUP method,larger bias(0.199)was produced by ST-GBLUP_combined than by ST-GBLUP_single.Likewise,unbiasedness was close to 1 with the STGBLUP_single method for growth and reproductive traits when predicting the traits of pigs from the Fujian population,while ST-GBLUP_combined yielded a 0.292 larger bias than ST-GBLUP_single(Table 3).Similarly,ST-BayesCπalso produced more bias for the combined population than for the single populations in most situations(Table 4).

    In addition,prediction accuracies obtained with training on the Beijing population and validation in the Fujian population and training on the Fujian population and validation in the Beijing population using STGBLUP and ST-BayesCπwere determined.Low prediction accuracies and large bias were observed in all scenarios(Additional file:Table S1),indicating that using one population to predict another is infeasible,even though the genetic backgrounds of the Beijing and Fujian populations were similar in this study.

    Comparison of the multi-trait model with the single-trait model

    Both single-trait and multi-trait models with GBLUP and BayesCπwere implemented to address the combined population in this study.Tables 3 and 4 show the accuracies and unbiasedness of genomic prediction obtained with the multi-trait and single-trait models in the same scenarios.In general,the multi-trait model showed significant superiority over the single-trait model.For the GBLUP method,as shown in Table 3,when predicting the traits of pigs from the Beijing population,the multi-trait GBLUP model(MT-GBLUP)yielded approximately 1.3%(1.2%)and 2.2%(0.8%)higher accuracies than ST-GBLUP_single(ST-GBLUP_combined)for the growth and reproductive traits,respectively.Furthermore,a small amount of bias was obtained with MTGBLUP in most scenarios.A similar trend was found in the Fujian population.As shown in Table 3,MT-GBLUP yielded 3.2%(2.2%)and 3.0%(1.1%)higher accuracies on average for the growth and reproductive traits,respectively,than ST-GBLUP_single(ST-GBLUP_combined).The unbiasedness was close to 1 when using the MTGBLUP methods for growth and reproductive traits in most situations.Again,the gains obtained by MTGBLUP in the Fujian population were larger than those in the Beijing population.

    For the multi-trait model with BayesCπ,as shown in Table 4,when predicting the traits of pigs from the Beijing population,MT-BayesCπyielded approximately 1.3%(1.4%)and 2.5%(2.5%)higher accuracies than STBayesCπ_single(ST-BayesCπ_combined)for the growth and reproductive traits,respectively.The largest gain over ST-BayesCπ_single and ST-BayesCπ_combined was 2.7% and 2.9% for NBA,respectively.When predicting the traits of pigs from the Fujian population,MT-BayesCπobtained higher accuracies for all traits(Table 4),yielding 2.7%(2.6%)and 2.4%(2.3%)gains on average over ST-BayesCπ_single(ST-BayesCπ_combined)for the growth and reproductive traits,respectively.MT-BayesCπalso produced bias as low as that of ST-BayesCπ_single in most scenarios.

    Comparison of the GBLUP and BayesCπmethods

    Figure 3 presents the accuracy of genomic prediction of growth and reproductive traits achieved by applying single-population ST-GBLUP(BayesCπ)and MTGBLUP(BayesCπ)methods.When predicting the traits of pigs from the Beijing population,GBLUP yielded approximately 1.0% higher accuracy than BayesCπregardless of whether the single-or multi-trait model was used for AGE and a bias close to 0 with the GBLUP method in this scenario,while there was no difference in accuracy between GBLUP and BayesCπfor BFT.For reproductive traits in the Beijing population,the accuracy of genomic prediction was increased by 1.3% on average for MT-BayesCπcompared to MT-GBLUP,and no difference in accuracy was found between ST-GBLUP and ST-BayesCπ.When predicting the traits of pigs from the Fujian population,GBLUP and BayesCπproduced similar prediction accuracies for AGE;e.g.,the prediction accuracies with ST-GBLUP_single and ST-BayesCπ_single were 0.245 and 0.243,respectively,while less bias was obtained by GBLUP in this scenario.There was no difference in prediction accuracy or unbiasedness between ST-GBLUP_single and ST-BayesCπ_single for BFT,while MT-GBLUP achieved a 1.3% higher accuracy than MT-BayesCπ.For reproductive traits in the Fujian population,there was no difference in prediction accuracy between GBLUP and BayesCπ,except that MTGBLUP_combined achieved a 1.9% higher accuracy than MT-BayesCπ_combined for TNB.Generally,GBLUP performed comparably to BayesCπin the single-trait and multi-trait models.

    Discussion

    In this study,we investigated the accuracy of genomic prediction in two Yorkshire pig populations with similar genetic backgrounds when G×E interactions were taken into account.Our results revealed G×E interactions between the two populations for the reference growth and reproductive traits.We further explored whether directly combining populations to enlarge the reference population is efficient in improving the accuracy of genomic prediction in the presence of G×E interactions.Our results demonstrated the superiority of the multi-trait model for genomic prediction in dealing with G×E interactions,which could better accommodate the G×E interactions.Higher prediction accuracies and lower bias were obtained with the multi-trait model than with the single-trait model when implementing either GBLUP or BayesCπin a given situation.

    G×E interactions play an important role in pig populations,and they should be considered in breeding programs to select the best animals in different environments[11,24].The detection of G×E interactions relies on a genetic correlation of 0.8 for one trait in different environments,the threshold suggested by Robertson[23].Accordingly,G×E interactions between two pig populations for growth and reproductive traits(ranging from 0.618 to 0.723)were observed.Similarly,Li et al.[25]reported on the across-country genetic correlations(ranging from 0.604 to 0.726)of dairy cattle,indicating that an important G×E interaction exists between Brazilian and Nordic(or Nordic and French)populations.Hay and Roberts[26]also reported the existence of genotype-by-nutritional environment interactions for growth and carcass traits in beef cattle with genetic correlations below the threshold value of 0.8.In this study,SNP markers were used to construct the relationships(G matrix)between the two populations in order to estimate the genetic correlations,as no linked pedigree was available.Compared with a pedigree used to estimate genetic correlations,the realized relationships among individuals are captured by marker information and achieve more accurate estimates of genetic correlations.However,the number of genotyped animals was not large enough to produce a lower standard error for genetic correlations compared to that obtained by estimating heritability.We therefore used SNP markers to assess the population structure of the two pig populations,and the similar patterns of LD on each chromosome between these two populations further showed that they had similar genetic backgrounds.However,the correlation ofrLDvalues of adjacent SNPs with a mean of 0.55 indicated the insufficient consistency in LD between these two pig populations,which could be due to G×E interactions.Genetic correlations are caused by the pleiotropic action of genes and linkage between genes affecting different traits[9].If the performance values of the same trait in different populations(environments)are regarded as different traits and are affected by different genes(that are of course pleiotropic),the phase of linkage of these genes and the consistency of their linkage(not LD)in the two populations will affect the strength of the genetic correlation and G×E interaction.

    The size of the reference population is one key factor in genomic prediction[4].Combining populations is an easy but practical way to improve the accuracy of genomic prediction,and its advantages in practice have been reported by EuroGenomics[8]and North American[2]consortiums.Similar results were shown in this study,in which higher prediction accuracies were obtained by the single-trait model for the combined population than for the single populations in most situations(Table 3),and such gains were greater for the population with a small reference size;e.g.,the prediction accuracy in the combined population was 1.1% higher than that in the single population for the growth trait BFT in the Fujian population but only 0.3% higher in the Beijing population,and the number of genotyped animals in the Beijing population was much larger than that in the Fujian population(Table 2),in accordance with the features of other studies[5,14].However,in some situations,accuracy was not increased or slightly decreased in the combined population;e.g.,the prediction accuracies for AGE in the Beijing population obtained with ST-BayesCπ_combined and ST-BayesCπ_single were 0.306 and 0.304,respectively(Table 4).Generally,in this study,the single-trait model exhibited only a slight gain in genomic prediction accuracy for most traits when using the combined population.In addition,the single-trait model produced more bias for the combined population than for the single populations(Tables 3 and 4).The lower improvement or slight decrease in the combined population compared with the expectation could be due to the weak consistency in LD between populations,as pointed out by Lund et al.[8].

    Many studies have shown that inconsistent LD between SNPs and causative variants across populations causes disadvantages compared with using a combined population in genomic prediction[5,6,27,28].In this study,the correlation ofrLDvalues of adjacent SNPs with a mean of 0.55 indicated the insufficient consistency in LD between the Beijing and Fujian populations.Weak consistency in LD will lead to bias in SNP effect estimates when populations are simply combined with a single-trait model,further resulting in fewer improvements in accuracy and larger bias.Therefore,the effect of increasing the consistency in LD or reducing its noise on GEBVs is worth investigating.On the one hand,the consistency in LD between populations could be affected by the density of chip SNPs.De Roos et al.[29]found that markers in LD with causative variants across diverged breeds would require~300,000 SNPs.Hoze et al.[30]reported that 2% higher prediction accuracies were acquired for combined small populations of dairy cattle breeds when using the 777 K panel compared to the 54 K panel,and our previous study showed that imputed whole-genome sequencing data hold the potential to increase the accuracy of genomic prediction for combined populations in pigs[7].

    On the other hand,reasonable models should be explored to take advantage of the correlation of LD between populations.The single-trait model ignores the inconsistency in LD between populations,while the multi-trait model treats populations as different environments and uses covariance to consider the correlation between them,which is caused by the consistency(no matter how high or low)in LD between populations.In this study,the multi-trait model performed best in all scenarios(Tables 3 and 4).Similar results were obtained in a Brazilian Holstein population by adding data from Nordic and French Holstein populations and using a multi-trait model[25].In addition,our results also showed that less bias was obtained with the multi-trait model compared to the single-trait model,indicating that the multi-trait model can not only use the combined information of populations to increase prediction accuracy but also reduce the impact of LD inconsistencies.

    The multi-trait model and reaction norm model are two prevalent methods for handling G×E interactions in genetic evaluations.Only a few studies have explicated the advantages and disadvantages of these two kinds of methods.Falconer et al.[9]reported that G×E interactions could be detected using a multi-trait model,in which the trait values in each environment are treated as genetically distinct traits.In a limited number of environments(e.g.,in this study,with two breeding farms),a multi-trait model could capture the G×E interaction between these environments.However,the computational demand of multi-trait models will increase rapidly with an increase in the number of environmental levels,as more(co)variance components will have to be estimated and convergence will become increasingly difficult.In contrast to the multitrait model,the reaction norm model can handle a large number of or continuous environmental levels(e.g.,farm,year and seasonal effects),and genetic parameters could be estimated at each environmental level[12].However,the reaction norm model captures only part of the G×E interaction due to the limited number of environmental levels.In this study,the genomic prediction accuracies were not improved by using the reaction norm model with farm,year and seasonal effects as environmental covariates due to the limited number of environments(results not shown).Therefore,the reaction norm model could not perform well in a limited number of environments.In addition,convergence of the reaction norm model could be hindered because of the complexity of the model when the data are less informative[11].Therefore,a multi-trait model is optimal for analyzing G×E interactions in a limited number of environments.Furthermore,the multitrait model is more convenient in practical breeding,and it can accommodate phenotypes that are not measured in exactly the same way,e.g.,different scales of deregressed proofs and the existence of G×E interactions in Chinese and Nordic Holsteins[25,31].Such phenotypes cannot be analyzed by a single-trait model for joint genomic prediction,while the use of a multi-trait model is plausible in such cases.

    In this study,GBLUP and BayesCπ,assuming a common(co)variance for all(GBLUP)or a proportion(BayesC π)of the SNP effects,were used in this study.As many studies have reported for real data[32,33],the average prediction accuracies were not significantly different between GBLUP and BayesCπin the single populations and multi-trait model in this study(Fig.3).However,BayesCπ performed worse than GBLUP for the combined population,and no improvements in accuracy were obtained by BayesCπfor most traits compared to that obtained with a single population.In addition,BayesCπproduced a larger bias(Table 4).Table 5 provides the standard error of the marker effect estimated by BayesCπin single and combined populations.The BayesCπmethod for the combined population produced the largest standard error of the marker effect in all situations,implying larger bias when estimating genomic breeding values.A larger standard error for the marker effect might have been caused by the inconsistency in LD between the Beijing and Fujian populations.In contrast,the minimum standard error of the marker effects estimated by the multi-trait BayesCπ model also indicated the superiority of multi-trait BayesC πin genomic prediction.Multi-trait BayesCπwas first proposed by Jia and Jannink[34],with the restrictive assumption that a locus simultaneously affects all the traits or none of them.In this study,multi-trait BayesCπ allowed a locus to affect any combination of traits,which is biologically realistic,especially in multi-trait analyses involving many traits[21].However,it might not be possible to apply multi-trait single-step GBLUP(BayesCπ)[35-37]in this study,as no pedigree link between the two pig populations was available and the BayesCπmethod was too computationally demanding.

    Table 5 Mean value of the standard error of marker effects estimated by the BayesCπmethod using all genotyped animals

    Conclusions

    G×E interactions constitute a potential source of inefficiency in animal breeding.Ignoring possible G×E interactions could lead to a reduction in genetic gains.Our results demonstrated that directly combining populations to enlarge the reference population is not efficient in improving the accuracy of genomic prediction in the presence of G×E interactions in two Yorkshire pig populations,while the multi-trait model performed better in a limited number of environments with weak G×E interactions.This study will be helpful when leveraging G×E interactions in practical genomic selection.

    Supplementary information

    Supplementary informationaccompanies this paper at https://doi.org/10.1186/s40104-020-00493-8.

    Additional file 1:Fig.S1.Principal component analysis(PCA)of the Beijing and Fujian Yorkshire pig populations and a British Yorkshire pig population.

    Additional file 2:Table S1.Accuracy and unbiasedness of genomic prediction of growth and reproductive traits obtained when using one population to predict another.

    Abbreviations

    G×E:Genotype-by-environment;GBLUP:Genomic best linear unbiased prediction;LD:Linkage disequilibrium;SNPs:Single nucleotide polymorphisms;Beijing and Fujian represent two Yorkshire pig populations from two Chinese pig breeding farms;AGE:Days to 100 kg;BFT:Backfat thickness at 100 kg;NBA:Number of piglets born alive;TNB:Total number of piglets born;yc:Corrected phenotypic values;GEBV:Genomic EBV;STGBLUP(BayesCπ):Single-trait GBLUP(BayesCπ);MT-GBLUP(BayesCπ):Multitrait GBLUP(BayesCπ);CV:Cross validation;ST-GBLUP(BayesCπ)_single:STGBLUP(BayesCπ)based on single population;ST-GBLUP(BayesC π)_combined:ST-GBLUP(BayesCπ)based on combined population;PCA:Principal component analysis;PC:Principal component;MCMC:Markov chain Monte Carlo

    Acknowledgements

    The authors are grateful to Beijing Liuma and Fujian Yongcheng Pig Breeding Co.,Ltd.for providing blood samples.

    Authors’contributions

    HLS and XDD designed the experiments.HLS performed statistical analysis and wrote the manuscript.QZ and XDD revised the manuscript.All authors read and approved the final manuscript.

    Funding

    This work was supported by grants from the earmarked fund for China Agriculture Research System(CARS-35),Modern Agriculture Science and Technology Key Project of Hebei Province(19226376D),the National Key Research and Development Project(SQ2019YFE00771),the National Natural Science Foundation of China(31671327).Major Project of Selection for New Livestock and Poultry Breeds of Zhejiang Province(2016C02054-5).

    Availability of data and materials

    The datasets used during the current study are available from the corresponding author on reasonable request.

    Ethics approval and consent to participate

    The experiments involving animals were approved by the Animal Welfare Committee of China Agricultural University.There was no use of human participants,data or tissues.

    Consent for publication

    Not applicable.

    Competing interests

    The authors declare that they have no competing interests.

    Author details

    1National Engineering Laboratory for Animal Breeding,Laboratory of Animal Genetics,Breeding and Reproduction,Ministry of Agriculture,College of Animal Science and Technology,China Agricultural University,Beijing 100193,China.2Shandong Provincial Key Laboratory of Animal Biotechnology and Disease Control and Prevention,Shandong Agricultural University,Taian 271001,China.

    Received:1 March 2020 Accepted:7 July 2020

    啦啦啦韩国在线观看视频| 中出人妻视频一区二区| 长腿黑丝高跟| 五月伊人婷婷丁香| av中文乱码字幕在线| 你懂的网址亚洲精品在线观看 | 亚洲电影在线观看av| 日韩欧美在线乱码| 亚洲欧美日韩高清在线视频| 色视频www国产| 禁无遮挡网站| 少妇猛男粗大的猛烈进出视频 | 午夜精品一区二区三区免费看| 午夜福利成人在线免费观看| 色视频www国产| 在线免费观看不下载黄p国产| 久久天躁狠狠躁夜夜2o2o| 偷拍熟女少妇极品色| 国产精品久久久久久精品电影| 亚洲成人精品中文字幕电影| 波多野结衣高清无吗| av.在线天堂| 波野结衣二区三区在线| 人人妻人人看人人澡| 最近最新中文字幕大全电影3| 人妻久久中文字幕网| 女的被弄到高潮叫床怎么办| 九九爱精品视频在线观看| 国产亚洲精品久久久com| 国内精品宾馆在线| 在现免费观看毛片| 亚洲四区av| 国产精品一及| 欧美一区二区精品小视频在线| 国产高清不卡午夜福利| 在线观看免费视频日本深夜| 婷婷亚洲欧美| 人人妻人人澡人人爽人人夜夜 | 日韩亚洲欧美综合| 俄罗斯特黄特色一大片| 97在线视频观看| 亚洲最大成人av| 69av精品久久久久久| 亚洲欧美日韩卡通动漫| 亚洲国产精品成人综合色| 我要搜黄色片| 一级毛片久久久久久久久女| 乱人视频在线观看| 亚洲精品粉嫩美女一区| 亚洲熟妇中文字幕五十中出| 99热只有精品国产| 日韩中字成人| 亚洲成人中文字幕在线播放| 亚洲成人精品中文字幕电影| 91av网一区二区| 日本色播在线视频| 午夜影院日韩av| 你懂的网址亚洲精品在线观看 | 一区二区三区免费毛片| 欧美成人一区二区免费高清观看| 成人午夜高清在线视频| 欧美zozozo另类| 成人无遮挡网站| 好男人在线观看高清免费视频| 人妻少妇偷人精品九色| 在线免费观看不下载黄p国产| 国产69精品久久久久777片| 一级av片app| 老司机午夜福利在线观看视频| 啦啦啦啦在线视频资源| 亚洲欧美精品综合久久99| 国产高清三级在线| .国产精品久久| 中文字幕久久专区| 亚洲不卡免费看| 国产在线精品亚洲第一网站| 亚洲av成人精品一区久久| 日韩在线高清观看一区二区三区| 精品不卡国产一区二区三区| 我的老师免费观看完整版| 蜜臀久久99精品久久宅男| 3wmmmm亚洲av在线观看| 亚洲一级一片aⅴ在线观看| 国产精品一及| 精品久久久久久久久av| 国产一区二区三区在线臀色熟女| 91久久精品国产一区二区成人| 一区二区三区免费毛片| 春色校园在线视频观看| 亚洲欧美日韩高清专用| 一进一出抽搐动态| 欧美色欧美亚洲另类二区| 国产单亲对白刺激| 欧美成人免费av一区二区三区| 在线播放国产精品三级| 亚洲欧美成人精品一区二区| h日本视频在线播放| 亚洲电影在线观看av| 天天躁夜夜躁狠狠久久av| eeuss影院久久| 国产不卡一卡二| 久久久久久久久大av| 日韩强制内射视频| 偷拍熟女少妇极品色| 99久久无色码亚洲精品果冻| 国产淫片久久久久久久久| 欧美中文日本在线观看视频| 校园春色视频在线观看| 国产高清三级在线| 中文资源天堂在线| 波野结衣二区三区在线| 久久热精品热| 国产一区二区三区av在线 | 婷婷色综合大香蕉| 日韩三级伦理在线观看| 激情 狠狠 欧美| 午夜老司机福利剧场| 在线看三级毛片| 九九爱精品视频在线观看| 久久久a久久爽久久v久久| 赤兔流量卡办理| 婷婷精品国产亚洲av| 免费观看在线日韩| 亚洲自拍偷在线| 久久精品国产亚洲av香蕉五月| a级毛片免费高清观看在线播放| 美女内射精品一级片tv| 日韩亚洲欧美综合| 国产中年淑女户外野战色| 欧洲精品卡2卡3卡4卡5卡区| 国产又黄又爽又无遮挡在线| 亚洲成av人片在线播放无| 又黄又爽又免费观看的视频| 欧美人与善性xxx| 欧美日韩一区二区视频在线观看视频在线 | 夜夜爽天天搞| 亚洲欧美清纯卡通| 麻豆成人午夜福利视频| 女人被狂操c到高潮| 成人亚洲精品av一区二区| 免费观看精品视频网站| 午夜日韩欧美国产| 成人亚洲欧美一区二区av| 亚洲激情五月婷婷啪啪| 男女下面进入的视频免费午夜| 国产av不卡久久| av视频在线观看入口| 久久久久九九精品影院| 一级av片app| 日本一本二区三区精品| 黄色配什么色好看| 亚洲激情五月婷婷啪啪| 美女大奶头视频| 国产精品三级大全| 天美传媒精品一区二区| 国产精品,欧美在线| 久久久精品94久久精品| a级毛片a级免费在线| 一区二区三区四区激情视频 | 最后的刺客免费高清国语| 乱系列少妇在线播放| 欧美高清性xxxxhd video| 成人精品一区二区免费| av福利片在线观看| 国产男靠女视频免费网站| 99热只有精品国产| 一夜夜www| 99在线视频只有这里精品首页| 亚洲精品乱码久久久v下载方式| av专区在线播放| 欧美日韩综合久久久久久| 亚洲激情五月婷婷啪啪| 免费不卡的大黄色大毛片视频在线观看 | 一区福利在线观看| 亚洲精品日韩在线中文字幕 | 我的老师免费观看完整版| 深夜a级毛片| 男女边吃奶边做爰视频| 久久天躁狠狠躁夜夜2o2o| 久99久视频精品免费| 国产精品永久免费网站| 亚洲人成网站在线观看播放| 国产视频一区二区在线看| 天天躁夜夜躁狠狠久久av| 精品一区二区免费观看| 99九九线精品视频在线观看视频| 国产精品久久久久久av不卡| 亚洲欧美日韩高清专用| 老女人水多毛片| 国产伦一二天堂av在线观看| 国产成人精品久久久久久| 热99在线观看视频| 久久久久久大精品| 午夜精品在线福利| 99久久成人亚洲精品观看| 精品国内亚洲2022精品成人| 熟妇人妻久久中文字幕3abv| 一个人看的www免费观看视频| 人人妻,人人澡人人爽秒播| 99热网站在线观看| 精品久久国产蜜桃| 天天躁日日操中文字幕| 夜夜夜夜夜久久久久| 免费人成在线观看视频色| 婷婷精品国产亚洲av| 91av网一区二区| 久久久精品94久久精品| 国产亚洲av嫩草精品影院| 校园春色视频在线观看| 内射极品少妇av片p| 全区人妻精品视频| 久久精品国产自在天天线| 十八禁网站免费在线| 亚洲七黄色美女视频| 亚洲四区av| 久久精品国产亚洲网站| 精品久久久久久久人妻蜜臀av| 99视频精品全部免费 在线| 欧美日韩精品成人综合77777| 久久久成人免费电影| 国产精品日韩av在线免费观看| 六月丁香七月| 99久国产av精品| 91在线观看av| 日本-黄色视频高清免费观看| 五月伊人婷婷丁香| 有码 亚洲区| 毛片女人毛片| 色av中文字幕| 国产免费男女视频| 国产一级毛片七仙女欲春2| 人人妻人人澡欧美一区二区| 久久精品国产清高在天天线| 国产成人影院久久av| 日本三级黄在线观看| 免费无遮挡裸体视频| eeuss影院久久| 国产亚洲91精品色在线| av卡一久久| 欧美丝袜亚洲另类| 欧美一级a爱片免费观看看| 欧美+日韩+精品| 国产极品精品免费视频能看的| 少妇丰满av| 久久久精品欧美日韩精品| 简卡轻食公司| 国产在线精品亚洲第一网站| 久久精品91蜜桃| 国产高清视频在线观看网站| 精品久久久久久久末码| 国产v大片淫在线免费观看| 青春草视频在线免费观看| 国产亚洲精品久久久久久毛片| 欧美日韩精品成人综合77777| 欧美性感艳星| 女人被狂操c到高潮| 国产成人freesex在线 | 国产黄色视频一区二区在线观看 | 91久久精品国产一区二区三区| 性色avwww在线观看| 校园春色视频在线观看| 婷婷精品国产亚洲av| 久久精品夜夜夜夜夜久久蜜豆| 久久综合国产亚洲精品| 热99re8久久精品国产| 99热这里只有是精品50| 一个人看的www免费观看视频| 国产高清激情床上av| 观看美女的网站| 成人鲁丝片一二三区免费| 18+在线观看网站| 男人舔奶头视频| 一级毛片电影观看 | 成人精品一区二区免费| 悠悠久久av| 亚洲精品色激情综合| 99在线视频只有这里精品首页| 国产色爽女视频免费观看| 精品熟女少妇av免费看| 成人国产麻豆网| av卡一久久| 精品乱码久久久久久99久播| 日韩一区二区视频免费看| 一个人看的www免费观看视频| 久久久久久伊人网av| 男人的好看免费观看在线视频| 亚洲天堂国产精品一区在线| 美女被艹到高潮喷水动态| 搡老熟女国产l中国老女人| 精品一区二区三区人妻视频| АⅤ资源中文在线天堂| 韩国av在线不卡| 欧美日韩国产亚洲二区| 欧美日韩在线观看h| 一区二区三区四区激情视频 | 国产真实乱freesex| 国产色婷婷99| 在线观看av片永久免费下载| 乱人视频在线观看| 国产视频内射| 欧美人与善性xxx| 搡老岳熟女国产| 岛国在线免费视频观看| 免费黄网站久久成人精品| 色综合色国产| 一区福利在线观看| 日日干狠狠操夜夜爽| 色尼玛亚洲综合影院| 大香蕉久久网| 久久精品国产亚洲av香蕉五月| 禁无遮挡网站| 欧美最黄视频在线播放免费| 成年女人看的毛片在线观看| 亚洲熟妇中文字幕五十中出| 别揉我奶头~嗯~啊~动态视频| 日日摸夜夜添夜夜添小说| 97碰自拍视频| 午夜亚洲福利在线播放| 色视频www国产| 久久久久国产精品人妻aⅴ院| 亚洲精品成人久久久久久| 国产黄色小视频在线观看| 亚洲美女搞黄在线观看 | 亚洲aⅴ乱码一区二区在线播放| 老司机福利观看| 91狼人影院| 久久精品国产自在天天线| 国模一区二区三区四区视频| 色尼玛亚洲综合影院| 插逼视频在线观看| 国产国拍精品亚洲av在线观看| 99在线视频只有这里精品首页| 欧美最黄视频在线播放免费| 国产av在哪里看| 成人午夜高清在线视频| 欧美国产日韩亚洲一区| 欧美区成人在线视频| 亚洲国产欧美人成| 久久精品国产99精品国产亚洲性色| 白带黄色成豆腐渣| 久久中文看片网| 国产精品美女特级片免费视频播放器| 国内精品久久久久精免费| 国产男靠女视频免费网站| 亚洲av熟女| 真实男女啪啪啪动态图| 国产爱豆传媒在线观看| 深爱激情五月婷婷| 日韩欧美国产在线观看| 国产白丝娇喘喷水9色精品| 国产精品一及| 午夜免费激情av| 可以在线观看毛片的网站| 嫩草影视91久久| 天美传媒精品一区二区| 尤物成人国产欧美一区二区三区| 成人毛片a级毛片在线播放| 嫩草影院精品99| 国产精品伦人一区二区| 亚洲成人av在线免费| 日本精品一区二区三区蜜桃| 成人av在线播放网站| 搡老熟女国产l中国老女人| 菩萨蛮人人尽说江南好唐韦庄 | 国产激情偷乱视频一区二区| 神马国产精品三级电影在线观看| 超碰av人人做人人爽久久| 成年av动漫网址| 寂寞人妻少妇视频99o| 99九九线精品视频在线观看视频| 晚上一个人看的免费电影| 成人亚洲欧美一区二区av| 亚洲国产欧美人成| 晚上一个人看的免费电影| 寂寞人妻少妇视频99o| 日本a在线网址| 在线播放国产精品三级| 成人特级黄色片久久久久久久| 欧美性感艳星| 十八禁国产超污无遮挡网站| 欧美又色又爽又黄视频| 99久国产av精品国产电影| 久久久久国产精品人妻aⅴ院| 少妇的逼好多水| 精品久久久久久成人av| av在线观看视频网站免费| 嫩草影院入口| 国产精品野战在线观看| 3wmmmm亚洲av在线观看| 精品久久久久久久末码| 国产精品人妻久久久久久| 久久精品国产清高在天天线| 国产精品嫩草影院av在线观看| 久久久欧美国产精品| 久久国产乱子免费精品| 99久国产av精品国产电影| 国产精品99久久久久久久久| 免费av毛片视频| 在线看三级毛片| 亚洲中文字幕日韩| 18禁在线播放成人免费| 观看免费一级毛片| 国产一区二区三区av在线 | 久久午夜亚洲精品久久| 久久鲁丝午夜福利片| 床上黄色一级片| 色尼玛亚洲综合影院| 国产男靠女视频免费网站| 国产精华一区二区三区| 亚洲欧美日韩东京热| 长腿黑丝高跟| 99久国产av精品| 此物有八面人人有两片| 老熟妇乱子伦视频在线观看| 国产黄片美女视频| 国产一区二区亚洲精品在线观看| 午夜精品国产一区二区电影 | 老司机影院成人| 国产毛片a区久久久久| 女生性感内裤真人,穿戴方法视频| 男人狂女人下面高潮的视频| 亚洲国产高清在线一区二区三| 深爱激情五月婷婷| 淫妇啪啪啪对白视频| 久久人人精品亚洲av| 免费观看的影片在线观看| 色5月婷婷丁香| 久久人人爽人人爽人人片va| 九九爱精品视频在线观看| 午夜福利在线在线| 国产亚洲精品久久久com| 亚洲,欧美,日韩| 免费高清视频大片| 精品日产1卡2卡| 国产爱豆传媒在线观看| 高清毛片免费观看视频网站| 不卡视频在线观看欧美| 欧美一区二区国产精品久久精品| 日本三级黄在线观看| 亚洲真实伦在线观看| 国产欧美日韩一区二区精品| .国产精品久久| 亚洲av一区综合| 欧美最黄视频在线播放免费| 国产高清激情床上av| aaaaa片日本免费| 精品一区二区三区视频在线观看免费| 日日干狠狠操夜夜爽| 日韩欧美在线乱码| 国产一区亚洲一区在线观看| 欧美成人精品欧美一级黄| 国产精品综合久久久久久久免费| 18禁在线无遮挡免费观看视频 | 国产精品爽爽va在线观看网站| 国产色爽女视频免费观看| 午夜福利18| 少妇人妻精品综合一区二区 | 久久久欧美国产精品| 国产精品人妻久久久影院| 蜜桃久久精品国产亚洲av| 亚洲精品在线观看二区| 亚洲最大成人手机在线| 又黄又爽又免费观看的视频| 欧美一区二区亚洲| a级一级毛片免费在线观看| 国产精品久久电影中文字幕| 日韩三级伦理在线观看| 国内精品宾馆在线| 精品人妻一区二区三区麻豆 | 嫩草影视91久久| 春色校园在线视频观看| 国产精品一区二区三区四区免费观看 | 九九热线精品视视频播放| 久久久久精品国产欧美久久久| 日本色播在线视频| 可以在线观看的亚洲视频| 午夜福利在线观看免费完整高清在 | 日本五十路高清| 在线免费十八禁| 人妻夜夜爽99麻豆av| 一级毛片我不卡| 国产黄片美女视频| 亚洲欧美日韩卡通动漫| 国产成人freesex在线 | 亚洲无线观看免费| 麻豆国产97在线/欧美| 天堂网av新在线| 久久久久久久久久黄片| 一本精品99久久精品77| avwww免费| 色av中文字幕| 日本五十路高清| 好男人在线观看高清免费视频| 熟女电影av网| 高清毛片免费观看视频网站| 变态另类成人亚洲欧美熟女| 少妇裸体淫交视频免费看高清| 联通29元200g的流量卡| 国产av不卡久久| 天堂网av新在线| 丰满人妻一区二区三区视频av| 久久精品国产亚洲av天美| 精华霜和精华液先用哪个| 国产单亲对白刺激| 日本与韩国留学比较| 国产亚洲欧美98| 97碰自拍视频| 国产成人aa在线观看| 国产爱豆传媒在线观看| 99在线视频只有这里精品首页| 日本-黄色视频高清免费观看| 久久午夜福利片| 亚洲专区国产一区二区| 免费人成视频x8x8入口观看| 久久久久国产网址| 中国国产av一级| 在线国产一区二区在线| 国产一区二区在线av高清观看| 人妻丰满熟妇av一区二区三区| 亚洲熟妇中文字幕五十中出| 淫妇啪啪啪对白视频| 精品午夜福利在线看| 黄色欧美视频在线观看| 亚洲18禁久久av| 中文字幕久久专区| 亚洲精品色激情综合| 99热这里只有是精品50| 亚州av有码| 亚洲人成网站高清观看| 亚洲自偷自拍三级| 99热这里只有是精品50| 国产精品国产三级国产av玫瑰| 国产久久久一区二区三区| 精品国产三级普通话版| 精品日产1卡2卡| 深夜a级毛片| 极品教师在线视频| av国产免费在线观看| 一级av片app| eeuss影院久久| 精品久久久久久成人av| 国内揄拍国产精品人妻在线| 精品午夜福利视频在线观看一区| 国产精品福利在线免费观看| 国产亚洲欧美98| 国产精品无大码| 国产单亲对白刺激| 午夜精品一区二区三区免费看| 搡女人真爽免费视频火全软件 | 天天一区二区日本电影三级| 两个人视频免费观看高清| 黄色一级大片看看| 91在线精品国自产拍蜜月| 国产伦一二天堂av在线观看| 免费看美女性在线毛片视频| 国产精品99久久久久久久久| 在线观看美女被高潮喷水网站| av在线亚洲专区| 在线观看美女被高潮喷水网站| 成人高潮视频无遮挡免费网站| 免费人成视频x8x8入口观看| 波野结衣二区三区在线| 亚洲精品久久国产高清桃花| 国产高清不卡午夜福利| 夜夜爽天天搞| 在线观看免费视频日本深夜| 亚洲婷婷狠狠爱综合网| 亚洲专区国产一区二区| 男人的好看免费观看在线视频| 午夜福利在线观看免费完整高清在 | 国产一区二区在线观看日韩| 九九爱精品视频在线观看| 国产成人91sexporn| 男女那种视频在线观看| 男女边吃奶边做爰视频| 91狼人影院| 一区二区三区四区激情视频 | 在线a可以看的网站| av卡一久久| 日本黄色视频三级网站网址| 欧美色视频一区免费| 精品人妻偷拍中文字幕| 午夜激情福利司机影院| 亚洲成人精品中文字幕电影| 成熟少妇高潮喷水视频| 日韩精品有码人妻一区| 日本在线视频免费播放| 亚洲自偷自拍三级| 看片在线看免费视频| 噜噜噜噜噜久久久久久91| 国产免费男女视频| 少妇猛男粗大的猛烈进出视频 | 国产单亲对白刺激| 国产不卡一卡二| 亚洲图色成人| 亚洲国产欧洲综合997久久,| 亚洲人成网站在线观看播放| 免费观看人在逋| 三级经典国产精品| 春色校园在线视频观看| 天天躁夜夜躁狠狠久久av| 啦啦啦啦在线视频资源| 男女做爰动态图高潮gif福利片| 十八禁国产超污无遮挡网站| 有码 亚洲区| 在线播放无遮挡| 嫩草影院新地址| 亚洲国产精品合色在线| 18禁裸乳无遮挡免费网站照片| av在线蜜桃| 亚洲成a人片在线一区二区| 久久精品国产自在天天线| 老女人水多毛片| 国产成年人精品一区二区| 成人精品一区二区免费|