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

基于轉錄組挖掘球藥隔重樓螺甾烷型重樓皂苷生物合成相關基因

2023-10-19 10:51:36曲別阿香曲別軍長俸明康何曉勇馮景秋盛華春黃艷菲閻新佳張紹山
中草藥 2023年20期
關鍵詞:途徑生物

曲別阿香,曲別軍長,俸明康,何 斌,何曉勇,李 娟,馮景秋,盛華春,黃艷菲,閻新佳,劉 圓*,張紹山*

? 藥材與資源?

基于轉錄組挖掘球藥隔重樓螺甾烷型重樓皂苷生物合成相關基因

曲別阿香1, 2, 3,曲別軍長1, 2, 3,俸明康1, 2, 3,何 斌1, 2, 3,何曉勇1, 2, 3,李 娟2, 3,馮景秋2, 3,盛華春2, 3,黃艷菲1, 3, 4,閻新佳2, 3,劉 圓2, 3*,張紹山2, 3*

1. 西南民族大學藥學院,四川 成都 610041 2. 四川省羌彝藥用資源保護與利用技術工程實驗室,四川 成都 610225 3. 青藏高原民族藥用資源保護與利用國家民委重點實驗室,四川 成都 610225

以球藥隔重樓根莖、莖和葉為材料,探究螺甾烷型重樓皂苷的生物合成基因。采用HPLC分析6個螺甾烷型重樓皂苷含量,進一步利用Illumina HiSeq X Ten平臺進行轉錄組測序,結合定量和轉錄組數(shù)據進行生信分析。6個螺甾烷型重樓皂苷在球藥隔重樓根莖、莖和葉中的含量具有顯著差異。Illumina轉錄組測序共獲得61 755條非冗余unigenes,其中30 263條(49.0%)被成功注釋。催化生成薯蕷皂苷元的31個關鍵基因均被成功鑒定,且基因表達模式與皂苷的積累水平總體一致,其中、、、、、、、、、、、、、、和等基因與重樓皂苷的積累呈顯著正相關。共獲得了61條尿苷二磷酸糖基轉移酶()基因,其中30條潛在參與重樓皂苷的糖基化修飾。球藥隔重樓根莖、莖和葉中皂苷合成基因表達模式與皂苷積累規(guī)律存在一致性,鑒定的基因可指導該類活性產物的生物合成途徑解析。

球藥隔重樓;轉錄組;螺甾烷型重樓皂苷;生物合成途徑;尿苷二磷酸糖基轉移酶

螺甾烷型重樓皂苷是著名中藥材“重樓”的主要活性成分,具有顯著的抗腫瘤、抗菌、止血等藥理活性,該類皂苷也是歷版《中國藥典》重樓項下的質控成分[1-3]。據螺甾烷型皂苷元C-17位上有無-OH,又分為薯蕷皂苷和偏諾皂苷,而偏諾皂苷又可在其C-23、C-24和C-27位加上-OH;除苷元母核的差異外,各皂苷結構的差異主要來自于苷元母核上連接的寡糖基[3]。進一步的藥理研究表明,螺甾烷型重樓皂苷元的結構修飾,對其抗腫瘤、抗菌、止血等藥理活性有顯著影響,如糖基化后形成的重樓皂苷I、II、Ⅲ、VII、H,抗菌活性顯著增強,其中重樓皂苷H的止血活性亦顯著增強[3-5]。遺憾的是,自然界中,螺甾烷型重樓皂苷在重樓屬植物中含量極低,外加重樓屬植物生長緩慢,限制了該類化合物的開發(fā)利用[6-7]。因此,研究該類化合物在重樓屬植物中的合成代謝,對該藥材的藥效形成機制和品質提升等具有重要意義。

螺甾烷型重樓皂苷的生物合成途徑復雜。首先通過細胞質中的甲羥戊酸途徑(mevalonate)和質體中的甲基赤蘚糖醇-4-磷酸途徑(methylerythritol 4-phosphate)合成法尼基焦磷酸(farnesyl pyrophosphate,F(xiàn)PP)[8-9]。FPP在角鯊烯合成酶(SQS)[10]、角鯊烯單氧化酶(SQE)[11]、環(huán)阿屯醇合成酶(CAS)[12]、甾醇-4α-甲基氧化酶(SMO1-3和SMO2-2)[13-14]、環(huán)桉烯醇環(huán)異構酶(CPI-5)[14-15]、甾醇-14α-去甲基酶(CYP51G)[14-16]、甾醇-C14-還原酶(C14-R-2)[14,17]、膽甾烯醇δ-異構酶(8,7SI-4)[14,18]、甾醇C-5(6)去飽和酶(C5-SD1)[14,19]、7-脫氫膽甾烯醇還原酶(7-DR1)[14,20]和δ(24)-甾醇還原酶(SSR1-3)[14,21]的催化下,形成膽甾烯醇。目前,從重樓中鑒定了3種不同的基因組合(CYP90G4/ CYP94D108、CYP90G4/CYP94D109、CYP90G4/ CYP94D72A616),可催化膽甾烯醇形成螺甾烷型皂苷元[22]。皂苷元最終在糖基轉移酶(Uridine diphosphate glycosyltransferases,UGTs)的糖基化下,主要在C3-OH上連接不同種類(主要為葡萄糖基、鼠李糖基、阿拉伯糖基)、不同數(shù)量(主要為單糖、二糖、三糖和四糖)和不同位置(主要為C-2、C-3和C-4位的-OH)的糖基,從而形成不同的螺甾烷型重樓皂苷,糖基化步驟中報道了葡糖基化皂苷元C3-OH的UGT73CR1、UGT80A40和UGT80A41,2′--鼠李糖基化的UGT73CE1,以及6′--葡萄糖基化的UGT91AH1和UGT91AH2[12,23],而其他糖基化步驟的UGTs尚無研究報道。

球藥隔重樓Franch.為重樓屬L.多年生草本植物,為川渝地區(qū)珍稀品種,相較重樓屬其他物種,其同時富含薯蕷型和偏諾型重樓皂苷,是探究螺甾烷型重樓皂苷生物合成相關基因的理想材料[24-25]。目前,關于螺甾烷型重樓皂苷的分子水平研究主要基于資源豐富的Smith(Franch.) Hand.Mazz.或七葉一枝花Smithvar.(Franch.) Hara[12,26]2022年,Hou等[27]利用Illumina NovaSeq 6000平臺對七葉一枝花、云南重樓、球藥隔重樓、狹葉重樓var.Franch.和毛重樓H. Léveillé的根莖和葉分別做了轉錄組測序,以云南重樓為參照,在葉中檢測到了14 815個球藥隔重樓特異性表達的基因,是其他3組中最多的,且在根莖中也檢測到了6030個特異性表達基因,這說明球藥隔重樓相比其他重樓品種,有著自身獨特的遺傳背景。

