祁連山黃參葉片轉錄組測序及生物信息學分析

張春梅，閆 芳，宋 海，張喜峰1,，陳 葉1,

(1. 河西學院 農(nóng)業(yè)與生態(tài)工程學院，甘肅 張掖 734000；2. 甘肅省河西走廊特色資源利用重點實驗室，甘肅 張掖 734000;3. 河西學院 生態(tài)與綠洲農(nóng)業(yè)研究院，甘肅 張掖 734000)

黃參Sphallerocarpusracills，傘形科Apiaceae迷果芹屬Sphallerocarpus的單種植物[1]，在我國零星分布于西北、東北、華北地區(qū)，尤其分布于祁連山、焉支山(甘肅張掖市境內(nèi))，像張掖市山丹縣這樣大面積形成群落優(yōu)勢的情況很罕見[2]。研究黃參對祁連山區(qū)水土保持、遏制草場沙漠化和維持生物多樣性具有重要的生態(tài)和經(jīng)濟意義?！侗静菥V目》中記載，黃參具有補氣養(yǎng)血、滋補肝腎、通經(jīng)活絡等功效，其肉質根富含人體必需的16種氨基酸，被譽為“小人參”[3]，作為營養(yǎng)豐富、經(jīng)濟、藥用價值極高的天然珍品，黃參成為甘肅特產(chǎn)、西部開發(fā)交易會指定產(chǎn)品，在食品、醫(yī)藥、化工等領域具有極大的應用潛力與開發(fā)前景。由于生態(tài)環(huán)境惡化及當?shù)匕傩盏穆訆Z式經(jīng)營，野生黃參種群數(shù)量驟減[3]。

絕大部分非模式生物缺乏基因組數(shù)據(jù)，因此，獲得轉錄組學信息尤為重要[4-7]。轉錄組學是獲取基因序列的首選方式[8-9]，非常有利于研究無參考基因組的非模式植物[10-11]，研究轉錄組學有助于揭示生物體的基因表達、研究結構變異、新基因以及植物的優(yōu)良性狀及功能基因的定位[12-15]。

目前，瀕危野生黃參種質的功能基因組、基因序列信息及遺傳背景匱乏，對黃參的搶救性保護及其遺傳結構研究迫切需要基因組資源。關于黃參轉錄組信息尚未見報道，本研究采用BGISEQ-500平臺，對黃參幼嫩葉片進行轉錄組測序，研究黃參轉錄組，并結合生物信息學對獲得的Unigene(轉錄組數(shù)據(jù)庫序列)開展功能注釋、代謝通路和EST-SSR分析，獲得的轉錄組信息將為今后黃參分子標記的開發(fā)和關鍵基因的克隆以及功能分析等提供科學數(shù)據(jù)，為黃參基因組水平的研究奠定基礎。

1 材料和方法

1.1 實驗材料

材料采自甘肅省張掖市山丹縣軍馬場。野生黃參，于2018年5月進行單株取樣，采集當年剛生長出的幼嫩葉片，迅速放入液氮速凍后，保存于-80 ℃冰箱中。

1.2 主要試劑及儀器

1.2.1 主要儀器

超級工作臺(蘇州凈化工作臺設備有限公司)；多樣品組織研磨儀(上海凈信實業(yè)發(fā)展有限公司)；SANYO制冰機(濟南金茂科創(chuàng)科技有限公司)；冷凍離心機(CL-21R，Thermo，美國)；超聲破碎儀(FS-150，Ultrasonic processor，中國)；震蕩儀(Shaker，Thermo，美國)；金屬浴(GL-150，其林貝爾，中國)；分析天平(BSA224S，Sartorius，德國)；Bio-Rad凝膠成像系統(tǒng)(麥克森公司)；NanoDrop 2000微量分光亮度計(Thermo Scientific, USA)；Agilent 2100 Bioanalyzer(美國安捷倫公司)；NovaseqTM6000(美國Illumina公司)。

1.2.2 主要試劑

TRIzol?Reagent購自Invitrogen公司；氯仿、異丙醇、乙醇等化學試劑購自天津市大茂化學試劑廠，分析純；NovaSeqTM6000 v1.5試劑盒(美國Illumina公司)。

1.3 實驗方法及步驟

1.3.1 黃參葉片RNA 提取及轉錄組測序

