段本耀　張學(xué)杰　樊守金　張洛艷

(山東師范大學(xué)，濟(jì)南　250014)

低溫脅迫是目前威脅全球作物產(chǎn)量的重要非生物脅迫因子之一，嚴(yán)重影響植物的正常生長(zhǎng)發(fā)育，降低農(nóng)作物的產(chǎn)量，是農(nóng)業(yè)生產(chǎn)中常見的限制因子。研究表明，每年低溫使全球農(nóng)作物的損失高達(dá)數(shù)千億美元[1～2]。在全球環(huán)境氣候變化的背景下，培育具有較高產(chǎn)量的耐受低溫脅迫的作物是解決目前日益嚴(yán)重的糧食危機(jī)的有效手段之一。模式植物如擬南芥低溫脅迫應(yīng)答的生理、分子和遺傳響應(yīng)機(jī)制對(duì)研發(fā)低溫耐受作物有重要意義，同時(shí)還有助于挖掘植物響應(yīng)非生物脅迫的進(jìn)化分子基礎(chǔ)。

低溫脅迫對(duì)于植物的生理水平傷害主要體現(xiàn)在水分運(yùn)輸和存貯、礦質(zhì)營(yíng)養(yǎng)、光合作用、呼吸作用和總代謝水平的改變[3～5]。植物進(jìn)化出一套復(fù)雜體系來響應(yīng)低溫脅迫，包括信號(hào)轉(zhuǎn)導(dǎo)、基因調(diào)控、蛋白修飾和代謝調(diào)節(jié)等生理和生化等途徑[3,6～7]?；蚣捌潢P(guān)聯(lián)的生物學(xué)過程共同構(gòu)成了上述生理和生化響應(yīng)機(jī)制的分子基礎(chǔ)。信號(hào)通路和調(diào)控途徑的參與基因，包括植物內(nèi)源激素途徑關(guān)聯(lián)基因、信號(hào)轉(zhuǎn)導(dǎo)基因、轉(zhuǎn)錄因子(transcriptional factors，TFs)和蛋白激酶(protein kinases，PKs)等調(diào)節(jié)基因由于通常位于脅迫響應(yīng)上游，會(huì)首先感受及傳導(dǎo)脅迫刺激信號(hào)，并調(diào)控下游單功能基因的表達(dá)，因此被認(rèn)為發(fā)揮更重要的作用。雖然目前的研究已經(jīng)闡明了一部分植物應(yīng)答低溫脅迫的分子機(jī)制，然而由于生物學(xué)性狀分子基礎(chǔ)的復(fù)雜性，仍有大量潛在冷脅迫應(yīng)答關(guān)聯(lián)的單功能基因和調(diào)節(jié)基因有待挖掘。

微陣列芯片測(cè)序是在組學(xué)水平挖掘生物學(xué)途徑的有效手段之一，目前已有數(shù)百個(gè)非生物脅迫研究采用微陣列芯片測(cè)序來挖掘脅迫應(yīng)答潛在的分子基礎(chǔ)[8]，目前，十多個(gè)擬南芥的低溫脅迫應(yīng)答測(cè)序分析已經(jīng)鑒定出了13 000余個(gè)潛在的冷脅迫應(yīng)答基因。然而，不同實(shí)驗(yàn)鑒定的差異表達(dá)基因(differently expressed genes,DEGs)卻有一定的差異性(heterogeneity)[9～10]。薈萃整合分析是指利用不同的統(tǒng)計(jì)學(xué)及信息學(xué)手段，通過降低不同實(shí)驗(yàn)數(shù)據(jù)間的不一致性挖掘多個(gè)實(shí)驗(yàn)間的共性的整合分析方法。目前主要報(bào)道了4種常用的微陣列芯片表達(dá)數(shù)據(jù)的整合方法：不同實(shí)驗(yàn)的表達(dá)數(shù)據(jù)歸一化后直接整合(directly merge after normalization)，P值整合(Pvalues)、有效樣本數(shù)目整合(combine effect sizes，ES)和排位整合(combine ranks)等。歸一化后直接整合的方法通常受限于要求同一批次的實(shí)驗(yàn)或者要求一致的測(cè)序平臺(tái)[11]；而P值整合和排位整合則要求測(cè)序?qū)嶒?yàn)提供包含P值的差異表達(dá)基因名單或者基因表達(dá)名單[11～13]。

有效樣本數(shù)目指的是兩組樣本間在通過數(shù)據(jù)集(dataset)標(biāo)準(zhǔn)差(s.d)歸一化后均值的差異[12]。NetworkAnalyst網(wǎng)絡(luò)工具是針對(duì)微陣列和二代測(cè)序所設(shè)計(jì)的薈萃分析綜合工具，其中批量效應(yīng)矯正分析、數(shù)據(jù)不一致度的檢測(cè)及有效樣本數(shù)目整合分析可以科學(xué)地評(píng)價(jià)和有效地降低不同表達(dá)數(shù)據(jù)間的不一致性，且僅需測(cè)序?qū)嶒?yàn)提供原始表達(dá)數(shù)據(jù)矩陣(series matrix)，因此被廣泛的應(yīng)用到人類疾病關(guān)聯(lián)的基因及microRNA的微陣列測(cè)序及二代測(cè)序數(shù)據(jù)分析中[14]。然而基于有效樣本數(shù)目整合的薈萃分析還未被應(yīng)用于植物應(yīng)答非生物脅迫的測(cè)序整合分析領(lǐng)域。