因此,本研究利用前期繁育的7年生球藥隔重樓,基于RNA-Seq技術對其根莖、莖、葉進行轉錄組分析,進一步結合根莖、莖、葉中螺甾烷型重樓皂苷的含量數(shù)據,鑒定該類皂苷生物合成途徑基因,為該類皂苷生物合成途徑解析和優(yōu)良重樓品種的分子育種提供科學支撐。

1 材料

1.1 植物和取樣

通常重樓皂苷在重樓中的積累于7~8年時豐度最高[28],故本研究選擇7年生無病害的球藥隔重樓,于2022年8月采自四川省阿壩藏族羌族自治州汶川縣水磨鎮(zhèn)西南民族大學藥材繁育試驗基地,由西南民族大學劉圓教授和張紹山博士鑒定為球藥隔重樓Franch.。將植株均分為3份,每份4株,取每份植株的根莖、莖、葉,按組織部位混勻后,一半40 ℃烘干,粉碎,過40目篩,干燥器保存,用于螺甾烷型重樓皂苷含量分析;一半裝于潔凈的離心管中,迅速用液氮固定,?80 ℃保存,用于轉錄組分析。

1.2 試藥

對照品重樓皂苷I(批號RP200523)、重樓皂苷II(批號RP190315)、重樓皂苷H(批號RP190528)、重樓皂苷V(批號RP191217)、重樓皂苷VI(批號RP220807)購自成都麥德生科技有限公司,重樓皂苷VII(批號MUST-18062201)購自成都曼斯特生物有限公司。色譜級甲醇和乙腈均購于賽默非世爾科技有限公司;水為屈臣氏蒸餾水;其余試劑為分析純。

2 方法

2.1 螺甾烷型重樓皂苷含量測定

2.1.1 色譜條件 Agilent 1260 Infinity II色譜儀;Diamonsil-C18色譜柱(250 mm×4.6 mm,5 μm);流動相乙腈(A)-水(B),梯度洗脫:0~20 min,30%~40% A;20~45 min,40%~60% A;體積流量1 mL/min;柱溫30 ℃;進樣體積10 μL;檢測波長203 nm。

2.1.2 對照品溶液的制備 分別精密稱取重樓皂苷I、II、V、VI、VII和H對照品適量,加色譜甲醇制成含重樓皂苷I 0.478 mg/mL、重樓苷II 0.285 mg/mL、重樓皂苷V 0.500 mg/mL、重樓皂苷VI 0.514 mg/mL、重樓皂苷VII 0.240 mg/mL和重樓皂苷H 0.375 mg/mL的混合對照品溶液,搖勻后經0.45 μm濾膜濾過。

2.1.3 供試品溶液的制備 取不同部位球藥隔重樓樣品粉末約0.5 g,精密稱定,按料液比1∶25精密加入甲醇,密封,稱定質量,超聲60 min(功率500 W,頻率40 kHz),放冷至室溫,再稱定質量,加甲醇補足減失的質量,搖勻后經0.45 μm濾膜濾過。

2.1.4 線性方程的繪制 取對照品適量,配制成重樓皂苷I 、重樓皂苷II、重樓樓甙V、重樓皂苷VI、重樓皂苷VII、重樓皂苷H的對照品儲備液,使其質量濃度分別為0.478、0.285、0.500、0.514、0.240、0.375 mg/mL,取混合對照品溶液 0.1、0.5、1.0、1.5、2.0、2.5mL于5 mL量瓶中加甲醇定容至刻度,按上述梯度洗脫條件,分別精密吸取供試品溶液各10 μL,注入液相色譜儀,測定,即得。以6種甾體皂苷的進樣濃度為橫坐標(),相應峰面積為縱坐標(),繪制標準曲線,計算回歸方程,對應的對照品重樓皂苷I、重樓皂苷II、重樓皂苷H、重樓皂苷V、重樓皂苷VI、重樓皂苷VII的線性方程分別為=3 333.6+72.454、=3 003.1+3.740 4、=3 051.7+5.760 4、=3 500.8+169.78、=3 575.3+26.1、=10 330+31.919;對應的線性范圍分別為0.015~0.478、0.009 5~0.285 0、0.012~0.375、0.016~0.500、0.017~0.514、0.008~0.240 mg/mL;2分別為0.999 9、0.999 9、1.000 0、0.994 3、0.999 9、0.999 9。結果表明各成分在相應范圍內呈良好的線性關系。

2.1.5 精密度試驗 精密稱取“2.1.2”項下混合對照品一份,按“2.1.1”項下色譜條件連續(xù)進樣6針,測得各指標成分色譜峰的峰面積,結果表明樣品中6種指標成分峰面積RSD值在0.98%~1.70%。

2.1.6 重復性試驗 精密稱取球藥隔重樓樣品粉末,按“2.1.3”項下方法制備6份供試品溶液,按“2.1.1”項下色譜條件進樣10 μL,測得各指標成分色譜峰的峰面積,結果表明樣品中6種指標成分質量分數(shù)的RSD值在1.01%~1.92%。

2.1.7 穩(wěn)定性試驗 精密吸取同一份球藥隔重樓樣品粉末供試品溶液,按“2.1.1”項下色譜條件分別于制樣后0、2、6、8、12、24 h進樣,測得各指標成分色譜峰的峰面積,結果表明樣品中6種指標成分峰面積的RSD值在1.05%~1.74 %,表明供試品溶液在制備后24 h內穩(wěn)定。

2.1.8 加樣回收率試驗 精密稱取球藥隔重樓樣品粉末 0.25 g,精密加入等量的對照品,按“2.1.3”項下供試品溶液制備方法平行制備6份樣品溶液,測定6個指標性成分的含量,計算加樣回收率,加樣回收率均在95.72%~110.20%,RSD為0.99%~1.82%,符合定量分析要求。

2.1.9 樣品含量測定 將“2.1.3”項下制得的9批(每個組織部位3個生物學重復)供試品溶液,按“2.1.1”項下色譜條件測定。

2.2 轉錄組測序和分析

2.2.1 文庫構建和測序 利用mirVana miRNA Isolation Kit(Cat. AM1561,InVItrogen,Thermo Fisher Scientific Inc.,美國)分別提取根莖、莖和葉的總RNA。RNA用瓊脂糖凝膠電泳檢測降解和污染情況,NanoDrop 2000 spectrophotometer(Thermo Fisher Scientific,Waltham,MA,美國)檢測RNA的純度和濃度,Agilent 2100 Bioanalyzer(Agilent Technologies,Santa Clara,CA,美國)評估RNA的完整性。符合高通量測序質量要求后,委托上海歐易生物醫(yī)學科技有限公司采用Illumina HiSeq X Ten platform(Illumina,美國)進行測序,文庫的構建采用TruSeq Stranded mRNA LTSample Prep Kit(Illumina,美國)試劑盒。

