兩種誘變處理下紫花苜蓿全基因組DNA甲基化圖譜分析

申曉慧

(宜春學(xué)院， 生命科學(xué)與資源環(huán)境學(xué)院， 江西 宜春 336000)

紫花苜蓿(MedicagosativaL.)具有產(chǎn)量高、適口性好、利用年限長(zhǎng)、營(yíng)養(yǎng)價(jià)值豐富和廣適應(yīng)性等特點(diǎn)，隨著我國(guó)畜牧業(yè)的快速發(fā)展，家畜數(shù)量急劇上升，對(duì)飼草數(shù)量和品質(zhì)要求越來(lái)越高[1]，為此，紫花苜蓿在畜牧業(yè)發(fā)展中起著舉足輕重的作用。然而我國(guó)一些極端地區(qū)如干旱、高寒、鹽堿地等特殊的環(huán)境條件給紫花苜蓿豐產(chǎn)性帶來(lái)了嚴(yán)重影響，這就導(dǎo)致生產(chǎn)中迫切需要更新紫花苜蓿種質(zhì)資源來(lái)應(yīng)對(duì)這種特殊環(huán)境條件。因此篩選、培育抗性品種成為育種工作中的關(guān)鍵。實(shí)踐證明，利用人工誘變遺傳變異是豐富作物種質(zhì)資源，創(chuàng)制物種新材料及新品種選育的重要手段之一[2]。然而，誘變后代能否保證優(yōu)良性狀能夠穩(wěn)定遺傳對(duì)新品種選育至關(guān)重要。表觀遺傳可以通過(guò)調(diào)控基因表達(dá)進(jìn)而影響植物生長(zhǎng)發(fā)育，而DNA甲基化作為表觀遺傳學(xué)研究的一種，其在調(diào)節(jié)植物生長(zhǎng)發(fā)育、植物應(yīng)對(duì)生物或非生物脅迫環(huán)境以及維持基因組穩(wěn)定等方面發(fā)揮著至關(guān)重要的作用[3-5]。

目前，隨著基因二代測(cè)序技術(shù)的不斷向前發(fā)展，帶動(dòng)著全基因組的DNA甲基化的飛快進(jìn)步，其研究成果在多學(xué)科、多領(lǐng)域中都發(fā)揮著重要的積極作用。如薛倩[6]利用BS-seq技術(shù)構(gòu)建了京海黃雞腿肌全基因組DNA甲基化圖譜，揭示了其DNA甲基化特征。蘇暢[7]利用重亞硫酸鹽測(cè)序技術(shù)構(gòu)建了白樺基因組甲基化圖譜，綜合分析了DNA甲基化對(duì)于基因表達(dá)的作用。趙倩等[8]研究表明，水稻在高劑量重離子輻射下，DNA去甲基化的多態(tài)性水平明顯高于DNA甲基化的多態(tài)性水平。還有一些在小麥[9]、中草藥[10]以及綠化植物長(zhǎng)春花[11]等植物方面也都有諸多研究。在我國(guó)誘變技術(shù)在紫花苜蓿中已有廣泛應(yīng)用，如韓微波等[12]通過(guò)零磁空間處理肇東紫花苜蓿種子，從中選育出了新品種‘農(nóng)菁10號(hào)’。然而，在眾多研究中通過(guò)誘變手段導(dǎo)致紫花苜蓿全基因組DNA甲基化水平的變化研究卻鮮有報(bào)道。基于此，作者在已有誘變對(duì)苜蓿生長(zhǎng)、生理變化的影響[13-14]的研究基礎(chǔ)上，利用BS-seq技術(shù)進(jìn)一步深入分析誘變處理前后紫花苜蓿全基因組甲基化圖譜變化情況，為紫花苜蓿突變體篩選及誘變條件下苜蓿表觀遺傳變化的分子機(jī)制研究提供科學(xué)依據(jù)。

1 材料與方法

1.1 試驗(yàn)材料

選取紫花苜蓿‘龍牧806’，該品種由黑龍江省畜牧研究所提供。

1.2 試驗(yàn)方法

1.2.1試驗(yàn)處理 甲基磺酸乙酯(EMS)處理[15]：首先將成熟飽滿的苜蓿種子用濃硫酸處理5 min，之后用蒸餾水反復(fù)沖洗。然后在濃度為100 mmol·L-1,pH值 7.0的4℃磷酸緩沖液中浸泡12 h。最后用磷酸緩沖液配制成濃度為0.4%的EMS溶液，將浸泡過(guò)的種子置于室溫黑暗條件下處理15 h，而后用蒸餾水反復(fù)沖洗，去除殘留種子表面的EMS溶液。每處理保留100粒種子，3次重復(fù)。

