許小涵，唐志強(qiáng)，劉 謙, 2，李 佳, 2，劉振華, 2，張永清, 2，蒲高斌, 2*

? 藥材與資源 ?

干旱脅迫下忍冬全基因組DNA甲基化和轉(zhuǎn)錄組分析

許小涵1，唐志強(qiáng)1，劉 謙1, 2，李 佳1, 2，劉振華1, 2，張永清1, 2，蒲高斌1, 2*

1. 山東中醫(yī)藥大學(xué)，山東 濟(jì)南 250355 2. 山東省中藥質(zhì)量控制與全產(chǎn)業(yè)鏈建設(shè)協(xié)同創(chuàng)新中心，山東 濟(jì)南 250355

分析干旱脅迫下忍冬全基因組DNA甲基化水平變化、模式以及與基因表達(dá)的關(guān)系，為進(jìn)一步探究忍冬響應(yīng)干旱脅迫的分子機(jī)制奠定基礎(chǔ)。采用全基因組重亞硫酸鹽甲基化測序法（whole-genome bisulfite sequencing，WGBS），測定干旱脅迫下忍冬全基因組DNA甲基化水平，并聯(lián)合轉(zhuǎn)錄組測序結(jié)果對差異甲基化基因相關(guān)的差異表達(dá)基因進(jìn)行分析。干旱脅迫下忍冬整體甲基化水平普遍上升；其中CG類型甲基化水平明顯高于CHG或CHH甲基化類型；對不同干旱脅迫處理下的忍冬進(jìn)行差異甲基化區(qū)域（differentially methylated regions，DMR）分析，發(fā)現(xiàn)1935、2158和1750個差異甲基化相關(guān)基因，基因組百科全書（Kyoto encyclopedia of genes and genomes，KEGG）富集分析表明，顯著富集的通路主要包括亞油酸代謝，氮代謝通路和次生代謝物的生物合成通路等。通過轉(zhuǎn)錄組測序技術(shù)，發(fā)現(xiàn)差異甲基化基因相關(guān)的差異表達(dá)基因主要富集于次生代謝產(chǎn)物合成相關(guān)通路，同時篩選出4個在干旱脅迫下基因表達(dá)量上升顯著的忍冬MYB（MYB，LjMYB）轉(zhuǎn)錄因子。干旱脅迫下忍冬甲基化水平上升，推測DNA甲基化調(diào)控基因表達(dá)是干旱影響金銀花有效成分生物合成的作用機(jī)制之一。

忍冬；干旱脅迫；DNA甲基化；全基因組重亞硫酸鹽甲基化測序；轉(zhuǎn)錄組測序

金銀花為忍冬科忍冬屬植物忍冬Thunb.的干燥花蕾或帶初開的花，具有清熱解毒、疏散風(fēng)熱之功效[1]，為常用大宗中藥材之一。目前，金銀花商品藥材主要來自人工栽培，其中山東省的種植面積近6.667萬公頃，為主要道地產(chǎn)區(qū)?，F(xiàn)代生物學(xué)認(rèn)為，道地藥材的本質(zhì)是藥用植物所擁有的基因型受環(huán)境因子誘導(dǎo)后表達(dá)的產(chǎn)物[2]。受環(huán)境影響，金銀花藥材質(zhì)量差異較大。如金銀花“九豐一號”由山東曲阜引種至廣西武鳴后，在相同采摘時期和干燥條件下，其木犀草苷含量下降了57.4%[3]。另有研究表明，水分[4]、光照[5]、溫度[6-7]等環(huán)境因子對金銀花次生代謝產(chǎn)物的生物合成均有顯著影響。然而，環(huán)境因子究竟是如何通過基因型修飾而發(fā)生作用的，特別是環(huán)境因子影響金銀花有效成分生物合成的作用機(jī)制，至今缺少直接的實驗研究的揭示和證明。

表觀遺傳學(xué)研究為揭示環(huán)境與藥材品質(zhì)的關(guān)系提供了契機(jī)。表觀遺傳學(xué)是以不改變基因DNA序列編碼的方式來改變遺傳基因表達(dá)的一種遺傳方式，主要涉及DNA甲基化作用的改變、RNA沉默、染色質(zhì)組蛋白的修飾作用、基因印記[8]。DNA甲基化是表觀遺傳學(xué)的重要研究內(nèi)容之一，與基因表達(dá)、基因組防御、細(xì)胞分化、染色質(zhì)失活和基因組印記的調(diào)節(jié)有關(guān)[9]。研究發(fā)現(xiàn)紅光和遠(yuǎn)紅光照射沉香后，其葫蘆素含量發(fā)生改變，通過全基因組重亞硫酸氫鹽測序發(fā)現(xiàn)紅光和遠(yuǎn)紅光條件下甲基化水平變化較大，紅光可能通過改變沉香中CHH和CHG的甲基化來調(diào)控基因的表達(dá)[10]。謝宜杰等[11]通過MASP技術(shù)測定廣藿香DNA甲基化水平，能夠在不同來源的樣品中鑒別出石牌產(chǎn)地的廣藿香。Li等[12]采用HPLC法測定紅豆杉不同傳代培養(yǎng)細(xì)胞中全基因甲基化水平，發(fā)現(xiàn)DNA總體甲基化水平高低與紫杉醇的量呈顯著的負(fù)相關(guān)關(guān)系，使用5-azaC處理后細(xì)胞中紫杉醇積累量顯著增加。Zha等[13]研究發(fā)現(xiàn)忍冬Thunb. 和其自然誘變產(chǎn)生的變種紅白忍冬Thunb. var.(Wats.) Bak.綠原酸含量在兩者之間差異顯著，通過測定DNA甲基化水平發(fā)現(xiàn)紅白忍冬苯丙氨酸裂解酶2基因側(cè)翼區(qū)域的啟動子?109～?279 bp處DNA甲基化程度要高于忍冬，且二者間CG位點的DNA甲基化程度顯著不同。

