王 笑, 劉 凱, 孫 尚, 薛 超, 龔志云

(揚州大學農(nóng)學院/ 江蘇省作物基因組學和分子育種重點實驗室/ 植物功能基因組學教育部重點實驗室/江蘇省作物遺傳生理重點實驗室/ 江蘇省糧食作物現(xiàn)代產(chǎn)業(yè)技術(shù)協(xié)同創(chuàng)新中心, 江蘇揚州 225009)

花藥培養(yǎng)是常用的植物組織培養(yǎng)技術(shù)之一,能產(chǎn)生具有特定親本特性的單倍體,然后通過自然或人工加倍獲得純合子二倍體,大大加快育種速度,因而廣泛應(yīng)用于培育各種水稻栽培品種[1]。但水稻花藥培養(yǎng)中常產(chǎn)生大量黃化苗,有時高達80%～90%[2]。這已成為花藥培養(yǎng)在水稻育種中進一步發(fā)展和應(yīng)用的重大障礙之一?；蛲蛔冊徽J為是導致幼苗黃化的主要原因,但其發(fā)生頻率如此之高并不符合基因突變的規(guī)律。因此,急需對水稻花藥培養(yǎng)中高頻黃化的分子和遺傳調(diào)控機制進行深入研究[3]。

愈傷組織的誘導和分化是水稻花藥培養(yǎng)的2個重要過程。脫分化及營養(yǎng)、光溫不足等不利因素均可能阻礙愈傷組織的形成。這些因素會刺激表觀遺傳調(diào)控的變化,因此黃化苗的出現(xiàn)可能與表觀遺傳調(diào)控有關(guān)[4]。盡管水稻基因組已完全測序,但關(guān)于水稻花藥培養(yǎng)表觀遺傳學的相關(guān)報道卻很少[5-6]。在植物生長過程中,重金屬、干旱、鹽、高溫等諸多不利因素經(jīng)常阻礙植株的正常生長。生物通常進化出相應(yīng)的調(diào)節(jié)機制(基因轉(zhuǎn)錄、翻譯調(diào)控或表觀遺傳學)來抵御這些不利因素[7-8]。DNA甲基化是一種不改變DNA序列的保守表觀遺傳學機制,由DNA甲基轉(zhuǎn)移酶催化基因組CpG二核苷酸中胞嘧啶第5碳位的甲基發(fā)生共價結(jié)合反應(yīng)而產(chǎn)生,通過改變?nèi)旧|(zhì)結(jié)構(gòu)、DNA構(gòu)象和DNA穩(wěn)定性導致基因表達水平的改變[9], 在植物基因組穩(wěn)定性、環(huán)境響應(yīng)和發(fā)育調(diào)控中具有重要意義[10]。以往研究[11]表明,水稻在組織培養(yǎng)或再生過程中易獲得或失去DNA甲基化,影響其正常生長發(fā)育。

本研究通過亞硫酸鹽測序(BS-Seq)和RNA測序(RNA-Seq)技術(shù),分別對水稻花藥培養(yǎng)中產(chǎn)生的綠苗和黃化苗進行轉(zhuǎn)錄組和甲基化組分析,確定它們之間的差異性; 利用SPAD葉綠素儀和分光光度計檢測黃化苗葉綠素含量,利用透射電鏡成像技術(shù)觀察葉綠體形態(tài)的變化; 通過對綠苗和黃化苗的差異DNA甲基化基因以及轉(zhuǎn)錄組水平中差異表達基因的分析,探究花藥培養(yǎng)過程中高頻黃化苗的產(chǎn)生可能與DNA甲基化的影響機制。

1 材料與方法

1.1 試驗材料

供試品種為粳稻品種日本晴。用濕潤的紗布包裹水稻幼穗,放入4 ℃冰箱預(yù)處理7 d。經(jīng)消毒、沖洗等處理后,將合適的穎花花藥均勻地撒在花藥培養(yǎng)基表面,置于(25±1) ℃的暗環(huán)境中培養(yǎng),培養(yǎng)3～4周后開始觀察花藥對誘導的反應(yīng)。然后將愈傷組織轉(zhuǎn)移到再生培養(yǎng)基上,在(25±1) ℃的人工光照下培養(yǎng)。最后將再生的綠苗和黃化苗轉(zhuǎn)移到生根培養(yǎng)基中形成根系。

1.2 光合色素含量的測定

取生長14 d的綠苗和黃化苗用于光合色素含量的測定。利用SPAD-502型葉綠素儀(日本柯尼卡美能達公司)測定葉片上、中、下部的透光系數(shù),重復(fù)3次,計算平均SPAD值。利用紫外/可見分光光度計(Lambda650型,美國鉑金埃爾默公司)測定5張葉片不同部位的色素含量。

1.3 葉綠體透射電鏡觀察

