田雪慧 ，郁繼華，頡建明

（1.楊凌職業(yè)技術(shù)學(xué)院生態(tài)環(huán)境工程分院，陜西 楊凌 712100；2.甘肅農(nóng)業(yè)大學(xué)園藝學(xué)院，甘肅 蘭州 730070）

【研究意義】百合（Lilium）是國家衛(wèi)生部首批公布的藥食同源植物，也是著名的觀賞植物［1］。食用百合是營養(yǎng)價值很高的特種蔬菜，其鱗莖含碳水化合物、蛋白質(zhì)、果膠、脂肪、鉀、磷以及鈣、鐵、百合甙等抗癌性植物堿和多種維生素，兼具保健功能，產(chǎn)品暢銷中國及東南亞地區(qū)，在世界范圍享有盛譽，市場前景十分廣闊。蘭州百合〔Lilium davidiiDuch.var.unicolor(Hong) Cotton〕微甜、口感好，是食用百合中的精品，生長在 103°54′ 54′′ E、36°06′ 12′′ N、海拔2 000多米的七里河區(qū)，是食用百合中的耐寒品種代表［2］。溫度是植物生長最重要的環(huán)境因子之一，通過對蘭州百合抗凍轉(zhuǎn)錄組的研究，篩選出百合抗凍的基因及抗凍相關(guān)通路，為后續(xù)蘭州百合SSR標(biāo)記開發(fā)、抗寒性相關(guān)基因的發(fā)掘及分子育種提供借鑒。

【前人研究進展】目前，對百合抗寒性研究集中于亞洲百合［3］、東方百合［4］、卷丹百合［5］等觀賞百合品種，對百合耐寒的研究主要集中在細胞膜透性和保護酶系統(tǒng)方面，相對電導(dǎo)率、脂膜氧化產(chǎn)物丙二醛、抗氧化酶類（SOD、CAT、POD）活性通常被作為植物耐寒性生理指標(biāo)［6-7］。HOSHI等［8］對百合在農(nóng)桿菌介導(dǎo)下轉(zhuǎn)錄組進行研究，建立了東方百合雜交品種的轉(zhuǎn)基因生產(chǎn)體系；杜方等［9］對百合不同器官的轉(zhuǎn)錄數(shù)據(jù)進行測序并開發(fā)了SSR標(biāo)記。

【本研究切入點】關(guān)于抗寒食用百合代表——蘭州百合的抗寒轉(zhuǎn)錄因子研究至今仍是空白。本研究運用先進的Illumina高通量測序平臺的轉(zhuǎn)錄組測序技術(shù)，分別對正常溫度條件下和零下超低溫條件下生長的蘭州百合苗期葉片的轉(zhuǎn)錄產(chǎn)物mRNA進行測序，比較常溫與零低溫下差異表達的基因及通路，了解蘭州百合超低溫脅迫下的基因表達狀況?！緮M解決的關(guān)鍵問題】進一步對差異表達基因進行GO、KEGG(Kyoto Encyclopedia of Genes and Genomes)注釋和富集分析，以期進一步篩選出與蘭州百合抗凍相關(guān)的基因［10］。

1 材料與方法

1.1 試驗材料

在蘭州市七里河區(qū)采集露地自然生長成熟的蘭州百合種球，貯藏于2℃冰箱110 d［11］后，種植于15 cm口徑的花盆中，置于20℃智能低溫人工氣候箱（RXZ-0288型）中培養(yǎng)，相同的水肥條件管理。將生長至15 cm左右的百合幼苗放在智能低溫人工氣候箱中進行-4℃低溫處理［12］，降至目標(biāo)溫度后維持24 h［13-14］，以室溫20℃作為對照，每個處理3次重復(fù)。

1.2 試驗方法

1.2.1 轉(zhuǎn)錄測序 分別取3份20℃與-4℃的蘭州百合葉片樣品，采用TRIzol法提取總RNA［15］，并對RNA濃度、純度、完整性進行檢測。樣品檢測合格后，用帶有Oligo (dT)的磁珠富集真核生物mRNA［16］。隨后加入fragmentation buffer將mRNA打斷成短片段，以打斷后的 mRNA為模板合成一鏈cDNA，再合成二鏈cDNA。然后純化、修復(fù)，再用AMPure XP beads進行片段大小選擇，最后進行PCR擴增，純化PCR產(chǎn)物，得到最終文庫。質(zhì)檢合格后使用Illumina 4000高通量測序平臺(HiSeq/MiSeq)進行測序［16］。由Illumina 4000高通量測序平臺測序所得的數(shù)據(jù)稱為 raw reads或raw data。然后進行轉(zhuǎn)錄結(jié)果組裝與拼接，并進行基因功能注釋與表達水平分析。

