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

稀杯盔形珊瑚轉錄組分析

2015-09-21 07:38:51劉金豆黃元佳
廣東海洋大學學報 2015年6期
關鍵詞:數據庫功能

劉金豆,劉 麗,黃元佳

(廣東海洋大學水產學院//南海水產經濟動物增養殖廣東省普通高校重點實驗室//南海生物資源開發與利用協同創新中心 廣東 湛江 524088)

稀杯盔形珊瑚轉錄組分析

劉金豆,劉麗,黃元佳

(廣東海洋大學水產學院//南海水產經濟動物增養殖廣東省普通高校重點實驗室//南海生物資源開發與利用協同創新中心 廣東 湛江 524088)

以高通量測序技術Illumina HiSeqTM2000對稀杯盔形珊瑚(Galaxea astreata)進行轉錄組測序分析,共獲得50 360 620條短序列 (reads)。利用Trinity軟件對所有reads從頭組裝后得到81 014條單基因簇(Unigenes),60 471條(74.58%)編碼蛋白框(Coding Sequences,CDs)。與N(rNon-redundant, 非冗余)、COG(Cluster of Orthologous Groups of proteins, 蛋白相鄰類的聚簇)、KEGG(Kyoto encyclopedia of genes and genomes, 京都基因與基因組百科全書)、Swissprot四大數據庫比對共獲得36 545條注釋基因。其中與COG數據庫比對獲得14 491條注釋基因,分為24個功能類別,其中參與一般功能預測類的Unigene數最多,有4 642條;KEGG分析比對獲得16 021條注釋基因,分成241類,包括代謝通路、鈣離子信號通路、MAPK信號通路等,在所有通路中參與代謝途徑的基因數最多,共有2 450條(15.29%);GO(Gene Ontology, 基因本體)功能分類將Unigene分為47個類別。

轉錄組;稀杯盔形珊瑚;Unigene;信號通路

珊瑚礁生態系統是全球最大的生態系統之一,在熱帶海洋環境中扮演著極為重要的角色。近年由于溫度變化、海洋酸化、海水污染等環境因素的影響,導致世界范圍內珊瑚礁大面積白化和死亡[1-3]。李淑等[4]在珊瑚礁白化研究進展中從細胞機制和光抑制機制兩個方面綜述了珊瑚白化的機制,但尚未研究清楚。基因層面上,由于珊瑚研究起步晚,缺乏珊瑚基因信息,使得分子層面上珊瑚的白化機制也十分模糊[5-6]。通過轉錄組測序技術獲得珊瑚基因信息和功能,了解珊瑚的轉錄組特征,將有助在基因水平上揭示珊瑚的白化機制。

轉錄組(transcriptome)是指細胞或組織所表達的所有RNA總和。它反映特定條件下細胞或組織內的基因表達情況和生物學代謝途徑的調節。轉錄組測序技術是最近發展起來的利用深度測序進行轉錄組分析的技術,該技術能提供全面的轉錄組信息,且具有靈敏度高、分辨率高、重復性好的特性。目前,該技術已經大規模應用。2009年Feng等[7]對靜水椎實螺(Lymnaea stagnalis)中樞神經系統進行了轉錄組分析,獲得7 712條EST序列,注釋得到了很多基因。同年,Eli Meyer[8]等人首次將高通量測序技術應用到珊瑚上,獲得了大量的珊瑚基因和注釋信息。2011年Wei等[9]利用該技術進行芝麻轉錄組研究,并對芝麻轉錄組整體特征進行了初步的分析,獲得86 222條Unigene序列,在很大程度上豐富了芝麻生物信息數據。除此之外,還有大量生物進行了轉錄組測序分析,如模式生物中的玉米(Zea mays L.)[10-11]、擬南芥(Arabidopsis thaliana (L.) Heynh)[12-13]等;非模式生物中的油菜(Brassica campestris L.)[14]、白鮭(Salmonidae)[15-16]、大比目魚(Hippoglossus stenolepis)[17]等。

稀杯盔形珊瑚(Galaxea astreata) 隸屬于枇杷珊瑚科,盔形珊瑚屬。本研究以廣東徐聞縣珊瑚礁國家級自然保護區一種造礁石珊瑚——稀杯盔形珊瑚為研究對象,進行轉錄組測序及數據分析。由于受到環境變化的影響,導致保護區內珊瑚礁數量逐年減少,稀杯盔形珊瑚即為受害者之一。自2009年以來,本課題組從分子分類[18-19]、性腺發育[20]以及功能基因克隆[21]等方面對稀杯盔形珊瑚進行了研究。本研究通過高通量測序技術對該珊瑚進行轉錄組測序,以期豐富該珊瑚的轉錄組特征和信息,以及為珊瑚礁的保護和修復提供理論依據。

1 材料和方法

1.1材料

稀杯盔形珊瑚(Galaxea astreata )樣本均采自廣東省湛江市徐聞珊瑚礁國家級自然保護區(20°10′36″N—20°27′00″N)。于室內循環水族箱暫養15 d左右,自然光照11~12 h·d–1, 海水水溫為(27±1.0),℃ 每周更換1/3 體積海水并清洗水循環系統。

1.2方法

1.2.1RNA的提取剪取珊瑚蟲體3~5個,采用Connolly et al[22]的Trizol法提取珊瑚總RNA, -80℃長久保存。部分RNA送廣州基迪奧生物科技有限公司進行檢測、建庫、組裝、測序等。

1.2.2RNA的質量檢測和文庫的構建稀杯盔形珊瑚總RNA經Aglilent 2100檢測,符合轉錄組RNA 檢測標準[質量濃度≥200 ng/μL,總量≥10~15 μg,OD260/280在1.8~2.2范圍,RNA Ratio (28S∶18S) ≥1 , RIN≥7.5]后,進行建庫。

樣品總RNA,用帶有Oligo(dT)的磁珠富集mRNA。加入fragmentation buffer將mRNA打斷成短片段,以mRNA為模板,用六堿基隨機引物(random hexamers)合成第一條cDNA鏈,然后加入緩沖液、dNTPs、RNase H和DNA polymerase I合成第二條cDNA鏈,在經過QiaQuick PCR試劑盒純化并加EB緩沖液洗脫之后做末端修復, 加poly(A)并連接測序接頭,再用瓊脂糖凝膠電泳進行片段大小選擇,進行PCR擴增,測序文庫建好后用Illumina HiSeqTM2000進行測序。

1.3數據處理

1.3.1初始數據處理測序產生的原始圖像數據經base calling轉化為序列數據(raw reads),過濾后得到clean reads。

1.3.2序列組裝使用短 reads 組裝軟件 Trinity做轉錄組從頭組裝。首先將具有一定長度overlap的reads連成更長的片段,從這些通過reads overlap關系得到的不含N的組裝片段得到組裝出來的 Unigene。

1.3.3功能注釋及編碼蛋白框(CDS)的預測通過blastx將Unigene序列比對到Nr (Non-redundant,非冗余)、COG(Cluster of Orthologous Groups of proteins,蛋白相鄰類的聚簇)、KEGG(Kyoto encyclopedia of genes and genomes,京都基因與基因組百科全書)、Swissprot四大蛋白質據庫(evalue<0.000 01)。獲取跟給定 Unigene 具有最高序列相似性的蛋白,從而得到該 Unigene 的蛋白功能注釋信息。對blast比對結果中rank最高的蛋白確定該Unigene的編碼區序列,然后根據標準密碼子表將編碼區序列翻譯成氨基酸序列,從而得到該Unigene編碼區的核酸序列(序列方向5'->3')和氨基酸序列。最后,與以上蛋白庫均比對不上的Unigene用軟件ESTScan預測其編碼區,獲取其編碼區的核酸序列(序列方向5'->3')和氨基酸序列。

2 結果與分析

測序共獲得50 360 620條clean reads,Q20(質量不低于20的堿基比例)達到96.44%,不確定的堿基比例為0,GC%值為49.65%。由此可見此次轉錄組測序結果較好,可以進行后續的數據組裝和分析。組裝拼接后獲得81014條unigene,N50(將所有 Unigene 從長到短排序,并依次累加長度。當累加片段長度達到總片段長度的50%時,對應片段的長度和數量,即為Unigene N50長度和數量)為1 096。這些Unigene中長度最長的為18 230 nt,最短的為201 nt,平均長度為751.55 nt。整體看來,本研究測序結果較好。

2.1Unigene長度分布

從稀杯盔形珊瑚的Unigene長度分布(圖1)中看出,所有Unigene的長度都大于等于200 nt。長度在200~299 nt的Unigene數目最多為19 869條,其次為長度大于或等于3 000 nt的Unigene數為1 419條。而長度在2 900~2999 nt的Unigene數目最少為149條。

圖1 稀杯盔形珊瑚的Unigene長度分布Fig.1 Length distribution of Galaxea astreata-Unigene

2.2Unigene的功能注釋

通過與Nr、COG、KEGG、Swissprot等四大數據庫進行比對共獲得36 545條注釋基因。Nr數據庫比對注釋獲得的Unigene數最多,為35 548條。占總Unigene數的43.88%。KEGG數據庫比對注釋的Unigene數最少,為16 021條,占總Unigene數的19.78%。COG、Swissprot數據庫比對注釋的Unigene數分別為14 491條和27 872條。未注釋的Unigene數為44 469條,占總Unigene的比率為54.89%,超過Unigene總數的一半以上。在Unigene與四大數據庫比對注釋時,有一些Unigene基因會在數據庫中重復出現,多次被注釋。稀杯盔形珊瑚的四大數據庫注釋維恩圖簡明清晰的展示出四大數據庫注釋的Unigene之間的關系(圖2)。由圖2可見,共同注釋的Unigene數為9 089條,占注釋的Unigene總數的24.87%。

圖2 稀杯盔形珊瑚的四大數據庫注釋維恩Fig.2 Venn of the four databases Galaxea astreata-unigene

2.3Unigene的功能分類

2.3.1Unigene的COG分類COG功能分類將Unigene分為24個類別,如圖3所示。從圖中看出Unigene的COG功能種類比較全面,幾乎包含了所有的生命活動。參與一般功能預測類的Unigene數最多,有4 642條,所占比率為32.03%。細胞核結構功能類別含有的Unigene只有9條,所占比率為0.62%,是所有類別中Unigene數最少的。其他的功能分類詳見圖3。

圖3 稀杯盔形珊瑚的Unigene COG功能分類Fig.3 COG Function Classification of Galaxea astreata-Unigene

2.3.2Unigene的KEGG分類與KEGG數據庫進行比對,獲得16 021條注釋Unigene,分為241條途徑,包括代謝通路、MAPK信號通路、鈣離子信號通路、蛋白質在內質網加工過程、細胞周期、RNA轉運等多種途徑。部分途徑列于下表1。從表1中可見,代謝途徑所占有的Unigene數最多,為2 450,占15.29%。其次為MAPK信號通路,共有Unigene數為657條,所占比率為4.1%。鈣離子信號通路所占有的Unigene數為610條,占3.81%,位居第三。

表1 稀杯盔形珊瑚Unigene的KEGG功能分類Table 1 KEGG function Classification of Galaxea astreata-Unigene

2.3.3Unigene的GO分類如圖4所示,GO功能分類將注釋的Unigene分為3類:生物過程、細胞組成和分子功能,共47個分支。其中生物過程包含23個不同的類別,是三大類別中所含類別最多的一類。所含類別如下:細胞生理過程、新陳代謝的過程、生物調節、生長、死亡等。細胞生理過程包含的Unigene數目最多為8 185條,其次為新陳代謝的過程,有7 192條Unigene。死亡類別的Unigene數最少,只有8條。

細胞組成和分子功能又分別分為13和11個類別。細胞組成中所含Unigene數最多的類別為細胞,擁有7 993條Unigene。細胞部分包含的Unigene數比細胞類只少兩條,位列第二。而突觸部分僅僅有35條Unigene,是細胞組成中Unigene數目最少的一類。分子功能中的催化活性所占有的Unigene數最多有7 359條。電子轉運活性所占Unigene數最少,只有10條。

2.4編碼蛋白框(CDs)的預測

通過blastx比對(evalue < 0.000 01)和ESTScan預測Unigene編碼區,得到其編碼框。兩者共得到編碼蛋白框60 471條(74.58%),大部分Unigene的蛋白編碼框已經得到預測。其中blast比對得到36 564條CDs,占總Unigene 數的比為45.13%,ESTScan預測得到的CDs為23 907條,占Unigene總數的29.51%。各自得到的CDs長度分布圖具體如圖5和圖6所示,由圖可知,所得CDs的長度在200~299 nt之間所占比例最大。

圖4 稀杯盔形珊瑚的Unigene GO功能分類Fig.4 GO class of Galaxea astreata-Unigene

3 討論

本研究通過Illumina HiSeqTM2000對稀杯盔形珊瑚(Galaxea astreata)進行轉錄組測序,共獲得50 360 620 reads,組裝后獲得81 014條Unigene。比對后得到36 545條注釋基因,未注釋的基因有44 469條,占總Unigene的百分比為54.89%。出現大量未能注釋的Unigene大致包括以下2個因素。

(1)珊瑚基因信息的缺乏。相對于其他物種如:擬南芥、斑馬魚而言,珊瑚的研究起步較晚,目前關于珊瑚基因組和轉錄組的研究很少,其數據庫信息不足,缺乏可供參考的基因組信息。這致使對基因進行功能注釋時往往找不到對應的注釋信息。

(2)測序技術的局限性。高通量測序技術雖然有很多的優勢,但是同時也存在一定的不足。由于珊瑚沒有參考基因組信息,只能進行重頭轉錄組測序分析。重頭組裝測序中序列長度較短的會影響后期數據質量。根據魏利斌等[23]的研究認為轉錄組序列越短,獲得注釋信息的可能性就越小。

從實驗結果Unigene長度分布圖和N50值的大小來看,Unigene長度都是在200 nt及其以上,N50值為1 096(Unigene N50越長,數量越少,說明組裝質量越好。),說明本實驗獲得的數據較好。同時,將NCBI數據庫中已有的10條稀杯盔形珊瑚的基因序列與本研究測序結果進blast比對發現,10條序列都能夠比對至本實驗的數據庫中,并且序列之間的相似程度很大,可靠性很高。如,將NCBI中稀杯盔形珊瑚的Fe基因序列(登錄號:KJ59 1050.1)與本研究測序的數據庫進行本地blast比對,找到了其對應的基因,并且結果顯示其S值為1 489 Bits,E值為0.0,表明兩序列具有高度同源性,并且可靠性很高,說明本研究測序數據的準確性較高。

圖6 ESTscan預測的CDs長度分布Fig.6 Length distribution of ESTscan.CDs

與Nr、COG、KEGG、Swissprot等四大數據庫進行比對共獲得36 545條注釋基因。GO功能分類中細胞生理過程包含的Unigene數目最多為8 185條,其次為新陳代謝的過程,有7 192條Unigene;死亡類別的Unigene數最少,僅有8條。KEGG數據庫比對共注釋16 021條Unigene,分為241條途徑,其中代謝途徑的Unigene數最多,有2 450條,所占比例為15.29%。鈣離子信號通路的Unigene數為610條,所占比例為3.81%,位居第三。GO功能分類與KEGG Pathway分析結果一致,發現珊瑚代謝相關的Unigene數目幾乎都是最多,表明珊瑚在生長生活過程中代謝活動非常旺盛。從KEGG Pathway分類結果看參與珊瑚鈣離子信號通路的Unigene數目很多,說明珊瑚鈣化在珊瑚生長過程中具有重要作用。Souter P等[24]人的研究表明珊瑚白化中珊瑚Ca2+平衡失調,這說明鈣化基因參與調解珊瑚的白化。本研究轉錄組測序找到了鈣離子信號通路,并發現了相關鈣化基因,這為研究其調控機制,揭示珊瑚白化機制提供可能,也為珊瑚的保護提供理論基礎和研究方向。

廣東省徐聞縣徐聞珊瑚礁國家級自然保護區給予本研究幫助和支持,特致謝忱!

[1]HARVELL C D,MERKEL S,ROSENBERG E,et al. Coral disease,environmental drivers, and the balance between coral and microbial associates[J]. Oceanography,2007,20(1):172-195.

[2]HOEGH-GULDBERG O,MUMBY P J,HOOTEN A J,et al. Coral reefs under rapid climate change and ocean acidification[J]. Science,2007,318(5857):1737-1742.

[3]HUGHES T P,BAIRD A H,BELLWOOD D R,et al. Climate change,human impacts,and the resilience of coral reefs[J]. Science,2003,301(5635):929-933.

[4]李淑,余克服. 珊瑚礁白化研究進展[J]. 生態學報,2007,27(5):2059-2069.

[5]BANUMS I B. A restoration genetics guide for coral reef conservation[J]. Mol Ecol,2008,17(12):2796-2811.

[6]DAY T,NAGEL L,VAN OPPEN M J,et al. Factors affecting the evolution of bleaching resist- ance in corals[J]. Am Nat,2008,171(2):E72-88.

[7]FENG Z P,ZHANG Z,VAN KESTEREN R E,et al. Transcriptome analysis of the central nervous system of the mollusc Lymnaea stagnalis BMC Genomics[J]. BMC Genomics, 2009,451(10):1471 -2164.

[8]MEYER E,AGLYAMOVA GV,WANG S,et al. Sequencing and de novo analysis of a coral larval transcriptome using 454 GSFlx. BMC Genomics,2009,10(1):219.

[9]WEI W L, QI X Q,WANG L H, et al. Characterization of the sesame (Sesamum indicum L.) global transcriptome using Illumina paired-end sequencing and development of EST-SSR markers[J]. BMC Genomics,2011,12:451.

[10]EMRICH S J,BARBAZUK W B,LI L,et al. Gene discovery and annotation using LCM- 454 transcriptome sequencing[J]. Genome Res,2008,17(1):69–73.

[11]OHTUS K,SMITH M B,EMRICH S J,et al. Global gene expression analysis of the shoot apical meristem of maize(Zea mays L.)[J]. Plant J,2007,52(3):391–404.

[12]JONES-RHOADES M W,BOREVITZ J O,PREUSS D. Genome-wide expression profiling of the Arabidopsis female gametophyte identifies families of small,secreted proteins [J]. PloS Genet,2007,3(10):1848–1861.

[13]WEBER A P M,WEBER K L,CARR K,et al. Sampling the arabidopsis transcriptome with massively parallel pyrosequencing [J]. Plant Physiol,2007,144(1):32-42.

[14]TRICK M,LONG Y,MENG J L,et al. Single nucleotide polymorphism(SNP)discovery in the polyploid Brassica napus using Solexa transcriptome sequencing [J]. Plant Biotechnol J,2009,7(4):334–346.

[15]RENAUT S,NOLTE A W,BERNATCHEZ L. Mining transcriptome sequences towards identifying adaptive single nucleotide polymorphisms in lake whitefish species pairs(Coregonus spp. Salmonidae)[J]. Mol Ecol,2010,19(S1):115–131.

[16]JEUKENS J,RENAUT S,ST-CYR J,et al. The transcriptomics of sympatric dwarf and normal lake whitefish(Coregonus clupeaformis spp,Salmonidae)divergence as revealed by next-generation sequencing[J]. Mol Ecol,2010,19(24):5- 389–5403.

[17]PEREIRO P,BALSEIRO P,ROMERO A,et al. Highthroughput sequence analysis of turbot(Scophthal musmaximus)transcriptome using 454-pyrosequencing for the discovery of antiviral immune genes[J]. PLoS ONE,2012,7(5):e35369.

[18]劉麗,李曉娜,陳育盛,等. 基于線粒體基因的石珊瑚分子系統學研究[J]. 海洋與湖沼,2012,43(4):814-820.

[19]劉麗,陳育盛,李曉娜,等. 基于線粒體Cyt b基因的10種石珊瑚的系統發育關系[J]. 廣東海洋大學學報,2011,31 (1):6-11.

[20]金磊. 盾形陀螺珊瑚和稀杯盔形珊瑚性腺發育與生長規律的研究[D]. 湛江:廣東海洋大學,2014.

[21]范程輝,劉麗,沈城,等. 稀杯盔形珊瑚銅鋅超氧化物歧化酶基因全長cDNA序列的克隆與分析[J]. 熱帶海洋學報,2015,34(1):83-89.

[22]CONNOLLY M A,CLAUSEN P A,LAZAR J G. Purification of RNA from animal cells using Trizol[J]. Cold Spring Harb Protoc,2006(1): 10-15. DOI:10. 1101/pdb.prot4104.

[23]魏利斌,苗紅梅,張海洋. 芝麻發育轉錄組分析[J]. 中國農業科學,2012,45(7):1246-1 256.

[24]SOUTER P,BAY L K,ANDREAKIS N,et a1. A multilocus, temperature stress-related gene expression profile assay in Acropora millepora,a dominant reef-building coral[J]. Mol Ecol Resour,2011,11(2):328-334.

(責任編輯:陳莊)

Transcriptome Analysis of Galaxea astreata

LIU Jin-Dou,LIU Li,HUANG Yuan-Jia
(Fisheries College,Guangdong Ocean University//Key Laboratory of Aquaculture in South China Sea for Aquatic Economic Animal of Guangdong Higher Education Institutes//South China Sea Bio-Resource Exploitation And Utilization Collaborative Innovation Center,Zhanjiang 524088,China)

High-throughput sequencing technology Illumina HiSeqTM2000 was used to sequence the Galaxea astreata transcriptome,and 50 360 620 reads were got in all.And then assembled software Trinity was used to denovo assembled the reads,getting 81014 Unigenes and 60 471 Coding Sequences(CDs). Blasted with the four database Nr(Non-redundant), COG(Cluster of Orthologous Groups of proteins),KEGG(Kyoto encyclopedia of genes and genomes), Swissprot,36 545 annotated Unigenes were gotten. Blasted with COG database,14 491 unig- enes were gained and according to the function it was divided into 24 categories. 4 642 Unigenes join in the general function prediction only category, which own the most numbers of Unigenes. Blasted with KEGG pathway,16 021 Unigenes were annotated, and it was divided into 241 categories according to different pathways,including metabolic pathway, calcium signaling pathway, and MAPK signaling pathway and so on. Of all the pathways, the number of genes join in the metabolic pathways was maximum in 2450 (15.29%). GO(Gene Ontology)functional classification divided all Unigenes into 47 categories. Through this study,a lot of Galaxea astreata coral Unigenes were gotten. It can help us find some new functional genes and to know about the growth of coraland adaptation to the environment changes.

transcriptome;Galaxea astreata;Unigene;signal pathway

S917.4

A

1673-9159(2015)06-0001-08

10.3969/j.issn.1673-9159.2015.06.001

2015-09-21

國家海洋公益性行業科研專項(201105012);廣東省自然科學基金項目(S2011010000269); 廣東省海洋漁業科技推廣專項(A201308E02)

劉金豆(1989—),男,碩士研究生,研究方向為海洋經濟動物發育生物學的研究。E-mail:735677865@qq.com

劉麗,女,教授,主要從事海洋經濟動物發育生物學的研究。E-mail:zjouliuli@163.com.

猜你喜歡
數據庫功能
也談詩的“功能”
中華詩詞(2022年6期)2022-12-31 06:41:24
關于非首都功能疏解的幾點思考
數據庫
財經(2017年15期)2017-07-03 22:40:49
數據庫
財經(2017年2期)2017-03-10 14:35:35
懷孕了,凝血功能怎么變?
媽媽寶寶(2017年2期)2017-02-21 01:21:24
“簡直”和“幾乎”的表達功能
數據庫
財經(2016年15期)2016-06-03 07:38:02
數據庫
財經(2016年3期)2016-03-07 07:44:46
數據庫
財經(2016年6期)2016-02-24 07:41:51
中西醫結合治療甲狀腺功能亢進癥31例
主站蜘蛛池模板: 婷婷综合在线观看丁香| 精品一区二区三区无码视频无码| www.狠狠| 中文字幕欧美成人免费| 中文字幕亚洲另类天堂| 午夜欧美在线| 米奇精品一区二区三区| 亚洲系列无码专区偷窥无码| 成人午夜久久| 国产成人综合亚洲网址| 成年人福利视频| 亚洲男人的天堂久久精品| 亚洲专区一区二区在线观看| 欧美伦理一区| 无码免费视频| 成人一区专区在线观看| 国产精品分类视频分类一区| 99尹人香蕉国产免费天天拍| 九九免费观看全部免费视频| 婷婷激情亚洲| 久久久久中文字幕精品视频| 中文字幕久久亚洲一区 | 国产成人夜色91| 91精品在线视频观看| 国产成人精品一区二区不卡| 中国精品自拍| 日韩大乳视频中文字幕 | 久无码久无码av无码| 亚洲精品欧美重口| 色135综合网| 国产精品女主播| 波多野结衣爽到高潮漏水大喷| 国产精品女在线观看| 日本黄色不卡视频| 色婷婷在线播放| 精品精品国产高清A毛片| 国产福利在线免费观看| 亚洲国产成人超福利久久精品| 97国产在线播放| 亚洲欧美日韩天堂| 久久伊人操| 9啪在线视频| 狠狠躁天天躁夜夜躁婷婷| 亚洲精品男人天堂| 日韩在线视频网站| 国产一级片网址| 欧洲欧美人成免费全部视频| 国产男女免费视频| 久久久久亚洲av成人网人人软件| 91探花在线观看国产最新| 婷婷午夜影院| 国产综合色在线视频播放线视| 69综合网| 国产精品99一区不卡| 国产18在线| 国产av色站网站| a级毛片视频免费观看| 成人福利在线视频| 亚洲日韩Av中文字幕无码| 国产欧美日韩va| 国产美女在线观看| 大香网伊人久久综合网2020| 人妻一区二区三区无码精品一区| 久久99热这里只有精品免费看| 国产日本欧美亚洲精品视| 免费毛片全部不收费的| 又粗又硬又大又爽免费视频播放| 国产乱肥老妇精品视频| 国产免费好大好硬视频| 国产精品成人AⅤ在线一二三四| 四虎影视库国产精品一区| 国产丝袜丝视频在线观看| 国产微拍一区二区三区四区| 中文精品久久久久国产网址| 91麻豆精品国产高清在线| 国产精品30p| 婷婷综合在线观看丁香| 国产91无码福利在线| 一级一级一片免费| 在线国产欧美| 18禁高潮出水呻吟娇喘蜜芽| 一级毛片在线播放免费观看|