可見，探明忍冬植株對干旱脅迫的應(yīng)答機(jī)制，對于穩(wěn)定金銀花藥材質(zhì)量、指導(dǎo)金銀花生產(chǎn)、發(fā)展金銀花生態(tài)種植具有重要意義。然而，相關(guān)研究尚未見報道。本研究利用重亞硫酸鹽測序法測定在甘露醇模擬干旱脅迫環(huán)境下忍冬葉片的DNA 甲基化水平，并分析其甲基化模式，闡述干旱脅迫對于忍冬全基因組DNA甲基化水平變化的影響，同時結(jié)合轉(zhuǎn)錄組數(shù)據(jù)分析干旱脅迫下忍冬DNA甲基化差異區(qū)域與表達(dá)水平差異基因之間的調(diào)控關(guān)系。研究結(jié)果將為今后忍冬新型植物分子抗旱育種及干旱脅迫下忍冬表觀遺傳變化的分子機(jī)制奠定基礎(chǔ)。

1 材料與儀器

1.1 材料

樣品采自經(jīng)山東中醫(yī)藥大學(xué)藥草園，由山東中醫(yī)藥大學(xué)蒲高斌教授鑒定為忍冬科忍冬屬植物忍冬Thunb.

1.2 儀器

CFX96 Touch Real-Time PCR儀（Bid-Rad 公司）；5424D型高速離心機(jī)（Eppendorf 公司）；甘露醇購自上海源葉生物科技有限公司；Hoagland’s 營養(yǎng)液購自青島海博生物技術(shù)有限公司；多糖多酚植物RNA提取試劑盒購自南京諾唯贊生物科技股份有限公司；SYBR Preminx EX TaqTM（ROX）、RT 反轉(zhuǎn)錄試劑盒購自TaKaRa公司。

2 方法

2.1 樣品處理

全基因組甲基化及轉(zhuǎn)錄組測定材料為3年生“華金2號”忍冬扦插苗。將園土與基質(zhì)按1∶1比例混合，基質(zhì)含草炭土、木質(zhì)泥炭、優(yōu)等椰殼粉、蛭石、高嶺石及有機(jī)物質(zhì)和生長所需的營養(yǎng)物質(zhì)，裝入內(nèi)徑為38 cm的花盆中。選擇生長情況良好、大小近似的忍冬苗移入花盆內(nèi)，每盆1株，植物在溫室條件中生長，進(jìn)行常規(guī)水分管理。向每盆扦插苗中澆灌2 L 0.2 mol/L的甘露醇用來進(jìn)行模擬干旱處理，在處理的0（CK）、3（DS-3d）、6（DS-6d）、9 d（DS-9d）分別取植株基本一致位置的幼葉，并立即用液氮冷凍，放入?80 ℃冰箱中備用。

MYB轉(zhuǎn)錄因子表達(dá)模式分析采用“華金2號”忍冬幼苗。培養(yǎng)條件為光照培養(yǎng)16 h、暗培養(yǎng)8 h，溫度為25 ℃。待長至第6片葉時，將材料分成2組，一組設(shè)置為對照，一組在培養(yǎng)液中加入甘露醇使其終濃度為0.2 mol/L，分別選取對照及干旱脅迫處理1、3、8、12、24 h植物材料的葉片，液氮冷凍后放入?80 ℃冰箱中保存?zhèn)溆谩?/p>

2.2 DNA的提取、質(zhì)檢、建庫及測序

利用植物DNA提取試劑盒提取不同干旱脅迫下樣品DNA，DNA樣品經(jīng)廣州基迪奧生物科技有限公司質(zhì)檢合格后構(gòu)建亞硫酸鹽（BS-seq）文庫構(gòu)建，測序步驟為（1）超聲處理DNA破碎為100～300 bp的片段。（2）DNA片段末端修復(fù)，并在3’端加上A堿基，然后將甲基化測序接頭連接到基因組片段上。（3）采用ZYMO EZ DNA Methylation- Gold kit進(jìn)行Bisulfite處理，Bisulfite處理會使未甲基化的胞嘧啶轉(zhuǎn)化為尿嘧啶。（4）脫鹽處理后切膠回收并進(jìn)行文庫片段大小選擇，PCR擴(kuò)增后重復(fù)大小選擇；構(gòu)建完成后，檢測質(zhì)量。最后使用Illumina Hi-SeqTM2500平臺測序。

2.3 測序數(shù)據(jù)過濾及比對分析

測序后得到原始數(shù)據(jù)，從原始數(shù)據(jù)去除含N比例大于10%、低質(zhì)量的reads后得到高質(zhì)量干凈讀數(shù)（clean data）。使用BSMAP（version：2.90）軟件將clean data比對到參考基因組序列上，選取唯一比對到的序列進(jìn)行下一步分析。參考基因組選用忍冬基因組，源自國家基因組數(shù)據(jù)中心（National Genomics Data Center，https://ngdc.cncb.ac.cn/ gwh/），登錄號為GWHAAZE00000000[14]。