1.2.2 實時熒光定量PCR驗證 隨機選取10個差異性表達的基因進行qRT-PCR分析驗證。以20 ℃與-4 ℃的蘭州百合葉片RNA為材料，采用HiScriptⅡQRT SuperMix for qPCR試劑盒反轉(zhuǎn)錄合成cDNA，設(shè)計引物（表1）進行qRT-PCR，以Actin基因作為內(nèi)參［17］。各處理均設(shè)3次重復(fù)，計算基因的相對表達量。

表1 定量PCR引物序列Table 1 Primer sequences for qRT-PCR

2 結(jié)果與分析

2.1 測序結(jié)果與Trinity拼接

蘭州百合20℃轉(zhuǎn)錄組共得到56 882 837 raw reads，56 456 634 clean reads，Q20值為94.56%，Q30值為87.32%，序列堿基GC含量為44.84%；-4℃轉(zhuǎn)錄組共得到58 092 923 raw reads，57 800 442 clean reads，Q20值為94.51%，Q30值為87.32%，序列堿基GC含量為44.54%。表明蘭州百合轉(zhuǎn)錄組測序質(zhì)量較高，可為后續(xù)的數(shù)據(jù)組裝與分析提供良好基礎(chǔ)。經(jīng)過Trinity 轉(zhuǎn)錄本拼接和CD-HIT去冗余序列處理結(jié)果見表2。Unigene序列長度＜300的最多（圖1），說明轉(zhuǎn)錄拼接較好，利于后續(xù)分析研究。

表2 處理前后序列數(shù)量、N50、N75、N90長度和總堿基數(shù)統(tǒng)計Table 2 Sequence number, length of N50, N75, N90 and total base number by different data processing

圖1 拼接后Unigene 長度分布Fig.1 Length distribution of all unigenes after splicing

2.2 基因功能注釋結(jié)果

對獲得的轉(zhuǎn)錄組序列進行數(shù)據(jù)庫比對注釋可知，注釋的Unigene數(shù)量共238 301個，Nr（NCBI 蛋白質(zhì)非冗余數(shù)據(jù)庫）注釋的Unigene數(shù)量為65 314個、占27.41%，KOG（真核生物直系同源基因數(shù)據(jù)庫）注釋的Unigene數(shù)量為42 138個、占17.68%，GO庫注釋的Unigene數(shù)量為35 244個、占14.79%，Swiss-Prot注釋的Unigene數(shù)量為46 337個、占19.44%，Pfam注釋的Unigene 數(shù)量為58 605個、占24.59%（表3）。

2.2.1 GO 分類 對Unigene進行GO注釋后，按照成功注釋到生物學(xué)過程（Biological Process，BP）、細胞組分（Cellular Component，CC）和分子功能（Molecular Function，MF）的基因進行統(tǒng)計。其中生物學(xué)過程中較多Unigene注釋有高分子代謝過程28 326個、初級代謝過程16 483個、細胞代謝過程38 883個；分子功能中較多Unigene注釋有轉(zhuǎn)移酶活性5 090個、水解酶活性5 153個、跨膜轉(zhuǎn)運功能3 069個；細胞組分中較多Unigene注釋有細胞壁8 663個、細胞膜7 832個、細胞器9 581個（圖2）。

表3 Unigenes 注釋結(jié)果統(tǒng)計Table 3 Annotation statistics for all Unigenes

圖2 GO分類Fig.2 GO function classification

2.2.2 KEGG分類 由表3可知，蘭州百合測序結(jié)果有26 291條Unigene注釋到KEGG數(shù)據(jù)庫，涉及5大類代謝通路。其中代謝作用中注釋4 914條、占40.66%；有機系統(tǒng)中注釋905條、占7.49%；細胞過程中注釋923條、占7.64%；遺傳信息處理中注釋3 609條、占29.86%；環(huán)境信息處理中注釋1 734條、占14.35%。

2.3 基因表達差異分析結(jié)果

2.3.1 差異基因GO富集分析 將所有差異表達基因映射到Gene Ontology 數(shù)據(jù)庫的各個功能分類中，計算每個分類的基因數(shù)目，與整體基因組背景相比，找出在差異表達基因中顯著富集的功能分類。

對上調(diào)的差異基因進行富集分析，結(jié)果（圖4）發(fā)現(xiàn)，在細胞成分中，差異基因富集效果較明顯的有細胞核，基因占比13.81%；細胞內(nèi)細胞器，基因占比16.06%；面上細胞器，基因占比16.06%。分子功能中差異基因富集效果較明顯的功能有鈣離子結(jié)合，基因占比3.52%；離子束縛，基因占比8.81%；陽離子綁定，基因占比8.81%；金屬離子結(jié)合，基因占比8.81%。生物過程中上調(diào)基因富集效果較明顯的功能有水響應(yīng)，基因占比1.27%；非生物刺激反應(yīng)，基因占比1.37%；磷酸化作用，基因占比8.62%。結(jié)果表明，含量增加的差異基因在細胞成分、分子功能、生物過程三方面主要調(diào)節(jié)物質(zhì)的合成代謝。

