菊葉香藜轉錄組單核苷酸多態(tài)性(SNP)信息挖掘及功能注釋分析

付蘇宏，雷 鳴，張勇群*，施 靜，郝豆豆

(1. 西藏自治區(qū)人民政府駐成都辦事處醫(yī)院，四川 成都 610000；2. 西藏大學，西藏 拉薩 850000)

【研究意義】菊葉香藜(Dysphaniaschraderiana)，具有強烈的氣味因此又稱之為臭菜，為藜科(Chenopodiaceae)刺藜屬一年生草本植物，廣泛分布于我國高海拔地區(qū)如西藏、內蒙古、甘肅、青海、陜西等省份。菊葉香藜被用作蒙藥(蒙文名稱：烏努日特-諾衣樂)[1]，具有平喘解痙、清熱解毒、止痛、止癢等功效[2]。此外，菊葉香藜富含植物精油，表現出良好的生物活性，對植物螨蟲、昆蟲、細菌均具有顯著的抑制活性[3-5]，菊葉香藜中最具有潛在天然藥物開發(fā)價值的是其中的活性化合物。此外，由于菊葉香藜對極端環(huán)境具有較好的耐受性，也有研究報道稱菊葉香藜可以改善其生存的生態(tài)環(huán)境[6]。然而，菊葉香藜的研究主要局限在傳統民族藥物方面，近幾年才逐漸開展植物學、植物化學以及生物活性研究，在遺傳學與基因組學方面的研究更是相當匱乏。單核苷酸多態(tài)性(single nucleotide polymorphism, SNP)指的是基因組DNA上的單個堿基發(fā)生變異且變異頻率不小于1 %所引起的DNA序列多態(tài)性[7]。SNP可以作為DNA分子標記，具有位點豐富、分布廣泛、高遺傳穩(wěn)定性、易于檢測和分型等優(yōu)點，現已成為最常用的分子標記之一，可廣泛用于構建高密度遺傳圖譜、植物分子輔助育種、品種鑒定、個體性狀遺傳等方面的研究[8-11]。因此，通過對菊葉香藜SNP分子標記的開發(fā)，可以為鑒定活性物質高產性狀以及極端環(huán)境耐受性狀的亞型提供參考依據。【前人研究進展】菊葉香藜不是模式植物，目前還未獲得其基因組數據以及其他近緣屬植物的基因組測序數據，NCBI數據庫中也尚無菊葉香藜的SNP信息登錄，因此大大限制了菊葉香藜SNP分子標記的開發(fā)。【本研究切入點】由于基因測序技術發(fā)展迅猛，非模式植物可以進行轉錄組測序，其轉錄組數據是開發(fā)SNP分子標記的良好資源，并且這些基于轉錄組SNP分子標記具有信息量大、通用性好的優(yōu)點[12]。【擬解決的關鍵問題】本研究利用菊葉香藜花和葉組織的轉錄組測序數據對其基因內部單核苷酸多態(tài)性(genic SNP)位點進行搜索，分析這些SNP的基本特征，并對SNP-unigene進行功能注釋，為菊葉香藜的SNP分子標記的開發(fā)及其功能研究奠定基礎。

1 材料與方法

1.1 試驗材料與測序數據

菊葉香藜的花和葉組織采集于西藏拉薩市(N29°38′，E91°10′)，采集后立即進行液氮速凍，然后送上海美吉生物醫(yī)藥科技有限公司于IlluminaHiseq4000平臺上進行轉錄組測序，花和葉組織分別獲得24 187 800和27 832 004條干凈數據(Clean reads)，所有的Clean reads通過組裝后共獲得40 142條Unigene (4.72 Gb)。菊葉香藜花和葉組織的轉錄組測序數據已上傳至NCBI數據庫，登錄號分別為SRX3145241和SRX3145242。

1.2 菊葉香藜轉錄組SNP位點檢測

