張 毅，吳萬億，劉霞宇，張 潔，唐銳敏，賈小云

(1.山西農(nóng)業(yè)大學(xué) 生命科學(xué)學(xué)院，山西 太谷 030801； 2.山西農(nóng)業(yè)大學(xué) 農(nóng)學(xué)院，山西 太谷 030801)

隨著全球人口的快速增加，人類對(duì)糧食的需求量日益提高，但全球土地資源有限，并且約有20%的土地正遭受鹽堿的危害。土地鹽堿化是一個(gè)全球性問題，造成大片耕地不能正常使用，據(jù)估計(jì)，到2050年，將有50%的耕地鹽堿化[1]。在鹽堿化土地上種植的植物會(huì)遭受離子毒害、滲透脅迫和氧化脅迫等，進(jìn)而抑制植物生長(zhǎng)甚至產(chǎn)量[2]。甘薯[Ipomoeabatatas(L.)Lam]是重要的糧食作物，其塊根富含淀粉、膳食纖維、維生素和花青素等，具有抗癌和保健作用[3-4]。甘薯對(duì)環(huán)境有較好的適應(yīng)快，但鹽堿脅迫仍然是影響甘薯正常生長(zhǎng)的主要因素。因此，挖掘甘薯鹽脅迫響應(yīng)基因?qū)榻馕龈适磉m應(yīng)鹽環(huán)境的機(jī)制和培育甘薯耐鹽新品種提供重要的理論依據(jù)。

隨著高通量測(cè)序技術(shù)的發(fā)展和測(cè)序價(jià)格的降低，海量的轉(zhuǎn)錄組數(shù)據(jù)產(chǎn)生[5]。傳統(tǒng)的2個(gè)樣本間轉(zhuǎn)錄組數(shù)據(jù)的差異分析已經(jīng)滿足不了需求，而且傳統(tǒng)的轉(zhuǎn)錄組數(shù)據(jù)差異分析只能對(duì)有功能注釋的基因進(jìn)行分析。加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析(Weighted gene co-expression network analysis，WGCNA)可以利用多個(gè)轉(zhuǎn)錄組數(shù)據(jù)構(gòu)建基因共表達(dá)網(wǎng)絡(luò)，通過軟閾值來鑒定表達(dá)模式相似的基因集合模塊(Module)，并特異性篩選出與性狀相關(guān)的共表達(dá)模塊，進(jìn)而挖掘出共表達(dá)模塊中的關(guān)鍵基因。表達(dá)模式相似的基因通常被認(rèn)為可能參與相似的生物學(xué)過程，WGCNA是研究表達(dá)模式相似基因集的重要工具[6]。WGCNA通過構(gòu)建無尺度分布的網(wǎng)絡(luò)，可以對(duì)參與生物學(xué)過程的沒有功能注釋的新基因進(jìn)行預(yù)測(cè)，并結(jié)合已有的基因功能注釋發(fā)現(xiàn)參與特定代謝途徑的共表達(dá)模塊、預(yù)測(cè)新基因的功能。WGCNA在植物學(xué)研究中廣泛應(yīng)用。前人利用WGCNA分析水稻不同時(shí)間點(diǎn)的鎘脅迫轉(zhuǎn)錄組數(shù)據(jù)，鑒定到了271個(gè)差異表達(dá)基因參與鎘的調(diào)控[7]；利用WGCNA分析了紅皮甘薯及其突變體，鑒定出2個(gè)調(diào)節(jié)甘薯皮色的基因模塊[8]；利用WGCNA分析了紅皮蘋果及其突變體黃皮蘋果的轉(zhuǎn)錄組數(shù)據(jù)，發(fā)現(xiàn)1個(gè)與花青素含量高度相關(guān)的共表達(dá)模塊，同時(shí)發(fā)現(xiàn)MdMYB10基因的啟動(dòng)子甲基化是導(dǎo)致果皮顏色差異的原因[9]。但是目前尚未見關(guān)于利用WGCNA鑒定甘薯耐鹽相關(guān)共表達(dá)網(wǎng)絡(luò)及核心基因的報(bào)道。為此，利用甘薯鹽脅迫轉(zhuǎn)錄組數(shù)據(jù)結(jié)合WGCNA方法，鑒定甘薯耐鹽相關(guān)共表達(dá)網(wǎng)絡(luò)及核心基因，為進(jìn)一步研究甘薯耐鹽分子機(jī)制及耐鹽甘薯新品種的培育奠定基礎(chǔ)。