2.4 甲基化水平及差異甲基化區(qū)域分析

對不同處理組之間進(jìn)行甲基化分析，主要是統(tǒng)計不同處理下樣品全基因組的甲基化水平和特定元件甲基化水平。使用methylkit軟件（version：1.7.10）進(jìn)行DNA差異甲基化區(qū)域（differentially methylated regions，DMR）分析，采用200 bp窗口在全基因組掃描，統(tǒng)計各個窗口內(nèi)特殊位點的平均DNA甲基化率，比較差異。

大學(xué)英語教學(xué)有其自身的特點：從教學(xué)內(nèi)容來看，大學(xué)英語教學(xué)內(nèi)容具有前沿性和實踐性；從教學(xué)管理來看，大學(xué)的管理模式與中學(xué)有很大的差異，學(xué)生的自由支配時間很多；從大學(xué)生學(xué)習(xí)過程的特點來看，學(xué)生可以根據(jù)自身未來發(fā)展的需要，制定學(xué)習(xí)計劃和目標(biāo)，在學(xué)習(xí)的過程中不斷鉆研，充分發(fā)揮學(xué)習(xí)的自主性，以彌補(bǔ)課堂教學(xué)中欠缺的在認(rèn)知風(fēng)格、學(xué)習(xí)策略、學(xué)習(xí)動機(jī)以及能力方面的差異，以實現(xiàn)自我職業(yè)發(fā)展的需求。

2.5 差異甲基化區(qū)域相關(guān)基因功能富集分析

得到全基因組DNA甲基化數(shù)據(jù)后，可以通過對DMR相關(guān)基因進(jìn)行基因本體（gene ontology，GO）和基因組百科全書（Kyoto encyclopedia of genes and genomes，KEGG）通路富集分析，了解顯著富集的生物學(xué)功能，從而推斷出干旱脅迫對忍冬基因表達(dá)的影響[15]。

2.6 轉(zhuǎn)錄組測序

轉(zhuǎn)錄組數(shù)據(jù)為本課題組前期實驗所得。利用植物RNA提取試劑盒提取樣品RNA，RNA樣本檢測合格后進(jìn)行普通轉(zhuǎn)錄組文庫構(gòu)建，文庫檢測合格后，利用Illumina NovaSeq6000進(jìn)行測序。對測序結(jié)果進(jìn)行過濾，去除掉低質(zhì)量讀數(shù)后得到高質(zhì)量的clean data，將得到的clean data與參考基因組進(jìn)行比對，參考基因組同甲基化測序參考基因組。利用featureCounts（1.5.0-p3）計算映射到每個基因的讀數(shù)，然后根據(jù)基因的長度計算每個基因的FPKM，并計算映射到該基因的讀數(shù)，以此來定量表示基因表達(dá)水平。基于基因表達(dá)定量分析得到數(shù)據(jù)，采用edgeR進(jìn)行基因差異表達(dá)顯著性分析，根據(jù)log2|(fold change)|＞1且值＜0.05，得到不同干旱脅迫處理下與對照組相比的差異表達(dá)基因（differentially expressed gene，DEG）。

2.7 全基因組DNA甲基化聯(lián)合轉(zhuǎn)錄組分析

根據(jù)不同處理組全基因組DNA甲基化數(shù)據(jù)中meth.diff值將差異甲基化基因分為甲基化水平上升組（A1）和甲基化水平下降組（A2）。根據(jù)不同處理組轉(zhuǎn)錄組數(shù)據(jù)中基因的FPKM值和log2FC值將差異表達(dá)基因分為表達(dá)水平上升組（B1）和表達(dá)水平下降組（B2），綜合分析數(shù)據(jù)。

2.8 干旱脅迫下忍冬MYB（L. japonica MYB，LjMYB）轉(zhuǎn)錄因子表達(dá)模式分析

從“2.7”項結(jié)果中，篩選出甲基化水平下降同時轉(zhuǎn)錄組水平上升的轉(zhuǎn)錄因子，利用qRT-PCR驗證干旱脅迫下MYB轉(zhuǎn)錄因子的表達(dá)情況，基因序列在網(wǎng)站[Home-Genome Warehouse（www.cncb.ac.cn）]中查找，利用Primer 3設(shè)計Real-Time PCR引物，內(nèi)參基因為金銀花Lj18S，所用引物鉑尚生物技術(shù)有限公司負(fù)責(zé)合成。將對照以及干旱脅迫處理1、3、8、12、24 h葉片加液氮研磨成粉末，利用多糖多酚植物RNA提取試劑盒（南京諾唯贊）提取RNA，采用反轉(zhuǎn)錄試劑盒（Takara）將RNA反轉(zhuǎn)錄為cDNA。以反轉(zhuǎn)錄后的cDNA為模板，利用SYBR試劑盒（Takara）進(jìn)行qRT-PCR。反應(yīng)體系25 μL：TB Green Premix Ex TaqⅡ（2×）12.5 μL，正向引物1 μL，反向引物1μL，cDNA模板2 μL，滅菌水8.5 μL。反應(yīng)程序為：95 ℃預(yù)變性3 min；95 ℃變性10 s，54 ℃退火30 s、72 ℃延伸1 min，39個循環(huán)，65～95 ℃過程中每隔0.5 ℃作溶解曲線，設(shè)3個重復(fù)。采用2?ΔΔCt法計算基因的相關(guān)表達(dá)水平。