以組裝好的轉錄本為模板序列，將原始序列與其進行比對，利用Samtools (http://samtools.sourceforge.net/)[13]和VarScan v.2.2.7 (http://varscan.sourceforge.net/)[14]軟件尋找候選的SNP位點。

1.3 SNP-unigene序列功能注釋分析

對含有SNP的unigene序列(SNP-unigene)進行功能注釋，探究SNP-unigene的生物學功能，注釋所選數據庫有基因本體數據庫(gene ontology database, GO)、蛋白直系同源簇數據庫(Clusters of Orthologous Groups of proteins, COG)、京都基因和基因組百科全書數據庫(kyoto encyclopedia of genes and genomes, KEGG)數據庫，期望值E設定為1e-5[15]。

2 結果與分析

2.1 菊葉香藜轉錄組SNP鑒定

使用Samtools軟件和VarScan v.2.2.7軟件對菊葉香藜花和葉的轉錄組序列進行SNP位點檢測，在花和葉轉錄組中分別鑒定到889個和673個SNP位點。菊葉香藜花轉錄組中的889個SNP位點中轉換(Transition)占有560個，顛換(Transversion)占有329個；6種核苷酸變異中屬于轉換的A/G和C/T發(fā)生頻率最高，所占比例均為31.50 %，屬于顛換的A/T、A/C、T/G和C/G的比例則分別為15.00 %、6.07 %、7.76 %和8.21 % (圖1A)；889個SNP中分布在密碼子的第3位(Third codon)的比例最高，占花轉錄組SNP總數的40.16 %，然后依次為3端非編碼區(qū)(3-UTR，15.52 %)、密碼子第1位(First codon，12.60 %)、密碼子第2位(Second codon，12.26 %)和5端非翻譯區(qū)(5-UTR，5.62 %)，有13.84 %分布位置無法確定(圖1B)。菊葉香藜葉轉錄組中SNP變異類型與分布位置趨勢與花組織轉錄組SNP類似，但數量較少(圖1)。

圖1 菊葉香藜轉錄組中genic SNP概況Fig.1 Summary of genic SNP within transcriptome ofD.schraderiana

2.2 SNP-unigene序列注釋

從菊葉香藜花和葉組織中鑒定出的SNP位點總共位于643條SNP-unigene上，為進一步探索這些SNP-unigene的生物學功能，將這643條SNP-unigene進行GO注釋、COG注釋和KEGG注釋，結果顯示，643條SNP-unigene中總共有440 (68.43 %)條SNP-unigene注釋到了GO數據庫、COG數據庫和KEGG數據庫。在有注釋信息的440條SNP-unigene中，注釋到KEGG數據庫的有232條，注釋到COG數據庫的有370條，注釋到GO數據庫的有343條；同時注釋到GO數據庫、COG數據庫和KEGG數據庫的總共有179條(40.7 %)；同時注釋到KEGG數據庫和COG數據庫197條(44.8 %)，同時注釋到COG數據庫和GO數據庫285條(64.8 %)，同時注釋到GO數據庫和KEGG數據庫201條(45.7 %)；只注釋到KEGG數據庫的有13條(64.8 %)，只注釋到COG數據庫的有67條(15.2 %)，只注釋到GO數據庫的有35條(8.0 %) (圖2)。

2.3 SNP-unigene序列GO分類

GO分類可分為3個方面：生物過程(Biological Process，BP)、分子功能(Molecular Function，MF)及細胞組分(Cellular Component，CC)，而在該3個大分支下又可分為更加細致的功能，通過對SNP-unigene進行GO注釋有助于理解基因背后所代表的生物學意義。使用blast2go軟件將643條SNP-unigene與GO數據庫比對，其中總共有342條SNP-unigene (53.19 %)注釋到GO條目中，分布于分子功能、細胞組分和生物學過程的level2 GO條目分別有397、477和756個。在生物過程中，結合活性(162條)和催化活性(168條)中所涉及的SNP-unigene較多；在細胞組分中，細胞(358條)中所涉及的SNP-unigene較多；在生物學過程中，代謝過程(379條)和細胞過程(165條)中所涉及的SNP-unigene較多(圖3)。

圖2 SNP-unigene注釋信息統計Venn圖Fig.2 Venn diagram of annotation information of SNP-unigene

2.4 SNP-unigene序列COG分類

將SNP-unigene與COG數據庫的比對，可以獲取序列的功能注釋、歸類以及蛋白進化等信息。COG注釋根據功能可分為25類(圖3，A-Z表示，除X)，總共有370條(57.54 %)SNP-unigene在COG

圖3 GO功能分類統計Fig.3 GO functional classifications

A：RNA加工和修飾；B：染色質結構和動力學；C：能量產生和轉換；D：細胞周期控制、細胞分裂和染色體分離；E：氨基酸轉運和新陳代謝；F：核酸轉運和代謝；G：碳運輸和新陳代謝；H：輔酶運輸和代謝；I：脂類轉運和代謝；J：翻譯，核糖體結構和合成；K：轉錄；L：復制,重組和修復；M：細胞壁膜核膜的合成；N：細胞機動性；O：翻譯后修飾，蛋白翻轉和分子伴侶；P：無機離子轉運和代謝；Q：二級代謝生物加工,轉運和分解代謝；R：一般預測功能組；S：未知功能；T：信號傳導機制；U：胞內的交換，分泌和膜泡輸送；V：防御機制；W：胞外結構；Y：核結構；Z：細胞骨架圖4 COG分類統計Fig.4 COG functional classification

圖5 KEGG注釋統計Fig.5 Summary of SNP-unigene with mainly KEGG annotation

數據庫中找到了相對應的注釋信息，得到385個COG注釋結果。通過對每一類功能的數量進行統計發(fā)現，除了“細胞機動性”(N)和“核結構”(Y)兩個功能沒有注釋有SNP-unigene以外，其余的23個COG功能均注釋有SNP-unigene，由此可見，被注釋的370條SNP-unigene功能種類比較全面，覆蓋了大多數(23/25)的生命活動過程及功能。在注釋有SNP-unigene的COG功能分類中，“翻譯，核糖體結構和合成”是數量最多的一大類，包括65個SNP-unigene；然后依次為“翻譯后修飾，蛋白翻轉和分子伴侶”、“一般預測功能組”、“轉錄”，包含的SNP-unigene分別為32、31、27條(圖4)。

2.5 SNP-unigene序列KEGG分析

KEGG數據庫可以將基因組的信息與功能信息聯系起來，更加系統地分析和破譯基因的功能。將unigene與KEGG數據庫進行比對，可以獲得與unigene相對應的KO編號，根據KO編號則可進一步獲得該unigene可能參與的具體生物學通路。對643條SNP-unigene進行KEGG通路注釋，其中總共有232條SNP-unigene注釋到241條KEGG代謝通路中。將這232條SNP-unigene根據所參與的KEGG代謝通路劃分為5個分支：代謝(A, Metabolism)、遺傳信息處理(B, Genetic Information Processing)、環(huán)境信息處理(C, Environmental Information Processing)、細胞過程(D, Cellular Processes)、有機系統(E, Organismal Systems)，其五大類的KEGG通路統計情況如圖5所示。

篩選注釋基因比例(占所有KEGG注釋unigene)大于1 %的71條KEGG通路列于表1。從表1可知，注釋到“代謝途徑”的SNP-unigene最多，可達到78條，占所有KEGG注釋SNP-unigene的33.62 %，然后依次為“核糖體”、“次級代謝產物的生物合成”、“微生物在不同環(huán)境中的代謝”、“嘌呤代謝”以及“RNA運輸”，分別注釋有48、36、19、13和13條SNP-unigene。通過KEGG的注釋分析，可以著重研究這些具有代表性的生物學進程，探索菊葉香藜中這些生物學進程所涉及的功能。

表1 主要KEGG pathway注釋統計

續(xù)表1 Continued table 1

序號No.通路IDPathway ID通路PathwaySNP-unigene數量No. of SNP-unigene比例(%)Ratio 序號No.通路IDPathway ID通路PathwaySNP-unigene數量No. of SNP-unigene比例(%)Ratio20path:ko05169EB病毒感染Epstein-Barr virus infection62.5956path:ko02010ABC轉運蛋白ABC transporters31.2921path:ko05205癌癥中的蛋白聚糖Proteoglycans in cancer52.1657path:ko04010MAPK信號通路MAPK signaling pathway31.2922path:ko00480谷胱甘肽代謝Glutathione metabolism52.1658path:ko04510焦點黏連Focal adhesion31.2923path:ko05145弓形體病Toxoplasmosis52.1659path:ko00900萜骨架生物合成Terpenoid backbone biosynthesis31.2924path:ko00240嘧啶代謝Pyrimidine metabolism52.1660path:ko04066HIF-1信號通路HIF-1 signaling pathway31.2925path:ko04144內吞作用Endocytosis52.1661path:ko05200癌癥中的信號通路Pathways in cancer31.2926path:ko04141內質網蛋白加工Protein processing in endoplasmic reticulum52.1662path:ko00620丙酮酸鹽代謝Pyruvate metabolism31.2927path:ko04712晝夜節(jié)律-植物Circadian rhythm - plant41.7263path:ko04113減數分裂-酵母Meiosis - yeast31.2928path:ko05204化學致癌性Chemical carcinogenesis41.7264path:ko04722神經營養(yǎng)蛋白信號通路Neurotrophin signaling pathway31.2929path:ko00982藥物代謝-細胞色素P450Drug metabolism - cytochrome P45041.7265path:ko00360苯丙氨酸代謝Phenylalanine metabolism31.2930path:ko00980細胞色素P450的異種生物代謝Metabolism of xenobiotics by cytochrome P45041.7266path:ko00520氨基糖和核苷酸糖代謝Amino sugar and nucleotide sugar metabolism31.2931path:ko00250丙氨酸,天冬氨酸和谷氨酸代謝Alanine, aspartate and glutamate metabolism41.7267path:ko03018RNA降解RNA degradation31.2932path:ko04921催產素信號通路Oxytocin signaling pathway41.7268path:ko05166HTLV-I 感染HTLV-I infection31.2933path:ko00680甲烷代謝Methane metabolism41.7269path:ko04110細胞周期Cell cycle31.2934path:ko00630乙醛酸和二羧酸代謝Glyoxylate and dicarboxylate metabolism41.7270path:ko00010糖酵解/葡萄糖異生作用Glycolysis/Gluconeogenesis31.2935path:ko05034酒精中毒Alcoholism41.7271path:ko04626植物-病原體相互作用Plant-pathogen interaction31.2936path:ko05168單純性皰疹感染Herpes simplex infection41.72

3 討 論

隨著新一代測序技術的成熟以及成本的降低，對于非模式生物來說，轉錄組序列是開發(fā)SNP分子標記最有效的策略之一[16]。目前，已經有大量的非模式生物通過轉錄組序列進行了SNP分子標記的研究和開發(fā)，在松蘿鳳梨葉片轉錄組數據中總共挖掘到了69 570個基因內部SNP，發(fā)生頻率為1/579 bp[17]；同時，在龍眼轉錄組數據中挖掘并驗證了60個SNP分子標記并進行了基因型鑒定[18]；通過對盾葉薯蕷轉錄組數據進行挖掘，得到了124 692個SNP位點，每條unigene中平均含有3.36個SNP位點[19]。大量研究表明，挖掘轉錄組數據中的SNP信息是非模式生物開發(fā)SNP分子標記可行且有效的方法。

本次研究中，菊葉香藜SNP位點處堿基變異類型主要為C/T和A/G，以轉換為主，花組織和葉組織中轉換和顛換的比值分別為1.70和1.68，均遠遠大于理論值的0.5。生物中SNP位點堿基變異的轉換和顛換的比值往往會大于理論值的0.5，這種現象稱為“轉換偏差”[20]，而這一現象也恰恰說明了SNP位點的堿基突變并不是隨機發(fā)生的，可能與環(huán)境選擇和進化有關。除了有少量的SNP分布位置不確定以外，位于非編碼區(qū)的SNP在花和葉轉錄組序列中分別占21.15 %和21.10 %，而位于編碼區(qū)的SNP在花和葉轉錄組序列中則分別占到了67.27 %和67.46 %。由此可見，菊葉香藜的SNP位點大多位于編碼區(qū)內，此類SNP位點處的堿基變異可能會引起氨基酸序列的變化，進一步引起菊葉香藜生物性狀的變化，因此，開發(fā)此類的SNP標記具有重要意義。

挖掘菊葉香藜轉錄組數據獲得了大量的SNP信息，通過對這些SNP-unigene進行功能注釋則有利于探索菊葉香藜在哪些生物和代謝途徑中可能存在個體性狀差異，從而將這些SNP和其表型關聯起來進行更深入的研究。在菊葉香藜SNP-unigene功能注釋中，菊葉香藜SNP-unigene涉及較多的功能主要與代謝、核糖體、次生代謝產物的生物合成相關，可針對這些SNP-unigene設計引物進行SNP驗證并做表型關聯研究，所以針對這類功能的SNP分子標記進行開發(fā)研究，可以為以后分子標記輔助育種等工作奠定堅實的基礎。

4 結 論

本研究挖掘獲得了菊葉香藜轉錄組測序數據中的SNP信息，并獲取了SNP-unigene的GO功能注釋和KEGG通路注釋信息，本研究成果將積極地促進菊葉香藜的遺傳圖譜構建、遺傳育種、遺傳多樣性分析等方面的研究。