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

陸地棉Nudix 基因家族的全基因組鑒定及表達分析

2021-04-14 06:56:58竇玲玲孫亞如趙琴田瑞潔康洋洋朱怡然楊蕾蕾王彩虹馮宇王文博肖光輝
棉花學報 2021年2期
關鍵詞:分析

竇玲玲,孫亞如,趙琴,田瑞潔,康洋洋,朱怡然,楊蕾蕾,王彩虹,馮宇,王文博,肖光輝*

(1. 咸陽師范學院化學與化工學院,陜西 咸陽712000;2. 河南大學生命科學學院,河南 開封475004;3. 陜西師范大學生命科學學院,西安710119)

Nudix 是一類廣泛存在于真核和原核生物(細菌、古細菌和病毒)中,水解RNA 帽子結構和包括核苷糖類、二核苷多聚磷酸鹽、三磷酸核苷在內的多種有機焦磷酸鹽的水解酶[1],在DNA 損傷修復和逆境脅迫生理方面發揮著重要作用。 Nudix 水解酶包含典型的Nudix 基序“GX5EX7REUXEEXGU”, 其中X 代表任意氨基酸,U 代表疏水性氨基酸。

Nudix 水解酶參與修復鳥嘌呤核苷酸(Guanosine triphosphate,GTP)的氧化損傷。 GTP在活性氧的作用下易氧化形成8- 氧-7,8 二氫鳥嘌呤(8-oxo-7,8-dihydroguanine, 8-oxo-GTP),而DNA 復制時8-oxo-GTP 可以將G:C 顛換突變為T:A。 研究表明,大腸桿菌Nudix 水解酶(MutT)可以水解8-oxo-GTP, 防止DNA 復制過程中顛換的產生以減少核酸損傷和突變[2]。 另外,Nudix水解酶在DNA 修復中的作用在原核生物和哺乳動物的研究中也有報道[3]。

Nudix 水解酶參與植物的生物和非生物脅迫應答反應。 在擬南芥中,AtNudix6 正向調控水楊酸(Salicylic acid,SA)誘導的病原菌響應基因NPR1(Non-expressor of Pathogenesis Related gene 1)的表達,從而提高植物對病原菌的防御作用[4-5]。另外,AtNudix7 在受到臭氧和病原菌脅迫時快速上調表達[6]。 Nudix 還參與了包括黃素代謝、輔酶A 分解代謝等在內的細胞代謝過程[7-8]。 由此表明,Nudix 水解酶參與調控植物逆境脅迫應答反應。

近年來,隨著測序技術的發展和進步,愈來愈多的物種完成了基因組的測序工作。 基于全基因組的Nudix 基因家族分析也在多個物種中有見報道,包括擬南芥[7],大腸桿菌[9],釀酒酵母[10]和人類[11]等。 但棉花中Nudix 基因家族還未被挖掘,因此,我們開展了Nudix 基因家族在陸地棉中的系統進化和轉錄組表達分析,初步揭示了該基因家族在陸地棉纖維發育過程中的作用。

1 材料與方法

1.1 陸地棉Nudix 基因家族的鑒定

從CottonFGD 數據庫[12]下載陸地棉(Gossypium hirsutum acc.TM-1,ZJU_v2.1)和雷蒙德氏棉(G. raimondii,JGI_v2.1)和亞洲棉(G. arboreum,CRI_v3.0) 基因組序列。 從擬南芥數據庫(The arabidopsis information resource,TAIR)下載擬南芥AtNudix 蛋白序列信息[13]。從Pfam 數據庫[14]下載Nudix 蛋白的種子文件PF00293,并利用HMMER 3.0 軟件的hmmsearch 搜索程序鑒定含有Nudix 保守結構域的陸地棉蛋白序列。 將所有獲得的陸地棉Nudix 蛋白序列提交NCBI 網站的保守結構域數據庫 (Conserved domain database,CDD)進行保守結構域驗證。

1.2 基因染色體位置、結構和蛋白序列的理化性質分析