3 結(jié)果與分析

3.1 亞硫酸氫鹽測序總結(jié)

采用全基因組重亞硫酸氫鹽甲基化測序測定忍冬在不同干旱脅迫處理下DNA甲基化水平，利用Illumina Hi-Seq 2500測序平臺測序，測序數(shù)據(jù)處理及與參考基因組的比對結(jié)果見表1。在不同處理時間下，過濾測序得到的原始數(shù)據(jù)，分別得到304 537 910、273 790 354、185 030 316、207 805 832個clean reads。使用BSMAP比對數(shù)據(jù)與基因組序列，比對率分別為83.60%、83.31%、83.66%、84.39%，表明該方法和結(jié)果具有較高的可靠性和準(zhǔn)確性。此外，基因組的有效覆蓋率在70%以上，不同甲基化類型之間甲基化水平存在明顯差異如圖1所示。

表1 重亞硫酸氫鹽測序結(jié)果

Table 1 Results of bisulfite sequencing

樣品名總序列讀數(shù)比對序列讀數(shù)比對率/%亞硫酸鹽處理轉(zhuǎn)化率/%比對深度覆蓋率/% CK304 537 910254 578 51183.6099.4942.2573.50 DS-3d273 790 354228 083 36583.3199.5037.8572.49 DS-6d185 030 316154 798 26583.6699.5825.6971.07 DS-9d207 805 832175 362 37684.3999.5229.10 72.67

圖1 不同序列環(huán)境甲基化占比

3.2 忍冬DNA甲基化的全基因組模式

3.2.1 忍冬全基因組DNA甲基化分析 全基因組平均甲基化水平由有效覆蓋C堿基的甲基化水平之和與有效覆蓋堿基總數(shù)之比決定，甲基化胞嘧啶的百分比根據(jù)局部序列環(huán)境和外部處理而變化。在CK、DS-3d、DS-6d、DS-9d 4種環(huán)境下，忍冬全基因組mC水平分別為23.60%、26.97%、28.37%、25.06%。與CK相比在干旱處理3、6 d時，全基因組甲基化水平會隨著干旱時間的延長不斷上升。9 d時，雖然甲基化水平與3、6 d相比有所下降，但仍然高于CK組。

CK組處理的基因組中mCG、mCHG、mCHH的甲基化水平分別為79.84%、50.95%、11.42%（表2），這反映了忍冬基因組中甲基化水平的百分比。同時，通過分析mCG、mCHG和mCHH序列相對于總mC位點出現(xiàn)的比例，可以反映出mC位點在3種序列環(huán)境中的分布。如圖1所示，甲基化胞嘧啶最常出現(xiàn)在mCHH位點（60.01%），在mCG、mCHG序列出現(xiàn)頻率較低，分別為23.64%和16.31%。

表2 甲基化水平統(tǒng)計

Table 2 Statistical table of methylation level

樣品名C/%CG/%CHG/% CHH/% CK23.6079.8450.9511.42 DS-3d26.9780.3852.4615.46 DS-6d28.3781.3954.2016.56 DS-9d25.0680.5351.7512.79

對于每一區(qū)域來說，在3種不同序列環(huán)境中CG位點的平均甲基化水最高，而CHH序列環(huán)境則呈現(xiàn)出最低的甲基化水平。具體來說，在Up2k區(qū)域發(fā)生甲基化水平較高即啟動子區(qū)域甲基化水平較高，其次是內(nèi)含子和外顯子區(qū)域，在CG序列環(huán)境下，內(nèi)含子區(qū)域CG平均甲基化水平最高，其次為外顯子區(qū)域與啟動子區(qū)域。在CHG序列環(huán)境下，啟動子區(qū)域CHG平均甲基化水平最高，遠(yuǎn)高于外顯子和內(nèi)含子區(qū)域。此外，CHH序列環(huán)境下各基因組區(qū)域的DNA甲基化水平模式與CHG序列環(huán)境基本類似。干旱脅迫后，各基因組區(qū)域DNA甲基化水平模式基本未變，但在甲基化率方面有所變化，啟動子區(qū)域中，CHH類型甲基化率有所上升；外顯子區(qū)域甲基化率基本未變，內(nèi)含子區(qū)域各甲基化類型甲基化水平均有小幅上升（圖2）。

3.3 胞嘧啶DNA甲基化序列偏好性分析

在真核生物中，全基因組范圍內(nèi)甲基化位點附近堿基的序列特征，對了解甲基化發(fā)生的序列偏好有重要作用。為了揭示序列環(huán)境和忍冬甲基化偏好之間是否存在關(guān)系，使用Logo Plots的工具對CHG堿基進(jìn)行統(tǒng)計，得到全基因組中CHG與mCHG堿基的附近9 bp的堿基構(gòu)成，通過CHG到mCHG的比較，探索甲基化位點或其附近區(qū)域的序列信息。在所有甲基化的mC位點中，mC位點周圍的序列偏好隨著CG背景和非CG背景而變化（圖3），在對稱CG環(huán)境中，mC經(jīng)常出現(xiàn)在TCGA序列中，而mCHG背景下mC位點的甲基化最經(jīng)常出現(xiàn)在CAG序列中，而mCHH范圍內(nèi)的甲基化最常發(fā)生在CAA序列中。因此，推斷mC位點的大多數(shù)甲基化發(fā)生在CAN（N代表A、T、C、G）。

