杜文麗 陳中釤 許端祥 徐同偉 高 山 溫慶放

(1福州市蔬菜科學研究所，福建 福州 350111；2福建省農(nóng)業(yè)科學院作物所，福建 福州 350013)

苦瓜(Momordica charantiaL.)又名涼瓜，起源于東印度，現(xiàn)主要種植于熱帶、亞熱帶和溫帶區(qū)域，加勒比海和南美地區(qū)也有栽培。近年來，人們逐步認識到苦瓜的營養(yǎng)價值及保健功能。我國南方地區(qū)生產(chǎn)的苦瓜除供本地和北運外，還遠銷東南亞地區(qū)，己成為南方部分地區(qū)出口創(chuàng)匯的主要蔬菜之一[1-2]?？喙蠈傧矞厥卟?，在實際生產(chǎn)中，苦瓜苗期常常遭受早春寒流或秋季氣溫驟降的侵襲，10℃以下的冷害導致生長緩慢甚至形成僵苗，嚴重制約苦瓜種植范圍的擴大和上市時節(jié)的調(diào)節(jié)，生產(chǎn)上亟需選育耐低溫且綜合性狀好的新品種[3]。

與植物耐低溫相關(guān)的研究逐漸由細胞學、生理生化轉(zhuǎn)向分子生物學，且植物耐低溫分子機制的研究已取得了長足進步，現(xiàn)已有100 余種轉(zhuǎn)基因耐低溫植物，共涉及50 多個物種[4-5]。目前，低溫脅迫涉及一系列復雜的信號傳導網(wǎng)絡，誘導植物體內(nèi)相關(guān)基因的表達，但對于植物將低溫脅迫信號傳遞到細胞核內(nèi)激活轉(zhuǎn)錄因子，以及激活下游基因表達使機體出現(xiàn)耐低溫特性的機制仍不清楚。前人研究發(fā)現(xiàn)CBF1 通過結(jié)合CRT/DRE(c-repeat/dehydration responsive element)順式作用元件，調(diào)控耐低溫特性基因表達來參與植物中低溫脅迫下的信號傳導，這種通過轉(zhuǎn)錄因子所調(diào)控的信號路在學術(shù)界已被廣泛驗證[6-7]。有研究學者先后在擬南芥中獲得CBF2、CBF3 和CBF4 基因轉(zhuǎn)入擬南芥中，發(fā)現(xiàn)擬南芥植株耐低溫脅迫的能力得到顯著提升，且導致許多與低溫脅迫相關(guān)的生理生化指標發(fā)生變化，如脯氨酸和可溶性糖含量增加等[8]。周麗霞等[9]發(fā)現(xiàn)WRKY7、WRKY1 等基因?qū)τ妥氐哪偷蜏匦园l(fā)揮調(diào)控作用。李佳[10]研究發(fā)現(xiàn)馬鈴薯在4℃處理后 St2RMYB80 基因表達量下調(diào)，St4RMYB1 和St3RMYB2 基因表達量上調(diào)。

第二代測序技術(shù)已廣泛應用于轉(zhuǎn)錄組測序(RNAseq)領(lǐng)域，通過比對植物生物和非生物脅迫轉(zhuǎn)錄組基因reads 數(shù)獲得植物體大量基因信息，能夠在分子水平上了解基因表達量的變化，揭示特定的生物學過程和分子機制[11-13]。Chen 等[13]通過Illumina nextseq 2500 測序平臺對低溫脅迫下的甘蔗葉片進行研究，發(fā)現(xiàn)2 583 個基因上調(diào)表達，3 302 個基因下調(diào)表達。黃超[11]利用Illumina Hiseq-2000 測序平臺對4℃低溫處理0.8 和32 h 的甜瓜葉片進行測序，篩選出19個包括編碼WRKY、b HLH、NAC、AP2/ERF等轉(zhuǎn)錄因子在內(nèi)可能與耐低溫具有潛在相關(guān)的候選基因。此外，轉(zhuǎn)錄組測序技術(shù)已應用于大白菜、花椰菜、黃瓜等蔬菜作物研究中[14-17]。目前，關(guān)于苦瓜在低溫脅迫方面的研究還限于表型、生理生化指標鑒定和遺傳規(guī)律等方面[18-19]，通過RNA-seq 揭示苦瓜的耐低溫機制方面研究鮮有報道，在一定程度上約束了苦瓜育種進程。