根據陸地棉基因組注釋文件(General feature format,GFF),利用MapChart 軟件[15]對GhNudix基因的染色體位置進行可視化處理。根據GFF 文件分析基因“外顯子-內含子” 結構特征, 并用GSDS 2.0 在線軟件[16]進行基因結構可視化處理。

通過ExPASy 在線軟件對GhNudix 蛋白序列的等電點(Isoelectric point,pI),相對分子質量(Relative molecular mass,Mr), 親水性平均系數(Grand average of hydropathicity,GRAVY) 等理化性質進行預測[17]。 根據CDD 保守結構域分析結果,利用WebLogo 在線軟件[18]繪制陸地棉GhNudix 蛋白保守結構域序列情況。 利用SWISS-MODEL 軟 件[19]展 示GhNudix 保 守 結 構域的空間結構特征。

1.3 蛋白序列的多重比對和系統進化樹構建

使用ClustalX 2.0 對GhNudix、GaNudix、Gr-Nudix 和AtNudix 蛋白序列進行多重序列比對[20],并使用MEGA 7.0 軟件[21]的鄰近連接法(Neighbour-joining method)構建無根系統發育樹,Bootstrap 設置為1 000 次。利用MEME 在線軟件[22]分析GhNudix 蛋白序列保守基序, 參數設置如下:最大發現數目為5, 基序最長為50 個氨基酸;功能域分布類型設置為zoops, 即功能域在每條序列中只出現1 次。

1.4 陸地棉GhNudix 基因復制分析

使用MCScanX 軟件的duplicate_gene_classifier 軟件包對GhNudix,GrNudix 和GaNudix 的復制類型和共線性區段進行分析[23],利用Circos軟件[24]對GhNudix 基因的共線性進行可視化處理。 利用KaKs_Calculator 2.0 軟件分析陸地棉GhNudix 復制基因對的同義替換率(Synonymous,Ks) 和 非 同 義 替 換 率(Non-synonymous,Ka)。 Ka/Ks<1 表示純化選擇,Ka/Ks=1 表示中性選 擇,Ka/Ks>1 表 示 正 向 選 擇[25]。 利 用 公 式“t=Ks/2r” 計算基因的分歧時間, 其中r=2.6×10-9,代表中性替換率[26]。

1.5 陸地棉GhNudix 基因的啟動子及表達分析

利用本地BLAST 軟件調取基因起始密碼子上游2 000 bp 的序列作為陸地棉GhNudix 基因的啟動子序列。 利用在線軟件Plant cis-acting regulatory element (Plant CARE)數 據 庫[27]分 析GhNudix 基因啟動子的順式作用元件。

為了進一步分析GhNudix 基因的表達情況,從National Center for Biotechnology Information(NCBI)的序列閱讀檔案(Sequence read archive,SRA)數據庫下載了陸地棉開花當天(0 DPA,Day post anthesis)、5 DPA、10 DPA、20 DPA 和25 DPA的纖維和不同組織(根、莖、葉、萼片、花瓣、雌蕊、雄蕊、花托和0 DPA 胚珠)轉錄組數據[28],項目登錄號為SRA180756。 采用FPKM (Fragments per kilobase of transcript per million mapped reads,每百萬片段中來自某基因每千堿基長度的片段數)法對表達reads 進行歸一化處理[29],進一步利用R 語言的pheatmap 軟件包[30]繪制基因表達量熱圖。

2 結果與分析

2.1 陸地棉GhNudix 基因家族的鑒定