Up2k為基因上游2 kb，Down2k為基因下游2 kb，5’UTR為5’非翻譯區(qū)，3’UTR為3’非翻譯區(qū)，Exon為外顯子，CDS為編碼序列，Intron為內(nèi)含子

3.4 忍冬基因組的DMR及相關(guān)基因富集分析

3.4.1 DNA甲基化的DMR分析 將不同干旱處理組與對照組進(jìn)行比較，分別檢測出5926、6978和5172個DMR（圖4）。對其進(jìn)行分析發(fā)現(xiàn)，在DS-3d組中，有3953個DMR甲基化水平下降，有1973個DMR甲基化水平上升。在DS-6d組中，有2068個DMR甲基化水平下降，有4910個DMR甲基化水平上升。在DS-9d組中，有3224個DMR甲基化水平下降，有1948個DMR甲基化水平上升。與對照組相比，處理組甲基化水平下降DMR數(shù)量均大于甲基化水平上升DMR數(shù)量。

3.4.2 DMR相關(guān)基因富集分析 將前期得到的DMR相關(guān)基因進(jìn)行GO和KEGG富集分析，挖掘出與忍冬在干旱脅迫條件下化學(xué)成分變化相關(guān)的生物學(xué)過程和pathway調(diào)控通路，進(jìn)而探討差異甲基化基因在忍冬響應(yīng)干旱脅迫中的調(diào)控作用。

GO富集分析結(jié)果顯示在DS-3d、DS-6d、DS-9d處理下，分別有478、557和464個顯著富集的GO條目（＜0.05），分別來自于193、2158和1750個DMR相關(guān)基因，稱為DMG（差異性甲基化相關(guān)基因）。將3組比較組的差異甲基化基因富集到GO的3類中，其二級注釋如圖5所示，在細(xì)胞組成項下，DMG主要富集到細(xì)胞，細(xì)胞組分和細(xì)胞器。在生物過程項下，DMG主要富集到細(xì)胞過程，代謝過程和單一生物過程。在分子功能項下，DMG主要富集到整合，催化反應(yīng)和運(yùn)輸反應(yīng)，表明甲基化介導(dǎo)忍冬中這些類型的途徑。

X軸表示甲基化胞嘧啶位于第4位的堿基位置，Y軸表示堿基的熵

1-DS-3d組與CK相比 2-DS-6d組與CK相比 3-DS-9d組與CK相比

KEGG富集分析顯示（圖6），在DS-3d、DS-6d、DS-9d處理下分別有7、15和12條通路顯著富集（＜0.05）。DS-3d組與CK相比，中顯著富集的pathway包含亞油酸代謝、甘油酯代謝、組氨酸代謝、萜類生物合成、苯丙氨酸、酪氨酸和色氨酸生物合成以及次生代謝物的生物合成等7條通路，涉及到387個相關(guān)基因。DS-6d組與CK相比，顯著富集的pathway包括次生代謝物的生物合成、谷胱甘肽代謝、脂肪酸代謝、異喹啉生物堿生物合成、淀粉與蔗糖的代謝、脂肪酸生物合成以及ABC轉(zhuǎn)運(yùn)黃酮類生物合成等15條通路，涉及468個相關(guān)基因。DS-9d組與CK相比，顯著富集的途徑包括氮代謝、酪氨酸代謝、異喹啉生物堿生物合成、亞油酸代謝、過氧化物酶、二苯乙烯、二芳基庚烷和姜酚的生物合成和次生代謝物的生物合成等12條通路，涉及393個相關(guān)基因，這些途徑可能與忍冬對干旱的響應(yīng)密切相關(guān)。

3.5 轉(zhuǎn)錄組數(shù)據(jù)分析

對照組及不同干旱脅迫處理組測序結(jié)果經(jīng)過濾后，分別得到44 402 319、44 655 626、43 873 937、42 504 058個clean reads，與參考基因組的比對率分別為93.28%、91.98%、92.36%、92.43%。通過與對照組相比，在DS-3d、DS-6d、DS-9d分別鑒定出8531、4171和6712個差異表達(dá)基因。其中在DS-3d時，有4126個基因表達(dá)上升，4405個基因表達(dá)下降。在DS-6d時，有2849個基因表達(dá)上升，1322個基因表達(dá)下降。在DS-9d時，有3727個基因表達(dá)上升，2985個基因表達(dá)下降。

3.6 全基因組DNA甲基化與轉(zhuǎn)錄組聯(lián)合分析

為了探討干旱脅迫期間觀察到的mC甲基化變化是否會改變基因表達(dá)，綜合分析差異甲基化基因和差異表達(dá)基因的數(shù)據(jù)。在DS-3d時有82個高表達(dá)基因其甲基化水平升高以及264個高表達(dá)基因其甲基化水平降低；在DS-6d時有56個高表達(dá)基因其甲基化水平升高以及341個高表達(dá)基因其甲基化水平降低（圖7）；在DS-9d時有96個高表達(dá)基因其甲基化水平升高以及237個高表達(dá)基因其甲基化水平降低；前期KEGG富集分析顯示，這些基因廣泛參與次級代謝產(chǎn)物的生物合成、苯丙烷代謝、類黃酮生物合成、醌和其他萜醌生物合成、異喹啉生物堿生物合成、MAPK信號通路-植物、激素信號轉(zhuǎn)導(dǎo)、倍半萜類和三萜類生物合成、萜類骨架生物合成以及二萜類生物合成等，由此推測干旱脅迫下忍冬可通過體內(nèi)DNA甲基化水平的改變調(diào)控次生代謝產(chǎn)物合成途徑中相關(guān)基因的表達(dá)。

