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

秦艽WRKY轉錄因子家族生物信息學分析

2022-12-08 07:19:32何懿菡尹洋洋孫曉春黃文靜岳正剛
中草藥 2022年23期
關鍵詞:分析

何懿菡,尹洋洋,胡 偉,李 鉑,孫曉春,王 楠,黃文靜,岳正剛

? 藥材與資源?

秦艽WRKY轉錄因子家族生物信息學分析

何懿菡,尹洋洋,胡 偉,李 鉑,孫曉春,王 楠,黃文靜,岳正剛*

陜西中醫藥大學陜西省中藥資源產業化省部共建協同創新中心秦藥資源研究開發國家重點實驗室(培育),陜西 咸陽 712046

利用生信技術分析秦艽WRKY家族,為進一步研究GmWRKY轉錄因子對秦艽次生代謝產物生物合成的分子調控機制提供信息。基于秦艽轉錄組注釋結果挖掘WRKY家族成員,通過NCBI-tBlastx和SMART保守結構域預測進一步篩選得到41條WRKY序列并進行生信分析。基于getorf預測開放閱讀框(open reading frame,ORF)區,對41條序列的ORF片段進行MEME保守基序分析和ClustalW比對,采用MEGAX構建鄰接法(neighbor-joining,NJ)系統發育樹。利用Omicshare GO對41個WRKY成員進行富集分析,利用String進行WRKY成員蛋白互作網絡預測分析。通過轉錄組FPKM值對成員的基因表達情況進行熱圖分析,并對蛋白互作核心成員進行基因相對表達分析。41個秦艽WRKY成員中大部分成員的保守基序為“WRKYGQK”,僅GmWRKY27保守基序發生氨基酸突變,為“WRKYGKK”。依據WRKY特征結構域將41條序列的51個結構域分為3大類和8個亞類型,第I大類為氮端(N)和碳端(C),第II大類下設5個亞類型,以及第III大類;此外,系統進化樹可將51個WRKY結構域分成4大組,Group I [I(N)]、GroupII(IIa+IIb)、GroupIII [I(C)+IIc]、GroupIV(IId+IIe和III),該結果與經典的WRKY家族分類特點基本一致。GO分析發現秦艽WRKY家族成員廣泛參與代謝調控過程。當茉莉酸甲酯(jasmonic acid methyl ester,MeJA)誘導秦艽后,蛋白網絡互作的核心成員GmWRKY1和GmWRKY17表達顯著提高,推測這些WRKY成員與秦艽次生代謝產物的調控關系密切。為進一步開展秦艽WRKY轉錄因子功能研究奠定基礎。

秦艽;WRKY轉錄因子;生物信息學;轉錄組分析;次生代謝產物

秦艽Pall.為“秦藥”大宗道地中藥材,是龍膽科(Gentianaceae)龍膽屬L. 的一種多年生草本植物,以干燥根入藥,具有祛風濕、清溫熱、止痹痛、退虛熱之功效[1-2],其主要的次生代謝產物是以龍膽苦苷(gentiopicroside)為代表的裂環環烯醚萜苷類化合物[3],這類化合物是秦艽發揮臨床療效的物質基礎,也是評價秦艽藥材質量的重要指標[4]。

植物次生代謝產物也稱“植保素”,是植物為了適應環境脅迫而產生的物質[5]。當植物受到外界脅迫后,體內會有一系列信號連鎖反應,其中,轉錄調控是細胞內部調控網絡中的重要一環,而轉錄因子能夠通過協調轉錄酶與DNA模板的相互作用,來激活或抑制下游基因的表達[6]。WRKY是植物中重要的一類轉錄因子,可參與植物的生長發育[7],調控植物次生代謝產物的生物合成[8],以及響應各種生物和非生物脅迫[9]等生物進程。

