李冠嶸，何 好，朱國慶，陳詩雅，徐 陽，金淑梅

（東北鹽堿植被恢復與重建教育部重點實驗室，東北林業(yè)大學生命科學學院，哈爾濱150040）

0 引言

RNA-Seq轉錄組測序是分析基因表達水平、篩選差異表達基因、挖掘功能基因的重要手段，很多高等植物都進行了耐鹽相關的轉錄組方面的研究。對耐鹽棉花進行轉錄組測序，分析了與鹽脅迫相關的基因[1]；通過轉錄組測序揭示了鹽脅迫下柑橘根中基因表達的情況[2]；對在鹽堿脅迫下高粱的轉錄組測序鑒定了與ABA相關的差異表達基因，推測了參與ABA信號傳導的家族成員[3]；對煙草的轉錄組測序確定了分別響應NaCl和NaHCO3差異表達的miRNA，并確定了在NaCl和NaHCO3處理下的mRNA-miRNA相互作用[4]；通過轉錄組分析3種不同表型的陸地棉鹽堿脅迫處理后反應的分子機制[5]。

百合的轉錄組學主要是在發(fā)育方面，如利用微陣列技術研究了麝香百合生殖細胞中基因的表達情況[6]，通過Illumina HiSeqTM2000測序平臺研究了百合花粉萌發(fā)過程基因的表達[7]，在花色方面的研究主要有，為了研究與百合花顏色合成有關的調(diào)控信息，研究了東方百合和亞洲百合春花過程中的分子響應機制[8-10]；通過轉錄組測序，揭示了花色苷生物合成和運輸所需表達的基因[11]；Li等[12]根據(jù)百合轉錄組測序結果研究了3種植物激素的信號傳導途徑。利用Roche454測序平臺獲得了4種基因型百合葉片的轉錄組測序信息[13]，在逆境方面也有研究，建立了卷丹冷脅迫鹽分、激素壓力下的轉錄譜[14-15]。

近幾年，對細葉百合的轉錄組也有研究，對低溫處理0天、30天、60天、90天后的細葉百合鱗莖進行轉錄組測序，探究鱗莖休眠解除的分子機制[16]。對參與細葉百合休眠解除過程中的WRKY差異表達基因進行篩選，詳細解析了LpWRKY20基因在細葉百合休眠解除過程中的調(diào)控作用[17]?？寺〉玫組ADS-box家族中的LpSVP闡明LpSVP基因的作用機制[18]。篩選得到21個 NAC基 因 ，并 得 到 含 LpNAC5、LpNAC13和LpNAC21基因的陽性植株進行研究[19]。

細葉百合(Lilium pumilum DC.)是百合科百合屬的多年生草本植物，花美姿艷。不僅是重要的觀賞性花卉，也具有很高的食用、藥用等經(jīng)濟價值。大多數(shù)野生百合植株嬌小細弱，其抗逆性不強，這極大地限制了百合在土壤鹽堿化問題嚴重的地區(qū)的種植和應用[20]。細葉百合是在鹽堿化土壤中生存的少數(shù)幾種百合之一，是百合的優(yōu)良種質資源之一，具有抗寒、耐旱、耐鹽堿等優(yōu)良的特性。細葉百合采自于鹽堿土壤比較嚴重的黑龍江省大慶市地區(qū)，對于能在惡劣的鹽堿條件下生長，其基因的表達有什么特殊之處呢？為了探究這個問題，對細葉百合在碳酸鹽逆境下的轉錄組數(shù)據(jù)進行了詳細的分析，希望揭示細葉百合在鹽堿逆境中的基因表達方式，為合理利用鹽堿地和大面積種植百合提供理論基礎。

1 材料與方法

1.1 試驗材料與處理

以培養(yǎng)在組培瓶中的一年生細葉百合為材料，將長勢一致的細葉百合幼苗移植到1/2MS培養(yǎng)基(CK)和含有20 mmol/L NaHCO3（處理組）的MS培養(yǎng)基中，處理48 h后，將對照組(CK)和處理組(NaHCO3)各兩組分別取樣，液氮速凍后置于-80℃保存。試驗在東北林業(yè)大學東北鹽堿植被恢復與重建教育部重點實驗室

1.2 細葉百合葉片總RNA的提取、文庫構建、Illumina測序及數(shù)據(jù)質量控制

