李俊良

(哈爾濱工業(yè)大學)

0 前言

植物為適應環(huán)境的變化需要不斷地調整自己的生長及生命進程來應對這種變化。這種適應過程受到基因表達的動態(tài)調控，同時相關調控基因的表達水平通常也會隨著環(huán)境的變化而實時改變。與生理水平的研究進展相比，作物在基因水平上對鹽脅迫應答的分子機制研究可謂是突飛猛進[1]。表達譜測序是一種可以進行轉錄鑒定和表達定量的高通量測序方法，廣泛用于研究基因在轉錄水平的表達差異。本研究使用表達譜測序技術并結合WGCNA方法從轉錄水平分析甜菜對鹽脅迫應答的分子機理。

1 材料與方法

1.1 生物材料培育

實驗選用了本課題組選育的適應性強/品質性狀優(yōu)良特別是耐鹽性較強的栽培甜菜“O68”做為實驗材料。將催芽后的甜菜種子定植于定植棉中，定植棉放置在育苗盤上，定期加入適量水，待子葉完全展開后將定植棉插入定植籃后移入培養(yǎng)槽，培養(yǎng)槽內裝入以1/2濃度的改良的霍格蘭營養(yǎng)液為培養(yǎng)基質，水培法培育甜菜幼苗至三對真葉期。培育條件設置為光照(16 h；24℃)/黑暗(8 h；18℃)。培養(yǎng)基質更換周期為3天。

1.2 鹽脅迫處理與樣品收集

選取長勢一致的3對真葉期幼苗分為5組：(1)使用正常培養(yǎng)基質繼續(xù)培育72 h后取樣作為空白對照；(2)使用正常培養(yǎng)基質培育60 h后移入添加了300 mM NaCl的培養(yǎng)基質中繼續(xù)培育12 h后取樣作為鹽處理12 h組；(3)使用正常培養(yǎng)基質培育48 h后移入添加了300 mM NaCl的培養(yǎng)基質中繼續(xù)培育24 h后取樣作為鹽處理24 h組；(4)使用正常培養(yǎng)基質培育24 h后移入添加了300 mM NaCl的培養(yǎng)基質中繼續(xù)培育48 h后取樣作為鹽處理48 h組；(5)直接使用添加了300 mM NaCl的培養(yǎng)基質中培育72 h后取樣作為鹽處理72 h組。表達譜測序使用1-5組的葉片和根共10組樣品。全轉錄組測序和small RNA測序使用Ⅰ和Ⅱ的葉片和根作為樣品。將1-5組的葉片和根樣品混合后分別執(zhí)行葉片和根的降解組測序實驗。代謝組學使用1和2 的葉片樣品。

葉片的取樣部位為第3對真葉，根的取樣部位為自根尖向上起約6 cm長的區(qū)域。生理指標測定所需樣本經(jīng)PBS漂洗后直接使用；高通量測序所需樣本經(jīng)PBS漂洗后使用鋁箔包裹，置于液氮中冷凍30分鐘后保存于-80℃待用。

1.3 表達譜測序挖掘甜菜鹽脅迫應答基因

1.3.1 RNA的提取

稱取0.5 g葉片或根樣本，利用TRIZOL法(Invitrogen, CA, USA)提取每個樣品的總RNA。使用Bioanalyzer 2100和RNA 6000 Nano LabChip檢測RNA質量，以RIN值>7.0作為合格標準。為盡可能排除植株個體差異對結果的影響，制作高質量Total RNA預混池(分別提取6株相同處理葉片或根的RNA，質檢合格后進行混合)用于高通量測序文庫的構建。

1.3.2 表達譜測序文庫構建與測序

對于RNA-seq，取1 μg總RNA使用Poly-Toligo磁珠(Invitrogen)分離Poly(A)mRNA。純化后，使用二價陽離子(在高溫下)將mRNA裂解成小片段。然后依照mRNA Seq sample preparation kit(Illumina, San Diego，USA)方案，將打碎后的RNA片段進行反轉錄構建最終的cDNA文庫。文庫的平均插入長度為300 bp(±50 bp)。按前述構建10組RNA-seq文庫，每組設置3個生物學重復，總計30個cDNA文庫。然后在聯(lián)川生物的Illumina NovaseqTM6000平臺上進行雙端測序。

1.3.3 鹽脅迫應答表達基因的鑒定與表達模式分析