60CO-γ射線處理[15]：將成熟飽滿的苜蓿種子進(jìn)行150 Gy劑量的60CO-γ射線輻照處理，處理周期為7 d。以未處理種子作為對(duì)照，每處理保留100粒種子，3次試驗(yàn)重復(fù)。

1.2.2DNA提取、建庫(kù)及測(cè)序 將上述各重復(fù)處理的種子分別進(jìn)行混合，每100粒作為一個(gè)樣本，每個(gè)樣品3個(gè)生物學(xué)重復(fù)，進(jìn)行DNA提取，DNA的提取按照張延召等[16]的方法進(jìn)行。DNA樣品經(jīng)武漢康測(cè)科技有限公司質(zhì)檢合格后進(jìn)行重亞硫酸鹽(BS-seq)文庫(kù)構(gòu)建，主要包括基因組DNA提取(CTAB法)，片段化，利用EZ DNA Methylation-Gold kit進(jìn)行Bisulfite處理及基于單鏈的DNA文庫(kù)構(gòu)建等4個(gè)步驟。文庫(kù)質(zhì)量檢測(cè)后通過(guò)Illumina XTen (Illumina，San Diego，CA，USA)測(cè)序儀進(jìn)行雙末端測(cè)序，具體方法參照文獻(xiàn)[17]。試驗(yàn)選取樣本材料分別為對(duì)照組(A34-1A)，EMS處理組(A34-2A)及60CO-γ射線處理組(A34-3A)。

表1 單鏈DNA濃度展示表Table 1 Single-stranded DNA concentration display table

1.2.3測(cè)序數(shù)據(jù)的過(guò)濾、比對(duì)分析 測(cè)序結(jié)果經(jīng)過(guò)識(shí)別堿基后轉(zhuǎn)化為原始測(cè)序序列(Raw reads)，然后除去Raw reads中含接頭序列、低質(zhì)量、含N比例>10%的Reads后得到高質(zhì)量序列(Clean Reads)。使用Bismark(v0.17.0)軟件對(duì)Clean reads進(jìn)行比對(duì)。通過(guò)將reads和參考基因組分別進(jìn)行C>T和G>A的轉(zhuǎn)換之后，分別進(jìn)行比對(duì)，篩選出比對(duì)結(jié)果最好的組合方式[18]。參考基因組選用Medicagotruncatula，源自ftp://ftp.ensemblgenomes.org/pub/release-33/plants/fasta/medicago_truncatula/dna/Medicago_truncatula.MedtrA17_4.0.dna.toplevel.fa.gz。

1.2.4苜蓿組織差異甲基化分析 使用swDMR軟件(version 1.0.0)檢測(cè)DMR[19]。swDMR采用滑動(dòng)窗口(sliding window size:1 000 bp;step length:100 bp)檢測(cè)DMR。大于2個(gè)樣本時(shí)，采用ANOVA方法。

1.2.5DMR(Differentially methylated region)基因KEGG通路富集分析 通路富集分析以KEGG Pathway為單位[20]，應(yīng)用超幾何檢驗(yàn)，找出與整個(gè)苜蓿基因組背景相比，在DMR相關(guān)基因中顯著富集的通路。采用KOBAS[21]對(duì)DMR所在基因列表進(jìn)行基因通路富集分析，篩選出通路富集P-value(或校正后P-value)小于0.05以及查詢基因在此通路中的基因數(shù)目大于2的通路。

2 結(jié)果與分析

2.1 紫花苜蓿全基因組DNA甲基化測(cè)序數(shù)據(jù)分析

使用Illumina XTen測(cè)序平臺(tái)對(duì)供試樣品進(jìn)行測(cè)序，共得到58.11 G的Clean reads (表2)。使用軟件SOAPnuke(version 1.6.0)將A34-1A，A34-2A 及A34-2A樣品得到的Clean reads比對(duì)到苜蓿基因組序列上，比對(duì)率分別為 90.40%，88.72%和90.60%。

表2 不同處理苜蓿基因測(cè)序數(shù)據(jù)統(tǒng)計(jì)及比對(duì)到參考基因組結(jié)果Table 2 Statistical results of sequencing data and clean reads map to reference genome results

2.2 紫花苜蓿全基因組DNA甲基化C特征分析