采集生長14 d的綠苗和黃化苗葉片,切成約1 mm×1 mm的小片。樣品用25%戊二醛磷酸鹽緩沖液(pH 7.4)在4 ℃下預(yù)處理4 h, 用0.1 mol·l-1PBS清洗3次, 每次15 min。接著用鋨酸固定 4 h, 再用0.1 mol·l-1PBS清洗3次, 每次15 min。然后用梯度乙醇(分別為50%、70%、80%、90%、95%、100%)溶液、100%丙醇和100%丙酮與無水Na2SO4逐級脫水。丙醇與樹脂1∶1混合浸泡1 h, 丙酮與樹脂1∶2浸泡2 h, 純樹脂浸泡2 h。用LR White樹脂浸漬, 60 ℃聚合2 d。用Lecia EM uc6超顯微切片儀(德國徠卡公司)制備超薄切片,超薄切片染色正常,并用Tecnai Spirit (120 kV)透射電鏡(荷蘭飛利浦公司)觀察葉綠體結(jié)構(gòu)。

1.4 DNA提取和BS-seq文庫構(gòu)建

使用DNeasy Plant Mini (Qiagen 69104, 德國凱杰公司)從綠苗和黃化苗的葉片中提取總基因組DNA (gDNA)。使用瓊脂糖凝膠電泳和紫外分光光度儀評估基因組DNA質(zhì)量,使用Bioruptor (比利時Diagenode公司)將gDNA切成100～500 bp的片段。在末端修復(fù)、腺苷酸化和接頭序列連接(以防止亞硫酸鹽轉(zhuǎn)化)后,以ZYMO甲基化試劑盒(美國ZYMO Research公司)用亞硫酸氫鹽處理DNA片段。處理后DNA在自旋柱上純化,用于制備測序文庫。在這一過程中,亞硫酸鹽處理的單鏈DNA隨機引物使用聚合酶能讀取尿嘧啶核苷酸合成含特定序列標簽的DNA。然后用這些標簽通過PCR在原始DNA鏈的5′和3′端添加接頭序列。表觀基因組庫被稀釋并加載到cBot DNA集群生成系統(tǒng)。聚類生成完成后,將2個樣本轉(zhuǎn)移到Illumina HiSeq 2000平臺(北京諾禾致源公司)進行測序。所有步驟均按照制造商的說明進行。

1.5 BS-Seq數(shù)據(jù)分析

