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

    RTM-GWAS方法應(yīng)用于大豆RIL群體百粒重QTL檢測的功效

    2020-06-03 08:02:36潘麗媛賀建波趙晉銘王吳彬邢光南喻德躍張小燕李春燕陳受宜蓋鈞鎰
    中國農(nóng)業(yè)科學(xué) 2020年9期
    關(guān)鍵詞:表型變異遺傳

    潘麗媛,賀建波,趙晉銘,王吳彬,邢光南,喻德躍,張小燕,李春燕,陳受宜,蓋鈞鎰

    RTM-GWAS方法應(yīng)用于大豆RIL群體百粒重QTL檢測的功效

    潘麗媛1,賀建波1,趙晉銘1,王吳彬1,邢光南1,喻德躍1,張小燕3,李春燕3,陳受宜2,蓋鈞鎰1

    (1南京農(nóng)業(yè)大學(xué)大豆研究所/國家大豆改良中心/農(nóng)業(yè)部大豆生物學(xué)與遺傳育種重點實驗室/作物遺傳與種質(zhì)創(chuàng)新國家重點實驗室/江蘇省現(xiàn)代作物 生產(chǎn)協(xié)同創(chuàng)新中心,南京 210095;2中國科學(xué)院遺傳發(fā)育研究所/植物基因組學(xué)國家重點實驗室,北京 100101;3山東圣豐種業(yè)科技有限公司,山東嘉祥 272400)

    【】為全面解析大豆重組自交系群體中調(diào)控百粒重性狀的QTL體系,將限制性兩階段多位點全基因組關(guān)聯(lián)分析方法(RTM-GWAS)和不同定位方法進行比較、優(yōu)選,為后續(xù)候選基因體系探索及分子標(biāo)記輔助育種設(shè)計提供依據(jù)。利用以科豐1號和南農(nóng)1138-2為親本衍生的重組自交系群體NJRIKY的427個家系,通過由全基因組39 353個SNP構(gòu)建的3 683個SNPLDB標(biāo)記及3個環(huán)境下的百粒重表型數(shù)據(jù),選用復(fù)合區(qū)間作圖法(CIM)、基于混合線性模型的全基因組關(guān)聯(lián)分析方法(MLM-GWAS)和RTM-GWAS 3種方法檢測百粒重QTL,通過QTL數(shù)目和總的表型變異解釋率比較檢測功效,挑選最佳定位結(jié)果進行NJRIKY群體中的百粒重遺傳體系解析。通過候選基因體系的功能注釋,挖掘調(diào)控大豆百粒重的生物學(xué)途徑。科豐1號與南農(nóng)1138-2的百粒重差異較大,多環(huán)境平均數(shù)分別為9.0和17.9 g,遺傳變異系數(shù)為12.4%,遺傳率為85.4%,適用于百粒重性狀的遺傳解析。比較3種方法定位結(jié)果表明RTM-GWAS方法表現(xiàn)最佳,檢測QTL數(shù)目最多(57個),解釋表型變異最多(70.78%)。而CIM僅檢測到14個QTL,解釋了56.47%的表型變異,MLM-GWAS僅定位到6個QTL,解釋了18.47%的表型變異。RTM-GWAS共檢測到57個QTL,分布在19條染色體上,表型變異解釋率為0.03%—7.57%,其中41個QTL覆蓋了已報道的來自30個雙親群體的81個百粒重QTL,16個QTL為新發(fā)現(xiàn)位點,包含一個表型變異解釋率大于3%的大效應(yīng)位點。此外,檢測的57個QTL中有20個位點與環(huán)境存在互作效應(yīng)。這57個QTL構(gòu)成了影響NJRIKY群體百粒重性狀的遺傳體系。通過SNPLDB標(biāo)記與預(yù)測基因內(nèi)的SNP進行χ2檢驗,共篩選到36個候選基因,其中4個候選基因來自大效應(yīng)QTL,剩余32個候選基因來自小效應(yīng)QTL。通過GO注釋發(fā)現(xiàn)這些候選基因功能注釋豐富,其中13個候選基因與籽粒發(fā)育直接相關(guān),剩余的候選基因功能豐富,包含轉(zhuǎn)運、轉(zhuǎn)錄調(diào)節(jié)因子等,表明不同生物學(xué)途徑的基因共同調(diào)控NJRIKY群體中百粒重性狀的表達。3種定位方法中,高效的RTM-GWAS方法檢測到較為全面的NJRIKY群體的百粒重QTL,更適用于雙親RIL群體的QTL定位。不同功能的候選基因共同調(diào)控了復(fù)雜的百粒重性狀的表達。

    大豆;百粒重;QTL;重組自交家系;限制性兩階段多位點全基因組關(guān)聯(lián)分析

    0 引言

    【研究意義】大豆是人類植物蛋白和脂肪的主要食物來源[1],產(chǎn)量一直是大豆生產(chǎn)發(fā)展的關(guān)注點,而百粒重是大豆產(chǎn)量構(gòu)成的重要因子,一般與產(chǎn)量呈正相關(guān)[2]。百粒重作為數(shù)量性狀,具有復(fù)雜的遺傳基礎(chǔ),而現(xiàn)有的連鎖定位方法,僅能定位到少數(shù)幾個位點,限制了其遺傳體系的全面解析,需要高效的定位方法進行全基因組QTL檢測以全面解析其遺傳基礎(chǔ)?!厩叭搜芯窟M展】前人基于連鎖定位的復(fù)合區(qū)間作圖法(composite interval mapping,CIM)[3]對多個雙親群體進行了百粒重性狀的QTL定位分析。目前,SoyBase數(shù)據(jù)庫(https://www.soybase.org/,截至2017年3月)已經(jīng)收錄了基于40個不同群體、利用不同定位方法檢測到的250個大豆百粒重QTL。最近,F(xiàn)ujii等[4]利用一個包含181個家系的重組自交系群體(recombinant inbred line,RIL),在2個環(huán)境下利用CIM方法共檢測到5個百粒重QTL,分布于第12、13、17和20等4條染色體上。Zhang等[5]通過半野生豆花皮豆(Huapidou)和一個栽培豆齊黃26(Qihuang26)衍生的RIL群體,利用完備區(qū)間作圖法(inclusive composite interval mapping,ICIM)和CIM分別檢測到10個和5個百粒重QTL。由于連鎖定位往往僅能定位到少數(shù)幾個大效應(yīng)位點,并且置信區(qū)間較大,不利于進一步篩選候選基因?,F(xiàn)今,眾多研究者通過關(guān)聯(lián)分析的方法,利用自然群體進行大豆百粒重遺傳解析[6-13]。Copley等[14]利用86個早熟組(MG 000—MG 00)大豆材料在2個環(huán)境下僅檢測到5個百粒重QTL。Jing等[15]通過由地方種質(zhì)和育成品種構(gòu)成的包含185個大豆材料的自然群體檢測到20個與百粒重關(guān)聯(lián)的SNP位點,其中11個SNP連續(xù)分布在同一基因組區(qū)域。為了解決現(xiàn)有GWAS方法中遺傳率缺失(missing heritability)和自然群體中的復(fù)等位變異問題,He等[16]開發(fā)了一種新的限制性兩階段多位點關(guān)聯(lián)分析方法(restricted two-stage multi-locus genome- wide association study,RTM-GWAS),該方法通過遺傳率控制位點的總表型變異解釋率,利用二階段方法和多位點模型提高了檢測功效,并成功應(yīng)用于多個自然群體和多親本群體[17-22]?!颈狙芯壳腥朦c】利用連鎖定位方法檢測到的大豆百粒重QTL較少,表型變異解釋率丟失較多。而隨著分子標(biāo)記信息的加密,關(guān)聯(lián)分析的精度提高,并且雙親群體中的群體結(jié)構(gòu)簡單[23],前人已多次成功將關(guān)聯(lián)分析方法應(yīng)用于雙親群體中[24-26],選擇高效的關(guān)聯(lián)分析方法有助于大豆百粒重數(shù)量性狀遺傳基礎(chǔ)的全面解析,可為揭示其完整基因體系提供支持?!緮M解決的關(guān)鍵問題】本研究基于科豐1號×南農(nóng)1138-2衍生的重組自交系NJRIKY群體,利用3個環(huán)境下的大豆百粒重表型數(shù)據(jù)和3 683個SNPLDB標(biāo)記,使用CIM、MLM-GWAS和RTM-GWAS 3種定位方法檢測百粒重QTL,以QTL數(shù)目及總的表型變異解釋率判定其檢測功效,篩選最佳定位方法,進一步解析NJRIKY群體中百粒重的全基因組QTL體系及候選基因體系。

    1 材料與方法

    1.1 供試材料

    以栽培大豆科豐1號和南農(nóng)1138-2為母本和父本,衍生了一個含有427個家系的重組自交系(recombinant inbred line,RIL)群體。群體命名為NJRIKY,由兩部分組成;其中184個家系創(chuàng)制于1994年夏,F(xiàn)2世代采用單粒傳法繁代,F(xiàn)3-F7家系內(nèi)混合收獲播種,F(xiàn)7家系內(nèi)隨機收獲1株,種植F2:8世代,最終經(jīng)過群體校正后獲得;另外243個家系采用同樣方法獲得[27]。親本科豐1號來自黃淮海二熟制春夏作大豆品種生態(tài)區(qū),由中科院遺傳研究所選育,熟期組為MG II,黑色種皮;親本南農(nóng)1138-2來自長江中下游二熟制春夏作大豆品種生態(tài)區(qū),由南京農(nóng)業(yè)大學(xué)選育,熟期組為MG V,黃色種皮[28]。

    1.2 田間試驗

    NJRIKY群體及親本,2012—2013年夏季于江蘇省南京市南京農(nóng)業(yè)大學(xué)江浦試驗站(32.07°N, 118.62°E)和2013年夏季于山東省濟寧市圣豐試驗站(36.67°N,116.98°E)3個環(huán)境進行田間試驗。3個環(huán)境分別簡稱為12JP、13JP和13SF,代表2012年江浦試驗站、2013年江浦試驗站和2013年圣豐試驗站。田間試驗均采用立方格子設(shè)計,每個環(huán)境3次重復(fù),江浦試驗站采用穴播,兩穴一個小區(qū),穴距為0.7 m×0.8 m,每穴定苗6株;圣豐試驗站采用行播,行長2 m,行距0.5 m,常規(guī)田間管理,籽粒成熟后按小區(qū)分別收獲脫粒,并在35—40℃條件下烘干至恒重,取100個完整籽粒置于精度0.01 g天平稱量記錄。

    1.3 標(biāo)記分型

    利用427家系及2個親本的新鮮幼嫩葉片采用CTAB法[29]提取DNA,建庫并利用Illumina Hiseq2000進行雙端測序,通過與參考基因組Williams82 v1.1比對獲得SNP。按照缺失≤20%、錯誤率≤1%和雜合率<5%的標(biāo)準過濾獲得39 353個高質(zhì)量SNP。利用NPUTE[30]對SNP缺失基因型進行填補,用于SNPLDB(SNP linkage disequilibrium block)標(biāo)記構(gòu)建[25]。利用賀建波等[16]開發(fā)的RTM-GWAS軟件,根據(jù)雙親帶型劃分SNPLDB標(biāo)記,設(shè)置的最大窗口為200 kb使用默認設(shè)置的置信區(qū)間法定義區(qū)段,計算標(biāo)記間的LD(D’),若95%以上的SNP的D’值位于0.70—0.98,則說明具有強LD水平,則被劃分到一個單體型區(qū)塊中,這樣全基因組的SNP被劃分為SNP連鎖不平衡區(qū)塊(SNP linkage disequilibrium block,SNPLDB)。根據(jù)雙親帶型進行SNPLDB基因型分型,少量非親本單體型被與其最相似的親本單體型替換,最終獲得3 683個SNPLDB標(biāo)記用于遺傳圖譜構(gòu)建及QTL定位[25]。

    1.4 表型數(shù)據(jù)分析

    采用SAS軟件PROC GLM進行多環(huán)境下表型數(shù)據(jù)的聯(lián)合方差分析,同時利用PROC MIXED的REML方法估計方差組分,根據(jù)和分別估計多環(huán)境和單環(huán)境下的性狀遺傳率,遺傳變異系數(shù)(genotypic coefficient of variation,GCV)為100×σ/,其中為遺傳方差,為基因型與環(huán)境互作方差,為誤差方差,為環(huán)境的個數(shù),為一個環(huán)境內(nèi)的試驗重復(fù)數(shù),為整體平均數(shù)。

    1.5 QTL定位方法

    利用3種方法進行NJRIKY群體的百粒重QTL定位,分別為復(fù)合區(qū)間作圖法(CIM)[3]、混合線性模型全基因組關(guān)聯(lián)分析方法(mixed linear model GWAS或MLM-GWAS)[31]、限制性兩階段多位點全基因組關(guān)聯(lián)分析方法(RTM-GWAS)[16]。

    CIM方法應(yīng)用Windows QTL Cartographer V2.5軟件,表型數(shù)據(jù)為多環(huán)境下各家系百粒重均值,前景及背景選擇的顯著水平為0.05,步長為1 cM,并進行排列測驗1 000次,確定QTL檢測LOD閾值。

    MLM-GWAS方法[31]應(yīng)用TASSEL 3.0軟件[32],表型數(shù)據(jù)為多環(huán)境下各家系百粒重均值,參數(shù)設(shè)定為默認參數(shù),采用5個主成分(principal component,P)和個體間親緣關(guān)系(kinship,K)矩陣,即P+K模型的MLM進行關(guān)聯(lián)分析,并利用R軟件fdrtool包[33]進行假發(fā)現(xiàn)率(false discovery rate,F(xiàn)DR)多重測驗校正,顯著水平為0.05。

    RTM-GWAS方法利用限制性二階段多位點全基因組關(guān)聯(lián)分析方法(RTM-GWAS)軟件[16]進行QTL定位,表型數(shù)據(jù)為多環(huán)境小區(qū)原始觀測值,顯著水平設(shè)為0.05,分析模型同時考慮QTL與環(huán)境互作效應(yīng)。根據(jù)表型變異解釋率(phenotypic variation explained,PVE)大小將QTL分為兩類:大貢獻(large contribution,LC)主效位點(PVE≥3%);小貢獻(small contribution,SC)主效位點(PVE<3%)。

    1.6 候選基因體系注釋

    首先,檢測SNPLDB標(biāo)記(±100 kb)物理區(qū)間內(nèi)的注釋基因;然后利用卡平方(χ2)獨立性測驗,顯著水平設(shè)為0.05,檢測SNPLDB標(biāo)記與注釋基因內(nèi)SNP之間的相關(guān)性。當(dāng)這二者顯著相關(guān)時,則該SNP所在的注釋基因被確定為候選基因;同時,利用SoyBase數(shù)據(jù)庫(https://www. soybase.org)提取所有注釋基因的功能注釋(gene ontology,GO),優(yōu)先選擇與籽粒發(fā)育相關(guān)的注釋基因確定為候選基因。同時給出候選基因的功能注釋(GO)結(jié)果。

    2 結(jié)果

    2.1 多環(huán)境下NJRIKY群體百粒重變異

    從表1可以看出3個環(huán)境下NJRIKY群體的親本科豐1號和南農(nóng)1138-2的表型差異較大,多環(huán)境平均值分別為9.0和17.9 g,重組自交家系的百粒重變異范圍為8.6—17.1 g(圖1),遺傳變異系數(shù)為12.4%,變異范圍位于雙親之間。該群體百粒重的遺傳率較高,多環(huán)境聯(lián)合下為85.4%,單環(huán)境下為92.5%—93.9%。多環(huán)境聯(lián)合方差分析的結(jié)果表明家系間、環(huán)境間及家系與環(huán)境互作的差異均達到極顯著水平(表2)。

    Mean表示各家系3個環(huán)境百粒重平均數(shù)

    表1 多環(huán)境下NJRIKY群體百粒重的次數(shù)分布和描述統(tǒng)計

    Table 1 Frequency distribution and descriptive statistics of 100-seed weight of NJRIKY under multiple environments

    :遺傳變異系數(shù): genotypic coefficient of variation

    表2 NJRIKY群體多環(huán)境聯(lián)合方差分析結(jié)果

    2.2 不同定位方法下百粒重QTL檢測結(jié)果

    利用3 683個SNPLDB標(biāo)記及CIM、MLM-GWAS和RTM-GWAS 3種定位方法對NJRIKY群體的百粒重性狀進行QTL定位(表3),并通過QTL的數(shù)目和總的表型變異解釋率來比較各個方法的檢測功效,選取最佳的定位結(jié)果進行該群體的百粒重遺傳體系解析。

    表3 NJRIKY群體中利用不同定位方法檢測到的百粒重位點結(jié)果概要

    檢測的QTL:本研究檢測的QTL數(shù)目和染色體數(shù)目,57(19):57個QTL位于19條染色體上。Mapped QTL:所有定位到的QTL的匯總。大貢獻主效QTL和小貢獻主效QTL分別表示貢獻率≥3%與<3%的主效QTL。QTL×Env.表示QTL與環(huán)境互作效應(yīng)。已報道的QTL:SoyBase (http://www.soybase.org) 數(shù)據(jù)庫中前人報道的與始花期位于鄰近1 Mb位置的QTL數(shù)目,41(84)表示41個QTL與SoyBase (http://www.soybase.org) 中84個QTL位置一致或鄰近

    Detected QTLs: the number of detected QTLs and chromosomes in the present study, 57(19): 57 QTLs on 19 chromosomes. Mapped QTL: the sum of all major QTLs. LC major QTL: the QTL detected from RTM-GWAS with its PVE more than 3% is called large-contribution (LC) major QTLs; SC major QTL: small-contribution major QTL with PVE less than 3%. QTL × Env.: the interaction between QTL and environments. Reported QTLs: the number of QTLs reported in SoyBase (http://www.soybase.org) that is nearby the present detected QTL in the RTM-GWAS procedure according to the physical position within 1 Mb, 41 (84): 41 QTLs shared same confidence regions with 84 SoyBase QTLs

    在連鎖定位中較為普遍使用的CIM方法下共檢測到14個百粒重QTL,分布在8個連鎖群上,總表型變異解釋率(PVE)為56.47%(圖2和表4)。在普遍使用的關(guān)聯(lián)分析方法MLM-GWAS中僅檢測到6個QTL,位于3條染色體上,總PVE為18.47%(圖2和表4)。利用RTM-GWAS方法,檢測到57個百粒重QTL,分布在除Gm19以外的19條染色體上,共解釋70.78%的表型變異,結(jié)合方差分析估計得到的遺傳率為85.4%,剩余14.52%的表型變異來源于未定位到的微效QTL(圖2)。其中位于13條染色體上的20個百粒重QTL與環(huán)境存在互作效應(yīng),共解釋4.20%的表型變異(表3和表5)。

    遵循常規(guī),以位點前后1 Mb的范圍來歸并QTL,發(fā)現(xiàn)RTM-GWAS方法幾乎覆蓋了CIM(8個QTL/共14個QTL,PVE=40.85%/總PVE=56.47%)和MLM-GWAS(4個QTL/共6個QTL,PVE=12.96%/總PVE=18.47%)方法的結(jié)果,包括位于第11染色體上的效應(yīng)最大的位點Gm11_BLOCK_5228756_ 5269867,這個位點在3種方法中的表型變異解釋率均最大,分別為11.05%(CIM)、3.69%(MLM-GWAS)和7.57%(RTM-GWAS)。

    表4 NJRIKY群體中在CIM和MLM-GWAS方法下的百粒重QTL

    LOD/-lg:LOD為CIM方法中QTL檢測的LOD(logarithm of odds)得分,-lg為MLM-GWAS方法中QTL檢測的值的對數(shù)值。SoyBase QTL:SoyBase (http://soybase.org) 數(shù)據(jù)庫中前人報道的與始花期位于鄰近1 Mb位置的QTL名稱

    LOD/-lg: LOD is the logarithm of odds score in CIM and -lgis thevalue in log10 scale in MLM-GWAS. SoyBase QTL: the QTL names reported in SoyBase (http://soybase.org) that is nearby the present detected QTL in the RTM-GWAS procedure according to the physical position within 1 Mb

    表5 NJRIKY群體中與百粒重性狀關(guān)聯(lián)的QTL位點

    *表示該QTL具有QTL與環(huán)境互作效應(yīng)。標(biāo)記中Gm后的數(shù)字表示其所在染色體,末端數(shù)字為其對應(yīng)物理位置

    * indicates the QTL interacted with environments. The number after “Gm” means the chromosome number, the number at the end of a marker indicates the corresponding physical position

    以上3種定位方法的百粒重QTL結(jié)果表明,RTM-GWAS方法的檢測功效最高,檢測到更多的QTL(57 vs. 14和6),解釋了更多的表型變異(70.78% vs. 56.47%和18.47%),在NJRIKY群體中檢測到了較為全面的百粒重性狀的遺傳體系,所以將對RTM-GWAS方法檢測到的由57個QTL構(gòu)成的影響NJRIKY群體百粒重性狀的遺傳體系進行進一步的遺傳解析。這里要特別指出,RTM-GWAS還檢測到57個QTL中有20個是與環(huán)境存在互作(或因環(huán)境而波動)的QTL,這是其他2種定位方法所不及的。

    2.3 NJRIKY RIL群體中百粒重的QTL及候選基因體系

    利用3 683個SNPLDB標(biāo)記及RTM-GWAS方法對NJRIKY群體百粒重進行全基因QTL檢測,第一階段中,利用一般線性模型(general linear model,GLM)篩選到1 860個候選SNPLDB標(biāo)記;第二階段中,多位點模型最終檢測到57個與百粒重相關(guān)的QTL(圖3和表5)。

    a:染色體結(jié)構(gòu),異染色質(zhì)區(qū)域標(biāo)為淺紅色(單位為Mb)

    圖3 NJRIKY中百粒重關(guān)聯(lián)分析的曼哈頓圖和QQ圖

    57個百粒重QTL分布于19條染色體,表型變異解釋率變化范圍為7.57%—0.13%,其中35個QTL的PVE為0.13%—1.00%,17個為1.00%—3.00%,剩余5個為3.00%—7.57%,表明百粒重作為一個復(fù)雜的數(shù)量性狀,其所包含的位點貢獻率變化范圍廣,不僅僅包含表型變異解釋率很大的QTL,還包含更多小貢獻率的QTL。

    57個QTL構(gòu)成的調(diào)控NJRIKY群體百粒重的QTL體系中,第4染色體包含的QTL數(shù)目最多,共14個QTL,第20染色體包含的QTL最少,僅為1個QTL。、、、、分別解釋了7.57%、5.12%、3.72%和3.61%的表型變異(phenotypic variation,PV),為大貢獻主效QTL,共解釋了23.30%的PV,剩余52個小貢獻主效QTL共解釋了47.48%的PV。NJRIKY群體的QTL體系中41個QTL與SoyBase中收錄的包含52個親本的30個雙親群體中檢測到的81個QTL(共250個QTL)的位置一致(表3)。

    RTM-GWAS在3個環(huán)境下還檢測到20個QTL與環(huán)境存在互作效應(yīng),其中,3個位于第4、9和11染色體的大貢獻主效QTL和17個位于12條染色體上的小貢獻主效QTL(表5)。

    通過RTM-GWAS方法對NJRIKY群體的百粒重性狀進行較為全面的檢測后,基于高密度的遺傳標(biāo)記,進一步推斷其候選基因體系。通過SNPLDB標(biāo)記與預(yù)測基因內(nèi)的SNP進行χ2檢驗,共在57個QTL中檢測到36個候選基因解釋了54.12%的PV,其中4個候選基因來自4個大貢獻主效QTL,剩余32個候選基因來自32個小貢獻主效位點,還有20個QTL中未檢測到符合條件的候選基因(存在2個QTL的候選基因相同)(表6)。通過GO注釋可以看出13個候選基因與籽粒發(fā)育相關(guān),剩余的候選基因的功能變化豐富,包含了轉(zhuǎn)運、蛋白糖基化、轉(zhuǎn)錄調(diào)節(jié)因子、脫落酸響應(yīng)等,表明這些候選基因通過不同生物途徑共同調(diào)控百粒重性狀的表達。

    3 討論

    3.1 高效的RTM-GWAS方法適用于雙親群體的QTL定位

    隨著測序技術(shù)的快速發(fā)展,植物生物學(xué)家的研究重心逐漸從少數(shù)基因功能向挖掘調(diào)控目標(biāo)性狀的全面基因體系和基因網(wǎng)絡(luò)轉(zhuǎn)移[34]。前人結(jié)果表明,RTM-GWAS方法適用于自然群體和巢式關(guān)聯(lián)群體的數(shù)量性狀全基因組位點檢測[16-17, 19]。本研究在3種不同的QTL方法中,RTM-GWAS幾乎覆蓋了其他方法(CIM和MLM-GWAS)檢測到的所有QTL。RTM-GWAS檢測到較多的QTL(57個QTL)解釋了更多的表型變異(70.78%),而在CIM和MLM-GWAS檢測到的QTL僅解釋了56.47%和18.47%的表型變異。由于雙親群體中不存在復(fù)雜的群體結(jié)構(gòu)問題,將關(guān)聯(lián)分析的方法應(yīng)用到雙親群體中,提高了檢測精度,將大大加快雙親群體數(shù)量性狀的遺傳體系的全面解析。另外,利用RTM-GWAS方法,由于其通過遺傳率控制了位點的總表型變異解釋率,利用二階段和多位點模型提高了檢測功效,僅利用一個雙親群體就檢測到了41個QTL與已報道的52個大豆材料中的81個QTL位置一致,并且13個候選基因與籽粒發(fā)育相關(guān),表明RTM-GWAS方法在雙親群體RIL群體中QTL定位表現(xiàn)高效,適用于雙親群體的QTL定位。

    3.2 大豆百粒重性狀的遺傳體系

    大豆是人類重要的植物蛋白和油脂的食物來源,大豆百粒重與產(chǎn)量的直接相關(guān)[2],直接決定了大豆百粒重性狀的重要性。前人已經(jīng)通過不同類型的群體對大豆百粒重的遺傳體系進行了研究[18, 35],本研究利用RTM-GWAS方法在親本具有較大差異的NJRIKY群體中檢測到57個百粒重相關(guān)聯(lián)的QTL,其中16個QTL是該群體檢測到的新位點,的PVE達到5.12%。5個大效應(yīng)主效QTL中,的貢獻率最高為7.57%,其位置與SoyBase中的、和一致,其中和來源于科豐1號和南農(nóng)1138-2的184個家系后代群體[36-37],QTL區(qū)間較大(5 132 968—6 184 217 bp),群體擴大及標(biāo)記加密后,SNPLDB的位置為5 228 756-5 269 867 bp。的貢獻率為3.71%,與已報道的和位于相近位置。的PVE為3.60%,與已報道的3個QTL位置相近,其中的定位群體為NJRIKY群體的184個家系。的貢獻率為3.27%,與和[38-39]的位置相近,這兩個已報道位點的標(biāo)記區(qū)間極大,約7 Mb,可能是定位的群體大小及標(biāo)記密度的影響。其他小貢獻主效QTL中也不乏被多個群體多次定位到的位點,如的貢獻率在本群體中僅為0.68%,其與SoyBase中的7個雙親群體中的8個QTL位置相近。

    表6 RTM-GWAS方法中檢測到的與百粒重相關(guān)的候選基因體系

    加粗的候選基因表示該基因的GO注釋與籽粒發(fā)育相關(guān)

    The candidate genes in boldface mean their gene ontology description is associated with seed development

    NJRIKY群體中檢測到36個候選基因中,13個候選基因與籽粒發(fā)育相關(guān)。其中2個候選基因、分別來自2個大貢獻主效QTL——和,剩余11個候選基因來自11個小貢獻主效位點。如,(PVE=1.87%)的候選基因與胚發(fā)育相關(guān),其在擬南芥中的同源基因表達產(chǎn)物為TOR蛋白,直接調(diào)控了植物的胚生長發(fā)育和逆境響應(yīng)[40]。具有不同功能的候選基因,共同調(diào)控了復(fù)雜的百粒重性狀的表達,而現(xiàn)今的研究還具有局限性,包括材料的限制等,因此,綜合越來越多的群體中百粒重性狀的結(jié)果,結(jié)合高效的QTL定位方法,為全面解析基因體系及網(wǎng)絡(luò)提供了必要基礎(chǔ)。

    4 結(jié)論

    NJRIKY RIL群體中,雙親的百粒重差異較大,適用于百粒重性狀的遺傳解析。利用3種定位方法,在3個環(huán)境下,RTM-GWAS方法定位到的QTL更多,表型變異解釋率更高,其檢測功效優(yōu)于其他2種方法(CIM和MLM-GWAS)。通過RTM-GWAS定位到57個百粒重QTL構(gòu)成了NJRIKY群體的百粒重遺傳體系,表明高效的QTL定位方法將加快全面解析復(fù)雜的百粒重性狀的遺傳體系。

    [1] SMITH T J, CAMPER H M. Effect of seed size on soybean performance., 1970, 67(5): 681-684.

    [2] BURRIS J S, EDJE O T, WAHAB A H. Effects of seed size on seedling performance in soybeans: II. seedling growth and photosynthesis and field performance1., 1973, 13(2): 207.

    [3] ZENG Z. Precision mapping of quantitative trait loci., 1994, 136(4): 1457-1468.

    [4] FUJII K, SAYAMA T, TAKAGI K, TAKAGI K, KOSUGE K, OKANO K, KAGA A, ISHINOTO M. Identification and dissection of single seed weight QTLs by analysis of seed yield components in soybean., 2018, 68(2): 177-187.

    [5] ZHANG Y, LI W, LIN Y, ZHANG L, WANG C, XU R. Construction of a high-density genetic map and mapping of QTLs for soybean () agronomic and seed quality traits by specific length amplified fragment sequencing., 2018, 19(1): 641.

    [6] ZHOU Z, JIANG Y, WANG Z, GOU Z, LYU J, LI W, YU Y, SHU L, ZHAO Y, MA Y, FANG C, SHEN Y, LIU T, LI C, LI Q, WU M, WANG M, WU Y, DONG Y, WAN W, WANG X, DING Z, GAO Y, XIANG H, ZHU B, LEE S H, WANG W, TIAN Z. Resequencing 302 wild and cultivated accessions identifies genes related to domestication and improvement in soybean., 2015, 33(4): 408-414.

    [7] ZHANG H, HAO D, SITOE H M, YIN Z, HU Z, ZHANG G, YU D, SINGH R. Genetic dissection of the relationship between plant architecture and yield component traits in soybean () by association analysis across multiple environments., 2015, 134(5): 564-572.

    [8] CONTRERAS-SOTO R I, MORA F, DE OLIVEIRA M A, HIGASHI, W, SCAPIM C A, SCHUSTER I A. genome-wide association study for agronomic traits in soybean using SNP markers and SNP-based haplotype analysis., 2017, 12(2): e0171105.

    [9] HU Z, ZHANG D, ZHANG G, KAN G, HONG D, YU D. Association mapping of yield-related traits and SSR markers in wild soybean (Sieb. and Zucc.)., 2014, 63(5): 441-449.

    [10] ZHANG J, SONG Q, CREGAN P B, JIANG G. Genome-wide association study, genomic prediction and marker-assisted selection for seed weight in soybean ()., 2016, 129(1): 117-130.

    [11] SONAH H, O'DONOUGHUE L, COBER E, RAJCAN I, BELZILE F. Identification of loci governing eight agronomic traits using a GBS-GWAS approach and validation by QTL mapping in soya bean., 2015, 13(2): 211-221.

    [12] YAN L, HOFMANN N, LI S, FERREIRA M E, SONG B, JIANG G, REN S, QUIGLEY C, FICKUS E, GREGAN P, SONG Q. Identification of QTL with large effect on seed weight in a selective population of soybean with genome-wide association and fixation index analyses., 2017, 18(1): 529.

    [13] WANG J, CHU S, ZHANG H, ZHU Y, CHENG HAO, YU D. Development and application of a novel genome-wide SNP array reveals domestication history in soybean., 2016, 6: 20728.

    [14] COPLEY T R, DUCEPPE M O, O'DONOUGHUE L S. Identification of novel loci associated with maturity and yield traits in early maturity soybean plant introduction lines., 2018, 19(1): 167.

    [15] JING Y, ZHAO X, WANG J, TENG W, QIU L, HAN Y, LI W. Identification of the genomic region underlying seed weight per plant in soybean (L. Merr.) via high-throughput single- nucleotide polymorphisms and a genome-wide association study., 2018, 9: 1392.

    [16] HE J, MENG S, ZHAO T, XING G, YANG S, LI Y, GUANG R, LU J, WANG Y, XIA Q, YANG B, GAI J. An innovative procedure of genome-wide association analysis fits studies on germplasm population and plant breeding., 2017, 130(11): 2327-2343.

    [17] LI S, CAO Y, HE J, ZHAO T, GAI J. Detecting the QTL-allele system conferring flowering date in a nested association mapping population of soybean using a novel procedure., 2017, 130(11): 2297-2314.

    [18] ZHANG Y, HE J, WANG Y, XING G, ZHAO J, LI Y, YANG S, PALMER R G, ZHAO T, GAI J. Establishment of a 100-seed weight quantitative trait locus-allele matrix of the germplasm population for optimal recombination design in soybean breeding programmes., 2015, 66(20): 6311-6325.

    [19] MENG S, HE J, ZHAO T, XING G, LI Y, YANG S, LU J, WANG Y, GAI J. Detecting the QTL-allele system of seed isoflavone content in Chinese soybean landrace population for optimal cross design and gene system exploration., 2016, 129(8): 1557-1576.

    [20] ZHANG Y, HE J, WANG H, WANG H, MENG S, XING G, LI Y, YANG S, ZHAO J, ZHAO T, GAI J. Detecting the QTL-allele system of seed oil traits using multi-locus genome-wide association analysis for population characterization and optimal cross prediction in soybean., 2018, 9: 1793.

    [21] ZHANG Y, HE J, MENG S, LIU M, XING G, LI Y, YANG S, YANG J, ZHAO T, GAI J. Identifying QTL-allele system of seed protein content in Chinese soybean landraces for population differentiation studies and optimal cross predictions., 2018, 214(9),157.

    [22] KHAN M A, TONG F, WANG W, HE J, ZHAO T, GAI J. Analysis of QTL-allele system conferring drought tolerance at seedling stage in a nested association mapping population of soybean [(L.) Merr.] using a novel GWAS procedure., 2018, 248(4): 947-962.

    [23] MACKAY I, POWELL W. Methods for linkage disequilibrium mapping in crops., 2007, 12(2): 57-63.

    [24] MALOSETTI M, VAN EEUWIJK F A, BOER M P, CASAS A M, ELIA M, MORALEJO M, BHAT P R, RAMSAY L, MOLINA CONO J L. Gene and QTL detection in a three-way barley cross under selection by a mixed model with kinship information using SNPs., 2011, 122(8): 1605-1616.

    [25] PAN L, HE J, ZHAO T, XING G, WANG Y, YU D, CHEN S, GAI J. Efficient QTL detection of flowering date in a soybean RIL population using the novel restricted two-stage multi-locus GWAS procedure., 2018, 131(12): 2581-2599.

    [26] RIYAN C, LIM J E, SAMOCHA K E, SOKOLOFF G, ABNEY M, SKOL A D, PALMER A A. Genome-wide association studies and the problem of relatedness among advanced intercross lines and other highly recombinant populations., 2010, 185(3): 1033.

    [27] 王永軍, 喻德躍, 章元明, 陳受宜, 蓋鈞鎰. 重組自交系群體的檢測調(diào)整方法及其在大豆NJRIKY群體的應(yīng)用. 作物學(xué)報, 2004, 30(5): 413-418.

    WANY Y J, YU D Y, ZHANG Y M, CHEN S Y, GAI J Y. Method of evaluation and adjustment of recombinant inbred line population and its application to the soybean RIL population NJRIKY.,2004, 30(5): 413-418. (in Chinese)

    [28] 蓋鈞鎰, 邱家馴, 趙團結(jié). 大豆品種南農(nóng)493-1和南農(nóng)1138-2與其衍生新品種的親緣關(guān)系及其育種價值分析. 南京農(nóng)業(yè)大學(xué)學(xué)報, 1997, 20 (1): 1-8.

    GAI J Y, QIU J X, ZHAO T J. An analysis of genetic relationship of nannong493-1 and nannong 1138-2 with their derivative cultivars and their potential in future breeding., 1997, 20(1): 1-8. (in Chinese)

    [29] MURRAY M G, THOMPSON W F. Rapid isolation of high molecular weight plant DNA., 1980, 8(19): 4321-4325.

    [30] ROBERTS A, MCMILLAN L, WANG W, PARKER J, RUSYN I, THREADGILL D. Inferring missing genotypes in large SNP panels using fast nearest-neighbor searches over sliding windows., 2007, 23(13): 401-407.

    [31] YU J M, PRESSOIR G, BRIGGS W H, BI I V, YAMASAKI M, DOEBLEY J F, MCMULLEN M D, GAUT B S, NIELSEN D M, HOLLAND J B, KRESOVICH S, BUCKLER E S. A unified mixed-model method for association mapping that accounts for multiple levels of relatedness., 2006, 38(2): 203-208.

    [32] BRADBURY P J, ZHANG Z, KROON D E, CASSTEVENS T M, RAMDOSS V, BUCKER E S. TASSEL: Software for association mapping of complex traits in diverse samples., 2007, 23(19): 2633-2635.

    [33] STRIMMER K. fdrtool: A versatile R package for estimating local and tail area-based false discovery rates., 2008, 24(12): 1461-1462.

    [34] GAI J. Quantitative Inheritance//MALOY S, HUGHES K. eds.. 2nd ed. San Diego: Academic Press, 2013: 18-21.

    [35] LU X, XIONG Q, CHENG T, LI Q, LIU X, BI Y, LI W, ZHANG W, MA B, LAI Y, DU W, MAN W, CHEN S, ZHANG J. A PP2C-1 allele underlying a quantitative trait locus enhances soybean 100-seed weight., 2017, 10(5): 670-684.

    [36] GAI J, WANG Y, WU X, CHEN S. A comparative study on segregation analysis and QTL mapping of quantitative traits in plants-with a case in soybean., 2007, 1(1): 1-7.

    [37] ZHANG W K, WANG Y J, LUO G Z, ZHANG J, HE C, WU X, GAI J, CHEN S. QTL mapping of ten agronomic traits on the soybean (L. Merr.) genetic map and their association with EST markers., 2004, 108(6): 1131-1139.

    [38] HAN Y, LI D, ZHU D, LI H, LI X, TENG W, LI W. QTL analysis of soybean seed weight across multi-genetic backgrounds and environments., 2012, 125(4): 671-683.

    [39] CHEN Q, ZHANG Z, LIU C, XIN D, QIU H, SHAN D, SHAN C, HU G. QTL analysis of major agronomic traits in soybean., 2007, 6(4): 399-405.

    [40] WANG P, ZHAO Y, LI Z, HSU C, LIU C, FU L, HOU Y, DU Y, XIE S, ZHANG C, GAO J, CAO M, HUANG X, ZHU Y, TANG K, WANG X, TAO W, XIONG Y, ZHU J. Reciprocal regulation of the TOR kinase and ABA receptor balances plant growth and stress response., 2018, 69(1): 100-112.

    Detection power of RTM-GWAS applied to 100-seed weight QTL identification in a recombinant inbred lines population of soybean

    PAN LiYuan1, HE JianBo1, ZHAO JinMing1, WANG WuBin1, XING GuangNan1, YU DeYue1, ZHANG XiaoYan3, LI ChunYan3, CHEN ShouYi2, GAI JunYi1

    (1Soybean Research Institute, Nanjing Agricultural University /National Center for Soybean Improvement /Key Laboratory of Biology and Genetic Improvement of Soybean (General), Ministry of Agriculture/State Key Laboratory for Crop Genetics and Germplasm Enhancement/Jiangsu Collaborative Innovation Center for Modern Crop Production, Nanjing 210095;2Institute of Genetics and Developmental Biology, Chinese Academy of Sciences/State Key Laboratory of Plant Genomics, Beijing 100101;3Shandong Shofine Seed Technology Co. Ltd., Jiaxiang 272400, Shandong)

    【】To thoroughly dissect the QTL system conferring 100-seed weight in a recombinant inbred lines population, the restricted two-stage multi-locus genome-wide association analysis (RTM-GWAS) method was compared with other mapping methods for method optimization, which will provides basis for further exploration of candidate gene system and molecular marker-assisted design breeding. 【】 A recombinant inbred line population consisting of 427 lines derived from a cross between Kefeng-1 and NN1138-2 was tested for its 100-seed weight under three environments. A total of 3 683 SNPLDBs (SNP linkage disequilibrium blocks) composed of 39 353 SNPs were applied to QTL mapping using three different mapping procedures, including the composite interval mapping (CIM) method, the mixed linear model (MLM-GWAS) method and the RTM-GWAS method, and the best mapping procedure was selected for the analysis of the 100-seed weight genetic system in NJRIKY population through comparing their detection power, including the detected number of QTLs and total phenotypic variation explained. 【】The 100-seed weight of Kefeng-1 and NN1138-2 were 9.0 g and 17.9 g, respectively, showing significant difference. The genotypic coefficient of variation and heritability of the trait were 12.4% and 85.4%, respectively. These results indicated that the population was suitable for genetic analysis of 100-seed weight trait. The RTM-GWAS procedure performed the best with the largest number of QTLs (57) explaining the most phenotypic variation (PVE=70.78%), while a total of 14 and 6 QTLs contributing 56.47% and 18.47% phenotypic variation were detected using CIM and MLM-GWAS, respectively. The 57 QTLs detected from the RTM-GWAS distributed on 19 chromosomes, of which 41 QTLs overlapped with 81 QTLs identified from 30 bi-parental populations in the literature. Furthermore, the PVE of 57 QTLs ranged from 0.03% to 7.57%, of which 16 QTLs were novel ones, including one large contribution major QTL(PVE>3%). Furthermore, 20 QTLs had significant interaction effect with environment. A total of 36 candidate genes were annotated from 37 QTLs through χ2test between SNPLDB markers and SNPs harboring in the predicted genes, of which 4 candidate genes were from the large contribution QTLs and other 32 candidate genes were from the small contribution QTLs. These candidate genes were included in different biological processes, of which 13 candidate genes were grouped in seed development directly, and the remaining candidate genes were grouped into different functions, such as transport, transcriptional regulators, etc., indicating that these genes from different biological pathways regulate the expression of 100-seed weight trait in NJRIKY together. 【】Among the three different mapping procedures, RTM-GWAS procedure is the most powerful one which can provide a relatively thorough detection of 100-seed weight QTLs in NJRIKY population, therefore, it is more suitable for QTL mapping in bi-parental population such as RIL population. The candidate genes with various functions jointly regulated the complex expression of 100-seed weight trait.

    soybean [(L.) Merr.]; 100-seed weight; QTL (quantitative trait locus); recombinant inbred lines population; restricted two-stage multi-locus genome-wide association analysis

    10.3864/j.issn.0578-1752.2020.09.004

    2019-08-24;

    2020-01-02

    國家自然科學(xué)基金(31701447)、國家作物育種重點研發(fā)計劃(2017YFD0101500,2017YFD0102002)、長江學(xué)者和創(chuàng)新團隊發(fā)展計劃(PCSIRT_17R55)、教育部111項目(B08025)、中央高?;究蒲袠I(yè)務(wù)費項目(KYT201801)、農(nóng)業(yè)部國家大豆產(chǎn)業(yè)技術(shù)體系CARS-04、江蘇省優(yōu)勢學(xué)科建設(shè)工程專項、江蘇省JCIC-MCP項目

    潘麗媛,E-mail:panly89@126.com。通信作者賀建波,E-mail:hjbxyz@gmail.com。通信作者蓋鈞鎰,E-mail:sri@njau.edu.cn

    (責(zé)任編輯 李莉)

    猜你喜歡
    表型變異遺傳
    非遺傳承
    變異危機
    變異
    還有什么會遺傳?
    還有什么會遺傳
    還有什么會遺傳?
    建蘭、寒蘭花表型分析
    GABABR2基因遺傳變異與肥胖及代謝相關(guān)表型的關(guān)系
    變異的蚊子
    百科知識(2015年18期)2015-09-10 07:22:44
    慢性乙型肝炎患者HBV基因表型與血清學(xué)測定的臨床意義
    精品久久久久久久毛片微露脸 | 国产男人的电影天堂91| 9191精品国产免费久久| 晚上一个人看的免费电影| 午夜免费男女啪啪视频观看| 亚洲av日韩在线播放| 99精品久久久久人妻精品| www.熟女人妻精品国产| 9热在线视频观看99| 在线观看人妻少妇| 肉色欧美久久久久久久蜜桃| 十分钟在线观看高清视频www| 久久99热这里只频精品6学生| 午夜两性在线视频| 宅男免费午夜| 日本av免费视频播放| 国产精品免费大片| 9热在线视频观看99| 巨乳人妻的诱惑在线观看| 国产熟女欧美一区二区| 亚洲国产欧美在线一区| 亚洲一区中文字幕在线| 欧美黑人精品巨大| 色精品久久人妻99蜜桃| 国产人伦9x9x在线观看| 一级黄片播放器| avwww免费| 日本a在线网址| 十分钟在线观看高清视频www| av国产久精品久网站免费入址| 嫩草影视91久久| 亚洲国产日韩一区二区| 国产成人av教育| 啦啦啦在线免费观看视频4| 激情视频va一区二区三区| 亚洲精品一二三| 亚洲av国产av综合av卡| 天堂8中文在线网| 一边摸一边做爽爽视频免费| 国产欧美日韩一区二区三 | 国产精品一国产av| 国产99久久九九免费精品| 亚洲精品中文字幕在线视频| 男女床上黄色一级片免费看| 日韩伦理黄色片| 日日夜夜操网爽| 免费在线观看影片大全网站 | 免费久久久久久久精品成人欧美视频| 国产黄色视频一区二区在线观看| 18禁国产床啪视频网站| 最黄视频免费看| 丁香六月欧美| 精品亚洲成a人片在线观看| 波野结衣二区三区在线| 色婷婷av一区二区三区视频| 欧美性长视频在线观看| 在线观看免费视频网站a站| 成年美女黄网站色视频大全免费| 美女扒开内裤让男人捅视频| 五月天丁香电影| 亚洲精品一卡2卡三卡4卡5卡 | 久久人人97超碰香蕉20202| 9191精品国产免费久久| 日韩av在线免费看完整版不卡| 人人妻人人澡人人看| 午夜久久久在线观看| av国产精品久久久久影院| 好男人视频免费观看在线| 黄色视频在线播放观看不卡| 老司机靠b影院| 狠狠精品人妻久久久久久综合| 久久人人97超碰香蕉20202| 日本91视频免费播放| 赤兔流量卡办理| 国产免费现黄频在线看| 国产成人精品无人区| 欧美人与性动交α欧美精品济南到| av又黄又爽大尺度在线免费看| 欧美日韩视频高清一区二区三区二| 国产在视频线精品| 精品国产超薄肉色丝袜足j| 免费黄频网站在线观看国产| 欧美老熟妇乱子伦牲交| 成在线人永久免费视频| 久久性视频一级片| 国产97色在线日韩免费| 又黄又粗又硬又大视频| 日本欧美视频一区| 国产精品偷伦视频观看了| 人人妻人人添人人爽欧美一区卜| 亚洲图色成人| 99国产精品一区二区蜜桃av | 波多野结衣av一区二区av| 精品国产乱码久久久久久男人| 18禁裸乳无遮挡动漫免费视频| 日本欧美国产在线视频| 午夜福利乱码中文字幕| 一级,二级,三级黄色视频| 亚洲国产精品一区三区| 午夜免费成人在线视频| 天天躁夜夜躁狠狠久久av| 国产片内射在线| 日韩中文字幕欧美一区二区 | 一区二区三区乱码不卡18| 搡老岳熟女国产| 纵有疾风起免费观看全集完整版| 50天的宝宝边吃奶边哭怎么回事| 少妇 在线观看| 国产一区有黄有色的免费视频| 搡老岳熟女国产| av电影中文网址| 日本欧美国产在线视频| 久久久久久免费高清国产稀缺| 国产欧美日韩一区二区三区在线| 一区福利在线观看| 国产一区二区激情短视频 | 久久九九热精品免费| 一级毛片 在线播放| 日本av手机在线免费观看| 一本大道久久a久久精品| 欧美日韩亚洲高清精品| 欧美黑人精品巨大| 国产男女超爽视频在线观看| 一区二区三区精品91| 成年女人毛片免费观看观看9 | 丰满迷人的少妇在线观看| 五月天丁香电影| 精品人妻在线不人妻| 18在线观看网站| 黑人欧美特级aaaaaa片| 在线观看人妻少妇| 日韩av在线免费看完整版不卡| 精品一区二区三区av网在线观看 | 欧美乱码精品一区二区三区| 在线观看人妻少妇| 亚洲国产欧美一区二区综合| 久久99热这里只频精品6学生| 亚洲av在线观看美女高潮| 性色av乱码一区二区三区2| 免费av中文字幕在线| 啦啦啦啦在线视频资源| 国产av精品麻豆| 亚洲精品国产一区二区精华液| 精品少妇一区二区三区视频日本电影| 1024香蕉在线观看| 首页视频小说图片口味搜索 | 69精品国产乱码久久久| 国产1区2区3区精品| 国产精品一二三区在线看| 高清欧美精品videossex| 久久精品久久精品一区二区三区| 蜜桃国产av成人99| 成年动漫av网址| 最近最新中文字幕大全免费视频 | 亚洲中文字幕日韩| 成人午夜精彩视频在线观看| 高清黄色对白视频在线免费看| 成人国语在线视频| 国产免费视频播放在线视频| 亚洲图色成人| 国产有黄有色有爽视频| 18禁裸乳无遮挡动漫免费视频| 超碰成人久久| 精品亚洲成国产av| 在线看a的网站| 国产男人的电影天堂91| 王馨瑶露胸无遮挡在线观看| 亚洲国产毛片av蜜桃av| 亚洲av日韩在线播放| 啦啦啦在线观看免费高清www| 亚洲欧美中文字幕日韩二区| 精品国产一区二区三区四区第35| 国产亚洲精品久久久久5区| 只有这里有精品99| 黄色a级毛片大全视频| 亚洲三区欧美一区| 国产成人一区二区在线| 亚洲精品日韩在线中文字幕| 爱豆传媒免费全集在线观看| 99九九在线精品视频| 欧美成人精品欧美一级黄| 巨乳人妻的诱惑在线观看| 国产av精品麻豆| 最黄视频免费看| 青青草视频在线视频观看| 国产精品免费视频内射| 久久国产精品男人的天堂亚洲| 天天添夜夜摸| 精品国产乱码久久久久久小说| 久久人人爽av亚洲精品天堂| 久热爱精品视频在线9| 美女中出高潮动态图| 亚洲七黄色美女视频| 桃花免费在线播放| 91九色精品人成在线观看| 丝瓜视频免费看黄片| 视频区欧美日本亚洲| 国产麻豆69| 亚洲欧美精品综合一区二区三区| 19禁男女啪啪无遮挡网站| 成年女人毛片免费观看观看9 | 男女之事视频高清在线观看 | 精品视频人人做人人爽| 国产成人91sexporn| 亚洲午夜精品一区,二区,三区| 午夜视频精品福利| 国产成人精品久久二区二区免费| 下体分泌物呈黄色| 久久狼人影院| 夫妻性生交免费视频一级片| 视频在线观看一区二区三区| 国产亚洲欧美精品永久| 无限看片的www在线观看| 成人国产一区最新在线观看 | 美女扒开内裤让男人捅视频| 欧美人与性动交α欧美软件| 欧美+亚洲+日韩+国产| 久久精品人人爽人人爽视色| 99热国产这里只有精品6| av在线播放精品| 老汉色∧v一级毛片| 91精品三级在线观看| 亚洲精品久久久久久婷婷小说| 青青草视频在线视频观看| 欧美成人精品欧美一级黄| av有码第一页| 久久ye,这里只有精品| 好男人电影高清在线观看| 老鸭窝网址在线观看| 美女福利国产在线| 欧美人与善性xxx| 国产男人的电影天堂91| 涩涩av久久男人的天堂| 免费在线观看黄色视频的| 黄色视频在线播放观看不卡| 黄频高清免费视频| 99久久综合免费| 另类亚洲欧美激情| 又粗又硬又长又爽又黄的视频| 欧美成人精品欧美一级黄| 亚洲熟女精品中文字幕| 国产亚洲av高清不卡| 成在线人永久免费视频| 黄色毛片三级朝国网站| 亚洲,欧美,日韩| 国产成人欧美在线观看 | 18禁裸乳无遮挡动漫免费视频| 亚洲av综合色区一区| 久久99一区二区三区| 中文字幕av电影在线播放| 国产精品国产av在线观看| 日韩大片免费观看网站| 亚洲av在线观看美女高潮| 亚洲五月婷婷丁香| 亚洲精品成人av观看孕妇| 国产精品av久久久久免费| 色视频在线一区二区三区| 日韩精品免费视频一区二区三区| 一区二区av电影网| 亚洲精品国产av蜜桃| 最新的欧美精品一区二区| 又黄又粗又硬又大视频| 99精国产麻豆久久婷婷| 精品少妇内射三级| 日本欧美视频一区| 亚洲国产精品一区二区三区在线| 国产午夜精品一二区理论片| 亚洲黑人精品在线| 在线av久久热| 又黄又粗又硬又大视频| 91麻豆精品激情在线观看国产 | 欧美乱码精品一区二区三区| 午夜福利视频在线观看免费| 亚洲av日韩在线播放| 国产激情久久老熟女| 狂野欧美激情性bbbbbb| 精品少妇黑人巨大在线播放| 亚洲av男天堂| 中文字幕精品免费在线观看视频| 久久国产精品男人的天堂亚洲| 母亲3免费完整高清在线观看| 亚洲av电影在线观看一区二区三区| 亚洲欧美激情在线| 欧美黄色片欧美黄色片| 大香蕉久久成人网| 久久国产亚洲av麻豆专区| 国产深夜福利视频在线观看| 国产成人精品久久二区二区91| 精品一区在线观看国产| 久久久久久久国产电影| 亚洲一区二区三区欧美精品| 国产精品99久久99久久久不卡| 香蕉丝袜av| 高清不卡的av网站| 国产精品二区激情视频| 国产日韩一区二区三区精品不卡| 宅男免费午夜| 国产精品免费大片| 国产淫语在线视频| 国产男女内射视频| 男人添女人高潮全过程视频| 又粗又硬又长又爽又黄的视频| 女人高潮潮喷娇喘18禁视频| 精品一区二区三区四区五区乱码 | 永久免费av网站大全| 亚洲国产最新在线播放| 90打野战视频偷拍视频| 啦啦啦在线免费观看视频4| 亚洲,欧美精品.| 欧美国产精品一级二级三级| 国产精品香港三级国产av潘金莲 | 久久久欧美国产精品| 成年女人毛片免费观看观看9 | 国产国语露脸激情在线看| 天堂俺去俺来也www色官网| 亚洲av欧美aⅴ国产| 国语对白做爰xxxⅹ性视频网站| 国产一区二区激情短视频 | 色综合欧美亚洲国产小说| 国产成人av教育| 午夜福利视频精品| 午夜福利视频在线观看免费| 日韩av不卡免费在线播放| 黑人猛操日本美女一级片| 啦啦啦啦在线视频资源| 亚洲av电影在线进入| 2018国产大陆天天弄谢| 一边摸一边做爽爽视频免费| 波野结衣二区三区在线| 免费av中文字幕在线| 国产日韩欧美视频二区| 看十八女毛片水多多多| 亚洲国产av新网站| 欧美 亚洲 国产 日韩一| 国产精品香港三级国产av潘金莲 | 黄色怎么调成土黄色| 久久精品国产a三级三级三级| 亚洲欧美精品自产自拍| 在线观看国产h片| 亚洲av在线观看美女高潮| 母亲3免费完整高清在线观看| 人人妻人人添人人爽欧美一区卜| 久久久国产欧美日韩av| 天天躁夜夜躁狠狠久久av| 99国产精品99久久久久| 精品视频人人做人人爽| 午夜av观看不卡| 日本欧美国产在线视频| 欧美日本中文国产一区发布| 美女扒开内裤让男人捅视频| www.自偷自拍.com| 久久精品成人免费网站| 大片免费播放器 马上看| 91老司机精品| 在线 av 中文字幕| 久久av网站| 少妇的丰满在线观看| 亚洲中文字幕日韩| 午夜福利免费观看在线| 黑人猛操日本美女一级片| 国产欧美日韩一区二区三区在线| 亚洲国产精品一区二区三区在线| 久久综合国产亚洲精品| 国产日韩欧美在线精品| 国产精品成人在线| 亚洲国产欧美网| 欧美在线黄色| 欧美日韩视频精品一区| 99国产精品一区二区三区| kizo精华| 十分钟在线观看高清视频www| 日韩中文字幕欧美一区二区 | 国产精品一区二区免费欧美 | 国产黄频视频在线观看| 每晚都被弄得嗷嗷叫到高潮| 少妇人妻 视频| 亚洲av在线观看美女高潮| 亚洲伊人久久精品综合| 少妇精品久久久久久久| 2018国产大陆天天弄谢| 1024香蕉在线观看| 免费在线观看日本一区| 久久久国产精品麻豆| 王馨瑶露胸无遮挡在线观看| 国产伦人伦偷精品视频| 国产精品久久久久久精品电影小说| 久久久国产精品麻豆| 欧美久久黑人一区二区| 在线天堂中文资源库| 国产1区2区3区精品| 色94色欧美一区二区| 精品卡一卡二卡四卡免费| 国产亚洲欧美在线一区二区| 天天躁日日躁夜夜躁夜夜| videosex国产| 午夜福利影视在线免费观看| 男女下面插进去视频免费观看| 欧美日韩av久久| 天天躁夜夜躁狠狠久久av| 亚洲精品成人av观看孕妇| 大香蕉久久成人网| 天天躁夜夜躁狠狠躁躁| 国产精品亚洲av一区麻豆| 欧美在线黄色| 国产高清国产精品国产三级| 欧美日本中文国产一区发布| 欧美日韩国产mv在线观看视频| 最近中文字幕2019免费版| av电影中文网址| 亚洲av日韩在线播放| 97人妻天天添夜夜摸| 丝袜人妻中文字幕| 久久久久久久大尺度免费视频| 国产成人av教育| 91成人精品电影| 成人三级做爰电影| 欧美老熟妇乱子伦牲交| 啦啦啦视频在线资源免费观看| 夫妻性生交免费视频一级片| 在线 av 中文字幕| 嫁个100分男人电影在线观看 | 国产97色在线日韩免费| 一区福利在线观看| 精品一区二区三区av网在线观看 | av天堂在线播放| 狂野欧美激情性xxxx| h视频一区二区三区| www.熟女人妻精品国产| 欧美人与善性xxx| 国产免费视频播放在线视频| 精品国产一区二区三区久久久樱花| 日本欧美国产在线视频| 汤姆久久久久久久影院中文字幕| 每晚都被弄得嗷嗷叫到高潮| 最新在线观看一区二区三区 | 人人妻人人澡人人爽人人夜夜| 国产成人精品久久二区二区免费| 亚洲欧美一区二区三区黑人| 亚洲欧美日韩另类电影网站| 一区在线观看完整版| 亚洲av在线观看美女高潮| 久久久久国产一级毛片高清牌| h视频一区二区三区| 久久热在线av| 波野结衣二区三区在线| 国产精品久久久久久精品电影小说| 国产成人影院久久av| 美女视频免费永久观看网站| 一级毛片 在线播放| 美女中出高潮动态图| 91麻豆av在线| 80岁老熟妇乱子伦牲交| 久久精品久久精品一区二区三区| 国产av精品麻豆| 手机成人av网站| 又大又爽又粗| www.自偷自拍.com| 一边亲一边摸免费视频| 精品少妇黑人巨大在线播放| 女人爽到高潮嗷嗷叫在线视频| 纵有疾风起免费观看全集完整版| 妹子高潮喷水视频| 菩萨蛮人人尽说江南好唐韦庄| 欧美97在线视频| 成人18禁高潮啪啪吃奶动态图| 国产三级黄色录像| 久久精品亚洲熟妇少妇任你| 亚洲av国产av综合av卡| 自线自在国产av| 日韩中文字幕视频在线看片| 黄片小视频在线播放| 国产成人a∨麻豆精品| 国产免费福利视频在线观看| 在线观看www视频免费| 久久这里只有精品19| 日韩电影二区| av国产精品久久久久影院| 欧美激情高清一区二区三区| 国产日韩欧美视频二区| 汤姆久久久久久久影院中文字幕| 亚洲欧美激情在线| 超色免费av| 免费不卡黄色视频| 国产精品一二三区在线看| 中国美女看黄片| 免费日韩欧美在线观看| 男女之事视频高清在线观看 | 亚洲专区国产一区二区| 国产成人精品久久二区二区91| 中文字幕精品免费在线观看视频| 90打野战视频偷拍视频| 日日摸夜夜添夜夜爱| 亚洲国产日韩一区二区| 亚洲精品久久午夜乱码| 在线观看免费视频网站a站| 熟女少妇亚洲综合色aaa.| 日韩 欧美 亚洲 中文字幕| 一级毛片 在线播放| 精品人妻一区二区三区麻豆| 又大又黄又爽视频免费| 午夜精品国产一区二区电影| 黄色片一级片一级黄色片| 亚洲专区中文字幕在线| 女警被强在线播放| 自线自在国产av| 精品久久久久久电影网| 十八禁高潮呻吟视频| 99热网站在线观看| 久久久国产欧美日韩av| 一本色道久久久久久精品综合| 免费黄频网站在线观看国产| 操美女的视频在线观看| 亚洲专区中文字幕在线| 观看av在线不卡| 18禁观看日本| 极品人妻少妇av视频| 久久精品aⅴ一区二区三区四区| av在线播放精品| 在线看a的网站| 国产成人欧美| 国产精品人妻久久久影院| 亚洲成人国产一区在线观看 | 国产亚洲欧美在线一区二区| 丰满迷人的少妇在线观看| 久久国产精品男人的天堂亚洲| 精品人妻熟女毛片av久久网站| 亚洲情色 制服丝袜| 午夜免费观看性视频| 激情五月婷婷亚洲| 久久亚洲精品不卡| av网站在线播放免费| 美女中出高潮动态图| www.精华液| 成人午夜精彩视频在线观看| 亚洲精品国产色婷婷电影| 蜜桃在线观看..| 青春草亚洲视频在线观看| 欧美老熟妇乱子伦牲交| 日本av手机在线免费观看| 少妇精品久久久久久久| 中文字幕人妻熟女乱码| 久久人人97超碰香蕉20202| 女警被强在线播放| 黄网站色视频无遮挡免费观看| 纵有疾风起免费观看全集完整版| 777久久人妻少妇嫩草av网站| 国产免费福利视频在线观看| 女性被躁到高潮视频| 交换朋友夫妻互换小说| 国产精品一区二区精品视频观看| 免费看十八禁软件| 国产精品国产三级国产专区5o| 性少妇av在线| 一本一本久久a久久精品综合妖精| 夜夜骑夜夜射夜夜干| 亚洲精品久久久久久婷婷小说| 女人被躁到高潮嗷嗷叫费观| 国产亚洲午夜精品一区二区久久| 夫妻性生交免费视频一级片| 丰满饥渴人妻一区二区三| 国产97色在线日韩免费| 99re6热这里在线精品视频| 亚洲精品日本国产第一区| 久久精品久久久久久噜噜老黄| 午夜福利一区二区在线看| 丝袜美足系列| 免费高清在线观看视频在线观看| av在线老鸭窝| 亚洲,一卡二卡三卡| 日韩一本色道免费dvd| 国产日韩欧美亚洲二区| 亚洲欧美日韩高清在线视频 | 男女边摸边吃奶| 啦啦啦视频在线资源免费观看| 宅男免费午夜| 国产一区二区 视频在线| 久久精品国产亚洲av涩爱| 欧美精品av麻豆av| 国产精品二区激情视频| 国产1区2区3区精品| 欧美xxⅹ黑人| 九色亚洲精品在线播放| 欧美精品一区二区大全| 亚洲欧美日韩另类电影网站| 久9热在线精品视频| 在线观看人妻少妇| 国产成人啪精品午夜网站| 国产欧美日韩精品亚洲av| 人妻人人澡人人爽人人| 国产精品 国内视频| 下体分泌物呈黄色| 日本色播在线视频| 国产爽快片一区二区三区| 国产成人av教育| 女人高潮潮喷娇喘18禁视频| av网站在线播放免费| 欧美人与善性xxx| 国产免费视频播放在线视频| 色网站视频免费| 水蜜桃什么品种好| 别揉我奶头~嗯~啊~动态视频 | 热99国产精品久久久久久7| 久久人妻熟女aⅴ| 美女国产高潮福利片在线看| 大香蕉久久成人网| 久久精品国产a三级三级三级| 国产精品人妻久久久影院| 热99久久久久精品小说推荐| 真人做人爱边吃奶动态| 999精品在线视频| 日韩人妻精品一区2区三区| 国产一区二区 视频在线| videosex国产| 成年人午夜在线观看视频| 视频在线观看一区二区三区|