999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

甜玉米籽粒體積和粒重的全基因組關聯分析

2022-07-25 06:38:24肖穎妮李高科于永濤李光玉高穎珊胡建廣
中國農業大學學報 2022年7期
關鍵詞:關聯

肖穎妮 李高科 李 坤 于永濤 李光玉 李 武 高穎珊 胡建廣

(廣東省農業科學院 作物研究所/廣東省農作物遺傳改良重點實驗室,廣州 510640)

甜玉米是普通玉米中1個或幾個基因發生突變而形成的一種特用玉米,其營養豐富,適口性好,風味獨特,深受廣大消費者青睞。隨著人們生活水平的提高,甜玉米由于容易種植、經濟效益好、生長周期短而深受種植農戶的青睞,種植面積逐年增加,2019年,我國種植面積已達53萬hm。高產穩產一直都是作物育種重要的目標。粒重和籽粒大小作為產量的主要構成因子,與產量直接相關。甜玉米籽粒鮮重和干重顯著正相關(

r

=0.89,

P

=2.2×10),提高甜玉米干粒重不僅可以提高制種產量,也意味著鮮穗產量的提高,增加籽粒體積一定程度上也可以提高產量。因此,深入剖析甜玉米籽粒體積和粒重的遺傳基礎對提高甜玉米產量具有極其重要的意義。大量的籽粒突變體為克隆調控玉米籽粒粒重和籽粒大小的基因提供了寶貴的遺傳資源。通過突變體克隆影響籽粒粒重和籽粒大小的基因主要包含3類:第一類是胚乳中淀粉合成代謝相關基因(

sh2

sh1

bt1

bt2

su1

等),這些基因突變主要導致胚乳淀粉含量下降,從而使得粒重降低;第二類是影響胚發育的基因(

emp2

emp12

dek14

dek15

等),這些基因突變導致胚敗育,最終籽粒發育不良;第三類是基部胚乳轉移層區相關基因(

mn1

),這類基因突變主要導致籽粒變小,粒重也隨之降低。但是,突變體導致籽粒發生劇烈變化,極端的表型很難直接用于籽粒的遺傳改良。利用QTL定位以及全基因組關聯分析方法是解析數量性狀遺傳基礎的重要手段,截止目前,已經定位了上千個與玉米粒重和籽粒大小相關的QTL (www.maizegdb.org/qtl)或位點,并成功克隆了4個與粒重相關基因,分別為編碼己糖轉運蛋白的

ZmSWEET4c

、編碼細胞壁轉移酶的

ZmINCW1

、編碼BAM相關蛋白激酶的

ZmBAM1d

和編碼三角狀五肽重復(PPR)蛋白的

qKW9

隨著玉米基因組參考序列的測序完成以及測序技術的發展,全基因組關聯分析(genome-wide association study, GWAS)已成功解析玉米多個復雜數量性狀的遺傳基礎,包含籽粒品質性狀、植株抗逆性狀和株型性狀等。然而,對甜玉米籽粒體積和粒重遺傳基礎解析的研究卻鮮少報道。本研究以209份具有廣泛代表性的甜玉米自交系為關聯群體,結合全基因組重測序獲得980萬個SNP基因型,對甜玉米籽粒體積和粒重進行全基因組關聯分析,挖掘控制甜玉米籽粒體積和粒重的相關基因,以期為甜玉米高產育種提供參考。

1 材料與方法

1.1 試驗材料及田間試驗

供試甜玉米材料共包含209份甜玉米自交系,其中,國內來源158份,國外來源51份(表1),由廣東省農業科學院作物研究所提供。所有材料在廣州市天河區連續種植3年(2017、2018和2019年),采用完全隨機區組設計,每個材料種1 行,每行11株,行距0.5 m,株距0.25 m,按照當地栽培方法進行田間管理。

表1 本研究所使用的甜玉米自交系材料
Table 1 Sweet corn inbreds used in this study

