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

多組學分析揭示拉格酵母應答高濃度麥汁機制

2021-05-21 09:16:46吳卓凡呼子暄王金晶劉春鳳鈕成拓鄭飛云李崎
食品與發酵工業 2021年9期
關鍵詞:差異分析

吳卓凡,呼子暄,王金晶,劉春鳳,鈕成拓,鄭飛云,李崎*

1(工業生物技術教育部重點實驗室(江南大學),江蘇 無錫,214122)2(江南大學 生物工程學院,江蘇 無錫,214122)

在全球變暖和化石燃料走向枯竭的背景下,節能減排一直是食品與發酵工業關注的主要問題。在啤酒生產中,高濃釀造作為一種有效的方式,可以減少水資源消耗、廢棄物排放、發酵罐數量和人工成本[1]。高濃釀造技術主要應用于拉格啤酒生產,拉格啤酒占全球產量的90%以上[2]。在高濃釀造過程中,酵母細胞暴露于高濃度麥汁的環境中,與常濃釀造相比,酵母的遲滯期延長、發酵周期延長并且殘余麥芽三糖顯著增多[3]。

自高濃釀造技術出現以來,許多研究者試圖從多方面提高酵母發酵性能:提升酵母抗氧化能力以改善發酵[4];篩選抗葡萄糖阻遏菌株來縮短發酵遲滯期[5];通過增強麥芽糖轉運來提高發酵末期的糖利用[6];通過添加外源氨基酸來提高酵母活性等[7]。提升酵母對高濃度麥汁的耐受是其中一個重要方向,高濃度麥汁是一個復雜的環境,高濃度的葡萄糖、麥芽糖與α-氨基氮都會影響酵母細胞的基因調控方向。然而目前對高濃度麥汁脅迫酵母細胞的機制沒有完善的研究。因此,揭示拉格酵母應對高濃麥汁的應答機制,可以為提升高濃釀造技術水平奠定基礎。隨著測序技術的發展和生物信息學分析方法的優化[8],轉錄組和代謝組的交叉分析可以很好地揭示不同條件下菌株的表達與代謝差異[9]。轉錄組學分析可以篩選不同條件下的顯著差異基因與代謝通路[10]。代謝組學分析可以驗證基因表達調整后的代謝物變化,并找到核心代謝物。

本研究通過轉錄組和代謝組分析了拉格酵母M14對常濃與高濃麥汁的應答機制。通過組學關聯分析篩選出13個具有顯著差異的重要基因。初步分析拉格酵母對高濃麥汁的應答機制,為后續高濃釀造優良酵母的選育提供方向。

1 材料與方法

1.1 菌株與培養基

工業拉格酵母(SaccharomycespastorianusM14)保藏于本實驗室,經過全基因組測序[11](NCBI ID:MVPU00000000)。

活化培養基:葡萄糖20 g/L、蛋白胨20 g/L、酵母粉10 g/L、瓊脂粉20 g/L(用于斜面培養基配制)。

麥汁培養基:加拿大麥芽3.2 kg、美國卡斯卡特啤酒花4 g、水9 L、瓊脂粉20 g/L(用于平板培養基配制),參照《釀造酒工藝學》[12]中的方法,制備成24 °P麥汁,其他濃度麥汁使用純凈水稀釋。

1.2 試劑與儀器

PBS緩沖液,Biosharp公司;標準品葡萄糖、果糖、麥芽糖,國藥試劑有限公司;RNA提取試劑盒(UNIQ-10 column Trizol total RNA extraction kit),生工生物工程(上海)股份有限公司;逆轉錄試劑盒(NEBNext?UltraTMRNA Library Prep Kit),Illumina公司。

