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

虹鱒傳染性胰臟壞死性狀全基因組關(guān)聯(lián)分析

2022-10-19 06:35:36歐陽少琪趙云峰蔣麗楊潤清
關(guān)鍵詞:關(guān)聯(lián)分析

歐陽少琪,趙云峰,蔣麗*,楊潤清*

虹鱒傳染性胰臟壞死性狀全基因組關(guān)聯(lián)分析

歐陽少琪1,2,趙云峰2,蔣麗2*,楊潤清2*

1. 上海海洋大學, 水產(chǎn)科學國家級實驗教學示范中心, 上海 201306 2. 中國水產(chǎn)科學研究院, 農(nóng)業(yè)農(nóng)村部水生動物基因組學重點實驗室, 北京市漁業(yè)生物技術(shù)重點實驗室, 北京 100141

為了探究虹鱒傳染性胰臟壞死癥(Infectious Pancreatic Necrosis, IPN)抗性的遺傳基礎,本文采用三種二歧性狀的分析方法,對來自58個全同胞家族的749尾虹鱒(271尾抗性和478尾易感),注射傳染性胰臟壞死病毒后的存活情況進行全基因組關(guān)聯(lián)分析(Genome-Wide Association Study, GWAS)。關(guān)聯(lián)分析結(jié)果表明,共鑒定到兩個QTN與虹鱒IPN抗性顯著關(guān)聯(lián),位于四號染色體和對應的Scaffold上。以這些檢測到的QTNs作為探針,對比虹鱒全基因組信息,篩選得到四個候選基因,它們可能是影響虹鱒IPN抗性的重要功能基因。通過對比三種方法的檢測結(jié)果,得知GRAMMAR-Lambda法對檢測二歧性狀QTNs的統(tǒng)計能力略高于SAIGE和GMMAT法,具有更好的統(tǒng)計性能。這些研究結(jié)果為虹鱒抗傳染性胰臟壞死病毒的遺傳選育提供了分子標記,為虹鱒抗傳染性胰臟壞死病毒的研究提供了理論基礎。

虹鱒; 胰臟; 全基因組關(guān)聯(lián)分析

虹鱒系鮭形目鮭科類大馬哈魚屬的重要代表性魚類,是一種具有極高經(jīng)濟價值的世界性冷水養(yǎng)殖魚類。由傳染性胰臟壞死病毒(Infectious Pancreatic Necrosis Virus, IPNV)引起的傳染性胰臟壞死癥(Infectious Pancreatic Necrosis, IPN)最先于北美,后傳入歐洲。1984年,IPNV首次在我國出現(xiàn)。隨后,在我國的黑龍江、遼寧、山西、甘肅、臺灣等地陸續(xù)檢測出同樣的病原體。2016年,四川省某虹鱒養(yǎng)殖場因IPVN感染導致虹鱒大量死亡,這是在我國西南地區(qū)發(fā)現(xiàn)的第一例由IPNV引起的魚病[1]。時至今日,IPNV已成為制約大西洋鮭魚和虹鱒魚養(yǎng)殖業(yè)發(fā)展的主要病害之一,給我國的水產(chǎn)業(yè)造成了嚴重的損失。當感染IPNV后,魚體是否死亡受品種、年齡、環(huán)境、遺傳背景等的影響[2],現(xiàn)有的研究表明,在鮭科類中,可能存在對細菌[3,4]、寄生蟲[5]和病毒性疾病[6]免疫相關(guān)的顯著遺傳變異。這意味著有可能通過人工選育的方式,來提高魚類對各種疾病的抵抗力,從而作為魚類養(yǎng)殖的一種可行性疾病控制策略。傳統(tǒng)育種一般采用對群體中大量個體進行性狀觀測的方法來選育優(yōu)良后代,這一過程不僅耗時長、操作繁瑣,而且效率低?;谌蚪M關(guān)聯(lián)分析(Genome-Wide Association Study, GWAS)的分子輔助標記育種手段,以全基因組所有單核苷酸多態(tài)性為分子遺傳標記,進行相關(guān)性分析,篩選出與性狀相關(guān)聯(lián)的數(shù)量性狀核苷酸(Quantitative Trait Nucleotide, QTN),精確定位主效基因,從而極大縮短了選育時間。