編號ID名稱Name地理來源Original region編號ID名稱Name地理來源Original regionC11022中國C41新美中國C21132中國C42國區中國C3C5-1中國C44超甜201中國C4日超1日本C48航1132中國C5C4中國C49航日超日本C6794中國C50改良794中國C7日超2日本C51群2-7中國C9夏威夷美國C52群2白-3中國C10鮮美選系中國C53群2白-4中國C11華珍1中國C5408秋群1耐寒后代-2中國C12華珍3中國C5508秋群1耐寒后代-3中國C13265OR1美國C5608秋群1耐寒后代-4中國C15世珍中國C5708秋群1耐寒后代-5中國C16中農419中國C58王朝開放美國C17Z5阿根廷C59超引黃中國C18N16中國C60奧甜8706美國C20MH70美國C61金煌1號泰國C21庫普拉選系美國C62國區T12-1中國C22金珠蜜美國C63國區T12-2中國C23金珠蜜白美國C64GLAY10黃-2美國C24曾061487泰國C65曾引me153-1泰國C25南引3中國C66曾引me153-2泰國C26群1穗選1中國C67QX331開放選黃中國C27群1穗選2中國C68QX331開放選白中國C28群1穗選3中國C69熱群-3中國C30農甜88白中國C70群2-8中國C31農甜88黃中國C71群2-9中國C32225美國C72群2-10中國C33群1穗選4中國C73群2-11中國C35熱群穗選中國C7408秋群1耐寒后代-6中國C36群1穗選6中國C75MH70×誘選甜中國C37SSC016美國C76群2甜誘選甜白中國C38奧甜美國C77李美引M6美國C39田雜中國C78開放華珍-1中國C40穗甜202中國C79開放華珍-2中國

表1(續)

編號ID名稱Name地理來源Original region編號ID名稱Name地理來源Original regionC8108秋泰引親本泰國C118群2-1中國C8208秋泰引開放泰國C119群2-2中國C8308秋泰引組合-2泰國C120群2-3中國C8408秋泰引組合-3泰國C12108秋耐寒后代中國C8508秋泰引組合-4泰國C12208秋耐寒后代白-1中國C86HSC0803美國C12312春高引白-1中國C87HSC0803白美國C12412春高引-2白中國C88金銀08白美國C12512春高引-3白中國C89田蜜6號中國C12612春高引-4白中國C90黃金15-2美國C12712春高引-5白中國C91泰引超甜-2泰國C12812春高引-6-1中國C92泰引白泰國C12912春高引-6-2中國C93泰引6高泰國C131GLAY10黃-1美國C94HIBRIX51大粒泰國C132群1-3中國C96湖北引-5白-1中國C133群2-4中國C97湖北引-5白-2中國C134群2-5中國C98湖北引-10中國C135群2-6中國C9909引-21中國C13709秋國區T9中國C10010春國區T9中國C138群1甜誘白中國C10110春國區T9白中國C139MH70×誘中國C102華珍母本中國C140群2甜誘-1中國C10309國區T4中國C141群2甜誘白中國C104群1-1中國C142群1甜誘 中國C105群2白-1中國C143群2甜誘 混中國C106群2白-2中國C144群2甜誘 單穗中國C107德甜36白中國C145群2甜誘-2中國C108S103中國C146MH80開放美國C109大28-2中國C14708秋泰引組合-1泰國C11012春高引9-1中國C148黃金15-1美國C11112春高引9-2中國C149泰引超甜-4泰國C112豐蜜中國C15009引-2-1中國C113群1-2中國C15109引-2-2中國C114熱群-1中國C15209引-11中國C115熱群-2中國C15510秋柯引彩-4中國C116群1-1-1中國C15710秋柯引彩-7白中國C117群1-1-2中國C158庫普拉×誘美國

表1(續)

編號ID名稱Name地理來源Original region編號ID名稱Name地理來源Original regionC15911春GT3-1中國C202華寶8號中國C16011春GT3-2中國C203泰豐母泰國C16111春GT8-1中國C204泰豐父泰國C16211春GT8-2中國C2525#-7中國C16311春GT11中國C253A453海-4-2中國C164群1-1-3中國C256A003海中國C16508秋耐寒后代白-2中國C25704冬Z15-1中國C166SB903美國C25804冬Z19中國C167921似親本泰國C25904冬Z52-1中國C169群1-4中國C260HN-6中國C171華美甜似親本白-1中國C261xsj2號中國C172華美甜似親本白-2中國C264A126-1中國C173美國225美國C272ZHU2中國C175山西引超甜中國C274JFM中國C17612春高引-7中國C27507SHXJF中國C17712春高引-10中國C27610QXJFM中國C17812春高引-10白-1中國C277XJLN中國C17912春高引-10白-2中國C281737中國C18012春高引-13中國C2825#中國C18112春高引-17中國C2841167-廣西TF中國C182奧費蘭選系美國C285AB1000 白粒超甜中國C183SC1388泰國C28936-1 白粒超甜中國C18608秋群1耐寒后代-1中國C292A02中國C18708秋群1耐寒后代×誘選甜中國C294A04中國C189群2 葉披,多中國C307A17中國C193農寶白中國C308A18中國C194麥哥娜姆選系美國C311A21 白粒超甜中國C195華威2045泰引開放系中國C315A25中國C196新美208中國C317A27中國C197SC858泰國C318A28 白粒超甜中國C198Yes37泰國C320A30中國C199改良泰引組合6白泰國C323A33中國C200仲鮮甜3號中國C329A39 白粒超甜中國C201群1白中國