SW-CJ-2D雙人單面凈化工作臺,蘇州凈化有限公司;Compass CX電子天平,Ohaus公司;UV11II紫外-可見分光光度計,上海天美公司;HYL-X恒溫搖床,強樂實驗設備公司;Research移液器,德國Eppendorf公司;Hve-50全自動高壓蒸汽滅菌鍋,日本平山制作所株式會社;R5408冷凍離心機,德國Eppendorf公司;Nanodrop2000核酸定量儀,美國Thermo Fisher公司;DMA4500啤酒分析儀,Anton-Paar公司。

1.3 濃度梯度麥汁平板脅迫測定

挑取1環M14細胞于YPD 培養基中,28 ℃,180 r/min,活化12 h。以體積分數10%的接種量進行翻接,28 ℃,180 r/min,培養12 h至穩定期(OD600 nm> 1.2)。7 000 r/min下4 ℃冷凍離心3 min,PBS緩沖液洗滌1次后使用無菌水梯度稀釋至10-5,取10 μL菌液在8、12、16、20、24 °P麥汁平板上點板,倒置于11 ℃培養箱靜置培養48 h。

1.4 發酵性能驗證及指標測定

使用12 °P麥汁分三級對拉格酵母M14進行擴培。第一級25 ℃,180 r/min,擴培18 h;第二級18 ℃,120 r/min,擴培24 h;第三級11 ℃,靜置,擴培24 h。取酵母細胞在7 000 r/min下4 ℃冷凍離心3 min,接種于300 mL 12 °P與24 °P麥汁培養基中,接種量為106個/(mL·Brix)。置于11 ℃培養箱發酵,每天記錄失重并取1 mL樣本分析糖譜。設置3個生物學重復。

發酵度、殘留浸出物與酒精度測定采用DMA4500啤酒分析儀,每個樣本重復測定2次。

糖譜分析使用高效液相色譜測定,參考ASBC分析方法[13],每個樣本測定1次。

1.5 RNA提取與轉錄組測序

使用1.4小節中方法制備種子液,取酵母細胞在7 000 r/min下4 ℃冷凍離心3 min后迅速接種于預冷至11 ℃的1 000 mL 12 °P與24 °P麥汁培養基中,接種量為106個/(mL·Brix)。11 ℃靜置2 h后取酵母細胞,7 000 r/min下4 ℃冷凍離心3 min,使用預冷的PBS緩沖液清洗3次,去除麥汁雜質后存儲細胞于1.5 mL EP管內。使用液氮速凍,存于-80 ℃冰箱備用。按照試劑盒說明書對酵母總RNA進行提取,提取后的RNA使用Nanodrop 2000進行核酸定量并進行瓊脂糖電泳檢測,檢測合格的RNA交于北京諾禾致源科技股份有限公司進行Illumina測序。設置3個生物學重復。

高通量測序儀測得的圖像數據經CASAVA堿基識別轉化為序列數據,然后對原始數據進行過濾,去除接頭與低質量數據得到clean data。同時,對clean data進行Q20、Q30和GC含量計算。后續所有分析均是基于clean data進行的高質量分析。使用HISAT2構建參考基因組的索引并將配對末端clean reads與M14基因組比對。featureCounts(1.5.0-p3)用于計算映射到每個基因的讀數。然后根據基因的長度計算每個基因的FPKM,并計算映射到該基因的讀數。使用DESeq2軟件進行2個比較組合之間的差異表達分析[14](每個組3個生物學重復),設置校正后的P<0.05以及|log2fold change|≥1作為顯著差異表達的閾值。通過clusterProfiler(3.4.4)軟件實現差異表達基因的GO(gene ontology)富集分析與KEGG(Kyoto encyclopedia of genes and genomes)通路中差異表達基因的統計富集,其中修正了基因長度偏差。

1.6 胞內代謝物提取與代謝組檢測

代謝組學樣本與轉錄組學樣本取樣方法相同,采用正離子與負離子模式分別檢測[15],在本文的結果與分析中合并分析,設置6個生物學重復。