2.2.2數(shù)據質量控制與組裝 首先使用Trimmomatic軟件進行質控并去除接頭,在此基礎上過濾掉低質量堿基以及N堿基,最終得到高質量的clean reads。使用Trinity(version:2.4)軟件paired-end的拼接方法,將clean reads拼接得到transcript序列,根據序列相似性及長度,挑選出最長的1條作為unigene。最終,利用CD-HIT(identify=98%)軟件聚類去冗余得到非冗余的unigenes。

2.2.3 Unigenes功能注釋和定量分析 利用diamond軟件將unigenes的蛋白序列比對到NR、KOG、GO、Swiss-Prot、eggNOG、KEGG數(shù)據庫,以及利用HMMER軟件比對到Pfam數(shù)據庫,對unigenes進行功能注釋分析。使用軟件bowtie2獲取每個樣本中比對到unigenes上的reads數(shù),同時使用 eXpress軟件來計算unigenes在各樣本中的基因表達量(fragments per kilobase of exon model per million mapped fragments,F(xiàn)PKM)值。

2.2.4 差異基因分析 利用DESeq2 R package軟件包篩選各樣本間的差異基因(differentially expressed genes,),篩選差異的條件為<0.05且差異倍數(shù)(foldchange)>2,對差異表達基因做KEGG通路顯著性分析。

2.3 螺甾烷型重樓皂苷生物合成相關基因的篩選和分析

以unigene序列集為數(shù)據庫,利用已發(fā)表三萜甾體皂苷生物合成途徑基因序列,進行blast檢索,鑒定球藥隔重樓中的同源基因,采用pearson相關性算法,計算上述同源基因與螺甾烷型重樓皂苷積累水平的關系。對尚未解析皂苷母核糖基修飾的,首先依據UGT的隱馬爾可夫模型(HMM)文件(PF00201)進行基因家族鑒定和功能分析,結合上游已鑒定基因的表達模式和重樓皂苷定量分析結果,篩選參與皂苷母核糖基修飾的候選。運用TBtools繪制基因不同組織相對表達量熱圖;運用MEGA7軟件,選取neighbor-joining方法,bootstrap為1000個重復,構建UGTs的系統(tǒng)發(fā)育進化樹。

2.4 qRT-PCR

利用RNA提取試劑盒“FastPure? Universal Plant Total RNA Isolation kit(Vazyme,China)”從9份RNA-seq材料中提取RNA,質量經檢測符合要求后,利用HiScript?II Q RT SuperMix for qPCR (+gDNA wiper) KitRNA反轉錄試劑盒獲取cDNA,Primer6.0軟件設計引物,RT-PCR檢查引物特異性。qRT-PCR實驗采用2×Taq Pro Universal SYBR qPCR Master Mix kit(Vazyme)試劑盒,具體操作按說明書進行。qRT-PCR 體系為10.0 μL:2×Taq Pro Universal SYBR qPCR Master Mix 5.0 μL,上游引物 0.20 μL,下游引物 0.20 μL,cDNA 1.0 μL,ddH2O 3.6 μL,配制過程在冰上完成。PCR 擴增程序為95 ℃、30 s;95 ℃、5 s,53 ℃、30 s,共40個循環(huán);95 ℃、15 s,60 ℃、60 s,95 ℃、15 s。以為內參基因,每個樣品設置3個生物學重復和3個技術重復。

3 結果

3.1 不同組織部位6個螺甾烷型重樓皂苷含量測定

經對建立的高效液相色譜(High-performance liquid chromatography,HPLC)方法學考察,結果表明各方法學均符合要求,方法可用于相關重樓皂苷的定量分析。含量測定結果表明,球藥隔重樓中含有薯蕷皂苷(重樓皂苷I、II和V)和偏諾皂苷(重樓皂苷H、VI和VII);重樓皂苷I在莖中沒有檢測到;重樓皂苷II和V僅在根莖中檢測到;重樓皂苷H、VI和VII在3個組織部位均檢測到,重樓皂苷H在根莖中的質量分數(shù)超過了30.0 mg/g,除重樓皂苷H在莖和葉中差異不顯著外,其余皂苷在各組織部位間的含量差異顯著(<0.05);此外,除重樓皂苷VII在葉中含量最高外,其余皂苷均在根莖中積累水平最豐富(圖1)。以上結果表明球藥隔重樓的根莖、莖和葉是揭示螺甾烷型重樓皂苷生物合成通路的理想材料。

3.2 轉錄組測序和de novo組裝

球藥隔重樓根莖、莖和葉的9個cDNA文庫經Illumina HiSeq X Ten平臺測序,獲得了66.78 Gb raw reads(數(shù)據NCBI查詢號:PRJNA915171),質控后得到64.15 Gb clean reads,所有文庫raw_bases中Q phred數(shù)值大于30(Q30)的堿基數(shù)占總堿基的百分比均超過93%,表明所有樣本文庫均獲得了高質量的raw reads。由于缺乏參考基因組,數(shù)據經拼接和去冗余后,最終獲得61 755條非冗余unigenes,平均長度為971.5 bp,50值為1399 bp,其中38 070條(61.65%)的長度大于500 bp,19 867條(32.19%)的長度大于1000 bp。BUSCO評估轉錄本完整性,結果表明本課題組獲得的unigenes集覆蓋了BUSCO庫中的87.9%,其中完全匹配且匹配上一個(complete and single-copy BUSCOs,S)的比例為77.2%。Raw reads質量和BUSCO評估結果表明,本研究獲得的轉錄組數(shù)據是可靠和有效的。

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

3.3 轉錄本注釋和功能歸類

61 755條非冗余unigenes中29 837條(48.32%)、21 328條(34.54%)、6467條(10.47%)、17 726條(28.70%)、26 812條(43.42%)、18 557條(30.05%)、16 641條(26.95%)分別注釋到NR、Swissprot、KEGG、KOG、eggNOG、GO和Pfam庫。經去冗余分析,共30 263條(49.0%)unigenes至少在1個數(shù)據庫中被成功注釋到,3234條(5.23%)unigenes在7個數(shù)據庫中均注釋到。目前,重樓屬植物所在目百合目(Liliales)中沒有一個物種進行過基因組測序,因此本研究中超過50%的unigenes不能被成功注釋是可以理解的,這些未注釋的基因可能涉及未知的功能。

