NaCl脅迫下黑果枸杞轉錄組測序分析

馬彥軍 段慧榮 魏佳 Richard John Tiika 單立山 馬瑞

（1. 甘肅農業大學林學院，蘭州 730070；2. 中國農業科學院蘭州畜牧與獸藥研究所，蘭州 730050）

土壤鹽漬化是全球性的生態問題，不僅導致生態環境惡化，而且也是植物生長的主要非生物脅迫因素［1］。我國鹽漬化危害嚴重，各類鹽堿地面積約9 913萬hm2［2］。如何高效利用這些鹽堿化土地已成為一個重要的課題和研究方向［3］。利用鹽生植物進行鹽堿地生物改良是鹽堿地恢復和改良的重要途徑，而選擇適宜的種植材料是實施鹽堿地生物恢復和改良的關鍵。目前，適宜于鹽堿地改良的植物種類比較有限，不同種源的同一植物種耐鹽性差異也較大，因此，有必要通過選育和培育抗鹽品種來提高植物的抗鹽性［4］。為此，在分子生物學層面上闡明植物的耐鹽機制，獲得植物的耐鹽相關基因，通過現代生物技術實現耐鹽基因轉化，并以此提高植物耐鹽性已成為目前耐鹽植物研究的熱點［5-6］。

黑果枸杞（Lycium ruthenicumMurr.）為茄科（Solanceae）枸杞屬（LyciumL.）多年生灌木，分布于新疆、青海、甘肅、寧夏和內蒙古等地［7］。常以灌叢狀生于鹽堿荒地、鹽化沙地、鹽湖岸邊、道路兩側集水區等各種鹽漬化土壤或荒漠環境，是我國荒漠區特有的耐鹽抗旱同時具有很高經濟價值及營養價值的野生植物［8-9］。黑果枸杞具有抗旱、抗寒、耐鹽堿、根蘗性強、耐土壤貧瘠及營養價值高等諸多優點，是鹽堿地治理的先鋒樹種，具有廣泛的經濟前景和開發潛力［10-11］。目前有關黑果枸杞抗鹽性研究主要集中在生理生化和解剖結構等方面［12-14］，對其分子水平的研究報道則相對較少［15-17］。隨著高通量測序技術的快速發展，轉錄組測序成為研究基因差異表達的重要手段［18-20］。本研究利用Illumina測序平臺，進行不同濃度鹽脅迫處理下的黑果枸杞根和葉的轉錄組測序，篩選鹽脅迫下葉和根中的差異表達基因，旨為揭示黑果枸杞響應鹽脅迫的分子機制提供基礎資料。

1 材料與方法

1.1 材料

以甘肅農業大學林學院實驗室保存的黑果枸杞組培苗為材料。

1.2 方法

1.2.1 NaCl 處理 2018年10月，選取生長健壯、長勢一致的黑果枸杞組培苗，開蓋煉苗5 d后，清洗根部殘留的培養基，轉移到塑料盒（口徑10 cm，高10 cm，帶有通氣裝置）。塑料盒中裝有等量含有不同濃度NaCl的1/2 Hoagland 營養液進行組培苗鹽脅迫處理，每個處理1株，3次重復。1/2 Hoagland營養液中添加的NaCl 濃度為：0（CK）mmol/L、50 mmol/L和250 mmol/L。上述處理分別在脅迫0、1和12 h時取根和葉（含嫩莖）各0.1 g用于轉錄組和基因表達譜的分析。

1.2.2 總RNA的提取、文庫構建和轉錄組測序 用mirVanaTMmiRNA ISOlation Kit（Ambion-1561）試劑盒提取和純化黑果枸杞葉片和根的總RNA。分別用Nanodrop 2000與Agilent2100 Bioanalyzer 進行總RNA質量和純度檢測。分別取NaCl脅迫處理組和對照組的黑果枸杞根和葉樣品的RNA，送樣至上海歐易生物醫學科技有限公司進行轉錄組測序。