取100 mg液氮研磨的組織樣本,置于EP管中,加入500 μL含1 L/m3甲酸的800 L/m3甲醇水溶液,渦旋振蕩,冰浴5 min,4 ℃、15 000 r/min離心10 min,取一定量的上清液使用超純水稀釋至甲醇含量為600 L/m3,并置于帶有0.22 μm濾膜的離心管中4 ℃、15 000 r/min離心10 min,收集濾液,進樣LC-MS進行分析。從每個實驗樣本中取等體積樣本混勻作為QC樣本??瞻讟颖緸楹? L/m3甲酸的600 L/m3甲醇水溶液代替實驗樣本,前處理過程與實驗樣本相同。

表1 色譜梯度洗脫程序Table 1 Chromatographic gradient elution program

超高效液相色譜使用Hyperil Gold column(C18)色譜柱,色譜條件:柱溫40 ℃、流速0.2 mL/min;正模式:流動相A,1 L/m3甲酸、流動相B,甲醇;負模式:流動相A,5 mmol/L醋酸銨,pH 9.0、流動相B,甲醇。

質譜掃描范圍選擇m/z100~1 500;ESI源的設置如下:Spray Voltage:3.2 kV;Sheath gas flow rate:35 arb;Aux Gasflow rate:10 arb;Capillary Temp:320 ℃。MS/MS二級掃描為data-dependent scans。

將下機數據文件導入Compound Discoverer 3.0軟件中,進行保留時間、質荷比等參數的簡單篩選,然后對不同樣品根據保留時間偏差0.2 min和質量偏差5 mg/L進行峰對齊,使鑒定更準確,隨后根據設置的質量偏差5 mg/L、信號強度偏差30%、信噪比3、最小信號強度100 000、加和離子等信息進行峰提取,同時對峰面積進行定量,再整合目標離子,然后通過分子離子峰和碎片離子進行分子式的預測并與mzCloud和Chemspider數據庫進行比對,用blank樣本去除背景離子,并對定量結果進行歸一化,最后得到數據的鑒定和定量結果。采用PLS-DA模型第一主成分的變量投影重要度(variable importance in the projection,VIP)值,VIP值表示不同分組中代謝物差異的貢獻率;并結合T-test的P值來尋找差異性表達代謝物,設置VIP>1.0、校正后的P<0.05以及|log2foldchange|≥1作為顯著差異代謝物的閾值。

1.7 數據分析

使用SPSS 20.0軟件進行統計學顯著性分析。P<0.05視為顯著差異。使用GraphPad Prism 8進行圖片繪制。

2 結果與分析

2.1 濃度梯度麥汁平板脅迫測定及發酵實驗

為了研究高濃度麥汁對拉格酵母的脅迫水平,通過濃度梯度麥汁平板進行了簡單的酵母生長情況與耐受性測定。如圖1所示,拉格酵母M14在8 °P與12 °P麥汁平板上基本沒有生長抑制,在16 °P麥汁平板上10-5濃度的細胞生長量顯著降低,而在20 °P與24 °P麥汁平板上10-5濃度的酵母基本無法生長。表明高濃度麥汁給酵母細胞生長帶來了較大的壓力。

圖1 M14菌株的高濃度麥汁脅迫平板實驗Fig.1 Plate experiment of strain M14 under high gravity wort stress

通過小瓶發酵實驗比較常濃與高濃釀造發酵指標,并跟蹤發酵過程的失重與糖降,有助于了解常濃度與高濃度麥汁造成的發酵性能與糖利用差異。由圖2可知,M14的發酵度由12 °P中的58.56%下降到24 °P中的56.01%(圖2-a);12 °P麥汁中酵母起酵較快,同時發酵5 d主酵接近結束,麥芽糖被完全利用(圖2-c);24 °P麥汁中酵母起酵較慢,同時7 d未能完全消耗麥芽糖(圖2-d);實驗使用的麥汁糖譜符合ASBC中相關比例標準[13](圖2-b)。