18557條注釋到了GO數(shù)據庫的unigenes中,注釋到“metabolic process”和“catalytic actiVIty”的條數(shù)占比相對高于其他GO子類別(圖2-A),這表明球藥隔重樓中相當一部分基因參與物質代謝。6467條注釋到了KEGG數(shù)據庫的unigenes中,4066條(62.87%)注釋為“Metabolism”,其中歸類為“Metabolism of terpenoids and polyketides”的條數(shù)為213,進一步分析發(fā)現(xiàn),70條注釋為參與萜類化合物骨架的生成,39條注釋為參與三萜的生物合成,這些基因可能與螺甾烷型重樓皂苷的生物合成相關(圖2-B)。

3.4 差異基因分析

根莖莖組有14 056個,其中5344個上調,8712個下調;根莖葉組有18 724個,其中8259個上調,10 465個下調;葉莖組有10806個,其中5026個上調,5780個下調(圖3-A)。分別對上調和下調的進行KEGG富集分析,結果表明,“photosynthesis”“porphyrin and chlorophyll metabolism”“photosynthesis-antenna proteins”等植物光合作用、生長發(fā)育相關途徑的基因,主要在葉和莖中富集;“flavonoid biosynthesis”“flavone and flavonol biosynthesis”“carotenoid- biosynthesis pathway”等次生代謝途徑基因,也主要在葉和莖中富集(圖3-B、C)。而與萜類化合物上游代謝密切相關的“terpenoid backbone biosynthesis”的,主要在根莖中富集,與三萜重樓甾體皂苷下游代謝密切相關“steroid biosynthesis”的,也主要在根莖中富集(圖3-B、C),這與重樓皂苷主要在根莖中積累一致(圖1)。這些結果表明,球藥隔重樓中螺甾烷型重樓皂苷生物合成密切相關的基因可能屬于這些。

圖2 Unigenes的GO (A) 和KEGG (B) 注釋

3.5 螺甾烷型重樓皂苷生物合成途徑已鑒定基因分析

基于KEGG注釋和蛋白編碼區(qū)(ORF)完整性結果,篩選到26條unigenes參與植物三萜類化合物共同前體角鯊烯的生成。其中8條基因編碼MVA途徑6個關鍵酶[8-9],包括2條乙酰輔酶A轉乙酰酶(acetyl-CoA acetyltransferase,AACT)、2條甲羥戊二酰輔酶A合成酶(hydroxy-methylglutaryl-CoA synthase,HMGCS)、1條羥甲基戊二酰輔酶a還原酶(hydroxymethylglutaryl-CoA reductase,HMGCR)、1條甲羥戊酸激酶(mevalonate kinase,MVK)、1條磷酸甲羥戊酸激酶(phosphomevalonate kinase,mvaK2)和1條mevalonate diphosphate decarboxylase(MVD),這6個關鍵酶最終催化形成異戊烯基焦磷酸(isopentenyl pyrophosphate,IPP);10條基因編碼MEP途徑7個關鍵酶[8-9],包括3條1-脫氧--酮糖5-磷酸合成酶(1-deoxy--xylulose 5-phosphate synthase,DXS)、1條1-脫氧--酮糖5-磷酸還原異構酶(1-deoxy--xylulose 5-phosphate reductoisomerase, DXR)、1條2--甲基--赤蘚糖醇4-磷酸胞苷基轉移酶(2--methyl--erythritol 4-phosphate cytidylyl- transferase,ispD)、2條4-二磷酸胞苷-2-甲基--赤蘚糖醇激酶(4-diphosphocytidyl-2--methyl--erythritol kinase,ispE)、1條2--甲基--赤蘚糖醇2,4-環(huán)二磷酸合酶(2--methyl--erythritol 2,4-cyclodiphosphate synthase,ispF)、1條()-4-羥基-3-甲基丁-2-烯基二磷酸合成酶[()-4-hydroxy-3-methylbut-2- enyl-diphosphate synthase,ispG]、1條()-4-羥基-3-甲基丁-2-烯基二磷酸還原酶[()-4-hydroxy-3- methylbut-2-enyl diphosphate reductase,ispH],這7個關鍵酶最終催化形成IPP和二甲基丙烯基二磷酸(dimethylpropyl diphosphate,DMAPP);此外,還包括了1條IDI,該基因可使上述MEP和MVA途徑生成的IPP與DMAPP間相互轉化,以及5條法呢基焦磷酸合酶(farnesyl diphosphate synthase,F(xiàn)PPS),F(xiàn)PPS可利用IPP與DMAPP為底物,生成倍半萜的共同前體FPP,最終2分子的FPP在SQS的催化下,形成植物三萜類化合物共同前體角鯊烯,本研究注釋到了1條基因(圖4)。基因表達量結果表明,MVA途徑基因和基因在根莖中的表達量顯著高于莖和葉,而MEP途徑基因在葉中高度表達,而螺甾烷型重樓皂苷主要在根莖中富集,這表明球藥隔重樓中螺甾烷型重樓皂苷生物合成所需的IPP和DMAPP可能主要來自MVA途徑,這與三萜類化合物生物合成主要由MVA途徑參與的結果一致[29]。

A-DEGs統(tǒng)計柱狀圖 B-表達上調DEGs的KEGG注釋 C-表達下調DEGs的KEGG注釋

基于已發(fā)表催化形成螺甾烷型甾體母核的文獻報道(表1),利用同源基因blast的檢索方法,去除不完整ORF及同源性≤80%的unigenes,最終篩選到17條unigenes參與后續(xù)螺甾烷型重樓皂苷母核“薯蕷皂苷元”的生成,其中2條SQE、1條CAS、1條SMO1-3、1條CPI-5、1條CYP51G、1條C14-R-2、1條8,7SI-4、1條C5-SD1、2條7-DR1、1條SSR1-3、2條CYP90G4、1條CYP94D108、1條CYP94D109和1條CYP72A616,在這些基因的共同催化下,最終生成了薯蕷皂苷元(圖4)。而參與皂苷元側鏈寡糖基的催化基因,報道了葡糖基化皂苷元C3-OH的UGT73CR1、UGT80A40和UGT80A41,2′--鼠李糖基化的UGT73CE1,以及6′--葡萄糖基化的UGT91AH1和UGT91AH2,本研究分別同源blast到了1條相關基因的同源序列。基因表達量結果表明,下游基因主要在根莖中高表達,這與螺甾烷型重樓皂苷主要在根莖中富集的結果一致;而末端的、、、、、和基因,在根莖中的表達量相對較低,可能是這些基因除參與螺甾烷型重樓皂苷合成外,還有其他更重要的生物學功能。此外,本研究在球藥隔重樓中檢索到的同源基因的編碼區(qū)氨基酸序列,與在云南重樓中驗證的同源基因高度同源,表明螺甾烷型重樓皂苷生物合成途徑基因在重樓屬植物中高度保守。