本研究以福州市蔬菜科學研究所選育的較耐低溫苦瓜自交系43 為研究對象，依托百邁客Illumina HiSeq2500 平臺對26℃(對照)和模擬福建早春夜溫8℃氣候[20]處理12 h 的兩個樣品葉片進行轉(zhuǎn)錄組測序分析，再通過生理指標分析和分子手段對苦瓜響應低溫脅迫的分子機制進行探討，同時進一步挖掘苦瓜耐低溫相關(guān)的候選基因，以期為苦瓜耐低溫基因工程提供一定的理論支撐。

1 材料與方法

1.1 試驗材料

供試材料為耐低溫苦瓜自交系43，由福州市蔬菜科學研究所苦瓜課題組選育。選取飽滿健壯的苦瓜種子采用溫湯浸種法催芽，待種子露白后轉(zhuǎn)移至營養(yǎng)缽中(椰糠＋草木灰)中，隨后放置于26℃培養(yǎng)箱(立德泰勒科學儀器有限公司，型號:LT-ACC)中培養(yǎng)，光強8 000 lx，相對濕度70%)。約2 周后苦瓜苗生長至四葉一心時，將培養(yǎng)箱溫度調(diào)至8℃，選取生長良好、發(fā)育狀態(tài)一致的苦瓜苗進行低溫處理，處理前澆透水。觀察每個時期的植株形態(tài)，其他培養(yǎng)條件同上[20]。處理0(對照)、4、8、12 h 后剪取幼苗真葉測定相關(guān)生理指標。剪取處理后剪取12 h 幼苗真葉用錫箔紙包好，以0 h 為對照，立即放置于液氮中速凍后送北京百邁客生物科技有限公司進行RNA-seq。

1.2 生理指標測定方法

硫代巴比妥酸(thiobarbituric acid，TAB)法用于測定丙二醛(malondialdehyde，MDA)含量。酸性茚三酮法用于測定游離脯氨酸(proline，Pro)含量。電解質(zhì)滲出率(rate of electrolyte efflux，REC)的測定參考文獻[21]。

1.3 RNA-seq 方法

總RNA 按照提取后經(jīng)消化處理后及時對RNA 樣品的純度、濃度和完整性進行檢測，保證使用高質(zhì)量的RNA 進行轉(zhuǎn)錄組測序。依托Illumina HiSeq2500 高通量測序平臺進行上機測試工作，PE150 為測序讀設置值。利用FASTX 軟件對獲得的數(shù)據(jù)進行質(zhì)量控制，即去除Raw Data 中的接頭和低質(zhì)量序列后得到Clean Data，再使用Trinity 軟件進行組裝最終獲得Unigene序列。利用BLAST 程序(設置E-value= 10-5)使用NR、COG、KOG 及KEGG 等數(shù)據(jù)庫對Unigene 進行注釋及和分類?；虮磉_衡量指標采用FPKM 法[22]，使用DESeq[23]對樣品間差異表達進行分析，以樣品間表達量符合錯誤發(fā)現(xiàn)率數(shù)≥2 且≤0.5 作為篩選標準獲得差異表達基因。最后對差異表達的基因進行GO 功能注釋、COG 功能注釋和富集分析。

1.4 實時熒光定量PCR 檢測

參照植物RNA 提取試劑盒(RNAprep Pure Plant Kit，天根生化科技有限公司)獲得苦瓜葉片的總RNA，然后參照[Prime ScriptTM RT reagent Kit，寶生物工程(大連)有限公司]試劑盒說明反轉(zhuǎn)錄合成cDNA。用Primer 5.0 軟件設計與低溫逆境相關(guān)的苦瓜差異表達基因引物，以McCYP(GeneBank 登錄號:HQ171897)作為內(nèi)參基因，候選基因的引物序列見表1，引物由北京六合華大基因科技股份有限公司合成。實時熒光定量PCR(quantitative real-tirne PCR，qRTPCR)采用ABI 7500 Fast Real-Time PCR System(美國)進行PCR 擴增，每次試驗設3 次重復。反應體系參照[Super Real Pre Mix Plus (SYBR Green)，天根生化科技有限公司)]試劑盒說明書進行，數(shù)據(jù)采用2-ΔΔCt方法計算[24]。