在本研究中，我們通過薈萃整合分析擬南芥在低溫脅迫條件和標(biāo)準(zhǔn)條件下的微陣列表達(dá)譜芯片數(shù)據(jù)挖掘出了擬南芥中潛在的低溫脅迫響應(yīng)基因。由于所選擇表達(dá)譜數(shù)據(jù)間的差異性，我們選擇NetworkAnalyst網(wǎng)站的批量效應(yīng)矯正分析和有效樣本數(shù)目整合方法來降低實(shí)驗(yàn)間的不一致性并整合表達(dá)數(shù)據(jù)。本研究重點(diǎn)關(guān)注挖掘擬南芥不同測(cè)序?qū)嶒?yàn)間共有的響應(yīng)冷脅迫的分子基礎(chǔ)，并基于此來探索植物響應(yīng)冷脅迫的分子機(jī)制，并為抗冷作物的研發(fā)提供一定的理論基礎(chǔ)。

1　材料與方法

1.1　文獻(xiàn)數(shù)據(jù)庫挖掘

本研究采用文獻(xiàn)搜索策略在不同數(shù)據(jù)庫中搜索擬南芥的低溫脅迫響應(yīng)表達(dá)芯片數(shù)據(jù)。首先，采用下列關(guān)鍵詞在Gene Expression Omnibus數(shù)據(jù)庫(GEO,www.ncbi.nlm.nih.gov/geo/)、PubMed數(shù)據(jù)庫(https://www.ncbi.nlm.nih.gov/pubmed)和ArrayExpress數(shù)據(jù)庫(www.ebi.ac.uk/arrayexpress)中進(jìn)行低溫脅迫關(guān)聯(lián)的文獻(xiàn)及表達(dá)芯片數(shù)據(jù)搜索：(((((((((cold damage)OR cold injury) OR cold stress)OR cold harm) OR chilling damage) OR chilling injury) OR chilling stress)) AND ((((sequencing) OR expression) OR profile) OR profiling)) AND ((Arabidopsisthaliana) ORArabidopsis)。數(shù)據(jù)庫搜索的截止時(shí)間為2017年4月30日。

1.2　表達(dá)芯片數(shù)據(jù)篩選

將搜索到的文獻(xiàn)或者表達(dá)芯片數(shù)據(jù)進(jìn)行篩選，保留的低溫脅迫實(shí)驗(yàn)均符合下列標(biāo)準(zhǔn)：實(shí)驗(yàn)材料為擬南芥；同時(shí)包含低溫脅迫處理樣本和對(duì)照組樣本；基因表達(dá)測(cè)定實(shí)驗(yàn)；實(shí)驗(yàn)組與對(duì)照組間有可比較的實(shí)驗(yàn)條件；需提供可下載的表達(dá)芯片數(shù)據(jù)。具有下列任一情況的實(shí)驗(yàn)被排除于本次分析：實(shí)驗(yàn)材料非擬南芥；沒有對(duì)照組實(shí)驗(yàn)樣本；研究材料為非野生型植物材料；采用熒光定量PCR來測(cè)定基因的表達(dá)量而未進(jìn)行測(cè)序分析；綜述類研究；生物信息類研究；薈萃分析類研究。在篩選后的研究中摘錄下列信息：GEO關(guān)聯(lián)號(hào)(accession number)；文獻(xiàn)ID(reference ID)；測(cè)定平臺(tái)；脅迫處理組及對(duì)照組樣本數(shù)；研究的組織部位；低溫脅迫水平(溫度)；低溫脅迫處理時(shí)間(表1)。

注：*.相關(guān)信息未在原文中列出；UP.上調(diào)基因數(shù)(TAIR ID)；DOWN.下調(diào)基因數(shù)(TAIR)

Note: *.The information was not specified in the original research; UP. Up-regulated genes’ number(TAIR ID); DOWN. Down-regulated genes’ number(TAIR ID)

1.3　表達(dá)芯片數(shù)據(jù)預(yù)處理

在GEO數(shù)據(jù)庫下載篩選出的測(cè)序?qū)嶒?yàn)的表達(dá)數(shù)據(jù)矩陣，并對(duì)下載的數(shù)據(jù)進(jìn)行下列預(yù)處理：首先，基于不同測(cè)序平臺(tái)提供的ID轉(zhuǎn)換文件，分別將每個(gè)實(shí)驗(yàn)中的探針I(yè)D轉(zhuǎn)換成TAIR基因組ID。刪除一個(gè)探針I(yè)D對(duì)應(yīng)多個(gè)TAID基因組ID所對(duì)應(yīng)的表達(dá)數(shù)據(jù)。當(dāng)多個(gè)探針I(yè)D對(duì)應(yīng)一個(gè)TAIR基因組ID時(shí)，采用取均值的方法來處理他們的表達(dá)值。隨后檢查每個(gè)表達(dá)數(shù)據(jù)矩陣是否已進(jìn)行表達(dá)值的log2轉(zhuǎn)換。將經(jīng)過預(yù)處理的表達(dá)數(shù)據(jù)上傳到NetworkAnalyst網(wǎng)站(http://www.networkanalyst.ca)[12～13]進(jìn)行個(gè)體和批量效應(yīng)調(diào)整分析(individual and batch effect adjustment analysis)。首先對(duì)每個(gè)表達(dá)譜數(shù)據(jù)進(jìn)行個(gè)體分析，主要包含以下步驟：數(shù)據(jù)處理(pre-processing)、基因注釋(gene annotation)、歸一化分析(quantile normalization)和差異表達(dá)基因篩選分析。數(shù)據(jù)處理過程中通過檢查基因在樣本中的表達(dá)量分布刪除表達(dá)量顯著異常的基因。差異表達(dá)基因篩選分析中對(duì)照組和處理組樣本依據(jù)原實(shí)驗(yàn)中的設(shè)計(jì)。

1.4　表達(dá)芯片數(shù)據(jù)個(gè)體和批量效應(yīng)矯正分析及薈萃整合分析

利用Limma算法中的t-test分析篩選單個(gè)表達(dá)譜數(shù)據(jù)中低溫脅迫樣本和對(duì)照樣本間的差異表達(dá)基因，界限值(cutoff)為FDR(false discovery rate)矯正后的Pvalue<0.05[12～13]。隨后對(duì)經(jīng)過個(gè)體分析后的表達(dá)數(shù)據(jù)集進(jìn)行基于經(jīng)驗(yàn)貝葉斯法(empirical Bayes approach)批量效應(yīng)矯正分析，主要目的為降低不同實(shí)驗(yàn)間基因的表達(dá)率(the expression ratios)差異，使不同實(shí)驗(yàn)數(shù)據(jù)適于后續(xù)整合分析。經(jīng)過批量效應(yīng)矯正分析的8個(gè)表達(dá)譜數(shù)據(jù)被整合成一個(gè)數(shù)據(jù)集，通過數(shù)據(jù)完整度評(píng)價(jià)后，利用Cochrans’ Q test分析，評(píng)價(jià)表達(dá)數(shù)據(jù)集中不同實(shí)驗(yàn)間不一致度，并基于此分析結(jié)果，選擇random effect model模型(REM)進(jìn)行有效樣本數(shù)目整合分析，篩選整合后數(shù)據(jù)集中低溫脅迫處理組和對(duì)照組中的差異表達(dá)基因，界限值為Pvalue<0.05。采用R(version 3.2.4)平臺(tái)pheatmap軟件包進(jìn)行前25個(gè)上調(diào)和下調(diào)差異表達(dá)基因的熱圖展示。

1.5　功能富集分析

利用R平臺(tái)topGO軟件包[15]對(duì)整合后數(shù)據(jù)集中在低溫脅迫樣本和對(duì)照樣本中的差異表達(dá)基因進(jìn)行功能富集分析。在基因功能本體論數(shù)據(jù)庫(Gene ontology,http://amigo.geneontology.org)中下載擬南芥“響應(yīng)低溫脅迫(response to cold)”(GO:0009409)類目中的具有實(shí)驗(yàn)證據(jù)(experimental evidences,EXP,IDA,IPI,IMP,IGI and IEP)的注釋基因，并將上述基因作為已進(jìn)行功能鑒定的低溫脅迫響應(yīng)基因。不同家族的轉(zhuǎn)錄因子、轉(zhuǎn)錄調(diào)節(jié)因子和蛋白激酶下載于iTAK數(shù)據(jù)庫(http://bioinfo.bti.cornell.edu/cgi-bin/itak/index.cgi)。植物激素代謝途徑的參與基因下載自基因本體論數(shù)據(jù)庫。

2　結(jié)果與分析

2.1　表達(dá)芯片數(shù)據(jù)選擇

本研究通過文本篩選利用關(guān)鍵詞在GEO數(shù)據(jù)庫、PubMed數(shù)據(jù)庫和ArrayExpress數(shù)據(jù)庫分別初步篩選出1 166篇關(guān)聯(lián)文獻(xiàn)、298個(gè)測(cè)序?qū)嶒?yàn)和695個(gè)測(cè)序?qū)嶒?yàn)。通過去除三個(gè)數(shù)據(jù)庫間的冗余及非研究性實(shí)驗(yàn)，一共獲得了693個(gè)實(shí)驗(yàn)。進(jìn)一步通過材料與方法中所述的篩選條件進(jìn)行精細(xì)篩選，進(jìn)一步從59個(gè)實(shí)驗(yàn)中篩選出8組完全符合條件的含低溫脅迫處理樣本和對(duì)照樣本的表達(dá)譜測(cè)序?qū)嶒?yàn)進(jìn)行薈萃整合分析(GEO關(guān)聯(lián)號(hào)：GSE3326,GSE5620,GSE5621,GSE39090,GSE43818,GSE43819,GSE55906,GSE63186)[16～20]，其中GSE5621按照實(shí)驗(yàn)處理部位將實(shí)驗(yàn)劃分成幼苗(Shoots)和根(Roots)兩組實(shí)驗(yàn)文獻(xiàn)(表1，圖1)。篩選出的8組表達(dá)譜測(cè)序?qū)嶒?yàn)主要基于2個(gè)測(cè)序平臺(tái)：AffymetrixArabidopsisATH1 Genome Array(GPL198)和Agilent-021169Arabidopsis4 Oligo Microarray(V4)(GPL9020)；兩個(gè)平臺(tái)所含的基因數(shù)目分別為20 163個(gè)和28 401個(gè)。8組實(shí)驗(yàn)中處理組樣本數(shù)目范圍為3個(gè)到12個(gè)；對(duì)照組的樣本數(shù)目范圍為從2個(gè)到18個(gè)。8組實(shí)驗(yàn)中的植物材料有：整個(gè)地上部分、幼苗、嫩芽、葉和根。不同實(shí)驗(yàn)中低溫脅迫的處理范圍為0℃到4℃；低溫脅迫的處理時(shí)間范圍為0.5～504 h。