隨著近年來算法的發(fā)展,線性混合模型(Linear mixed model, LMM)已成為GWAS廣泛采用的模型之一,用于校正群體分層和親緣相關(guān)性帶來的干擾。LMM假設性狀具有恒定的殘差方差,而二歧性狀不符合該假設,故LMM并不適用于分析二歧性狀。廣義線性混合模型(Generalized Linear Mixed Model, GLMM),考慮隨機多基因效應,提高了疾病性狀的定位能力。近年來,出現(xiàn)了一些針對二歧性狀的方法,試圖解決某些由基因型類別或低頻率事件導致的計算和分離問題。Chen H等[7]提出了GMMAT(Generalized Mixed Model Association Test),使用logistic混合模型,并拓展了score檢驗,比極大似然法更快。在GLMM的基礎上,Zhou X等[8]提出了SAIGE(Scalable and Accurate Implementation of Generalized Mixed Model),將BOLT-LMM正態(tài)分布性狀擴展到了二元疾病性狀。Yang R等[9]將預先估計的基因組育種值視為logit回歸中的已知預測因子,然后使用基因組控制值來校正關(guān)聯(lián)統(tǒng)計檢驗量,從而將GRAMMAR-Lambda法拓展到復雜群體的二元疾病分析。

目前,已經(jīng)在鮭魚和鱒魚中找到了與各種疾病抗性相關(guān)的QTL。根據(jù)兩起IPVN疫情和攻毒實驗結(jié)果,在大西洋鮭魚基因組中存在著一個抗IPVN的主效QTL[10,11]。隨后,在大西洋鮭魚中還發(fā)現(xiàn)了一個可以與病毒結(jié)合并通過了標記輔助選擇的抗IPNV候選基因[12]。已有研究表明,虹鱒種群對IPNV的抗性與收獲時魚的體重沒有明確相關(guān)性。

本研究中,將虹鱒在攻毒處理后是否存活作為二歧性狀,采用GMMAT、SAIGE、GRAMMAR-Lambda法進行全基因關(guān)聯(lián)分析,挖掘與IPNV抗性相關(guān)的候選基因,為虹鱒IPNV抗性遺傳改良提供參考。

1 材料與方法

1.1 實驗群體與表型鑒定

本實驗虹鱒由ATC 巴塔哥尼亞研究中心(智利,蒙特港)提供的58條雌魚和20條雄魚繁殖而來,包括58個全同胞家系。相關(guān)飼養(yǎng)條件和更多群體管理請參閱Flores-Mara[13],Neto[14]和Rodríguez[15]等人的文章,有關(guān)攻毒實驗的詳細分步方案參見Yoshida等[16]人的文章。實驗期間死亡個體的表型值記為1,實驗結(jié)束時仍然存活的個體表型值記為0。

1.2 基因型與質(zhì)量控制

