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

2個銀葉真蘚HSP70序列的生物信息學相關分析

2016-01-15 01:37:42王明強,高貝,張道遠
生物信息學 2015年1期

2個銀葉真蘚HSP70序列的生物信息學相關分析

王明強1,2,高貝1,張道遠1*

(1.中國科學院干旱區生物地理與生物資源重點實驗室,中國科學院新疆生態與地理研究所,烏魯木齊 830011;

2.中國科學院大學,北京 100049)

摘要:從銀葉真蘚(Bryum argenteum)轉錄組數據庫出發,使用Pfam數據庫提供的HMM模型共得到33條長度大于200aa,注釋的熱休克蛋白BaHSP70;其中2條(BaHSP70-1, BaHSP70-2)具有完整ORF,NCBI核酸數據庫登錄號為KP087877和KP087878。使用生物信息學在線分析工具和軟件,對真蘚HSP70的兩條蛋白質序列從氨基酸組成、保守結構域、理化性質、疏水性/親水性、信號肽、蛋白質結構、模體的識別及同源性分析等方面進行了預測和分析。結果表明:2條BaHSP70s基因序列ORF全長分別為2 396 bp和2 356 bp,分別編碼649aa和650aa。序列模體分析表明BaHSP70s和其它報道的植物HSP70均含有4個相同的模體,并且各模體在蛋白質序列上順序一致。通過對2條BaHSP70s進行氨基酸多序列比對及基因樹分析,發現BaHSP70-1和BaHSP70-2雪蓮相似度最高,分別是91.2%和86.6%。本研究為進一步研究HSP70基因的克隆和功能驗證奠定了基礎。

關鍵詞:真蘚;熱激蛋白70;RNA-Seq;生物信息學

中圖分類號:Q811.4文獻標志碼:A

收稿日期:2014-12-18;修回日期:2015-01-06.

基金項目:國家自然科學基金(No.31371685) 。

作者簡介:劉巧紅,女,副研究員,研究方向:甜菜遺傳育種與植物分子標記;E-mail:Qiaohong941@hit.edu.cn.

doi:10.3969/j.issn.1672-5565.2015.01.02

Bioinformatics relevant analysis of two heat shock protein 70 sequence

fromBryumargenteum

WANG Mingqiang1,2, GAO Bei1, ZHANG Daoyuan1*

(1.KeyLaboratoryofBiogeographyandBioresourceinAridLand,XinjiangInstituteofEcologyand

Geography,ChineseAcademyofSciences,Urumqi830011,China;

2.UniversityofChineseAcademyofSciences,Beijing10049,China)

Abstract:Our study utilizes the RNA-seq technology which is based on the Hi-seq 2000 sequencing platform from BGI to acquire the RNA sequences data from Bryum argenteum dealt with dehydration and rehydration stress. Using the Hidden Markov Model provided by the Pfam database, we have identified 33 HSP70 sequences with the length longer than 200aa. Two of them (BaHSP70-1,BaHSP70-2) contain Open Reading Frame, the Accession No. in NCBI nucleotide database are KP087877 and KP087878, respectively. Applying web-based bioinformatics tools and softwares, we predicted and analyzed B. argenteum Heat Shock Protein 70 sequences from the following aspects: amino acid component, conserved domain, physicochemical property, hydrophobicity and hydrophily, singal peptide, protein structure, motif recognition, homologous analysis. The result showed the length of two possessed complete ORF HSP70s are 649aa and 650aa respectively. The length of corresponding anscripts in the transcriptome database are 2 396 bp and 2 356 bp. Sequence motif analysis indicates B. argenteum and other plant species have 4 identical motifs. Meanwhile, the order of the motifs arranged in protein sequence is the same in all selected species. Multiple sequence alignment and homologous analysis demonstrated that HSP70-1 and HSP70-2 are most similar to Saussurea up to the similarity of 91.2% and 86.6% respectively. This study laid fundation for the B.argenteum further research in gene clone and functional confirmation.

Keywords:Bryum argenteum; Heat shock protein 70; RNA-Seq; Bioinformatics

