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

栗雜交F1代群體遺傳結構及其農藝性狀關聯分析

2022-08-03 10:57:40江錫兵章平生張東北吳仁超吳聰連賴俊聲龔榜初
林業科學研究 2022年4期
關鍵詞:關聯分析

江錫兵,章平生,張東北,吳仁超,吳 劍,吳聰連,賴俊聲,龔榜初*

(1.中國林業科學研究院亞熱帶林業研究所,浙江 杭州 311400;2.浙江省慶元縣自然資源和規劃局,浙江 慶元 323800)

殼斗科(Fagaceae)栗屬(Castanea Mill.)植物共有7個種,自然分布于歐亞及北美大陸的溫帶廣闊地域。其中,板栗(C.mollissima Blume)、錐栗(C.henryi (Skan) Rehd.et Wils.)和茅栗(C.seguinii Dode)為我國特有種。栗屬中國特有種在栗屬植物研究和利用中占有重要地位,其品質優良,適應性和抗逆性強,尤其對栗疫病和墨水病具有較強的抵抗力,是進行食用栗品種改良的重要基因來源,對世界栗屬資源保護和可持續利用具有重要作用[1]。我國栗屬育種研究起始于上世紀50—60年代,通過常規育種方法選育了一大批應用于生產的優良品種等[2-3],目前傳統的選擇育種和雜交育種仍是栗屬品種選育和性狀遺傳改良的主要方法。

栗屬為多年生木本植物,其童期長,選育一個新品種或對某一性狀進行遺傳改良周期漫長,至少需要10~15 a,且花費巨大。同時,栗屬農藝性狀大多是受微效多基因控制的數量性狀,遺傳基礎復雜,且易受環境影響,導致傳統的育種方法效率不高。確定控制農藝性狀的自然等位變異,不僅有助于解析其遺傳基礎,且可以為分子標記輔助育種乃至定向育種提供有效的基因資源和分子標記,具有重要的理論意義和應用價值[4]。近年來,關聯分析(Association analysis)成為解析復雜數量性狀遺傳基礎的一個重要手段。關聯分析又稱連鎖不平衡作圖(Linkage disequilibrium mapping)或關聯作圖(Association mapping),是一種以連鎖不平衡為基礎,檢測群體內處于連鎖不平衡狀態的標記或候選基因的遺傳變異與目標性狀顯著關聯頻率的方法,目前已廣泛應用于農林作物。截至目前,關聯分析在林木中應用主要以自然群體為主,而以雜交群體為材料進行關聯分析鮮有報道。

簡單重復序列(SSR)廣泛分布于各類真核生物基因組中,為共顯性標記,具有多態性高、重復性好、便于檢測等特點[5]。SSR標記現已成為生物遺傳特性研究的一種重要分子標記,被廣泛應用于栗屬植物遺傳多樣性分析、品種鑒定、遺傳圖譜構建等研究[5-11]。本課題組以錐栗種內、板栗和錐栗種間共9個雜交組合F1代群體(235個單株)為材料,前期采用32對高多態性栗屬SSR標記對其遺傳多樣性、遺傳效應等進行了分析[12]。在此基礎上,本研究進一步對其進行群體遺傳結構和連鎖不平衡分析,并將SSR標記與F1代生長、枝條、葉片表型及其光合生理等農藝性狀表型數據進行關聯分析,以期獲得與農藝性狀相關的分子標記,為栗屬植物分子標記輔助選擇和高效育種奠定基礎。

1 材料與方法

1.1 試驗材料

供試材料為錐栗種內、錐栗和板栗種間共9個雜交組合235份雜交F1代單株及其親本(表1)。其中,4個親本YLZ 1、YLZ 2、YLZ 24和YLZ 26為浙江省審(認)定錐栗良種,YLZ 24早熟、大果高產,YLZ 26特別高產穩產,YLZ 1和YLZ 2耐貯藏、優質且適宜加工;YLZ 14和YLZ 15為錐栗優株無性系,均是晚熟、優質品系;‘魁栗’為浙江省板栗良種,果型特大,平均單果質量20~25 g。2011年通過人工控制授粉獲得雜交種實,2012年播種育苗,2013年春季建立雜交子代測定林,同時栽植所有親本1年生嫁接苗若干株,株行距4 m×4 m,每年進行精細撫育管理,235份子代單株及其親本生長發育狀況良好。

表1 栗雜交組合概況Table 1 Cross combination of chestnut