圖1　文本挖掘分析、研究篩選/薈萃分析及功能富集分析流程圖Fig.1　Flow chart of the literature screening,meta-analysis,enrichment strategy

2.2　個(gè)體及批量效應(yīng)矯正分析結(jié)果

8組實(shí)驗(yàn)分別進(jìn)行的個(gè)體矯正分析和差異表達(dá)基因篩選，一共篩選出12 384個(gè)在低溫脅迫處理樣本和對(duì)照樣本間差異表達(dá)的基因，其中含有7 957個(gè)上調(diào)表達(dá)基因和5 863個(gè)下調(diào)表達(dá)基因(表1)。預(yù)處理和歸一化的單個(gè)表達(dá)譜數(shù)據(jù)上傳到NetworkAnalyst網(wǎng)站進(jìn)行批量效應(yīng)矯正分析，以減少由于不同批次實(shí)驗(yàn)及其他非生物學(xué)原因造成的測(cè)序數(shù)據(jù)間的變異，利用主成分分析(Principal Component Analysis，PCA)來展示批量效應(yīng)矯正分析的效果。批量效應(yīng)矯正分析前的PCA結(jié)果顯示，8個(gè)實(shí)驗(yàn)分布離散，僅來源于同一個(gè)實(shí)驗(yàn)的樣本聚集在一起(圖2A)；在批量效應(yīng)矯正分析后的PCA結(jié)果中，可以看到8個(gè)實(shí)驗(yàn)的樣本聚集程度顯著提高(圖2B)。