實線箭頭為已驗證的步驟,虛線為未驗證步驟

表1 同源基因檢索所用基因清單及相關信息

3.6 螺甾烷型重樓皂苷與其生物合成途徑基因相關性分析

使用pearson相關性計算方法,對螺甾烷型重樓皂苷的積累水平與其生物合成途徑基因表達水平進行了關聯(lián)性分析,結果見圖5。由圖可知,、、、、、、、、、、、、、、和等17個已知途徑中的基因表達水平與螺甾烷型重樓皂苷的積累存在顯著的正相關性,而基因的表達水平似乎負調控螺甾烷型重樓皂苷的積累。

3.7 螺甾烷型重樓皂苷生物合成途徑中UGT的篩選與分析

螺甾烷型重樓皂苷的結構差異主要來自側鏈的寡糖基,寡糖基的形成可能需要多個基因的參與,但僅被研究報道[12]。因此,我們結合上述皂苷的定量數(shù)據、轉錄組數(shù)據和已鑒定合成途徑基因分析結果,篩選潛在參與螺甾烷型重樓皂苷糖基化的基因。首先,基于球藥隔重樓unigenes的CDS數(shù)據庫和UGT的HMM文件(PF00201),共檢索得到75條UGTs,其中61條有完整的ORF。在糖基轉移酶(GT)命名系統(tǒng)指導下[30],結合擬南芥和玉米中UGTs的分類[31-32],對61條PfUGTs進行了系統(tǒng)命名和分類,最終將PfUGTs分類到16個UGT家族,即UGT72、UGT73、UGT74、UGT75、UGT78、UGT80、UGT83、UGT84、UGT85、UGT87、UGT88、UGT89、UGT90、UGT91、UGT92和UGT93。其次,對61條在球藥隔重樓根莖、莖和葉中的表達模式進行分析,結果表明,的表達模式可分為根莖中相對高表達(I類)、莖中相對高表達(II類)、莖葉中均相對高表達(III類)和葉中相對高表達(IV類)4類,分別有15條、12條、8條和26條基因(圖5)。螺甾烷型重樓皂苷定量分析和皂苷母核合成途徑基因表達分析結果表明,螺甾烷型重樓皂苷主要在根莖中積累(圖1),同時皂苷母核合成途徑關鍵基因在根莖中相對高表達(圖4)。因此,推測參與皂苷母核糖基修飾的可能也在根莖中相對高表達,故I類中的15條可能參與了皂苷母核后期的糖基化修飾。此外,重樓皂苷I、H和VII在葉片中也有較高的積累水平,特別是重樓皂苷VII在葉中積累最高,且在重樓屬植物中功能驗證的、1、和基因在葉片中高表達,故認為在葉片中相對高表達且與、1、和存在相同表達模式的34條Ⅲ和Ⅳ類,也是潛在的候選基因。植物源UGT的研究表明,參與三萜類化合物糖基化的UGTs主要屬于A組的UGT91和UGT94、D組的UGT73、E組的UGT71、L組的UGT74和UGT75、P組的UGT720家族[32]。上述篩選到的49條PfUGTs(I、Ⅲ和Ⅳ類)中,有2條屬于UGT91家族,有24條屬于UGT73家族,有3條屬于UGT74家族,有1條屬于UGT75家族(圖6)。

圖5 螺甾烷型重樓皂苷與合成途徑基因表達相關性分析

進一步將61條PfUGTs與其他物種中已鑒定功能的UGTs構建混合系統(tǒng)發(fā)育樹,發(fā)現(xiàn)來自苦蕎麥L.和矮牽牛花L.能識別UDP-鼠李糖作為糖基供體的FeUGT79A8[33]和PhUGT79A1[34],來自黃芩L.能利用UDP-葡萄糖醛酸作為糖供體的UBGAT[35],以及來自擬南芥可利用UDP-木糖作為糖供體的AtUGT78D1[36],與PfUGTs成功構建了混合發(fā)育樹(圖7)。這表明與FeUGT79A8和PhUGT79A1高度聚合的A組PfUGT91家族基因,可能是球藥隔重樓中其他識別UDP-鼠李糖作為糖基供體的UGTs;與UBGAT高度聚合的E組PfUGT72和PfUGT88家族基因可能是識別UDP-葡萄糖醛酸的UGTs;與AtUGT78D1高度聚合的F組PfUGT78D1基因可能是識別UDP-木糖的UGTs。有趣的是,螺甾烷型重樓皂苷的糖基化十分復雜,包括葡萄糖基化、鼠李糖基化、木糖化等,故這些家族的基因也是值得關注的。

各組織表達量值為3個重復的平均值;紅星號的為重樓屬植物中功能驗證的基因

Expression values for each tissue are the average of three replicates; Red asterisks are functionally verifiedgenes in the plant in

圖6在根莖、莖和葉中的熱圖分析

Fig. 6 Heatmap ofin rhizome, stem and leaf tissues

紅色星號標識的是來自其他物種的UGT,綠色正方形表示bootstrap值超過50%

3.8 qRT-PCR驗證基因表達的可靠性

為了確認RNA-seq結果的可靠性,選取了螺甾烷型重樓皂苷合成途徑中的5個關鍵酶基因(、、、和),以及篩選的3個基因(、和)進行qRT-PCR定量分析。qRT-PCR的引物見表3。

實驗結果表明,8個基因的qRT-PCR定量結果和轉錄組測序獲得的FPKM值在球藥隔重樓中根莖、莖和葉中的表達趨勢相似(圖8),這表明轉錄組測序結果是可靠的。

表3 qRT-PCR引物序列及產物大小

圖8 qRT-PCR驗證轉錄組測序數(shù)據

4 討論

通常,植物合成特定的次生代謝產物幫助其適應逆境,如食草動物損傷和微生物感染可以刺激植物防御類化合物的積累[37]。三萜皂苷類次生產物具抗菌、殺蟲、抗草食、植物毒性等廣泛的生物活性[1,3,38]。相比其他重樓屬植物,球藥隔重樓更喜生長于易被微生物和蚊蟲感染的陰濕環(huán)境,富含螺甾烷型三萜皂苷類成分可能是其遺傳特性和環(huán)境互做的結果。參與螺甾烷型三萜皂苷生物合成的基因表達和該類成分的HPLC定量結果,反映了球藥隔重樓莖、葉和根莖在該類皂苷生物合成和積累中的不同作用。本研究中,螺甾烷型三萜皂苷總體上在根莖中積累水平比較高,這與云南重樓等其他同屬植物測定結果一致[6-7,26],但偏諾皂苷在球藥隔重樓中積累更加豐富,如其他重樓中很難檢測到的重樓皂苷VI[7,39]。此外,莖和葉中也積累了一定數(shù)量的三萜皂苷,特別是重樓皂苷VII。因此,球藥隔重樓是偏諾皂苷類成分來源的理想植物材料,同時,地上部分具有資源替代利用的潛力。