1.2 樣品采集與DNA提取

于2019年4月份采集235個子代單株及其親本幼嫩、無病蟲害的新鮮葉片,每個單株采集嫩葉3~5片,裝于事先做好標記的自封袋中,并將自封袋置于裝有干冰的泡沫箱中帶回實驗室。

采用北京艾德萊生物科技公司生產的DNA試劑盒提取235份子代及其親本基因組DNA,所提取的DNA經1%的瓊脂糖凝膠電泳檢測質量和紫外分光光度計測定濃度后,?20 ℃低溫下保存備用。

1.3 SSR標記引物篩選與毛細管電泳分析

從120對引物中篩選出條帶清晰、多態性高、重復性好的32對SSR引物,采用毛細管電泳儀對所有樣本的SSR擴增產物進行分析(圖1)。32對SSR引物具體信息、SSR擴增程序、毛細管電泳數據獲得等內容詳見本課題組前期研究[12]。

圖1 SSR 標記毛細管電泳Fig.1 Capillary electrophoresis of SSR markers

1.4 表型性狀測定

分別于2018—2019年(連續2 a)1、7、11月份對235份子代單株的枝條性狀、葉片表型及其光合生理指標、生長性狀進行調查測定,具體方法參照本課題組章平生等[13-14]前期研究。

1.5 數據分析

利用TASSEL 5.0軟件和R語言完成標記間的連鎖不平衡(LD)分析及關聯分析,STRUCTURE 2.3.4 軟件進行群體遺傳結構分析。

根據R2確定標記間的關聯程度,即假設有兩個連鎖的座位A和B,其等位基因分別為A、a和B、b,用πA、πa、πB、πb分別表示相應的等位基因頻率,用πAB、πaB、πAb和πab分別表示4種單倍型AB、aB、Ab和ab的頻率,則實際觀測到的單倍型頻率與期望單倍型頻率間的差異D及R2的計算公式為:

Dab=πAB?πAπB

R2=(Dab)2/πAπaπBπb

應用STRUCTURE 2.3.4軟件估測栗雜交子代的群體遺傳結構。依據Evanno等[15]的統計模型,以最大似然法和Delta K值對235份雜交子代進行亞群劃分。數據錄入后,初始階段MCMC的不作數迭代默認設置為10 000次,再將不作數迭代后的MCMC調為100 000次,群體的數目(K)設為2~13,以此估計可能的群體數目。計算Q參數,作為關聯分析的協變量。

采用TASSEL 5.0軟件中的一般線性模型(General linear model,GLM)和混合線性模型(Mixed linear model,MLM)對SSR標記和農藝性狀進行關聯分析。GLM以Q作為協變量進行回歸分析;MLM 采用Q+K方法分析。

2 結果與分析

2.1 雜交F1代群體遺傳結構分析

群體遺傳結構是影響連鎖不平衡(LD)水平的一個重要因素,亞群數量的增加會使連鎖不平衡的程度增強,可能導致多態性基因位點與目標性狀間產生偽關聯,造成假陽性結果的出現,因此,在進行關聯分析時對群體遺傳結構的分析是不可或缺的步驟。本研究結果顯示,當K=4時,圖像出現明顯的拐點(圖2),Delta K值也接近最大。據此,235個雜交F1代最終被劃分為4個亞群(圖3),平均混合度為0.120。

圖2 不同K值下Delta K值的變化Fig.2 Change of Delta K value under different K values

圖3 雜交F1代群體遺傳結構Fig.3 Genetic structure of F1 hybrid population

由圖3可知,235份雜交子代群體被劃分為4個亞群,而并未按照其雜交組合分為9個亞群,但每個亞群中子代分布呈現出一定的遺傳分化。其中,第1亞群主要包括組合C8(‘魁栗’×YLZ 15號)和C9(‘魁栗’בYLZ 1號’),且有少部分組合C2、C6和C7的子代混入其中;第2亞群主要為組合C2(‘YLZ 26號’×YLZ 15號)、C4(‘YLZ 24 號’×YLZ 15 號)和C1(‘YLZ 26 號’×YLZ 14號),同樣有少部分組合C3、C5和C7的子代混入其中;第3亞群主要包括組合C3(‘YLZ 24號’בYLZ 1號’)和C5(‘YLZ 1號’בYLZ 24號’),少部分組合C6和C7的子代混入其中;而第4亞群則主要為組合C6(‘YLZ 1號’בYLZ 2號’)和C7(YLZ 14號בYLZ 1號’),組合C9的個別子代混入其中。此結果與前期研究中UPGMA和PCA分類結果基本一致[12]。

