黃長兵，程培蕾，楊紹宗，張煥朝，姜正之，金立敏

(1.南京林業(yè)大學(xué) 南方現(xiàn)代林業(yè)協(xié)同創(chuàng)新中心/林學(xué)院，江蘇 南京 210037； 2.蘇州農(nóng)業(yè)職業(yè)技術(shù)學(xué)院 園藝科技學(xué)院，江蘇 蘇州 215008； 3.蘇州市華冠園創(chuàng)園藝科技有限公司，江蘇 蘇州 215557)

萱草(Hemerocallisfulva)又被稱為蘐草、忘憂等，屬于百合科萱草屬多年生植物，其根莖較短但其紡錘形肉質(zhì)根較厚[1-3]，在中國有著悠久的栽培歷史。萱草既可以作為園林植物作觀賞之用，也可以作為蔬菜食用，亦可以作為藥物治??；由于其品種多樣，色彩豐富，抗逆性強(qiáng)等因素，成為道路綠化、公園綠地、屋頂花園等生態(tài)城市建設(shè)廣泛使用的園林綠化材料[4]。作為萱草屬中較少幾種可食用植物，萱草擁有較高的經(jīng)濟(jì)價(jià)值。此外，萱草中富含大量次級代謝產(chǎn)物，這些代謝物對于人體的焦慮、腫脹等癥狀具有良好的治療效果[5-6]。目前，研究者們已經(jīng)對萱草的生理生化、分類鑒定、組織培養(yǎng)和育種等領(lǐng)域進(jìn)行了較多研究[7-11]，但有關(guān)于萱草逆境脅迫的報(bào)道和研究相對較少。低溫脅迫作為最常見的非生物脅迫因素之一，能夠?qū)е轮参锏纳L和發(fā)育出現(xiàn)延遲，同時(shí)也會(huì)影響作物的產(chǎn)量[12-14]。冷害是影響萱草觀賞價(jià)值和經(jīng)濟(jì)價(jià)值的重要因素，因此，研究萱草在低溫脅迫下基因的表達(dá)變化對于解決萱草的冷害具有重要意義。

轉(zhuǎn)錄組測序與生物信息學(xué)分析技術(shù)相結(jié)合可以大大促進(jìn)功能基因組學(xué)的研究，通過測序得到一定條件下細(xì)胞的整體轉(zhuǎn)錄信息，為從轉(zhuǎn)錄水平來發(fā)掘生物體的生長發(fā)育、次級代謝與調(diào)控機(jī)制[15-16]提供了便利。近些年，研究者們用轉(zhuǎn)錄組分析策略來研究植物抗低溫脅迫基因組學(xué)，已經(jīng)獲得擬南芥[17]、番茄[18]、紫花苜蓿[19]、黃瓜[20]和水稻[21]等眾多植物轉(zhuǎn)錄組數(shù)據(jù)，這些都為闡明植物抗低溫脅迫機(jī)制研究奠定了基礎(chǔ)。本研究利用高通量測序平臺Illumina HiSeq對低溫脅迫條件下2種萱草的根莖進(jìn)行轉(zhuǎn)錄組測序分析，期望揭示萱草根莖轉(zhuǎn)錄組在低溫脅迫條件下發(fā)生的表達(dá)差異性特征，為萱草冷應(yīng)激功能基因的鑒定、次生代謝途徑的解析及其低溫脅迫應(yīng)激調(diào)控機(jī)制的研究提供參考。

1 材料與方法

1.1 材料

萱草Jinyan 7(耐冷型品種，冬季地上部常綠)和Lucretius 4(冷敏感型品種，冬季地上部枯萎)由蘇州農(nóng)業(yè)職業(yè)技術(shù)學(xué)院相城科技園提供。在通風(fēng)溫室中(20 ℃，光周期：16 h光照8 h黑暗)，植物于塑料盆中生長至6片葉子，然后移入生長室進(jìn)行處理。將溫度從20 ℃逐漸降低至-4 ℃開始低溫處理，具體方法：讓植株20 ℃生長48 h，然后分別在10 ℃和4 ℃誘導(dǎo)48 h，最后將它們在-4 ℃處理8 h。收集在20 ℃(作為對照組，CK)和-4 ℃(作為處理組，T)生長的萱草根莖，在液氮中冷凍用于RNA提取和生理分析。收集的4個(gè)樣品分別命名為4CK(20 ℃時(shí)Lucretius 4的根莖)、4T(-4 ℃處理的Lucretius 4根莖)、7CK(20 ℃時(shí)Jinyan 7根莖)、7T(-4 ℃處理的Jinyan 7根莖)。每3株植物根莖作為1個(gè)生物學(xué)重復(fù)，共設(shè)置2組重復(fù)。