1.2.3 Unigene功能注釋和差異表達基因分析 使用 blastx［21］，取閾值 e＜1e-5，將 Unigene 分別注釋到NR、KOG和Swissprot數據庫?；赟wissprot結果，將Swissprot ID映射到GO term，獲得Unigene的GO注釋，最后將 Unigene比對到KEGG［22］數據庫獲得通路信息。Unigene的FPKM［21］和count使用 bowtie2［22］和 eXpress［23］軟件分析得到。通過eXpress軟件獲取落到各個樣本中Unigene 的 reads 數目，使用 DESeq［23］R package 的 estimateSizeFactors函數對數據進行標準化，并使用nbinomTest函數計算差異比較的 pvalue和foldchange值。挑選出P＜0.05，差異倍數大于2的差異Unigene，并進行差異 Unigene的 GO 和 KEGG［24］富集分析，以判定差異Unigene主要影響的生物學功能或通路。利用HMMER［24］軟件比對 Pfam 數據庫來進行 Unigene 的功能分析。

2 結果

2.1 總RNA 質量檢測和組裝結果分析

經檢測，樣品總RNA的A260/A280均大于1.8，28S/18S≥0.7，RIN≥7，表明所提RNA可滿足后繼實驗分析需求。

黑果枸杞幼苗在不同濃度NaCl脅迫和不同取樣時間下葉和根共30個樣本，通過轉錄組測序共獲得 222.49 G的原始數據，各樣本的 Q30 堿基比均在 95%以上，GC含量均在43%以上（表1）。說明測序結果準確度較高，可用于后續分析。利用Trinity 軟件對測序所得數據進行合并組裝，共獲得86 037個Unigene，其中長度在 1 kb 以上的Unigene有 31 070 條，這些Unigene可作為后續實驗研究的重點對象；Unigene 平均長度為 1 097.27 bp（表2）。

表1 有效數據評估統計

表2 組裝結果統計分析

表3 單基因序列注釋統計表

2.2 基因功能注釋

由表3可知，鹽脅迫下黑果枸杞根和葉的86 037個Unigene中有47 108個Unigene在不同數據庫中的得到了注釋，占總Unigened的54.76%，還有38 929個Unigene在這些數據庫中沒有得到注釋。其中注釋到NR數據庫的Unigene 最多，達到46 594個，占54.16%；注釋到Pfm數據庫的Unigene 最少，僅有123個，占0.14%；共有16個Unigene在所有數據庫中都得以注釋。

由表4可知，黑果枸杞轉錄組所得Unigene在NR數據庫注釋中與馬鈴薯（Solanum tuberosum）同源序列最多，為8 594個，占注釋Unigene 的18.48%，番茄（Solanum lycopersicum）同原序列最少，為2 223個，占注釋Unigened 4.78%，其他物種的5 555，占注釋Unigene 的11.95%。

圖4 黑果枸杞轉錄組Unigenes 的 NR注釋物種分布表

2.3 GO 注釋和分類

對獲得相應GO條目的Unigenes進行統計分析（表5），共有28 727個Unigenes獲得235 271 個GO注釋，平均每條8.19個。獲得的GO 數據庫注釋的Unigene 可分為細胞組分（Cellular component）、分子功能（Molecular function）和生物過程（Biological process）3 個GO 類別的 51個小組。細胞組分涉及 104 407個GO條目，分為 14 個功能組，其中細胞（Cell，23 813個）、細胞部分（Cell part，23 778個）、細胞器（Organelle，18 975）和細胞膜（Membrane，9 479個）涉及的Unigene 較多；分子功能涉及 37 337個GO 條目，分為 15個功能組，其中催化活性（Catalytic activity，14 615個）和結合（Binding，17 383 個）含 Unigene 較多；93 527個GO 條目歸屬于生物過程，分為 22個功能組，其中代謝過程（Metabolic process，16 042 個）、細胞進程（Cellular process，19 213個）和刺激響應（Response to stimulus，8 326個）涉及的 Unigene 較多。

2.4 KOG注釋和分類