熱激蛋白70( Heat Shock Protein 70, HSP70) 是生物機體應對高溫及其他脅迫環境時所產生的一類特定的應激蛋白,在細胞內的大量表達可以明顯改善細胞的生存能力,提高生物機體對環境脅迫或傷害的耐受性[1]。HSP70在原核生物和真核生物體內都存在,并存在于所有的細胞區室和器官中,序列非常保守;在細胞內參與多肽的折疊[2]、寡聚蛋白的組裝和跨膜運輸、蛋白質活性的調節,影響植物的生長發育,是十分重要且基本的細胞功能蛋白[3-4]。

銀葉真蘚(Bryumargenteum)隸屬于真蘚科真蘚屬,是蘚類植物的典型代表種[5]。它的形態結構相對簡單,但是可以生存在極端干旱、高溫、高鹽的環境中,是中國荒漠區生物土壤結皮層的優勢成分和群落演替過程中的主要先鋒種[6],也是生態系統穩定和退化生態系統恢復的重要指示植物之一[7]。近年來,關于銀葉真蘚的形態發育[8],系統分類[9],生理生化的研究已有報道[10],但有關其HSP70蛋白生物信息學分析未見報道。

前期研究過程中,項目組構建了銀葉真蘚干燥-復蘇過程的轉錄組數據庫。本研究從轉錄組數據庫出發,使用Pfam數據庫中提供HSP70的HMM模型,從已注釋的真蘚BaHSP70s中得到了2條全長的BaHSP70 序列,分別命名為BaHSP70-1(GenBank No. KP087877)和BaHSP70-2(GenBank No.KP087878)。利用生物信息學相關軟件分析了這2條基因的序列特征、編碼蛋白質的理化性質、蛋白質結構、相關的生物學功能預測等,并與其他報道的HSP70進行了多序列比對與基因樹分析[11]。研究有助于進一步理解BaHSP70s蛋白抗逆調控機制,并為后期BaHSP70s基因實驗克隆和功能驗證奠定基礎。

1材料與方法

1.1數據收集與預處理

本研究實驗材料銀葉真蘚取自于新疆古爾班通古特沙漠,通過無菌操作,在人工控制條件下進行組織培養獲得再生完整植株,提取銀葉真蘚在干旱-復水過程中的RNA構建RNA測序文庫[12];利用華大基因提供的Illumina二代測序RNA-seq技術[13],使用Trinity軟件對序列進行組裝[14],同時利用相關蛋白質數據庫對組裝的Contig進行注釋,從而得到銀葉真蘚轉錄組Unigene數據庫和蛋白質數據庫。綠藻植物、小立碗蘚、擬南芥、水稻、玉米等植物的HSP70基因從NCBI的PROTEIN數據庫和UniProtKB數據庫中下載得到。從Pfam數據庫中下載HSP70的HMM概率模型[15],文件命名為hsp70.hmm。

1.2真蘚HSP70基因的鑒定

使用HMMER-3.1軟件搜索HSP70的同源序列[16],其中e-value值設置為5,輸出格式為pfamtblout,輸入命令:hmmsearch-E5-pfamtblout HSP1 hsp70.hmm Protein > result,獲得具有HSP70結構域的序列;根據輸出的比對信息,提取比對結果具有全長的片段,再次輸入Pfam提供的在線工具進行驗證,如圖1所示,具有完整的結構域。同時使用NCBI提供的在線分析工具ORF Finder確認序列的起始密碼子、終止密碼子和長度;Blast兩兩比對結果證明序列間的一致性達83%,因此分別命名為BaHSP70-1和BaHSP70-2。

圖1 基于Pfam數據庫確認HSP70蛋白質結構域

1.3HSP70的生物信息學分析