首先使用Cutadapt[2]去除RNA-Seq數(shù)據(jù)原始數(shù)據(jù)中的的引物/適配器污染、低質量堿基和未確定堿基。然后利用FastQC(http://www.bioinformatics.babraham.ac.uk/projects/fastqc/)對序列質量進行驗證。使用Bowtie 2[3]和Hierarchical Indexing for Spliced Alignment of Transcripts (HISAT)[4]將clean reads映射到甜菜參考基因組RefBeet-1.2上(http://bvseq.boku.ac.at/Genome/Download/RefBeet-1.2/)。使用StringTie[5]計算各樣本的映射讀數(shù)。在最終轉錄組生成后利用StringTie和edgeR[6]計算樣本中各轉錄本的表達量。FPKM(fragments per kilobase of transcript per million mapped reads)值被用來來衡量基因的表達量。使用edgeR進行差異表達分析；使用T-test檢驗評估差異基因的顯著性；差異變動基因(differentially expression gene，DEG)的篩選閾值為p<0.05且|log2(fold change)|大于1。組間基因表達模式差異使用UpSetR[7]進行可視化展示。

1.3.4 差異表達基因的功能分析

將所有表達基因與NCBI非冗余蛋白數(shù)據(jù)庫(Nr)進行Blastx比對，E值篩選閾值設置為10-5。使用GO(Gene Ontology)數(shù)據(jù)庫[8]對所有表達基因進行GO注釋；計算每個term的基因數(shù)量，根據(jù)超幾何分布檢驗結果對各term進行富集分析。使用KEGG(Kyoto Encyclopedia of Genes and Genomes)[9]數(shù)據(jù)庫對所有基因進行注釋，根據(jù)超幾何分布檢驗結果對各pathway進行富集分析。使用STRING數(shù)據(jù)庫分析DEG編碼蛋白間的互作關系[10]。

2 結果與分析

2.1 鹽處理對甜菜表型的影響

在甜菜受到脅迫后通常會出現(xiàn)一段生長滯后，在其適應脅迫后生長逐漸恢復，但是這種恢復的生長活力往往低于正常條件[11]。甜菜幼苗在300 mM NaCl處理后其表型如圖1所示：鹽處理下葉片首先失水萎蔫(12 h)；而處理時間的延長至(24 h)時幼葉(第三對真葉)逐漸從失水狀態(tài)恢復正常；當處理時間的繼續(xù)延長至(72 h)時，第一對真葉表現(xiàn)出部分失綠(圖1a)。鹽處理下根系隨著處理時間的延長顏色逐漸加深，24 h時可以明顯觀察到根尖變黑，并且黑色區(qū)域隨處理時間的延長而逐步擴大(圖1b)。

(a)葉片；(b)根圖1 鹽脅迫對甜菜幼苗表型的影響Fig.1 Effects of salt stress on the phenotype of sugar beet seedlings

2.2 甜菜鹽脅迫應答基因的鑒定

2.2.1 甜菜鹽脅迫下的基因表達圖譜

實驗對28天苗齡的甜菜幼苗進行300 mM NaCl處理， 5個取樣時間點分別為0 h、12 h、24 h、48 h和72 h。每個時間點設計3個生物學重復，最終構建了30個轉錄組測序文庫。RNA-seq測序獲得了超過13.5億的raw reads，平均每個樣本4500萬個reads。經(jīng)過質量控制篩選后共獲得13.2億，入選率達到97.6%的有效數(shù)據(jù)。這些Valid Data在基因組上的平均定位率為92.1%約12.2億，其中66.8%屬于Unique Mapped reads。RNA-seq數(shù)據(jù)信息、有效數(shù)據(jù)、映射比例見表1。

表1 表達譜數(shù)據(jù)概覽Tab.1 Summary of Expression profile data

本實驗共鑒定出27,677條表達基因，將基因表達水平量化成片段，每千堿基轉錄本每百萬映射reads(FPKM)。使用小提琴圖展示樣本基因表達分布統(tǒng)計如圖2。從圖中可以明確看出有樣本的測序質量較好，且各生物學重復樣本間重復性較好?；赑earson相關系數(shù)法對所有樣本間的基因表達相關性進行評估見圖3，結果表明各生物學重復樣本間基因表達相關性極高(>0.99)；同一組織內，不同處理時間點的各樣本間基因表達相關性相對較高，且鹽處理組之間的相關性要明顯高于處理組與對照組(ck)之間的相關性。但葉片和根不同的組織間的基因表達相關性較低。

圖2 RNA-seq所有樣本的基因表達分布統(tǒng)計Fig.2 Distribution of gene expression among RNA-seq samples

圖3 樣本間相關性分析Fig.3 Correlation analysis among RNA-seq samples

2.2.2 甜菜鹽脅迫應答相關基因的鑒定

過濾掉所有樣本中低豐度表達基因(FPKM<1)后，對剩余基因在鹽脅迫下的表達模式進行評估。最終鑒定到11,992條差異表達基因(differentially expressed gene，DEG)，其中葉片和根中的DEG分別有7,391和8,729條。對不同鹽處理時間下DEG的統(tǒng)計分析表明，葉片和根中分別含有39.4%(2,909條)和37.1%(3,238條)，僅在單一時間點呈現(xiàn)顯著差異的DEG，表明某些基因對鹽脅迫的應答存在時空特異性。葉片和根中DEG數(shù)量最高的處理時間點是12 h，其次是72 h、24 h和48 h(見圖4)。DEG的染色體定位顯示DEG的分布無明顯的染色體偏好(見圖5)。

圖4 葉片(a)和根(b)各樣本間差異基因分布圖Fig.4 Distribution of DEGs among samplesin leaves (a) and roots (b)

2.3 甜菜鹽脅迫應答基因的功能分析

對DEG的GO分析獲得了5,847條GO term注釋，其中3,074條為生物過程、719條為細胞成分和2,054條為分子功能。細胞組排名前5的分依次是細胞、細胞區(qū)域、膜、細胞器和膜區(qū)域；生物進程排名前5的依次是細胞進程、新陳代謝進程、刺激應答、生物調控和生物進程調控；分子功能排名前5的依次是結合、催化活性、轉運活性、分子結構活性和轉錄調控活性(圖6)。

本研究利用RNA-seq技術，對鹽脅迫后甜菜葉片和根進行了基因表達模式研究，在4個鹽處理時間點鑒定出11,992條DEG。通過對這些DEG功能分類表明，其功能可分為三個大類：信號轉導相關基因、轉錄調控基因、功能基因。(Ⅰ)信號轉導相關基因：該類基因調控信號轉導途徑來傳遞脅迫信號并激活下游應答途徑，如蛋白激酶、激素信號相關和ROS信號相關基因等；(Ⅱ)轉錄調控基因：該類基因通過對轉錄水平、轉錄后水平及翻譯水平上的調控重構植物在鹽脅迫下的基因表達系統(tǒng)，如可變剪接因子、轉錄因子和一些核酸內切酶等；(Ⅲ)功能基因：它們編碼脅迫應答相關物質的合成，直接參與植物脅迫過程中生存和脅迫后恢復的一些特定過程，如光合酶、滲透調節(jié)物合成相關基因、和自由基清除相關基因等。

從上到下的順序為：紅色：生物進程；藍色：細胞組分；綠色：分子功能圖6 差異基因的GO注釋Fig.6 GO annotation of DEGs

2.3.1 甜菜鹽脅迫應答過程中的信號轉導

在植物接觸NaCl的初始反應中，刺激在數(shù)秒內轉化為胞質Ca2+濃度變化，構成一種特殊的鈣信號[12]；隨后Ca2+信號依次被鈣依賴蛋白激酶(CDPK)、鈣調蛋白(calmodulin，CaM)/ calmodulin-like 蛋白、鈣調磷酸酶B類蛋白(CBLs)和膜聯(lián)蛋白傳遞解碼，進而觸發(fā)下游的細胞反應[13]。

通過本研究鑒定出了300多條差異蛋白激酶編碼基因，囊括了前述所有類型蛋白激酶：如編碼CDPK 26的Bv9_215750_yqnm、編碼calmodulin-like protein 2的Bv2_035110_kyjk、編碼CBL相作絲氨酸/蘇氨酸蛋白激酶的Bv7_157470_urmk、以及編碼膜聯(lián)蛋白RJ4的Bv2_041310_rnso等。Na+/H+逆轉運蛋白(NHX)是植物鹽超敏感信號通路(SOS)中的關鍵組分，它是植物細胞膜上發(fā)現(xiàn)的唯一能夠把胞質中Na+排除到細胞外的蛋白。本研究中發(fā)現(xiàn)了鹽脅迫下顯著上調的兩個NHX基因：Bv1_006770_xkmw和Bv_07430_jkji，它們可能在鹽脅迫信號轉導過程中起到重要作用。

胞質Ca2+濃度的提高還會刺激ROS和ABA的產生，在細胞水平上可以激活鹽過度敏感(SOS)通路，并排除多余的Na+來進而調控離子穩(wěn)態(tài)[14]。其中ROS信號的傳播依賴于呼吸爆發(fā)氧化酶同源物的激活，該同源物由乙烯、ABA依賴蛋白激酶和磷脂酸(PA)聯(lián)合控制。本研究中在鹽脅迫過程中編碼呼吸爆發(fā)氧化酶同源蛋白A(Bv5_099740_zpio)顯著上調，可能在葉片ROS信號轉導過程中發(fā)揮至關重要的作用。細胞內ABA的快速累積依賴于雙加氧酶對類胡蘿卜素前體的裂解，在鹽脅迫下甜菜中的Bv_01350_texx和NCED5(Bv_004050_xzzo)兩個9-順式-環(huán)氧類胡蘿卜素雙加氧酶基因NCED1顯著上調。我們認為鹽脅迫下誘導甜菜ABA合成。Ca2+或ROS均可調節(jié)ABA的釋放，進而引導下游離子通道的打開或閉合、滲透保護劑的合成以及轉錄因子的表達等等。此外，一些ABA和乙烯代謝及信號轉導相關基因如編碼PP2C家族成員的PP2C8(Bv3_052290_rcph)、PP2C24(Bv7_173230_chqm)和PP2C56(Bv8_197440_jnqo)；編碼玉米黃質環(huán)氧酶的Bv_003210_fkxn、以及編碼乙烯生物合成途徑關鍵酶1-氨基環(huán)丙烷-1-羧酸氧化酶基因(ACO)的Bv8_185300_wndh和Bv1_004040_gcto等均在鹽脅迫下顯著上調。

2.3.2 甜菜鹽脅迫應答過程中轉錄調控

轉錄因子可以與靶基因5′ORF區(qū)域目標序列特異性性結合，從而調控靶標基因在特定的時間和空間以特定的強度表達。它通過調控一個復雜的網(wǎng)絡而改變細胞內的轉錄組，從而在植物適應環(huán)境過程中發(fā)揮核心作用。本研究鑒定的500多條轉錄因子(圖7)包含了來自DEG中的53個家族。排名Top5的家族依次是MYB、AP2-EREBP、bHLH、MYB-related和NAC。其中的39個家已經(jīng)族被植物轉錄因子數(shù)據(jù)庫Plant Transcription Factor Database 4.0(http://planttfdb.cbi.pku.edu.cn/index.php?sp=Bvu)中的甜菜子數(shù)據(jù)庫收錄。此外，本研究鑒定到14個甜菜中存在的轉錄因子家族還未被該數(shù)據(jù)庫收錄如mTERF、ABI3/VP1、LOB、OFP等。

圖7 差異基因中的轉錄因子分布Fig.7 Transcription factor distribution in DEGs

DEG對鹽脅迫的應答通?？梢苑譃閮煞N類型：ABA依賴/非ABA依賴型。非ABA依賴基因啟動子中通常包含脫水反應元件/C-重復序列(DRE/CRTs)，并受AP2-EREBP轉錄因子的調控。本研究共鑒定出73條AP2-EREBP家族轉錄因子，其中絕大部分在鹽脅迫下差異表達，對該家族成員的染色體定位分析發(fā)現(xiàn)它們廣泛分布于甜菜的9條染色體上(圖8)，可能與甜菜非ABA依賴型鹽脅迫應答基因的激活密切相關。例如，AP2-EREBP家族的兩個脫水反應元件結合蛋白DREB 1(Bv3_066590_ignp)和DREB 2(Bv2_044600_drzk)在鹽脅迫下顯著上調。另一類依賴于ABA介導的鹽脅迫應答基因，通常在其啟動子區(qū)域包含一個保守的ABA應答順式作用元件(ABRE)。ABRE可以被含有基礎亮氨酸拉鏈結構(basic leucine zipper structure, bZIP)的轉錄因子識別。圖7是本研究中發(fā)現(xiàn)的編碼bZIP家族成員的DEG。如：Bv3_050990_kqjn、Bv9_214630_fnpa、Bv1_013750_smoy等。

Dicer-like(DCL)蛋白是雙鏈RNA(dsRNA)特異性的核糖核酸內切酶，在雙鏈RNA或帶莖環(huán)結構的單鏈RNA(如miRNA前體)切割過程中發(fā)揮關鍵作用[15]。本研究發(fā)現(xiàn)了7個鹽脅迫下差異表達的甜菜DCL基因，如Dicer同源蛋白1(Bv1_001950_imjw)。Argonaute蛋白直接與小RNA結合，并在所有小RNA介導的基因沉默過程中發(fā)揮重要作用，此外有研究表明Argonaute蛋白還與染色質修飾和可變剪接相關[16]。本研究發(fā)現(xiàn)了6個鹽脅迫下差異表達的甜菜Argonaute基因。上述結果表明，在多種轉錄調控途徑之間存在相互作用，構成了鹽脅迫下復雜的轉錄調控網(wǎng)絡。

圖8 甜菜差異表達AP2/ERF家族成員染色體分布Fig.8 Chromosome distribution of DE AP2/ERF family members in sugar beet

2.3.3 甜菜鹽脅迫應答過程中的功能基因

植物適應鹽脅迫的第一步是解決環(huán)境中高鹽度引起的滲透脅迫，甜菜堿作為細胞內的滲透保護劑在這一過程中起著至關重要的作用。編碼甜菜堿合成的兩個關鍵酶基因膽堿單加氧酶(Bv6_146100_wdro)和甜菜堿醛脫氫酶(Bv5_116230_ntjn)在鹽脅迫下顯著上調。水通道蛋白也是一種響應鹽脅迫的功能蛋白。編碼甜菜水通道蛋白PIP2-1的基因Bv7_162410_cqnr在鹽脅迫下表現(xiàn)為葉片和根中均顯著上調。超氧化物歧化酶(SOD)和過氧化物酶(POD)經(jīng)常用來評估植株的抗逆性，編碼甜菜Gu/Zn-SOD基因Bv_08010_smke、編碼Fe-SOD基因Bv5_097450_sxsu和編碼Mn-SOD 基因Bv2_045060_jxsu在鹽脅迫下表達上調。根中編碼POD 66的Hub基因Bv3_066080_qmgn也在鹽脅迫下顯著上調。硫氧還蛋白(thioredoxin)和谷氧還蛋白(glutaredoxin)在維持植物細胞的氧化還原狀態(tài)中也起著重要作用[17]。本研究鑒定了16條編碼硫氧還蛋白的DEG和10條編碼谷氧還蛋白的DEG。包括葉中編碼F型硫氧還蛋白的Hub基因Bv4_075150_andd。另外，本研究發(fā)現(xiàn)了25條鹽脅迫下差異表達的甜菜GST基因，其中甜菜葉片中的Hub基因Bv4_095830_aore。ASA-GSH循環(huán)是植物清除ROS的一個重要的途徑。谷胱甘肽過氧化物酶(GPX)是該系統(tǒng)的關鍵酶。本研究中編碼GPX的基因Bv5u_121470_tcou的表達受鹽脅迫的顯著誘導。

葉綠體是對鹽脅迫最敏感的細胞器。參與光系統(tǒng)Ⅰ光合電子傳遞的質子梯度調節(jié)蛋白5編碼基因Bv3_048340_odeme、參與光反應調控的STN7蛋白激酶編碼基因Bv6_155590_ugjw等直接參與光合作用的基因在鹽脅迫下顯著上調。

2.4 甜菜鹽脅迫應答core基因的篩選與功能分析

部分DEG在甜菜葉片或根的表達明顯不同具有時空特異性，而有部分DEG在葉片和根中的所有鹽處理時間點下均呈現(xiàn)顯著差異表達。即這些鹽脅迫應答基因的表達不具有明顯的時空特異性，這類DEG可能構成了甜菜對鹽脅的基礎應答，本研究將中將這類DEG歸為core基因。我們通過交叉比較了葉片和根中各處理時間點之間的DEG，最終鑒定出244條core基因，其中206條在葉片和根之間的表達趨勢相同，38條表達趨勢相反(圖9)。

圖9 基礎應答基因表達模式分析Fig.9 Expression pattern analysis of core DEGs

對core基因的GO分析顯示：超過90%的core基因都可以注釋到氧化還原進程(GO:0055114)；超過80%的core基因都可以注釋到轉錄調控，DNA-模板(GO:0006355)；此外許多core基因被注釋到與脅迫應答直接相關的生物進程如缺水響應(GO:0009414)、脫落酸應答(GO:0009737)、氧化應激響應(GO:0006979)、鹽脅迫響應(GO:0009651)和滲透脅迫響應(GO:0006970)等(圖10a)。分子功能注釋顯示core基因主要執(zhí)行蛋白結合(GO:0005515)和金屬離子結合(GO:0046872)等功能。KEGG分析顯示core基因主要富集在淀粉和蔗糖代謝與植物激素信號轉導通路中(圖10b)。

圖10 基礎應答基因的GO注釋(a)和KEGG分析(b)Fig.10 GO and KEGG analysis of core DEGs

對所有core DEG的功能進行分析后發(fā)現(xiàn)，鹽脅迫下葉片和根中的蔗糖合成酶(Bv7_163460_jmqz)和β-淀粉酶1(Bv_005770_atmm)顯著上調。研究顯示許多植物在非生物脅迫下會積累蔗糖合成酶和β-淀粉酶1，糖代謝產物不僅在滲透調節(jié)中發(fā)揮重要作用，而且在光合作用受限時也發(fā)揮著供應能量和碳源的作用。植物在感知環(huán)境脅迫后，通常會犧牲生長并激活保護性脅迫反應。ABA在脅迫信號整合和下游反應激活過程中發(fā)揮重要作用。3個ABA相關core DEG(Bv3_052290_rcph、Bv4_087620_atdx和Bv7_173230_chqm)在鹽脅迫下顯著上調。Bv3_052290_rcph和Bv7_173230_chqm編碼兩個蛋白磷酸酶PP2C(protein phosphatase 2C)家族成員，可以作為負調控因子在ABA信號轉導中發(fā)揮重要作用[18]；另外一項研究顯示鹽脅迫下擬南芥PP2C也可以與AtHKT1相互作用，調節(jié)Na+的分布[19]。此外，本研究還發(fā)現(xiàn)苯丙素和黃酮醇相關core DEG在葉片和根中表達均上調。非生物脅迫條件下可能激活苯丙素和黃酮醇生物合成途徑，通過累積各種酚類化合物清除有害ROS[20]?？傊@些core DEG通過調控甜菜的蔗糖-淀粉代謝、調節(jié)植物激素信號通路、激活苯丙素和黃酮醇代謝通路來響應鹽脅迫。我們認為對這些代謝途徑的調控構成了甜菜對鹽脅迫的基礎應答。

2.5 甜菜響應鹽脅迫的hub基因鑒定與生物信息學分析

關鍵基因(hub gene)指在某些生物學進程中起至關重要作用的基因，在一些相關通路中，其他基因的調控往往會受到該基因的影響[21]。因此，hub gene 通常可以作為重要的作用靶點。目前可以通過共表達或蛋白相互作用構建互作網(wǎng)絡，然后根據(jù)網(wǎng)絡拓撲結構篩選關鍵基因。

2.5.1 甜菜DEG共表達網(wǎng)絡構建

本研究依據(jù)WGCNA算法分別利用葉片和根中DEG的表達值(FPKM值)構建共表達模塊。最終分別將葉片和根中的DEG劃歸為11個共表達模塊。無法劃歸任何模塊的基因均統(tǒng)一劃歸為灰色模塊，本研究中葉片的11個模塊不包含灰色模塊，而根中的11個模塊包含灰色模塊，該模塊擁有28個DEG，占根中DEG的0.32%。繼續(xù)使用WGCNA算法構建模塊與處理時間的對應關系(圖11)。

(a)葉片；(b)根圖11 模塊-樣本相關性分析Fig.11 Module-sample correlation analysis

通常與樣本具有較高Spearman相關系數(shù)的模塊也是與處理時間最相關的模塊。依據(jù)P值最終在葉片和根中分別篩選了5個和4個符合標準的模塊作為關鍵模塊進行進一步分析。其中葉片中的5個模塊分別是：對應0 h的包含2181條DEG的藍綠色模塊 (Module Turquoise)、對應12 h的包含1115條DEG的棕色模塊(Module Brown)、對應24 h的包含58條DEG的黃綠色模塊(Module Greenyellow)、對應48 h的包含121條DEG的紫色模塊(Module Purple)以及對應72 h的包含1638條DEG的藍色模塊(Module Blue)(圖11a)。根中的4個模塊分別是：對應0 h的包含2501條DEG的藍綠色模塊(Module Turquoise)、對應12 h包含1038條DEG的黃色模塊(Module Yellow)、對應24 h的包含980條DEG的綠色模塊(Module Green) 以及對應于72 h的包含2052條DEG的藍色模塊(Module Blue)(圖11b)。根中48 h無對應的顯著相關(P<0.1)模塊。

2.5.2 關鍵模塊基因表達模式分析

表達模式分析顯示關鍵模塊間的基因表達存在顯著差異，各模塊內基因均在模塊關聯(lián)的時間點處呈現(xiàn)顯著差異表達(圖12)。例如葉片Brown模塊內基因在鹽處理12 h時呈現(xiàn)顯著差異表達(圖12a)，根 Yellow模塊內基因在鹽處理12 h時顯著差異表達(圖12b)。總的來看，對照組關聯(lián)模塊(Turquoise)的大部分基因表達被鹽脅迫抑制，處理組關聯(lián)模塊的大部分基因表達受鹽脅迫誘導。

(a)葉片中的5個關鍵模塊；(b)根中的4個關鍵模塊圖12 每個關鍵模塊的基因表達模式Fig.12 The expression profiles of genes in each critical module

2.5.3 關鍵模塊基因功能分析

葉片中關鍵模塊基因的GO富集分析顯示(圖13)：鹽脅迫下Turquoise模塊內富集度較高的生物進程是翻譯(GO:0006412)和核糖體大亞基組合(GO:0000027)；富集度較高的的細胞組分是葉綠體基質(GO:0009570)和胞質大核糖體亞基(GO:0022625)；富集度較高的分子功能是核糖體的結構組成部分(GO:0003735)。Brown模塊內基因富集度較高的生物進程是乙烯激活信號通路負調控(GO:0010105)；富集度較高的細胞組分是細胞核(GO:0005634)；富集度較高的分子功能是鈣離子結合(GO:0005509)。Greenyellow模塊內基因主富集度較高的生物進程是RNA聚合酶Ⅱ對轉錄的調控(GO:0006357)和磷酸化(GO:0016310)；富集度較高的細胞組分是核仁(GO:0005730)；富集度較高的分子功能是RNA聚合酶Ⅱ調控的序列特異性DNA結合(GO:0000977)。Purple模塊內基因富集度較高的細胞組分是中央葉泡(GO:0042807)；富集度較高的分子功能是蛋白二聚活動(GO:0046983)。Blue模塊內基因富集度較高的生物進程是蛋白酶體介導的泛素依賴的蛋白降解過程(GO:0043161)和蛋白水解相關細胞蛋白分解代謝過程(GO:0051603)；富集度較高的細胞組分是細胞外空間(GO:0005615)和溶酶體(GO:0005764)；富集度較高的分子功能是半胱氨酸型肽鏈內切酶活性(GO:0004197)。

圖13 葉片關鍵模塊基因的GO富集分析Fig.13 Go enrichment analysis of critical module genes in leaves

根中關鍵模塊基因的GO富集分析顯示(圖14)：鹽脅迫下Turquoise模塊內富集度較高的生物進程是草酸代謝過程(GO:0033609)、細胞表面受體信號通路(GO:0007166)和防御應答(GO:0006952)；富集度較高的細胞組分是膜的整體構件(GO:0016021)、細胞壁(GO:0005618)和液泡膜(GO:0005774)；富集度較高的分子功能是ATP結合(GO:0005524)、蛋白絲氨酸/蘇氨酸激酶活性(GO:0004674)和草酸脫羧酶活性(GO:0046564)。Yellow模塊內富集度較高的生物進程是鹽脅迫應答(GO:0009651)和類黃酮生物合成(GO:0009813)；富集度較高的細胞組分是胞質(GO:0005829)；富集度較高的分子功能是槲皮素3-O-葡糖基轉移酶活性(GO:0080043)和槲皮素7-O-葡糖基轉移酶活性(GO:0080044)。Green模塊內富集度較高的生物進程是植物型次生細胞壁發(fā)生(GO:0009834)、木質素生物合成(GO:0009809)和脂肪酸生物合成(GO:0006633)；富集度較高的細胞組分是膜的整體構件(GO:0016021)和細胞外區(qū)域(GO:0005576)；富集度較高的分子功能是過氧化物酶活性(GO:0004601)、血紅素結合(GO:0020037)和作用于酯鍵的水解酶活性(GO:0016788)。Blue模塊內富集度較高的生物進程是DNA模板轉錄負調控(GO:0045892)，DNA模板轉錄正調控(GO:0045893)和脫落酸應答(GO:0009737)；富集度較高的細胞組分是細胞核(GO:0005634)；富集度較高的分子功能是核酸結合(GO:0003676)。

圖14 根關鍵模塊基因的GO富集分析Fig.14 Go enrichment analysis of critical module genes in roots

葉片中關鍵模塊的基因KEGG富集分析顯示(圖15a)：Turquoise模塊基因富集度較高的通路是核糖體(ko03010)和脂肪酸生物合成(ko00061)；Brown模塊基因富集度較高的通路是精氨酸和脯氨酸的代謝(ko00330)和內吞作用(ko04144)；Greenyellow模塊基因富集度較高的通路是光合作用生物中的碳固定(ko00710)和植物-病原互作(ko04626)；Blue模塊基因富集度較高的通路是植物激素信號轉導(ko04075)和其他類型的鄰聚糖生物合成(ko00514)。根中關鍵模塊的基因KEGG富集分析顯示(圖15b)：Turquoise模塊基因富集度較高的通路是植物-病原互作(ko04626)和RNA聚合酶(ko03020)；Yellow模塊基因富集度較高的通路是谷胱甘肽代謝(ko00480)和花生四烯酸代謝(ko00590)；Green模塊基因富集度較高的通路是苯丙素的生物合成(ko00940)和丙酮酸代謝(ko00620)；Blue模塊基因富集度較高的通路是RNA運輸(ko03013)和mRNA監(jiān)測通路(ko03015)。

總的來看，使用WGCNA法對DEG進行歸類后，各關鍵模塊基因的表達模式和功能均存在高度的特異性，符合我們的預期目標。有助于從大量DEG之間挖掘關鍵耐鹽基因，探尋甜菜應對鹽脅迫的分子機理。

2.5.4 鹽脅迫應答Hub基因的鑒定

本研究首先利用加權基因共表達網(wǎng)絡計算模塊內基因的連通度，通常連通度值越高的基因越可能是潛在的Hub基因。然而當連通度值極度相近的情況下(如葉片Turquoise模塊中Bv5_105240_scsk和Bv1_016890_iiqe的連通度值分別為1140.9和1140.7)，以連通度值作為唯一判定標準是不夠科學的。因此我們選擇模塊內連通度值前20 %的基因作為候選Hub基因，利用STRING構建這些基因編碼蛋白的互作網(wǎng)絡；然后利用CytoHubba中的11種算法分別評估各模塊內的top10 Hub基因；最后以候選基因在11種算法top10 Hub基因中的出現(xiàn)頻次聯(lián)合基因的連通度值為每個模塊篩選出6個最可能的Hub基因。

圖15 葉片(a)和根(b)中關鍵模塊基因的KEGG富集分析Fig.15 KEGG enrichment analysis of critical module genes inleaves (a) and roots (b)

2.5.4.1 Hub基因在葉片中的鑒定與功能分析

葉片中篩選出的Hub基因見表3-2。鹽脅迫下，Turquoise模塊內大部分基因的表達呈現(xiàn)下調趨勢。Hub基因Bv2_041930_qfdi編碼的膽色素原脫氨酶是植物葉綠素生物合成過程中的關鍵酶，甜菜葉片遭受鹽脅迫后葉綠素含量顯著降低可能與受到該Hub基因網(wǎng)絡的調控。Hub基因Bv1_000340_ezio編碼葉綠體30S 核糖體蛋白、Hub基因Bv2_036440_xjzx編碼40S核糖體蛋白、Hub基因Bv5_110020_fucg編碼mRNA結合伴侶ACD11，這些基因參與mRNA的結合及后續(xù)加工過程，它們的下調與鹽脅迫后生長發(fā)育受阻相符，這些Hub基網(wǎng)絡可能參與調控植株從生長向應激狀態(tài)的轉換。Hub基因Bv2_029170_pxtq編碼環(huán)阿屯醇-C-24-甲基轉移酶(SMT)，它是植物甾醇代謝途徑中的關鍵限速酶，從圖16可以看出，抑制該基因的表達可以提高植物膽甾醇的累積，阻礙了環(huán)阿屯醇前體流向豆甾醇合成分支，促進植物甾體類化合物流向其他合成分支方向(如油菜素甾醇)，對擬南芥SMT基因的早期研究支持了這種推斷[22]。

Brown模塊Hub基因Bv2_023810_guyf編碼DExH-box ATP依賴RNA解旋酶DExH 17，最新研究顯示該基因調控核酸-蛋白質相分離為細胞提供了各種RNA處理步驟的時空控制，通過引導RNA分子進入相分離的區(qū)室，調節(jié)RNA成熟狀態(tài)、核糖核蛋白復合體(RNP)的組成并最終調節(jié)RNA命運[23]。

表2 葉片關鍵模塊中鑒定的Hub基因Tab.2 The Hub genes detected in critical modules in leaves

本研究中該基因在鹽脅迫后顯著上調并在12 h時達到峰值，說明以它為核心的基因網(wǎng)絡可能通過參與早期鹽脅迫應答基因的加工處理為甜菜鹽脅迫應激過程做出貢獻。Hub基因Bv6_130450_umcu編碼葉綠體早期脫水響應蛋白。Hub基因Bv5_099740_zpio編碼呼吸爆發(fā)氧化酶同源蛋白A參與ROS信號的轉導。Hub基因Bv6_149930_qmfq編碼莽草酸激酶，莽草酸是芳香氨基酸、生物堿和吲哚衍生物等的前體，該酶通過催化莽草酸的磷酸化參與調控莽草酸代謝，已有研究顯示莽草酸途徑與水稻的鹽脅迫應答過程相關[24]。

綠色：甜菜中存在組件；黃色：DEG；紅色邊框：植物甾醇代謝圖16 鹽脅迫下甜菜類固醇生物合成通路分析Fig.16 Pathway analysis of steroid biosynthesis in sugar beet under salt stress

Greenyellow模塊受樣本總量限制只篩選3個Hub基因。Bv4_074970_andd編碼葉綠體中F型硫氧還蛋白，Bv4_095830_aore編碼GST蛋白，以他們?yōu)楹诵牡腍ub基因網(wǎng)絡可能為細胞的抗氧化過程做出巨大貢獻。

葉片Purple模塊Hub基因Bv1_016400_jjoj編碼葉綠體中的伴侶蛋白dnaJ 11；Hub基因Bv4_086990_pgms編碼泛素連接酶E2 26；以它們?yōu)楹诵牡幕蚓W(wǎng)絡通過調控蛋白質的加工參與鹽脅迫應答。兩個Hub基因編碼2種轉錄因子bHLH 35和鋅指蛋白22。Hub基因Bv_005490_qmph編碼磷酸代謝蛋白8；Hub基因Bv7_172340_kzsr編碼非特異性磷脂酶C1(NPC)，該酶可以催化磷脂酰膽堿(PC)水解產生甘油二酯(DAG)，DAG再次代謝后生成磷脂酸(PA)，PA是細胞內重要的第二信使，可以激活多種脅迫應答途徑。以它們?yōu)楹诵牡幕蚓W(wǎng)絡可能通過PLC-DGK/PA途徑[25]參與鹽脅迫信號的轉導與應答。