1.5 數(shù)據(jù)處理

采用SPSS 16 軟件對試驗數(shù)據(jù)進行統(tǒng)計，利用Excel2016.1.0 軟件作圖。

2 結(jié)果與分析

2.1 低溫脅迫下苦瓜的表型與生理生化指標變化

苦瓜經(jīng)低溫脅迫4 h 后，苦瓜的真葉就開始出現(xiàn)萎蔫失水的現(xiàn)象，這是植物受到低溫脅迫的典型反應；低溫脅迫8 h 后，葉柄軟化、向下彎曲，未成熟莖出現(xiàn)更嚴重的萎蔫現(xiàn)象(圖1-A)。但在低溫脅迫12 h 后，又出現(xiàn)明顯的恢復現(xiàn)象，表明苦瓜43 號自交系具有良好的低溫耐受能力，適于作為本研究的材料。

圖1 低溫脅迫前后苦瓜表型變化和生理生化變化Fig.1 Phenotypic and physiological changes in cold-stressed bitter gourd

由圖1-B 可知，在經(jīng)過低溫脅迫后MDA 含量先升高后降低，在低溫脅迫8 和12 h 后，比0 h 分別增加了96.80%和60.34%，峰值約為對照的2 倍，說明此時苦瓜的細胞膜受到損傷導致MDA 含量增加，但在低溫脅迫12 h 后損傷有所恢復。由圖1-C 可知，隨著低溫脅迫時間延長，Pro 含量整體呈現(xiàn)先升高后降低的趨勢，在低溫脅迫4、8 和12 h 時，Pro 含量相比于對照分別提高了114.47%、125.20%和87.07%，這與MDA 含量的變化情況相一致。由圖1-D 可知，REC 隨脅迫時間的延長而持續(xù)增加并保持在較高水平，在低溫脅迫4、8 和12 h 時，比對照分別增加了7.55、9.38 和12.47個百分點。

2.2 測序數(shù)據(jù)整體評估及差異表達基因的獲得

為深入探索苦瓜在轉(zhuǎn)錄水平上對低溫脅迫的應答，剪取8℃低溫條件下處理0 和12 h 的苦瓜葉片，采用Illumiana HiSeq2500 平臺進行轉(zhuǎn)錄組測序，共獲得1.77 Gb Clean Data，Q30 值均大于91.00%(表2)，GC含量均在48%的水平。Clean Data 經(jīng)過濾、拼接和組裝后獲得的Unigene 數(shù)量為89 757 條，N50 長度為1 239 bp，平均長度達到700.75 bp (表3)。Unigene序列數(shù)量最多的分布在200 ～300 bp 之間，在500 ～1 000 bp 之間次之(42.64%)，Unigene 長度達到1 000 bp 以上的占18.2%以上，共計16 329 條。以上結(jié)果表明RNA-seq 結(jié)果符合下一步生物信息學分析的要求。

表2 測序數(shù)據(jù)結(jié)果統(tǒng)計Table 2 Statistic of evaluating the sequencing data

表3 測序數(shù)據(jù)結(jié)果統(tǒng)計Table 3 Statistic of evaluating the sequencing data

將得到的Unigene 序列比對結(jié)果經(jīng)NR、KEGG等多個數(shù)據(jù)平臺得到對應的注釋信息(表4)。在預設的E 值范圍最終有48 662 個注釋信息的Unigene，占全部Unigenes 的54.2%。比對得到的Unigene 長度大于1 000 nt 的達到14 964 個，占所有得到注釋Unigene 的30.5%。但仍有41 095 個Unigenes 在數(shù)據(jù)庫中沒有得到注釋。通過對差異表達分析結(jié)果進行統(tǒng)計，在低溫脅迫處理后共獲得1 285 個差異表達基因，其中916 個基因上調(diào)、369個基因下調(diào)。