采集當年剛生長出的幼嫩葉片，液氮速凍并保存于-80 ℃冰箱中，以保證RNA的相態(tài)保持在剛采摘時的狀態(tài)。采用試劑盒法，Trizol試劑，TIANGEN提取方法參照試劑說明。葉肉細胞移入1.5 mL 離心管中，加入1 mL Trizol，混勻，室溫靜置5 min；加入0.2 mL氯仿，振蕩15 s， 靜置2 min；4 ℃12 000 g離心15 min，裂解液分層成水相和有機相，取上清；加入0.5 mL異丙醇，將管中液體輕輕混勻，室溫靜置10 min；4 ℃ 12 000 g 離心10 min，棄上清，水相轉移后，加入1 mL體積分數(shù)為75%乙醇，輕輕洗滌沉淀。4 ℃ 7 500 g離心5 min，短暫離心，吸去殘留液體，棄上清；晾干3～5 min；將獲得的RNA沉淀溶于30 μL的DEPC(焦炭酸二乙酯)水中。DEPC是RNA酶的強抑制劑，是一種潛在的致癌物質，操作中應在通風條件下進行，并避免接觸皮膚。

1.3.2 cDNA 文庫構建及轉錄組序列組裝

取一定量的RNA 樣品，使用 oligo(dT)磁珠富集mRNA。加入試劑盒中提供的打斷試劑，適溫反應一定時間后mRNA片段化。合成一鏈、二鏈cDNA。配制反應體系，使接頭與cDNA連接。PCR反應及產(chǎn)物回收、擴增。PCR產(chǎn)物變性，充分混勻，得到單鏈環(huán)形產(chǎn)物，PCR產(chǎn)物變性，即得到文庫。文庫質量委托陜西致研生物科技有限公司使用 Agilent 2100 Bioanalyzer 檢測，檢測合格后采用BGISEQ-500測序。利用Trinity軟件對reads進行序列組裝，使用BUSCO軟件對組裝序列進行質量評估，評估基因組裝完整性。

1.3.3 功能注釋及SSR檢測

利用生物信息學分析獲得的黃參Unigene。為獲得全面的基因功能信息，對組裝得到的Unigene進行7大功能數(shù)據(jù)庫注釋，包括NR(NCBI non redundant protein sequences)、NT(NCBI nucleotide sequences)、KOG/COG(clusters of orthologous groups of proteins/eu-Karyotic ortholog groups)、GO (gene ontology，基因本體)、KEGG(kyoto encyelopedia of genes and genomes)、SwissProt (reviewed protein sequence database)和Pfam(protein family)。使用MISA對Unigene進行檢測，軟件參數(shù)為1-12、2-6、3-5、4-5、5-4、6-4、100，150，其中：1-12代表單堿基重復至少12次才算SSR，100表示2個SSR之間的間隔堿基數(shù)大于100 bp，150表示SSR位點距離兩端側翼序列大于150 bp。雙堿基6次，三堿基5次，以此類推，重復單元最多有6個堿基，2個微衛(wèi)星之間的距離小于100 bp。隨機選取10對引物，委托陜西致研生物科技有限公司合成。

2 結果與分析

2.1 高通量測序、de novo組裝

測序共獲得70.24 M原始讀數(shù)，去除低質量的reads，得到68.48 M clean reads，最終獲得6.85 Gb。質量評估結果見表1，轉錄本和Unigene組裝統(tǒng)計結果見表2。由表可見，測序得到的黃參數(shù)據(jù)質量較高，可以滿足后續(xù)的生物信息學分析。其中，長度為200～300 bp的有25 343條Unigene(40.66%)，長度為>300～500 bp的有16 463條Unigene(26.42%)，長度為>500～3 000 bp的有15 276條Unigene(24.51%)，5 241條Unigene(10.01%)的長度大于3 000 bp，由此可知，隨著基因長度增加，基因數(shù)量下降。使用Transdecoder檢測出62 323個CDS(見表2)。

表1 黃參測序后的質量評估Tab. 1 Quality evaluation of-sequencing output data in Sphallerocarpus racills

表2 黃參測序后的轉錄本和Unigene組裝統(tǒng)計Tab. 2 Data assembly for transcript and unigene in transcriptome of Sphallerocarpus racills

2.2 黃參轉錄組基因總體注釋情況

對黃參轉錄組基因進行7大功能數(shù)據(jù)庫注釋，結果見表3。結果顯示，注釋成功的Unigene數(shù)目最多的是NR(e≤10-5)，有66 451條，占總Unigene的66.46%；NT(e≤10-5)有49 390條，占49.40%；Swissprot(e≤10-5)有48 281條，占48.29%；KEGG(e≤10-10)有51 479條，占51.49%； KOG(e≤10-3)有61 116條，占61.13%；Pfam(e≤0.01)有55 859條，占55.87%。GO(e≤10-3)數(shù)據(jù)庫注釋到的基因最少，有36 958條，占36.97%；比對結果顯示，在7大數(shù)據(jù)庫中均能成功注釋的Unigene有17 074條，占總Unigene的17.08%。