葉片Blue模塊Hub基因Bv5_118050_deeh編碼一種E3泛素蛋白U-box結構域的蛋白52。Hub基因Bv6_139660_wjii編碼亮氨酸拉鏈蛋白ATHB-40；hub基因Bv6_133290_wkjd編碼翻譯起始因子IF-2；以他們?yōu)楹诵牡幕蚓W(wǎng)絡可能通過調控轉錄和翻譯過程響應鹽脅迫應答。Hub基因Bv2_042060_ccft編碼膽堿激酶2，以該基因為核心的網(wǎng)絡最終也可能通過PLC-DGK/PA途徑參與鹽脅迫信號的轉導與應答。脅迫72 h后的甜菜老葉呈出萎蔫與失綠的性狀，說明長時間的鹽脅迫已經(jīng)對甜菜幼苗的存活造成一定影響，Bv2_033120_umqq編碼液泡加工酶(VPE)；以VPE為核心的基因網(wǎng)絡可能與細胞衰老相關。

2.5.4.2 Hub基因在根中的鑒定與功能分析

根中Hub基因見表3。磺化反應是機體對內源性和外源性化合物進行生物轉化的關鍵反應。Turquoise模塊Hub基因Bv2_042180_yfwq編碼胞質磺基轉移酶，它可以催化生物體內黃酮類化合物轉化成硫酸鹽而降低利用率[26]，因此以該基因為核心的網(wǎng)絡可能通過調控激素和黃酮類化合物的利用來影響甜菜的發(fā)育與抗逆。此外，一些次生代謝相關的Hub基因如編碼重金屬相關異戊烯基類黃酮植物蛋白的Bv_016430_xjdf和編碼7-去氧番木鱉酸葡萄糖轉移酶的Bv4_077620_xkxs也在受到脅迫后顯著下調。這些結果說明甜菜根細胞在遭受鹽脅迫后會調整次生代謝方向，調控植株從生長狀態(tài)向應激狀態(tài)轉換以適應環(huán)境改變。另外編碼包含L型凝集素結構的受體激酶、LRR受體類絲氨酸/蘇氨酸蛋白激酶和包含重復的錨蛋白NPR4的3個Hub基因都與信號轉導相關，這些Hub基因網(wǎng)絡可能通過調控甜菜根系在鹽脅迫下的信號轉導來參與鹽脅迫應答過程。