使用ProtParam在線分析工具檢索蛋白質序列的理化參數,包括分子質量,氨基酸組成,消光系數等;使用CDD分析保守結構域;使用ProtScale程序繪制親疏水性序列譜,反映蛋白質在水環境下的折疊情況;使用NetNES分析信號肽位點;使用NPSA在線分析工具提供的MRLC方法預測真蘚HSP70蛋白的二級結構;使用瑞士模型工作區(Swiss-Model Workspace)網絡綜合服務器對HSP70蛋白質結構做同源建模;使用ClustalW對獲取的HSP70蛋白序列進行多序列比對;除去比對后序列兩側的冗余序列,提取序列間相似性較高的序列區域,使用MEGA5.2構建NJ進化樹并進行Bootstrap檢測[17],使用工具和地址如表1所示。

2結果與分析

2.1銀葉真蘚HSP70基因序列分析

基于Pfam數據庫提供的HMM模型方法,從銀葉真蘚轉錄組注釋蛋白質庫中鑒定出2條具有全長的HSP70序列,分別命名為BaHSP70-1和BaHSP70-2,它們在Unigene數據庫中對應的的轉錄本核酸序列長度分別為2 396 bp和2 356 bp。用NCBI提供的在線工具ORF Finder預測其編碼蛋白質序列的長度為649和650個氨基酸,開放式閱讀框具有起始密碼子和終止密碼子,多序列比對分析表明,BaHSP70最保守的序列區段位于序列的第200位至250位,氨基酸序列為LDVTPLSLGLETLGGVMTKLIPRNTTIPTKKSQVFSTYAD

NQTGVLIQVY,該區段位于核酸結合區(Nucleotide Binding Domain, NBD),具有結合并水解 ATP 的活性中心,它與ATP的結合-水解過程調控著它對多肽的親和力。

2.2銀葉真蘚HSP70保守域分析

本研究使用NCBI提供的CDD數據庫對BaHSP70s序列進行保守域分析[18],結果表明,銀葉真蘚HSP70序列含有HSP70結構域,其中核酸結合區結構域較為的保守,屬于NBD sugar-kinase HSP70 action superfamily,結果如圖2所示;研究進一步表明基于HMM模型鑒定的完整基因屬于熱激蛋白HSP70家族。

表1 研究使用生物信息學在線分析工具

圖2 BaHSP70保守結構域分析

2.3銀葉真蘚HSP70蛋白理化性質分析

使用ProtParam對HSP70蛋白質序列進行理化性質分析,結果如表2所示:兩個蛋白質序列消光系數相同,BaHSP70-1不穩定指數高于BaHSP70-2,說明一般情況下BaHSP70-2相對較穩定;但是總體平均親水性BaHSP70-1相對較高。

2.4銀葉真蘚HSP70蛋白質親疏水性分析

蛋白質氨基酸組成是蛋白質折疊,發揮功能的主要驅動力;蛋白質折疊時會形成疏水內核和親水的表面,同時在潛在的跨膜區出現高疏水值區域。

使用ProtScale在線工具對BaHSP70蛋白質氨基酸序列的親水性和疏水性進行預測,選取信號最為顯著的Kyte& Doolittle氨基酸標度值,其它參數默認。

表2 BaHSP70一級結構理化性質分析

如圖3所示,橫坐標表示序列位置,縱坐標表示氨基酸的標度值。HSP70-1有11個低分值峰,HSP70-2有9個低分值峰,這些氨基酸的區域如圖方框標記,屬于高親水性區域。其表面往往富集親水性氨基酸,同時也是蛋白質進化中氨基酸插入的主要位點。HSP70-1和HSP70-2疏水性的波形圖不一樣也反映兩種蛋白在理化性質上有所差別。

圖3 真蘚HSP70蛋白質疏水性和親水性預測

2.5真蘚HSP70蛋白信號肽預測和分析

采用基于神經網絡模型(Neural Networks models)和隱馬可夫模型(Hidden Markov Models)的NetNES 1.1 server預測銀葉真蘚熱激蛋白70的富含亮氨酸的核輸出信號[19];結果如圖4所示,BaHSP70-1在序列第400位有富含亮氨酸的核輸出信號,BaHSP70-2在序列第310位和第400位出現富含亮氨酸的核輸出信號。此區段說明BaHSP70有核輸出信號,可以被某些特定的輸出蛋白所結合,引導HSP70蛋白經由核孔進入細胞質特定位置執行相關功能。