2010年,Pellicer等[40]用流式細胞儀分析日本重樓,發(fā)現(xiàn)其基因組大小高達約150 Gb,是最大的真核生物基因組之一,可見重樓屬植物在植物進化過程中占有舉足輕重的作用,巨大的基因組也制約了該屬植物基因組的獲取,因此,關于該屬植物基因層面的研究主要集中于轉錄組和葉綠體基因組等。Gao等[26]利用RNA-seq技術對云南重樓營養(yǎng)期和果期的葉和根莖分別進行轉錄組測序,拼接獲得了85 677條unigenes,N50為668 bp,僅34 020條(39.71%)unigenes被成功注釋。Yang等[41]同樣利用RNA-seq對華重樓花期的根、根莖、莖、葉和花進行轉錄組測序,獲得了57 537條unigenes,N50為664 bp,平均長度為992 bp,有33 449條(58.13%)被成功注釋。本研究獲得了球藥隔重樓的61755條非冗余unigenes,平均長度為971.5 bp,N50為1399 bp,30 263條(49.0%)unigenes被成功注釋,三萜皂苷生物合成途徑基因分析中,除基因未檢測到外,其余通路基因均同源blast到,以上結果表明本研究的RNA-seq轉錄組數(shù)據,用于研究重樓三萜皂苷合成途徑基因是科學可靠的,亦可為今后該品種的分子研究提供一定的數(shù)據支撐。

三萜重樓甾體皂苷下游合成密切相關的“steroid biosynthesis”的主要在根莖中富集,合成途徑下游關鍵基因也主要在根莖中表達;“terpenoid backbone biosynthesis”途徑和“steroid biosynthesis”途徑不僅是三萜重樓甾體皂苷生物合成的中間步驟,也是類胡蘿卜素生物合成不可缺少的步驟[42],而類胡蘿卜素生物合成相關基因主要在葉和莖中高表達。此外,油菜素甾醇代謝途徑(brassinosteroid biosynthesis)也需要消耗三萜甾體途徑的中間產物[43],而“Brassinosteroid biosynthesis”途徑基因主要在葉和莖中高表達。這些因素可能限制了葉片和莖生物合成重樓甾體皂苷的能力。三萜甾體皂苷生物合成途徑中的關鍵酶基因研究表明,MVA途徑中的HMGCR和甾體類化合物下游代謝總開關的CAS為代謝途徑上的限速酶[44];此外,最近在本氏煙草中異源合成膽甾烯醇的研究表明,SSR1-3、SMO1-3、CPI-5、CYP51G、SMO2-2和C5-SD1對膽甾烯醇的合成有顯著影響[14];而本研究中除上述8個關鍵酶外,發(fā)現(xiàn)IDI、8,7SI-4、CYP90G4、C14-R-2、CYP94D108、HMGCS、SQS、ispH和DXR也可能顯著影響著重樓三萜甾體皂苷的積累,這些關鍵酶可作為實現(xiàn)人工精準調控重樓皂苷積累的潛在靶標基因。Hua等[12]利用PacBio三代測序平臺對云南重樓進行了全長轉錄組測序,從39 875條unigenes中鑒定到了199條UGTs,系統(tǒng)發(fā)育分析將其分到11個不同的組(A~E、I、J、L、O、P和OG),其中D組的UGT73家族成員最多。本研究中,僅鑒定到61條UGTs,這可能是重樓屬基因在不同生長期的轉錄水平有巨大差異導致(本研究僅檢測了果期,而全長轉錄組檢測多個時期),也可能是因球藥隔重樓和云南重樓的遺傳背景差異所致。此外,本研究中61條UGTs劃分到了13個不同的組(A~G、I、J、L、M、O和OG),較云南重樓多了F、G和M組,缺失了P組,這表明球藥隔重樓基因較云南重樓發(fā)生了明顯的分化,可能是該植物三萜皂苷較同屬其他植物差異的原因之一,球藥隔重樓中D組UGT73家族的成員也是最多的,這與前人研究結果一致[12, 32]。已有研究表明,OG組的UGT成員主要負責修飾甾醇和脂類化合物[45-46],螺甾烷型重樓皂苷屬于三萜甾醇類皂苷,因此該組的UGT基因(即UGT80家族)也是值得重點關注的。

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

[1] Gupta D D, Mishra S, Verma S S,. Evaluation of antioxidant, anti-inflammatory and anticancer actiVIties of diosgenin enrichedrhizome extract of Indian Himalayan landraces [J]., 2021, 270: 113842.

[2] 中國藥典[S]. 一部. 2020: 108.

[3] Puwein A, Thomas S C. An overVIew of, a highly vulnerable medicinal herb of eastern Himalayan region for sustainable exploitation [J]., 2020, 10(1): 3-14.

[4] Ding Y G, Zhao Y L, Zhang J,. The traditional uses, phytochemistry, and pharmacological properties ofL. (Liliaceae): A review [J]., 2021, 278: 114293.

[5] 管鑫, 李若詩, 段寶忠, 等. 重樓屬植物化學成分、藥理作用研究進展及質量標志物預測分析 [J]. 中草藥, 2019, 50(19): 4838-4852.

[6] Man S L, Gao W Y, Zhang Y J,. Qualitative and quantitative determination of major saponins inandby HPLC-ELSD and HPLC-MS/MS [J]., 2010, 878(29): 2943-2948.

[7] 劉巨釗, 馬苗苗, 田曉黎, 等. 基于文獻計量學的重樓研究現(xiàn)狀及熱點分析[J]. 藥物評價研究, 2022, 45(11): 2347-2356.

[8] Tholl D. Biosynthesis and biological functions of terpenoids in plants [J]., 2015, 148: 63-106.

[9] Wang J F, Li S Y, Xiong Z Q,. Pathway mining-based integration of critical enzyme parts forbiosynthesis of steviolglycosides sweetener in[J]., 2016, 26(2): 258-261.

[10] Gao J X, Chen Y G, Li D S,. Cloning and functional characterization of a squalene synthase fromvar.[J]., 2021, 18(7): e2100342.