2.3　薈萃整合分析及差異表達(dá)基因篩選

將8個(gè)獨(dú)立的表達(dá)譜數(shù)據(jù)整合后數(shù)據(jù)集上載到NetworkAnalyst網(wǎng)站進(jìn)行Cochran’s Q test分析。由于Cochran’s Q test分析顯示不同實(shí)驗(yàn)間有顯著的雜合度(heterogeneities)(圖2C)，因此采用REM模型進(jìn)行有效樣本數(shù)目整合分析來篩選整合后表達(dá)數(shù)據(jù)集中在脅迫組和對(duì)照組中差異表達(dá)的基因。在Pvalue<0.05界限值下，一共篩選出5 880個(gè)差異表達(dá)基因，其中包含1 553個(gè)上調(diào)基因和4 327個(gè)下調(diào)基因。5 138個(gè)薈萃整合分析后篩選出的差異表達(dá)基因在單個(gè)表達(dá)譜數(shù)據(jù)的差異表達(dá)基因篩選分析時(shí)也被鑒定為差異表達(dá)基因，有742個(gè)基因?yàn)樗C萃整合分析特異的差異表達(dá)基因(圖2D)。在1 553個(gè)薈萃分析所篩選出的上調(diào)基因中，有86個(gè)基因是已知的低溫脅迫功能基因；在4 327個(gè)下調(diào)基因中有58個(gè)是已有功能鑒定的低溫脅迫關(guān)聯(lián)基因。表2列出了前20個(gè)差異表達(dá)基因的名單；圖2E展示50個(gè)顯著上調(diào)和下調(diào)的基因的表達(dá)式樣。已報(bào)道與響應(yīng)低溫脅迫關(guān)聯(lián)的ERF/AP2轉(zhuǎn)錄因子家族中的DREB亞家族成員C-repeat/DRE binding factor 2(CBF2,AT4G25470)和C-repeat/DRE binding factor 3(CBF3,AT4G25480)被鑒定為最顯著上調(diào)的差異表達(dá)基因(表2，圖2E)。Octicosapeptide/Phox/Bem1p家族蛋白(AT3G26510)和葉綠素酶1(chlorophyllase 1,AT1G19670)的基因被鑒定為最顯著下調(diào)的基因(表2，圖2E)。

表2　薈萃分析鑒定出的前10個(gè)上調(diào)基因及前10個(gè)下調(diào)基因

圖2　A.批量效應(yīng)矯正分析前主成分分析(Principal Component Analysis，PCA)結(jié)果；B.批量效應(yīng)矯正分析后主成分分析結(jié)果；C. Cochran’s Q分析的Q-Q plot圖；D.薈萃分析篩選出的差異表達(dá)基因(Meta-DEGs)和單個(gè)分析篩選出的差異表達(dá)基因(Individual-DEGs)的venn圖；D.前25個(gè)上調(diào)基因及前25個(gè)下調(diào)基因在8組實(shí)驗(yàn)中的基因表達(dá)熱圖Fig.2　A. The principal component analysis(PCA) for datasets distributions before batch effect adjustment analysis; B. The PCA for datasets distributions after batch effect adjustment analysis; C. The Q-Q plot of Cochran’s Q test; venn diagram of differentially expressed genes identified from the meta-analysis(Meta-DEGs) and those from each individual microarray analysis(Individual-DEGs); D. Gene expression heatmap of the top 25 up-regulated genes and top 25 down-regulated genes from meta-analysis