圖4 BaHSP70核輸出信號預測

2.6真蘚HSP70蛋白二級結構分析

使用NPSA在線分析工具提供的MRLC方法對BaHSP70-1和BaHSP70-2蛋白質的二級結構預測分析結果表明,HSP70蛋白二級結構主要由α螺旋,β折疊,和無規卷曲三種形式構成,其中α螺旋所占比例較高,N-段主要分布無規卷曲結構,C-端主要分布α螺旋,統計的比例如表3所示,二級結構分布如圖5所示。

表3 銀葉真蘚HSP70的二級結構預測分析

圖5 銀葉真蘚HSP70二級結構

2.7三維結構模型分析

使用瑞士模型工作區(Swiss-Model Workspace)網絡綜合服務器對HSP70蛋白質結構做同源建模[20],建模方式采用自動方式,同源模型建立主要分四個步驟:(1)鑒定模板結構,即首先將蛋白質序列通過Blast比對,發現數據庫中一致性較高的序列;根據服務器分析結果,HSP70-1和HSP70-2同時與數據庫中的Hsc71 kDa protein結構最相似,前者一致性達到55%,后者達54%;(2)目標序列和模板序列結構的比對,即已找到的模板序列為參考序列進行比對,對于氨基酸序列完全一致的區域,建立相同的三級結構并進行結構疊和;(3)建立模型,根據網絡綜合服務器提供的算法,構建蛋白質模型,對于比對出未知的蛋白質,首先找出主鏈結構,然后側鏈建模,最后利用能量計算的方法進行結構優化;(4)模型質量的評估,即對構建好的模型用全局發和局部法進行質量評估。

本次模型評估使用全局法中的QMEANscore4記分值估算。構建的蛋白質三維結構如圖6所示,對于HSP70-1,可信度值為0.52,可信度較高,模板序列方法使用X-ray,分辨率為1.90 ?,比對的范圍在第7~559個氨基酸,覆蓋度達到0.85;對于HSP70-2,可信度值為0.34,模型可參考,模板序列方法使用X-ray,分辨率為1.90 ?,比對的范圍在第7~559個氨基酸,覆蓋度達到0.85。

圖6 HSP70-1和HSP70-2蛋白質的三維結構模型

2.8功能預測

使用CBS的Protfun軟件預測真蘚HSP70基因編碼的功能[21],如表4所示,熱激蛋白70具有轉錄調控、生長因子、結構蛋白的可能性較高,在BaHSP70-1中可能性依次為0.442,0.456,0.436;在BaHSP70-2分別中可能性依次為0.396,0.372,0.465;這符合生物學上植物受到外界環境脅迫時,脅迫誘導產生HSP70蛋白,HSP70作為分子伴侶識別細胞內的信號來引導轉運新合成多肽的折疊以使其形成正確的構象發揮生物學功能的功能。

表4 真蘚HSP70 GO功能預測

2.9序列模體識別

使用MEME程序包對綠藻植物(Helicosporidiumsp)、銀葉真蘚、小立碗蘚(Physcomitrellapatens)、水稻(Oryzasativa)、擬南芥(Arabidopsisthaliana)、玉米(Zeamays)HSP70蛋白進行序列模體識別分析,模體分布參數設置為"Any number of repetitions",其它參數默認。結果如圖7所示,從序列中發現了5個motif,其中以上物種含有公共Motif5, Motif1, Motif2, Motif4,并且模體在HSP70基因上的分布順序是一致的,間隔區長度也相近;其中,銀葉真蘚的HSP70-1和HSP70-2含有相同的模體結構;真蘚HSP70蛋白的保守基序串聯形式與水稻所具有的保守基序最相近,包括所有發現的5種模體,并且模體間的順序完全一致。與綠藻和小立碗蘚比較發現,綠藻HSP70序列的C端和小立碗蘚序列的N端多含有Motif4,但是末端出現的Motif4保守度相對中間區段的Motif較低;與擬南芥和玉米比較發現,銀葉真蘚HSP70序列在C端多含有一個Motif3。從植物由低等進化到高等的進化的角度來看,說明了HSP70序列家族在關鍵的功能結構域上保持著較高的一致性,可能這些保守的、順序一致的模體在HSP70功能中發揮關鍵作用,但同時根據兩端的較低模體的保守度,又豐富了HSP70在進化上的復雜性和多樣性。