表3 黃參轉錄組基因注釋情況統(tǒng)計Tab. 3 Statistics of gene annotation of Sphallerocarpus racills transcriptome

2.2.1 NR功能注釋

NR庫注釋結果的物種分布統(tǒng)計結果見圖1。圖1顯示，匹配最多的物種是胡蘿卜Daucuscarotasubsp.sativus，注釋到的基因數(shù)量最多，共有45 239條，占比為68.08%，同源性最高；其余依次為大麥Hordeumvulgaresubsp.vulgare、向日葵Helianthusannuus、藍隱藻GuillardiathetaCCMP2712、輪藻Klebsormidiumnitens，分別有963、711、638、524條Unigene，占比分別為1.45%、1.07%、0.96%、0.85%，黃參與這4種植物的同源性均相對較低，同源序列都不足2%，剩下近30%分布于其他物種，共有18 333條Unigene(占27.59%)屬于其他序列。

圖1 根據(jù)NR注釋結果統(tǒng)計注釋不同物種的分布Fig. 1 Statistics of species distribution of Unigenes annotation in NR database for Sphallerocarpus racills

2.2.2 GO功能注釋分類

注釋到生物學過程大類(共15個亞類)的基因有18 317條，其中，占比最高的亞類為細胞過程(cellular process)，注釋數(shù)量為11 026，其次為生物調(diào)節(jié)過程(biological regulation)，注釋數(shù)量為3 710，有關細胞增殖、碳利用等的基因表達甚少。

注釋到細胞組分大類(共11個亞類)的基因有25 108條，其中，細胞(cell)和細胞膜組分(membrane part)占比最高，分別為11 062和10 341條，再次是細胞器部分(organelle part)，基因為3 512條，而定位于病毒核心、細胞器的基因幾乎未表達。

注釋到分子功能大類(共14個亞類)的基因為34 615條，其中，最具有代表性的是參與分子結合功能(binding)，注釋數(shù)量最多(18 119)，其次是與催化活性(catalytic activity)相關的基因，注釋數(shù)量為16 074，而與蛋白質標簽、分子載體活性、分子轉運活性、養(yǎng)分庫活動、毒素活性相關的基因幾乎未表達(見圖2)。

圖2 Unigene的GO功能分類Fig. 2 GO function classification of Unigene

2.2.3 代謝通路分析

KEGG注釋結果見圖3，注釋到51 479條Unigene，占總Unigene的51.49%。注釋成功的所有Unigene歸屬到5大類代謝通路中的20條通路。11條與代謝有關的通路，有15 162個Unigene，占29.45%：以全局和總覽圖代謝通路(global and overview maps)在所有代謝通路中所占比例最高(11 346個，22.04%)。其次為碳水化合物代謝通路(4 137個，8.04%)和脂類代謝通路(2 618個，5.09%)。與遺傳信息處理相關的代謝通路有4條，其中以翻譯過程(translation，4 341，8.43%)所占比例最高，其次是折疊、分類與降解(folding, sorting and degradation，3 851，7.49%)代謝通路。與環(huán)境信息處理相關的通路有2條——信號轉導(signal transduction)和膜運輸(membrane transport)，分別為2 489、446條。細胞過程和生物系統(tǒng)的相關通路最少，各有1條，占4.67%和3.30%。

圖3 黃參轉錄組Unigene的代謝通路功能分布統(tǒng)計Fig. 3 KEGG classification of Sphallerocarpus racills transcriptome

2.2.4 KOG注釋分類

KOG數(shù)據(jù)庫包含了7個完整基因組真核生物的直系同源家族蛋白質。黃參轉錄組有61 116條Unigeine獲得了注釋，占總Unigene的61.13%，分為26大類功能區(qū)，包括功能預測、碳水化合物運輸與代謝、翻譯后修飾、蛋白質運輸、分子加工、信號轉導代謝等蛋白質家族。在不同的功能分類中，基因數(shù)量存在明顯差異，一般功能預測類基因數(shù)量最多(11 205，18.33%)，其次是信號轉導機制(7 618，12.46%)、翻譯后修飾、蛋白質運輸(5 306，8.68%)、未知功能(5 019，8.21%)、轉錄(4 006，6.56%)、翻譯、核糖體結構和生物發(fā)生(3 694，6.04%)、胞內(nèi)轉運、分泌和囊泡運輸(3 616，5.92%)、RNA加工和修飾(3 412，5.59%)。除此之外，負責碳水化合物運輸與代謝有2 510個Unigene(4.11%)，負責脂類運輸和代謝有2 391個Unigene(3.91%)，同時1 984個Unigene負責能源生產(chǎn)與轉換。氨基酸轉運與代謝、細胞骨架、次生代謝產(chǎn)物生物合成與轉運、復制、重組與修復、細胞壁/膜/包膜生物發(fā)生、細胞分裂、核苷酸轉運和代謝、防御機制及輔酶轉運與代謝匹配數(shù)目較少(均小于3%)(圖4)。只有極少數(shù)Unigene負責細胞運動、細胞核結構(均小于0.13%)。由此可見，黃參中功能預測的基因最多，信號轉導、翻譯、修飾及蛋白質運輸參與的基因次之，細胞運動、細胞核結構最少。

