甘藍(lán)型油菜鹽脅迫響應(yīng)基因表達(dá)譜分析及共表達(dá)網(wǎng)絡(luò)的構(gòu)建

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

華中農(nóng)業(yè)大學(xué)作物遺傳改良全國重點(diǎn)實(shí)驗(yàn)室 / 國家油菜工程技術(shù)研究中心 / 洪山實(shí)驗(yàn)室, 湖北武漢 430070

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

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

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

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

利用WGCNA構(gòu)建基因共表達(dá)網(wǎng)絡(luò), 已經(jīng)被廣泛地應(yīng)用于多種植物的轉(zhuǎn)錄組研究以及核心基因篩選。本研究通過測(cè)得的90份RNA-seq數(shù)據(jù), 獲得了甘藍(lán)型油菜鹽脅迫處理前后的高分辨率時(shí)間動(dòng)態(tài)轉(zhuǎn)錄表達(dá)譜, 進(jìn)一步通過差異基因分析和WGCNA分別構(gòu)建油菜根系以及葉片組織響應(yīng)鹽脅迫的基因共表達(dá)網(wǎng)絡(luò), 并從中挖掘油菜響應(yīng)鹽脅迫的核心基因。研究結(jié)果可為甘藍(lán)型油菜的鹽脅迫共表達(dá)網(wǎng)絡(luò)研究以及油菜的耐鹽性分子育種提供理論支持。

1 材料與方法

1.1 植物材料和鹽脅迫處理

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

圖1 鹽脅迫樣品處理?xiàng)l件Fig. 1 Treatment conditions of salt stress samples

1.2 RNA-seq數(shù)據(jù)處理與分析

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

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

2 結(jié)果與分析

2.1 鹽脅迫轉(zhuǎn)錄組數(shù)據(jù)的產(chǎn)生及初步驗(yàn)證

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

圖2 鹽脅迫響應(yīng)marker基因的表達(dá)模式Fig. 2 Relative expression pattern of marker genes under salt stress

2.2 鹽脅迫的相關(guān)性分析

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

圖3 鹽脅迫的相關(guān)性分析Fig. 3 Correlation analysis of salt stress

2.3 差異基因分析

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

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

2.4 權(quán)重基因共表達(dá)網(wǎng)絡(luò)的構(gòu)建

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

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

2.5 早期及后期響應(yīng)模塊的鑒定

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

圖6 基因模塊與樣本時(shí)間點(diǎn)的相關(guān)性熱圖Fig. 6 Correlation heat map between gene modules and sample time point

圖7 早期及后期響應(yīng)模塊的基因表達(dá)譜Fig. 7 Gene expression profiles of early and late response modules

2.6 早期及后期響應(yīng)模塊的GO富集分析

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

(圖8)

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

2.7 核心基因的鑒定與分析

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

圖9 早期及后期響應(yīng)模塊的基因共表達(dá)網(wǎng)絡(luò)可視化Fig. 9 Co-expression network visualization of early and late response modules

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

2.8 核心轉(zhuǎn)錄因子在505份鹽脅迫油菜群體變異數(shù)據(jù)中的SNPs及單倍型分析

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

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

3 討論

3.1 甘藍(lán)型油菜響應(yīng)鹽脅迫的時(shí)間動(dòng)態(tài)轉(zhuǎn)錄表達(dá)譜分析

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

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

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

3.2 利用WGCNA構(gòu)建鹽脅迫響應(yīng)基因共表達(dá)網(wǎng)絡(luò)

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

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

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

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

附表請(qǐng)見網(wǎng)絡(luò)版: 1) 本刊網(wǎng)站http://zwxb.chinacrops.org/; 2) 中國知網(wǎng)http://www.cnki.net/; 3) 萬方數(shù)據(jù)http://c.wanfangdata.com.cn/Periodical-zuowxb.aspx。