圖7 HSP70序列模體分析

2.10同源性分析及基因樹構建

將銀葉真蘚HSP70蛋白質序列與膠球藻、團藻、扁藻等低等綠藻類植物和小立碗蘚,擬南芥、水稻、玉米、小麥等高等植物的HSP70蛋白進行多序列比對發現,從低等綠藻類的植物到高等維管植物的蛋白質序列有較高的同源性。使用ClustalW軟件進行多序列比對,選取比對結果中保守區段的蛋白質序列,使用MEGA5.2軟件Phylogeny提供的鄰位連接法(Neighbor-Joining)對同源的序列構建基因樹,對于參數選項設置:檢驗方法選擇“Bootstrap method”,重復次數為1 000次進行檢測,以確保基因樹分枝的穩定性,其它參數選擇默認值。

構建的基因樹如圖8所示:基因樹大致可以分為5支,高等植物節節草、小麥、菠菜、玉米等聚在一支(分支I);低等植物螺旋藻、膠球藻、團藻、原藻、植物病原藻、扁藻、愛爾蘭海藻、小立碗蘚等聚在一支(分支II);高等植物苜蓿、桑、玉米、節節草、小麥聚在一支(分支III);低等植物愛爾蘭海藻、團藻、植物病原藻、等較低等植物聚在一支(分支IV);高等植物蘋果、菠菜、桑、雪蓮、鐵皮石斛、銀葉真蘚聚在一支(分支V)。選取的5條擬南芥HSP70序列分別位于基因樹不同的分支上,在高等植物和低等植物類群分支中都有分布,說明不同HSP家族內部成員序列差別較大,分布廣泛;而根據下載到HSP70蛋白質序列的注釋信息,分析發現定位于不同亞細胞器的HSP70蛋白距離較近,其中,小立碗蘚、擬南芥、黃瓜的葉綠體HSP70在基因樹上位置較近(菱形標識),而小麥、節節草、原藻、擬南芥的線粒體HSP70在基因樹上位置較近(三角形標識)。此外,銀葉真蘚的兩條HSP70序列與雪蓮的親緣關系較近,聚為一支(圓形標識),同源性分別是91.2%和86.6%,其次是輪藻,同源性分別是76.5%和72.4%。

3結論

本研究以真蘚轉錄組數據庫為基礎,通過Pfam數據庫提供的HSP70隱馬可夫模型從數據庫中鑒定出2條含有完整ORF的HSP70序列:HSP70-1和HSP70-2,二者一致性達83%,全長含有649aa和650aa,其編碼的基因在Unigene數據庫中的長度分別為2 396 bp和2 356 bp。序列模體識別和解析研究發現真蘚和其它植物含有4個相同的模體,并且各模體在蛋白質序列上分布順序一致,可能在生命活動中發揮關鍵作用,N端和C端也發現有保守度比較低的模體,豐富了HSP70在進化上的復雜度和多樣性。真蘚HSP70蛋白屬于親水性蛋白,HSP70-1有11個親水區,HSP70-2有9個親水區。二級結構、三級結構預測發現包括α螺旋,β折疊,β轉角和卷曲螺旋結構。三級結構可信度高,在以后的研究中模型可參考。銀葉真蘚的兩條熱激蛋白序列BaHSP70-1和BaHSP70-2位于高等植物的進化支上,與雪蓮同源性最高,同源性分別是91.2%和86.6%。綜上所述,該研究獲得的真蘚兩條HSP70蛋白的序列和基因序列,并且認識了HSP70蛋白質的生物化學特性,蛋白質結構功能,以及在植物同源關系中的位置。