2.2 SSR標記間連鎖不平衡分析

利用TASSEL 5.0軟件對32對SSR標記進行連鎖不平衡(LD)分析,根據標記間R2來確定兩兩標記間的關聯程度,利用R語言ggplot和pheatmap軟件包對所得數據進行可視化分析(圖4)。可以看出,統計概率P<0.05時,32對SSR標記組成496個位點中存在一定的連鎖不平衡(圖4中對角線上面淡紅及紅色位點),其中,74個位點的連鎖不平衡水平較高,占總標記對的比例為14.92%,其余位點的連鎖不平衡水平相對較低;而P<0.01時,大多數位點間連鎖不平衡水平較低,且有37個標記對的R2值為0,表現出完全的連鎖平衡現象。LD分析結果表明32個SSR位點間連鎖不平衡總體水平較低。

圖4 32個 SSR 標記間 LD 的分布情況Fig.4 The distribution of LD among 32 SSR markers

2.3 農藝性狀與SSR標記的關聯分析

利用TASSEL 5.0軟件一般線性模型(GLM)和混合線性模型(MLM)對雜交子代涵蓋生長、枝條、葉片表型及其光合生理的25個農藝性狀(連續2 a數據平均值)與32個SSR標記進行關聯分析,并結合比較包法隆尼(Bonferroni threshold,P<1/1 048=9.5×10?4)閾值,對標記進行篩選,結果見表2。

表2 GML和MLM模型關聯標記及其對農藝性狀的解釋率Table 2 Explanation ratio of associated markers with agronomic traits in GML and MLM models %

·續表2·

GLM模型中,高達28個SSR標記位點與樹高等24個農藝性狀呈極顯著關聯,每個性狀可檢測到關聯的SSR位點數1~22個不等,解釋率范圍為4.65%~24.02%。不同性狀關聯的位點數不同,3個生長性狀中,與樹高相關聯的位點數最多,達16個,其中位點ICMA017s、PRD21關聯解釋率相對較高,分別能解釋12.06%、11.60%的樹高變異,而其它位點的解釋率均不足10.00%。3個枝條性狀中,關聯位點數最多的性狀為節間距,但其位點解釋率偏低,均在8.00%以下。在葉片表型及其光合生理性狀中,葉片寬度、葉面積及葉片干質量關聯的位點數最多,均達22個。其中,葉片寬度、葉面積均與位點CsCAT7表現出較高的關聯解釋率,分別為16.70%、18.73%;與葉片干質量關聯解釋率最高的位點是ICMA017s,能解釋15.50%的葉片干質量變異;部分葉片性狀關聯的位點較少,如類胡蘿卜素僅與位點CsCAT18相關聯,葉緣齒數僅與ICMA003相關聯,葉綠素a僅與CmTCR22和CsCAT18相關聯,葉脈分枝角僅與CsCAT26和ICMA010相關聯。此外,不同位點關聯的性狀數也不相同,其中位點CsCAT41僅能解釋4.67%的葉尖角變異,位點ICMA022僅能解釋6.57%的樹高變異,而位點CmTCR22與樹高、地徑等17個農藝性狀關聯,關聯的性狀最多,但解釋率均在10.00%以下。

MLM模型中,得出26個SSR位點與樹高等23個農藝性狀呈極顯著關聯,解釋率范圍為5.04%~24.02%。與GLM模型關聯分析的結果相比,部分性狀關聯的位點數以及部分位點關聯的性狀數均有所減少,但解釋率差異不明顯,其中位點CsCAT18不再與類胡蘿卜素和葉綠素a關聯,位點CsCAT26不再與葉脈分枝角關聯。