表3 根中從WGCNA核心模塊中鑒定的Hub基因Tab.3 The Hub genes detected in critical WGCNA modules in roots

Yellow模塊中，Hub基因Bv5_124870_ydpo編碼解毒蛋白 48，以該基因為核心的網(wǎng)絡可能與金屬離子解毒相關。此外，以Bv6_129310_ndoa為核心的Hub基因網(wǎng)絡可能通過調控依賴于ATP的離子運輸輔助解毒過程。

以Hub基因Bv1_009510_hfps為核心的網(wǎng)絡與蛋白泛素化過程相關。編碼液泡氨基酸轉運蛋白的Hub基因Bv2_046760_tquf在鹽脅迫12時表達上調了80倍，液泡占據(jù)植物細胞的大部分體積，它參與調節(jié)細胞代謝，控制細胞內物質累積和運輸，提高植物適應環(huán)境變化和生存的能力，該Hub基因網(wǎng)絡可能通過調控氨基酸的調配過程在鹽脅迫過程中發(fā)揮重要作用。Hub基因Bv3_067190_wskz編碼過氧化物酶體中脂肪酸β氧化多功能蛋白MFP 2，以該基因為核心的網(wǎng)絡可能通過調控能量代謝為甜菜根細胞鹽脅迫應答供能。編碼谷氨酸脫羧酶4的Hub基因Bv6_131780_uxzh在鹽脅迫12時表達上調了7倍，該酶催化谷氨酸鹽生成γ-氨基丁酸(GABA)，許多研究顯示植物在逆境脅迫下會大量累積GABA，首先GABA的合成過程會消耗大量氫離子從而維護細胞內的pH環(huán)境；其次生成的GABA可以作為小分子滲透調節(jié)物質；過量的GABA會抑制植物的生長，GABA含量的升高可以提高植物對氧化脅迫的耐受性；此外，迅速累積的GABA可能是一種內源性逆境信號分子[27-28]。因此該Hub基因網(wǎng)絡通過調控GABA相關代謝參與甜菜根細胞對鹽脅迫的應答。