圖8 不同植物HSP70蛋白質序列基因樹

參考文獻(References)

[1]BALLINGER D G, PARDUE M L. The control of protein-synthesis during heat-shock in drosophila cells involves altered polypeptide elongation rates[J]. Cell,1983,33(1):103-113.

[2]HORWICH A L. Molecular chaperones in cellular protein folding: The birth of a field[J]. Cell,2014,157(2):285-288.

[3]BUKAU B, HORWICH A L. The Hsp70 and Hsp60 chaperone machines[J]. Cell,1998,92(3):351-366.

[4]MAYER M P, BUKAU B. Hsp70 chaperones: Cellular functions and molecular mechanism[J]. Cellular and Molecular Life Sciences,2005,62(6):670-684.

[5]HUI R, LI X R, CHEN C Y,et al. Responses of photosynthetic properties and chloroplast ultrastructure of Bryum argenteum from a desert biological soil crust to elevated ultraviolet-B radiation[J]. Physiol Plantarum,2013,147(4):489-501.

[6]鄭云普, 趙建成, 張丙昌, 等. 荒漠生物結皮中藻類和苔蘚植物研究進展[J]. 植物學報,2009,3(3):371-380.

ZHENG Yunpu,ZHAO Jiancheng, ZHANG Bingchang,et al. Advances on ecological studies of algae and mosses in biological soil crust[J].Chinese Bulletin of Botany, 2009, 3(3):371-380.

[7]BELNAP J. Nitrogen fixation in biological soil crusts from southeast Utah, USA[J]. Biol Fert Soils, 2002, 35(2):128-135.

[8]李利博, 趙建成. 真蘚屬(Bryum Hedw.)蒴齒形態特征及其分類學意義[J]. 植物研究, 2009, 29(6):651-658.

LI Libo, ZHAO Jiancheng. Peristome morphology of bryum hedw, (Musci Bryaceae) and its taxonomic significance[J]. Bulletin of Botanical Research, 2009, 29(6):651-658.

[9]李利博, 趙建成, 曹娜. 河北省真蘚屬(Bryum Hedw.)植物分類學研究[J]. 河北師范大學學報(自然科學版), 2009, 29(6):800-809.

LI Libo, ZHAO Jiancheng, CAO Na. A taxonomical study on genusbryumhedwof Hebei, China[J]. Journal of Hebei Normal University(Natural Science Edition), 2009, 29(6):800-809.

[10]回嶸, 李新榮, 賈榮亮. 增強UV-B輻射對真蘚結皮生理特性的影響[J]. 生態學雜志,2012,31(1):38-43.

HUI Rong, LI Xinrong, JIA Rongliang. Effects of enhanced UV-B radiation on physiological characteristics of Bryum argenteum[J]. Chinese Journal of Ecology, 2012, 31(1):38-43.

[11]SAKUDO A, XUE G A, KAWASHITA N, et al. Structure of the Prion Protein and Its Gene: An Analysis Using Bioinformatics and Computer Simulation[J]. Curr Protein Pept Sc, 2010, 11(2):166-79.

[12]張道遠, 馬瑞, 楊紅蘭, 等. 荒漠生物結皮中齒肋赤蘚總RNA最佳提取方法的確立[J]. 生物技術, 2010, 20(3):52-54.

ZHANG Daoyuan, MA Rui,YANG Honglan,et al. A method established for total RNA extraction from syntrichia caninervis in a temperate desert[J]. Biotechnology, 2010, 20(3):52-54.

[13]WANG Z, GERSTEIN M, SNYDER M. RNA-Seq: a revolutionary tool for transcriptomics[J]. Nat Rev Genet, 2009,10(1):57-63.

[14]HAAS B J, PAPANICOLAOU A, YASSOUR M,et al. De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis[J]. Nature Protocols, 2013, 8(8):1494-1512.

[15]FINN R D, BATEMAN A, CLEMENTS J, et al. Pfam: the protein families database[J]. Nucleic Acids Research, 2014, 42(1):222-230.