綜合兩種模型關聯分析結果,發現15個SSR標記與樹高、地徑、冠徑3個生長性狀關聯,14個SSR標記與1年生枝條長度、1年生枝條粗度、節間距3個枝條性狀關聯,26個SSR標記與葉片長度等18個葉片表型及其光合生理性狀關聯。進一步對各性狀高度關聯的SSR標記進行篩選分析,最終得出ICMA017s等15個標記分別與樹高等13個農藝性狀高度關聯,標記對性狀的解釋率均在10.00%以上(表3)。其中,ICMA017s與3個生長性狀高度關聯,PRD21與樹高性狀高度關聯,表明這2個SSR標記與栗樹體生長高度相關;PRD26、CmTCR4分別與1年生枝條粗度、1年生枝條長度高度關聯,表明這2個標記與栗枝條生長發育高度相關;PRA86等15個標記分別與代表葉片大小、形狀、厚度以及質量的8個表型性狀高度關聯,表明這些標記與栗葉片生長發育高度相關,同時,葉片各性狀高度關聯的標記數量差異較大,其中葉片干質量高度關聯的標記數量達10個,其次為葉片寬度和葉面積,均為9個,而葉形指數關聯的標記數量最少,僅2個。此外,還存在同一標記同時與多個性狀高度關聯的現象,如ICMA017s與樹高、地徑、冠徑、葉片寬度、葉面積、葉片鮮質量及葉片干質量均存在高度關聯,PRD21與樹高、葉片長度、葉片寬度、葉形指數均高度關聯,CmTCR4與1年生枝條長度、葉片厚度、葉片鮮質量、葉片干質量及比葉質量均高度關聯,表明某些標記位點可能同時控制多個性狀,且某一性狀同時受多個標記位點控制。

表3 性狀高度關聯SSR標記統計Table 3 The traits with high associated SSR markers

3 討論

農藝性狀大多是受多基因控制的數量性狀,且各性狀間存在一定的關聯,通過常規育種技術對農藝性狀進行遺傳改良往往效果不顯著,且周期長、效率低。在DNA水平利用分子標記對農藝性狀進行早期選擇和輔助育種前景廣闊。本研究以涵蓋9個雜交組合的235份栗雜交子代混合群體為材料,在進行群體遺傳結構和連鎖不平衡(LD)分析的基礎上,利用32個高多態性SSR標記與25個農藝性狀進行關聯分析,以獲得與農藝性狀相關聯的SSR標記位點,從而為栗屬植物分子標記輔助育種特別是雜交子代的選擇和高效育種奠定基礎。

連鎖不平衡也稱為配子相不平衡(Gametic phase disequilibrium)、配子不平衡(Gametic disequilibrium)或等位基因關聯(Allelic association),是指群體內不同座位等位基因(可以是標記,也可是基因/QTL間與標記)間的非隨機關聯[16]。LD在染色體上的分布一般用LD衰減散點圖和LD配對檢測的矩陣圖來描述,前者可以觀測LD隨遺傳或物理距離的衰減速率,后者可以直接觀測同一染色體的基因座或基因的多態性位點間LD的線性排列[17]。本研究LD分析采用了后者,結果較為直觀,統計概率P<0.05時,32對SSR標記組成496個位點中存在一定的連鎖不平衡,其中,74個位點的連鎖不平衡水平較高,占總標記對的14.92%,其余位點的連鎖不平衡水平相對較低;而P<0.01時,32個SSR位點間連鎖不平衡總體水平較低。LD水平決定了關聯分析的解析精度,是開展關聯分析研究的理論基礎。LD的一個明顯特性是群體依賴性,選擇的群體不同,其LD水平顯著不同。當所用群體來源有限時,其LD將維持在一個較高的水平。而選用多樣性較高的群體則包括更多不同來源的研究個體,其LD水平一般較低;同時,群體混合可通過引進不同祖先來源和等位基因頻率的染色體而影響群體的LD水平[18]。本研究所使用的材料為涵蓋9個雜交組合235份雜交子代的混合群體,其親本涉及錐栗和板栗2個物種的7個優良品種或優株無性系,并通過雜交即基因的交換和重組產生了大量的遺傳變異類型,來源較為復雜多樣,且前期研究表明不同組合子代群體具有豐富的遺傳多樣性(Shannon’s多樣性指數范圍0.881 6~1.131 7)[12],其LD水平較低。

群體結構的存在和亞群內等位基因頻率的不均等分布將導致多態性位點和表型的假陽性關聯,分析群體遺傳結構對表型性狀的影響,對于防止假陽性關聯是非常必要的[19]。因此,群體結構分析是關聯分析的前提。本研究所用的栗不同雜交組合群體并非自然群體,在一定程度上不同組合間存在一定的親緣關系,而對供試群體進行結構分析可將群體分為適宜的幾個亞群,將群體結構信息作為協變量來進行關聯分析,可以提高結果的準確性,在一定程度上使得假陽性得到控制。群體遺傳結構分析表明,當K=4時,Delta K值接近最大。據此,235個雜交F1代群體最終被劃分為4個亞群,且每個亞群中子代分布也呈現出一定的遺傳分化,平均混合度為0.120。此結果與前期研究中UPGMA和PCA分類結果基本一致。