Green模塊Hub基因Bv2_044780_rzxx編碼3-酮脂酰-CoA合酶，與內質網(wǎng)中的脂肪酸鏈的伸長相關，該Hub基因網(wǎng)絡可能影響膜脂的不飽和度以及非雙分子層脂質的比例，從而改變細胞膜的流動性來調控甜菜根系對鹽脅迫的應答。Hub基因Bv1_016480_prnw編碼IN2-1同源蛋白，該蛋白屬于GST蛋白，以該Hub基因為核心的網(wǎng)絡可能在根細胞非酶促抗氧化防御過程中發(fā)揮重要作用。Hub基因Bv2_043550_pccg編碼ABC轉運蛋白G超家族成員32。木質素是對植物結構支撐和水分運輸十分重要的生物大分子，它也可以提高作物的抗逆性，編碼過氧化物酶66的Hub基因Bv3_066080_qmgn可以催化木質素的生物合成。有研究顯示鹽脅迫會抑制根系的生長并誘導根系木質化[28-30]。以該基因為核心的網(wǎng)絡可能通過調控根系木質素的合成響應鹽脅迫。

總的來看，我們篩選出的Hub基因在生物學上具有重要意義，有助于解析甜菜鹽脅迫應答的分子機理。

2.5.5 甜菜轉錄水平的鹽脅迫應答網(wǎng)絡構建