使用Palti等[17]開發(fā)的57K Affymetrix Axiom SNP芯片對749尾虹鱒的基因型進行分型,得到40150個單核苷酸多態(tài)性(Single Nucleotide Polymorphisms, SNP)。使用PLINK1.9軟件(http://www.cog-genomics.org/plink2/)對基因型數(shù)據(jù)進行嚴格的質(zhì)量控制,去除低于90%最小檢出頻率的個體,以及最小哈代溫伯格平衡檢驗值小于1×10-5、最小等位基因頻率小于0.01、最小檢出頻率低于95%的SNPs。

1.3 全基因關(guān)聯(lián)分析

其中是虹鱒在實驗期間死亡率的均值,為標記和群體結(jié)構(gòu)以外的固定效應值矩陣,包括有年齡、性別和當前檢驗SNP的遺傳效應,是扣除當前檢驗SNP后的隨機多基因效應值行向量,和分別是和的指示變量矩陣。

本研究利用GRAMMAR-Lambda v1.0(URL:https://github.com/RunKingProgram/BinaryGRAMMAR-Lambda/),以實驗期間虹鱒攻毒處理后是否存活為二歧性狀,進行全基因組關(guān)聯(lián)分析,設置基因組顯著水平(5% Bonferroni校正閾值,1.325697×10-6)。簡單來說,GRAMMAR-Lambda處理二歧性狀時可大致分為三步,第一步,通過全基因組標記來計算親緣關(guān)系矩陣;第二步,構(gòu)建GBLUP方程估計基因組育種值,而不是多基因效應,然后將基因組育種值作為偏移項構(gòu)建線性模型;第三步,執(zhí)行PLINK并使用檢驗統(tǒng)計量的平均值或中位數(shù)校正統(tǒng)計量。同時,使用SAIGE v0.29.5(URL:https://github.com/weizhouUMICH/SAIGE/)和GMMAT v1.3.2(URL:https://github.com/hanchenphd/GMMAT/)作為對比。結(jié)果圖由R3.6.1繪制而成。

1.4 數(shù)據(jù)來源

表型數(shù)據(jù)和基因型數(shù)據(jù)由在線數(shù)據(jù)庫Figshare(URL:https://figshare.com/articles/Untitled_Item/7725668/)提供。

2 結(jié)果與分析

2.1 表型和基因型數(shù)據(jù)分析

來自58個家系的749尾虹鱒,平均個體數(shù)在10~18不等,實驗期間,死亡數(shù)共計271,存活數(shù)為478,累積死亡率達到36.18%,標準差為0.48,標準誤為0.02。實驗結(jié)束時參與實驗的所有虹鱒的平均體重為2.17,標準差為0.74,標準誤為0.03。共有749個樣本的37716個SNP通過了質(zhì)量控制。每個染色體上的SNP平均有1257個,介于591和1796之間。SNP密度分布見圖1。

圖 1 虹鱒IPN性狀的SNP密度圖

Fig.1 The SNP density map of IPN trait in

2.2 GWAS結(jié)果分析

在對GWAS結(jié)果的優(yōu)劣進行評估時,通常根據(jù)QQ圖和基因組控制(Genomic Control, GC)值或膨脹系數(shù)來判斷分析結(jié)果的群體分層和假陽性影響。圖2顯示,749個個體,37716個SNP虹鱒IPNV存活數(shù)據(jù)進行GMMAT(a)、SAIGE(b)、GRAMMR-Lambda(c)分析得到的基因定位結(jié)果。比較不同分析結(jié)果下的QQ圖,可知觀測值和期望值均基本一致,散點分布與趨勢線吻合度較高,絕大多數(shù)位點的實際觀測P值貼近理論線,分別計算三種方法下的GC值(每個標記的卡方統(tǒng)計量的平均數(shù)),GMMAT和GRAMMAR-Lambda的GC值分別為0.996和1,表明分析結(jié)果的群體分層和假陽性影響得到了有效的控制,SAIGE的GC值為1.103,檢測結(jié)果有顯著的假陽性。比較不同方法的曼哈頓圖,GRAMMAR-Lambda超過基因組顯著水平的QTNs有2個,位于14號染色體和Scaffold上(表1),GMMAT和SAIGE沒有檢測到超過基因組顯著水平的QTNs。

圖 2 虹鱒IPNV抗性性狀GWAS 曼哈頓和QQ圖

(a)GMMAT; (b)SAIGE; (c)GRAMMAR-Lambda

表 1 虹鱒IPN抗性顯著關(guān)聯(lián)的QTNs信息

注: *表示未在虹鱒全基因組中進行功能注釋的新基因。

Note: * Represent novel gene which has not found any gene function annotation in genome of.

3 討 論

水產(chǎn)動物相關(guān)病害都是多基因控制的復雜性狀,探索與疾病相關(guān)聯(lián)的QTNs、候選基因,對解析水產(chǎn)動物的疾病遺傳機制,利用分子輔助育種提高遺傳選育效率有著重要意義。隨著二代測序技術(shù)的快速發(fā)展,測序的成本也相應的降低,有關(guān)水產(chǎn)動物疾病性狀的GWAS研究逐漸增多。盡管GLMM對離散性狀具有可解釋和可預測的特性,但由于求解和計算上的困難,無法有效應用于復雜疾病的GWAS。基于GLMM的關(guān)聯(lián)方法,GMMAT和SAIGE,在一定程度上簡化了對疾病性狀的全基因組混合模型關(guān)聯(lián)分析,但它們更適用于處理親緣關(guān)系簡單的群體,如人類。當研究群體的親緣關(guān)系復雜,估計出的遺傳力并不準確,此時更適合使用GRAMMAR-Lambda法。GRAMMAR-Lambda法在執(zhí)行二歧性狀的關(guān)聯(lián)分析時,并不需要估計遺傳力,在處理大樣本群體時,還可通過使用少量標記來計算GRM達到全標記的效果,極大的簡化了計算復雜度。

此前在大西洋鮭魚中找到的兩個IPNV抗性的QTL,其中最顯著的一個位于26號染色體上,解釋了29%的表型變異和83%的遺傳變異[10,18],該QTL在挪威的一個獨立的大西洋鮭種群中得以證實[11]。2001年,在虹鱒中發(fā)現(xiàn)2個與IPNV抗性相關(guān)的QTL,分別位于3號和22號連鎖群中,每個QTL可以解釋約17%的表型變異[19],隨后的一項研究也證實存在著同樣的QTL[20],同時還發(fā)現(xiàn)了在12號連鎖群中存在另一個顯著的QTL。根據(jù)連鎖圖譜得知,這些重要標記分別位于14號染色體和16號染色體上[21]。

根據(jù)GWAS分析結(jié)果,在虹鱒全基因組上尋找高度關(guān)聯(lián)的QTN附近的基因,共找到四個已知的候選基因。其中基因在神經(jīng)元發(fā)育過程中起主要調(diào)節(jié)作用,此外還在介導細胞黏附和Gaq蛋白偶聯(lián)受體誘導的信號通路中發(fā)揮著重要作用[22];是一種印跡基因,其編碼的Grb10蛋白,通過拮抗胰島素信號,在哺乳動物細細胞培養(yǎng)中起到生長抑制因子的作用[23];在以斑馬魚為實驗動物模型時,發(fā)現(xiàn)基因在調(diào)節(jié)腺酸酯環(huán)酶和磷氧酶的活性過程中起著主要作用,同時與ATP二磷酸裂解酶(環(huán)化)、cAMP生成肽活性有著密切關(guān)聯(lián)[24];是一種核苷酸末端轉(zhuǎn)移酶基因,在mRNA鳥苷酸化過程中,起保護mRNA免受快速去腺苷酸化的作用,此外還參與調(diào)節(jié)NTP聚合酶、聚腺苷酸合成酶、聚腺苷酸聚合酶的活性,在生物體轉(zhuǎn)錄過程中起重要作用[25]。這些新位點的發(fā)現(xiàn)對后續(xù)虹鱒傳染性胰臟壞死癥抗性基因的發(fā)掘以及進一步的功能驗證、分子輔助選育工作具有重要意義,同時定位到的這些候選基因可能在虹鱒參與抗IPNV過程中發(fā)揮著重要作用,它們的篩選定位極大地提高了對虹鱒抗IPNV疾病遺傳基礎的認識?;谝陨涎芯拷Y(jié)果,下一步研究工作將主要集中在虹鱒IPNV抗性基因的克隆與功能驗證,為虹鱒IPNV抗性遺傳改良提供理論基礎。

4 結(jié) 論

GRAMMAR-Lambda是一種更加實用且高效的二歧性狀關(guān)聯(lián)分析方法,相比于傳統(tǒng)的二歧性狀分析方法,如GMMAT法和SAIGE法,檢測效率更高,計算速度更快,特別是針對親緣關(guān)系復雜群體,無法估計遺傳力或者估計遺傳力有誤時,有更好的效果。在虹鱒中,共鑒定到了四個候選基因,位于14號染色體上,這些候選基因可能在虹鱒參與抗IPNV過程中發(fā)揮著重要作用。

[1] Zhu L, Wang X, Wang K,Outbreak of infectious pancreatic necrosis virus (IPNV) in farmed rainbow trout in China [J]. Acta tropica, 2017,170:63-69

[2] Dorson M, Touchy C. The influence of fish age and water temperature on mortalities of rainbow trout, Salmo gairdneri Richardson, caused by a European strain of infectious pancreatic necrosis virus [J]. Journal of Fish Diseases, 1981,4(3):213-221

[3] Leeds TD, Silverstein JT, Weber GM,Response to selection for bacterial cold water disease resistance in rainbow trout [J]. Journal of Animal Science, 2010,88(6):1936-1946

[4] Silverstein JT, Vallejo RL, Palti Y,Rainbow trout resistance to bacterial cold-water disease is moderately heritable and is not adversely correlated with growth [J]. Journal of Animal Science, 2009,87(3):860-867

[5] Lhorente JP, Gallardo JA, Villanueva B,Quantitative genetic basis for resistance tosea lice in a breeding population of() [J]. Aquaculture, 2012,324:55-59

[6] Guy DR, Bishop SC, Brotherstone S,Analysis of the incidence of infectious pancreatic necrosis mortality in pedigreed,r L., populations [J]. Journal of Fish Diseases, 2006,29(11):637-647

[7] Chen H, Wang C, Conomos MP,Control for population structure and relatedness for binary traits in genetic association studies via logistic mixed models [J]. Am J Hum Genet, 2016,98(4):653-666

[8] Zhou X, Stephens M. Genome-wide efficient mixed-model analysis for association studies [J]. Nat Genet, 2012,44(7):821-824

[9] Yang R, Bao J, Song Y,A Shortcut Approach for Large-scale Mixed Model Associations with Binary Traits [J/OL]. Research Square. [2021-03-26].https://doi.org/10.21203/rs.3.rs-312421/v1

[10] Houston RD, Haley CS, Hamilton A,Major quantitative trait loci affect resistance to infectious pancreatic necrosis in() [J]. Genetics, 2008,178(2):1109-1115

[11] Moen T, Baranski M, Sonesson AK,Confirmation and fine-mapping of a major QTL for resistance to infectious pancreatic necrosis in(): population-level associations between markers and trait [J]. J Bmc Genomics, 2009,10(1):1-14

[12] Moen T, Torgersen J, Santi N,Epithelial cadherin determines resistance to infectious pancreatic necrosis virus in[J]. Genetics, 2015,200(4):1313-1326

[13] Flores-Mara R, Rodríguez FH, Bangera R,Resistance against infectious pancreatic necrosis exhibits significant genetic variation and is not genetically correlated with harvest weight in rainbow trout () [J]. Aquaculture, 2017,479:155-160

[14] Reis Neto RV, Yoshida GM, Lhorente JP,Genome-wide association analysis for body weight identifies candidate genes related to development and metabolism in rainbow trout () [J]. Molecular Genetics and Genomics, 2019,294(3):563-571

[15] Rodríguez FH, Cáceres G, Lhorente JP,Genetic (co)variation in skin pigmentation patterns and growth in rainbow trout [J]. Animals, 2019,13(4):675-682

[16] Yoshida GM, Bangera R, Carvalheiro R,Genomic prediction accuracy for resistance againstin farmed rainbow trout [J]. G3: Genes Genomes Genetics, 2018,8(2):719-726

[17] Palti Y, Gao G, Liu S,The development and characterization of a 57K single nucleotide polymorphism array for rainbow trout [J]. Molecular Ecology Resources, 2015,15(3):662-672

[18] Houston RD, Davey JW, Bishop SC,Characterisation of QTL-linked and genome-wide restriction site-associated DNA (RAD) markers in farmed[J]. BMC Genomics, 2012,13(1):1-15

[19] Ozaki A, Sakamoto T, Khoo S,Quantitative trait loci (QTLs) associated with resistance/susceptibility to infectious pancreatic necrosis virus (IPNV) in rainbow trout () [J]. Molecular Genetics and Genomics, 2001,265(1):23-31

[20] Ozaki A, Khoo S. K, Yoshiura Y,Identification of additional Quantitative Trait Loci (QTL) responsible for susceptibility to infectious pancreatic necrosis virus in rainbow trout [J]. Fish Pathology, 2007,42(3):131-140

[21] Hu ZL, Park CA, Reecy JM. Development progress and current status of the animal QTLdb [J]. Nucleic Acids Research, 2016,44(D1):D827-D833

[22] Schmidt S, Debant A. Function and regulation of the Rho guanine nucleotide exchange factor Trio [J]. Small GTPases, 2014,5(4):e983880

[23] Vecchione A, Marchese A, Henry P,The Grb10/Nedd4 complex regulates ligand-induced ubiquitination and stability of the insulin-like growth factor 1 receptor[J]. Molecular and Cellular Biology, 2003,23(9):3363-3372

[24] Zheng Y, Yuan JL, Meng SL,Testicular transcriptome alterations in zebrafish () exposure to 17β- estradiol [J]. Chenosphere, 2019,218:14-25

[25] Lim J, Kim D, Lee Young Suk,Mixed tailing by TENT4A and TENT4B shields mRNA from rapid deadenylation [J]. Science, 2018,361(6403):701-704

Genome-wide Association Analysis of Infections Pancreatic Necrosis Traits in

OU-YANG Shao-qi1,2, ZHAO Yun-feng2, JIANG Li2*, YANG Run-qing2*

1.,201306,2.,100141,

In order to explore the genetic basis of resistance to infectious pancreatic necrosis (IPN) in rainbow trout, a genome wide association study (GWAS) was conducted on 749 rainbow trout (271 resistant and 478 susceptible) from 58 full sib families after injection of infectious pancreatic necrosis virus (IPN) by using three methods for analyzing binary traits. The results showed that two SNP markers were significantly correlated with the anti-infectious pancreatic necrosis virus of, located on chromosomes 14 and the corresponding Scaffolds, respectively. Using these detected QTNs as probes, we compared the whole genome information of, and four candidate genes nearby were found, which may be functional genes affecting the resistance to IPNV in. By comparing the three methods, GRAMMAR-Lambda has a slightly higher statistical ability than SAIGE and GMMAT in detecting QTNs of binary traits and has better statistical performance. These results provide molecular markers for genetic selection of anti-infectious pancreatic necrosis virus inand the molecular basis for further research of the mechanism of infectious pancreatic necrosis virus in.

; pancreas; genome wide association analysis

S917.4

A

1000-2324(2022)04-0618-06

10.3969/j.issn.1000-2324.2022.04.018

2022-04-11

2022-04-20

動物全基因混合模型關(guān)聯(lián)分析多功能優(yōu)化解析策略(32072726)

歐陽少琪(1997-),男,碩士研究生,研究方向:數(shù)量遺傳學. E-mail:ouyangsq7@163.com

Author for correspondence. E-mail:runqingyang@cafs.ac.cn; jiangl@cafs.ac.cn.

猜你喜歡
關(guān)聯(lián)分析
不懼于新,不困于形——一道函數(shù)“關(guān)聯(lián)”題的剖析與拓展
“苦”的關(guān)聯(lián)
當代陜西(2021年17期)2021-11-06 03:21:36
隱蔽失效適航要求符合性驗證分析
“一帶一路”遞進,關(guān)聯(lián)民生更緊
當代陜西(2019年15期)2019-09-02 01:52:00
電力系統(tǒng)不平衡分析
電子制作(2018年18期)2018-11-14 01:48:24
奇趣搭配
智趣
讀者(2017年5期)2017-02-15 18:04:18
電力系統(tǒng)及其自動化發(fā)展趨勢分析
中西醫(yī)結(jié)合治療抑郁癥100例分析
在線教育與MOOC的比較分析
主站蜘蛛池模板: 国产欧美日韩免费| 午夜国产精品视频黄| 女人爽到高潮免费视频大全| 欧美A级V片在线观看| 久久国产成人精品国产成人亚洲 | 中文字幕第4页| 亚洲国产91人成在线| 久久窝窝国产精品午夜看片| 国产成人麻豆精品| 九色综合伊人久久富二代| 黄色网址免费在线| 人与鲁专区| 毛片网站观看| 日本高清在线看免费观看| 亚洲另类国产欧美一区二区| 777国产精品永久免费观看| 亚洲视频欧美不卡| 国产成人综合在线观看| 国内a级毛片| 色综合激情网| 亚洲天堂在线免费| 国产精品黑色丝袜的老师| 亚洲人成在线精品| 91美女视频在线| 欧美成一级| 国产日韩欧美中文| 国产精品主播| 日韩精品无码不卡无码| 国产亚洲精品无码专| 91人妻在线视频| 亚洲国产天堂久久综合226114| 亚洲无码视频一区二区三区| 一级毛片免费不卡在线| 青青青国产在线播放| 波多野结衣视频网站| 中文一级毛片| 精品无码人妻一区二区| 国产jizzjizz视频| 国产黄色片在线看| 久久99国产综合精品1| 亚洲成aⅴ人片在线影院八| 国产成人综合久久| 一本久道久久综合多人| 亚洲娇小与黑人巨大交| 亚洲av片在线免费观看| 五月婷婷综合色| 毛片网站观看| 好紧好深好大乳无码中文字幕| 亚洲狠狠婷婷综合久久久久| 91九色国产在线| 国产一区二区三区夜色| 99热这里只有免费国产精品| 777国产精品永久免费观看| 在线观看免费黄色网址| 亚洲无码视频图片| 欧美福利在线| 亚洲成人网在线播放| 四虎成人精品| 1级黄色毛片| 亚洲91在线精品| 国产精品无码久久久久AV| 欧美性精品不卡在线观看| 欧美另类图片视频无弹跳第一页| 亚洲三级成人| 午夜影院a级片| 久久伊人操| AV在线麻免费观看网站| 国产1区2区在线观看| 成人福利在线观看| 亚洲欧洲日本在线| 免费毛片a| 天天操天天噜| 亚洲高清无码久久久| 91香蕉国产亚洲一二三区| 国产打屁股免费区网站| 亚洲Va中文字幕久久一区 | 国产69囗曝护士吞精在线视频 | 国产精品毛片在线直播完整版| 69国产精品视频免费| 国产麻豆91网在线看| 亚洲福利视频网址| 亚洲娇小与黑人巨大交|