關聯分析又稱連鎖不平衡作圖,是一種以基因間連鎖不平衡為基礎,鑒定某一群體內目標性狀與標記間關系的分析方法。與連鎖作圖分析相比,關聯分析所用群體具有更為廣泛的遺傳基礎,可同時對同一基因座的多個等位基因進行分析,而絕大部分常規QTL作圖所用群體通常為2親本雜交重組后代,其基因座一般只涉及2個等位基因;同時,關聯分析作圖定位更為準確和具有更高的分辨率,可實現對QTL的精細定位,甚至直接定位到基因本身[20]。因此,關聯分析與QTL作圖分析既互為補充,又比連鎖作圖分析有著更為廣泛的應用。曹珂等[21]以分布于桃(Prunus persica L.)8個連鎖群上的53個SSR標記與104份桃地方品種單果質量及物候期性狀表型數據進行關聯分析,得到27個與桃單果質量及6個物候期性狀關聯的QTLs。作者等[11]以來自山東等10個省份95個板栗地方品種為試材,將17對SSR標記與板栗地方品種農藝性狀進行關聯分析,經假陽性檢驗發現葉柄長度和淀粉含量分別與標記CsCAT 5和CsCAT 22呈顯著關聯。張琳等[22]以篩選出的16對高多態性SSR標記對36份產油量較好、觀賞性較佳的牡丹(Paeonia suffruticosa Andr.)材料進行關聯分析,發現3個SSR位點與單株種子產量、聚合蓇葖果數等9個性狀極顯著關聯(P<0.01)。關聯分析所使用材料既可以是自然群體,也可以是多親本雜交混合群體。如對測交群體進行關聯分析,可以直接定位F1雜種當代控制目標性狀的遺傳位點或區域,進而直接反映出F1雜種當代的遺傳效應,是對傳統連鎖分析和關聯分析群體定位結果的有益補充[23]。姚俊修從16個鵝掌楸(Liriodendron chinense (Hemsl.) Sarg.)種間雜交組合中抽取430份子代個體,在測定其樹高和胸徑兩個表型性狀的基礎上,利用101個SSR標記檢測每個個體的基因型;然后利用Tassel 3.0軟件分析不同標記間的連鎖不平衡(LD值),Structure 2.3軟件對子代群體進行分群;在群體結構矯正基礎上,將表型與分子標記進行關聯分析,發現6個標記與樹高性狀相關聯,7個標記與胸徑性狀關聯[24]。本研究中,綜合兩種模型關聯分析結果,發現15個SSR標記與樹高等3個生長性狀關聯,14個SSR標記與1年生枝條長度等3個枝條性狀關聯,26個SSR標記與葉片長度等18個葉片表型及其光合生理性狀關聯。進一步對各性狀高度關聯的SSR標記進行篩選分析,最終得出ICMA017s等15個SSR標記分別與樹高等13個農藝性狀高度關聯,標記對性狀的解釋率均在10.00%以上。

同時,本研究發現存在同一標記同時與多個性狀高度關聯的現象,如ICMA017s與樹高、地徑、冠徑、葉片寬度、葉面積、葉片鮮質量及葉片干質量7個性狀均存在高度關聯,CmTCR4與1年生枝條長度、葉片厚度、葉片鮮質量、葉片干質量及比葉質量等5個性狀均高度關聯,PRD21與樹高、葉片長度、葉片寬度、葉形指數4個性狀均高度關聯;還存在同一性狀與多個標記高度關聯的現象,如葉片干質量高度關聯的SSR標記數量多達10個,其次為葉片寬度和葉面積,均為9個,葉片鮮質量、葉片長度、比葉質量和樹高性狀高度關聯的標記分別為5、3、3和2個。表明某些標記位點可能同時控制多個性狀,且某一性狀同時受多個標記位點控制。此現象在柳樹(Salix babylonica L.)[25]、日本落葉松(Larix kaempferi(Lamb.)Carr.)[26]、豇豆(Vigna unguiculata (L.)Walp.)[27]的研究中也有發現。說明植物中存在一個基因控制多個農藝性狀或者影響這些農藝性狀的部分基因存在緊密連鎖的現象,從而導致標記位點變異具有多效性[28]。據此推測生長、枝條、葉片表型及其光合生理等農藝性狀的這種相關性可能是由于基因的一因多效表達或者是連鎖不平衡產生的。