2.4　差異表達(dá)基因的功能富集分析結(jié)果

功能富集分析用來挖掘差異表達(dá)基因集成員參與的高頻(overrepresented)生物學(xué)途徑和基因功能類目。在Pvalue<0.05的界限值下，1 553個(gè)上調(diào)基因和4 327個(gè)下調(diào)基因一共富集到了148個(gè)和193個(gè)生物學(xué)過程類目。表3展示了前15個(gè)上調(diào)的和前15個(gè)下調(diào)的生物學(xué)過程。其中顯著上調(diào)的生物學(xué)過程有：“響應(yīng)水缺乏”(GO:0009414,response to water deprivation)、“響應(yīng)冷脅迫”(GO:0009409,response to cold)、“DNA為模板的轉(zhuǎn)錄調(diào)控”(GO:0006355,regulation of transcription,DNA-templated)和“響應(yīng)脫落酸”(GO:0009737,response to abscisic acid)等。其中顯著下調(diào)的生物學(xué)過程有：“氧化還原過程”(GO:0055114,oxidation-reduction process)、“磷酸化作用”(GO:0016310,phosphorylation)、“甾醇生物合成過程”(GO:0016126,sterol biosynthetic process)和“芥子油苷的生物合成過程”(GO:0019761,glucosinolate biosynthetic process)等。

表3　功能富集分析鑒定出的前15個(gè)上調(diào)及前15個(gè)下調(diào)生物學(xué)過程

2.5　轉(zhuǎn)錄因子、激酶及植物激素關(guān)聯(lián)基因的差異表達(dá)式樣

本研究重點(diǎn)分析了轉(zhuǎn)錄因子、激酶及植物激素關(guān)聯(lián)的調(diào)控基因在整合后的表達(dá)樣本中的差異表達(dá)式樣，以發(fā)掘潛在的響應(yīng)冷脅迫的調(diào)控基因。在轉(zhuǎn)錄因子的家族成員中，本研究篩選出了26個(gè)上調(diào)和10個(gè)下調(diào)的AP2/ERF家族的轉(zhuǎn)錄因子(表4)，推測(cè)其可能與擬南芥響應(yīng)低溫脅迫有關(guān)，例如：已報(bào)道參與應(yīng)答滲透脅迫、水脅迫和熱脅迫的DREB亞家族A-2成員DRE結(jié)合蛋白2A(DREB2A，AT5G05410)和DRE結(jié)合蛋白2B(DREB2B，AT3G11020)均在冷脅迫樣本中上調(diào)。在bZIP家族中，本研究發(fā)現(xiàn)有5個(gè)上調(diào)和10個(gè)下調(diào)成員可能與響應(yīng)低溫脅迫有關(guān)(表4)，例如：與鹽脅迫和滲透脅迫響應(yīng)關(guān)聯(lián)的bZIP家族成員AT2G31370和參與調(diào)控干旱和鹽脅迫過程中鈣離子依賴蛋白激酶(Ca2+-dependent protein kinase genes)關(guān)聯(lián)的基因表達(dá)的bZIP家族成員BZIP60(AT1G42990)。

表4部分轉(zhuǎn)錄因子家族、激酶家族及激素代謝途徑中所包含的上調(diào)、下調(diào)及功能驗(yàn)證基因數(shù)目

Table4GenenumbersofDEGsincludedintheTFfamilies,PKsfamiliesandplanthormonepathways

上調(diào)基因數(shù)目Up-regulatedgenenumber下調(diào)基因數(shù)目Down-regulatedgenenumber功能驗(yàn)證基因數(shù)目Functionalverifiedgenenumber轉(zhuǎn)錄因子家族TFfamiliesTF-AP2/ERF32812TF-MYB23245TF-C2H21685TF-C2C21491TF-WRKY1353激酶家族PKsfamiliesCAMK_CAMKL-CHK1822CAMK_CDPK780RLK-Pelle_RLCK-VIIa-2730RLK-Pelle_LRR-III4120RLK-Pelle_DLSV390植物激素代謝途徑PlanthormonepathwaysResponsetoabscisicacid946871Responsetokarrikin28158Responsetoethylene433416Responsetosalicylicacid222913Responsetojasmonicacid222516

通過激酶家族成員的表達(dá)差異分析，在CAMK_CDPK家族中篩選出7個(gè)上調(diào)成員和8個(gè)下調(diào)成員，例如：參與快速防御信號(hào)傳播和溫度波動(dòng)下的轉(zhuǎn)導(dǎo)途徑的鈣調(diào)蛋白激酶5(CPK5，AT4G35310)和參與干旱和鹽脅迫響應(yīng)的鈣調(diào)蛋白激酶6(CPK6，AT2G17290)均在冷脅迫處理樣本中上調(diào)表達(dá)。促分裂素原活化蛋白激酶(mitogen-activated protein kinases，MAPK)是植物響應(yīng)非生物脅迫信號(hào)傳導(dǎo)通路中的另一個(gè)重要成員。本研究中篩選出CMGC_MAPK家族中上調(diào)成員3個(gè)和下調(diào)成員4個(gè)，例如MAP激酶11(MAP kinase 11，MPK11，AT1G01560)和MAP激酶5(MAP kinase 5，MPK5，AT4G11330)均在冷脅迫樣本中上調(diào)表達(dá)。