1.2 RNA提取與文庫構(gòu)建

采用TRNzol Universal Reagent植物總RNA提取試劑盒，分別提取冷處理前后耐冷型和冷敏感型萱草根莖總RNA，Agilent 2100 Bioanalyzer檢測總RNA的質(zhì)量，NanoDrop TM2000檢測總RNA的濃度。利用Oligo(dT)磁珠富集每個(gè)萱草文庫的mRNA，加入適量打斷劑高溫下片斷化mRNA，并以片段化的mRNA為模板合成cDNA，選擇合適的片段進(jìn)行PCR擴(kuò)增，構(gòu)建文庫。

1.3 轉(zhuǎn)錄組測序與組裝

將構(gòu)建的萱草根莖轉(zhuǎn)錄組文庫提交Illumina HiSeq平臺進(jìn)行高通量測序。首先對Illumina測序的原始reads進(jìn)行處理，去除接頭污染、未知堿基N含量偏高(>5%)和低質(zhì)量reads，然后將獲得的高質(zhì)量clean reads利用Trinity[22]進(jìn)行denovo組裝分析，最后通過Tgicl對其進(jìn)行聚類和去冗余處理，得到unigenes。

1.4 轉(zhuǎn)錄組功能注釋

利用Blastn對得到的unigene序列進(jìn)行Nt數(shù)據(jù)庫注釋，通過Blastx將unigenes序列與以下蛋白數(shù)據(jù)庫進(jìn)行比對(E-value<1×10-5，E-value為比對的期望值)：Nr(NCBI nonredundant protein)、SwissPort、KOG(clusters of orthologous groups)和KEGG(kyoto encyclopedia of genes and genomes)，得到和相應(yīng)unigenes序列相似性最高的蛋白序列，從而得到unigenes的對應(yīng)注釋信息。使用Blast2GO軟件，同時(shí)結(jié)合Nr注釋得到unigenes的GO(gene ontology)功能注釋。

1.5 蛋白編碼框和轉(zhuǎn)錄因子預(yù)測