a-12 °P與24 °P麥汁的發酵度;b-麥汁的糖譜;c-12 °P下葡萄糖、果糖、麥芽糖與失重變化;d-24 °P下葡萄糖、果糖、麥芽糖與失重變化圖2 M14菌株的常濃與高濃發酵比較Fig.2 Comparison of normal and high gravity fermentation of strain M14

綜上所述,高濃度麥汁對拉格酵母有著較強脅迫,酵母在高濃釀造下遲滯期延長、生長速度下降并且糖分利用不完全,導致發酵周期延長和發酵度下降。因此本研究采用轉錄組與代謝組手段,進一步分析發酵早期代謝物與表達情況,從分子水平解釋高濃度麥汁脅迫下的拉格酵母應答機制。

2.2 轉錄組分析

2.2.1 樣本設置、主成分分析與基因表達熱圖

起酵階段的應答對拉格酵母發酵性能至關重要。為了探究發酵早期拉格酵母對高濃麥汁的應答機制,轉錄組與代謝組樣本的取樣時機設置在起酵階段。接種前的酵母參照工業釀造工藝流程進行逐級降溫、擴培,確保其接入高濃與常濃發酵培養基時的狀態一致、活力良好,并于接種后2 h取樣。

轉錄組樣本上機測序后,經過原始數據過濾、測序錯誤率檢查、GC含量分布檢查,獲得后續分析使用的clean data。使用經過注釋的自身參考基因組進行比對,共注釋到6 598個基因,然后進行定量分析。采用主成分分析評估組間差異及組內樣本重復情況,主成分分析采用線性代數的計算方法,對所有基因變量進行降維處理及主成分提取。對6個樣本的基因表達值(FPKM)進行主成分分析(圖3-b),組間樣本分散,組內樣本聚集,表明兩組間具有顯著差異且測序情況良好,說明拉格酵母M14對高濃麥汁產生了顯著的應答調控。

同時繪制了表達基因聚類熱圖(圖3-a)。通過主流的層次聚類對基因的FPKM值進行聚類分析,熱圖中的基因越靠近其表達模式越相近。因為經過歸一化處理,熱圖中的顏色只能橫向比較(同一基因在不同樣本中的表達情況),不能縱向比較(同一樣本不同基因的表達情況)。由圖可知M14在應答12 °P與24 °P麥汁的情況下,組內的表達模式完全一致,組間的總體表達模式接近,但是存在局部的差異(圖3-a中紅框),有上調的差異區域也有下調的差異區域。這些顯著差異的區域是拉格酵母應答高濃麥汁處理的重要基因與途徑。

a-表達熱圖;b-主成分分析;c-差異基因火山圖圖3 M14菌株應答常濃與高濃麥汁的表達情況Fig.3 Gene expression of strain M14 response to normal and high gravity wort

通過宏觀分析,驗證了實驗設計正確和平行性良好,且說明了常濃度與高濃度麥汁處理的拉格酵母差異顯著,需要進一步進行顯著基因差異分析與通路富集分析。

2.2.2 差異表達基因分析

經過測序結果質檢與宏觀分析,樣本平行性良好,故在差異表達基因與顯著差異表達基因的篩選中設置padj閾值為0.05與0.01。如表2所示,經過2組樣本間差異表達基因的統計與篩選,一共獲得了2 312個差異表達基因,占總基因數的35.04%,其中1 137個上調,1 175個下調,下調基因略多于上調基因(圖3-c)。設置|log2fold change| ≥ 1作為顯著差異基因篩選的閾值,得到191個顯著差異基因,占總差異表達基因的8.26%,其中38個基因顯著上調,153個基因顯著下調,下調基因多于上調基因。顯著差異基因中最大上調倍數為7.15倍,最大下調倍數為11.53倍。

表2 轉錄組基因統計Table 2 Transcriptomic gene statistics