為了進一步分析黑果枸杞根和葉轉錄組 Unigene的功能，進行了 KOG 功能分類分析，共有25 246個Unigenes獲得 27 912個KOG 注釋，平均每條1.1個。分類結果如表6所示，共獲得25個不同的功能分類。一般功能預測（General function prediction only）的Unigene為7 022條，是最大的功能類群；其次是翻譯后修飾、蛋白質翻轉和分子伴侶（Post translation almodification，protein turnover，chaperones）2 755 條，信號轉導類機制（Signal transduction mechanisms）次之，有2 349 條，最少的是細胞活性（Cell motility）功能類別，僅12條。

2.5 KEGG代謝通路分析

注釋到KEGG的15 699個Unigenes獲得 15 844個KEGG 注釋。參與的代謝通路可歸為4大類別、24個子類。由表7可知，4種代謝通路大類中，與代謝（Metabolism）相關的通路獲得8 143個注釋，遺傳信息處理（Genetic information processing）3 893個，細胞過程（Cellular processes）與環境信息處理（Environmental information processing）分別是1 955個和1 853個。進一步細分為24個子類代謝途徑，其中，翻譯（Translation）獲得注釋最多，為1 983個；其次為信號傳導（Signal transduction）和碳水化合物代謝（Carbohydrate metabolism），分別為1 777個和1 512個。

共有21 215個 Unigenes 歸入 211條代謝通路，靠前的有 14 個代謝 通 路（表 8）。核糖體（Ribosome）代謝途徑注釋到的Unigenes數量最多，有1 241個Unigenes；其次為碳代謝（Carbon metabolism）723、氨基酸的生物合成（Biosynthesis of amino acids）546和內質網蛋白處理（Protein processing in endoplasmic reticulum）490 等途徑。

2.6 Unigene 差異表達分析

由表9可以看出，黑果枸杞在不同濃度NaCl和不同處理時間下與CK相比，葉片和根隨著NaCl濃度增大和處理時間延長，上調基因和下調基因數呈增加趨勢；葉片上調基因數（7 514）小于下調基因數（9 032），而根的上調基因數（12 347）大于下調基因數（11 559）。

2.7 SSR信息分析

利用軟件 MISA 對黑果枸杞鹽脅迫下轉錄組測序所獲得的Unigene進行 SSR 預測，結果見表10。由表10可知，共有28 325個SSR位點，單核苷酸SSR最多，為19 960，占70.47%；5核苷酸SSR最少，為0.11%。重復單元重復出現的次數大于11次以上最多，為9 257，占32.68%；重復單元重復出現7次的最少，為486，占1.72%。

表5 黑果枸杞轉錄組Unigene GO功能分類表

圖6 黑果枸杞轉錄組Unigenes的KOG功能分布表

3 討論

對于未完成全基因組測序的物種來說，采用轉錄組測序技術可獲得大量的轉錄本信息，這可為植物基因表達的綜合分析提供合理可靠的數據資源［25］。因此，對黑果枸杞在鹽脅下的轉錄組信息進行系統分析，為全面了解黑果枸杞抗鹽分子機制和挖掘新的抗鹽基因提供基礎資料。本研究在沒有參考基因組測序情況下，對不同濃度NaCl溶液脅迫不同時間時的黑果枸杞組培苗的葉片和根進行轉錄組測序，共獲得86 037條Unigenes。將這些Unigene在 NR、Swissprot、KEGG、KOG、eggNOG、GO 及Pfam數據庫中進行注釋，共有47 108個Unigene得到了注釋，占總Unigened的54.76%，還有38 929個Unigene沒有得到注釋。這一結果在目前已測序的許多植物當中都存在［26］，即并不是所有組裝得到的Unigene都能在已知的有關基因數據庫得到注釋。這是由于大量非模式植物的基因組測序工作沒有進行，缺少相關基礎研究數據，需要對這些沒有注釋的基因進一步分析研究，來確定這些基因在植物生長發育過程中的作用，從而豐富基因數據庫。