1.2 表型調查及統計分析

所有材料自交授粉,收獲后自然晾干至含水量13%以下,每個材料選取5個均勻一致的果穗混合脫粒,取其中300粒裝袋,再從中隨機選取30粒干籽粒進行后續表型鑒定。將籽粒去麩后,使用分析天平稱量30粒籽粒重量(

w

),計算得單粒重(

w

/30),g;采用排酒精法對籽粒體積進行測量,取30粒干籽粒倒入裝有酒精的滴定管中,讀取加入籽粒前后液面刻度變化的差值,即為30粒籽粒體積(

V

),計算得單粒體積(

V

/30), mL。利用R軟件對3年表型進行統計分析和相關性分析,并計算廣義遺傳力和最佳線性無偏預測值(best linear unbiased prediction, BLUP)。

廣義遺傳力計算公式為:

(1)

式(1)中:代表遺傳變異,代表誤差,

e

代表環境個數。

BLUP計算公式為:

y

=

μ

+

g

+

e

+

ε

(2)

式(2)中:

y

代表第

i

個家系的表型值,

μ

代表所有環境下的均值,

g

代表遺傳效應,

e

代表環境效應,

ε

代表隨機誤差。

1.3 基因型分析

取每個自交系3~5棵植株的幼嫩葉片,采用改良的CTAB方法提取DNA,利用illumina HiSeq 2500平臺對209份材料進行深度約11.5×的全基因組重測序,用BWA軟件將測序序列比對到玉米‘B73’參考基因組(AGPv4),而后用SAMTOOLS軟件對SNP進行分析,最后用GATK軟件對SNP進行提取,保留缺失率≤20%以及最小等位基因頻率≥5%的SNP,最終獲得980萬個覆蓋全基因組的高質量的SNP(數據未發表)用于后續分析。

1.4 全基因組關聯分析及候選基因提名

將基因型和表型導入TASSEL5.0中,采用一般線性模型(General Linear Model, GLM)進行全基因組關聯分析。由于非獨立的SNP之間有較強的LD距離,為了進一步更準確的確定閾值,利用PLINK軟件(window size 50, step size 50,

R

≥0.2)計算得到767 964個獨立標記數,因此,將關聯分析的顯著性閾值設為

P

<1×10(1/767 964)。根據之前對甜玉米自交系LD的計算,甜玉米的衰減距離為500 kb,因此,當500 kb區間內存在多個顯著SNP時,合并到1個最顯著的SNP,查找其上下游500 kb內的基因,根據基因注釋并結合基因表達譜數據,找尋與其他植物籽粒體積或粒重有關基因為候選基因,如無相關基因,提名離顯著位點最近且檢測到表達量有差異的基因為候選基因。針對每個表型,所有顯著位點解釋的表型變異用R軟件中的

lm

進行線性回歸分析計算而得。

1.5 候選基因表達量與表型的相關分析

所有材料自交授粉,采摘授粉后15 d(15 DAP)的果穗,分別選取每個自交系的3~4穗,每穗取1~2粒籽粒,混合提取籽粒RNA構建總RNA文庫,利用Illumina Hiseq 2500平臺對文庫進行150 bp的雙端轉錄組測序,每個樣品獲得了大約24.76 M高質量的Reads,利用R軟件對每個樣品中所有基因的FPKM進行計算,最終獲得了27 133個基因在不同自交系中的表達譜,利用R軟件對候選基因在不同自交系中的表達量和其對應的表型進行相關分析。

2 結果與分析

2.1 甜玉米籽粒體積和粒重表型分析

由圖1可知,籽粒體積變異范圍為0.06~0.15 mL,變異倍數為2.5倍,單粒重變異范圍為0.08~0.17 g,變異倍數為2.1倍(表2),表明這2個性狀在群體中都具有廣泛的表型變異,并且這2個性狀在群體中都呈正態分布(圖1)。此外,籽粒體積和粒重呈顯著正相關(

r

=0.92)(圖1)。利用3年的表型進行遺傳力分析,籽粒體積和粒重的遺傳力分別為0.70和0.79(表2),表明2個性狀的表型變異主要受遺傳因素控制,受環境因素影響較小。