WRKY轉錄因子調控藥用植物次生代謝產物生物合成是當前研究的熱點之一。研究發現,過表達可促進丹參毛狀根中丹參酮含量的積累,其主要通過激活丹參酮合成途徑酶關鍵基因來提高丹參酮的合成[10]。人參基因可響應病原菌的誘導,同時過表達基因可顯著提高人參皂苷含量[11]。過表達基因可將短小蛇根草中喜樹堿的含量提高3倍,干涉基因則顯著降低喜樹堿的含量[12]。過表達長春花基因可以顯著提高代謝途徑的關鍵酶基因色氨酸脫羧酶和轉錄抑制因子ZCT的表達,并抑制轉錄激活因子ORCA2、ORCA3和CrMYC2,從而使長春花轉基因毛狀根中的蛇根堿含量顯著增加[13]。這些研究都表明WRKY轉錄因子在調控植物代謝產物積累中發揮著重要的作用。

本研究建立在秦艽轉錄組數據的基礎上,對WRKY轉錄因子家族成員保守結構域、系統發育關系、基因表達水平、轉錄因子蛋白互作等方面進行了分析,為后期研究WRKY轉錄因子調控秦艽次生代謝產物的機制奠定基礎。

1 材料

樣品采自甘肅慶陽正寧縣由陜西中醫藥大學李鉑教授鑒定為秦艽Pall.的種子,于溫室撒播種子給予16 h 光照,8 h 黑暗處理,覆膜出苗后培養一年生用于試驗。多通道熒光定量PCR儀qTOWER 2.0(德國耶拿公司)。

2 方法

2.1 秦艽WRKY轉錄因子家族成員的鑒定