以陸地棉基因組(ZJU_v2.1)為參考序列,利用HMMER 3.0 軟件搜索PF00293 蛋白序列種子文件,通過進一步分析保守結構域,最終在陸地棉基因組中鑒定到76 個GhNudix 水解酶家族成員。陸地棉GhNudix 水解酶家族的蛋白序列長度范 圍 是60 (GhNudix5) 到776 個 氨 基 酸(Gh-Nudix65);蛋白質分子量為6.46~139.61 kDa。根據pI 分析,64 個GhNudix 蛋白的pI<7.0 (平均值5.75), 為酸性蛋白質;12 個GhNudix 蛋白的pI>7.0 (平均值8.32), 為堿性蛋白。 所有Gh-Nudix 蛋白的親水性平均系數(GRAVY)都小于0, 表明陸地棉GhNudix 蛋白全部都是親水性蛋白。 從亞洲棉和雷蒙德氏棉中分別鑒定出來35個GaNudix 和36 個GrNudix 蛋白, 其理化參數詳見附表1。

2.2 陸地棉GhNudix 基因的染色體分布

根據陸地棉基因組GFF 注釋文件分析Gh-Nudix 基因在染色體上的相對位置, 從染色體At_01 到Dt_13, 依次將該基因家族成員命名為GhNudix1 到GhNudix75; 還有一個GhNudix 基因位于scaffold (scaffold6-1_subseq_1_118545_obj_D02),命名為GhNudix76。 75 個GhNudix 基因在陸地棉26 條染色體上的分布是不均衡的,其中At 亞基因組含有43 個GhNudix 基因,Dt亞基因組上有33 個GhNudix 基因(圖1)。 由此說明,GhNudix 基因家族在At 亞組和Dt 亞組經歷了非對稱進化。

2.3 陸地棉Nudix 蛋白的多重序列比對和保守結構域分析

據文獻報道, 典型的 Nudix 基序為“GX5EX7REUXEEXGU”序列[7]。通過利用ClustalX 2.0 進行多重序列比對, 大部分陸地棉GhNudix蛋白序列中具有典型的Nudix 基序結構(圖2),而部分蛋白Nudix 基序的“GU”轉變為了“AU”,這可能是陸地棉GhNudix 基因在長期進化過程中的突變引起的。利用SWISS-MODEL 在線軟件分析Nudix 的基序, 該基序在空間結構上形成1個α-螺旋結構。

2.4 陸地棉GhNudix 基因結構和基序分析

圖1 GhNudix 基因在陸地棉染色體上的分布Fig. 1 Chromosomal distribution of GhNudix genes in G. hirsutum

圖2 陸地棉GhNudix 蛋白的Nudix 基序的空間構象及WebLogo 統計Fig. 2 Spatial conformation and WebLogo statistics of the typical motif of Nudix in G. hirsutum

通過比較GhNudix 基因的外顯子-內含子可以看出(圖3),基因外顯子最少的只有1 個(GhNudix5), 而最多有14 個 (GhNudix28 和GhNudix65)。 通過MEME 在線軟件分析陸地棉GhNudix 蛋白序列的保守motif 分析,共發現了5個保守的motif(圖3),其中motif 1 在所有的陸地棉Nudix 蛋白序列中都存在。

2.5 GhNudix 蛋白的系統發育分析

為分析陸地棉Nudix 基因家族的進化關系,將76 個GhNudix、35 個GaNudix、36 個GrNudix和26 個AtNudix 的全長蛋白序列導入MEGA 7.0 軟件,采用鄰接法構建系統發育樹,并根據擬南芥Nudix 蛋白家族的分類關系對陸地棉、亞洲棉和雷蒙德氏棉Nudix 蛋白進行分類(圖4)。 聚類結果顯示, 陸地棉、 亞洲棉和雷蒙德氏棉的Nudix 蛋白被分為7 個亞組,分別命名為Ⅰ、Ⅱ、Ⅲ、Ⅳ、Ⅴ、Ⅵ和Ⅶ亞組,分別含有18、22、8、6、4、8和10 個GhNudix 家族成員。

2.6 陸地棉GhNudix 基因的復制事件

圖3 陸地棉GhNudix 基因的系統進化樹、基因結構和蛋白序列保守基序Fig. 3 Phylogenetic tree, gene structure and conserved motif of GhNudix protein sequences in G. hirsutum