(a)籽粒粒重表型分布;(b)籽粒體積表型分布;(c)粒重和體積的相關性分析(a) Distribution of kernel weight; (b) Distribution of kernel volume; (c) The correlation coefficients between kernel weight and kernel volume圖1 甜玉米群體籽粒體積和粒重的表型分析Fig.1 Phenotypic analysis of kernel volume and kernel weight in the sweet corn panel

表2 甜玉米籽粒體積和粒重的表型分析和遺傳力分析
Table 2 Phenotypic variation and heritability of two traits in sweet corn kernels

性狀Trait自交系個數No. of lines最小值Min最大值Max平均值±標準差Mean±SD廣義遺傳力H2籽粒體積/mLKV2090.060.150.10±0.010.70單粒重/gKW2090.080.170.12±0.020.79

2.2 籽粒體積和粒重全基因組關聯分析

由表3可知,利用GLM模型,當閾值為

P

<1×10時,共檢測到15個顯著SNP位點,其中,籽粒體積檢測到9個位點,單個位點解釋的表型變異范圍為10.6%~13.6%,總共解釋55.7%的籽粒體積表型變異;籽粒粒重檢測到10個位點,單個位點解釋的表型變異范圍為9.3%~12.6%,總共解釋52.1%的粒重表型變異。由表4可知,在這15個位點中,有4個位點和2個性狀均顯著關聯,表現“一因多效” 性(圖2),這與籽粒體積和粒重呈顯著正相關的結果是一致的(圖1)。

(a)籽粒體積的曼哈頓圖;(b)籽粒體積的Q-Q圖;(c)單粒重的曼哈頓圖;(d)單粒重的Q-Q圖(a) Manhattan plot of kernel volume; (b) Quantile-quantile plot of kernel volume; (c) Manhattan plot of single kernel weight; (d) Quantile-quantile plot of single kernel weight水平虛線代表閾值(1×10-6)。The dashed horizontal line depicts the Bonferroni significance threshold (1×10-6).圖2 甜玉米籽粒體積和粒重全基因組關聯分析Fig.2 Genome-wide association study of kernel volume and kernel weight in sweet corn

表3 甜玉米籽粒性狀顯著位點
Table 3 Summary of significant SNPs identified for kernel traits in sweet corn

性狀Trait顯著位點數目Significant SNPs單個位點解釋表型變異/%Variation explained by each SNP所有位點解釋表型變異/%Variation explained by all SNPs籽粒體積/mLKV910.6~13.655.7單粒重/gKW109.3~12.652.1

2.3 與普通玉米研究位點比較

由表4可知,3個與籽粒體積關聯的位點均位于普通玉米籽粒大小相關性狀的QTL區間內。其中,S1_204753078位點與籽粒體積和粒重均顯著關聯, S1_288353444與籽粒體積顯著關聯,這2個位點均落在控制粒寬的QTL區間內,S9_24601539與籽粒體積顯著關聯,該位點落在粒長的QTL區間內;此外,S1_291060760與粒重顯著關聯,該位點落在普通玉米千粒重的QTL區間內。綜上,通過檢測的位點和普通玉米籽粒有關性狀位點進行比較分析,說明利用全基因組關聯分析方法挖掘與甜玉米籽粒性狀顯著關聯位點方法可靠。

2.4 候選基因預測和功能分析

基于玉米參考基因組數據,對顯著關聯位點上下游500 kb范圍內的基因進行預測分析,結合甜玉米籽粒15 DAP的表達譜數據,最后挖掘到15個候選基因。大部分候選基因都是轉錄因子等調控基因,見表4,比如控制籽粒體積的Zm00001 d034277和Zm00001 d011730分別為NAC轉錄因子和乙醇脫氫酶轉錄因子,控制粒重的Zm00001 d030577是bZIP轉錄因子,此外,控制籽粒體積和粒重的Zm00001 d034369和Zm00001 d011639分別是環指蛋白和脫落酸不敏感轉錄因子。

1號染色體上的位點S1_204753078與籽粒體積、粒重均顯著關聯,該位點位于基因

Br2

(Zm00001 d031871)第3個內含子上,

Br2