通過軟件Transdecoder(https://transdecoder.github.io)來識別unigenes的候選編碼區(qū)域，具體過程主要包括通過TransDecoder.LongOrfs軟件從unigenes中將最長的開放閱讀框(open reading frame，ORF)提取出來，接著將其與SwissProt數(shù)據(jù)庫進(jìn)行Blast比對，其結(jié)果經(jīng)Hmmscan搜索得到Pfam蛋白的同源序列，從而實(shí)現(xiàn)蛋白編碼框(coding sequence，CDS)的預(yù)測。通過getorf軟件來檢測unigenes的ORF，接著利用Hmmsearch檢索將ORF與PlantTFDB數(shù)據(jù)庫中轉(zhuǎn)錄因子蛋白結(jié)構(gòu)域進(jìn)行比對，最后根據(jù)PlantTFDB中關(guān)于轉(zhuǎn)錄因子家族特征的描述對unigenes進(jìn)行鑒定。

1.6 Unigenes基因表達(dá)水平分析

通過MISA(http://pgrc.ipk-gatersleben.de/misa)軟件來檢測萱草轉(zhuǎn)錄組的unigenes，搜索和驗(yàn)證簡單重復(fù)序列(simple sequence repeats， SSR)，并進(jìn)行相關(guān)統(tǒng)計(jì)分析。利用軟件Bowtie2將clean reads與unigenes進(jìn)行比對，然后通過RSEM軟件計(jì)算得到各個(gè)樣品的基因表達(dá)水平，并參照Saddhe等[23]的差異基因檢測方法，通過軟件PossionDis [fold change≥2.00和錯(cuò)誤發(fā)現(xiàn)率(FDR)≤0.001]得到不同樣本之間的差異基因。

1.7 差異表達(dá)基因聚類、GO功能、Pathway功能分析

差異基因的層次聚類分析采用R軟件中的pheatmap函數(shù)分析法。依據(jù)GO功能的官方分類對差異基因進(jìn)行功能注釋，根據(jù)KEGG的官方分類對差異基因的生物通路進(jìn)行注釋分析，差異基因富集分析則利用R軟件中的phyper函數(shù)分析法。

1.8 差異表達(dá)基因的qRT-PCR驗(yàn)證

通過對轉(zhuǎn)錄組中差異表達(dá)基因的進(jìn)一步篩選，隨機(jī)選擇8個(gè)顯著差異表達(dá)的基因，進(jìn)行qRT-PCR驗(yàn)證。qRT-PCR檢測使用Takara公司的SYBR?Premix Ex TaqTM試劑，在熒光定量儀中進(jìn)行。內(nèi)參基因?yàn)?8S rRNA，正反向引物序列分別為：5′-CCATAAACGATGCCGACC-3′；5′-TTGCGACCATACTCCCCC-3′。PCR反應(yīng)條件為：95 ℃ 5 min；95 ℃ 5 s，60 ℃ 30 s，72 ℃ 10 s，40個(gè)循環(huán)。qRT-PCR中使用的引物序列如表1所示。對每個(gè)基因和樣本進(jìn)行3個(gè)生物重復(fù)和3個(gè)技術(shù)重復(fù)。用2-ΔΔCt法計(jì)算差異基因相對表達(dá)水平。

2 結(jié)果與分析

2.1 萱草轉(zhuǎn)錄組組裝與質(zhì)量分析

從表2可以看出，4個(gè)樣品4CK、4T、7CK、7T分別得到47 738 570、50 110 130、50 918 552、38 540 808條raw reads，過濾后分別產(chǎn)生46 171 660、47 791 962、47 411 290、35 880 734條高質(zhì)量clean reads，各包含6 858 685 119、7 099 131 160、7 047 333 003、5 342 671 114個(gè)核苷酸信息，Q20(質(zhì)量值≥20的堿基所占百分比)分別為98.61%、98.63%、98.72%、98.82%，Q30(質(zhì)量值≥30的堿基所占百分比)分別為95.65%、95.7%、95.97%、96.23%，GC含量分別為48.83%、48.89%、48.57%、48.48%，表明clean reads的質(zhì)量符合實(shí)驗(yàn)要求，可用于后續(xù)分析。樣品4CK、4T、7CK、7T經(jīng)Trinity組裝后各獲得55 348、51 400、5 6962、49 915個(gè)unigene，平均長度分別為796、824、763、753 nt，N50(按轉(zhuǎn)錄本長度從大到小排序后逐個(gè)累加至所有轉(zhuǎn)錄本總長度的50%)分別為1 294、1 337、1 225、1 197 nt。

表2 有效數(shù)據(jù)評估統(tǒng)計(jì)Table 2 Valid data evaluation statistics

Unigenes的長度分布圖(圖1)顯示，Unigenes長度在300～>3 000 nt；長度為300 nt的unigenes數(shù)量最多，共25 011條；400 nt的unigenes數(shù)量次之，共11 341條；其他長度的unigenes均低于10 000條，大于3 000 nt的unigenes共2 194條。

2.2 萱草轉(zhuǎn)錄組unigenes的功能注釋

為獲得組裝好的unigenes的表達(dá)和功能注釋，將unigene序列與Nr、Nt、Swissprot、KEGG、KOG等數(shù)據(jù)庫進(jìn)行比對。結(jié)果表明，通過Nt數(shù)據(jù)庫可對10 453條unigenes(12.00%)進(jìn)行成功比對并進(jìn)行注釋，Nr數(shù)據(jù)庫可對38 447條unigenes(44.15%)進(jìn)行成功比對并進(jìn)行注釋，在SwissPort、KOG、KEGG、Pfam、GO數(shù)據(jù)庫中分別有28 095(32.27%)、41 418(47.57%)、27 835(31.97%)、36 767(42.23%)、15 933(18.30%)條unigenes獲得成功比對和注釋。

對Nr數(shù)據(jù)庫的比對注釋結(jié)果進(jìn)行分析，圖2是unigenes注釋同源基因的物種分布圖，可以發(fā)現(xiàn)，在相似序列較高的物種中，蘆筍(Asparagusofficinalis)所占比例最高，其數(shù)量可以達(dá)到24 115條(45.71%)；其次為油棕(Elaeisguineensis)、海棗(Phoenixdactylifera)和鳳梨(Ananascomosus)，其數(shù)量分別為5 695條(10.79%)、3 673條(6.96%)和1 638條(3.1%)；注釋到其他物種的unigenes共17 633條，占比為33.43%。

2.3 CDS和轉(zhuǎn)錄因子分析

通過BLAST比對分析，從4組萱草樣品轉(zhuǎn)錄組的全部unigenes中共獲得42 133個(gè)CDS 序列，CDS的長度全部大于200 nt，其中有16 337條CDS長度超過1 000 nt，3 923條CDS長度超過2 000 nt(圖3)。通過轉(zhuǎn)錄因子預(yù)測分析可以發(fā)現(xiàn)，這些unigenes分屬于56個(gè)家族(圖4)，而其中Homeobox、zf-C2H2、bHLH、HMG、NGFIB-like、MYB、zf-GATA、ZBTB、TF_bZIP與CG-1類占主體，這說明眾多轉(zhuǎn)錄因子參與萱草根莖的生理代謝。

圖3 CDS長度分布圖Fig.3 Distribution diagram of coding sequence length

圖4 轉(zhuǎn)錄因子家族分類圖Fig.4 Classification of transcription factor families

2.4 差異表達(dá)基因統(tǒng)計(jì)、GO和KEGG分析

表3 Unigenes注釋結(jié)果統(tǒng)計(jì)Table 3 Annotation results of unigene sequences

2.4.1 差異表達(dá)基因統(tǒng)計(jì)

對冷敏感型和耐冷型萱草冷處理后的差異基因數(shù)量進(jìn)行統(tǒng)計(jì)發(fā)現(xiàn)，與對照相比，冷敏感型萱草根莖在低溫脅迫下，總共檢測到6 153個(gè)基因表達(dá)變化，其中，上調(diào)表達(dá)基因有1 849個(gè)，下調(diào)表達(dá)基因共檢測到4 304個(gè)。而耐冷萱草根莖在低溫脅迫下，共有14 813個(gè)差異表達(dá)基因，其中，有6786個(gè)基因上調(diào)表達(dá)，8 027個(gè)基因表達(dá)量下調(diào)。對照組中，與冷敏感型萱草相比，耐冷萱草有10 282個(gè)基因上調(diào)表達(dá)，12 234個(gè)基因下調(diào)表達(dá)；而在低溫處理?xiàng)l件下，與冷敏感型萱草相比，耐冷萱草上調(diào)表達(dá)基因有9 224個(gè)，下調(diào)表達(dá)基因有10 614個(gè)(圖5)。

