邢建宏，陳 偉

（1.三明學(xué)院資源與化工學(xué)院/ 福建省資源環(huán)境監(jiān)測與可持續(xù)經(jīng)營利用重點實驗室，福建 三明 365004；2.福建農(nóng)林大學(xué)生命科學(xué)學(xué)院，福建 福州 350012）

近年來土壤鹽堿化問題日益突出，對農(nóng)業(yè)可持續(xù)發(fā)展造成重大威脅[1].因此，深入研究植物的耐鹽機理，挖掘耐鹽基因資源對改變?nèi)遮厙乐氐耐寥利}堿化問題有重要意義.紅樹植物生長在潮間海岸帶，受海水的周期性浸漬，形成了一套有別于淡水植物或陸生植物的耐鹽機制[2]，值得深入研究，并挖掘其中有價值的基因資源.Ilumina 高通量測序技術(shù)的出現(xiàn)，為從基因組和轉(zhuǎn)錄組層面深入解析植物耐鹽分子機理提供了有力手段.近年來研究者利用高通量測序技術(shù)對紅樹（Rhizophora apiculata）[3]和木欖（Bruguiera gymnorrhiza）[4]進行轉(zhuǎn)錄組分析，結(jié)果表明抗逆基因的表達可能是提高其抗極端環(huán)境能力的關(guān)鍵因素.對泌鹽紅樹植物海欖雌（Avicennia officinalis）鹽脅迫下根系的RNA-seq 分析表明，ABA 和乙烯介導(dǎo)的植物激素信號通路在其響應(yīng)鹽脅迫中發(fā)揮重要作用[5].這些研究為從轉(zhuǎn)錄組層面揭示紅樹植物耐鹽機理提供了前期探索，但對非泌鹽植物根系的耐鹽分子調(diào)控機制還未見系統(tǒng)和深入的研究.秋茄（Kandelia candel）是紅樹植物中典型的非泌鹽植物[6]，其耐鹽能力通過自身基因表達調(diào)控實現(xiàn)，蘊藏著豐富的耐鹽基因資源，因此從分子水平解析秋茄的耐鹽機理可為揭示木本植物耐鹽機制提供新的證據(jù).本研究利用Illumina HiSeqTM2000 平臺對不同濃度NaCl 處理下秋茄幼苗根系轉(zhuǎn)錄組進行測序，并建立轉(zhuǎn)錄組數(shù)據(jù)庫，為今后從分子層面深入解析其耐鹽機制提供基礎(chǔ)數(shù)據(jù).

1 材料與方法

1.1 材料

供試秋茄胚軸采自福建省漳州市漳江口紅樹林國家自然保護區(qū)（23°55′N,117°26′E）.挑選大小與成熟度相近、無機械損傷和病蟲害的胚軸種植于45 cm×35 cm×25 cm 塑料盆的干凈河沙中，每盆20 株，每盆澆灌1L Hoagland 營養(yǎng)液，每天傍晚用自來水補充散失的水分.

1.2 方法

1.2.1 NaCl 處理 待幼苗長至4 葉后（60 d），用含不同鹽濃度（0、200、400、600 mmol·L-1NaCl）的Hoagland 營養(yǎng)液澆灌幼苗根部，為避免原有水分對處理濃度的影響，處理前2 d 先放干培養(yǎng)容器中的溶液，停止水分補充.為避免高濃度鹽處理引起鹽激效應(yīng)，400、600 mmol·L-1NaCl 處理均從200 mmol·L-1開始增加，每隔1 d 施加1 次鹽處理，并在同一時間達到相應(yīng)處理濃度，3 d 后取樣.每個處理選取20 株幼苗，剪取所有植株的根，迅速放入液氮中速凍20 min 后，置于-80 ℃冰箱保存?zhèn)溆?每個處理2 次生物學(xué)重復(fù).