分別為葉片和根構建了由core基因、Hub基因和重要鹽脅迫應答基因組成的轉錄水平鹽脅迫應答網(wǎng)絡。同時利用STRING數(shù)據(jù)庫對這些基因間編碼蛋白間的互作關系進行注釋。在葉片中(圖17)，鹽脅迫造成甜菜葉片生長遲滯，并導致葉綠素合成受阻，因此葉片細胞下調了大量核糖體編碼基因的表達，以削減mRNA后續(xù)加工過程；同時迅速調整次生代謝物的合成流向于有利于脅迫應答的分支。在受到鹽脅迫的早期(12 h和24 h)，ATP依賴RNA解旋酶基因表達顯著上調，調控脅迫應答相關mRNA的緊急加工；同時一些ROS信號和ROS清除相關基因的表達顯著上調。在鹽脅迫的中期(48 h)，一些轉錄因子如bHLH35和鋅指蛋白22編碼基因顯著上調，同時NPC等基因的上調可能激活PLC-DGK/PA途徑調控鹽脅迫應答過程。在鹽脅迫的晚期(72 h)，VPE等基因的上調可能與葉片細胞的衰老相關。值得注意的是，葉片中許多Hub基因均定位于葉綠體，反映出該細胞器在鹽脅迫應答過程中的重要作用。此外，鹽脅迫后多種泛素蛋白編碼基因在各時間點特異性高表達，說明蛋白修飾過程在鹽脅迫應答過程中發(fā)揮重要作用。