圖4 黃參轉錄組Unigene的KOG注釋分類Fig. 4 KOG classification of Sphallerocarpus racills

2.3 基因轉錄因子分析

植物轉錄因子(transcription factor，TF)也稱為反式作用因子，是在轉錄過程中發(fā)揮重要作用的蛋白質，與應答生物和非生物脅迫密切相關。與逆境脅迫相關的轉錄因子主要有MYB類、bZIP類、WRKY類、AP2/EREBP類和NAC類5個大家族[16]。對黃參轉錄因子家族進行分類統(tǒng)計，結果見圖5，共預測到2 370個編碼轉錄因子的Unigene，分布在57個轉錄因子家族中。由圖5可知，C2H2屬于最大家族，Unigene數(shù)量為268個，占總Unigene數(shù)量的11.31%；其次是MYB、WRKY、BHLH、C3H、AP2-EREBP、NAC等轉錄因子家族較多，Unigene數(shù)量分別為220、154、145、130、121、110，分別占9.28%、6.49%、6.10%、5.49%、5.06%、4.64%。轉錄因子的分析可從功能基因組的水平上為進一步開展黃參研究提供數(shù)據(jù)支持。

圖5 黃參轉錄因子家族分布Fig. 5 Transcription factor family distributionin Sphallerocarpus racills

2.4 黃參轉錄組SSR檢測分析

對黃參轉錄組Unigene的簡單序列重復SSR(simple sequence repeat)進行檢測，結果顯示17 308個SSR分布于13 256個Unigene中。二堿基重復的數(shù)量最多(6 721，38.83%)，其次是三堿基重復(6 302，23.21%)，其余為單堿基重復(3 378，19.52%)、五堿基重復(363，2.10%)，最少的是四堿基和六堿基重復，均為272，占1.57%(見表4)。二堿基SSR中，重復頻率最高是TA、AT和TC，最低是GC；三堿基重復總共60種，發(fā)生頻率最高的是CAA、TGTTTG，最低的是CGA和CGG；四堿基重復頻率最高是AAAT、CACT和TTTG，最低是ACAG。重復次數(shù)越多，多態(tài)性越高，多態(tài)性位點較多的是二、三堿基重復，去除SSR位點靠前或靠后的序列，將序列輸入Primer3引物設計軟件，對其中的10對引物進行設計，結果見表5。

表4 黃參SSR分析結果統(tǒng)計Tab. 4 Summary of simple sequence repeat (SSR) in Sphallerocarpus racills

表5 黃參SSR部分引物(10對)設計表Tab. 5 Information of partial primers of development in Sphallerocarpus racills (10 pairs)

3 討論與結論

高通量測序深受研究者歡迎，并越來越多地應用第二代測序技術來解決生物學問題。例如，在基因組水平上對還沒有參考序列的物種進行從頭測序，獲得該物種的參考序列。該技術已廣泛應用到植物特殊功能基因的挖掘與鑒定[17]，為后續(xù)研究和分子育種奠定基礎。本研究采用BGISEQ-500平臺對黃參葉片進行轉錄組測序，獲得了野生黃參功能基因組信息。全長非嵌合序列中存在大量的冗余序列，將其聚類到一起進行去冗余，得到新的一致性序列，然后將非全長序列比對到一致性序列上進行校正，最終得到準確度大于99%的高質量序列，即為轉錄本。本研究最終得到總有效堿基數(shù)6.85 Gb，得到99 981個Unigene，轉錄本總長度113 850 816 bp，平均長度1 138 bp，N50的長度1 874 bp，GC含量39.93%。N50是評價組裝序列完整性的重要指標，N50越長，代表組裝的完整性越好。本研究結果顯示，黃參Q20(堿基正確識別率達99%)序列占96.42%，Q30(堿基正確識別率達99.9%)高質量序列占92.09%，堿基錯誤率為0.01%，低于1%，表明所獲得的黃參轉錄組序列質量較高。