圖3 KEGG 代謝通路Fig.3 KEGG metabolic pathway

對下調(diào)的差異基因富集分析，結(jié)果（圖4）發(fā)現(xiàn)，差異基因主要富集在生物過程中。其中富集較為明顯的是光合作用、光收獲，基因占比1.65%；光合作用、光反應(yīng)，基因占比1.65%。結(jié)果表明，含量相對減少的差異基因主要調(diào)節(jié)生物過程中的光合作用過程。

圖4 差異基因（20℃ Vs -4℃）的GO富集結(jié)果Fig.4 GO enrichment result of differential genes (20℃ Vs -4℃)

2.3.2 差異基因KEGG富集分析 差異基因KEGG富集結(jié)果顯示，常溫（20℃）與超低溫（-4℃）下差異表達基因共有5 848個，其中上調(diào)的差異基因有3 478個，下調(diào)的差異基因有2 370個（圖5）。經(jīng)分析后發(fā)現(xiàn)，與蛋白激酶、蛋白磷酸酶、碳代謝以及ABA等相關(guān)的基因顯著變化；此外，吲哚-3-甘油磷酸合成酶、蛋白磷酸酶、已糖激酶、鈣結(jié)合蛋白、葉綠素a/b結(jié)合蛋白等相關(guān)基因顯著上調(diào)或下調(diào)。

圖5 差異表達基因（20℃ Vs -4℃）火山圖Fig.5 Differential expressed genes（20℃ Vs -4℃）

Pathway顯著性富集分析以KEGG Pathway為單位，應(yīng)用超幾何檢驗，找出差異基因相對于所有注釋的基因顯著富集的Pathway。處理后含量增加的差異基因顯著富集的通路有植物病原菌互作、氨基酸和核苷酸糖代謝、植物晝夜節(jié)律、亞油酸代謝等（圖6）。如上調(diào)的差異基因c156798_g2、c163371_g2 (calcium-dependent protein kinase，EC:2.7.11.1)、c166423_g2、c145559_g1 (respiratory burst oxidase，EC:1.6.3.- 1.11.1)、c162064_g1、c164545_g2、c166855_g1(3.1027)c145934_g1(8.2812) c101051_g1(4.9356)、c144017_g1(Inf) c144017_g2(Inf) c159000_g1(2.9702)c140553_g1(9.9651) c170503_g1(Inf)、c140013_g1(3.7728)、c154815_g1(Inf) c173122_g1(7.8456)c168962_g1(Inf) c163915_g1(Inf)、c156536_g2(2.9685) c156536_g1(1.2021)（calmodulin）、c147877_g2（mitogen-activated protein kinase kinase 4/5，EC:2.7.12.2）主要調(diào)節(jié)植物病原互作，病原互作相關(guān)基因與植物的抗性息息相關(guān)；上調(diào)的差異基因c154320_g1、c125968_g1(5.6972)c171764_g2(5.6213)、c172816_g3(5.9018)、c237479_g1(3.5347) c161991_g1(2.2688)、c172816_g4(2.2231)、c113669_g1(5.4648) c161991_g2(1.9362)、c106411_g1調(diào)節(jié)氨基酸和核苷酸糖代謝；上調(diào)的差異基因c151010_g1、c156146_g2(2.1536) c171563_g1(2.3117)、c156146_g1、c167910_g3、c152959_g1、c131089_g1(Inf)c144654_g1(6.9885)主要調(diào)節(jié)植物晝夜節(jié)律。

處理后下調(diào)的差異基因顯著富集的通路有光合作用、光合作用-天線蛋白、卟啉和葉綠素代謝、光合生物體中的碳固定、脂肪酸生產(chǎn)等。如下調(diào)的差異基因c237057_g1、c13699_g1、c237057_g1、c76799_g1、c158296_g1、c152833_g1、c152881_g1、c158839_g1、c134148_g1、c140011_g1、c153193_g1、c131998_g1、c48634_g1、c71692_g1、c42744_g1、c153649_g1調(diào)節(jié)光合作用；下調(diào)的差異基因c149517_g1、c139230_g1、c152435_g1、c149445_g1、c153860_g1、c115377_g1、c164577_g1、c155535_g1、c174985_g1、c157047_g1調(diào)節(jié)光合作用-天線蛋白；下調(diào)的差異 基 因 c167743_g1、c167947_g1、c167743_g1、c71809_g1、c165450_g1c168519_g1、c169028_g1c133188_g1、c123480_g1、c185147_g2、c157598_g1、c157598_g1、c155686_g1、c166557_g2(-1.4097) c166557_g1(-1.5029)調(diào)節(jié)卟啉調(diào)節(jié)葉綠素代謝；下調(diào)的差異基因c104889_g2(-1.5635)c174574_g3(-3.7395)、c104889_g1、c71483_g1、c159323_g1、c172966_g1、c170442_g1、c162112_g2、c154502_g4、c43883_g1、c198353_g1、c168133_g3、c164307_g1調(diào)節(jié)光合生物體中的碳固定；下調(diào)的差異基因c134603_g2、c164323_g1、c166224_g1(-1.1149) c150017_g1(-1.3389)、c163278_g2、c148702_g1、c164323_g2、c134778_g1、c117604_g1調(diào)節(jié)脂肪酸伸長（圖7）。