表4 測序數(shù)據(jù)結(jié)果統(tǒng)計Table 4 Statistic of evaluating the sequencing data

2.3 差異表達基因GO 分類

對差異表達基因進行GO 聚類分析聚類，結(jié)果有1 495 條Unigenes 映射到GO 不同功能節(jié)點上，顯著富集主要包括在細胞組分、分子功能、生物過程中的50個亞類(圖2)。其中最大的主類群是參與生物學過程的差異表達基因，共有749 個；歸入細胞組分和分子功能類的分別為347 和399 個；細胞部分、細胞器和膜結(jié)構(gòu)在細胞組分中是排名前三的亞類，分別為123、62 和46 個差異表達基因。在分子功能主類群，大多數(shù)的差異表達基因歸屬于催化活性(177 個)和結(jié)合活性(149個)。在生物學過程類群，由191 條Unigene 組成的代謝進程是數(shù)量最多的亞類，其次是細胞過程(174 個)和單一生物過程(134 個)。生物過程分組中共富集到14 個GO 條目，差異基因所占比例較大的有代謝過程、細胞過程、單個有機體過程、應激反應、生物調(diào)節(jié)。在細胞組分分組中共富集到8 個GO 條目，差異表達基因所占比例較高的有細胞部分細胞、膜、細胞器、膜結(jié)構(gòu)。在分子功能分組中共富集11 個GO 條目，催化活性、結(jié)合活性占較大比例。這些大量積累的基因表明，苦瓜在遭受低溫脅迫過程中進行著復雜的新陳代謝和酶促反應。

2.4 差異表達基因的COG 分類

圖2 苦瓜轉(zhuǎn)錄組GO 功能分類Fig.2 Gene ontology classification of Momordica charantia L.transcriptome

將轉(zhuǎn)錄組獲得的差異表達基因與COG 數(shù)據(jù)庫進行比對，對其所編碼的蛋白質(zhì)進行直系同源分類，進一步預測基因功能并進行分類統(tǒng)計(表5)。共得到了239 個COG 功能注釋的信息，按其功能劃分為24 類，其中有3 類功能類型在本研究中無Unigene 比對。其中，僅一般功能預測(16.74%)和轉(zhuǎn)錄(10.04%)類所占的比例最大，隨后依次是次生代謝物合成，運輸和代謝(9.62%)及重復，重組和修飾(8.79%)。

2.5 轉(zhuǎn)錄組差異表達基因代謝途徑分析

將差異基因比對到KEGG 數(shù)據(jù)庫獲得基因注釋結(jié)果以進一步探索基因功能(表6)。結(jié)果顯示共有237 個差異基因歸入94 種代謝路徑。其中參與核糖體代謝途徑的差異表達基因有19 個；其次有18 個差異表達Unigene 與苯丙素生物合成相關(guān)；與苯丙氨酸代謝路徑相關(guān)的差異表達基因占第三位；此外，參與植物激素信號傳導以及淀粉和糖類代謝途徑的差異表達基因均為11 個。與脯氨酸合成相關(guān)的精氨酸和脯氨酸代謝通路與低溫脅迫相關(guān)的氧化磷酸化(oxidative phosphorylation)通路放入差異表達基因分別獲得1 個和2 個。代謝途徑包含碳水化合物、信號傳導、光合等涉及生命活動的各個階段，表明這些代謝途徑在苦瓜遭受低溫脅迫時可能發(fā)揮著重要作用。

2.6 苦瓜耐低溫相關(guān)差異表達候選基因

為了挖掘參與耐低溫脅迫信號網(wǎng)絡相關(guān)基因，在差異表達基因中篩選到了包括WRKY、ERF/AP2(乙烯應答)、ZAT轉(zhuǎn)錄因子等可能參與苦瓜耐低溫脅迫的10 個重要基因(表7)。10 個候選基因的序列和相似性分析見附表1。

2.7 qRT-PCR 分析

