NaCl脅迫下黑果枸杞轉(zhuǎn)錄組測(cè)序分析

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

（1. 甘肅農(nóng)業(yè)大學(xué)林學(xué)院，蘭州 730070；2. 中國(guó)農(nóng)業(yè)科學(xué)院蘭州畜牧與獸藥研究所，蘭州 730050）

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

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

1 材料與方法

1.1 材料

以甘肅農(nóng)業(yè)大學(xué)林學(xué)院實(shí)驗(yàn)室保存的黑果枸杞組培苗為材料。

1.2 方法

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

1.2.2 總RNA的提取、文庫(kù)構(gòu)建和轉(zhuǎn)錄組測(cè)序 用mirVanaTMmiRNA ISOlation Kit（Ambion-1561）試劑盒提取和純化黑果枸杞葉片和根的總RNA。分別用Nanodrop 2000與Agilent2100 Bioanalyzer 進(jìn)行總RNA質(zhì)量和純度檢測(cè)。分別取NaCl脅迫處理組和對(duì)照組的黑果枸杞根和葉樣品的RNA，送樣至上海歐易生物醫(yī)學(xué)科技有限公司進(jìn)行轉(zhuǎn)錄組測(cè)序。

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

2 結(jié)果

2.1 總RNA 質(zhì)量檢測(cè)和組裝結(jié)果分析

經(jīng)檢測(cè)，樣品總RNA的A260/A280均大于1.8，28S/18S≥0.7，RIN≥7，表明所提RNA可滿足后繼實(shí)驗(yàn)分析需求。

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

表1 有效數(shù)據(jù)評(píng)估統(tǒng)計(jì)

表2 組裝結(jié)果統(tǒng)計(jì)分析

表3 單基因序列注釋統(tǒng)計(jì)表

2.2 基因功能注釋

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

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

圖4 黑果枸杞轉(zhuǎn)錄組Unigenes 的 NR注釋物種分布表

2.3 GO 注釋和分類

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

2.4 KOG注釋和分類

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

2.5 KEGG代謝通路分析

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

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

2.6 Unigene 差異表達(dá)分析

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

2.7 SSR信息分析

利用軟件 MISA 對(duì)黑果枸杞鹽脅迫下轉(zhuǎn)錄組測(cè)序所獲得的Unigene進(jìn)行 SSR 預(yù)測(cè)，結(jié)果見(jiàn)表10。由表10可知，共有28 325個(gè)SSR位點(diǎn)，單核苷酸SSR最多，為19 960，占70.47%；5核苷酸SSR最少，為0.11%。重復(fù)單元重復(fù)出現(xiàn)的次數(shù)大于11次以上最多，為9 257，占32.68%；重復(fù)單元重復(fù)出現(xiàn)7次的最少，為486，占1.72%。

表5 黑果枸杞轉(zhuǎn)錄組Unigene GO功能分類表

圖6 黑果枸杞轉(zhuǎn)錄組Unigenes的KOG功能分布表

3 討論

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

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

表7 黑果枸杞轉(zhuǎn)錄組Unigenes的KEGG功能注釋分布統(tǒng)計(jì)表

表8 黑果枸杞轉(zhuǎn)錄組Unigenes Unigenes數(shù)量最多的14個(gè)代謝通路

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

表9 黑果枸杞轉(zhuǎn)錄組差異表達(dá)基因分析

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

SSR是以特異引物PCR為基礎(chǔ)的分子標(biāo)記技術(shù)［33］，因其多態(tài)性高、共顯性遺傳、覆蓋面廣等優(yōu)點(diǎn)［34］，廣泛運(yùn)用于遺傳圖譜構(gòu)建、目標(biāo)基因的標(biāo)定、遺傳多樣性分析和分子標(biāo)記輔助育種等方面［35］。已對(duì)許多物種進(jìn)行了基于轉(zhuǎn)錄組SSR信息的分析研究，如馬尾松（Pinus massoniana）［36］、藍(lán)靛果忍冬（Lonicera caerulea）［37］、山桐子（Idesia polycarpa）［38］等。這些研究結(jié)果表明轉(zhuǎn)錄組中具有豐富的SSR位點(diǎn)信息，單核苷酸、二核苷酸和三核苷酸為優(yōu)勢(shì)重復(fù)類型在不同物種中不同。本研究利用MISA 對(duì)黑果枸杞鹽脅迫下轉(zhuǎn)錄組測(cè)序所獲得的Unigene進(jìn)行 SSR預(yù)測(cè)，共發(fā)現(xiàn)28 325個(gè)SSR位點(diǎn)，優(yōu)勢(shì)重復(fù)基序?yàn)閱魏塑账幔伎?SSR 位點(diǎn)的70.47%，這一研究結(jié)果與尹躍等［39］的結(jié)果一致，即黑果枸杞轉(zhuǎn)錄組的 SSR主要類型為單核苷酸重復(fù)基序（74.33%）。

4 結(jié)論

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