圖5 上調(diào)與下調(diào)差異基因數(shù)量統(tǒng)計(jì)Fig.5 Numbers of up-regulated and down-regulated differentially expressed genes

2.4.2 差異表達(dá)基因GO富集分析

為了進(jìn)一步了解低溫脅迫條件下萱草的抗逆機(jī)制，對差異表達(dá)基因進(jìn)行GO富集分析。GO分析結(jié)果表明(圖6)，這些差異表達(dá)基因分為3大類：分子功能(molecular function)、細(xì)胞組分(cellular component)和生物學(xué)過程(biological process)，每個(gè)大類對應(yīng)的功能組數(shù)量分別為16、19、29。差異基因大部分富集在生物學(xué)過程中，而參與分子功能的差異基因相對較少。在生物學(xué)過程中，富集較明顯的包括細(xì)胞過程(7 216)、代謝過程(6 652)、生物調(diào)節(jié)(2 138)、生物過程調(diào)節(jié)(1 973)和定位(1 700)；參與細(xì)胞組分的差異基因主要富集在細(xì)胞(6 705)、細(xì)胞組分(6 614)、膜(5 386)、細(xì)胞器(4 947)和膜部分(4 934)；而參與分子功能的差異基因富集類別主要有催化活性(8 314)、結(jié)合(7 916)和轉(zhuǎn)運(yùn)活性(1 916)。

圖6 差異基因GO功能分類圖Fig.6 Gene ontology classification of differentially expressed genes

2.4.3 差異表達(dá)基因KEGG分析

根據(jù)KEGG代謝通路對差異基因進(jìn)行功能注釋，發(fā)現(xiàn)參與的代謝通路包含6個(gè)分支：細(xì)胞過程(cellular processes)、環(huán)境信息處理(environmental information processing)、遺傳信息處理(genetic information processing)、人類疾病(human diseases)、代謝(metabolism)和有機(jī)系統(tǒng)(organismal systems)，具體的差異基因數(shù)量見圖7。在代謝途徑中，占據(jù)前5的類別分別為全局概述圖(5 527)、碳水化合物代謝(1 694)、核苷酸代謝(1 673)、氨基酸代謝(1 311)和脂質(zhì)代謝(1 183)；在遺傳與信息處理類別中，翻譯(2 612)、轉(zhuǎn)錄(2 165)和折疊分揀與降解(1 781)差異基因顯著富集，而信號轉(zhuǎn)導(dǎo)(5 416)在環(huán)境信息處理通路中差異基因富集較多。

圖7 差異表達(dá)基因KEGG分類Fig.7 KEGG classification of differentially expressed genes

2.4.4 萱草低溫脅迫相關(guān)代謝通路分析