植物在受到鹽脅迫時會產生大量轉導信號，有離子信號轉導、滲透信號轉導、解毒類信號轉導等；同時植物體內蔗糖和淀粉的分解代謝過程加強，分解的糖類不僅為植物生長發育提供能量，而且分解產生的小分子單糖提高了細胞中可溶性糖的含量和細胞的滲透勢，以抵抗高濃度鹽離子帶來的滲透脅迫［27］。本研究在KEGG代謝途徑注釋中翻譯（Translation）獲得注釋最多，其次為信號傳導（Signal transduction）和碳水化合物代謝（Carbohydrate metabolism）。這是因為當黑果枸杞受到鹽脅迫時，自身的碳水化合物進行代謝來產生足夠的能量使植物能夠進行各種生理反應，植物需要合成大量的酶來催化體內的各種生理生化反應，同時將受到的刺激信號進行傳導，以便植物能產生應激反應來應對不良環境造成的影響。這一結果與濱柃［28］、遼寧堿蓬［29］、海濱錦葵［30］等植物在鹽脅迫下轉錄組分析結果基本一致，這也說明不同植物在鹽脅迫有相似的抗鹽機制。

表7 黑果枸杞轉錄組Unigenes的KEGG功能注釋分布統計表

表8 黑果枸杞轉錄組Unigenes Unigenes數量最多的14個代謝通路

在DEGs分析中，葉的上調基因數（7 514）小于下調基因數（9 032），而根的上調基因數（12 347）大于下調基因數（11 559）。這是由于在鹽脅迫下，植物根系最早感受逆境脅迫信號，并產生相應的生理反應，繼而影響地上部生長［31］。張倩倩［32］在對中國石竹研究過程也發現，在鹽脅迫過程中根中的DEGs比葉中高。這也說明當黑果枸杞受到鹽脅迫時根部就可能會誘導更多與抗鹽有關基因表達，使得黑果枸杞能夠在鹽脅迫下正常生長，而葉片中可能與抗鹽關系不大的基因被抑制，使得黑果枸杞與抗鹽關系不明顯的性狀得以抑制，把更多能量用于與抗鹽有關的性狀生長。需要進一步的對這些差異基因進行分析，篩選出與黑果枸杞抗鹽有關的基因，從分子生物學角度來解釋黑果枸杞抗鹽機制。

表9 黑果枸杞轉錄組差異表達基因分析

表10 EST-SSR的類型、數量及分布頻率

SSR是以特異引物PCR為基礎的分子標記技術［33］，因其多態性高、共顯性遺傳、覆蓋面廣等優點［34］，廣泛運用于遺傳圖譜構建、目標基因的標定、遺傳多樣性分析和分子標記輔助育種等方面［35］。已對許多物種進行了基于轉錄組SSR信息的分析研究，如馬尾松（Pinus massoniana）［36］、藍靛果忍冬（Lonicera caerulea）［37］、山桐子（Idesia polycarpa）［38］等。這些研究結果表明轉錄組中具有豐富的SSR位點信息，單核苷酸、二核苷酸和三核苷酸為優勢重復類型在不同物種中不同。本研究利用MISA 對黑果枸杞鹽脅迫下轉錄組測序所獲得的Unigene進行 SSR預測，共發現28 325個SSR位點，優勢重復基序為單核苷酸，占總 SSR 位點的70.47%，這一研究結果與尹躍等［39］的結果一致，即黑果枸杞轉錄組的 SSR主要類型為單核苷酸重復基序（74.33%）。

4 結論

本研究利用Illumina HiSeq X Ten 測序儀對不同濃度NaCl溶液脅迫下的黑果枸杞轉錄組進行測序，轉錄組測序共產生 222.49 Gb 原始數據，拼接出 Unigenes 86 037條，注釋到 7大功能數據庫（GO、KEGG、KOG、NR、Pfam、Swiss-Prot和 egg NOG）上的 Unigenes 總數為46 594個，還有38 929個Unigenes在這些數據庫中沒有得到注釋。通過 GO分類和 KEGG Pathway 富集性分析，將分別歸于 51個 GO 類別和 211條代謝途徑。差異表達基因分析顯示，黑果枸杞葉片葉片中的上調基因數（7 514）小于下調基因數（9 032），而根中的上調基因數（12 347）大于下調基因數（11 559）。在黑果枸杞鹽脅迫下轉錄組中發現28 325個SSR位點，最多的為單核苷酸 SSR，占 70.47%，重復單元重復出現的次數大于11次以上最多，占32.68%。