在植物激素關(guān)聯(lián)基因中，本研究鑒定出了與ABA途徑關(guān)聯(lián)的94個(gè)上調(diào)基因和68個(gè)下調(diào)基因(表4)，推測(cè)這些基因可能參與擬南芥寒冷脅迫信號(hào)響應(yīng)，如：脫落酸應(yīng)答元件結(jié)合因子1(ABF1,AT1G49720)及plant U-box 9(PUB9,AT3G07360)等。本研究中篩選出38個(gè)上調(diào)的和33個(gè)下調(diào)的乙烯途徑關(guān)聯(lián)的基因(表4)，推測(cè)為參與冷脅迫的潛在基因，例如：ERF家族蛋白38(ERF38，AT2G35700)和乙烯應(yīng)答元件結(jié)合因子2(ERF2，AT5G47220)。此外，在與冷脅迫關(guān)聯(lián)的茉莉酸和赤霉素等激素途徑中分別有18個(gè)上調(diào)、25個(gè)下調(diào)及12個(gè)上調(diào)、22個(gè)下調(diào)的基因被推測(cè)為冷脅迫潛在參與基因(表4)。

3　討論

挖掘低溫脅迫響應(yīng)基因或者耐受基因是促進(jìn)低溫脅迫分子機(jī)制研究和培育抗低溫脅迫作物的有效途徑之一。雖然目前已報(bào)道了近400個(gè)擬南芥低溫脅迫應(yīng)答基因，然而由于植物對(duì)低溫脅迫的響應(yīng)是一個(gè)非常復(fù)雜的體系，涉及到一系列的細(xì)胞內(nèi)信號(hào)轉(zhuǎn)導(dǎo)、基因表達(dá)調(diào)控、蛋白翻譯修飾和生理代謝調(diào)節(jié)等生化和細(xì)胞途徑[3～5]，挖掘植物響應(yīng)低溫脅迫的分子機(jī)制仍是亟待解決的關(guān)鍵問題之一。近年來,隨著分子生物學(xué)、基因工程和基因組學(xué)研究方法的快速發(fā)展,與植物抗性生理學(xué)及生態(tài)學(xué)結(jié)合的植物響應(yīng)逆境的分子機(jī)制研究已經(jīng)進(jìn)入了基因組信息時(shí)代。伴隨著基因芯片測(cè)序及高通量測(cè)序方法在植物抗逆領(lǐng)域的廣泛應(yīng)用，積累了大量與植物抗逆關(guān)聯(lián)的遺傳及表達(dá)信息資源,如何科學(xué)的整合和利用這些資源成為關(guān)鍵問題?；蛐酒C萃分析可以有效的利用上述資源，對(duì)不同實(shí)驗(yàn)的表達(dá)譜數(shù)據(jù)進(jìn)行二次分析，降低不同實(shí)驗(yàn)間的不一致性，從而獲得更加可靠的結(jié)果，挖掘出與植物抗逆性狀關(guān)聯(lián)的差異表達(dá)基因的深層次信息，并且可以對(duì)抗逆性狀關(guān)聯(lián)的分子機(jī)制作出科學(xué)的預(yù)測(cè)。近年來，整合薈萃分析在基因表達(dá)譜分析、基因診斷、藥物篩選及序列分析等諸多領(lǐng)域已呈現(xiàn)出廣闊的應(yīng)用前景[6,10～11]，但在植物基因響應(yīng)寒冷脅迫的表達(dá)數(shù)據(jù)整合方面仍有所欠缺。

本研究基于NetworkAnalyst網(wǎng)站的批量效應(yīng)校正和基于REM模型的有效樣本數(shù)目整合分析有效地降低了不同實(shí)驗(yàn)間的不一致性并進(jìn)行了科學(xué)的薈萃整合分析，并篩選出1 553個(gè)上調(diào)基因和4 327個(gè)下調(diào)基因作為不同測(cè)序?qū)嶒?yàn)響應(yīng)低溫脅迫的共有潛在基因。已有的研究表明，植物激素關(guān)聯(lián)的脅迫感受途徑、離子和激酶關(guān)聯(lián)的信號(hào)通路及轉(zhuǎn)錄因子參與的轉(zhuǎn)錄調(diào)控過程是擬南芥響應(yīng)寒冷脅迫的關(guān)鍵途徑[3,6,21]。本研究發(fā)現(xiàn)，上述三個(gè)途徑仍有潛在基因可能參與擬南芥寒冷脅迫信號(hào)響應(yīng)。例如，在植物冷脅迫應(yīng)答過程中，脅迫響應(yīng)基因表達(dá)的通路被分為ABA依賴型和ABA非依賴型[3～5]。在本研究中，鑒定出了94個(gè)上調(diào)和68個(gè)下調(diào)基因可能參與擬南芥寒冷脅迫信號(hào)響應(yīng)，如：ABA信號(hào)通路及蛋白泛素化途徑關(guān)聯(lián)的基因等。此外，乙烯應(yīng)答因子ERF是植物啟動(dòng)冷脅迫關(guān)聯(lián)基因表達(dá)的重要轉(zhuǎn)錄因子，主要通過其保守的DNA結(jié)合域與乙烯應(yīng)答元件GCCbox相互作用，以啟動(dòng)脅迫應(yīng)答基因的表達(dá)。本研究通過整合分析預(yù)測(cè)出的數(shù)十個(gè)上調(diào)及下調(diào)的ERF家族轉(zhuǎn)錄因子，推測(cè)此家族仍有潛在的冷脅迫應(yīng)答參與基因待通過功能性實(shí)驗(yàn)進(jìn)行驗(yàn)證。