在根中(圖18)，遭受鹽脅迫后根系生長受到抑制，許多受體激酶基因表達顯著下調，次生代謝方向迅速發(fā)生改變。在受到鹽脅迫的早期(12 h和24 h)，GABA等滲透調節(jié)物質迅速累積，多種轉運蛋白基因表達顯著上調促進根細胞內物質的重新調配，次生代謝流向木質素生物合成代謝，誘導根系木質化；脂肪酸分解代謝加劇為細胞提供能量；此外，谷胱甘肽相關基因也在這一時期顯著上調，有助于根細胞的抗氧化防御。在鹽脅迫的晚期(72 h)，一些轉錄因子如RAP2-4和ERF54、轉錄延伸因子AFF4和鈣調蛋白編碼基因顯著上調。這些基因通過調控轉錄過程參與鹽脅迫應答。

正方形代表與時間點對應的葉中關鍵模塊；大圓點代表hub基因，core基因聚集在中心區(qū)域，紅色：表達上調，綠色：表達下調；紅線代表STRING數(shù)據(jù)庫注釋的互作關系圖17 葉片轉錄水平鹽脅迫應答網(wǎng)絡Fig.17 Transcriptional network of leaves in responses to salt stress

菱形代表與時間點對應的葉中關鍵模塊；大圓點代表hub基因，core基因聚集在中心區(qū)域，紅色：表達上調，綠色：表達下調；紅線代表STRING數(shù)據(jù)庫注釋的互作關系圖18 根轉錄水平鹽脅迫應答網(wǎng)絡Fig.18 Transcriptional network of roots in responses to salt stress