對差異基因的代謝通路進(jìn)行分析發(fā)現(xiàn)，差異基因富集程度比較明顯的途徑主要包括：異黃酮生物合成(isoflavonoid biosynthesis)、牛磺酸和亞牛磺酸代謝(taurine and hypotaurine metabolism)、光合天線作用蛋白(photosynthesis-antenna proteins)、二苯基庚酮和姜酚生物合成(stilbenoid， diarylheptanoid and gingerol biosynthesis)、植物晝夜節(jié)律(circadian rhythm-plant)、不飽和脂肪酸的生物合成(biosynthesis of unsaturated fatty acids)等，具體的代謝通路與差異基因數(shù)量見表4。

表4 與萱草低溫脅迫作用相關(guān)的16條代謝通路Table 4 Sixteen metabolic pathways related to cold stress of Hemerocallis fulva

2.5 差異基因參與對低溫的反應(yīng)

通過對差異基因的進(jìn)一步篩選發(fā)現(xiàn)，Ca2+信號通路相關(guān)基因表達(dá)發(fā)生明顯差異(表5)。與對照組相比，低溫處理后，冷敏感型萱草中CBL互作絲氨酸/蘇氨酸蛋白激酶3基因顯著上調(diào)；

表5 鈣離子信號通路相關(guān)差異基因Table 5 Differentially expressed genes involved in Ca2+ signalling pathway in each comparison

在耐冷型萱草中多種Ca2+信號通路相關(guān)基因明顯差異表達(dá)，鈣依賴性蛋白激酶基因CAMK2A、鈣調(diào)蛋白基因(Calm5、Calm3)表達(dá)顯著上調(diào)，而鈣調(diào)蛋白和肌聯(lián)蛋白結(jié)合RhoGEF基因、CBL互作絲氨酸/蘇氨酸蛋白激酶14基因則顯著下調(diào)。對照組中，與耐冷型萱草相比，冷敏感型萱草鈣調(diào)蛋白基因(Calm5、Calm3)、CBL互作絲氨酸/蘇氨酸蛋白激酶14基因顯著上調(diào)，鈣依賴性蛋白激酶基因CAMK2A顯著下調(diào)。在低溫脅迫下，與耐冷型萱草相比，冷敏感型萱草鈣調(diào)蛋白基因(Calm5、Calm3)、鈣依賴性蛋白激酶基因CAMK2A與對照組表達(dá)趨勢相同。CL7725.Contig5_All基因

編碼了鈣調(diào)蛋白和肌聯(lián)蛋白結(jié)合RhoGEF(calmodulin and titin-interacting RhoGEF)。低溫脅迫下，與對照組相比，耐冷型萱草中CL7725.Contig5_All基因顯著上調(diào)；與耐冷型萱草相比，冷敏感型萱草的CL7725.Contig5_All基因顯著下調(diào)；與對照組相比，冷敏感型萱草的CL7725.Contig5_All基因表達(dá)無顯著差異。對照組中，耐冷型萱草與冷敏感型萱草的CL7725.Contig5_All基因表達(dá)也沒有顯著差異。Unigene6945_All基因編碼鈣調(diào)素結(jié)合轉(zhuǎn)錄激活子3(calmodulin-binding transcription activator 3)，低溫脅迫抑制了其在耐冷萱草中的表達(dá)，在對照組中，該基因在2種萱草中的表達(dá)沒有差異。

此外，也篩選到MAPK信號通路相關(guān)差異基因(表6)，低溫處理前后比對結(jié)果發(fā)現(xiàn)，在冷敏感型萱草(4CK-VS-4T)中，促分裂原活化的蛋白激酶基因MAPKKKK3顯著上調(diào)，MAPK1-like、MAPKKKK5顯著下調(diào)。在耐冷型萱草中(7CK-VS-7T)，促分裂原活化的蛋白激酶基因MLK4-like、MAPKKKK3、MAPK1、MAPK4顯著上調(diào)，MAPK1-like、MAPK7-like、MAPKKKK5顯著下調(diào)。對照組中，與冷敏感型萱草相比，MAPK8互作蛋白2基因(MAPK8IP2)、MAPKKKK3基因、MAPK7-like基因上調(diào)表達(dá)，MAPK4、MAPKKKK5基因下調(diào)表達(dá)。低溫處理后，耐冷型萱草與冷敏感型萱草的MAPKKKK3、MAPK8IP2基因表達(dá)上調(diào)，MAPK7-like、MLK4-like、MAPKKKK5表達(dá)下調(diào)。進(jìn)一步分析發(fā)現(xiàn)，低溫可以抑制MAPK1基因(CL4304.Contig1_All)的表達(dá)，常溫下，2種萱草中沒有檢測到該基因的表達(dá)差異，而在低溫下，冷敏感萱草中，其下調(diào)倍數(shù)為-1.36，在耐冷型萱草中，其表達(dá)被抑制的更加明顯，下調(diào)倍數(shù)為-3.90。