A-DS-3d組與CK相比 B-DS-6d組與CK相比 C-DS-9d組與CK相比；紅色代表生物過程類，綠色代表細(xì)胞組成類，藍(lán)色代表分子功能類

A-DS-3d組與CK相比 B-DS-6d組與CK相比 C-DS-9d組與CK相比

A1-甲基化水平上升DMR A2-甲基化水平下降DMR B1-表達(dá)水平上升DEG B2-表達(dá)水平下降DEG

3.7 干旱脅迫下LjMYB轉(zhuǎn)錄因子表達(dá)模式分析

在DS-3d、DS-6d、DS-9d 3個處理組轉(zhuǎn)錄組水平上升且甲基化水平下降的基因中共篩選出10個MYB轉(zhuǎn)錄因子，分別命名為～。利用qRT-PCR對10個轉(zhuǎn)錄因子在對照以及0.2 mol/L甘露醇處理不同時間下的表達(dá)模式進(jìn)行分析。

qRT-PCR結(jié)果如圖8所示，干旱脅迫下轉(zhuǎn)錄因子的表達(dá)各不相同，基因表達(dá)水平整體呈上升趨勢。其中和基因表達(dá)水平在脅迫處理3 h時差異顯著，分別升高13.19倍和7.73倍，和在脅迫處理24 h時差異顯著，分別升高9.13倍和7.64倍。在干旱脅迫8 h時，基因表達(dá)水平達(dá)到最高，而、在脅迫處理24 h表達(dá)量最高。、、在不同脅迫時間下表達(dá)量無顯著差異。結(jié)果表明在干旱脅迫處理不同時間下轉(zhuǎn)錄因子的變化有所不同，是一個動態(tài)的變化，推測在干旱脅迫下植物可通過改變MYB轉(zhuǎn)錄因子的甲基化水平來調(diào)控體內(nèi)成分變化以應(yīng)對脅迫。

4 討論

全基因組重亞硫酸鹽甲基化測序法（whole- genome bisulfite sequencing，WGBS）技術(shù)作為基因組甲基化測定的一種方法，因具有精確度高、結(jié)果可靠及單堿基分辨率等優(yōu)點成為DNA甲基化檢測的“金標(biāo)準(zhǔn)”[16]。目前在已經(jīng)進(jìn)行基因組DNA甲基化分析的34種植物中，mCG在全基因組范圍內(nèi)具有最高水平的DNA甲基化，擬南芥mCG在已測物種中最低，為30.5%，甜菜 mCG最高，為92.5%。CHG位點甲基化水平變化范圍從9.3%的鹽水水芹到81.2%的甜菜。葡萄中mCHH含量最低，僅為1.1%，甜菜中mCHH含量最高，為18.8%[17-18]。本研究利用WGBS技術(shù)測定忍冬在不同干旱脅迫條件下DNA甲基化水平，CG位點甲基化水平從79.84%～81.39%，CHG位點甲基化水平從50.95%～54.20%，CHH位點甲基化水平從11.42%～16.56%。與其他物種相比，忍冬DNA甲基化水平在變化范圍之內(nèi)，且忍冬基因組各序列環(huán)境具有較高的甲基化水平，甲基化水平隨著干旱脅迫條件的改變而變化。

*P＜0.05，**P＜0.01，***P＜0.001

有研究表明DNA甲基化在植物干旱脅迫響應(yīng)中具有重要作用，Wang等[19]研究結(jié)果表明干旱誘導(dǎo)的DNA甲基化變化對水稻植株的干旱脅迫響應(yīng)有相當(dāng)大的影響，干旱誘導(dǎo)的DMRs在DK151植物中檢測到更多，即在干旱條件下，甲基組在耐旱基因型中比在干旱敏感基因型中更穩(wěn)定；Li等[20]采用WGBS法測定桑樹在干旱脅迫下的甲基化變化發(fā)現(xiàn)，干旱脅迫后桑樹甲基化水平普遍高于對照。本研究發(fā)現(xiàn)在干旱脅迫處理后，忍冬基因組甲基化水平普遍上調(diào)，且隨著干旱脅迫時間的延長，甲基化水平呈上升趨勢。進(jìn)一步分析發(fā)現(xiàn)，在整體甲基化水平上升的同時，仍有部分基因甲基化水平下降，結(jié)合轉(zhuǎn)錄組數(shù)據(jù)進(jìn)行分析，發(fā)現(xiàn)甲基化水平下降的差異甲基化基因其基因表達(dá)水平大部分也發(fā)生變化。富集分析顯示這些基因主要富集于次生代謝產(chǎn)物合成相關(guān)通路，與忍冬中類黃酮類、萜類成分的合成密切相關(guān)。由此推測DNA甲基化調(diào)控基因表達(dá)是干旱影響金銀花有效成分生物合成的作用機(jī)制之一。近年來，越來越多的研究表明,合成和積累次生代謝產(chǎn)物是藥用植物最重要的防御環(huán)境脅迫的策略，環(huán)境脅迫導(dǎo)致藥用植物次生代謝產(chǎn)物的積累，次生代謝產(chǎn)物通常也是藥用植物的主要藥用成分[21-22]，通過本研究將為進(jìn)一步探究環(huán)境影響金銀花質(zhì)量的分子機(jī)制奠定基礎(chǔ)。