大量基因顯著下調可能是高濃度麥汁的脅迫造成的,但也有一定量的基因上調,這可能是拉格酵母M14為了應對麥汁脅迫做出的應答。通過顯著差異基因分析,目標基因被鎖定到191個以內,后續將通過GO、KEGG富集分析與代謝組數據支撐找到重要差異基因。

2.2.3 差異表達基因的GO富集分析

GO是描述基因功能的綜合性數據庫,基因被分為生物過程(biological process,BP)、細胞組成(cellular component,CC)和分子功能(molecular function,MF)3個部分。在3個部分中各選擇10個以內富集到的顯著差異通路進行柱狀圖繪制,各部分從左往右顯著性從高到低如圖4所示。高濃度與常濃度的差異基因主要富集到生物過程(BP)的碳水化合物代謝過程(GO:0005975)與氧化還原過程(GO:0055114);富集到細胞組成(CC)的膜的固有成分(GO:0031224)、膜的組成部分(GO:0016021)、蛋白酶體核心復合物(GO:0005839)與蛋白酶體復合體(GO:0000502);富集到細分子功能(MF)的作用于供體的氧化還原酶活性(GO:0016614)、維生素結合(GO:0019842)、作用于受體氧化還原酶活性 (GO:0016616)與輔因子結合(GO:0048037)。

圖4 差異基因的GO富集Fig.4 GO enrichment of DEGs

3個部分富集到的基因功能高度相關,可以將富集到的基因功能概括為碳代謝、氧化還原平衡、膜組成與蛋白酶體4個類別。其中碳代謝與氧化還原平衡與高濃脅迫直接相關,高濃麥汁中的高糖脅迫會直接影響能量代謝;細胞膜和酵母抗逆一直息息相關[16],膜功能的狀態決定了酵母生理活性的水平;蛋白酶體功能涉及細胞周期控制、應激反應、基因轉錄與信號轉導多方面[17],可能參與高濃脅迫下的細胞宏觀調控。上述結果能夠進一步縮小基因范圍,加快基因挖掘進程。

2.2.4 差異表達基因的KEGG富集分析

KEGG綜合性數據庫整合了基因組、化學和系統功能信息。選取最顯著的18個通路繪制散點圖(圖5),點的大小代表了富集到的基因數量,展示在下方的通路顯著性較高。

圖5 差異基因的KEGG富集Fig.5 KEGG enrichment of DEGs

具有顯著差異且注釋到基因數量較多的通路為:碳代謝、氨基酸代謝與蛋白酶體。其中氨基酸代謝普遍上調(40/66),可能由于脅迫刺激細胞內蛋白合成。碳代謝普遍下調(43/69),可能由于高濃脅迫壓力,導致糖轉運與EMP途徑無法快速響應,相比常濃度麥汁更晚進入穩定耗糖階段,也可能是由于細胞資源被用于應對脅迫壓力,影響能量代謝表達水平。KEGG富集與GO富集分析結果基本一致。

2.3 代謝組分析

2.3.1 差異代謝物分析

本研究采用了與轉錄組一致的樣本進行代謝組檢測,設置6個生物學平行。下機數據經過圖譜處理與數據庫搜索,對代謝物進行定性與定量,總共分辨到326種代謝物。質控數據顯示代謝組數據可靠、準確。主成分分析顯示2組間樣本分散且組內樣本聚集,表明2組間存在顯著差異代謝物且代謝物含量穩定。

設置顯著差異代謝物篩選的閾值為VIP>1.0,|log2fold change|≥1且P<0.05,最終篩選到30個顯著差異代謝物,占總代謝物數量的9.20%,其中18個代謝物顯著上調,12個代謝物顯著下調。顯著差異代謝物中最大上調倍數為2.75倍,最大下調倍數為2.31倍。代謝物調整具有滯后性,作為對轉錄組數據的支撐,篩選到的顯著代謝物可以論證相關基因通路的表達差異,但是未能表現出顯著差異的代謝物不能說明相關基因通路無意義。上述數據表明酵母在短時間內便對高濃麥汁應答,基因表達差異引起了部分代謝物含量變化。