基因編碼生長素轉運蛋白。通過對該顯著位點進行單倍型分析發現,基因型CC與基因型TT的植株相比,株高有一定程度的增加,穗位高、籽粒體積和粒重均有顯著增加,因此,

Br2

和甜玉米籽粒體積、粒重均顯著關聯,并且,大部分自交系都包含有利等位基因CC(n=129),見圖3。由圖4可知,不同玉米材料授粉15 d 后

Br2

表達量與籽粒體積和粒重均呈顯著正相關(

r

=0.37,

P

=1.18×10和

r

=0.36,

P

=1.28×10),說明該基因是通過調控基因的表達量從而引起表型發生變化的,亦即隨著該基因表達量的增加,籽粒體積和粒重均隨之增加。

圖3 Br2上顯著位點與不同農藝性狀的方差分析Fig.3 ANOVA based on the most significant associated SNP of Br2 with agronomic traits

圖4 Br2基因表達量與籽粒體積(a)和單粒重(b)的相關性分析Fig.4 Correlation of expression level of Br2 with kernel volume (a) and kernel weight(b)

3 討 論

3.1 甜玉米籽粒體積和粒重的全基因組關聯分析

盡管甜玉米育種更多關注的是其風味和營養品質,然而提高產量也是甜玉米育種研究的重要目標之一,籽粒體積和粒重是與甜玉米產量相關的最重要的農藝性狀。本研究對該關聯群體進行采收期甜玉米籽粒含水量的測定,發現籽粒的含水量平均為75%,變異范圍為61%~85%,籽粒干粒重和鮮重相關性極顯著(

r

=0.89,

P

=2.2×10),因此,甜玉米籽粒干重的增加,意味著籽粒鮮重也增加,有助于提高甜玉米的產量。

GWAS是解析復雜性狀的遺傳結構以及挖掘關鍵基因的重要手段,本研究連續3年對209份甜玉米自交系的干籽粒體積和粒重進行了測定,共檢測到了15個顯著位點,但單個位點的遺傳效應均不高(9.3%~13.6%),表明這2個性狀主要受微效多基因控制,符合產量是復雜數量性狀的特點。通過與徐運林等報道的甜玉米百粒重的關聯SNP進行比較,并未發現有共同檢測到的位點,這可能因為所使用的材料不同,也說明了粒重性狀的復雜性。

3.2 甜玉米籽粒體積和粒重的候選基因挖掘

本研究中檢測到

Br2

基因與籽粒體積、粒重均顯著關聯,該基因編碼的是一個生長素極性運輸蛋白PGP;生長素轉運蛋白和信號因子通過調節植物生長素的極性運輸對植物的生長發育起著重要作用,玉米

de

18

突變體中生長素合成受到影響,從而抑制了細胞的分裂,最終使種子粒重降低。在普通玉米中,

Br2

基因突變體主要通過減少植株縱向細胞分裂,導致植株變矮,穗位變低,并且不影響產量,此外,Li等發現在普通玉米自然變異群體中,

Br2

與籽粒長度顯著關聯。在甜玉米中,檢測到與籽粒體積、粒重最顯著關聯的位點位于基因

Br2

內含子上,這與Li等的研究結果類似,該位點可能并不是影響籽粒體積和粒重的功能位點,而是與功能位點存在連鎖不平衡(LD)而被檢測到的。

Br2

的表達量和籽粒體積以及粒重均顯著相關(

r

=0.37和

r

=0.36),這與該位點解釋的表型變異(13.6%和10.6%)的結果相符,說明該位點并不是主效位點,對籽粒體積和粒重的變異只起到微效作用,推測功能位點可能存在于影響基因表達水平的調控區(5’UTR或3’UTR),通過調控基因的表達水平,引起粒長的變化,最終導致籽粒體積和粒重都增加,因此,還需要通過候選基因重測序以及基因驗證等手段進一步來驗證。此外,在甜玉米中,該基因在調控籽粒產量(單粒重和籽粒體積)的同時,也調控穗位高,而對株高的影響相對較小(圖3),因此,通過合理選育,有望選育出籽粒體積和粒重均增加,而株高和穗位高相對合理的優良品種。本研究在9號染色體上發現1個與籽粒體積關聯的候選基因Zm00001 d045499,該基因位于Liu等報道的控制籽粒粒長的QTL區間內,編碼1個E3泛素蛋白連接酶,水稻中克隆的控制粒寬和粒重的

GW2

基因也是1個E3泛素蛋白連接酶,GW2蛋白缺失后,穎殼內細胞數量增加,從而導致籽粒變寬、變重。此外,在擬南芥中,