為了研究Nudix 基因在染色體區段上的復制關系,利用MCScanX 軟件分別分析了陸地棉、亞洲棉和雷蒙德氏棉Nudix 基因的進化情況。 分析結果顯示,65 個GhNudix 基因來源于片段重復(Segmental duplication),占GhNudix 基因家族成員的85.5%;剩余的GhNudix 基因在染色體上呈分散排列(圖5,附表2)。 我們進一步分析了陸地棉片段重復基因對的Ka、Ks和Ka/Ks值。 由于Ks值不易受進化選擇影響,因此常被用于估算進化的分歧時間[31-33]。 陸地棉Nudix 復制基因對的Ks值為0.032~0.160(附表3),由此推斷基因復制事件發生于0.62 百萬~31.45 百萬年前。 Ka/Ks值分析結果表明,絕大多數陸地棉片段重復基因對的Ka/Ks值小于1, 由此推斷大部分的陸地棉Nudix 基因在進化的過程中受到了純化選擇。 亞洲棉Nudix 基因家族成員擴增主要來自染色體的片段重復(54.3%);雷蒙德氏棉Nudix 基因家族成員擴增主要源于串聯重復(Tandem duplication),占總GrNudix 成員數目的52.8%,片段重復占比22.2%。 亞洲棉和雷蒙德氏棉的Nudix 復制基因對的Ka/Ks值均小于1, 說明亞洲棉和雷蒙德氏棉的Nudix 基因在進化過程中受到純化選擇。 綜上所述,棉花Nudix 基因家族成員的擴增主要來源于染色體片段重復,且在進化過程中主要受到了純化選擇的作用。

圖4 陸地棉、亞洲棉、雷蒙德氏棉和擬南芥的Nudix 蛋白的系統進化樹Fig. 4 Phylogenetic tree of Nudix protein sequences in G. hirsutum, G. arboreum, G. raimondii and A. thaliana

2.7 GhNudix 基因在陸地棉不同組織的表達分析

圖5 陸地棉、亞洲棉和雷蒙德氏棉Nudix 基因的共線性關系Fig. 5 The synteny relationships of Nudix genes among G. hirsutum, G. arboreum and G. raimondii

由于基因的表達和功能的執行具有組織特異性,為了分析陸地棉GhNudix 基因在不同組織中的表達情況,利用轉錄組數據分析了GhNudix基因在棉花根、莖、葉、萼片、花瓣、雌蕊、雄蕊、花托和胚珠(0 DPA)中的表達情況(圖6)。 根據GhNudix 基因在不同組織中的表達變化可以將其分為5 組:第一組,主要在花瓣、雌蕊和雄蕊中具有較高的表達量, 包括編號為11、13、14、16、20、23、31、39、48、49、58、60、65、67、68 和72 的Gh-Nudix 基因,主要集中在系統發育樹的Clade Ⅱ;第二組,主要在花托和0 DPA 胚珠中具有較高的表達量, 包括編號為4、5、19、26、27、29、44、45、54、62、63、66 和69 的GhNudix 基因, 主要來源于系統發育樹的Clade Ⅰ和Clade Ⅲ; 第三組,編號為1、7、8、46 和76 的GhNudix 基因在萼片中具有較高的表達, 主要集中在系統發育樹的Clade Ⅶ;第四組,編號為10、17、41 和52 的Gh-Nudix 基因在根中具有較高的表達, 主要集中在系統發育樹的Clade Ⅱ和Clade Ⅳ; 第五組在花托和莖的表達量較高, 包括編號為2、3、6、9、18、22、24、25、32、35、37、38、40、42、43、47、51、55、56、57、61、70、71、73、74 和75 的GhNudix 基因,主要集中于系統發育樹的Clade Ⅱ。 GhNudix 基因在棉花根、莖、葉、萼片、花瓣、雌蕊、雄蕊、花托和胚珠(0 DPA) 的轉錄組數據表達情況可以看出,GhNudix 基因的表達具有組織特異性, 其中Clade Ⅶ在萼片中特異性表達,Clade Ⅱ在根、花瓣、雌蕊和雄蕊中均有較高表達。