2.3.2 差異代謝物的KEGG富集分析

通過相同的方法對差異代謝物進行KEGG富集分析,選取最顯著的18個通路繪制散點圖(圖6)。其中有10個通路與氨基酸代謝或能量代謝相關,其中最顯著的5個通路內有3個為氨基酸代謝。如表3所示,天冬氨酸、絲氨酸、精氨酸與蛋氨酸含量顯著上調。說明胞內氨基酸和氨?;衔镌诟邼舛塞溨碳は潞铣杉铀?與轉錄組結果一致。而細胞內糖類化合物還沒有顯著變化,說明了碳代謝普遍下調不是因為胞內糖分過量酵解使得ATP積累導致的反饋抑制,而可能與麥汁脅迫壓力相關。

2.4 重要差異基因篩選與分析

轉錄組分析篩選到顯著差異基因并富集到能量代謝、氨基酸代謝與蛋白酶體,代謝組分析驗證了氨基酸代謝途徑的基因上調,并揭示了碳代謝途徑下調的原因。對代謝物和差異基因進行關聯分析,并結合酵母SGD數據庫(saccharomyces genome database)與

圖6 差異代謝物的KEGG富集Fig.6 KEGG enrichment of differential metabolites

表3 顯著差異氨基酸Table 3 Amino acids with significantly different content

現有研究文獻[18-20],篩選獲得了13個重要差異基因。如表4所示,共確定了13個重要差異基因,其中8個基因與碳代謝相關,2個基因與己糖轉運相關,3個基因與氨基酸的生物合成相關。

表4 重要差異基因Table 4 Important DEGs

代謝組中顯著上調的氨基酸與轉錄組中顯著上調的氨基酸相關基因(AAT1、SER2與ARG4)完全一致,除顯著差異基因之外,氨基酸生物合成途徑普遍上調,表明拉格酵母在應答高濃度麥汁時胞內蛋白質與酶合成活躍。氨基酸的積累對細胞抗逆性能有著積極作用[21],比如帶電氨基酸如精氨酸對細胞耐低溫有著明顯效果[22]。氨基酸的加速合成與積累是拉格酵母細胞適應外界環境的應答策略,有助于細胞維持生長代謝平衡。

高濃度麥汁對拉格酵母細胞最明顯的壓力是高濃度的葡萄糖與麥芽糖。EMP途徑的基因全線下調,包括關鍵的HXK1與PGK1,同一代謝通路上下游的PGM2、GPM1與GPM2也顯著下調。EMP途徑的顯著下調可能是拉格酵母在高濃麥汁中起酵緩慢的主要原因。受到高濃度麥汁脅迫,己糖激酶HXK1與磷酸甘油酸激酶PGK1等關鍵基因表達下調,糖酵解速度放緩導致細胞內ATP生產減少,不利于細胞其他代謝活動開展。同時細胞內海藻糖合成相關基因顯著上調(TSL1、TPS3與TPS2),海藻糖積累對細胞代謝平衡有重要作用[23]。

同時在重要通路之外觀測到了HXT家族基因的顯著變化(HXT1與HXT2),HXT家族編碼己糖轉運蛋白。已有研究表明HXT1的表達和發酵初期的極端高糖相關[24],所以HXT1可能是拉格酵母應答高濃葡萄糖的主要基因。HXT2同時具有低濃與高濃的底物特異性親和力[25],與HXT1的誘導系統不同,HXT2可能受到環境因素如滲透壓、ATP缺乏等調控。

3 結論