MYB轉(zhuǎn)錄因子是目前研究較多的轉(zhuǎn)錄因子家族，廣泛參與調(diào)控植物的生長發(fā)育、抗逆脅迫、次級代謝產(chǎn)物的積累[23]。研究發(fā)現(xiàn)MYB轉(zhuǎn)錄因子的甲基化水平影響功能基因的表達(dá)進(jìn)而對活性成分產(chǎn)生影響。通過基因組重測序結(jié)合轉(zhuǎn)錄組分析定量性狀定位，發(fā)現(xiàn)紅肉蘿卜L. RsMYB1啟動子區(qū)域的DNA高甲基化抑制基因表達(dá)，使花青素的生物合成收到抑制，導(dǎo)致了白色果肉表型[24]。Li等[25]研究發(fā)現(xiàn)MYB轉(zhuǎn)錄因子MdLUX和MdPCL啟動子區(qū)mCG環(huán)境的甲基化水平與其mRNA水平和花青素積累呈負(fù)相關(guān)，MdLUX和MdPCL-like通過啟動子區(qū)DNA低甲基化和結(jié)構(gòu)類黃酮基因的激活促進(jìn)蘋果果皮花青素生物合成。Tang等[26]對異色菊花（YP）無性繁殖產(chǎn)生的黃花（YP-Y）和粉色花（YP-P）展開研究發(fā)現(xiàn)，CmMYB6啟動子的甲基化水平是花色產(chǎn)生變異的關(guān)鍵，其甲基化程度影響著基因的表達(dá)和花青素的生物合成，對CmMYB6啟動子進(jìn)行去甲基化處理后，花朵顏色由黃色恢復(fù)為粉紅色。這些研究顯示MYB轉(zhuǎn)錄因子的甲基化水平對于植物體內(nèi)類黃酮成分的生物合成有著重要作用。