2.8 GhNudix 基因在棉纖維發育過程中的表達分析

圖6 GhNudix 基因在陸地棉不同組織的表達熱圖Fig. 6 Expression heatmap of GhNudix genes in different cotton tissues

為了進一步分析陸地棉GhNudix 基因在棉纖維發育過程中的作用,利用不同棉纖維發育階段的轉錄組數據分析了GhNudix 基因的表達情況(圖7)。棉纖維從胚珠表皮細胞發育而來,迅速伸長形成單細胞表皮毛[34]。 棉纖維發育可以分為4 個相互重疊的過程,包括起始期、伸長期、次級細胞壁加厚期和成熟期[35]。 0~5 DPA 是棉纖維發育的起始期; 棉纖維發育起始期至20 DPA是纖維發育的快速伸長階段[36],15~40 DPA 是棉纖維細胞壁的纖維素合成時期,該時期主要是促進次級細胞壁的加厚[37]。 根據GhNudix 基因在棉纖維發育過程中的表達情況分為3 個亞組(圖7)。 第一組,44 個GhNudix 基因在棉纖維發育的起始階段(0 DPA 和5 DPA)具有較高的表達量,主要集中在系統發育樹的Clade Ⅰ; 第二組,編號為1、4、9、13、20、22、41、47、55、58、69、71 和74 的GhNudix 基因在棉纖維發育的10 DPA 和20 DPA 表達量較高, 該階段是棉纖維快速伸長期,這些基因主要集中在系統發育樹的Clade Ⅱ;第 三 組, 編 號 為6、37、39、43、51 和68 的Gh-Nudix 基因在20 DPA 和25 DPA 棉纖維發育的階段表達量較高,該階段是棉纖維細胞的次級細胞壁加厚期,這些基因主要集中在系統發育樹的Clade Ⅱ。 由此推斷,GhNudix 基因特別是CladeⅠ在棉纖維發育的起始階段具有重要作用,CladeⅡ亞組的GhNudix 基因在棉纖維的伸長階段和細胞壁次生壁加厚過程中發揮重要作用。

2.9 GhNudix 基因啟動子激素相關的順式作用元件分析

圖7 GhNudix 基因在棉纖維發育過程中的表達熱圖Fig. 7 Expression heatmap of GhNudix genes during cotton fiber developmental stages

生長素(Auxin,IAA)和赤霉素(Gibberellin acid,GA) 是促進棉纖維伸長生長最顯著的促進劑[38-39],因此我們分析了GhNudix 基因啟動子與激素相關的順式作用元件。 對76 個陸地棉Nudix 基因順式作用元件分析表明,36 個GhNudix 基因的啟動子包含1~3 個赤霉素響應元件GARE (Gibberellin-responsive element,序列為TCTGTTG),P-box (CCTTTTG) 和TATC-box(TATCCCA), 其 中13 個 屬 于 系 統 發 育 樹 的Clade Ⅰ,12 個聚類于系統發育樹的Clade Ⅱ;33個GhNudix 基因的啟動子包含1~2 個生長素響應元件,AuxRR (GGTCCAT) 和TGA-element(AACGAC), 其中8 個GhNudix 基因聚類于系統發育樹的Clade Ⅰ,16 個GhNudix 基因聚類于Clade Ⅱ(詳見附表4)。

3 討論