2.6 qPCR檢驗Hub基因在鹽脅迫下的表達模式

采用實時熒光定量PCR(qRT-PCR)法在葉片和根中分別隨機抽選15條Hub基因進行表達模式驗證。結果顯示，葉片中有11條Hub基因在兩種檢測方法間的R2值(回歸方程決定系數(shù))大于0.9；有1條Hub基因在兩種檢測方法間的R2值在0.8～0.9之間；有2條Hub基因在兩種檢測方法間的R2值在0.7～0.8之間，僅有1條Hub基因在兩種檢測方法間的R2值小于0.7(圖19a)。根中Hub基因的qRT-PCR定量結果與葉片相似，有9條Hub基因在兩種檢測方法間的R2值大于0.9；又3條Hub基因在兩種檢測方法間的R2值在0.8～0.9之間；有1條Hub基因在兩種檢測方法間的R2值在0.7～0.8之間，僅有2條Hub基因在兩種檢測方法間的R2值小于0.7(圖19b)。

(a)葉片；(b)根；x軸為qRT-PCR結果，y軸為RNA-SEQ結果圖19 Hub基因qRT-PCR和高通量測序數(shù)據(jù)比對Fig.19 Comparison betweenqRT-PCR and sequencing of Hub genes

總的來看，高通量測序檢測的Hub基因的表達與和qPCR驗證結果在趨勢上高度一致，說明我們的測序結果高度可信。個別Hub基因在兩種檢測方法間存在一定差異(R2值較低)可能是由于檢測原理間存在差異或目的基因具有特殊結構(如莖環(huán))等因素造成的。

3 結論

(1)RNA-seq在甜菜中共鑒定出了11,992條鹽脅迫應答基因，其中葉片和根中分別包含7391條和8729條。樣本間DEG比對分析，鑒定出244條core應答基因，我們認為這些基因通過調控甜菜的蔗糖-淀粉代謝、調節(jié)植物激素信號通路、激活苯丙素和黃酮醇代謝通路來響應鹽脅迫進而構成了甜菜對鹽脅迫的基礎應答。

(2)通過對Hub基因的表達模式進行分驗證表明及功能分析顯示這些Hub基因具有重要的生物學意義，有助于解析甜菜對鹽脅迫應答的分子機理。

(3)本研究利用core基因、Hub基因和重要鹽脅迫應答相關DEG構建了甜菜轉錄水平鹽脅迫應答網(wǎng)絡，通過對該網(wǎng)絡中DEG的功能分析，從轉錄層面論述了甜菜對鹽脅迫應答的分子機理?？傊狙芯渴状螢樘鸩藰嫿艘粋€較為全面的轉錄水平鹽脅迫應答網(wǎng)絡。