2次重復的對照組和處理組的細葉百合鱗莖cDNA文庫的構建以及測序實驗由蘇州金唯智生物科技有限公司完成。具體實驗操作流程如下：（1）提取到總RNA，mRNA的富集、cDNA鏈的合成。（2）純化、洗脫、末端修復后進行PCR擴增。（3）使用Illumina HiSeq高通量測序平臺對cDNA文庫進行測序，對測序結果原始圖像數(shù)據(jù)利用軟件bcl2fastq進行圖像堿基識別，初步質量分析后得到原始測序數(shù)據(jù)。分析測序數(shù)據(jù)進行質量評估，使用二代測序數(shù)據(jù)質量統(tǒng)計軟件Cutadapt對測序原始數(shù)據(jù)去除接頭以及低質量序列，得到后續(xù)信息分析用的Clean Data。

1.3 測序數(shù)據(jù)組裝及基因功能注釋

采用組裝軟件Trinity對樣品數(shù)據(jù)從頭組裝，組裝結果通過序列聚類做進一步序列拼接和去冗余處理，得到長的非冗余的unigene序列并與Nr、Swiss-Prot、GO、KEGG、COG數(shù)據(jù)庫進行比對，獲得mRNA的注釋信息。

1.4 序列比對以及差異表達基因(DEG)的鑒定

使用bowtie2軟件進行短reads的比對，將Clean data比對Unigene上，進行分析，并RSEM軟件使用FPKM方法計算基因的表達量。采用Bioconductor軟件包的DESeq2進行樣品組間的差異表達檢測分析，以差異基因表達變化2倍以上且FDR≤0.05為標準進行篩選，得到最終的差異表達基因，選取9個差異表達基因進行QPCR驗證。

1.5 實時熒光定量PCR驗證

為了驗證RNA-Seq中表達量數(shù)據(jù)的準確性，對細葉百合進行20 mmol/L NaHCO3處理24 h后提取細葉百合鱗莖RNA，將RNA反轉錄成cDNA作為模板，選取了9個（5個上調(diào)基因和4個下調(diào)基因）與鹽堿脅迫密切相關的基因(BEL1、ERD6、GRP、MCM、FRO8、CPC、PPT、LBP、PEX7)設計熒光定量PCR引物（表1）進行qRT-PCR驗證。UltraSYBR Mixture試劑進行實時熒光定量qRT-PCR檢測，反應體系為20 μL（2×SYBR Green Mix 10 μL；10 μmol/L 上下游引物各0.5 μL；cDNA 模板 1 μL；ddH2O 8 μL）；反應條件為：95℃預變性10 min，95℃ 30 s，55℃ 30 s，72℃ 1 min，40個循環(huán)。得到細葉百合測序的轉錄組中上調(diào)和下調(diào)基因的表達量并進行分析。

表1 qRT-PCR所用引物序列

2 結果與分析

2.1 測序結果統(tǒng)計分析

對4個(CK1、CK2、Salt1、Salt2)cDNA文庫進行的測序，通過質量控制后，共得到23.4 Gb Clean Data，各樣品Q30堿基百分比均不小于89.35%。

Clean Data經(jīng)組裝后得到115061條unigene，N50長度為939 nt，115061條unigene中有65.18%的長度在200～500之間，具體信息見表3與圖1。

圖1 Unigene不同長度的占比情況

表3 經(jīng)組裝后的unigene情況

2.2 基因功能注釋結果

將Unigene序列與Nr、COG、Swissprot、KEGG數(shù)據(jù)庫進行Blast比對后，共獲得56828個mRN的注釋信息。55433個Unigene被注釋到Nr數(shù)據(jù)庫中，26973個Unigene被注釋到COG數(shù)據(jù)庫中，23610Unigene被注釋到GO數(shù)據(jù)庫中，37434個Unigene被注釋到Swissprot數(shù)據(jù)庫中，13142個Unigene被注釋到KEGG數(shù)據(jù)庫中。

表2 樣品測序數(shù)據(jù)統(tǒng)計表

2.3 Nr注釋結果

共55433個Unigene被注釋到Nr數(shù)據(jù)庫中，其中，注釋到油棕(Elaeis guineensis)的Unigene最多，占總數(shù)的25%，其余注釋到海棗(Phoenix dectylifera)、小果野蕉(Musa acuminatasubsp.malaccensis)、鳳梨(Ananas comosus)、蓮花(Nelumbo nucitera)、葡萄(Vitis vinifera)的Unigene分別占總數(shù)的19%、9%、5%、2%、2%。具體信息見圖2。

圖2 Nr注釋分類圖

2.4 KOG注釋結果

共有26973條Unigene被注釋到KOG數(shù)據(jù)庫的26個分類中，其中數(shù)量最多的一類是一般功能類(R:General function prediction only，共3478條，12.89%)，其次是信號轉導機制（T:Signal transduction mechanisms，共 3279條，12.16%），翻譯后修飾（O:Posttranslational modificadification，共 3218 條 ，11.93%），具體信息見圖3。