1.2.2 總RNA 的提取及文庫構(gòu)建 采用十六烷基三甲基溴化銨（CTAB）-LiCl 沉淀法提取秋茄根系中總RNA[7].將4 種不同NaCl 濃度處理的秋茄幼苗根系單獨提取RNA，每個處理2 個生物學(xué)重復(fù)，準確測定所提取RNA 濃度，將所有樣品RNA 等量混合，構(gòu)建含有不同處理濃度的RNA 樣品池進行測序.將檢測合格（濃度、純度、完整性和片段大小符合建庫要求）的RNA 按照Illumina 使用說明書構(gòu)建cDNA 文庫.首先將RNA 通過片段緩沖液(Ambion，Austin，TX，USA)進行孵育并破碎；接著應(yīng)用隨機六聚體引物hexamer-primers (Illumina，San Diego，CA，USA) 和逆轉(zhuǎn)錄酶(Invitrogen，Carlsbad，CA，USA) 合成第一鏈cDNA；然后應(yīng)用DNA 聚合酶I(Invitrogen，Grand Island，NY，USA)合成第二鏈cDNA；將雙鏈cDNA 片段應(yīng)用QIAquick PCR extraction kit 進行純化處理，并連接Illumina 轉(zhuǎn)錄適配器(Illumina，San Diego，CA,USA)，隨后選擇合適的片段進行PCR 擴增.

1.2.3 Illumina 測序、組裝及注釋分析 轉(zhuǎn)錄組文庫質(zhì)量檢測合格后，應(yīng)用Illumina HiSeqTM2000(Illumina，San Diego，CA，USA)儀器按照Illumina 公司的標準操作步驟進行轉(zhuǎn)錄組測序.對原始reads（雙端序列）進行評估，剔除低質(zhì)量的reads（more than 20% of the bases showed a Q-value≤10）和含有超過5%未知序列“[N]”的reads，獲得clean reads.隨后采用Trinity（v2012-10-05）軟件對去雜后的clean reads進行轉(zhuǎn)錄本De novo 組裝，獲得Unigene 庫（Unigene Library），利用Nonredundant protein (Nr)、Kyoto Encyclopedia of Genes and Genomes (KEGG)、Clusters of Orthologous Groups（COG）、Gene Ontology (GO)、Swiss-Prot 等數(shù)據(jù)庫對每條Unigene 進行BLAST 比對，獲得每個數(shù)據(jù)庫的比對結(jié)果及功能注釋.

1.2.4 耐鹽相關(guān)基因的qRT-PCR 驗證 為證明本研究所得數(shù)據(jù)庫包含了秋茄耐鹽相關(guān)基因[8]，選取前人發(fā)現(xiàn)的6 個紅樹植物鹽脅迫相關(guān)基因，以本研究獲得的序列信息設(shè)計引物，對不同NaCl 濃度處理下的基因表達量進行qRT-PCR 驗證.所驗證基因信息及引物見表1，方法見參考文獻[8].

表1 不同濃度NaCl 處理下差異表達基因qRT-PCR 驗證引物

續(xù)表 1

2 結(jié)果與分析

2.1 總RNA 提取與質(zhì)量檢測

分別對4 種NaCl 濃度處理下的8 個樣品的總RNA 進行瓊脂糖凝膠電泳檢測，結(jié)果（圖1）顯示，電泳條帶明亮、清晰，無拖尾現(xiàn)象，28S 的亮度在18S 的兩倍以上，表明所提取的RNA 質(zhì)量較高，符合建庫要求.

圖1 瓊脂糖電泳初步檢測總RNA 質(zhì)量

2.2 測序結(jié)果評估與序列組裝

由表2 可知，對不同NaCl 濃度處理下秋茄根系樣品進行轉(zhuǎn)錄組測序總共得到了65.79 M 的clean reads，每個樣品的clean reads 均超過7.50 M；測序總堿基數(shù)為6.64 Gb，每個樣品的總堿基數(shù)均高于760 Mb；各樣品堿基Q30 均大于86%；各樣品的GC 含量大約45%，與一般物種中的GC 含量相近.上述結(jié)果表明，秋茄根系的轉(zhuǎn)錄組測序質(zhì)量較好.應(yīng)用Trinity 組裝平臺，共獲得112 243 條Transcripts 和61 970條Unigenes，而且Transcripts 與Unigenes 的N50 分別為2 361 和1 510 bp，這表明秋茄根系轉(zhuǎn)錄組片段組裝完整性較高.