本研究通過濃度梯度麥汁平板脅迫測定及發酵實驗,比較了常濃度與高濃度下的酵母特性。通過對拉格酵母M14在12 °P與24 °P麥汁處理2 h后的樣本進行轉錄組測序、代謝組檢測與關聯分析,共獲得191個顯著差異基因,其中38個上調,153個下調;共獲得30個顯著差異代謝物,其中18個上調,12個下調。通過GO富集與KEGG富集分析,結合SGD數據庫,篩選獲得13個重要差異基因(HXT1、HXT2、HXK1、PGM2、PGK1、GPM1、GPM2、TSL1、TPS3、TPS2、AAT1、SER2與ARG4),差異基因主要與糖酵解、海藻糖合成、氨基酸合成和己糖轉運相關。本研究初步揭示了拉格酵母應對高濃度麥汁的應答機制,為高濃釀造優良酵母的選育提供技術參考。

猜你喜歡
差異分析
相似與差異
音樂探索(2022年2期)2022-05-30 21:01:37
隱蔽失效適航要求符合性驗證分析
找句子差異
DL/T 868—2014與NB/T 47014—2011主要差異比較與分析
電力系統不平衡分析
電子制作(2018年18期)2018-11-14 01:48:24
生物為什么會有差異?
電力系統及其自動化發展趨勢分析
M1型、M2型巨噬細胞及腫瘤相關巨噬細胞中miR-146a表達的差異
中西醫結合治療抑郁癥100例分析
收入性別歧視的職位差異
主站蜘蛛池模板: 日韩欧美在线观看| 国产综合在线观看视频| 99这里只有精品免费视频| 国产一国产一有一级毛片视频| 国产一级毛片网站| 国产在线一区二区视频| 成年网址网站在线观看| 97超碰精品成人国产| 91视频青青草| 91破解版在线亚洲| 精品超清无码视频在线观看| 欧美性天天| 91精品啪在线观看国产| 国产精品尤物铁牛tv | 久久婷婷国产综合尤物精品| 好吊妞欧美视频免费| AV天堂资源福利在线观看| 亚洲大学生视频在线播放| 成年av福利永久免费观看| 九色91在线视频| 高h视频在线| 国产精品不卡永久免费| 欧美一级色视频| 不卡午夜视频| 四虎国产精品永久一区| 亚洲AⅤ无码国产精品| 四虎永久在线精品国产免费 | 欧洲欧美人成免费全部视频| 就去吻亚洲精品国产欧美| 亚洲二区视频| 91精品国产丝袜| 欧美日韩成人在线观看| 欧美亚洲香蕉| 国产精品久久国产精麻豆99网站| 久久鸭综合久久国产| 亚洲欧美不卡中文字幕| 国产玖玖玖精品视频| 亚洲熟妇AV日韩熟妇在线| 亚洲天堂精品视频| 欧美亚洲一区二区三区导航| 精品無碼一區在線觀看 | 日韩欧美一区在线观看| 国产一区三区二区中文在线| 呦视频在线一区二区三区| 久久五月视频| 亚洲91在线精品| 日本午夜网站| 亚洲国产成人超福利久久精品| 久久精品丝袜高跟鞋| 亚洲三级成人| 亚洲精品无码不卡在线播放| 日韩无码黄色| 又爽又大又光又色的午夜视频| 国产婬乱a一级毛片多女| 青青草国产免费国产| 亚洲视频无码| 狠狠v日韩v欧美v| 久久精品国产精品一区二区| 国产女人在线视频| 99热国产这里只有精品无卡顿" | 国产打屁股免费区网站| 国产精品无码作爱| 亚洲欧洲美色一区二区三区| 亚洲另类色| 天天爽免费视频| 狠狠做深爱婷婷综合一区| 97在线公开视频| 99九九成人免费视频精品| 98超碰在线观看| 日韩不卡免费视频| 中文字幕伦视频| 91无码人妻精品一区| 一级毛片免费观看久| 亚洲欧美精品日韩欧美| 国产亚洲美日韩AV中文字幕无码成人| 欧美一级高清视频在线播放| 波多野结衣一区二区三区88| 岛国精品一区免费视频在线观看| 暴力调教一区二区三区| 波多野结衣第一页| 色妞永久免费视频| 欧美综合一区二区三区|