C位點(diǎn)的DNA甲基化主要有3種類型：mCG，mCHG和mCHH。這 3種類型的數(shù)量及構(gòu)成比例反映了特定物種的全基因組甲基化的特征。本研究發(fā)現(xiàn)，A34_1A，A34_2A，A34_3A中分別有32.00%，35.79%和30.82%的C位點(diǎn)發(fā)生了甲基化，且DNA甲基化變化主要有mCG，mCHH和mCHG 3種類型，其中A34_2A的各位點(diǎn)甲基化比例高于A34_1A，A34_3A的各位點(diǎn)甲基化比例低于A34_1A (表3)。對(duì)甲基化C堿基中CG，CHG和CHH的分布比例結(jié)果表明，苜蓿誘變前后其mCHH位點(diǎn)所占比例最高，變化趨勢(shì)為先升后降，而mCG和mCHG位點(diǎn)出現(xiàn)頻率表現(xiàn)為先降后升的變化趨勢(shì)(圖1)。研究結(jié)果表明，不同誘變處理紫花苜蓿后，C位點(diǎn)發(fā)生甲基化呈現(xiàn)出不同變化趨勢(shì)，不同類型的mC模式均發(fā)生了變化，說(shuō)明DNA甲基化水平變化在一定范圍內(nèi)是保持相對(duì)穩(wěn)定的。

表3 苜蓿不同分布類型甲基化C的數(shù)量及比例Table 3 The number and proportion methylated C in alfalfa

圖1 甲基化C堿基中CG，CHG和CHH的分布比例Fig.1 The distribution of mCG,mCHG and mCHH in all methylcytosine

2.1.2紫花苜蓿全基因組DNA甲基化水平分析 由供試苜蓿全基因組C位點(diǎn)DNA甲基化水平分析結(jié)果(表4)可知，EMS處理后苜蓿DNA甲基化水平由15.67% 提高到17.43%，達(dá)到差異顯著水平(P<0.05)，而60CO-γ射線輻射處理后苜蓿DNA甲基化水平由15.67% 下降到12.85%(P<0.05)。此外，除EMS處理下mCG水平由52.45%上升到54.89%(P<0.05)，mCHG和mCHH水平及60CO-γ射線輻射處理下mCG，mCHG和mCHH水平均呈現(xiàn)下降趨勢(shì)，由此表明，DNA甲基化主要發(fā)生在CG序列，CHH序列在mC中所占比例最小，從而說(shuō)明不同誘變處理甲基化的變化與mCG的類型變化關(guān)系更密切。

表4 不同誘變處理前后C，CG，CHG，CHH甲基化水平Table 4 Methylation levels of C，CG，CHG and CHH in different mutagenesis treatments

2.2 紫花苜蓿全基因組甲基化在基因區(qū)的分布規(guī)律

不同基因組元件甲基化水平可以反映出基因組的甲基化模式，同時(shí)也有利于充分理解與認(rèn)識(shí)甲基化功能機(jī)制。圖2是紫花苜蓿誘變前后基因不同區(qū)域上發(fā)生甲基化水平情況，由圖可知，在各基因區(qū)域上，CG甲基化水平最高，CHG甲基化水平次之，CHH甲基化水平最低。此外還可以發(fā)現(xiàn)，EMS處理在紫花苜蓿許多基因區(qū)域的CG甲基化水平都有所上升，說(shuō)明EMS處理在誘變過(guò)程中存在多個(gè)位點(diǎn)的超甲基化過(guò)程；而60CO-γ射線處理在紫花苜蓿許多基因區(qū)域的CG甲基化水平都有所下降，說(shuō)明60CO-γ射線處理在誘變過(guò)程中存在多個(gè)位點(diǎn)的去甲基化過(guò)程。研究結(jié)果表明，紫花苜蓿誘變過(guò)程中內(nèi)含子區(qū)域發(fā)生甲基化水平均較高，其次是外顯子區(qū)和啟動(dòng)子區(qū)域，5′UTR和3′UTR甲基化程度最低。這可能是因?yàn)橥怙@子及內(nèi)含子參與基因表達(dá)調(diào)控及防止了其序列的異常轉(zhuǎn)錄所導(dǎo)致的。

圖2 不同基因區(qū)域上各類型C的甲基化水平Fig.2 Methylation level of various types of C in different gene regions

2.3 紫花苜蓿DNA甲基化組的差異甲基化及KEGG通路富集分析

2.3.1DNA甲基化組的DMR分析 將苜蓿對(duì)照組與EMS處理組(A34-1A VS.A34-2A)及對(duì)照組與60CO-γ輻射處理組(A34-1A VS. A34-3A)比較中分別檢測(cè)出的3 066和1 964個(gè)DMRs。將其比對(duì)到苜蓿參考基因組中，發(fā)現(xiàn)分別有1 212和517個(gè)差異甲基化基因被DMR覆蓋到，這些基因被稱為DMR相關(guān)基因或差異甲基化基因(DMG)，其中有105個(gè)基因在兩組比較中均存在差異甲基化區(qū)域(見圖3)。將各組比較檢測(cè)出的DMR比對(duì)到不同的基因組元件上(圖4)，發(fā)現(xiàn)DMR主要分布于啟動(dòng)子區(qū)、內(nèi)含子區(qū)和外顯子區(qū)，而在5′UTR和3′UTR區(qū)的分布較少。