[16]FINN R D, CLEMENTS J, EDDY S R. HMMER web server: interactive sequence similarity searching[J]. Nucleic Acids Research, 2011, 39:29-37.

[17]TAMURA K, PETERSON D, PETERSON N, et al. MEGA5: Molecular evolutionary genetics analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods[J]. Molecular Biology and Evolution, 2011,28(10):2731-2739.

[18]MARCHLER-BAUER A, LU S N, ANDERSON J B, et al. CDD: a conserved domain database for the functional annotation of proteins[J]. Nucleic Acids Research, 2011,39:225-229.

[19]lA Cour T, KIEMER L, MOLGAARD A,et al. Analysis and prediction of leucine-rich nuclear export signals[J]. Protein Engineering, Design & Selection : PEDS, 2004,17(6):527-536.

[20]ARNOLD K, BORDOLI L, KOPP J, et al. The Swiss-Model workspace: a web-based environment for protein structure homology modelling[J]. Bioinformatics, 2006, 22(2):195-201.

[21]JENSEN LJ, GUPTA R, STAERFELDT H H,et al. Prediction of human protein function according to Gene Ontology categories[J]. Bioinformatics, 2003, 19(5):635-642.

主站蜘蛛池模板: 亚洲精品动漫| 亚洲精品无码在线播放网站| 久久伊人久久亚洲综合| 久久精品人人做人人| 色老头综合网| 久久香蕉国产线看观看精品蕉| 亚洲精品在线观看91| 99热这里只有精品免费国产| 婷婷亚洲最大| 色婷婷狠狠干| 免费观看三级毛片| 在线观看无码av免费不卡网站| 国产视频一区二区在线观看| 国产精品人成在线播放| 国产一区二区三区在线精品专区| a级毛片一区二区免费视频| 91麻豆精品国产91久久久久| 精品1区2区3区| 亚洲中文字幕在线一区播放| 亚洲精品桃花岛av在线| 国产人碰人摸人爱免费视频| 亚洲欧洲日产国码无码av喷潮| 国产成+人+综合+亚洲欧美| 最近最新中文字幕在线第一页| 成人精品午夜福利在线播放| 国产激情在线视频| 高清无码一本到东京热| 亚洲欧美天堂网| 成人国产精品视频频| 国产女人18水真多毛片18精品 | 亚洲无码精品在线播放| jizz在线免费播放| 久久综合AV免费观看| 在线免费亚洲无码视频| 美女被操91视频| 国产精品一区二区无码免费看片| 欧美日韩成人在线观看 | 成人福利在线看| 精品欧美一区二区三区久久久| 国产精品香蕉| 国产成人麻豆精品| 好紧太爽了视频免费无码| 亚洲天堂精品视频| 欧美不卡在线视频| 第一页亚洲| 日日碰狠狠添天天爽| 真人高潮娇喘嗯啊在线观看| 91视频首页| 夜夜高潮夜夜爽国产伦精品| 成人小视频网| 真实国产精品vr专区| AV天堂资源福利在线观看| 精品一區二區久久久久久久網站| 亚洲欧美日韩另类在线一| 色九九视频| 国产乱子伦无码精品小说 | 91精品专区| 456亚洲人成高清在线| 国产麻豆永久视频| 亚洲男人的天堂在线观看| 97综合久久| 亚洲中文字幕在线观看| 为你提供最新久久精品久久综合| 精品国产免费观看| 一边摸一边做爽的视频17国产| 亚洲Av激情网五月天| 婷婷综合在线观看丁香| 亚洲视频色图| 亚洲国产欧美自拍| 国产精品福利一区二区久久| 日韩欧美国产综合| 女人18毛片久久| 国产毛片片精品天天看视频| 久久久久青草大香线综合精品 | 毛片免费高清免费| 77777亚洲午夜久久多人| 欧美成人午夜在线全部免费| 欧美三級片黃色三級片黃色1| 国产成人禁片在线观看| 免费观看精品视频999| 成年人久久黄色网站| 九九视频免费看|