2.4 qRT-PCR 驗證

為驗證RNA-seq測序結(jié)果的可信度，隨機選取10個差異表達基因，其中含上調(diào)的5條差異基因（c166605_g1、c199368_g1、c171769_g1、c164813_g2、c153721_g1）和下調(diào)的5條差異基 因（c76799_g1、c153860_g1、c149445_g1、c115377_g1、c155535_g1），以Actin為內(nèi)參進行qRT-PCR驗證。由圖8可知，10個基因在低溫脅迫下的表達程度不同，但表達趨勢與高通量測序結(jié)果基本一致，即表明測序結(jié)果真實可靠。

圖6 上調(diào)差異基因（20℃ Vs -4℃）參與的KEGG代謝通路富集結(jié)果Fig.6 Differential up-regulation genes（20℃ Vs -4℃）assignment to KEGG metabolic pathway

圖7 下調(diào)差異基因（20℃ Vs -4℃）參與的KEGG代謝通路富集結(jié)果Fig.7 Differential down-regulation genes（20℃ Vs -4℃）assignment to KEGG metabolic pathway

圖8 差異基因的qRT-PCR驗證Fig.8 Validation of DEGs by qRT-PCR

3 討論

前人關(guān)于百合抗凍分子機制研究較少，不良生長環(huán)境尤其是零下低溫脅迫是植物生長發(fā)育中常遇到的自然災(zāi)害。本研究采用先進的測序平臺Illumina 4000高通量測序平臺對常溫（20℃）與零下低溫（-4℃）處理的蘭州百合進行測序，確保測序數(shù)據(jù)的準(zhǔn)確性和序列的高度豐富性［18］。測序結(jié)果在NCBI數(shù)據(jù)庫中共注釋到了238 301條Unigene，但注釋到每個庫中的基因最大不超過27.41%，說明還有很多基因沒有被注釋到NCBI庫中，原因可能是數(shù)據(jù)庫信息不足、測序結(jié)果中有新基因和基因片段短等。因此，在后續(xù)的研究中，我們可以進一步比對測序結(jié)果中差異極顯著的基因片段，進一步深入發(fā)掘蘭州百合抗凍轉(zhuǎn)錄新基因。

根據(jù)KEGG富集分析結(jié)果，可以看出共有25條差異基因參與病原互作途徑，10條差異基因參與氨基酸和核苷酸糖代謝途徑，8條差異基因參與植物晝夜節(jié)律途徑，15條差異基因參與卟啉調(diào)節(jié)葉綠素代謝途徑，16條差異基因參與光合作用途徑，10條差異基因調(diào)節(jié)光合作用-天線蛋白途徑，13條差異基因調(diào)節(jié)光合生物體中的碳固定途徑，9條差異基因調(diào)節(jié)脂肪酸伸長途徑。因此我們可以深入研究以上通路以進一步揭示蘭州百合抗凍的分子機制。另外也可以將蘭州百合相關(guān)的代謝產(chǎn)物作為研究對象進行深入探究。

4 結(jié)論

本研究對蘭州百合轉(zhuǎn)錄組測序后通過KEGG通路分析，結(jié)果檢測到24 465個差異表達的基因，篩選后獲得8 558條差異表達基因，共占差異表達基因總數(shù)的10.27%。經(jīng)富集分析后發(fā)現(xiàn)，與蛋白激酶、蛋白磷酸酶、碳代謝以及ABA等相關(guān)的基因可以作為研究蘭州百合寒冷響應(yīng)機制的主要研究對象。qRT-PCR分析結(jié)果表明，隨機選出的10個差異性表達基因的表達趨勢與高通量測序結(jié)果相一致。此外，還候選了吲哚-3-甘油磷酸合成酶、蛋白磷酸酶、已糖激酶、鈣結(jié)合蛋白、葉綠素a/b結(jié)合蛋白等基因作為與蘭州百合零下低溫響應(yīng)相關(guān)的應(yīng)答候選基因，為揭示蘭州百合抗凍分子機制奠定了基礎(chǔ)。