[11] Laranjeira S, Amorim-Silva V, Esteban A,.squalene epoxidase 3 (SQE3) complements SQE1 and is important for embryo development and bulk squalene epoxidase activity [J]., 2015, 8(7): 1090-1102.

[12] Hua X, Song W, Wang K Z,. Effective prediction of biosynthetic pathway genes involved in bioactive polyphyllins in[J]., 2022, 5(1): 50.

[13] Darnet S, Bard M, Rahier A. Functional identification of sterol-4alpha-methyl oxidase cDNAs fromby complementation of a yeast erg25 mutant lacking sterol-4alpha-methyl oxidation [J]., 2001, 508(1): 39-43.

[14] Yin X, Liu J, Kou C X,. Deciphering the network of cholesterol biosynthesis inlaid a base for efficient diosgenin production in plant chassis [J]., 2023, 76: 232-246.

[15] Lovato M A, Hart E A, Segura M J,. Functional cloning of ancDNA encoding cycloeucalenol cycloisomerase [J]., 2000, 275(18): 13394-13397.

[16] Kushiro M, Nakano T, Sato K,. Obtusifoliol 14α-demethylase (CYP51) antisenseshows slow growth and long life [J]., 2001, 285(1): 98-104.

[17] Schrick K, Mayer U, Horrichs A,. FACKEL is a sterol C-14 reductase required for organized cell division and expansion inembryogenesis [J]., 2000, 14(12): 1471-1484.

[18] Grebenok R J, Ohnmeiss T E, Yamamoto A,. Isolation and characterization of anC-8, 7 sterol isomerase: Functional and structural similarities to mammalian C-8, 7 sterol isomerase/emopamil-binding protein [J]., 1998, 38(5): 807-815.

[19] Husselstein T, Schaller H, Gachotte D,. Delta7-sterol-C5-desaturase: Molecular characterization and functional expression of wild-type and mutant alleles [J]., 1999, 39(5): 891-906.

[20] Lecain E, Chenivesse X, Spagnoli R,. Cloning by metabolic interference in yeast and enzymatic characterization ofsterol delta 7-reductase [J]., 1996, 271(18): 10866-10873.

[21] Choe S, Dilkes B P, Gregory B D,. Thedwarf1 mutant is defective in the conversion of 24-methylenecholesterol to campesterol in brassinosteroid biosynthesis [J]., 1999, 119(3): 897-907.

[22] Christ B, Xu C C, Xu M L,. Repeated evolution of cytochrome P450-mediated spiroketal steroid biosynthesis in plants [J]., 2019, 10(1): 3206.

[23] ChenY, YanQ, JiY,. Unraveling the serial glycosylation in the biosynthesis of steroidal saponins in the medicinal plantand their antifungal action [J]., 2023, 5: 33.

[24] Zhao W S, Gao W Y, Huang L Q,. Study on the chemical constituents of the[J]., 2011,23 (6) : 1017-1020.

[25] 汪瑤, 薛丹, 文飛燕, 等. HPLC測定球藥隔重樓中的5種重樓皂苷 [J]. 華西藥學雜志, 2015, 30(4): 469-470.

[26] Gao X Y, Zhang X, Chen W,. Transcriptome analysis ofvar.illuminates the biosynthesis and accumulation of steroidal saponins in rhizomes and leaves [J]., 2020, 178: 112460.

[27] Hou L X, Zhang F R, Yuan X C,. Comparative transcriptome analysis reveals key genes for polyphyllin difference in fivespecies [J]., 2022, 174(6): e13810.

[28] 尹顯梅, 張開元, 蔣桂華, 等. 華重樓皂苷類成分的動態(tài)分布規(guī)律對藥材質量的影響 [J]. 中草藥, 2017, 48(6): 1199-1204.

[29] Kempinski C, Chappell J. Engineering triterpene metabolism in the oilseed of[J]., 2019, 17(2): 386-396.

[30] MacKenzie P I, Owens I S, Burchell B,. The UDP glycosyltransferase gene superfamily: Recommended nomenclature update based on evolutionary divergence [J]., 1997, 7(4): 255-269.

[31] Li Y, Baldauf S, Lim E K,. Phylogenetic analysis of the UDP-glycosyltransferase multigene family of[J]., 2001, 276(6): 4338-4343.

[32] Wilson A E, Tian L. Phylogenomic analysis of UDP-dependent glycosyltransferases provides insights into the evolutionary landscape of glycosylation in plant metabolism [J]., 2019, 100(6): 1273-1288.

[33] Koja E, Ohata S, Maruyama Y,. Identification and characterization of a rhamnosyltransferase involved in rutin biosynthesis in(common buckwheat) [J]., 2018, 82(10): 1790-1802.

[34] Brugliera F, Holton T A, Stevenson T W,. Isolation and characterization of a cDNA clone corresponding to the Rt locus of[J]., 1994, 5(1): 81-92.

[35] Nagashima S, Hirotani M, Yoshikawa T. Purification and characterization of UDP-glucuronate: Baicalein 7-- glucuronosyltransferase fromGeorgi. cell suspension cultures [J]., 2000, 53(5): 533-538.

[36] Pandey R P, Malla S, Simkhada D,. Production of 3--xylosyl quercetin in[J]., 2013, 97(5): 1889-1901.

[37] Yang C Q, Fang X, Wu X M,. Transcriptional regulation of plant secondary metabolism [J]., 2012, 54(10): 703-712.

[38] WAller G R, Yamasaki K.[M]. Berlin: Springer Science & Business Media, 2013.

[39] 張紹山, 劉璇, 王景富, 等. UPLC法測定云南省不同地區(qū)云南重樓及多芽品系中7種甾體皂苷量及其指紋圖譜建立 [J]. 中草藥, 2016, 47(23): 4257-4263.

[40] Pellicer J, Fay M F, Leitch I J. The largest eukaryotic genome of them all? [J]., 2010, 164(1): 10-15.

[41] Yang Z Y, Yang L F, Liu C K,. Transcriptome analyses ofvar. chinensis,, andcharacterize their steroidal saponin biosynthesis pathway [J]., 2019, 135: 52-63.

[42] Sandmann G. Diversity and origin of carotenoid biosynthesis: Its history of coevolution towards plant photosynthesis [J]., 2021, 232(2): 479-493.

[43] Hou L X, Yuan X C, Li S,. Genome-wide identification ofgene family and expression patterns related to jasmonic acid treatment and steroidal saponin accumulation in[J]., 2021, 22(20): 10953.

[44] 尹艷, 關紅雨, 張夏楠. 甾體皂苷生物合成相關酶及基因研究進展 [J]. 天然產物研究與開發(fā), 2016, 28(8): 1332-1336.

[45] Paquette S, M?ller B L, Bak S. On the origin of family 1 plant glycosyltransferases [J]., 2003, 62(3): 399-413.