為了驗證RNA-seq 數(shù)據(jù)的準確性并進一步分析苦瓜候選基因的耐低溫機制，對挑選出的10 個候選基因利用qRT-PCR 觀測低溫脅迫下動態(tài)變化情況，結(jié)果如圖3 所示。WRKY33 的表達量在脅迫4 h 后急劇上升，8 h 時稍下調(diào)隨后又有所上升，但都始終高于對照(0 h)；MYB108 的表達量隨脅迫時間的延長呈逐步上升趨勢；ERF017 的表達量隨低溫脅迫時間的延長呈現(xiàn)先下降后上升的趨勢；ERF021 和ERF053 的表達量在受到低溫脅迫4 h 時急劇上升，且均在低溫脅迫處理8 h 時達到峰值，在低溫脅迫處理12 h 時稍有下降但仍高于對照；NAC25 的表達量在受到低溫脅迫后下降，低溫脅迫處理8 h 表達量上升后又下調(diào)；Ga2＋ATPase和Ga2＋ATPase2 的表達量隨著低溫脅迫時間的延長均呈現(xiàn)先上升后下降的趨勢；ZAT11 和ZAT12的表達量在受到低溫脅迫后先上調(diào)后下降，且峰值均出現(xiàn)在低溫脅迫8 h，低溫脅迫12 h 時仍高于對照。所有基因的表達變化與RNA-seq 所得的變化情況一致，說明轉(zhuǎn)錄組結(jié)果數(shù)據(jù)可靠。這些基因表達量的上調(diào)或下調(diào)也預示著它們在苦瓜耐低溫脅迫中可能發(fā)揮著重要作用。

表5 差異表達基因的COG 注釋分類情況Table 5 Classification of COG annotation about DEGs

3 討論

苦瓜經(jīng)低溫脅迫后出現(xiàn)萎蔫失水的現(xiàn)象，隨后又明顯恢復，表明該苦瓜材料具有良好的低溫耐受能力。MDA 是評價植物細胞膜系統(tǒng)損傷程度的重要指標，本研究中苦瓜受到一定程度的冷害后MDA 含量先升高后降低，說明植物損傷后逐漸得到恢復，這與在低溫脅迫下甜瓜不同品種間的耐低溫性與MDA 含量成顯著負相關(guān)的研究結(jié)果一致[10]。Pro 是植物細胞內(nèi)重要的滲透調(diào)節(jié)物質(zhì)，通過提高細胞持水能力從而提高植物的耐低溫特性[1，25]。本研究中Pro 含量先上升后下降，說明在低溫脅迫前期，苦瓜體內(nèi)能夠快速促進Pro含量的提升使植物細胞膜結(jié)構(gòu)免受損傷。MDA 和Pro的變化情況與苦瓜的表型變化一致，說明苦瓜在低溫脅迫時可能刺激相關(guān)通路基因的表達以保護細胞膜系統(tǒng)不受損害。電解質(zhì)參透率能反映植物在逆境條件下細胞膜透性程度[1]。本研究中苦瓜葉片的電解質(zhì)參透率在低溫脅迫的過程中始終維持較高水平，表明低溫造成的膜損傷并沒有在12 h 內(nèi)得到完全修復。

表6 苦瓜部分差異表達Unigene 的代謝通路Table 6 Some metabolic pathways involving DEG of Momordica charantia L

表7 差異表達中耐低溫脅迫候選基因Table 7 Low temperature candidate gene of the DGEs

圖3 低溫脅迫下苦瓜候選基因的相對表達量Fig.3 Relative expression of candidate gene under cold stress in Momordica charantia L.