首先對Illumina測序獲得的原始數(shù)據(jù)進行處理,過濾掉含有接頭序列、未知或低質(zhì)量堿基的reads, 然后用比對軟件Bowtie 2 (http://bowtiebio.sourceforge.net/bowtie2/index.shtml)將唯一的凈reads映射到水稻參考基因組(https://phytozome.jgi.doe.gov/pz/portal.html#info?alias=Org_Osativa)?；蚪M元素的鑒定采用Li等[12]方法進行。

利用BS-Seq數(shù)據(jù)檢測單個堿基的甲基化狀態(tài)、甲基化胞嘧啶位點的數(shù)量以及每個基因組原色的甲基化比例。甲基化水平(%)由每個甲基化胞嘧啶(mC)的讀數(shù)除以胞嘧啶的總讀數(shù)得到。差異甲基化區(qū)域(DMRs)包含至少5個CG (CHH或CHG)位點,甲基化水平發(fā)生2倍變化(Fisher精確檢驗,P≤0.05)。在多次試驗中,采用錯誤發(fā)現(xiàn)率(FDR)方法確定P值的閾值。使用R包進行基礎(chǔ)數(shù)據(jù)操作和統(tǒng)計分析。

當1個基因上至少有1個DMR時,該基因被認為是差異甲基化基因(DMG)。使用Blast2GO軟件獲取每個DMG的GO信息,利用WEGO網(wǎng)站(http://wego.genomics.org.cn/)上傳每個DMG的注釋信息,并根據(jù)功能對DMGs進行分類。功能意義的GO富集分析應(yīng)用超幾何測試將所有DMGs映射到GO數(shù)據(jù)庫中,在與全基因組相比的差異甲基化基因中搜索顯著富集的GO術(shù)語[13]。

與全基因組背景相比, KEGG通路富集分析可確定DMGs中顯著富集的代謝通路或信號轉(zhuǎn)導通路。計算公式和應(yīng)用程序與GO分析相同。隨后,以FDR (q值)≤0.05作為閾值來確定DMGs中最顯著富集的通路。

使用MapMan軟件(http://mapman.gabipb.org)描述與DMGs相關(guān)的生物途徑[14]。

為提供高或低DMGs的描述,使用Cluster軟件和JAVA TreeView軟件進行聚類分析。將一份包含2個樣本中DMGs及其甲基化比例的文件上傳到Gene Cluster 3.0。經(jīng)過數(shù)據(jù)過濾和日志轉(zhuǎn)換后,采用平均連鎖法對基因和陣列進行分層聚類。利用Cluster軟件生成的格式文件結(jié)合JAVA TreeView軟件,可直觀地觀察圖形化顯示的聚類結(jié)果。

1.6 RNA提取和轉(zhuǎn)錄組測序

使用TaKaRa RNAiso Plus (總RNA提取試劑,日本TaKaRa公司), 參照試劑盒的標準操作程序,從用于甲基組分析的同一綠苗和黃化苗中分離總RNA。使用NanoDrop ND-2000分光光度計(美國賽默飛世爾公司)和Agilent生物分析儀(2100, 美國安捷倫公司)檢測RNA濃度和質(zhì)量。使用RNA Clean XP Kit (Cat#A63987, 美國見克曼庫爾特公司)和RNase-Free DNase Set (Cat#79254, 德國凱杰公司)純化總RNA。

每個樣品使用5 μg RNA進行文庫構(gòu)建。文庫測序由上海伯豪生物公司完成。利用IlluminaHiSeq 2000/2500、MiSeq、Next Generation Sequencing (NGS)技術(shù)對綠苗和黃化苗的cDNA進行測序。移除接頭序列和低質(zhì)量reads, 執(zhí)行端測序。使用HISAT和Bowtie2工具將凈reads比對到水稻基因組注釋項目(http://rice.plantbiology.msu.edu/)?；虮磉_分析使用RESM軟件計算FPKM值。差異表達基因(DEGs)采用NOISeq軟件進行篩選,篩選標準為fold change≥2且P≤0.05。利用GO和KEGG工具對DEGs進行分析。

1.7 RNA提取和RT-qPCR分析基因表達

用TaKaRa RNAiso Plus (總RNA提取試劑)從生長14 d的綠苗和黃化苗幼葉中提取和純化總RNA。用TaKaRa PrimeScriptTM RT Reagent Kit合成第1鏈cDNA, 用不含RNase的ddH2O稀釋反轉(zhuǎn)錄產(chǎn)物。RT-qPCR分析使用TaKaRa SYBR Premis Ex Taq II (Perfect Real Time)試劑盒和ABI 7300實時定量PCR系統(tǒng)(美國Applied Biosystems公司)。泛素蛋白UBQ引物的UBQ-1作為內(nèi)參基因。采用循環(huán)閾值(Ct) 2-ΔΔCt方法計算每個樣本中所選基因的相對表達量。進行3次獨立的生物學重復(fù)。用于RT-qPCR的基因特異性引物見表1。

表1 RT-qPCR和重亞硫酸鹽測序試驗使用的引物Tab.1 Primers used in RT-qPCR experiments and Bisulfite-PCR analysis

1.8 DNA提取與重亞硫酸鹽測序

從生長14 d的綠苗和黃化苗幼葉中提取750 ng的DNA, 抽提的DNA用ZYMO甲基化試劑盒進行亞硫酸氫鹽處理。使用Free Bisulfite Primer Design Tool (https://www.zymoresearch.com/ pages/bisulfite-primer-seeker)設(shè)計所選基因序列的引物。簡并引物使用Ex Taq Hot Start (日本TaKaRa公司)從處理過的DNA中擴增目標片段。擴增產(chǎn)物使用Axygen凝膠提取試劑盒(美國Axygen公司)純化,所有片段連接到pMDTM18 T載體(日本TaKaRa公司)中,挑選至少6個獨立的單克隆由尚亞(杭州)生物技術(shù)公司測序。然后使用Kismeth (http://katahdin.mssm.edu/Kismeth)進行甲基化分析[15]。所用引物見表1。

2 結(jié)果與分析

2.1 水稻花藥培養(yǎng)中黃化苗的形成

采用日本晴(OryzasativaL. subspjaponicacv. Nipponbare)進行水稻花藥培養(yǎng)。花藥形成愈傷組織,愈傷組織再分化成芽和根,逐漸發(fā)育成幼苗。在這一過程中會形成大量的黃化苗(圖1.A、B), 2周后,綠苗和黃化苗被轉(zhuǎn)移到再生培養(yǎng)基中(圖1.C)。102根試管內(nèi)的1 126粒花藥形成293個愈傷組織,最終形成269株幼苗,其中綠苗為177株,黃化苗為92株,黃化苗率為34.20%。

圖1 綠苗和黃化苗的表型和葉綠體超微結(jié)構(gòu)分析*Fig.1 Phenotype and Chloroplast ultrastructure of green and etiolated seedlings* A. 花藥培養(yǎng)愈傷組織誘導; B. 愈傷組織再生形成綠苗和黃化苗; C. 綠苗和黃化苗; D. 綠苗葉綠體超微結(jié)構(gòu)(標尺=2 μm); E. 綠苗葉綠體超微結(jié)構(gòu)(標尺=0.5 μm); F. 黃化苗葉綠體超微結(jié)構(gòu)(標尺=2 μm); G. 黃化苗葉綠體超微結(jié)構(gòu)(標尺=0.3 μm); H. 綠苗和黃化苗SPAD葉綠素值(0～100無單位); I. 綠苗和黃化苗葉綠素含量為分光光度計的測定值。* A. induction of callus from cultured rice anthers; B. regenerated rice green and etiolated seedlings from callus; C. green and etiolated seedlings; D and E. ultrastructure of the chloroplast of green seedlings, bar=2 μm, bar=0.5 μm, respectively; F and G. ultrastructure of the chloroplast of etiolated seedlings, bar=2 μm, bar=0.5 μm, respectively; H. SPAD chlorophyll value (a unitless index from 0 to 100) in green and etiolated seedlings; I. chlorophyll content measured by spectrophotometer.

與綠苗相比,黃化苗的SPAD值顯著降低,說明相對葉綠素含量較低(圖1.H)。利用細胞學鑒定技術(shù),選擇二倍體的綠苗和黃化苗開展研究。植物光合色素主要分為葉綠素a、葉綠素b和胡蘿卜素等,位于葉綠體的類囊體膜上。黃化苗葉片上的葉綠素a、葉綠素b和胡蘿卜素含量顯著降低(圖1.I)。透射電鏡觀察發(fā)現(xiàn),綠苗綠葉中葉綠體結(jié)構(gòu)完整,可見類囊體結(jié)構(gòu)和明顯的片層(圖1.D、E), 但在黃化苗的幼葉中,類囊體片層不明顯,液泡占據(jù)細胞(圖1.F、G)。推測黃化苗的形成可能與葉綠體合成相關(guān)基因表達受到抑制有關(guān)。

2.2 綠苗和黃化苗的全基因組甲基化分析

對來自水稻花藥培養(yǎng)的綠苗和黃化苗樣品進行單堿基亞硫酸氫鹽測序,從每個樣本中過濾出大約7 780萬和6 070萬個凈reads。綠苗和黃化苗樣品的凈率均超過91%, 表明測序質(zhì)量較高,數(shù)據(jù)可靠,可用于高質(zhì)量的全基因組甲基化分析。

綠苗與黃化苗之間DNA甲基化水平的全基因組整體差異不顯著。在綠苗基因組中檢測到超過15億個胞嘧啶位點,在黃化苗基因組中檢測到超過18.4億個胞嘧啶位點,其中CG位點的mCs (Methylated C’s, 甲基化的C總數(shù)目)頻率最高, CHG和CHH位點的的mCs頻率最低(圖2)。通過比較,綠苗和黃化苗在CHH位點略有差異,黃化苗在CG和CHH位點的甲基化水平明顯高于綠苗(表2), 表明CG、CHG和CHH序列背景下的DNA甲基化受黃化的影響并不均勻。與綠苗相比,黃化苗在mCG、mCHG和mCHH位點檢測到的mCs數(shù)量增加(圖3.A-C), 表明在黃化苗形成過程中出現(xiàn)大量新的甲基化。綠苗和黃化苗之間的mC位點的比例在不同序列背景下也不同。結(jié)果發(fā)現(xiàn)mCHG位點的數(shù)量較mCHH位點的數(shù)量變化更多,而mCHH位點的數(shù)量較mCG位點的數(shù)量變化更多(圖3.A-C)。

圖2 綠苗和黃化苗的整個基因組中C、CpG、CHG和CHH的平均甲基化水平Fig.2 The average methylation levels of C, CpG, CHG and CHH in the whole genome of green and etiolated seedlings

圖3 綠苗和黃化苗的DNA甲基化水平*Fig.3 The DNA methylome of green and etiolated seedlings* A-C. 綠苗和黃化苗之間mCs在mCG、mCHG、mCHH位點的比較; D-F. 綠苗和黃化苗在每一背景下水稻基因組不同區(qū)域甲基化序列的百分比; R1. 啟動子, R2. 5′-非翻譯區(qū), R3. 外顯子, R4. 內(nèi)含子, R5. 3′-非翻譯區(qū), R6. 基因間區(qū), R7. 相關(guān)轉(zhuǎn)座因子。* A-C. comparison of mCs between green and etiolated seedlings for mCG, mCHG, mCHH sites; D-F. total percentage of methylated sequences in different regions of the rice genome for each context in green and etiolated seedlings; R1. promoter, R2. 5′-UTR, R3. exon, R4. intron, R5. 3′-UTR, R6. intergenic, R7. TE-related.

表2 綠苗和黃化苗整個基因組中C、CpG、CHG和CHH的平均甲基化水平Tab.2 The average methylation levels of C, CpG, CHG and CHH in the whole genome of green and etiolated seedlings

基于檢測到的甲基化位點,計算了3個序列背景(mCG、mCHG和mCHH)中綠苗和黃化苗的7個區(qū)域(啟動子、5′-UTR、外顯子、內(nèi)含子、3′-UTR、基因間和TE相關(guān)區(qū)域)的平均甲基化水平(圖3.D-F), 發(fā)現(xiàn)在黃化苗的基因間區(qū)中mCG位點的比例略高。有趣的是,黃化苗中mCHG位點啟動子區(qū)、內(nèi)含子區(qū)、3′-UTR和TE相關(guān)區(qū)域的相對比例分別高于綠苗,而黃化苗中mCHH位點啟動子區(qū)、內(nèi)含子區(qū)、3′-UTR 和基因間區(qū)相對比例均高于綠苗?？傮w而言,綠苗和黃化苗中mCG和mCHG的7個區(qū)域的水平高度一致,其中,基因間區(qū)甲基化比例最多,其次是外顯子區(qū)域。為更好地理解綠苗和黃化苗的基因和轉(zhuǎn)座因子中DNA甲基化的不同模式,將基因組劃分為幾個功能基因區(qū)域: 啟動子區(qū)、轉(zhuǎn)錄起始位點上游2 kb區(qū)域(TSS)、基因體、轉(zhuǎn)錄終止位點(TTS)下游2 kb處(圖4.A-D)。據(jù)分析,mCG的甲基化水平最高,mCHG次之,mCHH再次之。總體而言,綠苗和黃化苗的CG、CHG和CHH甲基化趨勢無明顯變化。

圖4 基因特征和染色體的DNA甲基化模式*Fig.4 DNA methylation pattern in gene features and chromosomesanalysis* A、C分別表示綠苗和黃化苗標準化甲基化水平, 這些區(qū)域是由轉(zhuǎn)錄起始位點(TSS)上游的2 kb, 以及轉(zhuǎn)錄區(qū)和轉(zhuǎn)錄終止位點(TTS)下游的2 kb構(gòu)成; B、D分別代表mCG、mCHG和mCHH背景下基因功能區(qū)域的DNA甲基化模式, a. 啟動子, b. 5′-UTR, c. 外顯子, d. 內(nèi)含子, e. 3′-UTR; E、F分別代表綠苗和黃化苗12條染色體上mCG、mCHG和mCHH背景下基因表達水平、轉(zhuǎn)錄因子(TEs)密度與甲基化水平的關(guān)系。* A and C represent the normalized methylation levels of green and etiolated seedlings, respectively, the plots were generated from 2 kb upstream of the transcription start site (TSS) to 2 kb downstream of the transcription termination site (TTS); B and D represent the DNA methylation patterns of gene functional regions in mCG, mCHG and mCHH environments in green and etiolated seedlings, respectively, a. promoter, b. 5′-UTR, c. exon, d. intron, e. 3′-UTR; E and F represent the circular plots of gene expression levels, transposable element (TE) density and methylation levels in mCG, mCHG and mCHH environments in rice green and etiolated seedlings, respectively.

為進一步確定12條染色體的DNA甲基化水平(mCG、mCHG和mCHH序列背景)轉(zhuǎn)座元件(TEs)和基因表達之間的相關(guān)性,分析了染色體位置關(guān)系、基因注釋和TEs信息,結(jié)果表明,在染色體水平上綠苗和黃化苗的DNA甲基化水平并沒有明顯的差異(圖4.E、F)。然而在不同序列背景下,mCG的甲基化水平最高,而mCHH的甲基化水平最低?；蜣D(zhuǎn)錄水平和DNA甲基化水平在每個染色體上的總體分布呈現(xiàn)出相反的模式(圖4.E、F)。在水稻基因組中, DNA甲基化可能與TE密度呈正相關(guān),與基因表達呈負相關(guān),但仍未被大規(guī)模鑒定無論是在基因組規(guī)模上還是在特定位點上, CHH甲基化的水平往往低于CG或CHG甲基化水平[16]。其中, 1號染色體的0～5 Mb區(qū)域甲基化水平顯著變化, 4號染色體的17～19 Mb區(qū)域甲基化水平顯著變化, 9號染色體9～12 Mb以及11號染色體的6～8和26～28 Mb區(qū)域類似且甲基化水平顯著變化。此外,每條染色體上的紅色線(mCG)水平明顯出現(xiàn)整體水平的上調(diào),而每條染色體上的綠色線(mCHH)未發(fā)生明顯變化(圖5、6)。

圖5 綠苗和黃化苗12條染色體基因組的DNA甲基化模式*Fig.5 Dynamics of DNA methylation in the genome of 12 chromosomes between green seedlings and etiolated seedlings* 黃化苗與綠苗的每條染色體上甲基化水平差異的比值; 紅色代表mCG, 藍線代表mCHG, 綠線代表mCHH。* tthe ratio of differences in methylation levels on each chromosome of etiolated seedlings to that of green seedlings; red linerepresents mCG, blue linerepresents mCHG, and green linerepresents mCHH.

圖6 綠苗和黃化苗的染色體上差異甲基化位點及區(qū)域的分布*Fig.6 Distribution of differential methylation on chromosomes among green and etiolated seedlings* A和C分別為綠苗和黃化苗每條染色體上不同甲基化位點的比例, 紅色代表高甲基化位點, 藍色代表低甲基化位點; B和D分別為綠苗和黃化苗中每條染色體上不同甲基化區(qū)域的比例, 紅色代表高甲基化區(qū)域, 藍色代表低甲基化區(qū)域。* A and B indicate the number of differentially methylated sites on each chromosome in green and etiolated seedlings, respectively, red color indicates hypermethylated sites and blue color indicates hypomethylated sites; B and D indicate the number of differentially methylated regions (DMRs) on each chromosome in green and etiolated seedlings, respectively, red color indicates hypermethylated regions and blue color indicates hypomethylated regions.

2.3 綠苗和黃化苗的差異甲基化區(qū)域和GO富集分析

通過對綠苗和黃化苗DNA甲基化序列的全基因組比較,本研究鑒定了1 058個DMRs。在功能區(qū)中,與基因間相關(guān)的DMRs占比最大(圖7.A)。從具有基因組特征的631個高甲基化DMRs和427個低甲基化DMRs中,鑒定出4 156個DMRs相關(guān)基因(圖7.B)。在綠苗和黃化苗的CG位點處檢測到2 769個DMRs相關(guān)基因,在CHG位點檢測到1 095個DMRs相關(guān)基因,在CHH位點檢測到292個DMRs相關(guān)基因,共獲得1 907個啟動子相關(guān)的DMGs和2 249個基因本體相關(guān)的DMGs, 這與早前的研究[17]結(jié)果一致,甲基化差異主要集中在CG位點,表明基因本體的DNA甲基化在整個基因組中受到廣泛調(diào)控,啟動子區(qū)域的甲基化相對次之。

圖7 綠苗和黃化苗甲基化差異區(qū)分布*Fig.7 Distribution of differentially methylated regions in green and etiolated seedlings* A. 基因組上高甲基化和低甲基化區(qū)域(DMRs)的分布; B. 綠苗和黃化苗的不同甲基化序列背景下相關(guān)基因的基因特征分布; C和D分別表示啟動子和基因體相關(guān)基因的GO分析。* A. distribution of hyper-and hypomethylated regions (DMRS) on the genome; B. the distribution of different methylated context-associated genes among the gene features in etiolated/green seedlings; C and D represent the Gene ontology (GO) analysis of genes associated with promoters and gene body, respectively.

為描述基因及其產(chǎn)物的性質(zhì),使用WEGO對DMGs進行功能分類。GO分析顯示,相對于綠苗,黃化苗參與多種分子功能和生物過程,包括有機物代謝過程、初級代謝過程、細胞代謝過程、催化活性、雜環(huán)化合物結(jié)合和細胞解剖實體(圖7.C), 表明與基因本體相關(guān)的甲基化以及與啟動子相關(guān)的甲基化各自占有一定比例。CHH序列背景的甲基化是基因本體內(nèi)的主要甲基化位點,而CG序列背景甲基化是基因本體內(nèi)的主要甲基化位點。通過對DMRs相關(guān)啟動子的GO分析,發(fā)現(xiàn)與參與有機物質(zhì)代謝過程、初級代謝過程、細胞代謝過程、細胞周期過程、結(jié)合、催化活性和雜環(huán)化合物結(jié)合等通路相關(guān)基因有關(guān)(圖7.D)。

2.4 綠苗和黃化苗的差異基因表達分析

有研究[18-19]顯示,基因表達與DNA甲基化水平呈負相關(guān)性。本研究基于RNA-Seq分析了綠苗和黃化苗之間的差異基因表達,測序數(shù)據(jù)表明測序結(jié)果良好、數(shù)據(jù)質(zhì)量合格,可用于基因組比對。綠苗和黃化苗樣品中唯一的比較reads分別為4 768萬和6 731萬。通過測序獲得了2個樣品的每百萬個作圖讀數(shù)的外顯子模型的每千堿基片段(FPKM)和倍數(shù)變化(FC)值。由于CG序列背景是基因本體內(nèi)的主要甲基化位點,因此,差異基因表達分析以CG序列背景進行,發(fā)現(xiàn)綠苗和黃化苗中大多數(shù)基因在CG位點的表達水平相似,僅有411個基因上調(diào), 646個基因下調(diào)(圖8.A)。

圖8 綠苗和黃化苗基因表達差異及GO富集分析*Fig.8 Differential gene expression among green and etiolated seedlings and GO enrichment analysis* A. CG序列背景下綠苗和黃化苗之間基因轉(zhuǎn)錄的火山圖, DEGs的閾值為P≤0.05且倍數(shù)>2, 紅點和藍點分別代表表達上調(diào)和下調(diào)的基因, 黑點代表無顯著表達的基因; B. DEGs的GO富集分析; C. DEGs的KEGG通路富集分析; D. 不同序列背景下甲基化變化的DEGs比較; E. 在一種或多種序列背景下甲基化差異的基因與DEGs之間的比較?！癎ene body”和“promoter”分別代表轉(zhuǎn)錄區(qū)和基因上游區(qū)域1 kb,“DEGs”代表綠苗和黃化苗RNA-Seq中的差異表達基因。* A. volcano plot indicating gene transcription between green and etiolated seedlings in CG context, DEGs were identified by using a threshold of P≤0.05 and fold change>2, red and blue dots represent genes with up-and down-regulated expression, respectively, whereas the black dots represent the genes with no significant expression; B. GO enrichment analysis of differentially expressed genes; C. KEGG pathway enrichment analysis of DEGs; D. comparison of DEGs according to methylation variations among sequence contexts; E. comparison of genes differentially methylated at one or more sequence contexts and differentially expressed between green and etiolated seedlings. ‘Gene body’ and ‘promoter’ represent the transcribed region and the 1 kb upstream region of genes, respectively; ‘DEGs’ represents genes that are differentially expressed genes in RNA-seq among green and etiolated seedlings.

利用GO分析綠苗和黃化苗在CG位點的1 057個DEGs, 結(jié)果表明,在生物過程方面主要參與多糖結(jié)合、蛋白水解調(diào)控、內(nèi)肽酶活性負調(diào)控、對低光強度刺激的響應(yīng)、光合作用、蛋白質(zhì)-色素細胞聯(lián)動、光采集在光系統(tǒng)I中的反應(yīng); 細胞組分方面主要與類囊體膜、質(zhì)體類囊體膜、光系統(tǒng)Ⅱ、質(zhì)體小球和光系統(tǒng)Ⅰ有關(guān); 分子功能方面主要與鐵離子結(jié)合、谷氨酰胺-氨連接酶活性、蛋白質(zhì)結(jié)構(gòu)域特異性結(jié)合、四吡咯結(jié)合和葉綠素結(jié)合有關(guān)(圖8.B)。KEGG富集分析發(fā)現(xiàn), DEGs在代謝途徑顯著富集,包括次生代謝物生物合成途徑、淀粉和蔗糖代謝、光合作用-天線蛋白、α-亞麻酸代謝、丙氨酸-天冬氨酸-谷氨酸代謝和氮代謝(圖8.C) GO與KEGG富集分析表明與光合作用相關(guān)的途徑均與蛋白質(zhì)表達的改變有關(guān)。

試驗鑒定出的644個基因在綠苗和黃化苗之間的甲基化和基因表達均表現(xiàn)出顯著差異(以下稱為共差異基因), 只有0.5% (644個基因中的3個)在所有3個序列背景中存在差異甲基化(圖8.D), 分別是編碼葉綠體酶基因LOC_Os06g01850(OsLFNR1, 鐵氧還蛋白-NADP+氧化還原酶)、基因LOC_Os01g41710(假定表達基因葉綠素A-B結(jié)合蛋白)和基因LOC_Os02g10390(假定表達基因葉綠素A-B 結(jié)合蛋白), 其中OsLFNR1在光合作用電子傳遞鏈的最后步驟中發(fā)揮功能[20]。這表明甲基化對葉綠素合成過程中的轉(zhuǎn)錄水平具有重要影響。此外,綠苗和黃化苗差異表達的644個基因中, 490個基因含有mCG位點, 122個基因含有mCHG位點, 74個基因含有mCHH位點。進一步比較顯示,在綠苗和黃化苗中,只有8個基因同時在基因本體和啟動子區(qū)域發(fā)生不同程度的甲基化(圖8.E), 表明基因本體和啟動子區(qū)域的甲基化在很大程度上是由獨立的機制分別調(diào)節(jié)的。

2.5 綠苗和黃化苗光合作用途徑的共差異基因

葉綠體生物合成在光合作用和植物生長中發(fā)揮著至關(guān)重要的作用,以響應(yīng)環(huán)境變化。為研究綠苗和黃化苗之間的差異基因表達,重點研究參與葉綠體生物合成相關(guān)途徑的共差異基因。根據(jù)GO和KEGG富集分析,確定了9個與綠苗和黃化苗光合作用途徑有潛在關(guān)系的基因,這些基因是甲基化上調(diào)而表達下調(diào)的基因(表3), 其中包括編碼葉綠素生物合成酶的基因OsPORA[21], 在水稻光合作用和葉綠體發(fā)育中發(fā)揮功能的基因OsTLP27[22], 編碼在光合作用電子傳遞鏈的最后步驟中發(fā)揮功能的一個重要的葉綠體酶基因OsLFNR1[22]參與光系統(tǒng)而影響光合作用能力的基因OsNPH1a[23], 通過影響水稻葉片光合率和糖代謝而參與葉綠素積累、葉綠體發(fā)育以及植株生長的基因OsAld-Y51[24], 具有轉(zhuǎn)氨酶作用的類Ⅰ和類Ⅱ結(jié)構(gòu)域含蛋白質(zhì)的表達基因LOC_Os07g01760, 以及假定的類囊體管腔 29.8 ku 蛋白質(zhì)的編碼基因LOC_Os12g08830和參與光系統(tǒng)Ⅱ組裝的基因LOC_Os07g05360。

表3 9個可能與光合作用相關(guān)的基因Tab.3 Nine genes potentially related to photosynthesis pathways

通過RT-qPCR分析發(fā)現(xiàn),上述9個DMGs中有8個基因的表達水平下降,這與轉(zhuǎn)錄組測序結(jié)果一致(圖9.A)。為驗證甲基化模式的結(jié)果,隨機選擇2個基因使用重亞硫酸鹽測序分析?；贐S-Seq結(jié)果,重點研究參與光合過程光系統(tǒng)II組裝的基因LOC_Os07g05360和參與結(jié)構(gòu)域含蛋白的表達基因LOC_Os07g01760, 結(jié)果顯示這2個基因的甲基化水平分別為18.145%和91.891%, 均處于高甲基化水平,重亞硫酸鹽測序的結(jié)果與甲基化組分析的結(jié)果一致(圖9.B), 表明測序結(jié)果可靠。

圖9 9個DMGs的RT-qPCR和重亞硫酸鹽測序驗證*Fig.9 Expression level of nine DMGs validated by RT-qPCR and bisulfite-PCR* A. 9個DMGs在綠苗和黃化苗的RT-qPCR分析,所有的RT-qPCR反應(yīng), 3個生物重復(fù); B. 綠苗和黃化苗LOC_Os07g05360和LOC_Os07g01760的DNA甲基化水平。* A. RT-qPCR analysis of nine DMGs in green and etiolated seedlings, all RT-qPCR reactions were performed with three biological replicates; B. DNA methylation levels of LOC_Os07g05360 and LOC_Os07g01760 in green and etiolated seedlings.

3 討論

表觀遺傳修飾在植物逆境響應(yīng)機制中起著重要作用。環(huán)境刺激如鹽和滲透脅迫可導致某些基因編碼區(qū)的去甲基化,從而激活其表達。在某些不利條件下,非生物脅迫也會導致DNA甲基化和基因表達模式的改變。有研究[25]表明,植物響應(yīng)高溫、低溫、干旱、鹽堿、重金屬等非生物脅迫過程中,通過全基因組或特定基因DNA甲基化水平變化,進而調(diào)控抗逆基因的表達以更好地適應(yīng)生存環(huán)境。此外,在轉(zhuǎn)基因過程中, DNA甲基化參與調(diào)節(jié)植物開花、花器官形成和植物基因沉默[26-27]。水稻花藥培養(yǎng)簡單、省時、適應(yīng)性廣泛,但高頻的黃化現(xiàn)象嚴重影響著水稻的育種進程。葉綠體生物合成是一個非常復(fù)雜的過程[28], 通過一系列協(xié)同反應(yīng),且有多種酶催化完成的[29]。葉綠體負責光合作用以及激素與代謝物的產(chǎn)生[30-31], 這些激素和代謝物是由存在于植物分生組織未成熟細胞中前質(zhì)體發(fā)展而來的。葉綠體的生物發(fā)生涉及質(zhì)體和核編碼基因[32]的協(xié)調(diào)功能。葉綠體的生物發(fā)生過程中伴隨著葉綠素、葉綠體蛋白、類胡蘿卜素和酯類的生物合成,以及光合蛋白復(fù)合體的組裝,包括捕光復(fù)合體(LHC)、光系統(tǒng)Ⅰ (PSⅠ)、光系統(tǒng)Ⅱ (PSⅡ)、細胞色素f/be (cytf/bs)和三磷酸腺苷(ATP)合酶[30,33]。水稻花藥培養(yǎng)的綠苗和黃化苗,經(jīng)轉(zhuǎn)錄組分析,發(fā)現(xiàn)通過RT-qPCR試驗鑒定出9個與光合途徑相關(guān)的表達下調(diào)的基因,同時經(jīng)重亞硫酸鹽測序分析,隨機鑒定其中2個DNA甲基化上調(diào)的基因。其中,LOC_Os07g053600[34]是光反應(yīng)和卡爾文循環(huán)中起作用的基因。LOC_Os07g01760[35]是一種氨基酸轉(zhuǎn)移酶,推測在多胺降解中起作用,這些基因在光合途徑中均起著至關(guān)重要的作用。

通過對水稻花藥培養(yǎng)獲得的綠苗和黃化苗觀察和細胞超微結(jié)構(gòu)分析發(fā)現(xiàn),黃化苗葉片黃化,葉綠體發(fā)育不完整,葉綠素含量極少,生長發(fā)育受阻等,造成水稻幼苗不能進行正常的光合作用,無法正常生長。同時,經(jīng)驗證的9個與光合作用相關(guān)的基因表達下調(diào),故推測DNA甲基化會影響葉綠體生物合成相關(guān)基因,進而引起葉綠體生物合成受阻,導致高頻黃化苗的產(chǎn)生。本研究結(jié)果說明, DNA甲基化與基因表達之間的關(guān)系比早前認識要復(fù)雜,而啟動子的甲基化通常與基因表達抑制有關(guān),有關(guān)基因甲基化的作用則更不確定,既有正向也有反向的調(diào)控關(guān)系。同時,本研究對水稻綠苗和黃化苗葉片差異表達基因和差異甲基化基因的功能注釋和富集分析有助于探討DNA甲基化在水稻葉片中各代謝通路的表達與調(diào)控,為進一步研究DNA甲基化對花藥培養(yǎng)中黃化苗的產(chǎn)生以及光合作用的影響機制提供理論基礎(chǔ)。