Nudix 水解酶廣泛存在于真核生物、細菌、古細菌和病毒中, 其主要作用是水解核苷二磷酸(Nucleoside diphosphate,NDP) 為 核 苷 一 磷 酸(Nucleoside monophosphate,NMP)[40]。 其中,大腸桿菌基因組有13 個Nudix 水解酶基因[9],釀酒酵母菌基因組有6 個Nudix 水解酶基因[10],人類基因組有24 個Nudix 水解酶基因[11],擬南芥基因組有27 個Nudix 水解酶基因[7]。 本研究從雷蒙德氏、亞洲棉和陸地棉中分別鑒定到了36、35 和76個Nudix 基因, 由此可以看出棉花基因組中Nudix 基因的數目遠高于其他物種。 大量的研究表明,基因復制是基因家族成員擴增的主要方式之一[41-43]。 通過分析陸地棉、雷蒙德氏棉和亞洲棉Nudix 基因復制可以看出, 陸地棉Nudix 基因的85.5%來源于片段重復,亞洲棉54.3%的GaNudix基因來源于片段重復,雷蒙德氏棉GrNudix 基因串聯重復和片段重復分別占52.8%和22.2%,由此證明片段重復是棉花Nudix 基因家族擴增的主要來源。

棉花基因組進化分析表明,亞洲棉和雷蒙德氏棉分別在115 百萬~146 百萬年和13 百萬~20 百萬年前經歷了兩輪的全基因組復制[44-45]。 在1 百萬~2 百萬年前,雷蒙德氏棉(DD,2n=2x=26)作為父本與母本亞洲棉(AA,2n=2x=26)雜交形成了異源四倍體的陸地棉 (AD1,2n=4x=52)[46]。 本文研究表明,亞洲棉Nudix 基因的片段重復發生時間是14.13 百萬~36.96 百萬年前,雷蒙德氏棉Nudix 基因片段重復發生的時間是22.29 百萬~76.10 百萬年前,從而推測雷蒙德氏棉和亞洲棉Nudix 基因家族的成員擴增于棉花基因組的兩輪全基因組復制事件[44,45]。 根據陸地棉Nudix 基因片段重復發生的時間范圍是0.62百萬~31.45 百萬年前,平均時間4.92 百萬年前,由此推斷造成陸地棉Nudix 基因家族成員擴增的片段重復事件主要發生在亞洲棉和雷蒙德氏棉雜交之前。

通過分析GhNudix 基因在根、莖、葉、萼片、花瓣、雌蕊、雄蕊、花托和胚珠(0 DPA)不同組織的轉錄組數據表明,GhNudix 基因在陸地棉中的表達具有組織特異性 (圖6)。 進一步分析Gh-Nudix 基因在棉纖維發育過程中的表達可以發現,在棉纖維發育的起始和伸長階段(0 DPA,5 DPA 和10 DPA), 聚類于系統發育樹的Clade I 和Clade II亞組的GhNudix 基因表達量較高(圖7)。生長素和赤霉素是2 種重要的促進細胞伸長的植物激素,生長素通過促進細胞膜H+轉移,使細胞壁酸化,從而促進細胞伸長生長[47];赤霉素通過誘導木葡聚糖內轉糖苷酶分泌而提高細胞壁的延展性促進細胞伸長生長[48]。 通過分析GhNudix 基因啟動子激素相關的順式作用元件分析發現,大部分的Clade Ⅰ和Clade Ⅱ亞組的基因啟動子含有赤霉素和生長素相關的順式作用元件(詳見附表4)。GhNudix 基因特別是系統發育樹Clade Ⅰ和Clade Ⅱ亞組的基因在棉纖維發育的起始和伸長階段上調表達,且大部分具有響應生長素和赤霉素的順式作用元件,由此推斷陸地棉Nudix 基因在棉纖維發育過程中具有重要作用。

4 結論

本文從陸地棉基因組共鑒定到76 個Gh-Nudix 基因。經聚類分析將其分為7 個亞組,基因復制分析表明,片段重復是導致陸地棉GhNudix基因家族成員擴增的主要原因。 根據轉錄組數據,GhNudix 基因的表達不僅有組織特異性還有時間特異性。 通過分析棉纖維發育不同階段的轉錄組數據可以看出,GhNudix 基因主要在纖維發育的起始和伸長階段高水平表達,其中以系統發育樹Clade Ⅰ和Clade Ⅱ亞組的成員為主; 啟動子分析結果表明,這2 個亞組的GhNudix 成員大部分含有生長素和赤霉素的順式作用元件,由此推測Clade Ⅰ和Clade Ⅱ亞組的GhNudix 基因參與陸地棉纖維發育的起始和伸長。本研究通過對陸地棉GhNudix 基因的全基因組鑒定、 進化和表達分析, 為深入解析GhNudix 基因的功能奠定了基礎。