1 材料和方法

1.1 數(shù)據(jù)獲取與分析

甘薯鹽脅迫轉(zhuǎn)錄組數(shù)據(jù)(SRP092215)是從EBI(European Bioinformatics Institute)的ENA(European Nucleotide Archive)數(shù)據(jù)庫獲得的[10]。該轉(zhuǎn)錄組數(shù)據(jù)是以栗子香和ND98兩個(gè)甘薯品種為材料，對(duì)生長(zhǎng)30 d的甘薯幼苗利用 200 mmol/L NaCl溶液進(jìn)行鹽脅迫處理，分別在處理后0、12、48 h取根系進(jìn)行轉(zhuǎn)錄組測(cè)序，每個(gè)處理2個(gè)重復(fù)，共12個(gè)樣本[11]。首先用Aspera軟件下載甘薯鹽脅迫轉(zhuǎn)錄組(SRP092215)數(shù)據(jù)的fastq文件，然后利用FastQC[12]軟件對(duì)原始數(shù)據(jù)進(jìn)行質(zhì)量評(píng)估處理，利用fastp軟件[13]去除接頭和低質(zhì)量堿基，利用Hisat2軟件[14]將過濾后的clean data與甘薯參考基因組[15](https://ipomoea-genome.org)比對(duì)，然后利用featureCount軟件[16]計(jì)算每個(gè)基因reads數(shù)，最后使用TPM(Transcripts per kilobase of exon model per million mapped reads)體現(xiàn)基因表達(dá)量。

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

利用R 3.61軟件中WGCNA 1.69軟件包[17]構(gòu)建加權(quán)基因共表達(dá)網(wǎng)絡(luò)，以標(biāo)準(zhǔn)化后的12個(gè)鹽脅迫轉(zhuǎn)錄組數(shù)據(jù)中的基因表達(dá)矩陣作為輸入數(shù)據(jù)，使用R軟件中g(shù)enefiliter 軟件包篩選在各個(gè)樣本間表達(dá)量變異最大的前50%基因構(gòu)建加權(quán)共表達(dá)網(wǎng)絡(luò)，利用WGCNA 1.69軟件中的pickSoftThreshold 函數(shù)計(jì)算共表達(dá)網(wǎng)絡(luò)的軟閾值，并利用powerEstimate函數(shù)估計(jì)最優(yōu)power值，最終選擇power=16對(duì)原始有尺度關(guān)系矩陣進(jìn)行冪處理，得到無尺度化鄰接矩陣。為更好評(píng)估基因間表達(dá)模式的相關(guān)性, 進(jìn)一步將鄰接矩陣轉(zhuǎn)化為拓?fù)渲丿B矩陣(Topological overlap matrix,TOM)，并利用拓?fù)湎喈惥仃?dissTOM=1-TOM)，采用動(dòng)態(tài)剪切算法進(jìn)行基因聚類及模塊劃分。使用WGCNA自動(dòng)網(wǎng)絡(luò)構(gòu)建函數(shù) blockwiseModules 構(gòu)建網(wǎng)絡(luò),其中，最小模塊基因數(shù)為30，相似模塊的合并閾值為0.25(mergeCutHeight=0.25)，其他參數(shù)按照默認(rèn)設(shè)置。

1.3 鹽脅迫特異性模塊的鑒定

