干旱脅迫下忍冬全基因組DNA甲基化和轉錄組分析

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

? 藥材與資源 ?

干旱脅迫下忍冬全基因組DNA甲基化和轉錄組分析

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

1. 山東中醫藥大學，山東 濟南 250355 2. 山東省中藥質量控制與全產業鏈建設協同創新中心，山東 濟南 250355

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

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

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

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

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

1 材料與儀器

1.1 材料

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

1.2 儀器

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

2 方法

2.1 樣品處理

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

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

2.2 DNA的提取、質檢、建庫及測序

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

2.3 測序數據過濾及比對分析

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

2.4 甲基化水平及差異甲基化區域分析

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

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

2.5 差異甲基化區域相關基因功能富集分析

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

2.6 轉錄組測序

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

2.7 全基因組DNA甲基化聯合轉錄組分析

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

2.8 干旱脅迫下忍冬MYB（L. japonica MYB，LjMYB）轉錄因子表達模式分析

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

3 結果與分析

3.1 亞硫酸氫鹽測序總結

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

表1 重亞硫酸氫鹽測序結果

Table 1 Results of bisulfite sequencing

樣品名總序列讀數比對序列讀數比對率/%亞硫酸鹽處理轉化率/%比對深度覆蓋率/% 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 不同序列環境甲基化占比

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

3.2.1 忍冬全基因組DNA甲基化分析 全基因組平均甲基化水平由有效覆蓋C堿基的甲基化水平之和與有效覆蓋堿基總數之比決定，甲基化胞嘧啶的百分比根據局部序列環境和外部處理而變化。在CK、DS-3d、DS-6d、DS-9d 4種環境下，忍冬全基因組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位點出現的比例，可以反映出mC位點在3種序列環境中的分布。如圖1所示，甲基化胞嘧啶最常出現在mCHH位點（60.01%），在mCG、mCHG序列出現頻率較低，分別為23.64%和16.31%。

表2 甲基化水平統計

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

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

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

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

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

3.4 忍冬基因組的DMR及相關基因富集分析

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

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

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

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個相關基因。DS-6d組與CK相比，顯著富集的pathway包括次生代謝物的生物合成、谷胱甘肽代謝、脂肪酸代謝、異喹啉生物堿生物合成、淀粉與蔗糖的代謝、脂肪酸生物合成以及ABC轉運黃酮類生物合成等15條通路，涉及468個相關基因。DS-9d組與CK相比，顯著富集的途徑包括氮代謝、酪氨酸代謝、異喹啉生物堿生物合成、亞油酸代謝、過氧化物酶、二苯乙烯、二芳基庚烷和姜酚的生物合成和次生代謝物的生物合成等12條通路，涉及393個相關基因，這些途徑可能與忍冬對干旱的響應密切相關。

3.5 轉錄組數據分析

對照組及不同干旱脅迫處理組測序結果經過濾后，分別得到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個差異表達基因。其中在DS-3d時，有4126個基因表達上升，4405個基因表達下降。在DS-6d時，有2849個基因表達上升，1322個基因表達下降。在DS-9d時，有3727個基因表達上升，2985個基因表達下降。

3.6 全基因組DNA甲基化與轉錄組聯合分析

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

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

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

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

3.7 干旱脅迫下LjMYB轉錄因子表達模式分析

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

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

4 討論

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

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

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

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

類黃酮類成分是忍冬中的重要活性成分，是忍冬次生代謝產物的一種，前期有研究表明，適當的干旱脅迫會使類黃酮類成分含量增加[27]，本研究通過組學分析加熒光定量PCR分析初步篩選出4個在干旱脅迫下表達量差異顯著的轉錄因子，并通過通路富集分析發現其與忍冬的次生代謝產物的生物合成相關，下一步將繼續深入探究這4個轉錄因子啟動子區域甲基化的變化以及不同干旱脅迫處理下忍冬中活性成分的變化，為進一步解析忍冬在干旱脅迫下類黃酮類成分變化分子機制提供基礎。

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

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

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

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

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

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

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

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

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

[9] 薛梅, 陳成彬, 陳力, 等. 半夏多倍體復合體基因組DNA甲基化狀態的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]. 中國科學: 生命科學, 2016, 46(4): 363-373.

[16] 邢朝斌, 張妍彤, 王卓, 等. DNA甲基化及重亞硫酸鹽測序法在藥用植物中的應用策略 [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]. 草地學報, 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]. 湖北民族大學學報: 自然科學版, 2022, 40(2): 129-134.

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

[23] 王曉童, 向麗, 鄔蘭, 等. 黃芩MYB轉錄因子家族全基因組鑒定與分析 [J]. 世界中醫藥, 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] 岳凱, 劉文瑜, 魏小紅. 干旱脅迫對不同品系藜麥內黃酮和抗氧化性的影響 [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

山東省重點研發計劃（鄉村振興科技創新提振行動計劃）項目（2022TZXD0036）；山東省農業良種工程項目（2021LZGC008）；山東省現代農業產業技術體系中草藥創新團隊項目（SDAIT-20）

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

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

[責任編輯 時圣明]