圖3 紫花苜蓿不同誘變處理差異甲基化基因Fig.3 Differentially methylated genes in each alfalfa mutagenesis comparison

圖4 DMR在不同基因組元件中的分布Fig.4 Distribution of DMRs in different genomic elements

2.3.2差異甲基化基因的KEGG通路富集分析 對(duì)所得到的DMR相關(guān)基因即差異甲基化基因進(jìn)行KEGG功能富集分析，挖掘出與苜蓿生長(zhǎng)及生理變化相關(guān)的生物學(xué)過(guò)程和pathway調(diào)控通路，進(jìn)而探討差異甲基化基因在紫花苜蓿生長(zhǎng)發(fā)育中的調(diào)控作用。對(duì)各組比較的差異甲基化基因進(jìn)行KEGG通路富集分析，結(jié)果表明，對(duì)照組和EMS處理組之間差異表達(dá)基因顯著富集于的Pathway(P<0.1)包括淀粉和蔗糖代謝途徑(Starch and sucrose metabolism)、氨基酸糖和核苷酸糖代謝(Amino sugar and nucleotide sugar metabolism)、氮代謝通路(Nitrogen metabolism)嘧啶代謝途徑(Pyrimidine metabolism)、精氨酸和脯氨酸代謝(Arginine and proline metabolism)等生理變化相關(guān)代謝及通路(圖5-1)，涉及到的基因共33個(gè)(表5)，其中氮代謝通路、淀粉和蔗糖代謝通路以及脯氨酸代謝通路是與苜蓿逆境生理息息相關(guān)的3個(gè)主要通路；對(duì)照組與60CO-γ輻射組的差異甲基化基因顯著富集于植物病原體互作(Plant-pathogen interaction)、鞘脂代謝通路(Sphingolipid metabolism)、代謝通路(Metabolic pathways)、次生代謝產(chǎn)物的生物合成(Biosynthesis of secondary metabolites)、生物合成植物激素信號(hào)轉(zhuǎn)導(dǎo)(Plant hormone signal transduction)、及氨基酸類生理代謝變化及通路(圖5-2)，涉及到的基因共9個(gè)(見表5)，這些通路參與調(diào)節(jié)細(xì)胞的生長(zhǎng)、分化、衰老和細(xì)胞程序性死亡等許多重要的信號(hào)轉(zhuǎn)導(dǎo)過(guò)程，使細(xì)胞產(chǎn)生各種不同的生物學(xué)功能。綜合通路富集分析結(jié)果可推測(cè)紫花苜蓿在受到誘變脅迫條件下，這些通路及相關(guān)基因在生理代謝過(guò)程中可能發(fā)揮著重要的調(diào)控作用。

圖5-1 A34-1A VS.A34-2A差異表達(dá)基因pathway富集分析Fig.5-1 The enrichment analysis of DGEs in the comparison of A34-1A VS.A34-2A

圖5-2 A34-1A VS.A34-3A差異表達(dá)基因pathway富集分析Fig.5-2 The enrichment analysis of DGEs in the comparison of A34-1A VS.A34-3A

表5 差異甲基化基因富集的與生理代謝相關(guān)通路及基因Table 5 The pathways and genes involved in physiological metabolism

3 討論

誘變能夠引起生物體基因組[22]、表達(dá)組[23]、蛋白質(zhì)組[24]和代謝組[25-26]的變化。其中DNA甲基化也被稱為基因組中的“第五堿基”，它是維持基因組穩(wěn)定性和調(diào)控基因表達(dá)的另一種方式，在植物各種脅迫中以及生長(zhǎng)發(fā)育過(guò)程中扮演著至關(guān)重要的角色[27]。有研究報(bào)道，DNA甲基化水平往往與物種基因組大小表現(xiàn)出一定的正相關(guān)關(guān)系[28-29]。前人研究CG位點(diǎn)的甲基化水平變化范圍從30.5%的擬南芥(Arabidopsis)到92.5%甜菜(Beet)；CHG甲基化水平的變化范圍從9.3%的鹽水水芹(Saltwater cress)到81.2%甜菜(Beet)，CHH甲基化水平變化范圍從1.1%葡萄(Grape)到18.8%甜菜(Beet)[30]。在本研究中，誘變前后紫花苜蓿基因組在C，CG，CHG和CHH位點(diǎn)的甲基化水平范圍分別為12.85%～15.67%，46.62%～54.89%，18.51%～22.68%和5.94%～8.35%。顯然與其他物種相比，紫花苜蓿DNA甲基化水平在變化范圍之內(nèi)，且紫花苜蓿基因組在CG和CHG序列環(huán)境下也存在相對(duì)較高的甲基化水平，且甲基化水平隨著不同誘變脅迫條件改變而變化。