DA2

編碼1個E3泛素蛋白連接酶,功能缺失突變體的種子體積由于細胞分化未被抑制而顯著增加。Zm00001 d045499在甜玉米中調控籽粒體積變化的作用機制還有待進一步驗證。

4 結 論

本研究結果表明,本群體(209份甜玉米自交系)包含廣泛的遺傳變異,結合覆蓋全基因組高密度的980萬個SNP標記,利用GLM模型通過全基因組關聯分析方法來挖掘控制甜玉米產量性狀的關聯位點,共檢測到15個顯著關聯位點,大部分候選基因是轉錄因子;其中,

Br2

與籽粒體積、粒重均顯著關聯,該基因主要通過提高表達量來增加甜玉米籽粒體積和粒重。

猜你喜歡
關聯
不懼于新,不困于形——一道函數“關聯”題的剖析與拓展
“苦”的關聯
當代陜西(2021年17期)2021-11-06 03:21:36
船山與宋學關聯的再探討
原道(2020年2期)2020-12-21 05:47:06
“一帶一路”遞進,關聯民生更緊
當代陜西(2019年15期)2019-09-02 01:52:00
新制度關聯、組織控制與社會組織的倡導行為
奇趣搭配
基于廣義關聯聚類圖的分層關聯多目標跟蹤
自動化學報(2017年1期)2017-03-11 17:31:17
智趣
讀者(2017年5期)2017-02-15 18:04:18
探討藏醫學與因明學之間的關聯
西藏科技(2016年5期)2016-09-26 12:16:39
GPS異常監測數據的關聯負選擇分步識別算法
主站蜘蛛池模板: 国产在线精彩视频二区| 熟妇无码人妻| 永久免费av网站可以直接看的 | 无码在线激情片| 国产精品第5页| 午夜视频在线观看免费网站| 青青草原国产| 57pao国产成视频免费播放| AV无码一区二区三区四区| jijzzizz老师出水喷水喷出| 尤物视频一区| 国产网站免费| 无码福利日韩神码福利片| 九九热这里只有国产精品| 91丝袜乱伦| 国产亚洲精品97AA片在线播放| 女人毛片a级大学毛片免费| 在线日韩日本国产亚洲| 五月婷婷亚洲综合| 91色综合综合热五月激情| 又猛又黄又爽无遮挡的视频网站| 久久青青草原亚洲av无码| 国产精品99r8在线观看| 老司国产精品视频| 亚洲动漫h| 亚洲国产在一区二区三区| 亚洲国产成熟视频在线多多 | 欧美在线导航| 亚洲另类第一页| 亚洲人成影院午夜网站| 亚洲天堂日韩av电影| 午夜视频免费试看| 亚洲黄色片免费看| 久久久久青草大香线综合精品| 欧美午夜一区| a天堂视频| 亚洲伊人久久精品影院| 国产精选自拍| www.亚洲一区二区三区| 国产在线精品网址你懂的| 人妻一本久道久久综合久久鬼色| 毛片网站免费在线观看| 国产一级裸网站| 亚洲侵犯无码网址在线观看| 久久久受www免费人成| 国产99久久亚洲综合精品西瓜tv| 欧美成人精品在线| 亚洲精品高清视频| 亚瑟天堂久久一区二区影院| 欧美激情,国产精品| 四虎影视8848永久精品| 九九视频免费看| 无码精品福利一区二区三区| 婷婷色在线视频| 黄色网站不卡无码| 东京热av无码电影一区二区| 日本国产精品| 日本免费新一区视频| 国产精品亚洲一区二区在线观看| 亚洲熟女中文字幕男人总站| 五月丁香在线视频| 欧美一区二区福利视频| 91在线精品麻豆欧美在线| 色噜噜狠狠狠综合曰曰曰| 久久91精品牛牛| 国产成人精品高清不卡在线| 一级一级特黄女人精品毛片| 2021国产v亚洲v天堂无码| 精品三级在线| 国产成本人片免费a∨短片| 在线播放真实国产乱子伦| 一级毛片免费不卡在线| 9久久伊人精品综合| 男人天堂伊人网| 国产sm重味一区二区三区| 国产原创自拍不卡第一页| 在线亚洲精品福利网址导航| 午夜人性色福利无码视频在线观看| 一本二本三本不卡无码| 久久无码免费束人妻| 热久久综合这里只有精品电影| 国产日本视频91|