對(duì)每個(gè)共表達(dá)模塊中的所有基因進(jìn)行主成分分析(Principle component analysis,PCA)，將主成分1(PC1)稱為此模塊的特征向量(Module eigengene,ME)，為了篩選鹽脅迫相關(guān)特異性模塊，計(jì)算各個(gè)模塊的ME值與不同性狀(此處為不同鹽脅迫時(shí)間)的相關(guān)系數(shù)r和P值。r>0表示正相關(guān)，r<0表示負(fù)相關(guān)。在本研究中，若模塊ME值和性狀相關(guān)系數(shù)|r|>0.70且P<0.05，則認(rèn)為該模塊是特異性模塊。

1.4 GO功能富集分析

利用TBtools (A tool kit for biologists integrating various biological data handling tools with a user-friendly interface)軟件[18]對(duì)每個(gè)共表達(dá)模塊中的基因進(jìn)行GO(Gene ontology)富集分析，閾值為P<0.05。

1.5 耐鹽特異性模塊核心基因的鑒定及基因互作網(wǎng)絡(luò)構(gòu)建

將共表達(dá)模塊中每個(gè)基因的表達(dá)量與模塊ME值進(jìn)行關(guān)聯(lián)，得到模塊關(guān)系值(Module membership，MM)，其可以體現(xiàn)基因在模塊中的重要性，MM值越高，基因在模塊中越重要。選取耐鹽正相關(guān)特異性模塊中MM值排名前5的基因?yàn)槟K中的核心基因，并利用Cytoscape軟件展示核心基因互作網(wǎng)絡(luò)。

1.6 轉(zhuǎn)錄因子分析