附表:

附表詳細內容參見http://journal.cricaas.com.cn

附表1 3 種棉花的Nudix 蛋白理化性質的詳細參數

Table S1 Detailed physicochemical parameters of the Nudix proteins from three cotton species

附表2 陸地棉、亞洲棉和雷蒙德氏棉Nudix 基因的復制事件

Table S2 Gene duplication events of Nudix genes in G.hirsutum,G.arboreum and G.raimondii

附表3 3 種棉花Nudix 片段重復基因對的KaKs值和分歧時間

Table S3 KaKsanalysis and divergence time of segmentally duplicated Nudix gene pairs of three cotton species

附表4 GhNudix 基因啟動子植物激素相關的順式元件

Table S4 Cis-elements related to plant hormones in the promoters of GhNudix genes

猜你喜歡
分析
禽大腸桿菌病的分析、診斷和防治
隱蔽失效適航要求符合性驗證分析
電力系統不平衡分析
電子制作(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例分析
偽造有價證券罪立法比較分析
主站蜘蛛池模板: 亚洲综合一区国产精品| 亚洲天堂视频在线观看| 国产成人精品亚洲77美色| 欧美劲爆第一页| 中文字幕在线一区二区在线| 日本高清有码人妻| 日本欧美一二三区色视频| 亚洲女同欧美在线| 最新日韩AV网址在线观看| 国产美女无遮挡免费视频| 亚洲av日韩综合一区尤物| 性网站在线观看| 亚洲乱码在线视频| 视频二区中文无码| 欧美翘臀一区二区三区| 国产精品一区二区在线播放| 亚洲人成人无码www| 欧美精品啪啪一区二区三区| 日韩国产一区二区三区无码| 黄色网址手机国内免费在线观看| 国产综合精品一区二区| 亚洲一级色| 国产精品无码一区二区桃花视频| 国产区免费| 欧美精品v| 免费a在线观看播放| 秘书高跟黑色丝袜国产91在线| 亚洲日本中文字幕乱码中文| 亚洲精品片911| 国产97公开成人免费视频| 青青草国产一区二区三区| 毛片网站在线播放| 国产精品成人第一区| 国模私拍一区二区| 亚洲国产精品人久久电影| 91在线精品麻豆欧美在线| 免费在线色| 东京热高清无码精品| 日韩国产 在线| 色婷婷亚洲综合五月| 免费国产好深啊好涨好硬视频| 日本精品视频| 国产精品视频导航| 亚洲欧美在线综合一区二区三区| 最新国产你懂的在线网址| 久久婷婷五月综合97色| 亚洲人成网18禁| 欧美日韩一区二区在线播放| 欧美专区日韩专区| 国产亚洲精品97在线观看| a色毛片免费视频| 亚洲欧美极品| 久久精品91麻豆| 91久久偷偷做嫩草影院电| 午夜精品福利影院| 亚洲综合专区| 国产精品妖精视频| 日本一区二区不卡视频| 亚洲第一综合天堂另类专| 国产精品手机在线播放| 欧美成人精品高清在线下载| 97视频在线精品国自产拍| 欧美成人午夜视频免看| 国产高清不卡| 国产欧美高清| 久久人妻xunleige无码| 精品一区二区三区自慰喷水| 91色综合综合热五月激情| 在线免费看黄的网站| 91福利国产成人精品导航| 少妇人妻无码首页| 五月婷婷伊人网| 伊人久久婷婷| 夜夜爽免费视频| 欧美有码在线观看| 99中文字幕亚洲一区二区| 久久久久国色AV免费观看性色| 久久99国产精品成人欧美| 亚洲国产系列| 高清色本在线www| 国产色伊人| 992Tv视频国产精品|