甘藍型油菜鹽脅迫響應基因表達譜分析及共表達網絡的構建

楊 闖 王 玲 全成滔 余良倩 戴 成 郭 亮 傅廷棟 馬朝芝

華中農業大學作物遺傳改良全國重點實驗室 / 國家油菜工程技術研究中心 / 洪山實驗室, 湖北武漢 430070

土壤鹽漬化是威脅全球農業生產的主要環境挑戰之一[1]。世界上大約有20%的灌溉農田會受到土壤鹽漬化的不利影響[2]。自然環境惡化、灌溉方式不當以及氣候變化等因素也正在進一步加劇土壤鹽漬化問題[3]。預計到2050年, 鹽漬化土壤面積將達到全球可利用耕地的50%以上[4]。近些年來, 我國的農業耕地也在不斷受到鹽漬化侵蝕, 鹽漬地總面積已達到3.69×107hm2, 給作物生產和農業可持續發展帶來了極大的挑戰[5-6]。油菜是世界四大油料作物之一, 具有重要的油用、飼用、菜用以及旅游觀光價值[7]。作為食用植物油和植物蛋白的主要來源, 油菜在我國農產品中占有重要地位, 對于保障我國油料供給安全尤為重要[8]。然而, 鹽脅迫目前是影響油菜生長發育的主要環境條件之一, 可能會造成油菜減產、品質下降甚至死亡。因此, 解析甘藍型油菜響應鹽脅迫的機制, 并從中挖掘油菜響應鹽脅迫的核心基因, 對于油菜的耐鹽性改良具有重要意義。

鹽脅迫對植物的毒害作用主要包括滲透脅迫和離子脅迫2種[9]。鹽脅迫的滲透作用會導致土壤溶液的水勢低于根細胞水勢, 從而抑制根系的吸水作用; 此外, 滲透脅迫還會導致葉片氣孔關閉從而使光合效率下降, 最終影響植物的正常生長發育。離子脅迫主要是由細胞中的Na+和Cl-過量積累導致的。Na+的毒性主要是通過抑制新陳代謝途徑中關鍵酶活性從而影響植物的生長發育; 同時,過量的Na+也會影響植物對氮、磷、鉀、鈣等其他礦質元素的吸收利用。由于植物對NO3-、SO42-以及Cl-的吸收都是通過相同的非選擇性陰離子通道, 因此過量的Cl-可能會導致氮素以及硫素的缺乏。

植物為了響應鹽脅迫, 也進化出了互相關聯的調節通路并引起細胞內的廣泛變化, 其中許多變化是適應性反應, 可導致抗逆性增加, 這也是植物在長期進化過程中保留下來的適應性機制[10]。鹽脅迫早期響應階段, 植物可以通過高滲感受器OSCA1以及低滲感受器MSL8感受滲透壓的變化[11-12], 同時植物還可以通過鈉離子感受器MOCA1感受鹽濃度變化[13]。植物感受到鹽脅迫后, 可通過滲透調節、組織特異性離子轉運、激素信號傳導、以及活性氧清除等機制來應對鹽脅迫, 其中也涉及到一系列基因的表達, 包括相關功能基因和調節基因, 如轉錄因子、蛋白激酶等。SOS信號通路是近20年來發現的典型的鹽脅迫響應通路,SOS1、SOS2以及SOS3編碼構成SOS信號通路的蛋白質, 在多種植物的鹽脅迫響應過程中發揮著重要功能[14]。AtNHX1和AtNHX2編碼液泡膜上的Na+/H+逆向轉運蛋白, 主要在根和葉中表達, 并選擇性地將Na+運輸到液泡, 其突變體對鹽脅迫更加敏感[15]。其中,轉錄因子作為關鍵的調控因子, 其介導的基因調控網絡在植物響應鹽脅迫時起著不可或缺的作用。擬南芥AtNIG1是植物中首先被鑒定出來參與耐鹽性的BHLH轉錄因子家族成員, 過表達AtNIG1表現出更高的耐鹽性[16]。玉米轉錄因子ZmWRKY114作為負調控因子, 通過ABA依賴性途徑參與鹽脅迫響應[17]。水稻轉錄因子基因OsDOF15通過介導乙烯生物合成途徑, 在鹽脅迫抑制初生根伸長的過程中起著關鍵作用[18]。棉花轉錄因子基因GhABF3可以通過降低細胞氧化程度顯著提高陸地棉的耐鹽性[19]。大豆轉錄因子基因GmNAC06可以通過控制毛狀根Na+/K+比值維持離子穩態, 從而增強大豆的耐鹽性[20]。這些研究結果表明, 相關轉錄因子的鑒定及功能解析對于理解植物響應鹽脅迫的分子機制至關重要。