表6 MAPK信號通路相關(guān)差異基因Table 6 Differentially expressed genes involved in MAPK cascades pathway in each comparison

對編碼熱激蛋白(表7)相關(guān)的基因進(jìn)行篩選發(fā)現(xiàn)，冷敏感型萱草(4CK-VS-4T)中，熱激蛋白基因HSP70-2、DNAJC6、HSP90顯著上調(diào)，而在耐冷型萱草中(7CK-VS-7T)，僅熱激蛋白基因HSP86上調(diào)表達(dá)，其他熱激蛋白基因HSP90-5、HSP90均顯著下調(diào)。不同的溫度條件下，對2種萱草進(jìn)行對比發(fā)現(xiàn)，對照條件下，與冷敏感型萱草相比，Hspa8基因上調(diào)表達(dá)，HSP70-2、HSP90-5基因下調(diào)表達(dá)，低溫處理后，2種萱草的對比中篩選到的熱激蛋白基因均下調(diào)表達(dá)。HSP90基因(Unigene6560_All)在2種萱草中表現(xiàn)出相反的表達(dá)趨勢，冷敏感萱草中，低溫處理后上調(diào)，而在耐冷型萱草中，與對照條件相比下調(diào)-1.61，與冷敏感萱草相比下調(diào)-2.61。

表7 編碼熱激蛋白的差異表達(dá)基因Table 7 Differentially expressed genes encoding HSPs protein in each comparison

對差異表達(dá)的轉(zhuǎn)錄因子(表8)進(jìn)行分析發(fā)現(xiàn)，在冷敏感型萱草(4CK-VS-4T)中，轉(zhuǎn)錄因子WRKY11、MYB102顯著下調(diào)表達(dá)，而在耐冷型萱草(7CK-VS-7T)中，轉(zhuǎn)錄因子MYB20、bHLH2上調(diào)表達(dá)，MYB102下調(diào)表達(dá)。與冷敏感萱草相比，對照組耐冷型萱草WRKY11上調(diào)表達(dá)，bHLH2、4CK-VS-4T，Lucretius 4低溫與對照相比；7CK-VS-7T，Jinyan 7低溫和對照相比；4CK-VS-7CK，對照條件下Jinyan 7與Lucretius 4相比；4T-VS-7T，低溫處理下Jinyan 7與Lucretius 4相比。下同。

表8 差異表達(dá)的轉(zhuǎn)錄因子Table 8 Differential expression of transcription factors

2.6 熒光定量PCR驗(yàn)證

為了驗(yàn)證轉(zhuǎn)錄組測序結(jié)果，隨機(jī)選擇8個(gè)表達(dá)水平差異顯著的基因進(jìn)行qRT-PCR試驗(yàn)，分別為：類硫氧還蛋白基因(thioredoxin-like)、絲氨酸/精氨酸重復(fù)基質(zhì)蛋白基因(serine/arginine repetitive matrix protein 3-like)、溶磷酰膽堿酰轉(zhuǎn)移酶1基因(Lysophos phatidylcholine acyltransferase 1)、1-磷脂酰肌醇 3-磷酸 5-激酶基因(1-phosphatidylinositol 3-phosphate 5-kinase)、粘液素-2基因(Mucin-2)、肌球蛋白VIIa基因(myosin-VIIa-like)、

含有LOB結(jié)構(gòu)域的蛋白38基因(LOB domain-containing protein 38-like)、多聚泛素-A基因(polyubiquitin-a)。結(jié)果顯示：低溫脅迫下，8個(gè)基因的qRT-PCR結(jié)果與轉(zhuǎn)錄組分析的趨勢一致(圖8)，只是具體表達(dá)倍數(shù)上存在差異，這些差異可能是由于轉(zhuǎn)錄組測序和qRT-PCR實(shí)驗(yàn)的檢測靈敏度，以及實(shí)驗(yàn)結(jié)果分析方法不同造成的，表明轉(zhuǎn)錄組測序結(jié)果可靠。

3 結(jié)論與討論