通過對低溫脅迫條件下苦瓜葉片進行RNA-seq分析，共獲得1 285 個差異表達基因，受低溫誘導表達的基因數(shù)目明顯多于低溫抑制的基因數(shù)目。差異表達基因GO 富集到的大量基因存在于代謝、催化活性和應激反應亞類中，表明苦瓜在低溫脅迫過程中進行著復雜的新陳代謝和應激反應。差異表達基因歸類到物質(zhì)代謝、植物激素信號傳導、光合作用相關(guān)蛋白及植物晝夜節(jié)律等通路，表明苦瓜可能通過改變信號傳導、環(huán)境應激及次級代謝物合成來參與低溫脅迫的響應，這與黃瓜、番茄、葡萄通過轉(zhuǎn)錄組測序技術(shù)研究低溫響應相關(guān)差異表達基因的主要通路結(jié)果類似[7，26-27]。研究表明WRKY、MYB、AP2/ERF、ZAT和NAC這幾類轉(zhuǎn)錄因子在調(diào)節(jié)植物低溫脅迫應答網(wǎng)絡中起著關(guān)鍵性的調(diào)控作用[28-30]。對篩選出表達量變化較為明顯的幾類轉(zhuǎn)錄因子進行qRT-PCR 驗證并對這些基因在苦瓜低溫脅迫過程進行動態(tài)跟蹤，發(fā)現(xiàn)WRKY33 基因表達量先上調(diào)后下降隨后又上升，MYB108 表達量持續(xù)上升、ZAT11 和ZAT12 上調(diào)表達后下調(diào)，與前人研究結(jié)果也相一致[31]。ERF021 和ERF053 上調(diào)表達，之后出現(xiàn)下降趨勢，ERF017 表達量隨著脅迫時間的延長先略下調(diào)后上調(diào)，這與Du 等[32]證實在冷脅迫后期大量ERF上調(diào)表達的研究結(jié)果一致。NAC25 表達量在低溫脅迫后劇烈下降，預示其可能負調(diào)控植物低溫應答。Ca2＋作為第二信使在植物感應外界逆境和信號傳遞過程中發(fā)揮重要的作用，Ca2＋ATPase的功能研究證明該家族基因響應各種逆境脅迫[33-34]，本研究Ca2＋ATPase基因表達受低溫誘導明顯，說明其可能在苦瓜應對低溫逆境時發(fā)揮重要作用。qRT-PCR 和轉(zhuǎn)錄組測序基因表達的上調(diào)和下調(diào)一致，具體的定量值有差異可能是機器或者不同批次處理植株中取樣引起的誤差，表明轉(zhuǎn)錄組測序的結(jié)果真實可信，基因表達和生理變化能夠相互印證。此外，發(fā)現(xiàn)脯氨酸合成、氧化磷酸化過程通路中3 個差異基因表達被激活，這與生理生化測試結(jié)果一致，可作為苦瓜耐低溫基因挖掘的補充資源。

由本研究結(jié)果推測在苦瓜耐低溫脅迫中Ca2＋作為信號分子激活相關(guān)轉(zhuǎn)錄因子的表達，被激活的轉(zhuǎn)錄因子再調(diào)控下游功能基因的表達及修飾，減弱了低溫脅迫對植物細胞膜產(chǎn)生的破壞，促使植物體內(nèi)失衡的活性氧和激素水平快速恢復，從而提高了植物的耐低溫性。本研究篩選出的候選基因可能參與苦瓜機體耐低溫脅迫過程，為下一步開展苦瓜耐低溫相關(guān)分子機制研究提供了一定的理論支持。

4 結(jié)論

本研究以耐低溫品種自交系43 號為材料對低溫冷脅迫下苦瓜進行形態(tài)和相關(guān)生理指標測定及分析，并對26、8℃處理下的苦瓜葉片進行RNA-seq 獲得了豐富的轉(zhuǎn)錄本信息，對差異表達基因的基因功能進行注釋，了解苦瓜葉片在低溫脅迫下生理變化過程，篩選出包括WRKY、NAC、AP2/ERF等轉(zhuǎn)錄因子家族候選基因，可能參與苦瓜機體耐低溫脅迫過程調(diào)控。10 個差異表達基因的qRT-PCR 動態(tài)分析結(jié)果與RNA-seq 變化一致，說明轉(zhuǎn)錄組數(shù)據(jù)可靠，基因表達變化和生理指標的變化也一致。本研究利用轉(zhuǎn)錄組測序技術(shù)對低溫脅迫下的苦瓜材料進行測序，通過生物信息學方法對苦瓜耐低溫性的多種代謝途徑進行分析，并從分子水平闡述了苦瓜低溫脅迫下的響應機制，為進一步研究苦瓜耐低溫脅迫分子生物學研究功能奠定了一定的理論基礎(chǔ)。

附表1: 10 個候選基因相似性比較結(jié)果Schedule 1: Similarity comparison results of 10 candidate genes