表2 組裝結(jié)果統(tǒng)計

續(xù)表2

采用Bowtie 軟件將測序得到的reads 與Unigene 庫進行比對，總比對效率大于85%，并且多個位置比對read 所占百分比在41%左右，證明測序數(shù)據(jù)的總體比對效果較好，可以用于各轉(zhuǎn)錄本表達量的確定.以比對到參考基因序列上的reads 在基因上的分布情況進行測序隨機性評估表明，本研究中樣品reads 在參考基因上對應(yīng)的分布曲線比較平緩（圖2），表明樣品的測序隨機性良好.依據(jù)測序堿基數(shù)量和發(fā)現(xiàn)基因數(shù)進行測序飽和度評價，結(jié)果顯示（圖3），本研究中總RNA 在當前的測序量條件下已達到飽和，表明測序結(jié)果可以用于全轉(zhuǎn)錄組分析.

圖2 cDNA 片段的隨機性檢驗圖

圖3 測序飽和度分析圖

2.3 Unigenes 基因功能注釋

經(jīng)詳細比對分析，本研究最終成功獲得26 804 個秋茄根系轉(zhuǎn)錄組Unigenes 序列的注釋信息.基因注釋的統(tǒng)計結(jié)果見表3.注釋結(jié)果顯示組裝成功的61 970 條Unigenes 序列有26 668 條可以在Nr 數(shù)據(jù)庫中獲得注釋，占總數(shù)的43.03%，17 761 條Unigenes 序列與Swissprot 數(shù)據(jù)庫中獲得注釋，占28.66%；22 143 條Unigenes 序列在GO 庫中獲得注釋，占35.73%；8 324 條在COG 數(shù)據(jù)庫中得到注釋，占13.43%；6 258 條可以注釋到KEGG 庫中，占10.10%.獲得成功注釋的基因為深入研究秋茄對鹽脅迫的應(yīng)答提供了參考依據(jù)，尚有35 166 條Unigenes 序列無法從現(xiàn)有數(shù)據(jù)庫中獲得注釋，獲得注釋的基因為將來進一步挖掘秋茄中的新基因提供了有價值的信息.

表3 Unigenes 序列注釋結(jié)果統(tǒng)計表

2.4 Unigenes 基因的COG 分類

為了進一步分析秋茄根系NaCl 處理下轉(zhuǎn)錄組序列的功能，本研究將COG 注釋到的8 324 條Unigenes 基因進行了功能分類，共獲得11 719 個功能注釋信息，可劃分為25 個功能類群，結(jié)果見圖4.

圖4 秋茄根系NaCl 處理下Unigenes 基因的COG 功能注釋分布圖

COG分類結(jié)果顯示，包含基因數(shù)量排在前10 位的類別分別為，一般功能預(yù)測基因（ceneral function prediction），達2 194 個，占26.36%；復(fù)制、重組與修復(fù)基因（replication,recombination and repair）1 022個，占12.28%；與轉(zhuǎn)錄（transcription）相關(guān)的功能基因1 008 個，占12.11%；與翻譯、核糖體結(jié)構(gòu)與生物合成（translation,ribosomal structure and biogenesis）相關(guān)的基因948 個，占11.39%；與信號轉(zhuǎn)導(dǎo)機制（signal transduction mechanisms）相關(guān)基因887 個，占10.66%；翻譯后修飾、蛋白轉(zhuǎn)運和分子伴侶（posttranslational modification,protein turnover,chaperones）基因872 個，占10.48%；碳水化合物運輸與代謝（carbohydrate transport and metabolism）基因654 個，占7.86%；氨基酸轉(zhuǎn)運和代謝（amino acid transport and metabolism）相關(guān)基因588 個，占7.06%；能量產(chǎn)生和轉(zhuǎn)換（energy production and conversion）相關(guān)基因467 個，占5.61%；無機離子轉(zhuǎn)運和代謝（inorganic ion transport and metabolism）相關(guān)基因406 個，占4.88%.從COG 分類結(jié)果來看，排在前10 位的基因類群大多數(shù)和植物應(yīng)答生物與非生物脅迫相關(guān)，這表明本研究所建立的秋茄NaCl 處理下的轉(zhuǎn)錄組庫包含了秋茄抵御逆境的相關(guān)基因.