低溫是許多植物所面臨的常見威脅，其不僅可誘導(dǎo)細(xì)胞脫水和膜損傷[23-24]，還會(huì)改變植物細(xì)胞的正常生理代謝過程，主要表現(xiàn)在2個(gè)方面：呼吸作用障礙[25]、活性氧的積累和去除[26-28]。萱草是一種在許多地區(qū)廣泛種植的營養(yǎng)植物，耐寒品種的選育對西北地區(qū)的營養(yǎng)植物生長具有重要意義。低溫脅迫通過抑制生長和發(fā)育嚴(yán)重限制了植物物種的地理分布，植物可以通過調(diào)節(jié)基因表達(dá)來適應(yīng)低溫脅迫?？沟蜏孛{迫的植物可以通過防御機(jī)制調(diào)節(jié)自身，包括生理反應(yīng)、代謝和信號轉(zhuǎn)導(dǎo)途徑的適應(yīng)，以及基因表達(dá)變化[29]。本研究對2種不同抗性萱草在低溫脅迫條件下的轉(zhuǎn)錄組進(jìn)行分析發(fā)現(xiàn)，在7個(gè)數(shù)據(jù)庫同時(shí)注釋的unigenes共有6 727條，至少有1種數(shù)據(jù)庫注釋成功的unigenes共有56 679條，其

qRT-PCR結(jié)果中的誤差線表示平均值的標(biāo)準(zhǔn)差(n=3)。The error line in results represents the standard deviation of the mean value (n=3).圖8 差異表達(dá)基因的qRT-PCR驗(yàn)證Fig.8 qRT-PCR validation for differentially expressed genes identified by RNA-Seq

比例達(dá)到65.09%；還有30 394條unigenes未獲得注釋，這一比例和多種植物的測定結(jié)果相類似[30-32]，說明萱草轉(zhuǎn)錄組中存在著大量特征序列，以及大量unigenes功能有待開發(fā)。與對照組相比，冷敏感型萱草根莖在低溫脅迫下，上調(diào)表達(dá)基因有1 849條，下調(diào)表達(dá)基因有4 304條；而耐冷萱草根莖在低溫脅迫下，有6 786條基因上調(diào)表達(dá)，8 027條基因下調(diào)表達(dá)，差異表達(dá)基因數(shù)量顯著多于冷敏感型萱草品種，猜測其抗凍性可能與迅速調(diào)控多種基因差異表達(dá)有關(guān)。差異表達(dá)基因的GO分類結(jié)果顯示，2種萱草根莖的DEG功能特性與生物過程、分子功能和細(xì)胞組分相關(guān)。2種萱草低溫脅迫和非低溫脅迫條件下差異表達(dá)基因均涉及KEGG代謝通路的五大分支。

次生代謝物在植物逆境脅迫下的合成和積累具有重要意義。已有研究發(fā)現(xiàn)，低溫等逆境脅迫能夠誘導(dǎo)類黃酮化合物的積累，提高植物抗逆性[33]；一些植物在低溫脅迫下體內(nèi)類黃酮的生物合成會(huì)發(fā)生顯著變化[34]。本研究也發(fā)現(xiàn)，低溫脅迫下萱草根莖的類黃酮合成基因表達(dá)水平顯著高于對照組，這說明低溫脅迫在一定程度上可能誘導(dǎo)了類黃酮的積累?；ㄇ嗨?、吲哚生物堿等物質(zhì)的合成也會(huì)受到低溫脅迫的影響。

低溫刺激植物體內(nèi)Ca2+信號通路的變化，Ca2+信號通路被三大類蛋白質(zhì)感知：CDPKs、CaMs和CBLs[35]。在水稻中，低溫脅迫下，鈣調(diào)蛋白OsMSR2基因在不同發(fā)育階段不同組織均顯著上調(diào)表達(dá)；擬南芥中外源表達(dá)該基因可以影響其耐鹽性和耐旱性[36]。在耐冷型萱草中，2個(gè)鈣調(diào)蛋白基因Calm5/Calm3均顯著上調(diào)表達(dá)，而在冷敏感型萱草中，沒有發(fā)現(xiàn)相應(yīng)基因的變化。在低溫處理的茶樹中，5個(gè)鈣調(diào)蛋白基因、2個(gè)CDPK基因和1個(gè)CBL基因參與了信號轉(zhuǎn)導(dǎo)[37]。本研究中，CDPK3基因表達(dá)在冷敏感型萱草中顯著上調(diào)，耐冷型萱草中CDPK14基因表達(dá)顯著下調(diào)。Ca2+信號通路相關(guān)基因表達(dá)在2種萱草中差異較顯著，其中calmodulin and titin-interacting RhoGEF(CL7725.Contig5_All)與calmodulin-binding transcription activator 3(Unigene6945_All)基因均在低溫脅迫的耐冷萱草中被顯著抑制，而在冷敏感萱草中無明顯差異，猜測這2個(gè)基因在Ca2+信號通路對萱草耐冷調(diào)控中發(fā)揮重要作用。MAPKs信號通路也影響植物非生物脅迫響應(yīng)，許多受體蛋白激酶(RLKs)負(fù)責(zé)感知外部環(huán)境的變化和相應(yīng)的信號轉(zhuǎn)導(dǎo)[38]。在冷敏感型萱草中，MAPKKKK3上調(diào)表達(dá)，MAPKKKK5下調(diào)表達(dá)；在耐冷型萱草中也發(fā)現(xiàn)這2個(gè)基因有相同的表達(dá)趨勢，且上下調(diào)程度更顯著。冷敏感型萱草中僅發(fā)現(xiàn)1個(gè)MAPK1基因下調(diào)表達(dá)，而耐冷型萱草中共檢測到3個(gè)差異表達(dá)的MAPK基因(MAPK1、MAPK4、MAPK7-like)和1個(gè)差異表達(dá)的MLK基因。這些發(fā)現(xiàn)為Ca2+在萱草低溫響應(yīng)中的關(guān)鍵作用提供了證據(jù)，同時(shí)進(jìn)一步證明MAPK通路影響萱草對低溫脅迫的響應(yīng)。MAPKs信號通路基因MAPK1(CL4304.Contig1_All)的表達(dá)受低溫抑制，在耐冷型萱草中的抑制程度明顯高于冷敏感萱草，表明MAPK1基因可能作為MAPKs信號通路的關(guān)鍵基因發(fā)揮作用。

轉(zhuǎn)錄因子在植物非生物脅迫應(yīng)答中發(fā)揮重要作用，在擬南芥中過量表達(dá)水稻MYB4(R2R3-型MYB轉(zhuǎn)錄因子)和OsMYB3R-2(R1R2R3型MYB轉(zhuǎn)錄因子)基因，可顯著提高擬南芥抗凍性[39]。冷敏感型萱草中，MYB家族轉(zhuǎn)錄因子MYB102顯著下調(diào)表達(dá)，該基因在耐冷型萱草中下調(diào)倍數(shù)較低。此外，另一個(gè)MYB家族轉(zhuǎn)錄因子MYB20僅在耐冷型萱草中上調(diào)表達(dá)，在冷敏感型萱草未發(fā)現(xiàn)該基因表達(dá)變化。研究表明，WRKY和bHLH轉(zhuǎn)錄因子家族與植物低溫脅迫密切相關(guān)[40-41]。bHLH轉(zhuǎn)錄因子家族基因bHLH2僅在耐冷型萱草低溫處理后上調(diào)表達(dá)，而WRKY轉(zhuǎn)錄因子家族基因WRKY11在冷敏感型萱草中下調(diào)表達(dá)；與冷敏感型萱草相比，耐冷萱草中的WRKY11在常溫和低溫條件表達(dá)水平均較高，WRKY11在2種萱草中的表達(dá)差異及其在應(yīng)對低溫脅迫中的作用機(jī)理需要進(jìn)一步探討。以上結(jié)果表明，MYB、WRKY、bHLH轉(zhuǎn)錄因子家族在冷敏感型萱草和耐冷型萱草低溫脅迫響應(yīng)中發(fā)揮了重要而不同的作用。除了蛋白激酶和轉(zhuǎn)錄因子，一些功能蛋白也參與了低溫脅迫過程。熱休克蛋白(HSPs)最初是通過參與對高溫的反應(yīng)而被識別的，越來越多的研究發(fā)現(xiàn)，許多熱休克蛋白也對低溫有反應(yīng)。在冷敏感型萱草和耐冷型萱草中均發(fā)現(xiàn)一些HSP蛋白基因呈現(xiàn)差異表達(dá)，但在2種萱草中篩選到的HSP家族成員不同。有趣的是，在2種萱草中均發(fā)生變化的HSP90，在冷敏感型萱草中上調(diào)表達(dá)，而在耐冷型萱草中則下調(diào)表達(dá)，推測HSP90可能是萱草響應(yīng)低溫脅迫的重要候選基因，其在萱草低溫脅迫中的作用還需要進(jìn)一步研究。

綜上所述，本研究利用高通量測序鑒定和比較了萱草Jinyan 7(耐冷型品種)和Lucretius 4(冷敏感型品種)的低溫響應(yīng)基因，發(fā)現(xiàn)在低溫脅迫下，Ca2+信號通路基因、MAPKs信號通路基因、轉(zhuǎn)錄因子和熱激蛋白基因在2種萱草中表現(xiàn)出不同的變化趨勢，這些關(guān)鍵路徑基因表達(dá)的差異可能導(dǎo)致了2種萱草在低溫脅迫中的不同反應(yīng)。本研究為萱草低溫響應(yīng)的分子機(jī)制提供了思路，同時(shí)，篩選出的差異基因也可能成為培育耐冷型萱草品種的候選基因。