除了上述調(diào)控途徑，本研究基于功能富集分析發(fā)現(xiàn)“細(xì)胞內(nèi)信號(hào)傳遞(intracellular signal transduction)”、“細(xì)胞內(nèi)水缺乏應(yīng)答(cellular response to water deprivation)”、“晝夜節(jié)律(circadian rhythm)”和“防御反應(yīng)(defense response)”等生物學(xué)過程在冷脅迫響應(yīng)過程中的功能仍有待進(jìn)一步挖掘。此外，由于脫水、高鹽及冷脅迫均會(huì)不同程度的影響滲透平衡等內(nèi)穩(wěn)態(tài)，植物響應(yīng)不同非生物脅迫轉(zhuǎn)錄因子調(diào)控機(jī)制有一定的相似性[7,22～23]，因此本研究挖掘的與其他脅迫關(guān)聯(lián)的功能基因及其參與的生物學(xué)過程可能參與低溫脅迫應(yīng)答，例如：“鹽脅迫應(yīng)答的負(fù)調(diào)控(negative regulation of response to salt stress)”“氧化脅迫應(yīng)答(response to oxidative stress)”和“干旱恢復(fù)(drought recovery)”等。雖然上述潛在的生物學(xué)途徑關(guān)聯(lián)的基因在植物冷脅迫響應(yīng)過程中的功能仍需功能實(shí)驗(yàn)進(jìn)行驗(yàn)證，但本研究通過薈萃分析為植物冷脅迫響應(yīng)的潛在功能基因和潛在分子機(jī)制的挖掘提供了更全面和可靠的候選研究對(duì)象。

雖然本研究通過薈萃分析挖掘出了參與擬南芥低溫脅迫響應(yīng)的潛在基因及其關(guān)聯(lián)的功能子過程，仍需要指出本研究的優(yōu)點(diǎn)及局限性。本研究的局限性主要是由所篩選的不同測(cè)序?qū)嶒?yàn)之間的差異性造成的：我們選擇的8個(gè)實(shí)驗(yàn)包括了5個(gè)組織部位，2個(gè)低溫脅迫條件，8個(gè)時(shí)間梯度和2個(gè)不同的測(cè)序平臺(tái)。首先，不同研究間由于單個(gè)實(shí)驗(yàn)操作技術(shù)、分子遺傳因素、環(huán)境條件變化等造成的不一致性是很難移除的。其次，由于基因在不同組織部位表達(dá)差異造成的不一致性也是很難評(píng)價(jià)和移除的。上述兩個(gè)因素構(gòu)成了本研究的主要局限性。然而，批量效應(yīng)校正方法被設(shè)計(jì)用來通過平衡基因的表達(dá)波動(dòng)來降低不同實(shí)驗(yàn)間的差異性，再進(jìn)一步結(jié)合基于REM模型的有效樣本數(shù)目整合分析，可以進(jìn)一步消減不同表達(dá)測(cè)序?qū)嶒?yàn)間基因表達(dá)值的不一致性[12～13]。此外，采用GO功能途徑基因功能富集分析可以在差異表達(dá)基因的基礎(chǔ)上查找與擬南芥冷脅迫響應(yīng)關(guān)聯(lián)的子過程，潛在的功能基因和子過程共同構(gòu)成了本研究預(yù)測(cè)的擬南芥響應(yīng)低溫脅迫的分子基礎(chǔ)。植物對(duì)低溫脅迫響應(yīng)是一個(gè)非常復(fù)雜的體系，包括信號(hào)轉(zhuǎn)導(dǎo)、基因調(diào)控、蛋白修飾和代謝調(diào)節(jié)等生理和生化途徑，基于本研究提供的潛在分子基礎(chǔ)，利用遺傳學(xué)、生物化學(xué)、基因組學(xué)、轉(zhuǎn)錄組學(xué)技術(shù)相結(jié)合，開展擬南芥及其他植物對(duì)低溫響應(yīng)的研究，將更深入地闡明植物應(yīng)答低溫脅迫的分子機(jī)制，并在此基礎(chǔ)上研發(fā)抗冷作物。

1.簡(jiǎn)令成,盧存福,李積宏,等.適宜低溫鍛煉提高冷敏感植物玉米和番茄的抗冷性及其生理基礎(chǔ)[J].作物學(xué)報(bào),2005,31(8):971-976.