圖3 KOG注釋結果

2.5 GO注釋結果

共有23610條Unigene被注釋到GO數(shù)據(jù)庫的3個大類與53個小類中，在分子功能(molecular function)大類中，注釋到催化活性(catalytic)與結合綁定(binding)分類的數(shù)量最多；在細胞組分(Cellular Component)大類中，注釋到細胞部分(cell part)與細胞器(organelle)分類的數(shù)量最多；在生物學過程(Biological process)大類中，注釋到代謝過程(metabolic process)及細胞生理過程(cellular process)中的Unigene數(shù)量最多。具體信息見圖4。

圖4 GO注釋結果

2.6 KEGG注釋結果

共13142條Unigene被注釋到KEGG數(shù)據(jù)庫的五大類中，其中注釋到代謝(Metabolism)這一大類的Unigene數(shù)量最多，占總數(shù)的38%，其余注釋到的大類中，有機體系(Organismal Systems)占20%，遺傳信息處理(Genetic Information Processing)占17%，環(huán)境信息處理(Environmental Information Processing)占14%，細胞過程(Cellular Processes)占11%。具體信息見圖5。

圖5 KEGG注釋結果

2.7 DEG的分析

經(jīng)過NaHCO3脅迫處理細葉百合植株24 h后，處理組和對照組鱗莖中的基因表達對比，共鑒定出390個差異表達基因，其中有155個基因上調(diào)，235個基因下調(diào)，具體信息見圖6。

圖6 差異表達基因DEG統(tǒng)計

2.8 熒光定量PCR驗證結果

qRT-PCR的結果顯示見圖7，BEL1基因與對照相比表達量降低了90%，ERD6基因與對照相比表達量降低了85%，GRP基因與對照相比表達量降低了43%，MCM基因與對照相比表達量降低了38%，PPT基因與對照相比表達量增加了47%，LBP基因與對照相比表達量增加了48%，F(xiàn)RO8基因與對照相比表達量增加了53%，CPC基因與對照相比表達量增加了176%，PEX7基因與對照相比表達量增加了191%。所選取的9個基因的表達結果與轉錄組的結果基本趨于一致，表明RNA-Seq的結果準確可信。

圖7 熒光定量PCR驗證轉錄組可靠性

3 討論與結論

百合的轉錄組在發(fā)育[6-7]和花色[8-11]方面研究比較多，對激素信號途徑的百合轉錄組也有研究[11-12]，在逆境方面研究了卷丹在冷脅迫下的基因表達，并且篩選出與鹽堿逆境有關的基因[14-15]，在低溫處理下研究了細葉百合鱗莖轉錄組[16]，本研究同樣是以細葉百合的鱗莖做為植物材料進行的轉錄組方面的研究，但對于鹽堿脅迫下百合尤其是堿性鹽脅迫下的基因轉錄組研究未見報道，本文的創(chuàng)新之處是在堿性鹽(NaHCO3)處理后轉錄組的研究，通過數(shù)據(jù)分析，力圖找到處理組和對照組中差異表達的基因，觀察差異基因的表達，為了探究細葉百合耐鹽堿的分子機理提供理論基礎。

在本研究轉錄組的差異數(shù)據(jù)中除了一些常見的與鹽堿有關的基因（脯氨酸轉運蛋白，過氧化物酶體生物合成酶）外，還出現(xiàn)幾個與糖轉運（蔗糖轉運蛋白）及金屬運輸（鐵還原氧化酶，藍銅蛋白）等基因，選取了9個與基因進行qRT-PCR驗證，BELL基因是參與調(diào)控植物的發(fā)育過程的轉錄因子，3% NaCl處理條件下一些基因的表達水平明顯上調(diào)，尤其NtBELL6基因的表達水平上調(diào)了3倍[21]，結果顯示BEL1基因與對照相比表達量降低了90%，BELL1的表達在NaHCO3處理后下調(diào)，是否在其他NaHCO3的濃度下，該基因的表達會有所不同，需要進一步的研究。ERD6蛋白屬于蔗糖轉運蛋白，發(fā)現(xiàn)AtERD6能夠快速地響應低溫誘導，在低溫處理1 h后表達顯著上調(diào)[22]，而茶樹中的CsERD6-7冷馴化過程中的表達顯著被抑制[23]。在轉錄組測序和qPCR驗證中，ERD6基因與對照相比表達量降低了85%，在NaHCO3處理下，基因的表達量下調(diào)，其作用機理還有待進一步的研究。GRP富含甘氨酸的RNA結合蛋白，鹽芥葉中的TsGRP表達量在300 mmol/L NaCl處理后表達下調(diào)，根中的表達量上調(diào)[24]，在轉錄組測序和qPCR驗證中，GRP基因與對照相比表達量降低了43%。是否是因為所選的脅迫濃度或者選取的器官不同，基因的表達量有所不同。還有待進一步研究。MCM編碼真核生物的轉錄因子家族，水分脅迫下大豆籽粒中發(fā)現(xiàn)被富集到DNA復制通路上的差異表達基因MCM[25]本轉錄組中MCM基因與對照相比表達量降低了38%，表達量降低的原因還有待進一步研究。