將黃參Unigene比對到7大功能數(shù)據(jù)庫進行注釋，其中：1)黃參Unigene比對到NR數(shù)據(jù)庫共有66 451條注釋成功，與胡蘿卜有較高同源性，而與其他物種的同源性較低。2)GO分析顯示，黃參有78 040條Unigene得到注釋，按功能分為生物過程、細胞組分、分子功能3大類，分別有15、11、14個亞類，最富集的通路主要是在生物學過程中。3)在KOG數(shù)據(jù)庫比對分析結果中，黃參共有61 116條Unigeine獲得了注釋。26個可能的功能大類中，黃參功能預測的基因最多，信號轉導、翻譯、修飾及蛋白質運輸參與的基因次之，細胞運動、細胞核結構最少。4)KEGG是系統(tǒng)分析基因功能、基因組信息的數(shù)據(jù)庫，是進行生物體內(nèi)代謝分析、代謝網(wǎng)絡研究的強有力工具，涵蓋了藥物開發(fā)(drug development)、細胞過程(cellular processes)、環(huán)境信息處理(environmental information processing)、遺傳信息處理(genetic information processing)、人類疾病，僅限動物(human diseases)、代謝(metabolism)、生物系統(tǒng)(organismal systems)等方面。本研究中黃參轉錄組數(shù)據(jù)注釋到KEGG數(shù)據(jù)庫的有富集在20條代謝通路中的51 479個Unigenes，與代謝相關的通路最多，第2位是碳水化合物代謝基因。該發(fā)現(xiàn)有助于揭示黃參淀粉及藥用物質合成途徑，功能基因的分析、代謝通路的注釋等，為后期開展黃參的代謝組學、功能基因組研究奠定了基礎，這部分內(nèi)容將在后續(xù)工作中做進一步研究。

轉錄因子也稱為反式作用因子，是指能夠與真核基因的順式作用元件發(fā)生特異性相互作用，并對基因的轉錄有激活或抑制作用的DNA結合蛋白。轉錄因子能調(diào)控多個與抗逆相關基因的表達，使植物的抗逆性得到改善，可為基礎研究及生產(chǎn)應用提供理論依據(jù)，因此，也逐漸成為植物抗逆機制研究的核心內(nèi)容。植物的抗逆性狀是多基因控制的數(shù)量性狀，多個轉錄因子家族均與植物的抗逆性有重要關系。目前已從高等植物中分離鑒定出數(shù)百種轉錄因子與植物抗逆性密切相關，可調(diào)控植物體感受干旱、高鹽、低溫和病原等信號的相關基因的表達[18]。C2H2型鋅指蛋白主要涉及生長發(fā)育和環(huán)境脅迫應答反應，調(diào)控植物抗逆境脅迫、抗病和生長發(fā)育方面的生物學功能[19]；MYB轉錄因子響應非生物脅迫逆境[20]；WRKY轉錄因子是近幾年研究比較熱的與植物脅迫應答相關的轉錄因子[21]，參與轉錄重編程的調(diào)控[22]，在多種生物和非生物脅迫及諸如水楊酸[23]、赤霉素[24]等植物激素信號轉導中起重要作用。目前，越來越多的研究表明WRKYs是ABA應答信號網(wǎng)絡的關鍵節(jié)點[24]。BHLH是第二大類轉錄因子，不僅影響植物生長發(fā)育，還參與調(diào)控信號轉導和激素合成[25]，結構高度保守，它和NAC一樣，是一類植物特有的轉錄因子，響應不同的生物脅迫和低溫、干旱及高鹽等非生物逆境[25]。生長在祁連山中的黃參抗逆性強，通過研究這些轉錄因子，將有利于進一步揭示這些轉錄因子如何調(diào)控黃參更好地適應逆境脅迫。

簡單序列重復SSR標記是目前最理想的分子標記，廣泛用于開展遺傳多樣性、基因定位、遺傳圖譜構建和比較基因組學研究[15]。本試驗的SSR分析結果為野生黃參的分子標記及遺傳學研究提供了一定的理論基礎。

本次轉錄組測序結果揭示，黃參具有豐富的基因表達，并通過生物信息學分析獲得了黃參基因的注釋信息及代謝通路，獲得的轉錄組信息將為后續(xù)黃參分子標記的開發(fā)和關鍵基因的克隆及功能分析等研究提供科學數(shù)據(jù)。