隨著高通量測序技術的快速發展, RNA-seq技術已經被廣泛地用于多種植物的研究中, 同時也為轉錄組學研究提供了大量的分析方法。權重基因共表達網絡分析(WGCNA)是最常用的分析RNA-seq數據的方法之一。相比于其他分析方法, WGCNA在產生具有生物學意義的結果方面也更加穩健[21]。該方法可以通過構建基因共表達網絡將上萬個基因劃分成十幾個具有相似表達模式的基因模塊, 并通過分析與性狀高度相關的關鍵模塊從中挖掘核心基因。Zheng等[22]利用324個RNA-Seq數據構建了玉米表皮蠟質合成相關基因的共表達網絡, 發現幾乎所有已知的有關蠟質合成的基因都在同一個模塊, 這為后續深入挖掘玉米表皮蠟質合成相關基因提供了重要的參考依據。李旭凱等[23]利用WGCNA預測出與多種脅迫高度相關的25個抗逆關鍵基因, 為水稻的多抗逆基因挖掘以及利用提供了參考依據。Ye等[24]利用鹽梯度處理的RNA-seq數據, 闡明了互花米草高耐鹽性的轉錄組動力學, 并通過WGCNA分析發現SaOST1等蛋白激酶編碼基因是耐鹽調控網絡中的核心基因。Zhao等[25]以花生為材料, 利用WGCNA篩選出一個與干旱脅迫生理數據密切相關的重要模塊, 并從中預測到9個抗旱候選基因。Yang等[26]基于剪接異構體的WGCNA分析, 篩選出與不同非生物脅迫相關的23個共表達模塊, 并發現許多重要的核心基因在單一應激或多種應激反應中起關鍵作用, 其中主要包括轉錄因子、RNA結合蛋白等。

利用WGCNA構建基因共表達網絡, 已經被廣泛地應用于多種植物的轉錄組研究以及核心基因篩選。本研究通過測得的90份RNA-seq數據, 獲得了甘藍型油菜鹽脅迫處理前后的高分辨率時間動態轉錄表達譜, 進一步通過差異基因分析和WGCNA分別構建油菜根系以及葉片組織響應鹽脅迫的基因共表達網絡, 并從中挖掘油菜響應鹽脅迫的核心基因。研究結果可為甘藍型油菜的鹽脅迫共表達網絡研究以及油菜的耐鹽性分子育種提供理論支持。

1 材料與方法

1.1 植物材料和鹽脅迫處理