2.5 Unigenes 基因GO 分類

對GO 注釋的2 2143 個Unigenes 基因進行功能顯著性富集分析，結(jié)果（圖5）表明，本研究獲得的秋茄根系Unigenes 基因涉及細胞組分（cellular component,CC）、分子功能（molecular function,MF）和生物學(xué)過程（biological process，BP）3 個大類中56 個小類.3 個大類中生物學(xué)過程獲得的注釋最多，占51.53%，包含：細胞過程（cellular process）、代謝過程（metabolic process）、刺激響應(yīng)（response to stimulus）、生物學(xué)調(diào)節(jié)（biological regulation）、發(fā)育過程（developmental process）和細胞組分組織與生物合成（cellular component organization or biogenesis）等24 個小類.細胞組分類獲得注釋占35.59%，包含：細胞部分（cell part）、細 胞（cell）、細 胞 器（organelle）、膜（membrane）、細 胞 器 部 分（organelle part）、大 分 子 復(fù) 合 體（macromolecular complex membrane part）等16 個小類.分子功能類獲得注釋較少，占12.88%，包含：核酸結(jié)合轉(zhuǎn)錄因子活性（nucleic acid binding transcription factor activity）、結(jié)合（binding）、催化活性（catalytic activity）、轉(zhuǎn)運體活性（transporter activity）、結(jié)構(gòu)分子活性（structural molecule activity）、分子傳感器活性（molecular transducer activity）和酶調(diào)節(jié)活性（enzyme regulator activity）等16 個小類.上述結(jié)果顯示，秋茄鹽脅迫下轉(zhuǎn)錄組中涉及生物學(xué)過程基因較多，表明秋茄應(yīng)答鹽脅迫的分子調(diào)控與細胞生物學(xué)過程關(guān)系密切.

圖5 秋茄根系NaCl 處理下Unigenes 基因的GO 分類圖

2.6 Unigenes 基因的代謝通路分析

對基因的KEGG pathway 分析可以幫助我們深入了解基因的生物功能及相互作用關(guān)系.對不同濃度NaCl 脅迫下秋茄轉(zhuǎn)錄組庫中Unigenes 基因的KEGG pathway 富集分析顯示，獲得注釋的6 258 條Unigenes 基因參與了116 條代謝途徑，按照獲得注釋的基因數(shù)多少，排在前10 位的通路如表4 所示.

表4 Unigenes 基因的KEGG pathway 富集結(jié)果統(tǒng)計

Unigenes 基因通路富集結(jié)果顯示核糖體、內(nèi)質(zhì)網(wǎng)蛋白加工、植物信號轉(zhuǎn)導(dǎo)、氧化磷酸化、RNA 轉(zhuǎn)運剪接、核酸代謝以及糖代謝等途徑為包含較多基因的通路，而這些通路已被證明與植物逆境應(yīng)答密切相關(guān).這表明，本研究建立的秋茄轉(zhuǎn)錄組庫包含了多個與抗逆相關(guān)的代謝途徑，可以作為后續(xù)秋茄耐鹽相關(guān)基因篩選的參考數(shù)據(jù)庫.

2.7 Unigenes 基因開放閱讀框預(yù)測

對Unigenes 基因開放閱讀框（ORF）進行預(yù)測，共獲得61 690 條序列，平均長度393.63 bp，其中長度大于2 000 bp 的有1 703 條，N50 為984 bp，與其他物種的預(yù)測結(jié)果基本一致.由于秋茄缺乏全基因組測序信息，轉(zhuǎn)錄組庫中尚有大量序列無法從現(xiàn)有數(shù)據(jù)庫中得到注釋，而這些序列可能為新的蛋白編碼序列.

2.8 耐鹽相關(guān)基因的qRT-PCR 分析