鐵還原氧化酶(ferric reduction oxidase，F(xiàn)RO)基因家族參與植物中廣泛存在的各種生物過程，AtFRO8為線粒體FRO家族基因成員，在不同發(fā)育階段參與植株體內(nèi)Fe3+螯合物的還原[26]，F(xiàn)RO8基因與對照相比表達量增加了13%；CPC是一種藍銅蛋白，影響生物的生理發(fā)育活動。擬南芥藍銅蛋白能催化離子的氧化，參與了在細胞膜區(qū)域的電子轉移反應來提高植物的逆境脅迫能力[27]本研究中CPC基因與對照相比表達量增加了22%；脯氨酸轉運蛋白（ProT或PPT）是典型的Na+依賴型亞氨基酸轉運蛋白類，在NaCl的脅迫下，馬鈴薯脯氨酸轉運蛋白2個基因在2 h就開始急劇響應，其表達量均是對照的19倍以上[28]，本轉錄組所得到的數(shù)據(jù)與文獻中的結果相一致，PPT基因與對照相比表達量增加了47%；BIP為熱激蛋白家族成員，是內(nèi)質網(wǎng)的一個主要的分子伴侶，大豆BIP能增強抗旱能力，減弱由干旱誘導的葉片衰老[29]本研究中BIP基因與對照相比表達量增加了48%，過氧化物酶體生物合成酶(PEX)是過氧化物酶體在生物體內(nèi)合成的關鍵酶，胡楊PePEX11能夠提高擬南芥在鹽脅迫下的抗氧化能力，提升耐鹽能力[30]PEX7基因與對照相比表達量增加了59%。這些在轉錄組中表達量上調(diào)的與文獻中基因在逆境處理后，表達量上調(diào)的結果是一致的，說明細葉百合，BELL、ERD、GRP、FRO、CPC、PPT、LBP、PEX7基因與鹽有一定的應答關系，為以后進一步深入研究這些基因的功能奠定一定的基礎。

通過對細葉百合鱗莖進行轉錄組測序，共獲得23.4 Gb的數(shù)據(jù)，通過組裝得到115061條unigene。將Unigene序列與Nr、COG、Swissprot、KEGG數(shù)據(jù)庫進行Blast比對后，共獲得56828個mRNA的注釋信息。其中，55433個Unigene被注釋到Nr數(shù)據(jù)庫中，其中，注釋到油棕的Unigene最多，占總數(shù)的25%，其余注釋到海棗、小果野蕉、鳳梨、蓮花、葡萄的Unigene分別占總數(shù)的19%、9%、5%、2%、2%。26973條Unigene被注釋到KOG數(shù)據(jù)庫的26個分類中，其中數(shù)量最多的一類是一般功能類，其次是信號轉導機制。有23610條Unigene被注釋到GO數(shù)據(jù)庫的3個大類與53個小類中，在分子功能(molecular function)大類中，注釋到催化活性與結合綁定分類的數(shù)量最多；在細胞組分(Cellular Component)大類中，注釋到細胞部分與細胞器分類的數(shù)量最多；在生物學過程(Biological process)大類中，注釋到代謝過程及細胞生理過程中的Unigene數(shù)量最多。13142條Unigene被注釋到KEGG數(shù)據(jù)庫的五大類中，其中注釋到代謝(Metabolism)這一大類的Unigene數(shù)量最多，占總數(shù)的38%，其余注釋到的大類中，有機體系(Organismal Systems)占20%，遺傳信息處理(Genetic Information Processing)占17%，環(huán)境信息處理(Environmental Information Processing)占14%，細胞過程(Cellular Processes)占11%。在兩組轉錄組數(shù)據(jù)分析中，共鑒定390個差異表達基因，其中有155個基因上調(diào)，235個基因下調(diào)。qRT-PCR結果顯示所選取的9個基因的表達結果與轉錄組的結果基本趨于一致，結果說明了RNA-Seq數(shù)據(jù)的可靠性。細葉百合在碳酸鹽(NaHCO3)逆境下的轉錄組數(shù)據(jù)對比，提供了許多的基因表達通路和表達量的差異，希望為細葉百合在鹽堿逆境中的基因表達方式研究提供理論基礎。