研究前期以50 μmol/L檸檬酸銨(ammonium citrate,AC)、200 μmol/L茉莉酸甲酯(jasmonic acid methyl ester,MeJA)為誘導子,水處理為對照分別處理秦艽幼苗6 h,送于深圳千年基因(Macrogen公司,中國)完成秦艽轉錄組測序和分析。本研究基于轉錄組基因注釋結果,初步篩得74條潛在的秦艽基因,逐條進行NCBI-tBlastx序列相似度比對,參數為默認,以判定序列確為WRKY家族。并進一步利用SMART數據庫(http://smart. embl.de/)比較序列家族特征,去掉WRKY結構域不完整的序列,最終篩選得到41條基因。

2.2 GmWRKY保守基序分析、蛋白多重比對以及系統進化樹的構建

首先利用getorf(http://emboss.bioinformatics.nl/ cgi-bin/emboss/getorf)在線分析41條序列的開放閱讀框(open reading frame,ORF)區,并選擇含有WRKY基序的ORF長片段。隨后,通過MEME在線軟件(https://meme-suite.org/meme/)對秦艽WRKY蛋白序列進行保守基序鑒定,參數選擇分析3個保守區。利用MEGAX[14]對秦艽WRKY家族序列的ORF區域進行ClustalW比對,選取保守區,采用鄰接法(neighbor-joining,NJ)進行系統進化樹的構建,Bootstrap檢測值為1000,其他使用參數均為默認值。

2.3 GmWRKY轉錄因子家族成員(GO)注釋及蛋白質網絡互作

利用Omicshare GO富集分析工具建立WRKY家族轉錄組注釋信息。利用String蛋白互作數據庫(http://string-db.org/)對秦艽WRKY蛋白進行互作網絡預測分析,選擇擬南芥為參考物種。并利用Cytoscape 3.8.0[15]對互作結果進行視圖優化。

2.4 GmWRKY轉錄因子家族成員的表達分析

基于轉錄組數據挖掘MeJA、AC和對照組的GmWRKY轉錄因子的FPKM值,用Omicshare熱圖分析工具對FPKM值進行聚類分析和可視化,分析類型選pearson。

對蛋白互作網絡核心成員和進行基因相對表達分析,為內參[16],引物信息見表1。用200 μmol/L MeJA和對照水分別噴施處理秦艽幼苗6 h,收集對照和處理組的葉,重復3組,分別提取RNA,對1%凝膠電泳和Nanodrop檢測合格的RNA用于后續的定量實驗。使用德國耶拿熒光定量梯度PCR儀qTOWER 2.0,反應條件為95 ℃、30 s預變性,95 ℃、5 s,60 ℃、30 s,40個循環。利用2?ΔΔCt計算基因相對表達水平[17]。

表1 實時熒光定量基因表達序列引物

3 結果

3.1 GmWRKY轉錄因子家族成員的鑒定

在秦艽轉錄組數據有74條基因被注釋為秦艽基因,序列最短為209 bp,最長為2671 bp。通過對74條基因逐條進行NCBI-tBlastx比對,初步判定序列均為WRKY家族成員,但因部分序列信息不完整,結構域缺失,結合在線軟件SMART分析結果,去除保守結構域不完整的序列,將明確具有WRKY保守結構域的41條序列列為秦艽WRKY家族成員,并依次命名為GmWRKY1~GmWRKY41,作為后續研究的對象。

3.2 GmWRKY家族成員蛋白質保守區分析

使用在線軟件MEME分析41條秦艽WRKY家族蛋白的保守基序,經分析比對后3個保守區Motif1、Motif2、Motif3的基序長度分別為29、29、50個氨基酸。其中,Motif 1和Motif 3是基因的七肽保守序列,Motif 2是秦艽基因的特征基序(圖1)。41條序列中,其中34條序列包含Motif1和Motif2,這34條序列中又有10條序列含Motif3。此外,GmWRKY7、24、29和31只含有1個保守區Motif1;GmWRKY16、26和30只含有1個保守區Motif2。

圖1 秦艽41個WRKY家族成員蛋白保守基序分析

3.3 蛋白質序列比對及系統進化分析

多重序列比對分析結果發現41個GmWRKY蛋白均含有完整的結構域核心序列(圖2)。大部分序列具有WRKY家族的保守基序“WRKYGQK”,僅GmWRKY27蛋白保守基序中的“Q”變成“K”。依據WRKY家族的結構特征,將41個WRKY蛋白分為3大類,其中10個蛋白為第I類,該類蛋白含有2個“WRKYGQK”序列,分別位于序列的N端和C端,以及C2H2型鋅指結構。第II類有23個蛋白,特點是含有1個WRKY結構域,鋅指結構為C2H2型,該大類可根據序列差異可進一步細化為5個小類,即IIa~IIe。第III類有6個WRKY轉錄因子,含有1個WRKY結構域,鋅指結構為C2HC型。GmWRKY7和GmWRKY16與其他家族成員之間的保守氨基酸位點存在較大差異,分類不明確。

基于MEGAX比對秦艽41個WRKY家族成員的51個WRKY保守域,并構建系統進化樹,其結果與保守結構域分類情況基本一致。系統進化樹可將秦艽的3大類WRKY家族成員分為4大組(圖3),Group1 [I (N)]、Group2 [IIc+I(C)]、Group3(IIa+IIb)、Group4(IId+IIe和III)。其中,Group2與Group1親緣關系較近,聚為一大支。Group3與Group4親緣關系近,被聚為另一大支。對于分類不明確的GmWRKY7和GmWRKY16分別聚類到了I(C)和I(N),與第I類親緣關系較近。

紅色方框表示WRKY核心結構域,紅色方框+ 表示鋅指結構;I類的2個WRKY結構域,分別標識為“N”和“C”端

圖3 WRKY家族的系統發育樹

3.4 GmWRKY轉錄因子家族成員GO注釋

基于轉錄組GO數據庫信息,對秦艽WRKY轉錄因子進行功能注釋。對于大部分秦艽基因在GO的二級分類的生物過程(biological process,BP)、細胞組分(cellular component,CC)和分子功能(molecular function,MF)3個大類中得到注釋(圖4)。在生物過程中細胞過程(cellular process)和代謝過程(metabolic process)富集了40條序列信息。在細胞組分中有40條可富集于細胞(cell)、細胞部分(cell part)和細胞器(organelle)。分子功能中有40條序列可富集于綁定(binding)和核酸結合轉錄因子活性(Nucleic acid binding transcription factor activity)。

圖4 秦艽WRKY家族成員的GO注釋信息

3.5 GmWRKY蛋白互作網絡預測

基于String數據庫,構建了秦艽41條WRKY轉錄因子家族成員的蛋白互作網絡,去除沒有相互關系作用的成員,秦艽WRKY家族中的11個成員間有蛋白質互作的關系(圖5)。可依據相互作用分為3塊互作網絡。其中,GmWRKY1、GmWRKY17、GmWRKY25和GmWRKY15為互作網絡中的4個主要中心節點,相互作用關系最強。這4個主要中心節點與GmWRKY31、GmWRKY13和GmWRKY28之間也存在著相互作用關系。GmWRKY4和GmWRKY16,GmWRKY39、GmWRKY19和GmWRKY41為2塊具有獨立互作的蛋白網絡。

圖5 秦艽WRKY家族成員的蛋白互作網絡

3.6 GmWRKY轉錄因子家族成員表達模式分析

基于秦艽轉錄組數據的FPKM值,對WRKY基因在MeJA和AC誘導后的表達水平差異進行熱圖分析,和對照(CK)相比,41條GmWRKY家族成員受誘導調控顯著,且不同成員對于MeJA和AC誘導響應強度不同(圖6)。基因表達聚類結果顯示,至,和對照相比,MeJA誘導后表達明顯降低。至在MeJA和AC誘導后普遍降低。至在MeJA和AC誘導后普遍升高。至在AC誘導后表達顯著升高。至在AC誘導后表達顯著升高,而對MeJA誘導響應不強。至,和對照相比,MeJA誘導后表達明顯升高而AC誘導后表達普遍降低。

圖6 基于轉錄組FPKM值分析秦艽GmWRKY基因表達差異

以秦艽WRKY家族成員蛋白互作網絡關系最強的4個節點GmWRKY1、GmWRKY17、GmWRKY25和GmWRKY15進行qRT-PCR基因表達分析(圖7)。當秦艽受到MeJA誘導后,葉中表達量顯著升高,其次是,和差異不顯著。

4 討論

WRKY轉錄因子家族因N端包含高度保守的WRKY結構域“WRKYGQK”而得名,該結構域能專一地與靶基因啟動子區中的W盒(TTGACC/T)序列結合,激活或抑制基因的表達,進而調控植物的生長代謝,參與響應應激等各種生理活動[18]。秦艽41條GmWRKY家族成員WRKY結構域保守,大部分為“WRKYGQK”,僅GmWRKY27突變為“WRKYGKK”。在丹參、水稻、擬南芥等物種中也發現了類似的突變體,這種突變的產生也可能與物種長期進化有關[19]。有研究發現突變基序中谷氨酰胺(glutamine,G)氨基酸的突變會降低保守基序與DNA結合的能力[20]。對于WRKY結構域突變體的進化事件、以及突變體參與植物代謝調控的機制還需要進一步的研究。

不同字母表示差異顯著P<0.05

Eulgem等[21]依據WRKY家族成員的結構特點將其分為3大類型I(N端+C端)、II(a、b、c、d、e)和III。本研究對GmWRKY蛋白家族成員進行多重比對、保守基序分析結果發現秦艽WRKY家族成員也符合該分類特點,這些結果表明,不同物種的WRKY結構保守,功能可能也會比較相似。此外,系統進化樹依據親緣關系將GmWRKY保守域分為2大支,4大組。II類型的WRKY成員不單獨分為一類,說明其并不是單系起源,且具有兩個WRKY結構域基因的C端可能是具有一個WRKY域編碼基因的祖先[22]。有研究發現許多來自不同系統發育群的基因參與植物響應同樣的脅迫應激,說明WRKY家族可廣泛參與植物的各種生理代謝調控網絡,對WRKY家族的起源、進化關系的了解可以為其功能的研究提供有價值的信息[22]。

隨著轉錄組測序技術的發展,越來越多藥用植物的WRKY轉錄因子家族被挖掘和研究。WRKY家族在陽春砂和丹參中被挖掘,并預測了參與次生代謝產物合成相關的WRKY轉錄因子[23-24]。在紅豆杉中挖掘了61條WRKY家族,并發現過表達WRKY成員能夠提高目標藥用成分的含量[25]。本研究基于秦艽轉錄組挖掘了41條WRKY轉錄因子,通過GO數據分析發現這些家族成員中有40條具有與核酸結合的轉錄因子活性,并廣泛參與代謝調控過程。基于轉錄組數據的基因表達熱圖顯示,WRKY成員對MeJA和AC的響應強度不同,推測不同成員在植物體內參與代謝調控的方式可能存在差異。此外,網絡互作中的核心成員在受到MeJA誘導后基因表達均有所升高,其中和升高最為顯著,推測這些WRKY成員與秦艽次生代謝產物的調控關系密切。對秦艽WRKY轉錄因子家族的分析將為進一步研究WRKY轉錄因子參與調控裂環環烯醚萜苷類代謝產物合成途徑奠定基礎。

利益沖突 所有作者均聲明不存在利益沖突

[1] 中國藥典 [S]. 四部. 2020: 45.

[2] 胡本祥, 彭亮, 楊冰月, 等. “秦藥”的現代研究概況 [J]. 中草藥, 2018, 49(21): 4949-4959.

[3] 陳千良, 石張燕, 涂光忠, 等. 陜西產秦艽的化學成分研究 [J]. 中國中藥雜志, 2005, 30(19): 1519-1522.

[4] 馬麗娜, 張鐵軍, 田成旺, 等. 大孔樹脂分離純化川西獐牙菜中環烯醚萜苷類和酮類成分的工藝研究 [J]. 中草藥, 2010, 41(2): 227-231.

[5] 郭蘭萍, 周良云, 康傳志, 等. 藥用植物適應環境脅迫的策略及道地藥材“擬境栽培” [J]. 中國中藥雜志, 2020, 45(9): 1969-1974.

[6] Lee T I, Young R A. Transcription of eukaryotic protein-coding genes [J]., 2000, 34: 77-137.

[7] Yang Y, Chi Y J, Wang Z,. Functional analysis of structurally related soybean GmWRKY58 and GmWRKY76 in plant growth and development [J]., 2016, 67(15): 4727-4742.

[8] Schluttenhofer C, Yuan L. Regulation of specialized metabolism by WRKY transcription factors [J]., 2015, 167(2): 295-306.

[9] Wani S H, Anand S, Singh B,. WRKY transcription factors and plant defense responses: Latest discoveries and future prospects [J]., 2021, 40(7): 1071-1085.

[10] Deng C P, Hao X L, Shi M,. Tanshinone production could be increased by the expression of SmWRKY2inhairy roots [J]., 2019, 284: 1-8.

[11] Yao L, Wang J, Sun J C,. A WRKY transcription factor, PgWRKY4X, positively regulates ginsenoside biosynthesis by activating squalene epoxidase transcription in[J]., 2020, 154: 112671.

[12] Hao X L, Xie C H, Ruan Q Y,. The transcription factor OpWRKY2positively regulates the biosynthesis of the anticancer drug camptothecin in[J]., 2021, 8(1): 7.

[13] Suttipanta N, Pattanaik S, Kulshrestha M,. The transcription factor CrWRKY1positively regulates the terpenoid indole alkaloid biosynthesis in[J]., 2011, 157(4): 2081-2093.

[14] Kumar S, Stecher G, Li M,. MEGA X: Molecular evolutionary genetics analysis across computing platforms [J]., 2018, 35(6): 1547-1549.

[15] Shannon P, Markiel A, Ozier O,. Cytoscape: a software environment for integrated models of biomolecular interaction networks [J]., 2003, 13(11): 2498-2504.

[16] He Y H, Yan H L, Hua W P,. Selection and validation of reference genes for quantitative real-time PCR in[J]., 2016, 7: 945.

[17] Erickson H S, Albert P S, Gillespie J W,. Assessment of normalization strategies for quantitative RT-PCR using microdissected tissue samples [J]., 2007, 87(9): 951-962.

[18] Bakshi M, Oelmüller R. WRKY transcription factors: Jack of many trades in plants [J]., 2014, 9(2): e27700.

[19] 田云, 盧向陽, 彭麗莎, 等. 植物WRKY轉錄因子結構特點及其生物學功能 [J]. 遺傳, 2006, 28(12): 1607-1612.

[20] Duan M R, Nan J, Liang Y H,. DNA binding mechanism revealed by high resolution crystal structure ofWRKY1protein [J]., 2007, 35(4): 1145-1154.

[21] Eulgem T, Rushton P J, Robatzek S,. The WRKY superfamily of plant transcription factors [J]., 2000, 5(5): 199-206.

[22] Zhang Y J, Wang L J. The WRKY transcription factor superfamily: Its origin in eukaryotes and expansion in plants [J]., 2005, 5: 1.

[23] He X Y, Wang H, Yang J F,. RNA sequencing onLour. induced by MeJA identifies the genes of WRKY and terpene synthases involved in terpene biosynthesis [J]., 2018, 61(2): 91-102.

[24] Yu H Z, Guo W L, Yang D F,. Transcriptional profiles offamily genes and their putative roles in the biosynthesis of tanshinone and phenolic acids in[J]., 2018, 19(6): E1593.

[25] Zhang M, Chen Y, Nie L,. Transcriptome-wide identification and screening of WRKY factors involved in the regulation of taxol biosynthesis in Taxus chinensis [J]., 2018, 8(1): 5197.

Bioinformatics analysis of WRKY transcription factor family of

HE Yi-han, YIN Yang-yang, HU Wei, LI Bo, SUN Xiao-chun, WANG Nan, HUANG Wen-jing, YUE Zheng-gang

State Key Laboratory of Research & Development of Characteristic Qin Medicine Resources (Cultivation), Shaanxi Collaborative Innovation Center of Chinese Medicinal Resources Industrialization, Shaanxi university of Chinese, Xianyang 712046, China

To analyze theWRKY family of(GmWRKY) using biotransformation techniques in order to provide information for further research on the molecular regulation mechanism of GmWRKY transcription factor on the biosynthesis of the secondary metabolites of.Based on the annotation results of WRKY family members, 41 WRKY sequences were further screened by NCBI-tBLastx and SMART domains for bioinformatics analysis. Based on getorf prediction of open reading frame (ORF) region, MEME conservative motif analysis and ClustalW alignment were performed on ORF fragments of 41 sequences, and the neighbor-joining (NJ) phylogenetic tree was constructed by MEGAX. The enrichment analysis of 41 WRKY members was performed using Omicshare GO. And prediction analysis of WRKY member protein interaction network was performed by String. The expression of WRKY members was analyzed by heat map based on the FPKM value of the transcriptome, and the relative gene expression of the core members of GmWRKY protein interaction was analyzed.Most of the 41 GmWRKY members had a conserved motif of "WRKYGQK", and only THE GmWRKY27 conserved motif had an amino acid mutation of "WRKYGKK". According to the WRKY characteristic domain, 51 domains of 41 sequences were divided into three categories and eight subtypes. The first category consisted of nitrogen terminal (N) and carbon terminal (C), the second category consisted of 5 subtypes and the third category. Phylogenetic tree can divide 51 WRKY domains into four groups: GroupI [I(N)], GroupII (IIa+IIb), GroupIII [I(C)+IIc], GroupIV (IId+IIe and III). These results were consistent with the classic WRKY family classification. GO analysis revealed that members of the GmWRKY family were extensively involved in metabolic regulation. The expression ofand, the core members of protein network interaction, were significantly increased whenwas induced by MeJA, suggesting that these WRKY members were closely related to the regulation of secondary metabolites of.This study can lay a foundation for further research on the function of GmWRKY transcription factors..

Pall.; WRKY transcription factors; bioinformatics analysis; transcriptome analysis; secondary metabolites

R286.12

A

0253 - 2670(2022)23 - 7499 - 08

10.7501/j.issn.0253-2670.2022.23.021

2022-05-06

國家自然科學基金項目(81903753);陜西中醫藥大學科學基金項目(2017PY06)

何懿菡,講師,研究方向為藥用植物資源及次生代謝產物調控。Tel: (029)38183207 E-mail: annaaid@126.com

通信作者:岳正剛,副教授,研究方向為中藥藥效物質基礎。Tel: (029)38183301 E-mail: liuxingjian1981@163.com

[責任編輯 時圣明]

猜你喜歡
分析
禽大腸桿菌病的分析、診斷和防治
隱蔽失效適航要求符合性驗證分析
電力系統不平衡分析
電子制作(2018年18期)2018-11-14 01:48:24
電力系統及其自動化發展趨勢分析
經濟危機下的均衡與非均衡分析
對計劃生育必要性以及其貫徹實施的分析
現代農業(2016年5期)2016-02-28 18:42:46
GB/T 7714-2015 與GB/T 7714-2005對比分析
出版與印刷(2016年3期)2016-02-02 01:20:11
中西醫結合治療抑郁癥100例分析
偽造有價證券罪立法比較分析
在線教育與MOOC的比較分析
主站蜘蛛池模板: 亚洲成人在线网| 91成人在线观看视频| 色有码无码视频| 亚洲三级视频在线观看| 五月天香蕉视频国产亚| 亚洲国产成人精品青青草原| 亚洲天堂网在线播放| 欧美69视频在线| 97狠狠操| 亚洲日本www| 亚洲第一成年网| 日本草草视频在线观看| 国产视频a| 四虎影视永久在线精品| 91精品啪在线观看国产| 在线a网站| 丁香婷婷激情网| 久久黄色毛片| 精品福利网| 日韩欧美色综合| 中文字幕亚洲乱码熟女1区2区| 成人在线观看一区| 性视频久久| 久久精品人妻中文系列| 国产色爱av资源综合区| AV不卡无码免费一区二区三区| 熟女成人国产精品视频| 精品视频第一页| 欧美精品在线观看视频| 国产女人18水真多毛片18精品 | 无码在线激情片| 欧美日韩在线第一页| 人人看人人鲁狠狠高清| 中文字幕精品一区二区三区视频| 91免费在线看| 国产成人精品一区二区| 亚洲高清中文字幕| 一级毛片在线免费看| 久久综合九色综合97网| 手机在线看片不卡中文字幕| 欧美综合成人| 中文字幕在线看视频一区二区三区| 免费一级全黄少妇性色生活片| 不卡无码网| 欧美日韩资源| 国产精品第一区| 国产精品v欧美| 午夜性刺激在线观看免费| 国产美女视频黄a视频全免费网站| 欧美成人一区午夜福利在线| 亚洲色欲色欲www网| 国产极品美女在线| 99在线视频免费观看| 波多野结衣无码AV在线| 国产亚洲精品无码专| 婷婷色婷婷| 亚洲一级毛片在线观播放| 国产在线第二页| 亚洲无码高清视频在线观看 | 久久婷婷综合色一区二区| 欧美全免费aaaaaa特黄在线| 国产伦片中文免费观看| 手机永久AV在线播放| 在线免费看片a| 欧美中文字幕一区二区三区| 久久香蕉国产线看观看精品蕉| 国产精品青青| 人人91人人澡人人妻人人爽| 日韩福利在线视频| 香蕉网久久| 性欧美在线| 91区国产福利在线观看午夜| 中文字幕 欧美日韩| 亚欧美国产综合| 91精品国产自产在线观看| 成人国产一区二区三区| 久久国产V一级毛多内射| 91美女视频在线| 久久综合干| 91小视频在线观看| 久久99国产视频| a级毛片免费看|