4 結論

供試雜交群體32個SSR位點間LD總體水平較低;根據Delta K值計算結果,235個栗雜交F1代混合群體被劃分為4個亞群,且每個亞群中子代分布呈現出一定的遺傳分化,平均混合度為0.120;綜合GLM和MLM兩種模型關聯分析結果,得出15個SSR標記與樹高等3個生長性狀關聯,14個SSR標記與1年生枝條長度等3個枝條性狀關聯,26個SSR標記與葉片長度等18個葉片表型及其光合生理性狀關聯。其中,ICMA017s等15個SSR標記分別與樹高等13個農藝性狀高度關聯,標記對性狀的解釋率均在10.00%以上,且存在同一標記與多個性狀、同一性狀與多個標記高度關聯的現象。研究結果為栗屬植物分子標記輔助育種特別是雜交子代的選擇和高效育種奠定了基礎。

猜你喜歡
關聯分析
不懼于新,不困于形——一道函數“關聯”題的剖析與拓展
“苦”的關聯
當代陜西(2021年17期)2021-11-06 03:21:36
隱蔽失效適航要求符合性驗證分析
“一帶一路”遞進,關聯民生更緊
當代陜西(2019年15期)2019-09-02 01:52:00
電力系統不平衡分析
電子制作(2018年18期)2018-11-14 01:48:24
奇趣搭配
智趣
讀者(2017年5期)2017-02-15 18:04:18
電力系統及其自動化發展趨勢分析
中西醫結合治療抑郁癥100例分析
在線教育與MOOC的比較分析
主站蜘蛛池模板: 婷婷六月激情综合一区| 99视频只有精品| 2020极品精品国产 | 55夜色66夜色国产精品视频| 欧美成人在线免费| 欧美.成人.综合在线 | 一本一本大道香蕉久在线播放| 国产成人综合在线观看| 亚洲欧美自拍中文| 国产高清不卡视频| 国产日韩精品欧美一区灰| 国产香蕉一区二区在线网站| 亚洲欧美日韩另类在线一| 99热亚洲精品6码| 国产xxxxx免费视频| 欧美综合激情| 午夜精品影院| 午夜天堂视频| 免费观看亚洲人成网站| 欧美精品啪啪| 91在线日韩在线播放| 国产打屁股免费区网站| 亚洲人成网18禁| 狠狠ⅴ日韩v欧美v天堂| 日本国产在线| 日韩在线视频网| 国产另类视频| 国产一区二区精品高清在线观看| 一本一道波多野结衣一区二区| 免费欧美一级| 全部免费特黄特色大片视频| 国产精品短篇二区| 国产一级裸网站| AV无码一区二区三区四区| 无码国产偷倩在线播放老年人 | 五月婷婷导航| 91精品最新国内在线播放| 国产精品成人观看视频国产| 亚洲av无码人妻| 国产精品人成在线播放| 国产迷奸在线看| 亚洲第一成年网| 国产原创第一页在线观看| 巨熟乳波霸若妻中文观看免费| 99re这里只有国产中文精品国产精品| 国产情侣一区| 亚洲男女在线| 午夜免费小视频| 国产网站免费观看| 国产综合另类小说色区色噜噜| 亚洲欧美在线精品一区二区| 波多野结衣一区二区三区四区| 免费xxxxx在线观看网站| 国产人成在线观看| 久久人体视频| 一级毛片不卡片免费观看| 午夜电影在线观看国产1区| 一区二区午夜| 亚洲人成人无码www| 精品三级网站| 青青青草国产| 凹凸国产分类在线观看| 国产Av无码精品色午夜| 91精品专区国产盗摄| 欧美成人国产| 国产欧美视频在线| av手机版在线播放| 91精品国产91久无码网站| 最新精品久久精品| A级毛片无码久久精品免费| 国产精品一老牛影视频| 欧美午夜理伦三级在线观看| 五月激情婷婷综合| 国产成人资源| 亚洲精品麻豆| 国产天天色| 99热这里只有精品5| 在线不卡免费视频| 亚洲第一区在线| 91亚洲精选| 自拍欧美亚洲| 男女精品视频|