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

    Joint association analysis method to dissect complex genetic architecture of multiple genetically related traits

    2020-10-21 10:01:26FengLinGuonQiTingXuXingyngLouYongoHongHimingXu
    The Crop Journal 2020年5期

    Feng Lin,Guon Qi,Ting Xu,Xingyng Lou,Yongo Hong*,Himing Xu,*

    aInstitute of Crop Science,College of Agriculture and Biotechnology,Zhejiang University,Hangzhou 310058,Zhejiang,China

    bDepartment of Mathematics,Zhejiang University,Hangzhou 310058,Zhejiang,China

    cDepartment of Biostatistics,Colleges of PHHP and the College of Medicine,University of Florida,Gainesville,FL 32610,USA

    dState Key Laboratory of Rice Biology and Zhejiang Key Laboratory of Super Rice Research,China National Rice Research Institute,Hangzhou 311401,Zhejiang,China

    ABSTRACT

    1.Introduction

    Most complex syndromes,such as asthma,metabolic syndrome,coronary heart disease,type 2 diabetes,consists of a large number of highly related,rather than independent,clinical or molecular phenotypes[1].The prevalent viewpoint suggests that phenotypic correlation is ascribed to 3 resources:pleiotropy,linkage disequilibrium and shared environment.When one or several same genetic variants affect genetic variation of multiple related traits,the phenotypes of an individual will tend to be correlated,known as pleiotropy of gene effect.Similarly,linkage disequilibrium and shared environment effect also introduce some level of correlation between traits.Taking these correlations into account is essential to genetic dissection of complex traits.

    Presently,the GWAS has been a routine way to explore complicate genetic architecture of quantitative traits;however,most GWAS to date have been conducted on each individual phenotype based on the single-locus additive association model.Less attention has been devoted to the fact that some surveyed phenotypes may be genetically correlated.Multiple studies have employed phenotypic data reduction methods to conduct association mapping;however,reduction methods will result in a considerable loss of power in detecting genetic variants[2–4].Joint analysis of multiple traits would increase power because intuitively,multiple trait measurements increase sample size relative to a single trait and exploit the covariances of other genetically related phenotypes;As a result,it will improve the ability to capture weak association signals which are most likely missed in separate analysis of trait.Hence,it is natural and imperative to develop new association analysis method for identifying genetic variants involved in genetic variation of multiple traits.

    Given the importance of joint analysis of multiple traits,there is a long history of multi-trait models in quantitative genetics.MultiPhen[5]uses ordinal regression to regress SNP genotypes on a group of phenotypes,and applies one omnibus statistic to test whether all regression weights in the model are together significantly different from zero.Unlike other methods,TATES[6]uncovers not only the genetic variants shared with multiple phenotypes as well as phenotype specific variants.Multivariate linear mixed models(mvLMMs)are considered powerful tools for testing associations between SNPs and multiple correlated phenotypes.Korte et.al[7]extended the linear mixed-model approach to carry out GWAS of correlated phenotypes,deriving a fully parameterized multi-trait mixed model(MTMM)that considers both the within-trait and between-trait variance components simultaneously for multiple traits.Zhou and Stephens[8]proposed an efficient algorithm in the genomewide efficient mixed model association(GEMMA)software for fitting mvLMMs,which improve computing speed,power and P-value calibration compared with existing methods.GEEGMDR[9]is an effective multivariate method to detecting gene-gene(epistasis)and gene-environment(GE)interactions underlying multiple complex traits.

    Moreover,all above mentioned multivariate methods achieve maximal power only in special circumstances and are based on simple additive model.It has been well documented that genetic variation of a quantitative trait is mostly controlled by a set of genes,which may function individually or jointly with additive or non-additive genetic effects[10,11].Non-additivity refers to interaction between alleles at same locus(dominance)or at different loci(epistasis),which has long been recognized fundamentally important to understand the structure and function of genetic pathways and the evolutionary dynamics of complex genetic system.Moreover,in contrast to Mendelian traits controlled by several major genes,the phenotypic variation of quantitative traits,also known as complex traits,depends not only on major genes but also many minor-effect alleles whose effects are sensitive to genetic,sexual,parent and environmental exposure[12].Thus,unravelling the causality of genetic variation and understanding the biological mechanism requires an analysis that jointly considers the genegene and gene-environment interactions.

    In this study,we developed a joint association analysis method for mapping common or phenotype-specific genetic causal variants for multiple related traits.Our mapping model,incorporates epistasis and GE interactions,can better reflect underlying biological processes.The effectiveness of the method was simulated,and demonstrated by the analysis of four correlated blood lipid traits in Multi-Ethnic Study of Atherosclerosis(MESA)Cohort[13]and two yield-related traits in a rice immortalized F2dataset[14].Finally,a fast parallelization computing software package for the proposed method,QTXNetwork-MT,was developed by combining computation of GPU(graphics processing unit,a type of hardware implementation of parallel computation adapted in many scientific computations)with CPU(central processing unit)[15],which can be downloaded from http://ibi.zju.edu.cn/software.

    2.Methods and materials

    2.1.Statistical model

    Considering a natural population consists of n individuals and the phenotypic performance of m correlated quantitative traits are measured in p different environments,the genetic variation of the traits are jointly controlled by s quantitative trait SNPs(QTSs)in which t pairs of QTSs are involved in epistatic interactions.A mixed linear model(MLM)for multiple traits association mapping can be established by extending the standard MLM for single trait as follows:

    For all individuals,the model can be expressed in vectormatrix format as follows:

    2.2.Mapping strategy

    Since all QTSs associated with the surveyed traits are unknown,we first used one/two-dimensional genome-wide screening to detect the QTS with significant individual effects or epistatic effects,then all detected QTS were used to construct the QTS full model for estimation of genetic effects of all QTS.

    One Dimensional Scanning:we conduct one dimensional whole genome scanning to search the significant loci using the following model

    where,ailkand dilkare the additive and dominance effects of the tested SNP in the k-th environment,respectively,μikis the mean of the i-th trait in the k-th environment;all other parameters have the same definitions as those in the Eq.(1).

    We also can express the above model in matrix and vector for all phenotypes as follows:

    Two-dimensional Scanning:In order to search for all possible candidate epistasis,we do two-dimensional scanning.In the scanning procedure,we include the detected significant SNPs in the model as cofactors.The mapping strategy is the same as one dimensional scanning.If we detected s significant SNPs in one-dimensional scanning,we can use the following model to test the significance of epistasis for each pair of SNPs:

    where,aailhk,adilhk,dailhk,and ddilhkare four epistatic effects(additive-additive,additive-dominance,dominance-additive and dominance-dominance effects)of the paired testing SNPs in the k-th environment for the i-th trait,respectively;all other parameters have same definitions as those in the model(1)and(3).The model(5)can also be expressed in matrixvector format for all phenotypes as the Eq.(4).Twodimensional scanning is conducted and the permutation testing is employed to determine the empirical threshold value of the Lambda statistic(Λ).All significant paired epistatic SNPs are distinguished.In order to remove the false positive QTSs,stepwise regression on all distinguished QTSs is performed,as a result,a QTS full model could be established based on all significant SNPs.

    2.3.Genetic effect estimation

    After obtaining the QTSs,the QTS full model could be constructed and then used to estimate all the parameters in the model(1).All the parameters in the model(1),as well as their standard errors and p-values of significance,were estimated based on the samples generated by the method of Markov chain Monte Carlo(MCMC)employing a conservative burn-in period of 20,000 cycles,a chain length of 200,000 and a thinning interval of 10 cycles and same prior distributions described by Yang et al.[16].We first generated a set of initial estimates of parameters including variances of random effects by the minimum norm quadric unbiased estimation(MINQUE)[17],the fixed effects by the ordinary least square estimation(OLSE)and the random effects by the adjusted unbiased prediction(AUP).And then these estimates are used as the prior values of parameters in the MCMC procedure for the mixed linear model implemented with Gibbs sampling.Parameter estimation and statistical inference are conducted based on the Gibbs samples[16,18,19].

    2.4.Analysis of human data and rice data

    2.4.1.Human data

    A set of real data in Multi-Ethnic Study of Atherosclerosis(MESA)was downloaded from the dbGaP of NCBI and exemplified the method.MESA is a prospective populationbased study focusing on investigating the prevalence,correlations,and progression of subclinical cardiovascular diseases(CVD),which consists of 6814 men and women aged 45–84 years[13]from four ethnic groups,38% of the cohort are European-American,28% African-American,22%Hispanic-American,and 12% Asian,predominantly of Chinese descent.The phenotype and genotype data are available on the database of Genotypes and Phenotypes(dbGaP)through dbGaP accession number phs000209.v13.p3.

    Single nucleotide polymorphisms(SNPs)were filtered according to the following criteria implemented by software PLINK[20]:(1)SNP call rate>95%;(2)subject call rate>95%;(3)Hardy Weinberg equilibrium P-value>1e-6;(4)polymorphic at least in one ethnic group(i.e.,no monomorphic SNP).After filtered,5018 individuals with 20,697 SNPs on the candidate genes of cardiovascular heart diseases were used to perform multi-trait association analysis.

    2.4.2.Rice data

    A rice immortalized F2(IF2)public dataset was downloaded from the website https://figshare.com/s/0773080c122d11e58b6306ec4bbcf141,which consists of 1619 markers and 278 samples[14].The traits were measured in two consecutive years.Here,the yield(Yield)and grain number per plant(Grain)traits were analyzed to illustrate our method.

    3.Results

    3.1.Simulation scenarios

    To investigate the efficiency and accuracy of the proposed method,we performed a Monte Carlo simulation study on the unbiasedness and robustness of our model.One natural population consisting of 500 subjects was simulated in two environments.100 equi-frequent biallelic SNPs were evenly distributed on 5 chromosomes,7 of which were randomly selected as simulated genetic variants(hereafter named by QTS)for two traits,wherein,5 QTSs were shared with both traits and 2 is trait-specific.The correlated phenotypes of two traits were generated by adding genetic effects to the 7 QTS.The heritability of each trait which could be accounted for the QTSs was set to 15%and 20%,respectively.The true values of the QTS position and genetic effects were presented in Tables 1 and 2.Allowing for variation in the phenotypic correlation,we also set the residual correlation to 0.40 between two traits.In order to compare the joint analysis and separate analysis of two traits,we also performed the single trait association analysis with the same simulated data.Each simulation was repeated 1000 times.The power(%)for a QTS was calculated as the percentage of this true locus correctly detected in the total number of simulations.

    3.2.Type one error rate and statistical power of detecting QTS

    Simulation results are summarized in Tables 1,2,and Figs.1–3.Fig.1 presents the Quantile-quantile plot of significance level versus Type I error rate drawn from the 1000 simulations in assumption of no gene-gene interaction and no residual correlation.All empirical Type I error rates were close to the nominal level.The simulation results indicate that this new method had the correct Type I error rates for multiple traits,supporting the validity of the new approach under the null model.Fig.2 presents the power of each QTS and paired epistatic QTSs detected by the joint analysis method and/or the separate analysis method,respectively.In general,the new approach had a marked advantage in improving power over the separate analyses.Fig.3 shows that the pleiotropic QTSs(Q3,Q4,Q5,and EQ1)could be detected with higher power at any significant level by the joint analysis than by the separate analyses of individual traits;especially for Q5,which contributed opposite genetic component effects to two traits,the power was substantially increased when the joint analysis was employed(Table 1 and Fig.3).In contrast with Q5,the pleiotropic Q4 had the same direction of component effects in additive,dominance and their interaction with environment on the two traits,and its power was also increased to over 95%in the joint analysis,but only 62.6%and 81.4%in two separate analyses.For Q3 and EQ1(Q6–Q7),although they had the component effects on the trait 1 and 2 in the same and the different directions,respectively(Table 1),they could be detected consistently with higher power by the joint analysis;the power were 99.9% at the significance level of 0.01 for Q3 and 94.5% at the significance level of 1×10-10for EQ1 in the joint analysis,but the counterparts were only 75.1%and 92.1%for Q3 and 74% and 85.8% for EQ1 in the individual analyses,respectively(Fig.3).On the other hand,for the trait-specific QTSs(Q1,Q2),although the separate analyses performed more efficiently than the joint analysis since the relatively simple model was employed in the analysis of individual trait,the loss of power in detecting trait-specific QTSs was minor when the joint analysis was applied and totally acceptable if the multiple testing issue in individual analyses was considered(Fig.3).These results suggested the proposed new method based on the joint analysis of multiple traits could effectively improve the detection of pleiotropic genes but might have an ignorable loss in power when the traits surveyed are genetically independent.

    Table 1–Mapping powers and estim ates of QTS m ain effects and their interaction effects with environments for two simulated traits(T1 and T2).

    Table 2–Mapping powers and estimates of epistasis QTS(EQ1)main effects and their interaction effects with environments for two simulated traits(T1 and T2).

    3.3.Estimation of QTS effects

    Tables 1 and 2 showed that both types of QTS,main effectonly QTSs(Q1,Q2,Q3,Q4,Q5)and epistatic QTSs EQ1(Q6-Q7)had bias(the difference between the estimated and the true value of parameter)of less than 5%of the true value for each of genetic main effects(A,D,AA,AD,DA,and DD),indicating that the genetic main effects could be unbiasedly estimated by our method.As for the components effects of gene by environment interaction,a small but acceptable bias would occasionally occur for some QTSs,such as Q2 with estimates of 0.56/-0.56 for the true DE effect of 0.63/-0.63,Q3 with estimates of 0.85/-0.86 for true AE of 0.94/-0.94,and Q4 with estimates of 0.51/-0.51 for true AE of 0.57/-0.57(Table 1),respectively;compared with the additive/dominance by environment interaction,relatively lager bias could be observed on the estimates of epistasis-environment interaction effects,such as the AAE,ADE and DAE of EQ1(Table 2).In general,the QTS effects could be unbiasedly or well estimated by our mixed linear model-based method.

    Fig.1–Quantile-quantile plot of significance level and type I error rate for joint and separate analysis.The type I error rate for main effect and epistatic effect is evaluated by Joint(T1,T2)and Separate(T1)in assumption of no gene-gene interaction and no residual correlation.Joint(T1,T2),the joint analysis of trait 1 and trait 2;Separate(T1),the separate analysis of trait.

    Fig.2–Statistical power at the significance level of 0.05.Separate(T1)and Separate(T2),the separate analysis of trait 1 and trait 2,respectively;Joint(T1,T2),the joint analysis of trait 1 and trait 2;Q1 and Q2,trait-specific QTSs,Q1 for trait 1,Q2 for trait 2;Q3,Q4,and Q5,pleiotropic main effect-only QTSs for trait 1 and trait 2;and EQ1 pleiotropic epistatic QTSs(Q6 and Q7)for trait 1 and trait 2.

    Fig.3–Statistical powers for all the simulated main effect-only QTSs and epistatic QTSs at various significance levels.Separate(T1)and Separate(T2),the separate analysis of trait 1 and trait 2,respectively;Joint(T1,T2),the joint analysis of trait 1 and trait 2;Q1 and Q2,trait-specific QTS,Q1 for trait 1,Q2 for trait 2;Q3,Q4,and Q5,pleiotropic main effect-only QTSs for trait 1 and trait 2;EQ1 pleiotropic epistatic QTSs(Q6 and Q7)for trait 1 and trait 2.

    3.4.Application to two genetic studies

    3.4.1.Human data

    To validate the effectiveness of the joint association mapping of multiple traits in dissection of complicated genetic architecture of complex traits,the new method was used to analyze four blood metabolites-related traits which are crucial clinic phenotypes of human cardiovascular heart disease[21],namely triglyceride(TG),low density protein(LDL),high density protein(HDL),and total cholesterol(TC).The genotype and phenotype data were collected from the Multi-Ethnic Study of Atherosclerosis(MESA)Cohort(see Methods for details)and downloaded from the dbGaP of NCBI.Covariance analysis was performed based on a mixed linear model with race as covariate,the results clearly showed there were significant genetic correlations between the four traits(Table 3).TC had a strong positive genetic correlation with other three traits,and the genetic correlation with HDL was the highest with a correlation coefficient of 0.90.HDL had a significant middle negative genetic correlation with TG(-0.40)and a significant but relatively small negative correlation with LDL(-0.05).

    We applied our methods to conduct the joint and the separate association analyses of four traits,respectively,the results are shown in Tables 4 and 5.Four main effect-only QTSs(rs3740767,rs10750097,rs1800588,and rs4775041)and one pair of epistatic QTSs(rs3799481–rs13254175)were detected to have significant association with the traits only by our new method but not by the separate analysis of individual traits(Table 4 and Fig.4A).Interestingly,rs3740767 is a variant which has been reported in the previous study on peripheral arterial disease[22];rs596985,rs2494283,rs301,rs5072,rs12679834,rs12287066,and rs3799481 were detected only by the single trait analysis,wherein,the first five QTSs were all associated with TG and the last two QTSs associated with HDL;rs12740374,rs13702,rs2266788,rs3764261 and rs7499892 were detected by both the joint or the separate analysis of four traits;generally,the joint analysis could identify more significant loci than the separate analysis,for example,the separate analysis method detected only 5,1,1,and 6 QTSs for HDL,LDL,TC,and TG,respectively,but accordingly,8,5,6,and 7 QTSs by the joint analysis(Fig.4C–F).Compared with the single-trait association analysis in ability to explore a complicated genetic architecture of complex traits,the multi-traits association analysis could provide more insights into gene pleiotropy and molecular mechanism of genetic dependency of related traits;forexample,the joint analysis revealed that rs3740767,rs3764261 and the paired epistatic QTSs rs3799481–rs13254175 are pleiotropic variants,they contributed genetic effects to all the four phenotypes of human peripheral arterial disease(Tables 4 and 5;Fig.4B),of which rs3764261 mainly took additive genetic effect on four traits,while rs3740767 mainly affected the traits via additive,dominance and additive-race interaction effects,and rs3799481–rs13254175 had epistatic effects and their interaction with races,but no genetic effects of individual locus.

    Table 3–Genetic correlations between four CVD_related traits.

    To understand regulation patterns and intensity of the detected significant QTSs for the four traits,we also estimated the main effects,gene-by-race interaction effects as well as epistatic effects(Tables 4 and 5).Overall,16 QTSs mostly had additive effects and interaction effects with races,and one pair of epistatic QTSs were detected by the joint analysis and separate analysis.It is noteworthy that rs12740374 on chromosome 1 was detected both in the joint and the separate analyses of TC and LDL,which compared with the other detected QTSs accounted for the largest phenotypic variation for TC and LDL,2.50% and 4.15% in the joint analysis and 3.11% and 5.38% in the separate analysis,respectively.Some QTSs(rs12740374,rs3740767)exhibited additive,dominance effects combined with additive,dominance by races interaction effects,while some QTSs(rs13702)had genetic main effects regardless of races,which indicates the QTSs are not sensitive to races.For examples,rs13702 on chromosome 8,for TG and HDL,only contributed additive effects of 10.73 and-1.73,respectively.In addition,some QTSs only took effects through interaction with race,which implies that theirexpression is affected by race,like rs7499892 on chromosome 16 for LDL,TC in the joint analysis expressed only in Hispanic with race-specific dominance effects of-2.75 and-3.06 respectively.Hence,we should pay attention on the difference of gene expression across races.Moreover,the joint analysis also detected one pair of epistatic QTSs rs3799481–rs13254175 on chromosomes 6 and 8,respectively,accounting for a relatively larger proportion(3.44% and 4.10%)of phenotypic variation for TG and HDL,respectively;notably,the epistasis took opposite effects on TG and HDL in the detected genetic components.rs13702 and rs2266788 had additive effects of 10.73,-19.49 on TG and of-1.73 and 1.00 on HDL respectively,which was in accordance with their negative correlation of-0.40(Table 3).

    Table 4–The detected significant QTSs and their effects by joint and separate analysis of four coronary heart disease related traits.

    Table 5–The detected paired QTSs with significant epistasis and the estimated effects by joint analyses of four coronary heart disease related traits.

    Fig.4–Venn diagram for the detected significant QTSs from joint analysis and separate analysis of four CVD traits.TG,triglyceride;TC,total cholesterol;LDL,low density lipoprotein;HDL,high density lipoprotein;TC+TG+LDL+HDL(Joint),TC(Joint),TG(Joint),HDL(Joint),LDL(Joint)indicate joint analysis of four traits;TC,TG,LDL,HDL indicate separate analysis of each trait;the number in non-overlapping area is the number of QTSs which could only be detected by one analysis method,in contrast,the number in overlapping area is the number of QTSs which could be detected by multiple analysis methods.

    Almost all of the detected QTSs are located in the annotated genes related with CVD traits,which demonstrates the effectiveness of our method.Totally,16 QTSs were identified in the regions of the nine annotated genes,which are CELSR2,ROR1,LPL,APOA1,APOA5,LIPC,CETP,PEX7,MSRA,respectively.QTS rs12740374 on chromosome 1 is located in noncoding region of gene CELSR2,as the causal gene at the 1p13 locus for LDL[23].rs301,rs13702,rs12679834 are located in the lipoprotein lipase(LPL)locus on chromosome 8 which encodes a crucial enzyme in lipoprotein metabolism whose function partially determines concentrations of circulating TG and HDL.Moreover,the LPL 30 UTR rs13702 has been identified to induce an allele-specific interaction with miR-410,a regulator of LPL in human[24].APOA1 encodes apolipoprotein A-I,which is the major protein component of high density lipoprotein(HDL)in plasma[25].The gene APOA5 encodes an apolipoprotein,a component of HDL that plays an important role in regulating TG,a major risk factor for coronary artery disease[26].The gene LIPC was not detected in the separate analysis,but showed highly significant association in joint analysis.The LIPC located on chromosome 15q22,has both triglyceride lipase and phospholipase activities and is involved in different processes of lipoprotein metabolism[27].CETP is a hydro-phobic glycoprotein which plays a major role in RCT(Reverse Cholesterol Transport)from tissues to the liver.By enabling transfer of cholesteryl esters from HDL to LDL,CETP enables remodeling of plasma lipoprotein[28].Generally,almost all of the identified genes are directly or indirectly involved in regulating the CVD traits.These findings provide greater insights into the genetic regulation and gene interactions in metabolic pathways of CVD traits.

    3.4.2.Rice data

    We applied our method to conduct the joint and the separate association analyses of Yield and Grain traits,respectively,the results were summarized in Tables S1 and S2.Two QTSs(BIN3,BIN1098)and one pair of epistatic QTSs(BIN189–BIN278)were detected to take significant genetic main effects on the traits only by the joint analysis not by the separate analysis of two traits.Almost all QTSs and epistatic QTSs of Yield or Grain which were detected by the separate analysis could also be detected by the joint analysis.Furthermore,the joint analysis revealed that BIN50,BIN436,BIN875,BIN1004,BIN1057,BIN1098 and the paired epistatic QTSs BIN189–BIN278 are pleiotropic variants,of which BIN875,BIN1004,BIN1057 were also found to be pleiotropic by the separate analysis.

    In terms of genetic effects,the detected QTSs except BIN1004 and BIN1057 took only genetic main effects in additive or dominance on two traits.In total,only two pairs of epistatic QTS(BIN189–BIN278,Bin594–Bin1468)were detected and contribute only epistasis main effects to Yield or Grain,wherein,the pair of BIN189 and BIN278 effect both on Yield and Grain,but the pair of Bin594 and Bin1468 only on Grain.Two QTSs(BIN1004,BIN1057)exhibited additive-year interaction effects,which indicates their expressions are sensitive to the environment,for example,BIN1004 could increase additive effect of 1.21 in the first environment but reduce 1.22 in the second year on the Yield according to the joint analyses.The pleiotropic epistatic QTSs BIN189–BIN278 accounted for a relatively larger proportion(19.05% and 10.05%)of phenotypic variation for Yield and Grain,respectively.

    4.Discussion

    Most human diseases such as diabetes,hypertension and obesity,are usually measured by several clinic traits between which there are some certain genetic correlation;on the other hand,same human disease symptom may arise from several pathogenesis etiology connected with different causal phenotypes.Further,ever-growing evidence has pointed to the view that gene by environment interactions are also important genetic buffer of phenotypic variations[29–31].Evidently,it is very imperative to take into account of these complicate genetic features in GWAS of multiple traits.This study proposed a new genome-wide association analysis method for joint analysis of multiple quantitative traits,which is built on a multivariate mixed linear model with inclusion of each genetic component effects of all QTSs,paired interacting QTSs and their interaction with environment,and the covariance information between traits is exploited in testing effect of QTS.Extensive simulation showed that the new method outperforms the separate analysis of individual traits not only in identifying more minor-effect genes as well as in revealing pleiotropic genes and molecular mechanism of complicate genetic correlation of traits.In addition,the new method can provide unbiased estimates of almost all kinds of genetic component effects of single gene and paired epistatic genes,although the simulation study observed a little bit over or lower estimates for very fewer interaction components effects of single or paired epistatic genes with environments.Thus,the proposed multi-traits GWAS method provides an effective solution to dissect complicate genetic architecture of a group of genetically correlated phenotypes.

    The application of the proposed approach to the MEAS data revealed that the seven detected pleiotropic genes CELSR2,LPL,APOA5,LIPC,CETP,PEX7,MSRA take different effects on the four human blood metabolites-related traits,most of which have been reported in the previous GWAS studies involved in human blood metabolites-related traits except for PEX7 and MSRA.PEX7 has been reported involved in many metabolic processes such as fatty acid oxidation,biosynthesis of ether lipids and free radical detoxification,and the function of MSRA is to repair oxidative damage to proteins by enzymatic reduction of methionine sulfoxide,associated with waist circumference.Hence,the paired epistatic genes PEX7 and MSRA may be involved in a same fatty acid metabolic pathway to regulate the four human blood metabolites-related traits,thus controlling the obesity[33].The QTS rs3740767,a pleiotropy locus of all the four traits,was only reported to be associated with peripheral arterial disease which is caused by the disorders of the blood lipid traits.Thus,it is biologically plausible for variant rs3740767 to be involved in regulating the blood lipid traits and thereby involved in peripheral arterial disease.Compared with the separate analysis of each trait,the multi-traits analysis detected more minor heritability QTSs,like rs1800588,the heritabilities for TG,HDL,and TC were respectively 0.55,0.89,and 0.17,which may make a contribution to unveiling the problem of“missing heritability”.

    Our method provides a new approach of association analysis for dissecting complex genetic architecture of a group of correlated traits.However,what is the best strategy to analyze a set of data with multiple traits depends on the molecular mechanism of genetic correlation of surveyed traits.If a variant is of trait-specific,then the separate analysis on the trait will generally be more powerful to detect.On the other hand,if a variant affects multiple traits,then the joint analysis of all traits will be more powerful.Hence,the joint and the separate association analysis of multiple traits are not exclusive but complementary to each other in GWAS.In addition,it should be noted that the Wilks’Lambda statistic in our method is designed to test the significance of effects on all surveyed traits for one individual or pair of QTSs,not on one specific trait,only in the step of parameter estimation of QTS via Gibbs sampling based on the final QTS full model,each component effect of QTS on each trait will be tested by t-statistic and then the trait-specific or pleiotropic effects of QTS is inferred.Different from the proposed method based on multivariate mixed linear model,the univariate conditional model of one trait on another or several other correlated traits can also be used to dissect genetic architecture of a group of related traits.If a QTS could be detected and similar strength of effect was estimated by both unconditional model and conditional model of one trait,we could conclude that the QTS is trait-specific and only take effect on the trait;on the contrary,if a QTS could be detected by unconditional model but not by conditional model,then no trait-specific effect of the QTS on the trait could be concluded;third,the power of pleiotropic QTS will be reduced when conditional model was applied since partial effects of the QTS has been removed by the conditional trait in the model.

    Onto implement the proposed method,we developed a software in C++programming language which is able to analyze multiple correlated traits,which can detect dominance,epistasis and GE interactions,and estimate the corresponding variance component and genetic effects.In addition,to deal with tremendous matrix computation and speed up the analysis of real data,especially for detecting epistasis for whole genome sequence data,a hybrid parallel computing technology of CPU and GPU was employed in the software QTXNetwork-MT,which can improve the computation efficiency substantially,so that the new method could be applied in association analysis of big data.Finally,the QTXNetwork-MT has been available at the website http://ibi.zju.edu.cn/software.

    5.Conclusions

    We proposed a multi-trait association analysis method based on multivariate mixed linear model including additive,dominance,epistasis effects of QTS and their interaction effects with environments.Wilks statistic is employed to detect the QTSs with significant individual or di-genic interaction effects.Based on all detected QTSs,a QTS full model is built and used to estimate the genetic component effects of QTSs.Simulation study and real data analysis illustrated that the new method could achieve higher statistical power in detecting middle-to minor-effect QTS compared with the conventional separate analysis of each trait,and be more effective to identify pleiotropic genes in correlated phenotypes,which may be involved simultaneously in several related biological pathways.

    Declaration of competing interest

    The authors declare that there are no conflicts of interest.

    Acknowledgments

    The study was supported by the National Key Research and Development Program of China(2016YFC1303300),National Natural Science Foundation of China(31671570,31871707),National Science Foundation(DMS2002865)and the 111 Project(BP2018021).The funders had no role in study design and data analysis.The authors thank the investigators,the staff,and the participants of the MESA(The Multi-Ethnic Study of Atherosclerosis)for their valuable contributions.A full list of participating MESA investigators and institutions can be found at http://www.mesa-nhlbi.org.

    CRediT authorship contribution statement

    Feng Lin conducted simulation,program development and wrote the manuscript.Guoan Qi checked the program codes.Ting Xu took part in the design of the software QTXNetowk-MT.Xiangyang Lou checked the method and improved the manuscript.Yongbo Hong analyzed the example data and took part in the writing of the manuscript.Haiming Xu developed the model and analysis method,and improved the manuscript.

    Appendix A. Supplementary data

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

    久久亚洲真实| 亚洲成人久久爱视频| 在线观看午夜福利视频| 91久久精品国产一区二区成人 | 免费看十八禁软件| 久久精品影院6| 窝窝影院91人妻| 国产色爽女视频免费观看| 最后的刺客免费高清国语| 亚洲欧美日韩卡通动漫| 久久久色成人| 国产中年淑女户外野战色| 欧美成人免费av一区二区三区| 国产精品永久免费网站| 国产aⅴ精品一区二区三区波| 亚洲av成人精品一区久久| 国产成人a区在线观看| 日韩国内少妇激情av| 高潮久久久久久久久久久不卡| 免费看日本二区| 色视频www国产| 国产v大片淫在线免费观看| 别揉我奶头~嗯~啊~动态视频| 99在线视频只有这里精品首页| 日韩亚洲欧美综合| 给我免费播放毛片高清在线观看| 黄片大片在线免费观看| 中亚洲国语对白在线视频| 免费搜索国产男女视频| 日本 欧美在线| 十八禁人妻一区二区| 国产精品久久电影中文字幕| 在线观看66精品国产| 午夜福利欧美成人| 中亚洲国语对白在线视频| 国产成人啪精品午夜网站| 亚洲一区二区三区不卡视频| 夜夜看夜夜爽夜夜摸| 老熟妇乱子伦视频在线观看| 两性午夜刺激爽爽歪歪视频在线观看| av天堂中文字幕网| 动漫黄色视频在线观看| 久久国产精品影院| 日韩中文字幕欧美一区二区| 久久久久国产精品人妻aⅴ院| 色综合婷婷激情| 两性午夜刺激爽爽歪歪视频在线观看| 国产激情欧美一区二区| 脱女人内裤的视频| a级一级毛片免费在线观看| 亚洲中文日韩欧美视频| 国内精品一区二区在线观看| 两个人的视频大全免费| 丰满乱子伦码专区| 美女 人体艺术 gogo| 国产激情偷乱视频一区二区| 人人妻人人澡欧美一区二区| 乱人视频在线观看| 亚洲欧美激情综合另类| 日韩国内少妇激情av| 18禁黄网站禁片免费观看直播| a级毛片a级免费在线| 国产美女午夜福利| 中文资源天堂在线| 国产真人三级小视频在线观看| 51午夜福利影视在线观看| av片东京热男人的天堂| www日本在线高清视频| 91av网一区二区| 欧美日韩乱码在线| svipshipincom国产片| 国语自产精品视频在线第100页| 色噜噜av男人的天堂激情| 69av精品久久久久久| 精品日产1卡2卡| 国内精品一区二区在线观看| 日韩欧美三级三区| 男插女下体视频免费在线播放| 国产在线精品亚洲第一网站| 精品电影一区二区在线| 18禁在线播放成人免费| 色视频www国产| 搞女人的毛片| 九九久久精品国产亚洲av麻豆| 亚洲精华国产精华精| 日本在线视频免费播放| 999久久久精品免费观看国产| 日本精品一区二区三区蜜桃| 狂野欧美白嫩少妇大欣赏| 一级毛片高清免费大全| 婷婷精品国产亚洲av在线| 亚洲性夜色夜夜综合| av片东京热男人的天堂| 亚洲精品一卡2卡三卡4卡5卡| 少妇人妻一区二区三区视频| 九九久久精品国产亚洲av麻豆| 美女黄网站色视频| 国产精品一区二区三区四区久久| 99精品欧美一区二区三区四区| 天堂√8在线中文| 丰满乱子伦码专区| 日韩av在线大香蕉| 99精品久久久久人妻精品| 欧美丝袜亚洲另类 | 激情在线观看视频在线高清| 我的老师免费观看完整版| 一级a爱片免费观看的视频| 亚洲精品粉嫩美女一区| 日韩有码中文字幕| 国产伦一二天堂av在线观看| 一a级毛片在线观看| 欧美zozozo另类| 亚洲国产精品久久男人天堂| 午夜福利高清视频| 叶爱在线成人免费视频播放| svipshipincom国产片| 亚洲 欧美 日韩 在线 免费| 亚洲av电影在线进入| 国产97色在线日韩免费| 网址你懂的国产日韩在线| 久久久成人免费电影| 美女cb高潮喷水在线观看| 脱女人内裤的视频| 美女免费视频网站| 成年人黄色毛片网站| 国产色婷婷99| 精品久久久久久,| 两个人视频免费观看高清| 尤物成人国产欧美一区二区三区| 波多野结衣高清无吗| 国产黄色小视频在线观看| 成人亚洲精品av一区二区| 日本成人三级电影网站| 国产精品三级大全| 欧美中文日本在线观看视频| 中文亚洲av片在线观看爽| 天堂动漫精品| 午夜激情福利司机影院| 老司机在亚洲福利影院| 女生性感内裤真人,穿戴方法视频| 免费在线观看亚洲国产| 亚洲av二区三区四区| 久久中文看片网| 午夜精品久久久久久毛片777| 18+在线观看网站| av国产免费在线观看| 色精品久久人妻99蜜桃| 亚洲在线自拍视频| 成人三级黄色视频| 毛片女人毛片| 日韩欧美在线乱码| 久久久成人免费电影| 精品一区二区三区视频在线 | 亚洲七黄色美女视频| 99久久成人亚洲精品观看| 国产一级毛片七仙女欲春2| 国产一区二区亚洲精品在线观看| 校园春色视频在线观看| 悠悠久久av| 精品久久久久久久毛片微露脸| 两人在一起打扑克的视频| 国产精品国产高清国产av| www.999成人在线观看| 亚洲专区国产一区二区| 国产精品99久久99久久久不卡| 免费观看精品视频网站| 久久精品亚洲精品国产色婷小说| 国产乱人伦免费视频| 婷婷精品国产亚洲av在线| 高清日韩中文字幕在线| www.www免费av| 精品无人区乱码1区二区| 国产精品 欧美亚洲| 啦啦啦韩国在线观看视频| 国产老妇女一区| 久久精品国产自在天天线| 欧美激情在线99| 国产精品久久电影中文字幕| 变态另类成人亚洲欧美熟女| 一a级毛片在线观看| 观看美女的网站| 久久草成人影院| 色老头精品视频在线观看| 黄色女人牲交| 成年人黄色毛片网站| 亚洲中文日韩欧美视频| 欧美av亚洲av综合av国产av| 欧美性猛交黑人性爽| bbb黄色大片| 国产熟女xx| 色噜噜av男人的天堂激情| 免费观看人在逋| 国产淫片久久久久久久久 | 精品久久久久久久毛片微露脸| 亚洲成av人片免费观看| 啦啦啦观看免费观看视频高清| 午夜福利高清视频| 久久久久久国产a免费观看| 别揉我奶头~嗯~啊~动态视频| 中文字幕高清在线视频| 免费看a级黄色片| 岛国在线免费视频观看| 久久久久精品国产欧美久久久| 国产亚洲精品一区二区www| 色哟哟哟哟哟哟| 国产极品精品免费视频能看的| 亚洲美女视频黄频| 亚洲精品成人久久久久久| 久久中文看片网| 99在线视频只有这里精品首页| 小蜜桃在线观看免费完整版高清| 悠悠久久av| 久久久成人免费电影| 久久久久国内视频| 国产精品永久免费网站| 岛国在线观看网站| 人人妻,人人澡人人爽秒播| 女人十人毛片免费观看3o分钟| 国产 一区 欧美 日韩| 亚洲精品亚洲一区二区| 18禁美女被吸乳视频| 又爽又黄无遮挡网站| 精品久久久久久久人妻蜜臀av| 免费看a级黄色片| 嫩草影院入口| 欧美大码av| 欧美+日韩+精品| 女人高潮潮喷娇喘18禁视频| 69av精品久久久久久| 亚洲在线观看片| 亚洲中文字幕日韩| 国产成人系列免费观看| 在线观看美女被高潮喷水网站 | 国产成人影院久久av| 丰满人妻熟妇乱又伦精品不卡| 在线观看66精品国产| 亚洲精品美女久久久久99蜜臀| 国产日本99.免费观看| 亚洲va日本ⅴa欧美va伊人久久| 久久久久久九九精品二区国产| 舔av片在线| 亚洲国产精品久久男人天堂| 国产伦一二天堂av在线观看| 国产精品一区二区三区四区免费观看 | 老司机在亚洲福利影院| 一a级毛片在线观看| 午夜免费激情av| 看免费av毛片| 婷婷亚洲欧美| 国产欧美日韩精品亚洲av| 嫁个100分男人电影在线观看| 天堂影院成人在线观看| av专区在线播放| 手机成人av网站| 九九热线精品视视频播放| 中文在线观看免费www的网站| 无遮挡黄片免费观看| 麻豆成人午夜福利视频| 国产极品精品免费视频能看的| 最近最新中文字幕大全免费视频| 香蕉久久夜色| 99热只有精品国产| 欧美国产日韩亚洲一区| 啦啦啦免费观看视频1| 黄色片一级片一级黄色片| 变态另类成人亚洲欧美熟女| 免费人成视频x8x8入口观看| 成人永久免费在线观看视频| 19禁男女啪啪无遮挡网站| 国产精品永久免费网站| 男女下面进入的视频免费午夜| 亚洲人成伊人成综合网2020| 精品电影一区二区在线| av女优亚洲男人天堂| 国产精品98久久久久久宅男小说| 亚洲欧美日韩卡通动漫| 91在线观看av| 综合色av麻豆| 亚洲在线自拍视频| 97超视频在线观看视频| 亚洲人成网站在线播| 午夜福利在线在线| 亚洲成人久久爱视频| 全区人妻精品视频| 欧美在线一区亚洲| 老司机深夜福利视频在线观看| 少妇的逼好多水| 成人18禁在线播放| 一a级毛片在线观看| 两人在一起打扑克的视频| 人妻夜夜爽99麻豆av| 天天躁日日操中文字幕| 亚洲av五月六月丁香网| 丝袜美腿在线中文| 国产毛片a区久久久久| 99国产综合亚洲精品| 欧美日韩中文字幕国产精品一区二区三区| 国产综合懂色| 99久久精品热视频| 老司机福利观看| 神马国产精品三级电影在线观看| 一级作爱视频免费观看| 国产一级毛片七仙女欲春2| 两个人视频免费观看高清| 91九色精品人成在线观看| 亚洲精品粉嫩美女一区| 可以在线观看毛片的网站| 色吧在线观看| av黄色大香蕉| 高清日韩中文字幕在线| 亚洲第一电影网av| 精品一区二区三区人妻视频| 久久九九热精品免费| 黄色成人免费大全| 欧美xxxx黑人xx丫x性爽| 久久久久国产精品人妻aⅴ院| 免费看美女性在线毛片视频| 国产精品亚洲美女久久久| 欧美日韩中文字幕国产精品一区二区三区| 国产伦精品一区二区三区四那| 男女床上黄色一级片免费看| 精品一区二区三区av网在线观看| 88av欧美| 亚洲欧美日韩高清在线视频| 长腿黑丝高跟| 日韩大尺度精品在线看网址| 丰满人妻一区二区三区视频av | 午夜福利在线在线| 免费在线观看亚洲国产| av中文乱码字幕在线| 国产麻豆成人av免费视频| 1000部很黄的大片| 成人性生交大片免费视频hd| 麻豆久久精品国产亚洲av| 嫁个100分男人电影在线观看| 国产精品1区2区在线观看.| av视频在线观看入口| 欧美+日韩+精品| 国产精品精品国产色婷婷| 国产伦精品一区二区三区四那| 十八禁网站免费在线| 一本精品99久久精品77| 高清日韩中文字幕在线| 两个人看的免费小视频| 美女黄网站色视频| 夜夜躁狠狠躁天天躁| 97超视频在线观看视频| 久久久精品大字幕| 少妇的逼好多水| 怎么达到女性高潮| АⅤ资源中文在线天堂| 脱女人内裤的视频| 亚洲av日韩精品久久久久久密| 人妻丰满熟妇av一区二区三区| 国产91精品成人一区二区三区| 亚洲国产中文字幕在线视频| 国产精品香港三级国产av潘金莲| 亚洲国产精品久久男人天堂| 一进一出好大好爽视频| 中文字幕人妻熟人妻熟丝袜美 | 可以在线观看毛片的网站| 亚洲精品美女久久久久99蜜臀| 男女那种视频在线观看| 久久中文看片网| 午夜久久久久精精品| 欧美午夜高清在线| 最新美女视频免费是黄的| 欧美午夜高清在线| 精品一区二区三区av网在线观看| 国产av不卡久久| 国产成人av激情在线播放| 99精品欧美一区二区三区四区| 国产精品一区二区免费欧美| 中文字幕人妻丝袜一区二区| 免费观看精品视频网站| 久久久久久久午夜电影| 中文资源天堂在线| 亚洲欧美一区二区三区黑人| 欧美激情久久久久久爽电影| 免费人成视频x8x8入口观看| 女同久久另类99精品国产91| tocl精华| 亚洲精品亚洲一区二区| 欧美丝袜亚洲另类 | 国产99白浆流出| 久久午夜亚洲精品久久| 三级国产精品欧美在线观看| 18禁黄网站禁片午夜丰满| 伊人久久大香线蕉亚洲五| 欧美日韩福利视频一区二区| 欧美+亚洲+日韩+国产| 亚洲成人久久爱视频| 又粗又爽又猛毛片免费看| 可以在线观看的亚洲视频| 午夜福利高清视频| 桃色一区二区三区在线观看| 在线视频色国产色| 亚洲精品乱码久久久v下载方式 | 午夜视频国产福利| 叶爱在线成人免费视频播放| 岛国在线观看网站| 免费看十八禁软件| av在线蜜桃| 亚洲成人久久爱视频| 97碰自拍视频| 美女高潮的动态| 18禁黄网站禁片免费观看直播| 亚洲av电影不卡..在线观看| 精品国产美女av久久久久小说| 高清在线国产一区| www.999成人在线观看| 国产av麻豆久久久久久久| 九九热线精品视视频播放| 在线观看舔阴道视频| 免费观看人在逋| 成人高潮视频无遮挡免费网站| 欧美成人性av电影在线观看| 日本成人三级电影网站| 麻豆国产av国片精品| 中文字幕av成人在线电影| 美女 人体艺术 gogo| 亚洲乱码一区二区免费版| 亚洲欧美精品综合久久99| 少妇丰满av| 法律面前人人平等表现在哪些方面| 国产蜜桃级精品一区二区三区| 99久久成人亚洲精品观看| 亚洲,欧美精品.| 天堂网av新在线| 久久久久久国产a免费观看| 国产亚洲精品av在线| 搡老妇女老女人老熟妇| 在线观看66精品国产| 毛片女人毛片| 91字幕亚洲| 成人鲁丝片一二三区免费| 在线观看一区二区三区| 黑人欧美特级aaaaaa片| 九九热线精品视视频播放| 国产精品亚洲美女久久久| 国产黄色小视频在线观看| 天天一区二区日本电影三级| www.www免费av| 精品一区二区三区视频在线观看免费| 热99re8久久精品国产| 少妇人妻精品综合一区二区 | 91av网一区二区| 国产成人啪精品午夜网站| 国产精品免费一区二区三区在线| 怎么达到女性高潮| 听说在线观看完整版免费高清| 国产精品乱码一区二三区的特点| 日韩高清综合在线| 午夜福利高清视频| 久久草成人影院| 99久久精品热视频| 啦啦啦韩国在线观看视频| 舔av片在线| 欧美激情久久久久久爽电影| 久久性视频一级片| 久久精品人妻少妇| 成熟少妇高潮喷水视频| 亚洲av免费在线观看| 天堂av国产一区二区熟女人妻| 久久久精品欧美日韩精品| 在线十欧美十亚洲十日本专区| 欧美xxxx黑人xx丫x性爽| 悠悠久久av| 桃色一区二区三区在线观看| 搡女人真爽免费视频火全软件 | 国产一区二区在线观看日韩 | 久久精品国产自在天天线| 国产aⅴ精品一区二区三区波| 麻豆国产97在线/欧美| 国产成+人综合+亚洲专区| 亚洲国产精品成人综合色| 午夜两性在线视频| 99热6这里只有精品| 美女大奶头视频| 久久伊人香网站| 99精品欧美一区二区三区四区| 18+在线观看网站| 欧美在线一区亚洲| 12—13女人毛片做爰片一| 成人特级av手机在线观看| 少妇熟女aⅴ在线视频| 俺也久久电影网| 又粗又爽又猛毛片免费看| 99久久综合精品五月天人人| 美女高潮喷水抽搐中文字幕| 人人妻人人看人人澡| 51午夜福利影视在线观看| 国产亚洲精品久久久久久毛片| av女优亚洲男人天堂| 国产真实伦视频高清在线观看 | 国产伦精品一区二区三区视频9 | 国产私拍福利视频在线观看| 中文字幕高清在线视频| 夜夜爽天天搞| 精品人妻偷拍中文字幕| 亚洲欧美精品综合久久99| 97超视频在线观看视频| 在线观看舔阴道视频| 精品久久久久久久末码| 91字幕亚洲| 欧美日韩亚洲国产一区二区在线观看| 88av欧美| 色综合婷婷激情| 久久亚洲真实| 日韩人妻高清精品专区| 国内精品久久久久久久电影| 欧美+亚洲+日韩+国产| 日韩有码中文字幕| 国产三级中文精品| 在线观看免费午夜福利视频| 亚洲专区中文字幕在线| 国产精品久久电影中文字幕| 亚洲精品国产精品久久久不卡| 天天一区二区日本电影三级| 欧美成狂野欧美在线观看| 99国产极品粉嫩在线观看| 日本在线视频免费播放| 午夜福利免费观看在线| av福利片在线观看| 中文字幕精品亚洲无线码一区| 国内少妇人妻偷人精品xxx网站| 成人性生交大片免费视频hd| 中文字幕av成人在线电影| 国产精品美女特级片免费视频播放器| av国产免费在线观看| 夜夜躁狠狠躁天天躁| 国产伦人伦偷精品视频| 不卡一级毛片| 99久国产av精品| 亚洲av成人av| 日韩精品中文字幕看吧| 久99久视频精品免费| 午夜a级毛片| 天堂影院成人在线观看| 欧美黑人欧美精品刺激| 午夜免费观看网址| 人妻久久中文字幕网| 男女下面进入的视频免费午夜| 少妇的丰满在线观看| 久久人人精品亚洲av| 亚洲av不卡在线观看| 国内揄拍国产精品人妻在线| 一本精品99久久精品77| 美女高潮喷水抽搐中文字幕| 日日摸夜夜添夜夜添小说| 亚洲av日韩精品久久久久久密| 中文资源天堂在线| 丝袜美腿在线中文| 熟女电影av网| 给我免费播放毛片高清在线观看| 精品电影一区二区在线| 91麻豆av在线| 久久精品国产清高在天天线| 亚洲在线自拍视频| 90打野战视频偷拍视频| 岛国在线观看网站| 男人的好看免费观看在线视频| 国产精品综合久久久久久久免费| av片东京热男人的天堂| eeuss影院久久| av黄色大香蕉| 女生性感内裤真人,穿戴方法视频| 日日摸夜夜添夜夜添小说| 日本免费a在线| 亚洲色图av天堂| 亚洲最大成人手机在线| 此物有八面人人有两片| 欧美3d第一页| av欧美777| 天堂网av新在线| 最近视频中文字幕2019在线8| 色老头精品视频在线观看| 国产精华一区二区三区| 欧美一区二区精品小视频在线| 首页视频小说图片口味搜索| 婷婷六月久久综合丁香| 怎么达到女性高潮| 国产精品精品国产色婷婷| 在线十欧美十亚洲十日本专区| 制服丝袜大香蕉在线| 国产欧美日韩精品一区二区| 制服丝袜大香蕉在线| av福利片在线观看| 国产亚洲av嫩草精品影院| 国内少妇人妻偷人精品xxx网站| 久久久久性生活片| 久久精品夜夜夜夜夜久久蜜豆| 精品一区二区三区av网在线观看| 90打野战视频偷拍视频| 在线十欧美十亚洲十日本专区| 免费观看的影片在线观看| 大型黄色视频在线免费观看| 成年免费大片在线观看| 国产高清激情床上av| 淫妇啪啪啪对白视频| 久久久久久国产a免费观看| 香蕉久久夜色| 日韩欧美精品免费久久 | 亚洲无线观看免费| 一夜夜www| 国产亚洲精品久久久久久毛片| 国产精品美女特级片免费视频播放器| 欧美av亚洲av综合av国产av| 日韩欧美在线二视频| 日本精品一区二区三区蜜桃| 好看av亚洲va欧美ⅴa在| 欧美激情在线99| 悠悠久久av| 国产v大片淫在线免费观看| 免费电影在线观看免费观看|