此外，本研究發(fā)現(xiàn)，C位點(diǎn)的DNA甲基化主要為mCG，mCHH和mCHG 3種類型，其中EMS誘變后苜蓿甲基化水平上升，而60CO-γ輻射處理后苜蓿的甲基化水平下降，且甲基化C堿基中CG，CHG和CHH的分布比例因誘變處理不同而不同，其中CG類型甲基化胞嘧啶所占比例最高，這與Shi等[31]研究發(fā)現(xiàn)，成熟期水稻DNA甲基化在重離子輻射處理下發(fā)生在CG位點(diǎn)多于CNG位點(diǎn)的結(jié)論相一致。此外，也有研究表明CNG位點(diǎn)甲基化變化在植物生長(zhǎng)發(fā)育和抵御逆境過(guò)程中發(fā)揮著重要作用[32-34]。進(jìn)一步表明研究誘變后甲基化水平變化對(duì)突變體篩選工作具有重要作用。

不同基因組元件甲基化水平能夠反應(yīng)基因組的甲基化模式，有利于分析甲基化功能機(jī)制研究。本實(shí)驗(yàn)通過(guò)計(jì)算各元件內(nèi)胞嘧啶位點(diǎn)的平均甲基化水平，分析各功能元件的甲基化水平差異。結(jié)果發(fā)現(xiàn)CG，CHG，CHH各序列環(huán)境下基因組各元件甲基化水平均發(fā)生了不同程度的變化，其中外顯子及內(nèi)含子區(qū)域發(fā)生甲基化水平均很高。這可能是因?yàn)橥怙@子及內(nèi)含子在參與基因表達(dá)調(diào)控為了阻止其序列的異常轉(zhuǎn)錄所導(dǎo)致的。這與趙振利研究的白花泡桐叢枝病發(fā)生過(guò)程中全基因組DNA甲基化差異分析的研究結(jié)果相似[35]。

本研究進(jìn)行對(duì)照與EMS處理及對(duì)照與60Co-γ射線處理(A34-1A VS.A34-2A和A34-1A VS. A34-3A)的DNA甲基化進(jìn)行比較，分別檢測(cè)出的3066和1964個(gè)DMR。將其比對(duì)到參考基因組中，發(fā)現(xiàn)分別有1212和517個(gè)差異甲基化基因被DMR覆蓋到，且主要分布在啟動(dòng)子區(qū)，說(shuō)明大部分的甲基化差異主要發(fā)生在啟動(dòng)子區(qū)。將這些基因進(jìn)行KEGG功能富集分析，發(fā)現(xiàn)這些基因顯著富集于氮代謝途徑、嘧啶代謝途徑、多種氨基酸代謝途徑及糖代謝等生理變化相關(guān)的信號(hào)通路中(P< 0.05)。推測(cè)這些通路及基因在紫花苜蓿生理代謝過(guò)程中發(fā)揮著重要的調(diào)控作用。進(jìn)一步說(shuō)明，DNA甲基化在一定程度上調(diào)控基因的表達(dá)，是受外界因素調(diào)控的眾多基因的公共開關(guān)，通過(guò)改變生理代謝變化，從而促進(jìn)或者抑制作物生長(zhǎng)過(guò)程中某些能量物質(zhì)的代謝與轉(zhuǎn)換，進(jìn)而影響植物的生活力，協(xié)助植物渡過(guò)逆境脅迫[36]。

4 結(jié)論

本研究利用BS-seq技術(shù)對(duì)EMS及60CO-γ兩種誘變處理下紫花苜蓿全基因組DNA甲基化的變化進(jìn)行了研究，分析了DNA甲基化的水平和模式變化情況，發(fā)現(xiàn)EMS誘變后苜蓿甲基化水平上升，60CO-γ輻射處理后苜蓿的甲基化水平下降，對(duì)兩組處理進(jìn)行甲基化區(qū)域的差異研究及KEGG通路富集分析，發(fā)現(xiàn)了3個(gè)與逆境生理相關(guān)的信號(hào)通路，涉及到的相關(guān)基因共有42個(gè)。