Jian L C,Lu C F,Li J H,et al.Increment of chilling tolerance and its physiological basis in chilling-sensitive corn sprouts and tomato seedlings after cold-hardening at optimum temperatures[J].Acta Agronomica Sinica,2005,31(8):971-976.

2.Lukatkin A S.Contribution of oxidative stress to the development of cold-induced damage to leaves of chilling-sensitive plants:2.The activity of antioxidant enzymes during plant chilling[J].Russian Journal of Plant Physiology,2002,49(6):782-788.

3.烏鳳章,王賀新,徐國(guó)輝,等.木本植物低溫脅迫生理及分子機(jī)制研究進(jìn)展[J].林業(yè)科學(xué),2015,51(7):116-128.

Wu F Z,Wang H X,Xu G H,et al.Research progress on the physiological and molecular mechanisms of woody plants under low temperature stress[J].Scientia Silvae Sinicae,2015,51(7):116-128.

4.Bloom A J,Zwieniecki M A,Passioura J B,et al.Water relations under root chilling in a sensitive and tolerant tomato species[J].Plant,Cell & Environment,2004,27(8):971-979.

5.劉勇,馮美利,曹紅星,等.低溫脅迫對(duì)油棕葉片養(yǎng)分含量變化的影響[J].熱帶農(nóng)業(yè)科學(xué),2014,34(6):16-19.

Liu Y,Feng M L,Cao H X,et al.Effect of low temperature on the leaf nutrient content of oil palm[J].Chinese Journal of Tropical Agriculture,2014,34(6):16-19.

6.Song B H,Mitchell-Olds T.Evolutionary and ecological genomics of non-model plants[J].Journal of Systematics and Evolution,2011,49(1):17-24.

7.Catalá R,Medina J,Salinas J.Integration of low temperature and light signaling during cold acclimation response inArabidopsis[J].Proceedings of the National Academy of Sciences of the United States of America,2011,108(39):16475-16480.

8.Atkinson N J,Lilley C J,Urwin P E.Identification of genes involved in the response ofArabidopsisto simultaneous biotic and abiotic stresses[J].Plant Physiology,2013,162(4):2028-2041.

9.Yu G W,Zang W D,Yang X Q,et al.Meta-analysis ofArabidopsisthalianaunder abscisic acid and salt stress[J].Journal of Medicinal Plants Research,2011,5(24):5889-5893.

10.Zhang L Y,Zhang X J,Fan S J.Meta-analysis of salt-related gene expression profiles identifies common signatures of salt stress responses inArabidopsis[J].Plant Systematics and Evolution,2017,303(6):757-774.

11.Tseng G C,Ghosh D,Feingold E.Comprehensive literature review and statistical considerations for microarray meta-analysis[J].Nucleic Acids Research,2012,40(9):3785-3799.

12.Xia J G,Benner M J,Hancock R E.NetworkAnalyst-integrative approaches for protein-protein interaction network analysis and visual exploration[J].Nucleic Acids Research,2014,42(W1):W167-W174.

13.Xia J G,Gill E E,Hancock R E.NetworkAnalyst for statistical,visual and network-based meta-analysis of gene expression data[J].Nature Protocols,2015,10(6):823-844.

14.Jha P K,Vijay A,Sahu A,et al.Comprehensive Gene expression meta-analysis and integrated bioinformatic approaches reveal shared signatures between thrombosis and myeloproliferative disorders[J].Scientific Reports,2016,6:37099.

15.Alexa A,Rahnenfuhrer J.topGO:enrichment analysis for gene ontology.[http://www.bioconductor.org/packages/2.11/bioc/html/topGO.html].

16.Lee B H,Henderson D A,Zhu J K.TheArabidopsiscold-responsive transcriptome and its regulation by ICE1[J].The Plant Cell,2005,17(11):3155-3175.

17.Kilian J,Whitehead D,Horak J,et al.The AtGenExpress global stress expression data set:protocols,evaluation and model data analysis of UV-B light,drought and cold stress responses[J].The Plant Journal,2007,50(2):347-363.

18.Guan Q M,Wu J M,Zhang Y Y,et al.A DEAD box RNA helicase is critical for pre-mRNA splicing,cold-responsive gene regulation,and cold tolerance inArabidopsis[J].The Plant Cell,2013,25(1):342-356.

19.Kim Y,Park S,Gilmour S J,et al.Roles of CAMTA transcription factors and salicylic acid in configuring the low-temperature transcriptome and freezing tolerance ofArabidopsis[J].The Plant Journal,2013,75(3):364-376.

20.Chan Z L,Wang Y P,Cao M J,et al.RDM4 modulates cold stress resistance inArabidopsispartially through the CBF-mediated pathway[J].New Phytologist,2016,209(4):1527-1539.

21.Gorsuch P A,Sargeant A W,Penfield S D,et al.Systemic low temperature signaling inArabidopsis[J].Plant & Cell Physiology,2010, 51(9) :1488-1498.

22.Agarwal P K,Shukla P S,Gupta K,et al.Bioengineering for salinity tolerance in plants:state of the art[J].Molecular Biotechnology,2013,54(1):102-123.

23.Munns R.Genes and salt tolerance:bringing them together[J].New Phytologist,2005,167(3):645-663.