類黃酮類成分是忍冬中的重要活性成分，是忍冬次生代謝產(chǎn)物的一種，前期有研究表明，適當(dāng)?shù)母珊得{迫會使類黃酮類成分含量增加[27]，本研究通過組學(xué)分析加熒光定量PCR分析初步篩選出4個在干旱脅迫下表達(dá)量差異顯著的轉(zhuǎn)錄因子，并通過通路富集分析發(fā)現(xiàn)其與忍冬的次生代謝產(chǎn)物的生物合成相關(guān)，下一步將繼續(xù)深入探究這4個轉(zhuǎn)錄因子啟動子區(qū)域甲基化的變化以及不同干旱脅迫處理下忍冬中活性成分的變化，為進(jìn)一步解析忍冬在干旱脅迫下類黃酮類成分變化分子機(jī)制提供基礎(chǔ)。

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

[1] 中國藥典[S]. 一部. 2020: 230.

[2] 黃璐琦, 郭蘭萍, 胡娟, 等. 道地藥材形成的分子機(jī)制及其遺傳基礎(chǔ) [J]. 中國中藥雜志, 2008, 33(20): 2303-2308.

[3] 師鳳華, 趙祥升, 莫喬程, 等. 不同產(chǎn)地不同種質(zhì)金銀花中綠原酸和木犀草苷的含量分析 [J]. 中藥材, 2014, 37(7): 1145-1148.

[4] 徐迎春, 張佳寶, 蔣其鰲, 等. 水分脅迫對忍冬生長及金銀花質(zhì)量的影響 [J]. 中藥材, 2006, 29(5): 420-423.

[5] 孫希芳, 王斌, 張永清, 等. 光照對金銀花藥材化學(xué)成分的影響 [J]. 山東農(nóng)業(yè)科學(xué), 2014, 46(10): 63-67.

[6] 王玲娜, 張金, 劉紅燕, 等. 金銀花質(zhì)量影響因素分析 [J]. 中國醫(yī)藥導(dǎo)報, 2014, 11(25): 159-161.

[7] 張萍, 蒲高斌. 氣候因子對金銀花綠原酸含量影響的研究 [J]. 山東農(nóng)業(yè)科學(xué), 2015, 47(9): 77-79.

[8] 呂芳, 蘇幼紅, 張富春, 等. 植物活性成分對表觀遺傳調(diào)節(jié)的研究概況 [J]. 中草藥, 2008, 39(10): 1580-1583.

[9] 薛梅, 陳成彬, 陳力, 等. 半夏多倍體復(fù)合體基因組DNA甲基化狀態(tài)的MSAP分析 [J]. 中草藥, 2008, 39(11): 1713-1716.

[10] Kuo T C, Chen C H, Chen S H,. The effect of red light and far-red light conditions on secondary metabolism in agarwood [J]., 2015, 15: 139.

[11] 謝宜杰, 劉曉瑩, 羅可可, 等. 醇型和酮型廣藿香DNA甲基化的MSAP分析 [J]. 中草藥, 2020, 51(20): 5293-5301.

[12] Li L Q, Li X L, Fu C H,. Sustainable use of Taxus media cell cultures through minimal growth conservation and manipulation of genome methylation [J]., 2013, 48(3): 525-531.

[13] Zha L P, Liu S, Liu J,. DNA methylation influences chlorogenic acid biosynthesis in Lonicera japonica by mediating LjbZIP8 to regulate phenylalanine ammonia-lyase 2 expression [J]., 2017, 8: 1178.

[14] Pu X D, Li Z, Tian Y,. The honeysuckle genome provides insight into the molecular mechanism of carotenoid metabolism underlying dynamic flower coloration [J]., 2020, 227(3): 930-943.

[15] 王瀟, 尹天舒, 李柏逸, 等. 基因功能富集分析的研究進(jìn)展 [J]. 中國科學(xué): 生命科學(xué), 2016, 46(4): 363-373.

[16] 邢朝斌, 張妍彤, 王卓, 等. DNA甲基化及重亞硫酸鹽測序法在藥用植物中的應(yīng)用策略 [J]. 中草藥, 2017, 48(24): 5063-5069.

[17] Niederhuth C E, Bewick A J, Ji L X,. Widespread natural variation of DNA methylation within angiosperms [J]., 2016, 17(1): 194.

[18] 申曉慧. 兩種誘變處理下紫花苜蓿全基因組DNA甲基化圖譜分析 [J]. 草地學(xué)報, 2022, 30(1): 46-54.

[19] Wang W S, Qin Q, Sun F,. Genome-wide differences in DNA methylation changes in two contrasting rice genotypes in response to drought conditions [J]., 2016, 7: 1675.

[20] Li R X, Hu F, Li B,. Whole genome bisulfite sequencing methylome analysis of mulberry () reveals epigenome modifications in response to drought stress [J]., 2020, 10(1): 8013.

[21] 蘇琦, 張珂, 黃淑琪, 等. 逆境脅迫對藥用植物藥效成分積累的影響 [J]. 湖北民族大學(xué)學(xué)報: 自然科學(xué)版, 2022, 40(2): 129-134.

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

[23] 王曉童, 向麗, 鄔蘭, 等. 黃芩MYB轉(zhuǎn)錄因子家族全基因組鑒定與分析 [J]. 世界中醫(yī)藥, 2022, 17(13): 1802-1812.

[24] Wang Q B, Wang Y P, Sun H H,. Transposon-induced methylation of the RsMYB1 promoter disturbs anthocyanin accumulation in red-fleshed radish [J]., 2020, 71(9): 2537-2550.

[25] Li W F, Ning G X, Zuo C W,. MYB_ SH[AL]QKY[RF]transcription factors MdLUX and MdPCL-like promote anthocyanin accumulation through DNA hypomethylation and MdF3H activation in apple [J]., 2021, 41(5): 836-848.

[26] Tang M W, Xue W J, Li X Q,. Mitotically heritable epigenetic modifications of CmMYB6 control anthocyanin biosynthesis in chrysanthemum [J]., 2022, 236(3): 1075-1088.

[27] 岳凱, 劉文瑜, 魏小紅. 干旱脅迫對不同品系藜麥內(nèi)黃酮和抗氧化性的影響 [J]. 分子植物育種, 2019, 17(3): 956-962.

Whole genome DNA methylation and transcriptome analysis ofin response to drought stress

XU Xiao-han1, TANG Zhi-qiang1, LIU Qian1, 2, LI Jia1, 2, LIU Zhen-hua1, 2, ZHANG Yong-qing1, 2, PU Gao-bin1, 2

1. Shandong University of Traditional Chinese Medicine, Jinan 250355, China 2. Shandong Provincial Collaborative Innovation Center for Quality Control and Construction of the Whole Industrial Chain of Traditional Chinese Medicine, Jinan 250355, China

To analyze the changes and patterns of genome-wide DNA methylation level and its relationship with gene expression inunder drought stress, and lay a foundation for further exploring the molecular mechanism ofin response to drought stress.The whole-genome bisulfite sequencing (WGBS) was used to determine the DNA methylation level ofunder drought stress, and the differentially expressed genes related to differentially methylated genes were analyzed in combination with transcriptome sequencing results.The overall methylation level ofunder drought stress generally increased; The methylation level of CG type was significantly higher than that of CHG or CHH type. Differentially methylated regions (DMR) analysis ofunder different drought stress treatments revealed 1935, 2158 and 1750 differentially methylated genes. Kyoto encyclopedia of genes and genomes (KEGG) enrichment analysis showed that the significantly enriched pathways included linoleic acid metabolism, nitrogen metabolism and biosynthesis of secondary metabolites. Through transcriptome sequencing technology, it was found that differentially expressed genes related to differentially methylated genes were mainly enriched in secondary metabolite synthesis-related pathways, and fourtranscription factors with significant increase in gene expression under drought stress were screened.Under drought stress, the methylation level ofincreased, suggesting that DNA methylation regulation of gene expression is one of the mechanisms of drought affecting the biosynthesis of effective components of.

Thunb.; drought stress; DNA methylation; whole genome bisulfite sequencing; transcriptome sequencing

R286.12

A

0253 - 2670(2023)16 - 5339 - 11

10.7501/j.issn.0253-2670.2023.16.023

2023-02-03

山東省重點研發(fā)計劃（鄉(xiāng)村振興科技創(chuàng)新提振行動計劃）項目（2022TZXD0036）；山東省農(nóng)業(yè)良種工程項目（2021LZGC008）；山東省現(xiàn)代農(nóng)業(yè)產(chǎn)業(yè)技術(shù)體系中草藥創(chuàng)新團(tuán)隊項目（SDAIT-20）

許小涵（1999—），女，碩士研究生，研究方向為中藥分子生物學(xué)。Tel: 15194181753 E-mail: xuxiaohan0719@163.com

蒲高斌（1979—），男，博士，教授，主要從事中藥資源與分子生藥學(xué)研究。E-mail: gbpu@163.com

[責(zé)任編輯 時圣明]