亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        普通小麥主要農(nóng)藝性狀的全基因組關聯(lián)分析

        2019-09-10 02:31:02翟俊鵬李海霞畢惠惠周思遠羅肖艷陳樹林程西永許海霞
        作物學報 2019年10期
        關鍵詞:農(nóng)藝染色體基因組

        翟俊鵬 李海霞 畢惠惠 周思遠 羅肖艷 陳樹林 程西永 許海霞

        普通小麥主要農(nóng)藝性狀的全基因組關聯(lián)分析

        翟俊鵬 李海霞 畢惠惠 周思遠 羅肖艷 陳樹林 程西永*許海霞*

        河南農(nóng)業(yè)大學 / 省部共建小麥玉米作物學國家重點實驗室 / 河南省糧食作物協(xié)同創(chuàng)新中心, 河南鄭州 450046

        為解析小麥復雜農(nóng)藝性狀的遺傳機制, 本研究以150份小麥品種(系)為自然群體, 在4個環(huán)境條件下測定了9個主要農(nóng)藝性狀, 利用小麥35K SNP芯片, 結合5種關聯(lián)模型(Q、PCA、K、PCA+K、Q+K), 進行全基因組關聯(lián)分析。結果表明, 全基因組多態(tài)性信息量PIC的范圍為0.0950~0.5000, 最小等位基因頻率MAF值為0.0500~0.5000; 群體結構分析和PCA分析均表明參試材料可分為兩個亞群; 連鎖不平衡分析發(fā)現(xiàn)A基因組、B基因組、D基因組和全基因組的LD衰減距離分別為4.7、8、11和6 Mb。9個性狀共檢測到652個顯著的關聯(lián)位點(≤0.001), 其中21個SNP在2個或2個以上的環(huán)境中被重復檢測到, 分布在1A(1)、1B(4)、2A(3)、2D(2)、3A(1)、5A(1)、5B(5)、6A(1)、6B(2)和7D(3)染色體上; 1個SNP標記的物理位置未知, 3個SNP標記同時與2個性狀顯著關聯(lián); 單個SNP 的表型貢獻率為7.67%~18.79%。8個優(yōu)勢等位變異在供試群體中所占比例較低, 篩選出14個可能與小麥農(nóng)藝性狀相關的候選基因, 其中、和可能在植物抵御生物與非生物脅迫中起作用,和可能與植物激素的合成和響應有關,可能與植物細胞壁的增強有關,可能參與葉綠體發(fā)育, 另外7個候選基因的功能未知。

        小麥; 農(nóng)藝性狀; 連鎖不平衡; 全基因組關聯(lián)分析; 候選基因

        小麥是中國主要糧食作物之一, 其產(chǎn)量占全國糧食總產(chǎn)量的22%左右, 小麥生產(chǎn)對中國糧食安全具有重要意義。小麥產(chǎn)量構成因子包括穗數(shù)、穗粒數(shù)和千粒重[1-2], 此外, 株高[3]、穗下節(jié)長和穗長[4]、粒寬[5]、可育小穗數(shù)[6]等其他農(nóng)藝性狀對小麥產(chǎn)量形成也有重要作用。育種工作者經(jīng)過長期的品種改良, 培育出了一大批廣適多抗的高產(chǎn)品種[7], 但傳統(tǒng)育種方法盲目性大, 效率較低, 導致品種同質(zhì)化嚴重[8]。分子標記輔助選擇(marker-assisted selection, MAS)有助于提高選擇效率和可預測性, 進而加速育種進程[9]。

        小麥多數(shù)性狀是多基因控制的數(shù)量性狀[10-11], 受基因和環(huán)境的共同調(diào)控。QTL定位和全基因組關聯(lián)分析(genome wide association study, GWAS)是研究數(shù)量性狀的主要方法, Ramya等[12]利用復合區(qū)間作圖法, 檢測到25個與小麥籽粒性狀有關的QTL, 單個位點的表型貢獻率為4.15%~15.53%。Ren等[13]在兩種灌溉條件下對小麥株高性狀進行研究, 共發(fā)現(xiàn)17個QTL, 其中12個在灌溉條件下被檢測到, 3個在干旱條件下被檢測到, 僅有2個在2種環(huán)境中被重復檢測到, 說明基因的表達受環(huán)境條件的影響較大。張坤普等[14]在5D染色體上和區(qū)間發(fā)現(xiàn)與穗粒數(shù)、總小穗數(shù)、小穗著生密度和可育小穗數(shù)相關的QTL, 遺傳效應值分別為11.67%、13.83%、12.26%和10.22%。GWAS亦稱關聯(lián)作圖(association mapping), 是以連鎖不平衡(linkage disequilibrium, LD)為基礎, 以自然群體為研究對象, 用于鑒定基因和定位基因的分析方法, 現(xiàn)已廣泛應用于小麥[15]、玉米[16]、大豆[17]、花生[18]、水稻[19]等農(nóng)作物。傳統(tǒng)的QTL定位精度較低, 只能將目標基因定位在5~20 cM之間, 相較于雙親連鎖分析, 全基因組關聯(lián)分析考慮了更加復雜的遺傳背景[20], 在解析小麥數(shù)量性狀遺傳機制方面的應用越來越廣泛[21]。朱玉磊等[22]利用181對SSR標記與264份小麥材料進行穗發(fā)芽抗性的 GWAS 研究, 發(fā)現(xiàn)8個重復關聯(lián)的標記。Wang等[23]對105份小麥材料進行GWAS分析, 共檢測到6個與穗數(shù)相關聯(lián)的位點, 分布在1D(1)、3B(2)、4A(1)和4B(2)染色體上, 表型貢獻率為12.18%~16.01%。Mwadzingeni等[24]以93個普通小麥為研究對象, 在干旱脅迫和非脅迫的環(huán)境條件下對部分農(nóng)藝性狀進行GWAS分析,檢測到20個與穗長相關聯(lián)的標記, 9個與穗粒數(shù)顯著關聯(lián)的標記。Liu等[25]運用GLM和MLM兩種模型對小麥23個農(nóng)藝性狀進行GWAS, 發(fā)現(xiàn)了25個與重要農(nóng)藝性狀相關的候選基因。

        目前利用SSR標記對小麥農(nóng)藝性狀進行關聯(lián)分析的研究較多, 與SSR標記相比, SNP標記在基因組中分布更廣, 數(shù)量更多, 定位精度更高, 是目前最具潛力的分子標記[26]。本研究以150份小麥品種(系)構成的自然群體為材料, 利用均勻分布于小麥21條染色體上的35,542個SNP標記, 對穗下節(jié)長、穗粒數(shù)、千粒重、株高、穗長、可育小穗數(shù)、抽穗期、每平方米穗數(shù)、粒寬共9個農(nóng)藝性狀進行全基因組關聯(lián)分析, 篩選優(yōu)勢等位變異, 發(fā)掘相關候選基因, 為基因克隆和分子標記輔助選擇提供理論依據(jù)。

        1 材料與方法

        1.1 試驗材料與表型鑒定

        選取150份小麥品種(系)作為自然群體(表1), 其中含黃淮南片國家小麥區(qū)域試驗品種(系) 35份, 河南省小麥區(qū)域試驗品種(系) 13份, 河南省小麥聯(lián)合體品種比較試驗品種(系) 30份, 本課題組選育出的綜合性狀表現(xiàn)優(yōu)異的高代品系72份, 該群體一定程度上反映了黃淮南片麥區(qū)高產(chǎn)育種目標。2016—2017年度將參試材料分別種植于河南鹿邑華冠種業(yè)選種基地(E1: 17-鹿邑)、河南農(nóng)業(yè)大學鄭州科教園區(qū)(E2: 17-鄭州); 2017—2018年度分別種植于河南農(nóng)業(yè)大學原陽科教園區(qū)(E3: 18-原陽)、河南農(nóng)業(yè)大學鄭州科教園區(qū)(E4: 18-鄭州)。采用隨機區(qū)組設計, 每份材料種植6行, 行長2 m, 行距20 cm, 2次重復, 基本苗18~22萬株 hm–2, 機器條播。在試驗田四周設置保護行, 試驗田土壤肥力、灌溉等田間管理基本一致。

        試驗小區(qū)中50%以上個體抽出1/3穗(不含芒)時定為抽穗期, 灌漿后期從每個小區(qū)隨機選取10株測量穗下節(jié)長、株高、穗長和可育小穗數(shù), 蠟熟后期調(diào)查每平方米穗數(shù), 每個材料取30穗脫粒, 調(diào)查穗粒數(shù)和千粒重, 采用Perten SKCS 4100單粒谷物特性分析儀測定粒寬。

        表1 供試小麥材料

        1.2 表型數(shù)據(jù)的統(tǒng)計分析

        利用SPSS 22.0軟件對表型數(shù)據(jù)進行正態(tài)性檢驗, 利用QTL IciMapping軟件進行描述性統(tǒng)計分析和方差分析[27], 廣義遺傳力2=2g/[2g+(2g/)+ (2e/)], 其中2g為基因型方差,2e為殘差方差,2g為基因與環(huán)境互作方差,為地點數(shù),為重復數(shù)。

        1.3 DNA提取及SNP分型

        小麥生長至三葉期時取10株幼苗, 將葉片混合后利用改良的CTAB法提取DNA[28], 用1.0%瓊脂糖凝膠電泳檢驗DNA質(zhì)量和濃度。利用博奧晶典生物技術有限公司與中國科學院遺傳與發(fā)育生物學研究所聯(lián)合開發(fā)的小麥基因組35K SNP芯片對參試材料進行基因分型, 共檢測到覆蓋小麥全基因組的35,143個SNP標記, Allen等[29]經(jīng)適應性評估發(fā)現(xiàn), 這套芯片非常適合于六倍體小麥的基因分型。采用Tassel 5.0軟件對標記信息去冗余, 剔除最小等位基因頻率(minor allele frequency, MAF)小于5 %和缺失率(miss rate)大于25%的SNP, 最終得到9552個多態(tài)性高、穩(wěn)定性好的優(yōu)質(zhì)SNP標記用于后續(xù)關聯(lián)分析。多態(tài)性信息量(polymorphism information content): PIC = 1?∑2(其中P表示第位點的第個等位變異出現(xiàn)的頻率)[30]。

        1.4 群體結構與主成分分析

        利用Structure 2.3.4進行群體結構分析, 將最佳群體數(shù)目設置為1~12, 每個值運算5次, 開始時的不作數(shù)迭代(length of burnin period)設置為10,000, 不作數(shù)迭代后的MCMC (number of MCMC reps after burnin)設置為100,000, 依據(jù)Evanno等[31]的方法選取最大Δ對應的值作為最佳亞群數(shù)目。利用Tassel 5.0 軟件對參試材料進行基于SNP芯片數(shù)據(jù)的主成分分析。

        1.5 連鎖不平衡分析

        利用Tassel 5.0軟件進行連鎖不平衡分析, 設置LD衰減窗口大小為100, 即計算每個標記前后100對標記之間的連鎖不平衡程度, 當標記間≤0.001時, 認為存在顯著的LD[32]。

        1.6 關聯(lián)分析與候選基因篩選

        利用GLM (general linear model)一般線性模型(Q、PCA)和MLM (mixed linear model)混合線性模型(K、PCA+K、Q+K)對各性狀進行模型修正, 采用RStudio 1.1.453 (https://www.rstudio.com/)軟件包lme4中l(wèi)mer函數(shù)[33]計算每個性狀的最佳線性無偏預測值(best linear unbiased prediction, BLUP), 去除環(huán)境效應, 依據(jù) BLUP值的關聯(lián)結果作出5種關聯(lián)模型的Quantile-Quantile圖(Q-Q Plot), 考慮到實際?lg ()值與預期?lg ()值的偏差, 以實際?lg ()值最接近?lg ()期望值的模型作為性狀的最優(yōu)模型[34],基于最優(yōu)模型分別對各性狀在4個環(huán)境中的表型值和BLUP值進行GWAS, 以BLUP值的關聯(lián)結果進行各性狀的曼哈頓圖(Manhattan Plot)作圖。

        由RStudio 1.1.453軟件完成Q-Q Plot、Manhattan Plot作圖。

        以LD距離作為候選基因的預測區(qū)間, 將多個環(huán)境中穩(wěn)定出現(xiàn)的SNP標記的延伸序列在小麥TGACv1.0數(shù)據(jù)庫(http://plants.ensembl.org/hmmer/ index.html)和NCBI網(wǎng)站(https://www.ncbi.nlm.nih. gov/)進行Blast比對, 對候選位點進行功能注釋。

        2 結果與分析

        2.1 農(nóng)藝性狀表型數(shù)據(jù)的統(tǒng)計分析

        從表2可以看出, 9個農(nóng)藝性狀的廣義遺傳力均高于60%, 說明這些性狀的表型主要受基因調(diào)控, 受環(huán)境影響較小。方差分析顯示所有性狀在基因型間和亞群間均存在顯著或極顯著差異, 表明研究的性狀在該群體中具有豐富的遺傳變異, 但不同性狀的變異程度具有明顯差異, 其中每平方米穗數(shù)(20.09%)和穗粒數(shù)(18.25%)的變異系數(shù)較大, 變異程度較高; 抽穗期(3.29%)和粒寬(4.61%)的變異系數(shù)較小, 變異程度較低。正態(tài)性檢驗結果表明, 除抽穗期外, 其他性狀的峰度和偏度均接近0, 基本符合正態(tài)分布, 說明目標性狀大部分是多基因控制的數(shù)量性狀。

        2.2 SNP標記分析

        全基因組分析中PIC的變異范圍為0.0950~ 0.5000, 平均值為0.3215, 其中3B染色體最高(0.3729), 4D染色體最低(0.2590), 3個染色體組中PIC表現(xiàn)為A (0.3317)>B (0.3240)>D (0.3100) (表3)。95.38% (9111/9552)的SNP標記具有物理位置信息, A、B、D染色體組分別檢測到2747、3708和2656個SNP標記, 分別占總標記數(shù)的28.76%、38.82%和27.80%。全基因組最小等位基因頻率 MAF 平均值為0.2329, 變異范圍為0.0500~0.5000, 49.11%標記的MAF大于0.2, 16.25% 標記的MAF大于0.4, 物理圖譜長度為14,541.4 Mb, 標記密度為1.70 Mb/marker。

        2.3 群體結構分析和主成分分析

        利用Structure 2.3.4軟件對群體遺傳結構進行分析, 當=2時, Δ值最大(圖1-A), 由圖1-B可知, 參試材料被劃分為2個亞群(sub-population), 其中亞群I (Sub-Pop I)包含39份材料, 亞群II (Sub-Pop II)包含111份材料; 基于基因型的PCA分析顯示整個群體可分為2個群組, 分類結果與Structure 結果高度吻合(圖1-C)。方差分析發(fā)現(xiàn), 亞群I中品種(系)的穗下節(jié)長、穗粒數(shù)、千粒重、穗長、可育小穗數(shù)、粒寬和抽穗期均顯著或極顯著低于亞群II, 而株高和每平方米穗數(shù)極顯著高于亞群II (圖1-D和表2), 說明2個亞群在表型性狀和遺傳背景上均有較大差異。

        2.4 連鎖不平衡分析

        考慮連鎖群(linkage group)間連鎖不平衡背景的影響, 對連鎖群間的2值進行平方根轉換, 以大于此分布95%的參數(shù)值為閾值, 用來截取同一連鎖群內(nèi)LD的衰減距離[35], A基因組、B基因組、D基因組和全基因組的2閾值分別為0.2846、0.3162、0.2763和0.2817, 從連鎖不平衡的衰減圖可以看出, A基因組、B基因組、D基因組和全基因組的LD距離分別為4.7、8、11和6 Mb (圖2)。依據(jù)全基因組的LD衰減距離, 將多環(huán)境中穩(wěn)定出現(xiàn)的SNP標記在物理圖譜上前后6 Mb的區(qū)間認定為一個候選位點, 推測可能有控制相關性狀的QTL存在。

        表3 SNP標記信息

        2.5 最優(yōu)模型選擇與全基因組關聯(lián)分析

        以Q-Q plot中實際?lg ()值最接近?lg ()期望值為判斷標準, 確定各性狀的最優(yōu)關聯(lián)模型, 從圖3可以看出, 穗粒數(shù)、株高、可育小穗數(shù)和粒寬的最優(yōu)關聯(lián)模型是Q+K模型, 千粒重和抽穗期的最優(yōu)關聯(lián)模型為PCA+K模型, 穗下節(jié)長、穗長和每平方穗數(shù)的最優(yōu)關聯(lián)模型為K模型。

        針對不同性狀選用其最優(yōu)關聯(lián)模型進行GWAS分析, 共檢測到652個顯著SNP位點(≤0.001), 其中21個SNP在2個或2個以上的環(huán)境中被檢測到, 屬于穩(wěn)定的SNP標記(表4)。檢測到與穗下節(jié)長穩(wěn)定關聯(lián)的標記4個, 分別位于2D(1)、6A(1)和6B(2)染色體上, 其中和在3個環(huán)境中被重復檢測到,在2個環(huán)境中被重復檢測到, 可解釋7.68%~10.51%的表型變異。檢測到與穗粒數(shù)穩(wěn)定關聯(lián)的標記 4個, 位于5B(3)和7D(1)染色體上, 均在2個環(huán)境中被重復檢測到,2范圍為7.67%~10.97%。7D染色體上的和1A染色體上的分別與千粒重和株高穩(wěn)定關聯(lián), 均能在3個環(huán)境中穩(wěn)定遺傳, 分別解釋9.16%~10.84%、8.18%~12.41%的表型變異。檢測到與穗長穩(wěn)定關聯(lián)的標記2個, 位于2D(1)和3A(1)染色體上, 均在2個環(huán)境中被重復檢測到,2范圍為9.89%~14.27%。檢測到與可育小穗數(shù)相關聯(lián)的穩(wěn)定標記2個, 位于5B染色體上, 均在3個環(huán)境中被重復檢測到,2范圍為7.69%~11.50%。檢測到與抽穗期相關聯(lián)的穩(wěn)定標記10個, 分別位于1B(4)、2A(3)、5A(1)和7D(1)染色體上,和在2個環(huán)境中被重復檢測到,和在3個環(huán)境中被重復檢測到, 可解釋8.07%~18.79%的表型變異, 其中在物理圖譜上沒有位置信息。沒有檢測到與每平方米穗數(shù)和粒寬穩(wěn)定關聯(lián)的SNP。

        圖1 群體結構、主成分分析和亞群表型

        1-A:D折線圖; 1-B: 群體結構; 1-C: PCA分析; 1-D: 亞群表型; 各性狀詳細名稱見表2。

        1-A: line chart ofD; 1-B: structure plot; 1-C: analysis of PCA; 1-D: sub-population phenotype. Name of each trait is given in Table 2.

        圖2 連鎖不平衡衰減圖

        21個穩(wěn)定關聯(lián)的位點中, 位于5B上的和同時與可育小穗數(shù)和穗粒數(shù)相關聯(lián), 位于2D上的同時與穗長和穗下節(jié)長相關聯(lián), 說明這3個位點為一因多效位點。

        各性狀詳細名稱見表 2; 在2個或2個以上環(huán)境中被重復檢測到的SNP用綠色高亮顯示。

        Name of each trait is given in Table 2; A SNP that detected in two or more environments with green highlighted.

        表4 小麥農(nóng)藝性狀穩(wěn)定關聯(lián)的位點信息

        (續(xù)表4)

        各性狀詳細名稱見表2。Name of each trait is given in Table 2.

        2.6 穩(wěn)定遺傳的SNP效應分析

        對穩(wěn)定出現(xiàn)的標記進行單標記效應分析, 共發(fā)現(xiàn)14個優(yōu)異等位變異(表5), 其中可使穗下節(jié)長度降低1.94 cm, 在群體中出現(xiàn)的頻率占75.84%;、和可使穗粒數(shù)增加5粒左右, 在群體中出現(xiàn)的頻率均在93%以上;和均可使可育小穗數(shù)增加1.25個, 在群體中所占比例為93.24%和93.29%。

        檢測到8個優(yōu)異等位變異在供試群體中所占比例較低, 其中和分別使穗下節(jié)長度降低2.79、2.58和2.58 cm, 在群體中所占比例分別為8.90%、10.74%和10.74%;和分別使穗長增加0.63 cm和0.69 cm, 在群體中占比24.16%和14.05%; 另外、分別增加穗粒數(shù)3.22個、提高千粒重4.91 g、降低株高3.42 cm, 在群體中出現(xiàn)的頻率分別為11.41%、8.16%和32.39%。

        檢測到與抽穗期相關聯(lián)的等位變異10個, 其中、、、、、、和分別延長抽穗期3.69、3.71、2.64、2.65、3.34、3.19、2.22和2.88 d, 在群體中占比分別為93.29%、93.79%、88.51%、87.92%、93.29%、93.79%、89.58%和93.75%。和分別使抽穗期提前1.88 d和3.95 d, 在群體中出現(xiàn)的頻率為75.00%和94.02%。

        表5 穩(wěn)定關聯(lián)的SNP及其等位變異的表型效應

        (續(xù)表5)

        各性狀詳細名稱見表 2。在供試群體中所占比例較小的優(yōu)異等位變異用粗體和下畫線標注。

        Name of each trait is given in Table 2. The small proportion of excellent alleles in the test population is marked in bold and underlined.

        2.7 候選基因分析

        利用21個穩(wěn)定表達的SNP的延伸序列在 TGACv1.0數(shù)據(jù)庫和NCBI網(wǎng)站中Blast, 共發(fā)現(xiàn)10個候選位點, 14個候選基因, 其中7個候選基因有功能注釋(表6)。其中編碼7TM區(qū)域的鈣依賴通道蛋白(Calcium-dependent channel),編碼NB-ARC結構域的植物抗性蛋白,編碼MATE蛋白,編碼植物生長素響應因子(Auxin response factor),編碼Rab GTPase 11家族蛋白,編碼植物損傷誘導蛋白WI12,編碼一種谷氨酰胺轉移酶II類蛋白(Glutamine amidotransferases class-II (GATase), 另外7個候選基因的功能未知。

        表6 小麥農(nóng)藝性狀候選基因及其功能注釋

        各性狀詳細名稱見表 2。Name of each trait is given in Table 2.

        3 討論

        3.1 模型修正與關聯(lián)分析

        Yu等[36]提出的將親緣關系作為協(xié)變量納入分析的混合線性模型被證明是一種有效的關聯(lián)模型, 即使是在較小的群體中也能表現(xiàn)出顯著的關聯(lián)能力,與一般線性模型相比, 減少了假陽性的概率, 但是盲目使用 Q+K 模型可能會導致群體結構的過度修正, 出現(xiàn)假陰性的錯誤[23,36], 故有必要針對不同的性狀進行不同模型的比較。本研究中穗粒數(shù)、株高、可育小穗數(shù)和粒寬適用Q+K模型, 千粒重和抽穗期適用 PCA+K模型, 穗下節(jié)長、穗長和每平方穗數(shù)適用 K 模型。

        3.2 小麥農(nóng)藝性狀的GWAS研究

        小麥基因組重復序列多, 基因組龐大, 即使厘摩范圍內(nèi)也包含上百萬個堿基對[37], 隨著高通量測序技術的發(fā)展, 高密度SNP分型芯片為小麥復雜農(nóng)藝性狀的研究提供了更多的變異信息。本研究利用35K SNP芯片對小麥9個主要農(nóng)藝性狀進行GWAS分析, 檢測到4個與穗下節(jié)長穩(wěn)定關聯(lián)的SNP, 位于2D(1)、6A(1)和6B(2)染色體上。Turuspekov等[21]和Jaiswal等[38]分別在2D和6B染色體上發(fā)現(xiàn)了與穗下節(jié)長顯著關聯(lián)的位點, 而位于6A染色體上的未見報道, 可能是一個新的位點。Jaiswal等[38]在5B染色體長臂上發(fā)現(xiàn)與穗粒數(shù)關聯(lián)的位點(114 cM), 與本研究檢測到的位點距離較近。Shi等[39]在7D染色體上發(fā)現(xiàn)4個與穗粒數(shù)關聯(lián)的位點(、、、), 其中與本研究發(fā)現(xiàn)的位點在物理圖譜上僅距離3 Mb。Zanke等[40]在7D染色體上發(fā)現(xiàn)2個與千粒重關聯(lián)的位點, 其中與本研究發(fā)現(xiàn)的距離5.9 Mb。在1A染色體長臂上檢測到1個與株高穩(wěn)定關聯(lián)的SNP。Sukumaran等[20,41]和Jaiswal等[38]分別在相似位置檢測到與株高關聯(lián)的位點(116~117 cM, 126 cM)。檢測到與穗長穩(wěn)定關聯(lián)的2個SNP分別位于2D和3A上。Mwadzingeni等[24]在相同染色體上也檢測到穗長的關聯(lián)位點。檢測到2個與可育小穗數(shù)穩(wěn)定關聯(lián)的SNP, 位于5B上。Turuspekov等[21]在5B染色體上也檢測到與可育小穗數(shù)關聯(lián)的位點, 但與本研究發(fā)現(xiàn)的位點距離較遠。檢測到10個與抽穗期穩(wěn)定關聯(lián)的SNP, 其中9個具有物理位置信息, 位于1B(4)、2A(3)、5A(1)和7D(1)染色體上。Ain等[42]、Sukumaran等[20]、Mwadzingeni等[24]、Jaiswal等[38]和Zanke等[43]分別在1B、5A和7D染色體上也檢測到與抽穗期關聯(lián)的位點, 其中Sukumaran等[20]發(fā)現(xiàn)的位點與本研究發(fā)現(xiàn)的距離較近。已知控制小麥光周期的基因位于2AS上, 本研究也在2A染色體短臂上發(fā)現(xiàn)3個與抽穗期穩(wěn)定關聯(lián)的SNP, 其中與相距7.16 Mb,可能與有一定聯(lián)系。

        3.3 等位變異及其應用

        本研究共發(fā)現(xiàn)14個優(yōu)異等位變異, 其中6個在供試群體中所占比例相對較高, 8個在供試群體中所占比例較低, 例如增加穗粒數(shù)的等位變異, 提高千粒重的等位變異, 這些等位變異效應較大, 可能與不利基因緊密連鎖, 也可能是新的優(yōu)異等位變異, 育種工作者可以結合分子標記在早期世代選擇中加以利用。

        發(fā)現(xiàn)10個與抽穗期相關的等位變異, 其中8個延長了抽穗期, 2個推遲了抽穗期。與其他性狀不同, 抽穗期等位變異的應用受地區(qū)環(huán)境的影響, 在適宜春性和弱春性品種種植的地區(qū), 注重對抽穗較早的材料保留選擇, 而在適宜冬性和半冬性品種種植的地區(qū), 應優(yōu)先選擇抽穗較晚的材料。對于春季凍害較重的地區(qū), 抽穗過早的品種更易受到春季凍害的影響, 而春季凍害較輕的地區(qū), 抽穗早的品種在后期籽粒灌漿方面更占優(yōu)勢[44]。

        3.4 候選基因分析

        對21個穩(wěn)定表達的SNP的延伸序列在 TGACv1.0數(shù)據(jù)庫和NCBI網(wǎng)站中Blast, 其中7個有功能注釋, 位于5B染色體上穗粒數(shù)和可育小穗數(shù)的候選基因編碼一種7TM區(qū)域的鈣依賴通道蛋白(calcium-dependent channel)。Pei等[45]發(fā)現(xiàn)植物中Ca2+通道的激活會介導脫落酸的合成, 小麥籽粒灌漿后期, 干物質(zhì)增多, 結合水含量上升, 自由水含量減少, 代謝減弱, 脫落酸開始增多, 推測可能參與植物中脫落酸的合成, 從而影響籽粒的結實和成熟。位于5B染色體上穗粒數(shù)的候選基因編碼NB-ARC結構域的植物抗性蛋白, Van等[46]發(fā)現(xiàn)NB-ARC結構域由NB、ARC1和ARC2三個結構域組成, 其結合狀態(tài)可調(diào)節(jié)植物抗性蛋白的活性, 在植物中參與病原體的識別和先天免疫反應的后續(xù)激活, 位于7D染色體上穗粒數(shù)的候選基因編碼一種MATE蛋白, 基因注釋顯示MATE是植物中多基因家族, 其與植物抗病性和滲透脅迫下鐵離子的穩(wěn)定有關, 2個候選基因與植物的抗性和免疫反應有關, 推測其通過提高小麥抵御脅迫的能力來增加穗粒數(shù)。位于7D染色體上千粒重的候選基因編碼植物生長素響應因子(Auxin response factor), 推測其通過調(diào)節(jié)植物中生長素的響應機制來影響粒重。位于1B染色體上抽穗期的候選基因編碼Rab GTPase 11家族蛋白, 參與膜轉運調(diào)節(jié)。一些RAB11成員在植物特有的發(fā)育過程中發(fā)揮作用, 包括細胞分裂和囊泡運輸[47-48], 但大多數(shù)的RAB11成員的分子和生理功能仍不清楚。Asaoka等[49]在擬南芥中對 RAB11 亞群進行了亞細胞定位和動力學研究, 發(fā)現(xiàn)其家族成員RABA1在擬南芥中介導跨高爾基體和質(zhì)膜之間的運輸, 在鹽脅迫耐受性上起關鍵作用。位于2A染色體上抽穗期的候選基因, 編碼植物損傷誘導蛋白WI12, 基因注釋顯示它優(yōu)先在細胞壁中積累, 研究表明它對植物細胞壁有增強作用[50]。位于5A染色體上抽穗期的候選基因編碼一種谷氨酰胺轉移酶II類蛋白[Glutamine amidotransferases class-II (GATase)], 蘇延萍等[51]通過擬南芥突變體試驗和轉錄本驗證, 發(fā)現(xiàn)編碼此類蛋白的基因GATL12在擬南芥葉綠體發(fā)育中發(fā)揮重要作用, 推測通過參與小麥中葉綠體的發(fā)育來調(diào)節(jié)植株的光感反應, 進而影響小麥的抽穗期。

        4 結論

        共檢測到21個與小麥主要農(nóng)藝性狀穩(wěn)定關聯(lián)的SNP標記, 有3個SNP與多個性狀相關聯(lián), 8個優(yōu)勢等位變異在供試材料中所占比例較低。篩選出10個候選位點, 包含14個可能與小麥農(nóng)藝性狀相關的候選基因, 其中3個可能在植物抵御生物與非生物脅迫過程中起作用, 2個可能與植物激素的合成和響應有關, 1個可能與植物細胞壁的增強有關, 1個可能參與葉綠體的發(fā)育, 另外7個功能未知。

        [1] Deng S, Wu X, Wu Y, Zhou R, Wang H, Jia J, Liu S. Characterization and precise mapping of a QTL increasing spike number with pleiotropic effects in wheat., 2011, 122: 281–289.

        [2] Zhang H, Chen J, Li R, Deng Z, Zhang K, Liu B, Tian J. Conditional QTL mapping of three yield components in common wheat (L.)., 2016, 4: 220–228.

        [3] 丁安明, 崔法, 李君, 趙春華, 王秀芹, 王洪剛. 小麥單株產(chǎn)量與株高的QTL分析. 中國農(nóng)業(yè)科學, 2011, 44: 2857–2867. Ding A M, Cui F, Li J, Zhao C H, Wang X Q, Wang H G. QTL analysis on grain yield per plant and plant height in wheat., 2011, 44: 2857–2867 (in Chinese with English abstract).

        [4] 金立橋. 冬小麥麥穗光合生理特性及其對籽粒灌漿的影響. 山東農(nóng)業(yè)大學博士學位論文, 山東泰安, 2017. Jin L Q. Photosynthetic Characteristics of Ear and Contribution of Ear Photosynthesis to Grain Filling in Winter Wheat. PhD Dissertation of Shandong Agricultural University, Tai,an, Shandong, China, 2017 (in Chinese with English abstract)

        [5] Gan Y, Stobbe E H. Seedling vigor and grain yield of ‘Roblin’ wheat affected by seed size., 1996, 88.

        [6] Rodrigues O, Teixeira M C C, Costenaro E R, Damo R. Influence of vernalization and photoperiod on the duration of stem elongation and spikelet fertility in wheat., 2014, 05: 1547–1557.

        [7] Dubcovsky J, Dvorak J. Genome plasticity a key factor in the success of polyploid wheat under domestication., 2007, 316: 1862–1866

        [8] Holland J B. Genetic architecture of complex traits in plants., 2007, 10: 156–161.

        [9] Collard B C Y, Mackill D J. Marker-assisted selection: an approach for precision plant breeding in the twenty-first century., 2008, 363: 557–572.

        [10] Coleman R K, Gill G S, Rebetzke G J. Identification of quantitative trait loci for traits conferring weed competitiveness in wheat (L.)., 2001, 52: 1235–1246.

        [11] Nelson J C, Andreescu C, Breseghello F, Finney P L, Gualberto D G, Bergman C, Pe?a R J, Perretant M R, Leroy P, Qualset C O, Sorrells M E. Quantitative trait locus analysis of wheat quality traits., 2006, 149: 145–159.

        [12] Ramya P, Chaubal A, Kulkarni K, Gupta L, Kadoo N, Dhaliwal H S, Chhuneja P, Lagu M, Gupt V. QTL mapping of 1000-kernel weight, kernel length, and kernel width in bread wheat (L.)., 2010, 51: 421–429.

        [13] Ren Y, Huang L, Shao M, Sun L, Zhao K, Wang J, Xu X, Feng W, Wang J, Yan L, Wang S, Wang L. QTL mapping for plant height of wheat under different irrigation modes., 2017: 15–19.

        [14] 張坤普, 徐憲斌, 田紀春. 小麥籽粒產(chǎn)量及穗部相關性狀的QTL定位. 作物學報, 2009, 35: 270–278. Zhang K P, Xu X B, Tian J C. QTL mapping for grain yield and spike related traits in common wheat., 2009, 35: 270–278 (in Chinese with English abstract).

        [15] Zegeye H, Rasheed A, Makdis F, Badebo A, Ogbonnaya F C. Genome-wide association mapping for seedling and adult plant resistance to stripe rust in synthetic hexaploid wheat., 2014, 9: e105593.

        [16] Ding J, Ali F, Chen G, Chen G, Li H, Mahuku G, Yang N, Narro L, Magorokosho C, Makumbi D, Yan J. Genome-wide association mapping reveals novel sources of resistance to northern corn leaf blight in maize., 2015, 15: 206.

        [17] Hwang E Y, Song Q, Jia G, Specht J E, Hyten D L, Costa J, Cregan P B. A genome-wide association study of seed protein and oil content in soybean., 2014, 15: 1.

        [18] Zhang X, Zhang J, He X, Wang Y, Ma X, Yin D. Genome-wide association study of major agronomic traits related to domestication in peanut., 2017, 8: 1611.

        [19] Huang X, Wei X, Sang T, Zhao Q, Feng Q, Zhao Y, Li C, Zhu C, Lu T, Zhang Z, Li M, Fan D, Guo Y, Wang A, Wang L, Deng L, Li W, Lu Y, Weng Q, Liu K, Huang T, Zhou T, Jing Y, Li W, Lin Z, Buckler E S, Qian Q, Zhang Q, Li J, Han B. Genome-wide association studies of 14 agronomic traits in rice landraces., 2010, 42: 961.

        [20] Sukumaran S, Dreisigacker S, Lopes M, Chavez P, Reynolds M P. Genome-wide association study for grain yield and related traits in an elite spring wheat population grown in temperate irrigated environments., 2015, 128: 353–363.

        [21] Turuspekov Y, Baibulatova A, Yermekbayev K, Tokhetova L, Chudinov V, Sereda G, Ganal M, Griffiths S, Abugalieva S. GWAS for plant growth stages and yield components in spring wheat (L.) harvested in three regions of Kazakhstan., 2017, 17: 51–61.

        [22] 朱玉磊, 王升星, 趙良俠, 張德新, 胡建幫, 曹雪連, 楊亞杰, 常成, 馬傳喜, 張海萍. 以關聯(lián)分析發(fā)掘小麥整穗發(fā)芽抗性基因分子標記. 作物學報, 2014, 40: 1725–1732. Zhu Y L, Wang S X, Zhao L X, Zhang D X, Hu J B, Cao X L, Yang Y J, Chang C, Ma C X, Zhang H P. Exploring molecular markers of preharvest sprouting resistance gene using wheat intact spikes by association analysis., 2014, 40: 1725–1732 (in Chinese with English abstract).

        [23] Wang S X, Zhu Y L, Zhang D X, Shao H, Liu P, Hu J B, Zhang H, Zhang H P, Chang C, Lu J, Xia X C, Sun G L, Ma C X. Genome-wide association study for grain yield and related traits in elite wheat varieties and advanced lines using SNP markers., 2017, 12: e0188662.

        [24] Mwadzingeni L, Shimelis H, Rees D J G, Tsilo T J. Genome-wide association analysis of agronomic traits in wheat under drought-stressed and non-stressed conditions., 2017, 12: e0171692.

        [25] Liu Y, Lin Y, Gao S, Li Z, Ma J, Deng M, Chen G, Wei Y, Zheng Y. A genome-wide association study of 23 agronomic traits in Chinese wheat landraces., 2017, 91: 861.

        [26] Würschum T, Langer S M, Longin C F, Korzun V, Akhunov E, Ebmeyer E, Schachschneider R, Schacht J, Kazman E, Reif J C. Population structure, genetic diversity and linkage disequilibrium in elite winter wheat assessed with SNP and SSR markers., 2013, 126: 1477–1486.

        [27] Meng L, Li H, Zhang L, Wang J. QTL IciMapping: integrated software for genetic linkage map construction and quantitative trait locus mapping in biparental populations., 2015, 3: 269–283.

        [28] 陳昆松, 李方, 徐昌杰, 張上隆, 傅承新. 改良CTAB法用于多年生植物組織基因組DNA的大量提取. 遺傳, 2004, 26: 529–531. Chen K S, Li F, Xu C J, Zhang S L, Fu C X. An efficient macro-method of genomic DNA isolation fromleaves., 2004, 26: 529–531 (in Chinese with English abstract).

        [29] Allen A M, Winfield M O, Burridge A J, Downie R C, Benbow H R, Barker G L, Wilkinson P A, Coghill J, Waterfall C, Davassi A, Scopes G, Pirani A, Webster T, Brew F, Bloor C, Griffiths S, Bentley A R, Alda M, Jack P, Phillips A L, Edwards K J. Characterization of a wheat breeders’ array suitable for high-throughput SNP genotyping of global accessions of hexaploid bread wheat ()., 2017, 15: 390–401.

        [30] Botstein D, White R, Skolnick M, Davis R. Construction of a genetic linkage map in man using restriction fragment length polymorphisms., 1980, 32: 314–331.

        [31] Evanno G, Regnaut S, Goudet J. Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study., 2005, 14: 2611–2620.

        [32] Neumann K, Kobiljski B, Den?i? S, Varshney R K, B?rner A. Genome-wide association mapping: a case study in bread wheat (L.)., 2011, 27: 37–58.

        [33] Kuznetsova A, Brockhoff P B, Christensen R H B. lmerTest: Tests for random and fixed effects for linear mixed effect models (lmer objects of lme4 package). 2013. https://rdrr.io/cran/lmerTest/man/lmerTest-package.html.

        [34] Sukumaran S, Xiang W, Bean S R, Pedersen J F. Association mapping for grain quality in a diverse sorghum collection., 2012, 5: 126–135.

        [35] Breseghello F, Sorrells M E. Association mapping of kernel size and milling quality in wheat (L.) cultivars., 2006, 172: 1165–1177.

        [36] Yu J, 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: 203–208.

        [37] 要燕杰. 我國冬小麥品種(系)遺傳多樣性及部分產(chǎn)量性狀和品質(zhì)性狀位點的關聯(lián)分析. 西北農(nóng)林科技大學碩士學位論文, 陜西楊凌, 2014. Yao Y J. Genetic Diversity and Association Analysis of Loci for Yield and Quality Traits in Wheat. MS Thesis of Northwest A&F University, Yangling, Shaansi, China, 2014 (in Chinese with English abstract).

        [38] Jaiswal V, Gahlaut V, Meher P K, Mir R R, Jaiswal J P, Rao A R, Balyan H S, Gupta P K. Genome wide single locus single trait, multi-locus and multi-trait association mapping for some important agronomic traits in common wheat (.L.)., 2016, 11: e0159343.

        [39] Shi W, Hao C, Zhang Y, Cheng J, Zhang Z, Liu J, Yi X, Cheng X, Sun D, Xu Y, Zhang X, Cheng S, Guo P, Guo J. A combined association mapping and linkage analysis of kernel number per spike in common wheat (L.)., 2017, 8: 1412.

        [40] Zanke C, Ling J, Plieske J, Kollers S, Ebmeyer E, Korzun V, Argillier O, Stiewe G, Hinze M, Neumann F, Eichhorn A, Polley A, Jaenecke C, Ganal M W, R?der M S. Analysis of main effect QTL for thousand grain weight in European winter wheat (L.) by genome-wide association mapping., 2015, 6: 644.

        [41] Sukumaran S, Reynolds M P, Sansaloni C. Genome-wide association analyses identify QTL hotspots for yield and component traits in durum wheat grown under yield potential, drought, and heat stress environments., 2018, 9: 81.

        [42] Ain Q, Rasheed A, Anwar A, Mahmood T, Imtiaz M, Mahmood T, Xia X, He Z, Quraishi U M. Genome-wide association for grain yield under rainfed conditions in historical wheat cultivars from Pakistan., 2015, 6: 743.

        [43] Zanke C, Ling J, Plieske J, Kollers S, Ebmeyer E, Korzun V, Argillier O, Stiewe G, Hinze M, Beier S, Ganal M W, R?der M S. Genetic architecture of main effect QTL for heading date in European winter wheat., 2014, 5: 217.

        [44] 靖華, 亢秀麗, 馬愛平, 崔歡虎, 王娟玲, 劉建華. 晉南旱垣春季低溫對不同播種期小麥凍害的影響. 中國農(nóng)學通報, 2011, 27(9): 76–80 Jing H, Kang X L, Ma A P, Cui H H, Wang J L, Liu J H. Effect of spring low temperature on different sowing date winter wheat frozen injury on the arid area of southern Shanxi province., 2011, 27(9): 76–80 (in Chinese with English abstract).

        [45] Pei Z M, Murata Y, Benning G, Thomine S, Klüsener B, Allen G J, Grill E, Schroeder J I. Calcium channels activated by hydrogen peroxide mediate abscisic acid signalling in guard cells., 2000, 406: 731–734.

        [46] Van Ooijen G, Mayr G, Kasiem M M, Albrecht M, Cornelissen B J, Takken F. Structure-function analysis of the NB-ARC domain of plant disease resistance proteins., 2008, 59: 1383–1397.

        [47] Hong M J, Lee Y M, Son Y S, Im C H, Yi Y B, Rim Y G, Bahk J D, Heo J B. Rice Rab11 is required for JA-mediated defense signaling., 2013, 434: 797–802.

        [48] Graaf B H J D, Cheung A Y, Andreyeva T, Levasseur K, Kieliszewski M, Wu H M. Rab11 GTPase-regulated membrane trafficking is crucial for tip-focused pollen tube growth in tobacco., 2005, 17: 2564–2579.

        [49] Asaoka R, Uemura T, Ito J, Fujimoto M, Ito E, Ueda T, Nakano A. Arabidopsis RABA1 GTPases are involved in transport between the trans-Golgi network and the plasma membrane, and are required for salinity stress tolerance., 2013, 73: 240–249.

        [50] Keates S E, Kostman T A, Anderson J D, Bailey B A. Altered gene expression in three plant species in response to treatment with Nep1, a fungal protein that causes necrosis., 2003, 132: 1610–1622.

        [51] 蘇延萍, 路小鐸, 沈頌東, 張春義. 擬南芥基因影響葉綠體的形成. 植物學報, 2011, 46: 379–385. Su Y P, Lu X D, Shen S D, Zhang C Y.is essential for chloroplast biogenesis in Arabidopsis., 2011, 46: 379–385 (in Chinese with English abstract).

        Genome-wide association study for main agronomic traits in common wheat

        ZHAI Jun-Peng, LI Hai-Xia, BI Hui-Hui, ZHOU Si-Yuan, LUO Xiao-Yan, CHEN Shu-Lin, CHENG Xi-Yong*, and XU Hai-Xia*

        Henan Agricultural University / National Key Laboratory of Wheat and Maize Crop Science / Collaborative Innovation Center of Henan Grain Crops, Zhengzhou 450046, Henan, China

        To illustrate the genetic mechanism of complex agronomic traits in wheat, we investigated nine agronomic traits using 150 wheat cultivars (lines) from China across four environments. Genome-wide association analysis was performed using wheat 35K genotyping assay with five association models (Q, PCA, K, PCA+K, Q+K). The results revealed that the polymorphic information content (PIC) of values was between 0.0950 and 0.5000, and the minimum allele frequency (MAF) was between 0.0500 and 0.5000. Both the population structure analysis and the PCA analysis showed that the tested materials could be divided into two sub-populations. Linkage disequilibrium analysis found that the LD decay distances of the A, B, D, and the whole genome were approximately 4.7, 8, 11, and 6 Mb, respectively. A total of 652 significant (≤ 0.001) marker-trait associations (MTAs) were detected, thereinto, 21 SNPs could be detected on chromosomes 1A(1), 1B(4), 2A(3), 2D(2), 3A(1), 5A(1), 5B(5), 6A(1), 6B(2), and 7D(3) in two or more environments. Three SNPs were significantly associated with two traits and the physical position of one SNP was unknown. Single SNPs could explain 7.67 % to 18.79 % phenotypic variation. It was found that eight favorable allelic variations accounted for a low proportion in the tested population. Fourteen candidate genes that may be related to agronomic traits of wheat were identified. Among them,,, andmay play important roles in plants resistance to biotic and abiotic stress.andmay be related to the hormones synthesis and response in plants.may enhance the cell wall formation of plants.may be involved in chloroplast development. The function of the other seven candidate genes is unknown.

        wheat; agronomic trait; linkage disequilibrium; genome-wide association study; candidate gene

        本研究由國家重點研發(fā)計劃項目(2017YFD0100706), 國家轉基因生物新品種培育重大專項(2016ZX08002003- 004)和國家重點基礎研究發(fā)展計劃(973計劃)項目(2014CB138105)資助。

        This study was supported by the National Key Research and Development Program of China (2017YFD0100706), the National Major Project for Developing New GM Crops (2016ZX08002003-004), and the National Key Basic Research Program (973 Program) (2014CB138105).

        許海霞, E-mail: hauxhx@henau.edu.cn; 程西永, E-mail: xyc634@163.com

        E-mail: 597452880@qq.com

        2019-01-06;

        2019-04-15;

        2019-05-14.

        10.3724/SP.J.1006.2019.91002

        URL: http://kns.cnki.net/kcms/detail/11.1809.S.20190511.1229.002.html

        猜你喜歡
        農(nóng)藝染色體基因組
        牛參考基因組中發(fā)現(xiàn)被忽視基因
        農(nóng)機需要農(nóng)藝“標準”,農(nóng)藝需要農(nóng)機“靈活”——2021國際農(nóng)機展不容錯過的果蔬茶機械化論壇
        落葉果樹(2021年6期)2021-02-12 01:29:20
        水稻主要農(nóng)藝性狀與產(chǎn)量的相關性分析
        多一條X染色體,壽命會更長
        科學之謎(2019年3期)2019-03-28 10:29:44
        為什么男性要有一條X染色體?
        科學之謎(2018年8期)2018-09-29 11:06:46
        14份蔞蒿種質(zhì)資源主要農(nóng)藝性狀及營養(yǎng)成分評價
        中國蔬菜(2016年8期)2017-01-15 14:23:37
        能忍的人壽命長
        90團舉辦初級農(nóng)藝工培訓班
        再論高等植物染色體雜交
        基因組DNA甲基化及組蛋白甲基化
        遺傳(2014年3期)2014-02-28 20:58:49
        又爽又黄无遮挡高潮视频网站| 变态调教一区二区三区女同| 亚洲一区二区三区影院| 亚洲av无码不卡久久| 久久国产亚洲高清观看5388| 看全色黄大黄大色免费久久| 蜜桃免费一区二区三区| 四川丰满妇女毛片四川话| 中文字幕亚洲乱码熟女一区二区 | 91超碰在线观看免费| 日本一二三区在线视频观看| 色综合久久中文字幕综合网| www插插插无码视频网站| 人妻熟妇乱系列| 日本熟妇免费一区二区三区| 77777亚洲午夜久久多喷| 人妻av无码系列一区二区三区| 亚洲AV无码成人精品区网页| 国产av三级精品车模| 激情综合色综合啪啪开心| 亚洲中文字幕无码久久| 精品人妻少妇一区二区中文字幕| 国产一区在线视频不卡| 92午夜少妇极品福利无码电影| 国产精品va在线观看无码| 亚洲欧美日本人成在线观看| 亚洲精品国产av日韩专区| 97在线观看播放| 另类亚洲欧美精品久久不卡 | 色狠狠色噜噜av天堂一区| 国产精品亚洲欧美云霸高清| 偷柏自拍亚洲综合在线| 日日碰日日摸日日澡视频播放 | 噜噜综合亚洲av中文无码| 日韩精品人妻系列无码专区免费| 久久精品国产精品亚洲艾| 99久久国内精品成人免费| 色橹橹欧美在线观看视频高清 | 国产精品成人久久一区二区| 午夜视频在线观看一区二区小| 亚洲欧洲∨国产一区二区三区|