選取甘藍型油菜ZS11作為試驗材料, 采用水培法進行培養, 培養條件為: 光照16 h、黑暗8 h, 溫度24℃左右, 相對濕度為75%。培養27 d后, 選取長勢一致的幼苗,使用200 mmol L–1氯化鈉模擬鹽脅迫處理, 分別在鹽脅迫處理后0、0.25、0.5、1、3、6、12、24 h進行取樣, 同時各個時間點均用正常培養條件下的幼苗作為對照組, 取樣時分離根部組織和葉片組織, 立即將其包入錫箔紙并投入–80℃的液氮中儲存, 每個處理3個生物學重復。取樣完畢后送公司進行RNA提取以及轉錄組測序。鹽脅迫處理條件如圖1所示, 其中0 h的樣本組未顯示。本研究所使用的原始數據尚未公開發表, 但單個基因受到鹽脅迫處理后的轉錄表達量已上傳并可通過BnIR[27]網站(http://yanglab.hzau.edu.cn/BnIR/eFP_single_gene)查詢。

圖1 鹽脅迫樣品處理條件Fig. 1 Treatment conditions of salt stress samples

1.2 RNA-seq數據處理與分析

1.2.1 RNA-seq數據的質控與定量 使用FastQC和MultiQC軟件[28]對RNA-seq數據進行質控, 使用Trim Galore過濾掉低質量的reads片段, 獲得clean reads用于后續分析。使用Hisat2[29]將其批量比對到ZS11參考基因組, 統計所有樣本的比對率。使用Samtools[30]將sam文件轉換成bam文件, 使用FeatureCounts[31]對回帖到參考基因上的CDS序列進行計數, 獲得原始counts值。根據counts值以及CDS序列長度, 使用TBtools[32]計算出基因的TPM值。同時使用原始counts值, 利用R包DESeq2[33]進行標準化, 對鹽脅迫根系以及葉片組織在各個時間點的處理組樣本與其對照組樣本進行兩兩差異基因分析,最終獲得組織在各個處理時間點的差異基因以及上下調基因數目。篩選差異基因的閾值設置為|log2(Fold Change)|>1 & FDR<0.05。

1.2.2 利用WGCNA構建基因共表達網絡 WGCNA(weighted gene co-expression network analysis), 即權重基因共表達網絡分析, 可以通過構建共表達網絡尋找協同表達的基因模塊, 探究基因模塊與樣本性狀之間的相關性, 從中挖掘核心基因, 經常適用于非生物逆境脅迫不同時間點應答、病原菌侵染后不同時間點應答等研究方向[34]。本研究考慮到低表達基因一般不具備生物學意義, 將Deseq2分析得到的鹽脅迫組織響應差異基因賦予TPM值,進一步篩選得到至少在3個樣本中TPM值大于1的組織響應差異表達基因用于WGCNA的輸入值。最小模塊的基因數不低于100個, 其余參數按照默認設置。使用R包Cluster Profiler[35]對模塊進行GO富集分析, 閾值設置為FDR<0.05。本研究定義模塊MEs與樣本時間點之間的相關性絕對值大于0.6,P值小于0.05為顯著相關; 以權重值(weight)大于0.2作為基因間的有效連接, 連通度前10%的基因為核心基因。使用Cytoscape軟件對關鍵模塊的共表達網絡進行可視化分析[36]。

2 結果與分析

2.1 鹽脅迫轉錄組數據的產生及初步驗證

本研究于鹽脅迫處理0、0.25、0.5、1、3、6、12、24 h對油菜葉片組織和根系組織進行轉錄組測序, 獲得90份RNA-seq數據。通過質檢和過濾, 共獲得52.78億條clean reads。將其比對到ZS11參考基因組, 發現所有樣本的比對率都超過95%, 表明本研究所獲得的RNA-seq數據具備較高質量。為了驗證鹽脅迫處理的有效性, 利用擬南芥已報道參與鹽脅迫響應的marker基因, 觀察其對應的油菜同源基因在鹽脅迫處理前后的表達模式。SOS1、WRKY33通常在鹽脅迫早期響應階段被激活[37-38], 油菜根系組織中的BnaC09.SOS1以及BnaC04.WRKY33在鹽脅迫1 h內表達量顯著上升(圖2-A, B)。HKT1作為Na+運載體, 可以裝載葉片中的鈉離子到韌皮部再循環到根部[39],油菜葉片組織中的BnaA09.HKT1在鹽脅迫處理3 h時出現顯著差異表達(圖2-C)。ABI1可以負調控ABA促進氣孔關閉從而影響植物耐鹽性[40], 油菜葉片以及根系組織中的BnaA03.ABI1在鹽脅迫處理3 h時表達量顯著上升(圖2-D)。以上結果均與擬南芥中的報道結果類似, 這表明本研究對甘藍型油菜幼苗進行的鹽脅迫處理是有效的,也進一步驗證了RNA-seq數據的可靠性。

圖2 鹽脅迫響應marker基因的表達模式Fig. 2 Relative expression pattern of marker genes under salt stress

2.2 鹽脅迫的相關性分析

為了探究樣本處理之間的相關性, 本研究進行了主坐標分析以及層次聚類分析。從鹽脅迫的主坐標分析結果來看, 第1主坐標主要解釋的是組織表達的差異性, 其貢獻度較大, 約占75.08%, 表明后續在利用WGCNA構建基因表達網絡時, 應將脅迫的根系和葉片組織分開分析,避免組織差異性覆蓋脅迫處理差異性; 第2主坐標主要解釋的是對照組以及鹽脅迫處理之間的差異性, 鹽脅迫處理3 h、6 h、12 h、24 h的樣本表達與對照組形成明顯差異,這也表明脅迫處理的有效性和數據的良好性(圖3-A)。從層次聚類結果來看, 根和葉由于組織表達的差異性較大聚為2類, 這與主坐標分析結果相對一致。另外, 無論是根系組織還是葉片組織, 1 h處理前后的組織樣本均形成了明顯的聚類差異。根據組織響應鹽脅迫的樣本聚類結果進行初步劃分, 可以將鹽脅迫處理1 h前后劃分為早期響應階段以及后期響應階段(圖3-B)。

圖3 鹽脅迫的相關性分析Fig. 3 Correlation analysis of salt stress

2.3 差異基因分析

本研究利用DESeq2包進行差異基因分析, 獲得了油菜根系以及葉片組織在鹽脅迫各個處理時間點的差異基因以及上下調差異基因數目, 并分析了差異基因數目隨鹽脅迫處理時間的變化趨勢。鹽脅迫處理1 h內, 葉片組織在0.25 h、0.5 h和1 h的差異基因數目始終小于根系組織, 表明根系組織在早期階段的響應程度大于葉片組織; 鹽脅迫處理1 h后, 葉片組織在任意一個處理時間點的差異基因數目都大于根系組織, 表明葉片組織在后期階段的響應程度大于根系組織(圖4-A)。其中, 根系組織在鹽脅迫處理3 h時達到最高的響應程度, 此時差異基因為10,883個; 而葉片組織在處理24 h時達到最高的響應程度, 此時差異基因數目為18,414個。進一步分析鹽脅迫上下調差異基因數目的變化趨勢, 發現鹽脅迫處理24 h內, 根系組織中的上調基因數目始終大于下調基因, 而葉片組織僅在處理3 h內的上調基因數目大于下調基因, 而隨著鹽脅迫時間的增長,葉片組織中的上調基因數開始降低并逐漸少于下調基因(圖4-B)。表明油菜幼苗對鹽脅迫的響應程度存在組織差異性。鹽脅迫處理12 h內, 根系組織和葉片組織的響應程度均呈現出先升高后下降的變化趨勢, 變化拐點均出現在處理3 h附近, 表明油菜幼苗對鹽脅迫的響應程度也存在組織相似性。本研究定義至少在任意一個時間點存在顯著差異的基因作為鹽脅迫組織響應差異基因集, 通過去重取并集, 得到根系組織以及葉片組織響應差異基因分別是20,462個和29,334個, 表明全程鹽脅迫處理24 h內, 葉片組織整體上的響應程度比根系更劇烈。

圖4 鹽脅迫處理24 h內的差異基因數目及變化趨勢Fig. 4 Number and trend of differential genes under salt stress treatment for 24 hours

2.4 權重基因共表達網絡的構建

為了構建油菜幼苗響應鹽脅迫的基因共表達網絡,并從中挖掘核心基因, 基于鹽脅迫處理24 h的轉錄組表達矩陣用于WGCNA分析。同時, 為了防止組織差異性過大影響分析結果, 在保證不低于WGCNA所要求的樣本數的情況下, 分別構建了根系組織以及葉片組織響應鹽脅迫的基因共表達網絡。過濾低表達基因, 篩選26,839個葉片組織差異表達基因, 使用動態樹切割方法, 選取最佳軟閾值β=14構建葉片組織響應鹽脅迫的基因共表達網絡(圖5-A), 共識別出14個不同的共表達模塊, 最大和最小模塊分別是turquoise和tan模塊, 分別含有8651個和260個基因(附表1)。過濾低表達基因, 篩選18,952個根系組織響應差異表達基因, 選取最佳軟閾值β=18構建根系組織響應鹽脅迫的基因共表達網絡(圖5-B), 共識別出11個不同的共表達模塊, 最大和最小模塊分別是turquoise以及greenyellow模塊, 分別含有4058個和234個基因(附表2)。

圖5 最佳軟閾值的選擇Fig. 5 Selection of optimal soft threshold (power)

2.5 早期及后期響應模塊的鑒定

在葉片組織響應鹽脅迫的共表達網絡中, 通過計算模塊與樣本處理時間點之間的相關性發現, tan模塊、salmon模塊分別與鹽脅迫早期響應階段0.5 h、1 h顯著正相關; red模塊、cyan模塊、turquiose模塊以及green模塊分別與鹽脅迫后期響應階段3 h、6 h、12 h、24 h顯著正相關(圖6-A)。在根系組織響應鹽脅迫的共表達網絡中,通過計算模塊與樣本處理時間點之間的相關性發現, yellow模塊與鹽脅迫早期響應階段0.25 h以及0.5 h顯著正相關; blue模塊和turquiose模塊分別與鹽脅迫處理1 h、3 h顯著正相關; black模塊與鹽脅迫處理12 h顯著正相關;purple模塊以及red模塊均與鹽脅迫處理24 h顯著正相關(圖6-B)。綜合考慮模塊與樣本處理時間點的相關性、模塊內的轉錄因子所占比例以及模塊的GO富集結果, 選取葉片組織中的tan模塊和根系組織中的yellow模塊作為組織早期響應鹽脅迫的關鍵模塊, 葉片組織中的green模塊和根系組織的red模塊作為組織后期響應鹽脅迫的關鍵模塊。葉片組織中的早期響應tan模塊包含260個基因, 其中有75個轉錄因子, 該模塊內的絕大多數基因在鹽脅迫處理0.5 h的表達量顯著高于其他時間點(圖7-A)。根系組織中的早期響應yellow模塊包含2322個基因, 其中含有276個轉錄因子, 該模塊內的絕大多數基因在鹽脅迫處理0.25 h以及0.5 h顯著表達(圖7-B)。葉片中的后期響應green模塊包含1654個基因, 其中含有138個轉錄因子;根系中的后期響應red模塊包含1238個基因, 含有99個轉錄因子, 2個模塊內的絕大多數基因均在鹽脅迫處理24 h時的表達量顯著高于其他時間點(圖7-C, D)。

圖6 基因模塊與樣本時間點的相關性熱圖Fig. 6 Correlation heat map between gene modules and sample time point

圖7 早期及后期響應模塊的基因表達譜Fig. 7 Gene expression profiles of early and late response modules

2.6 早期及后期響應模塊的GO富集分析

為了分析組織響應鹽脅迫的共表達模塊特征, 本研究分別對早期響應階段的tan模塊和yellow模塊, 后期響應階段的green和red模塊進行GO富集分析。早期響應階段, 根系最先受到鹽脅迫的刺激, 根系組織中的yellow模塊富集到了對生物刺激反應的調控、對外部刺激反應的調控、磷酸化信號轉導、對鹽的響應等生物學過程(圖8-B)。葉片中的tan模塊富集到與葉片衰老、鹽脅迫響應的正向調控、自噬調節等生物學過程(圖8-A)。通過對比分析發現, 早期響應模塊均富集到了脫落酸、乙烯、茉莉酸等激素信號通路, 這可能表明油菜幼苗早期響應鹽脅迫就會有多種激素參與全身器官的系統性反應。早期響應階段,葉片中還富集到了對細菌的防御反應和對昆蟲的防御反應, 根系也富集到了對細菌防御反應的調控以及對細菌來源分子的響應, 這可能表明油菜幼苗早期響應鹽脅迫可能與早期響應細菌等生物脅迫有相似之處。除此之外,葉片中也富集到了抗毒素合成通路, 根系中富集到了類黃酮生物合成過程的調控。這些早期響應階段的生物學通路富集, 可能表明油菜幼苗在早期響應鹽脅迫時也會通過合成抗毒素等物質來響應鹽脅迫。鹽脅迫后期響應階段,green和red模塊顯著富集到了對過氧化氫的響應、細胞對鹽脅迫的反應、對脫水反應的調控等鹽脅迫響應相關通路(圖8-C, D)。其中, 葉片中的green模塊也富集到了與早期響應相似的生物學通路, 如脫落酸激活信號通路;green模塊也特異富集到了有機羥基化合物合成、肌醇生物合成等通路(圖8-C)。根系中的red模塊特異富集到了有機酸轉運、植物型次生細胞壁的發生、苯丙烷生物合成過程、硫代葡萄糖苷代謝過程、木質素代謝過程等通路(圖8-D)。總體來看, 無論是早期響應模塊, 還是后期響應模塊, 都富集到了鹽脅迫響應的相關通路, 表明利用WGCNA可以構建具有生物學意義的共表達模塊, 這些模塊可成為本研究的重點關注對象。

(圖8)

圖8 早期及后期響應模塊的GO富集分析Fig. 8 GO enrichment of early and late response modules

2.7 核心基因的鑒定與分析

為了識別關鍵模塊中的核心基因,利用cytoscape進行共表達網絡可視化。以權重值大于0.2作為閾值, 篩選連通度最大的前10%基因作為核心基因, 早期響應tan模塊篩選到23個核心基因, 包含15個核心轉錄因子(圖9-A)。其中連通度最高的3個核心基因分別為BnaA04G0009600ZS(AGP20)、BnaC04G0304800ZS(CAMBP25)、BnaA10G 0150300ZS(RING1)。BnaA04G0009600ZS(AGP20)編碼調節細胞壁擴張的阿拉伯半乳糖蛋白, 其表達主要由鹽脅迫誘導[41]。早期響應yellow模塊篩選到221個核心基因,包含26個核心轉錄因子 (圖9-B)。其中連通度最高的3個核心基因分別為BnaC09G0584800ZS(NHL3)、BnaC08G 0479500ZS(WAKL8)、BnaA09G0091000ZS(SUPA)。BnaA09G 0091000ZS對應的擬南芥基因SUPA定位于過氧化物酶體中, 鹽脅迫處理15 min內迅速上調表達, 過表達可導致擬南芥中活性氧類物質積累增加以及抗逆性改變[42]。后期響應green模塊篩選到核心基因137個, 包含14個核心轉錄因子(圖9-C)。其中連通度最高的3個核心基因分別為BnaA03G0371900ZS(AOX1A)、BnaC07G0256700ZS(PBL19)、BnaA02G0027300ZS(DGK1)。BnaA03G 0371900ZS對應的擬南芥基因AOX1A編碼線粒體電子傳遞鏈的替代氧化酶, 會限制脯氨酸誘導的氧化應激, 從而利于擬南芥鹽度恢復[43]。后期響應red模塊篩選到核心基因115個, 包含12個核心轉錄因(圖9-D)。其中連通度最高的3個核心基因分別為BnaC05G0203500ZS(PME6)、BnaA02G0180600ZS(WRKY57)、BnaC06G0284600ZS(KTI1)。BnaC05G0203500ZS對應的擬南芥PME6編碼果膠甲酯酶,在保衛細胞中高度表達, 并且是氣孔功能所必需的, 基因恢復可以挽救保衛細胞壁果膠甲酯化狀態、氣孔功能和植物生長[44], 該基因也在油菜鹽脅迫后期顯著差異表達,可能通過氣孔調節應對鹽脅迫。

圖9 早期及后期響應模塊的基因共表達網絡可視化Fig. 9 Co-expression network visualization of early and late response modules

轉錄因子作為關鍵調控因子, 其介導的基因調控網絡經常在植物響應鹽脅迫時發揮著重要作用。本研究著重關注核心轉錄因子, 通過雙向blast將核心基因比對至擬南芥, 并借助TAIR數據庫注釋核心轉錄因子的功能(附表3), 發現油菜同源基因WRKY33、DDF1、AZF2、ARR1、ZFHD1、DREB2B等是擬南芥中已知參與編碼鹽脅迫響應的核心轉錄因子。其中, 擬南芥中的同源基因WRKY33對應油菜中的多個拷貝BnaA04G0247200ZS、BnaC03G 0217300ZS、BnaA03G0185200ZS、BnaC04G0562300ZS, 4個拷貝均在早期響應tan模塊中占據較高連通度, 并參與多種非生物脅迫響應, 尤其是對鹽脅迫的響應[37]。Yellow模塊中的早期響應基因BnaA06G0081500ZS、BnaC05G 0100700ZS在擬南芥中的同源基因DDF1編碼ERF/AP2轉錄因子家族DREB亞家族成員, 過表達可增強對高水平鹽的耐受性[45]。Yellow模塊中的早期響應基因BnaC05G 0398500ZS在擬南芥中的同源基因AZF2編碼鋅指蛋白轉錄因子, 響應ABA、高鹽和輕度脫水反應[46]。Green模塊中的后期響應基因BnaA01G0191600ZS在擬南芥中的同源基因ARR1屬于細胞分裂素信號成分擬南芥響應調節因子1, ARR1蛋白可以通過MPK3/6負調控作用保持穩定從而增強耐鹽性[47]。Red模塊中的后期響應基因BnaA07G 0311300ZS在擬南芥中的同源基因ZFHD1編碼鋅指轉錄因子家族成員, 可由干旱、高鹽和脫落酸誘導表達[48]。無論是早期響應模塊還是后期響應模塊, 均有已知的核心轉錄因子參與擬南芥的鹽脅迫響應, 這表明本研究所篩選的關鍵模塊以及核心轉錄因子具有較高的可靠性。據此初步推測, 這些核心基因可能也在甘藍型油菜的耐鹽性過程中發揮著重要作用。

2.8 核心轉錄因子在505份鹽脅迫油菜群體變異數據中的SNPs及單倍型分析

進一步利用已發表的505份鹽脅迫處理油菜群體變異數據庫, 分析核心轉錄因子在甘藍型油菜群體變異中的SNPs及單倍型情況(http://yanglab.hzau.edu.cn/BnIR/single_locus)[49]。結果表明, 多個核心轉錄因子在505份油菜群體中存在較為豐富的SNPs變異以及單倍型類型。其中,BnaC04G0015700ZS(BnWRKY46)是根系組織早期響應核心基因, 其在鹽脅迫處理0.5 h時可上調表達約6倍。BnaA07G0271900ZS(BnWRKY57)是根系組織后期響應核心基因, 其在鹽脅迫后期處理24 h可達到最大上調倍數約7倍。BnWRKY46和BnEWRKY57均由2個內含子和3個外顯子組成, 經在群體內序列變異分析發現, 這2個核心基因均呈現出較為豐富的SNPs變異(圖10-A, C)。單倍型分析發現,BnWRKY46可初步劃分為3種單倍型, 不同單倍型在鹽脅迫處理后的根長變短, 其中hap.B和hap.C型在鹽脅迫處理前后的根長相對值顯著小于hap.A型(圖10-B)。BnWRKY57可劃分為3種單倍型, 不同單倍型的根長在鹽脅迫處理后顯著縮短, 其中hap.C型在鹽脅迫處理前后的根長相對值顯著大于hap.A和hap.B型(圖10-D)。表明本研究所挖掘到的核心轉錄因子中很有可能存在參與鹽脅迫響應的關鍵候選基因, 深入鑒定鹽脅迫相關的極端單倍型也將有利于為耐鹽型油菜提供優質的種質資源。

圖10 核心基因BnWRKY46和BnWRKY57在505份群體數據中的SNPs及單倍型分析Fig. 10 SNPs and haplotype analysis of hub genes BnWRKY46 and BnWRKY57 in 505 population data

3 討論

3.1 甘藍型油菜響應鹽脅迫的時間動態轉錄表達譜分析

油菜作為世界上重要的油料作物, 在生長過程中經常受到鹽脅迫的影響, 可能會導致減產、品質下降甚至死亡。而甘藍型油菜被認為是中等耐鹽作物, 因此解析甘藍型油菜響應鹽脅迫的調控網絡, 對于油菜的耐鹽性改良以及鹽堿地的開發和利用都具有重要意義。本研究利用轉錄組測序獲得了甘藍型油菜幼苗在鹽脅迫處理前后的90份RNA-seq數據, 將其比對到ZS11參考基因組, 發現所有樣本的比對率都超過95%, 鹽脅迫響應marker基因的表達模式也與擬南芥報道結果類似, 表明鹽脅迫處理是有效的, 數據質量是可靠的。

由于差異基因數目可以在一定程度上反映脅迫響應程度, 本研究發現油菜幼苗對鹽脅迫的響應程度存在組織差異性和相似性。全程鹽脅迫處理24 h內, 葉片組織整體上的響應程度比根系更劇烈。根系組織24 h內響應鹽脅迫的差異基因動態變化趨勢與2007年發表的擬南芥鹽脅迫24 h內的轉錄動態極為類似, 均呈現出先上升后下降的趨勢。并且, 根系組織在24 h內的任意處理時間點的上調基因數目始終大于下調基因數目, 這一點也與擬南芥中的研究結果相對一致[50]。然而, 不同之處在于, 油菜中的組織響應差異基因數目遠高于擬南芥, 這可能是由于異源四倍體的油菜經歷過自然加倍等事件具有更多的同源基因。綜上所述, 本研究獲得了甘藍型油菜響應鹽脅迫的高分辨率時間動態轉錄表達譜, 可為后續研究人員提供可靠的數據參考。

研究人員曾在擬南芥中發現非生物逆境脅迫早期反應與后期反應的轉錄差異, 并認為早期反應可能是非特異性的, 主要是在脅迫處理1 h內, 而基因轉錄的特異性應激反應是在1～3 h后被檢測出來[51]。根據鹽脅迫樣本層次聚類結果, 本研究也發現油菜響應鹽脅迫存在明顯的早期響應和后期響應的聚類差異。對于全程鹽脅迫處理24 h的樣本, 處理1 h前后的根系以及葉片組織均形成了明顯的聚類差異。

3.2 利用WGCNA構建鹽脅迫響應基因共表達網絡

WGCNA分析經常被用于各種植物的轉錄組數據分析, 但多數限制于樣本數目或處理時間點有限而采用多脅迫或多組織混合分析的方式。本研究利用鹽脅迫處理24 h內的多時間點樣本數據, 可以更好地反映基因響應鹽脅迫的表達模式, 從而最大程度地排除其他干擾因素,將不同表達模式的基因歸類到相應模塊。同時, 為了避免組織差異性覆蓋脅迫處理的差異性, 本研究分別利用根系以及葉片組織響應差異表達基因進行WGCNA分析,構建組織響應鹽脅迫的基因共表達網絡。通過分析模塊與樣本時間點的相關性可以獲得組織在不同處理時間點響應鹽脅迫的顯著模塊, 并著重關注早期響應階段的tan模塊和yellow模塊以及后期響應階段的green和red模塊。

已有研究表明, 擬南芥中的自噬水平在鹽脅迫處理30 min內可達到峰值, 并且這些鹽誘導的PCD事件是由Na+特異引起的[52]。本研究在葉片組織中的tan模塊顯著富集到了自噬反應的調控通路, 該結果與擬南芥中的報道相對一致。油菜在鹽脅迫處理0.5 h左右, 體內自噬調節通路就會顯著富集, 可能表明此時油菜幼苗的葉片不單受到鹽處理的滲透作用, 鹽脅迫的離子作用已經可以通過根系傳導到葉片, 并引起葉片中鹽脅迫特異的離子脅迫反應。植物可以通過多種串擾途徑協調各種激素的合成、信號傳導以及新陳代謝來應對鹽脅迫, 從而建立有效的防御系統[53]。本研究在早期響應模塊中顯著富集到了脫落酸、乙烯、茉莉酸、水楊酸等激素信號通路, 這可能表明油菜在鹽脅迫1 h內就會有多種激素信號參與早期脅迫反應的調控。除激活多種激素信號之外, 鹽脅迫早期響應階段的根系和葉片組織中也顯著富集到了植物抗毒素合成、植保素代謝、類黃酮合成以及對細菌的防御反應等通路。已有研究表明, 植物抗毒素主要包括倍半萜、類黃酮和植保素等物質, 其中植保素是擬南芥中一種主要的植物抗毒素, 可以由多種植物病原體誘導產生[54]。這可能表明鹽脅迫早期響應階段可能和病原菌等生物脅迫具有一定的相似性, 同時也會引起機體合成抗毒素等物質進行初級的防御反應。

鹽脅迫造成的生理干旱會引起脫落酸積累從而調控氣孔關閉。本研究在早期響應模塊和后期響應模塊均富集到由不同共表達基因參與的脫落酸信號激活通路。這可能表明鹽脅迫全程響應過程中, 會有不同的共表達基因在不同的響應階段參與調控相似的生物學通路。與早期響應階段相比, 后期響應階段也具有特異性。后期響應階段的green模塊特異富集到了肌醇合成相關通路。已有研究表明, 原本是受到鹽脅迫影響的植物, 在肌醇生產保持不間斷時會表現出正常的生長情況, 推測肌醇代謝途徑的選擇性控制可能是提高耐鹽性的方法之一[55]。后期響應階段的red模塊特異富集到了植物型次生細胞壁的發生以及木質素代謝過程的通路。已有研究表明, 植物可以通過半纖維素和木質素沉積來加強次生細胞壁來增加細胞壁厚度, 從而維持鹽脅迫的穩態[56]。Red模塊也特異富集到了有機酸轉運以及苯丙烷合成通路。苯丙烷代謝作為植物重要的次級代謝途徑之一, 其代謝產物如木質素以及有機酸等, 在調控植物適應性生長的過程中發揮著重要功能[57]。以上結果表明本研究所挑選的關鍵模塊具有較高的可靠性, 同時也說明利用WGCNA可以構建出具有生物學意義的共表達模塊和通路, 這為下一步挖掘核心基因奠定了基礎。

在篩選核心基因時, 本研究發現BnaA04G0247200ZS由于具有較高的連通度而被篩選出來, 其對應擬南芥同源基因WRKY33。已有研究表明, 過表達WRKY33可增強擬南芥的耐鹽性, 該基因可通過控制根細胞凋亡屏障的形成, 賦予根組織的耐鹽能力[58]。油菜中WRKY33可通過增強植保素合成基因的表達增強油菜對菌核病的抗性[59]。本研究分析發現,WRKY33在油菜ZS11中有7個同源拷貝,其中4個拷貝均在鹽脅迫早期響應階段的同一個tan模塊中被篩選到, 這說明WGCNA對于挖掘具有多拷貝的核心基因, 也具有較高的適用性和可靠性, 同時推測這些核心轉錄因子可能也在油菜響應鹽脅迫的過程中發揮著重要功能。進一步結合505份鹽脅迫處理的油菜群體變異數據庫, 發現一些核心轉錄因子在甘藍型油菜群體中存在較為豐富的SNPs變異以及單倍型類型。其中, 早期響應基因BnaA08G0015900ZS在鹽脅迫處理0.5 h時可上調表達約18倍。其在擬南芥中的同源基因ERF8參與ABA以及免疫信號傳導, 并介導植物對病原體的防御反應[60]。單倍型分析發現, 該基因可初步劃分為5種單倍型, 且不同單倍型在鹽脅迫處理前后的根長方面存在極顯著差異。后期響應模塊中的核心轉錄因子基因BnaA05G0156300ZS(NFYA5)在3 h前與對照組基本無差異, 3 h表達量顯著上升并在24 h達到最大差異倍數。該基因有3種單倍型, 不同單倍型在幼苗期的植株總干重比(鹽脅迫/對照組)方面存在極顯著差異。如BnaA07G0311300ZS(ZFHD1)也是后期響應基因, 有2種單倍型, 不同單倍型在幼苗期的植株總干重比(鹽脅迫/對照組)方面存在極顯著差異。如BnaC01G0510700ZS(NAC047), 有2種單倍型, 不同單倍型在幼苗期的根長比(鹽脅迫/對照組)方面存在極顯著差異。這些基因在模式植物擬南芥中有所報道, 但在油菜中報道還相對較少, 可為后續甘藍型油菜耐鹽性改良提供重要的候選基因資源。

附表請見網絡版: 1) 本刊網站http://zwxb.chinacrops.org/; 2) 中國知網http://www.cnki.net/; 3) 萬方數據http://c.wanfangdata.com.cn/Periodical-zuowxb.aspx。