將各個(gè)模塊中的蛋白質(zhì)序列提交到plantTFDB數(shù)據(jù)庫[19](http://planttfdb.gao-lab.org/)進(jìn)行轉(zhuǎn)錄因子分析預(yù)測(cè)，獲得各個(gè)模塊中預(yù)測(cè)的轉(zhuǎn)錄因子。

2 結(jié)果與分析

2.1 甘薯加權(quán)共表達(dá)網(wǎng)絡(luò)的構(gòu)建

本研究過濾掉變異較小的基因，最終選擇32 147個(gè)基因用于構(gòu)建加權(quán)基因共表達(dá)網(wǎng)絡(luò)。使用WGCNA軟件包中powerEstimate函數(shù)估計(jì)最適power值為16，當(dāng)power值為16時(shí)，無尺度網(wǎng)絡(luò)擬合指數(shù)R2>0.85，平均連通性趨于0(圖1)。該power值可以獲得符合要求的無尺度網(wǎng)絡(luò)。采用動(dòng)態(tài)剪切法劃分模塊，合并相似度大于75%的模塊，最終構(gòu)建了20個(gè)共表達(dá)模塊，不同顏色表示不同模塊(圖2)。各模塊含基因數(shù)目如圖3，其中，Turquoise模塊所含基因數(shù)最多，為10 383個(gè)；Royalblue模塊中基因數(shù)最少，為124;其他模塊的基因數(shù)介于242～6 745；Grey模塊中的基因是一組沒有分配到其他模塊中的基因集。

a.不同軟閾值對(duì)應(yīng)的無尺度網(wǎng)絡(luò)擬合指數(shù)R2，直線代表R2=0.85; b.不同軟閾值對(duì)應(yīng)的平均連通性

上部為基因聚類樹，下部是按樹的分枝切割得到的模塊，相同的模塊用同一種顏色表示

圖3 模塊中基因數(shù)目分布

2.2 甘薯鹽脅迫相關(guān)特異性模塊的鑒定

由圖4可知，本研究共鑒定到6個(gè)與鹽脅迫特異相關(guān)的基因共表達(dá)模塊(|r|>0.70，P<0.05)。其中，Green模塊(r=0.71，P=0.01)和Red模塊(r=0.85，P=5e-04)與鹽脅迫48 h正相關(guān)；Black模塊(r=0.83，P=9e-04)和Yellow模塊(r=0.79，P=0.002)與鹽脅迫12 h正相關(guān)；Midnightblue模塊(r=-0.74,P=0.006)和Magenta模塊(r=-0.73,P=0.007)與鹽脅迫48 h負(fù)相關(guān)。為了驗(yàn)證模塊劃分的可靠性，進(jìn)一步選擇與鹽脅迫正相關(guān)的4個(gè)特異性模塊(Green模塊、Black模塊、Yellow模塊和Red模塊)作為目標(biāo)模塊做深入分析。

括號(hào)外數(shù)據(jù)、括號(hào)內(nèi)數(shù)據(jù)分別表示模塊和性狀之間的相關(guān)系數(shù)、P值

2.3 甘薯鹽脅迫正相關(guān)特異性模塊GO富集分析

由圖5可知，選取4個(gè)與甘薯鹽脅迫正相關(guān)的特異性模塊(Green模塊、Black模塊、Yellow模塊和Red模塊)進(jìn)行GO富集分析，發(fā)現(xiàn)其均可以富集到生物學(xué)過程、分子功能和細(xì)胞組分。Green模塊富集到139個(gè)生物學(xué)過程、26個(gè)細(xì)胞組分、41個(gè)分子功能，主要富集到干旱響應(yīng)(GO:0009414)、葉綠體(GO:0009507)、非生物刺激的反應(yīng)(GO:0009628)、高滲響應(yīng)(GO:0006972)等；Black模塊富集到84個(gè)生物學(xué)過程、51個(gè)細(xì)胞組分、24個(gè)分子功能，主要富集到跨膜運(yùn)輸相關(guān)通路，如活性離子跨膜轉(zhuǎn)運(yùn)活性(GO:0022853)、質(zhì)子跨膜運(yùn)輸(GO:1902600)、跨膜運(yùn)輸(GO:0055085 )等；Yellow模塊富集到265個(gè)生物學(xué)過程、24個(gè)細(xì)胞組分、53個(gè)分子功能，主要富集到信號(hào)傳遞(GO:0023052)、質(zhì)膜(GO:0005886)、脫落酸激活的信號(hào)通路(GO:0009738)、離子結(jié)合(GO:0043167)等；Red模塊富集到63個(gè)生物學(xué)過程、25個(gè)細(xì)胞組分、22個(gè)分子功能，主要富集到抗氧化活性(GO:0010035)、光合作用膜(GO:0034357)、水楊酸刺激的細(xì)胞反應(yīng)(GO:0071446)等。

圖5 甘薯耐鹽相關(guān)特異性模塊GO富集分析

2.4 甘薯鹽脅迫顯著正相關(guān)共表達(dá)模塊中核心基因的鑒定及基因互作網(wǎng)絡(luò)構(gòu)建

本研究選取與甘薯鹽脅迫顯著正相關(guān)的共表達(dá)模塊中MM 值排名前5的核心基因，并將其與擬南芥蛋白質(zhì)序列進(jìn)行比對(duì)，預(yù)測(cè)核心基因的功能(表1)。核心基因主要有HAK5(High affinity K+transporter 5)、TLP5(Tubby like protein 5)、NAGS2(N-acetyl-L-glutamate synthase 2)、bHLH115(Basic helix-loop-helix 115)、SK11(Shaggy-related kinase 11)、ABCG28(ABC transporter G family member 28)、DUF699(GNAT acetyltransferase)等，其中部分核心基因功能已經(jīng)明確，如Green模塊中核心基因bHLH115是一個(gè)參與植物缺鐵響應(yīng)的正調(diào)控轉(zhuǎn)錄因子基因，與參與響應(yīng)植物逆境脅迫的基因共表達(dá);Yellow模塊中的核心基因SK11可以提高植物的耐鹽性；Red模塊中的核心基因ABCG28有物質(zhì)轉(zhuǎn)運(yùn)功能。進(jìn)一步探尋與核心基因互作關(guān)系強(qiáng)的基因，利用Cytoscape軟件展示互作網(wǎng)絡(luò)(圖6)，發(fā)現(xiàn)與核心基因互作的基因也有相應(yīng)的功能并且已經(jīng)被報(bào)道，如Green模塊中的核心基因TLP5(g42340)與NAC[NAM(no apical meristem)，ATAF1(Arabidopsistranscription activation factor 1)，ATAF2，CUC2(cup-shaped cotyledon 2)]基因家族NAC2(g28394)基因互作，核心基因g48306與HsfB2a(Heat shock transcription factor B2a)(g20568)互作；Black模塊中的核心基因g14417與FAR1(Far-red-impaired response 1)(g54849)基因存在互作；Yellow模塊中核心基因g23125與DREB2C(Dehydration responsive element binding protein 2C)(g1860)基因互作;Red模塊中核心基因g11018與ERF1(Ethylene response factor 1)(g46237)基因互作。

表1 甘薯鹽脅迫顯著正相關(guān)共表達(dá)模塊中的核心基因及其功能

灰色圓形代表核心基因，三角形代表轉(zhuǎn)錄因子，灰色三角形代表既是核心基因又是轉(zhuǎn)錄因子，圓形代表其他基因

3 結(jié)論與討論

植物在長(zhǎng)期適應(yīng)鹽脅迫的過程中形成了應(yīng)對(duì)鹽脅迫的分子機(jī)制。本研究利用甘薯鹽脅迫的轉(zhuǎn)錄組數(shù)據(jù)，構(gòu)建了加權(quán)基因共表達(dá)網(wǎng)絡(luò)，挖掘出與鹽脅迫顯著相關(guān)的6個(gè)特異性模塊，并對(duì)其中的4個(gè)與鹽脅迫顯著正相關(guān)的基因模塊進(jìn)行GO富集分析，發(fā)現(xiàn)其可以顯著富集到與耐鹽相關(guān)的具有生物學(xué)意義的條目，主要包括跨膜運(yùn)輸、水楊酸刺激的細(xì)胞反應(yīng)、非生物刺激的反應(yīng)及高滲響應(yīng)等條目。其中，Green模塊可以顯著富集到非生物刺激的反應(yīng)(GO:0009628)和干旱響應(yīng)(GO:0009414)等GO條目。研究表明，植物在鹽脅迫過程中由于滲透脅迫引起細(xì)胞水勢(shì)過低，植物無法吸收足夠的水分滿足其正常生命活動(dòng)，進(jìn)而引起植物生理干旱[20]。Black模塊可以顯著富集到活性離子跨膜轉(zhuǎn)運(yùn)活性(GO:0022853)和質(zhì)子跨膜運(yùn)輸(GO:1902600)等相關(guān)通路。研究表明，植物在響應(yīng)鹽脅迫的過程中進(jìn)行滲透調(diào)節(jié)，植物生物膜上的質(zhì)子泵通過水解ATP(5′-adenylate triphosphate)產(chǎn)生能量，定向運(yùn)輸氫離子形成跨膜電勢(shì)梯度，從而可以使得鹽離子進(jìn)行跨膜運(yùn)輸[21]。植物的質(zhì)子泵與離子通道等決定了其對(duì)離子的選擇性吸收。植物在鹽脅迫下，質(zhì)子泵提供能量，把鈉離子排出細(xì)胞外，使得植物體內(nèi)鈉離子含量減少，進(jìn)而使細(xì)胞處于穩(wěn)態(tài)和代謝正常[22]。Red模塊可以顯著富集到水楊酸刺激的細(xì)胞反應(yīng)(GO:0071446)、抗氧化活性(GO:0010035)等條目。Yellow模塊可以顯著富集到脫落酸激活的信號(hào)通路(GO:0009738)。研究表明，脫落酸、茉莉酸、水楊酸等植物激素在植物響應(yīng)鹽脅迫中起著重要的作用，其中脫落酸最關(guān)鍵，其可以參與調(diào)控耐鹽相關(guān)基因的表達(dá)、氣孔關(guān)閉和離子平衡等[23]。

本研究選擇模塊中MM值排名前5的基因作為核心基因進(jìn)行研究，模塊中核心基因?yàn)镠AK5、TLP5、NAGS2、bHLH115、TBL19、SK11、DUF699等。由于核心基因在網(wǎng)絡(luò)中具有較高的連通性，所以其可能參與重要的生物學(xué)過程[24]。通過對(duì)目標(biāo)模塊中的核心基因進(jìn)行功能分析，結(jié)果表明，其主要功能為跨膜運(yùn)輸、蛋白質(zhì)磷酸化、轉(zhuǎn)錄調(diào)控等，然而發(fā)現(xiàn)大部分核心基因的功能在甘薯中鮮有報(bào)道，但是它們?cè)跀M南芥、番茄等植物中的同源基因已經(jīng)被報(bào)道過可響應(yīng)鹽脅迫。如Green模塊中的核心基因g34447在擬南芥中的同源基因bHLH115在維持體內(nèi)鐵穩(wěn)態(tài)中起著關(guān)鍵作用[25],g1032在番茄中的同源基因NAGS1可提高擬南芥的耐鹽性和抗旱性[26]；Black模塊中的核心基因g57672在擬南芥中的同源基因HAK5，屬于高親和K+轉(zhuǎn)運(yùn)體HAK，可以調(diào)節(jié)K+離子濃度，進(jìn)而提高植物的耐鹽性[27]；Yellow模塊中的核心基因g23125在擬南芥中的同源基因SK11在鹽脅迫下被激活，隨后又激活G6PD基因(Glucose-6-phosphate dehydrogenase)，提高了耐鹽性[28]；Red模塊中的核心基因g33775在擬南芥中的同源基因SRO5與P5CDH(Pyrroline-5-carboxylate dehydrogenase)共同控制活性氧(ROS)產(chǎn)生和響應(yīng)脅迫反應(yīng)[29]。植物耐鹽性不是由單個(gè)基因決定的，而是多個(gè)抗逆基因共同作用的結(jié)果。與核心基因互作的基因也可能參與重要的生物學(xué)過程。本研究通過對(duì)與核心基因互作的基因進(jìn)行分析，發(fā)現(xiàn)部分基因功能已經(jīng)被報(bào)道。如與Green模塊中核心基因g42340互作的g28394基因在擬南芥中的同源基因NAC2響應(yīng)多種非生物脅迫和灰霉病菌侵染[30]。與Green模塊中核心基因g48306互作的g20568基因在辣椒中的同源基因HsfB2a參與辣椒對(duì)青枯病菌的免疫和對(duì)高溫高濕的耐受性[31]。與Black模塊中核心基因g14417互作的g54849基因在擬南芥中的同源基因FAR1通過整合葉綠素生物合成和水楊酸信號(hào)通路，在調(diào)節(jié)植物免疫中發(fā)揮重要作用[32]。與Yellow模塊中核心基因g23125互作的g1860基因在擬南芥中的同源基因DREB2C在氧化應(yīng)激耐受性方面發(fā)揮重要作用[33]。與Yellow模塊中核心基因g56260互作的g58671基因在番茄中的同源基因ZFP1(Zinc-finger protein 1)可以提高番茄幼苗對(duì)非生物脅迫的耐受能力[34]。與Red模塊中核心基因g11018互作的g46237基因在擬南芥中的同源基因ERF1通過整合茉莉酸、乙烯和脫落酸信號(hào)調(diào)控脅迫相關(guān)基因，在高鹽、干旱和高溫脅迫中發(fā)揮積極作用[35]。

本研究側(cè)重分析了與鹽脅迫顯著正相關(guān)的模塊，其他模塊中也可能有一些重要基因參與甘薯的耐鹽脅迫。本研究所分析的核心基因在甘薯中的功能尚未明確，后續(xù)需要通過分子生物學(xué)試驗(yàn)進(jìn)一步驗(yàn)證。