依據(jù)本研究獲得的相關(guān)基因數(shù)據(jù)庫，對前人在紅樹植物中發(fā)現(xiàn)的鹽脅迫相關(guān)的6 個基因在不同NaCl 濃度下的表達差異進行了qRT-PCR 分析，結(jié)果（表5）表明，上述6 個基因在不同NaCl 濃度處理下差異表達，表明上述基因可能參與了秋茄對高鹽環(huán)境的響應(yīng)，也表明本研究構(gòu)建的數(shù)據(jù)庫序列可以用于秋茄耐鹽基因差異表達研究.

表5 本研究檢測到的部分與植物耐鹽相關(guān)的基因qRT-PCR 分析結(jié)果

3 結(jié)論與討論

測序深度是保證轉(zhuǎn)錄組庫質(zhì)量的關(guān)鍵因素.當測序量達到一定標準時，檢測到的表達基因數(shù)量就會逐漸趨近飽和.本研究共獲得33 246 108 條的reads，總堿基數(shù)達到了6.71 G，達到飽和測序量的200倍，表明測序深度很高，利用此數(shù)據(jù)構(gòu)建的轉(zhuǎn)錄組庫結(jié)果可靠，能保證低豐度基因的檢測.

選擇合適的組裝軟件是獲得高質(zhì)量轉(zhuǎn)錄組序列的關(guān)鍵.常見的幾種組裝軟件中Trinity 是進行無參考基因組非模式生物轉(zhuǎn)錄組從頭組裝的最佳軟件[9].本研究利用Trinity 組裝共得到112 243 條Transcripts 和61 970 條Unigenes，Transcripts 與Unigenes 的N50 分別為2 361 bp 和1 510 bp，有13 642 條Unigenes序列的長度大于1 000 bp，平均長度為790.75 bp，這高于前人獲得的同為紅樹植物且無參考基因組的海桑（Sonneratia caseolaris）（589 bp）[10]，表明本研究獲得秋茄轉(zhuǎn)錄組序列總數(shù)、平均長度、以及片段組裝完整性都比較好.

對測序結(jié)果的準確注釋是開展基因功能及代謝通路分析的前提.本研究借助Nr、Swissprot、GO、COG和KEGG 數(shù)據(jù)庫對秋茄根系Unigenes 基因進行了注釋，共獲得了43.25%的注釋信息，尚有56.75%的序列沒有獲得注釋信息，這主要原因可能是因為組裝出的短序列難以達到具有統(tǒng)計意義的匹配閥值，無法進行同源性比對，也可能是因為未比對上的序列可能為秋茄特有基因或非編碼序列，這需要后續(xù)工作進一步挖掘.

GO 富集分析、COG 功能分類和KEGG pathway 三者相互結(jié)合成為從新物種中挖掘功能基因的重要手段[11].本研究通過對Unigenes 基因的GO 富集、COG 功能分類和KEGG 通路分析發(fā)現(xiàn)，秋茄在鹽脅迫下啟動了大量與逆境脅迫相關(guān)的基因表達，這些基因參與了多種代謝調(diào)控途徑，這預(yù)示著秋茄響應(yīng)鹽脅迫的過程相當復(fù)雜，涉及到多條代謝通路交叉調(diào)控，需要多個基因組成調(diào)控網(wǎng)絡(luò)來調(diào)控.

綜上所述，通過對不同NaCl 濃度處理下秋茄根系轉(zhuǎn)錄組庫測序，共得到112 243 條Transcripts 和61 970 條Unigenes，Transcripts 與Unigenes 的N50 分別為2 361 bp 和1 510 bp，表明本研究建立的轉(zhuǎn)錄組庫片段的組裝完整性.利用Nr、SwissProt、GO、COG 和KEGG 數(shù)據(jù)庫對Unigenes 基因進行注釋，有26 804 個Unigenes 基因獲得注釋信息.對Unigenes 基因的生物信息分析表明，這些抗逆相關(guān)基因是秋茄能夠適應(yīng)濱海潮間帶高鹽環(huán)境的分子基礎(chǔ)，對這些基因功能及所參與的代謝路徑的深入研究將為木本植物耐鹽分子機理的解析及未來作物抗鹽基因工程育種奠定理論基礎(chǔ).

[責(zé)任編輯 楊玉玲]