轉錄組測序研究蘭州百合抗凍關鍵基因及途徑

田雪慧 ，郁繼華，頡建明

（1.楊凌職業技術學院生態環境工程分院，陜西 楊凌 712100；2.甘肅農業大學園藝學院，甘肅 蘭州 730070）

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

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

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

1 材料與方法

1.1 試驗材料

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

1.2 試驗方法

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

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

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

2 結果與分析

2.1 測序結果與Trinity拼接

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

表2 處理前后序列數量、N50、N75、N90長度和總堿基數統計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 基因功能注釋結果

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

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

表3 Unigenes 注釋結果統計Table 3 Annotation statistics for all Unigenes

圖2 GO分類Fig.2 GO function classification

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

2.3 基因表達差異分析結果

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

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

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

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

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

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

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

Pathway顯著性富集分析以KEGG Pathway為單位，應用超幾何檢驗，找出差異基因相對于所有注釋的基因顯著富集的Pathway。處理后含量增加的差異基因顯著富集的通路有植物病原菌互作、氨基酸和核苷酸糖代謝、植物晝夜節律、亞油酸代謝等（圖6）。如上調的差異基因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）主要調節植物病原互作，病原互作相關基因與植物的抗性息息相關；上調的差異基因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調節氨基酸和核苷酸糖代謝；上調的差異基因c151010_g1、c156146_g2(2.1536) c171563_g1(2.3117)、c156146_g1、c167910_g3、c152959_g1、c131089_g1(Inf)c144654_g1(6.9885)主要調節植物晝夜節律。

處理后下調的差異基因顯著富集的通路有光合作用、光合作用-天線蛋白、卟啉和葉綠素代謝、光合生物體中的碳固定、脂肪酸生產等。如下調的差異基因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調節光合作用；下調的差異基因c149517_g1、c139230_g1、c152435_g1、c149445_g1、c153860_g1、c115377_g1、c164577_g1、c155535_g1、c174985_g1、c157047_g1調節光合作用-天線蛋白；下調的差異 基 因 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)調節卟啉調節葉綠素代謝；下調的差異基因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調節光合生物體中的碳固定；下調的差異基因c134603_g2、c164323_g1、c166224_g1(-1.1149) c150017_g1(-1.3389)、c163278_g2、c148702_g1、c164323_g2、c134778_g1、c117604_g1調節脂肪酸伸長（圖7）。

2.4 qRT-PCR 驗證

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

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

圖7 下調差異基因（20℃ Vs -4℃）參與的KEGG代謝通路富集結果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 討論

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

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

4 結論

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