[46] Stucky D F, Arpin J C, Schrick K. Functional diversification of two UGT80 enzymes required for steryl glucoside synthesis in[J]., 2015, 66: 189-201.

Identification of genes involved in biosynthesis of spirostane-type polyphyllin inbased on transcriptome analysis

QUBIE A-xiang1, 2, 3, QUBIE Jun-zhang1, 2, 3, FENG Ming-kang1, 2, 3, HE Bin1, 2, 3, HE Xiao-yong1, 2, 3, LI Juan2, 3, FENG Jing-qiu2, 3, SHENG Hua-chun2, 3, HUANG Yan-fei1, 3, 4, YAN Xin-jia2, 3, LIU Yuan2, 3, ZHANG Shao-shan2, 3

1. College of Pharmacy, Southwest Minzu University, Chengdu 610041, China 2. Sichuan Technology and Engineering Laboratory of Qiang-Yi Medicinal Resources Protection and Utilization, Chengdu 610225, China 3. Key Laboratory of Tibetan Plateau Ethnic Medicinal Resources Protection and Utilization of National Ethnic Affairs Commission, Chengdu 610225, China

To investigate the biosynthetic genes of spirostane-type polyphyllins of.Based on the rhizome, stem and leaf materials of, the contents of six spirostane-type polyphyllins were analyzed by HPLC. The transcriptome sequencing was performed using Illumina HiSeq X Ten platform. The quantitative data and transcriptome data were combined for bioinformatics analysis.The content of six spirostane-type polyphyllins had significant differences in the rhizome, stem and leaf of. A total of 61755 non-redundant unigenes were obtained by Illumina sequencing, of which 30263 (49.0%) were successfully annotated. Thirty-one key genes in the biosynthesis pathway of diosgenin had been successfully identified, and the gene expression pattern was generally consistent with the accumulation level of spirostane-type polyphyllins. The accumulation level of spirostane-type polyphyllins were positively correlated with,,,,,,,,,,,,,,and. A total of 61 uridine diphosphate glycosyltransferase () genes were obtained, of which 30 were potentially involved in the glycosylation modification of spirostane-type polyphyllins.There was a consistency between the expression pattern of genes involved in the biosynthesis pathway of polyphyllins and the accumulation level of polyphyllins in the rhizome, stems and leaves of. The genes identified in this study can help to elucidate the biosynthetic pathway of this kinds of active products.

Franch.; transcriptome; spirostane-type polyphyllins; biosynthesis pathway; uridine diphosphate glycosyltransferase

R286.12

A

0253 - 2670(2023)20 - 6798 - 15

10.7501/j.issn.0253-2670.2023.20.023

2023-02-03

四川省自然科學基金資助項目(22NSFSC3136);西南民族大學2023年研究生創(chuàng)新型項目資助碩士重點項目(ZD2023459)

曲別阿香(1999—),女,在讀碩士,主要從事民族藥品質評價及品質形成分子機理研究。E-mail: 1814910468@qq.com

通信作者:張紹山,助理研究員,主要從事民族藥資源品質評價、生態(tài)栽培及次生代謝產物生物調控研究。E-mail: shaoshanzhang_smu@qq.com

劉 圓,教授,主要從事民族藥資源品質評價研究。E-mail: 499769896@qq.com

[責任編輯 時圣明]

猜你喜歡
途徑生物
生物多樣性
天天愛科學(2022年9期)2022-09-15 01:12:54
生物多樣性
天天愛科學(2022年4期)2022-05-23 12:41:48
上上生物
當代水產(2022年3期)2022-04-26 14:26:56
發(fā)現(xiàn)不明生物
科學大眾(2021年9期)2021-07-16 07:02:54
史上“最黑暗”的生物
軍事文摘(2020年20期)2020-11-28 11:42:50
第12話 完美生物
航空世界(2020年10期)2020-01-19 14:36:20
構造等腰三角形的途徑
多種途徑理解集合語言
減少運算量的途徑
醫(yī)保基金“可持續(xù)”的三條途徑
主站蜘蛛池模板: 国产精品亚洲va在线观看| 国产产在线精品亚洲aavv| 5388国产亚洲欧美在线观看| 精品色综合| 亚洲视频无码| 午夜少妇精品视频小电影| 日韩经典精品无码一区二区| 毛片免费试看| 美女国产在线| 97久久精品人人做人人爽| 在线看国产精品| 亚洲成综合人影院在院播放| 欧美色综合久久| 三级欧美在线| 99er这里只有精品| 色综合网址| 亚洲男人天堂2020| 国产sm重味一区二区三区| 亚洲中文字幕97久久精品少妇| 久久久久无码国产精品不卡| 国产本道久久一区二区三区| 成人欧美日韩| 亚洲人成成无码网WWW| 日本一本正道综合久久dvd | 中文字幕在线永久在线视频2020| 996免费视频国产在线播放| 四虎精品黑人视频| 欧美高清国产| 久久精品66| 亚洲人成网站色7777| 激情影院内射美女| 噜噜噜久久| 高清久久精品亚洲日韩Av| 91po国产在线精品免费观看| 日本影院一区| 午夜综合网| 福利视频一区| 国产一级毛片yw| 狠狠色香婷婷久久亚洲精品| 性网站在线观看| 日韩精品一区二区三区中文无码| 538精品在线观看| 波多野结衣爽到高潮漏水大喷| 综合色亚洲| 99伊人精品| 国产欧美亚洲精品第3页在线| 色悠久久综合| 真实国产乱子伦视频| 特级毛片免费视频| 欧美精品黑人粗大| 国产成人在线无码免费视频| 99视频在线免费| 国产系列在线| 日韩欧美高清视频| 亚洲资源站av无码网址| 美女无遮挡免费网站| 噜噜噜久久| 日韩AV无码一区| 欧美 亚洲 日韩 国产| 九月婷婷亚洲综合在线| 最近最新中文字幕免费的一页| 国产清纯在线一区二区WWW| 日韩成人高清无码| 国产一区二区精品高清在线观看| 国产成人免费高清AⅤ| 91麻豆国产视频| 扒开粉嫩的小缝隙喷白浆视频| 91久久精品国产| 国产激情无码一区二区APP | 国产视频一二三区| 精品久久久久成人码免费动漫| 精品撒尿视频一区二区三区| 亚洲精品成人片在线播放| 精品少妇人妻av无码久久| 欧美中文一区| 美女高潮全身流白浆福利区| 国产97公开成人免费视频| 免费日韩在线视频| 在线一级毛片| 亚洲AⅤ波多系列中文字幕| 免费毛片在线| 成人小视频在线观看免费|