王暉，謝巖，高玉軍，李季生，王彬彬，高妍夏

(承德醫(yī)學(xué)院蠶業(yè)研究所/河北省高校特產(chǎn)蠶桑應(yīng)用技術(shù)研發(fā)中心，河北 承德 067000)

桑樹是一種分布于世界各地的商品樹種，其經(jīng)濟(jì)效益主要來(lái)自于桑葉，而桑葉主要用于養(yǎng)蠶業(yè)。由于養(yǎng)蠶成本的增加和“東桑西移”戰(zhàn)略的實(shí)施，東中部地區(qū)養(yǎng)蠶業(yè)規(guī)模銳減，導(dǎo)致桑樹產(chǎn)業(yè)亟待轉(zhuǎn)型。桑葚作為桑樹的果實(shí)，含有礦質(zhì)元素、花青素、風(fēng)味物質(zhì)等多種有益成分。近些年桑葚作為一種新型水果受到了極大的關(guān)注，整個(gè)產(chǎn)業(yè)得以迅速發(fā)展，果用桑樹新品種的選育工作也取得了一定進(jìn)展。目前通過(guò)傳統(tǒng)育種手段選育出的果桑新品種主要有大10、白玉王、珍珠白、大馬牙等[1]。隨著多種分子標(biāo)記如RAPD、AFLP、SSR、SNP/Indel的不斷涌現(xiàn)，利用分子標(biāo)記進(jìn)行新品種的選育逐漸成為趨勢(shì)。相比其它類型的分子標(biāo)記，SNP/Indel作為一種新型分子標(biāo)記，其優(yōu)勢(shì)在于基因組及轉(zhuǎn)錄組數(shù)據(jù)中存在著大量該標(biāo)記位點(diǎn)。從茄子[2]、西葫蘆[3]、向日葵[4]、武昌魚[5]、碧桃[7]、松蘿鳳梨[8]、辣椒[9]、陸地棉[10]、大菱鲆[11]、花葉海棠[12]等多種動(dòng)植物的基因組或轉(zhuǎn)錄組數(shù)據(jù)庫(kù)，已成功挖掘出眾多SNP/Indel位點(diǎn)。SNP/Indel標(biāo)記可應(yīng)用于品種選育、基因型分型、種間親緣關(guān)系分析等。Vivek等[6]用桑樹査爾酮合酶基因作為候選基因開(kāi)發(fā)出SNP分子標(biāo)記，但只分析了査爾酮合酶活性與分子標(biāo)記分離方式的相關(guān)性。目前桑樹的基因組測(cè)序已經(jīng)完成，但針對(duì)桑樹SNP/Indel標(biāo)記的開(kāi)發(fā)并不充分，更缺乏適宜于果桑的SNP/Indel標(biāo)記。本研究擬利用RNA-seq技術(shù)對(duì)不同發(fā)育時(shí)期的桑葚進(jìn)行測(cè)序，獲得桑葚SNP/Indel標(biāo)記并進(jìn)行生物信息學(xué)分析，將含有SNP/Indel標(biāo)記的基因?qū)氲蕉鄠€(gè)數(shù)據(jù)庫(kù)中進(jìn)行功能注釋分析，獲得有重要價(jià)值的SNP/Indel分子標(biāo)記，以期為果桑種質(zhì)資源評(píng)價(jià)、品種選育、鑒定等工作的開(kāi)展提供參考依據(jù)。

1 材料與方法

1.1 材料

分別于綠果期、紅果期、黑果期收集果桑品種“安葚”的果實(shí)，清理表面后迅速投入液氮中，-80 ℃保存。轉(zhuǎn)錄組測(cè)序工作由北京康普森生物公司完成。

1.2 數(shù)據(jù)分析方法

將Clean reads通過(guò)Trinity軟件進(jìn)行拼接，獲得Unigene，將其在GO、COG、KEGG數(shù)據(jù)庫(kù)進(jìn)行比對(duì)，獲得其注釋信息。把各個(gè)樣品的數(shù)據(jù)與已組裝好的Unigene進(jìn)行比對(duì)，統(tǒng)計(jì)每個(gè)位點(diǎn)的堿基分布情況，然后結(jié)合GATK和SAMtools從頭挖掘SNPs和INDELs位點(diǎn)。按照質(zhì)量值≥20、測(cè)序深度≥2、SNP間距≥5、50 bp內(nèi)SNP數(shù)<5的過(guò)濾原則，統(tǒng)計(jì)基因的SNP分布情況。Excel 2010進(jìn)行數(shù)據(jù)分析，Sigmaplot 12.5軟件作圖。

2 結(jié)果與分析

2.1 桑葚轉(zhuǎn)錄組數(shù)據(jù)質(zhì)量評(píng)估

3個(gè)時(shí)期桑葚的轉(zhuǎn)錄組數(shù)據(jù)共獲得51 895條Unigene，各時(shí)期桑葚均含有特定表達(dá)的Unigene。桑葚轉(zhuǎn)錄組的序列數(shù)及堿基(A、T、C、G)數(shù)目由高到低依次為綠果期、黑果期、紅果期，序列平均長(zhǎng)度為300 bp。3個(gè)時(shí)期總GC含量差值較小，且總AT含量均高于總GC含量；測(cè)序質(zhì)量值Q20均大于97%、Q30均大于93%(表1)。表明測(cè)序結(jié)果良好，可以進(jìn)行后續(xù)的SNP/Indel數(shù)據(jù)挖掘、分析工作。

表1 測(cè)序產(chǎn)生數(shù)據(jù)質(zhì)量評(píng)估

2.2 桑葚轉(zhuǎn)錄組SNP/Indel類型分析

綠果期13 086條Unigene含有35 024個(gè)SNP/Indel位點(diǎn)，平均每條Unigene 2.68個(gè)SNP/Indel位點(diǎn)；紅果期11 001條Unigene含有28 856個(gè)SNP/Indel位點(diǎn)，平均每條Unigene 2.62個(gè)SNP/Indel位點(diǎn)；黑果期10 666條Unigene含有26 956個(gè)SNP/Indel位點(diǎn)，平均每條Unigene 2.53個(gè)SNP/Indel位點(diǎn)。綠果期桑葚轉(zhuǎn)錄組的SNP類型中，轉(zhuǎn)換、顛換分別占61.75%、38.25%；紅果期分別為61.37%、38.63%；黑果期分別為61.23%、38.77%。綠果期與黑果期轉(zhuǎn)換類型中A/G型最多，顛換類型中A/T型最多。紅果期轉(zhuǎn)換類型C/T型最多，顛換類型中A/T型最多；3個(gè)時(shí)期轉(zhuǎn)換與顛換類型之比均在1.6左右(表2)。

表2 SNP類型統(tǒng)計(jì)

桑葚3個(gè)時(shí)期轉(zhuǎn)錄組基因序列中，每1 000 bp的基因序列均以包含1、2、3個(gè)SNP位點(diǎn)的基因數(shù)目最多，分別占各自時(shí)期含有SNP位點(diǎn)的基因總數(shù)目的86.99%、88.51%、87.84%(圖1)。

圖1 綠果期(a)、紅果期(b)、黑果期(c)桑葚轉(zhuǎn)錄組SNP位點(diǎn)密度分布頻率

桑葚轉(zhuǎn)錄組序列中插入/缺失片段從1到10 bp，Indel數(shù)量逐漸減少，僅在6 bp時(shí)小幅增加；插入/缺失片段長(zhǎng)度以1、2、3 bp為主；綠果期3種長(zhǎng)度類型的插入、缺失片段數(shù)目占總插入、缺失片段數(shù)目的76.53%、57.98%；紅果期為71.07%、58.15%；黑果期為77.31%、61.21%。桑葚轉(zhuǎn)錄組序列中大于10 bp的缺失突變數(shù)量遠(yuǎn)大于插入突變數(shù)量(圖2)。

圖2 Indel類型統(tǒng)計(jì)

2.3 含有SNP位點(diǎn)Unigene GO功能注釋

BLAST結(jié)果顯示，3個(gè)時(shí)期桑葚轉(zhuǎn)錄組中含有SNP/Indel位點(diǎn)的Unigene共有28345條序列注釋到GO數(shù)據(jù)庫(kù)。在注釋后的3個(gè)主要通路中，占比最高的是細(xì)胞組分(50.55%)，其次為分子功能(26.64%)，最后為生物過(guò)程(22.81%)。細(xì)胞組分通路分為12個(gè)小類，細(xì)胞器部分和細(xì)胞通路的基因數(shù)目最多；分子功能通路分為7個(gè)小類，催化活性和結(jié)合通路的基因數(shù)目最多；生物過(guò)程通路分為21個(gè)小類，細(xì)胞過(guò)程和生殖過(guò)程通路的基因數(shù)目最多(圖3)。

圖3 桑葚轉(zhuǎn)錄組序列中SNP/Indel基因的GO功能類別：生物過(guò)程(a)、細(xì)胞組分(b)、分子功能(c)

2.4 含有SNP位點(diǎn)Unigene KOG注釋

將含有SNP/Indel位點(diǎn)的序列在KOG數(shù)據(jù)庫(kù)中進(jìn)行功能注釋，共6 299條Unigene得到注釋，分為4個(gè)大類通路：細(xì)胞進(jìn)程與信號(hào)、信息存儲(chǔ)與加工、代謝及表征不明顯。可進(jìn)一步分為23個(gè)通路。功能預(yù)測(cè)通路1501條基因得到注釋，所有亞類中最多；其次為翻譯后修飾、蛋白轉(zhuǎn)換、分子伴侶通路606條、信號(hào)轉(zhuǎn)導(dǎo)機(jī)制通路504條等；最少為細(xì)胞核結(jié)構(gòu)通路，5條(圖4)。

A—細(xì)胞進(jìn)程與信號(hào)；B—信息存儲(chǔ)與加工；C—代謝；D—表征不明顯。圖4 桑葚轉(zhuǎn)錄組序列中SNP/Indel基因的KOG功能通路

2.5 含有SNP位點(diǎn)Unigene KEGG注釋

通過(guò)比對(duì)分析，共有5737條序列在KEGG數(shù)據(jù)庫(kù)中得到注釋，共包括5個(gè)第一層級(jí)通路：細(xì)胞過(guò)程、環(huán)境信息處理、遺傳信息處理、代謝、有機(jī)系統(tǒng)。5個(gè)第一層級(jí)通路包含19個(gè)第二層級(jí)通路，主要以碳水化合物代謝、翻譯為主(圖5)。

A—細(xì)胞過(guò)程；B—環(huán)境信息處理；C—遺傳信息處理；D—代謝；E—有機(jī)系統(tǒng)。圖5 桑葚轉(zhuǎn)錄組序列中SNP/Indel基因的KEGG功能類別

再進(jìn)一步可分為129個(gè)通路，其中注釋到核糖體、碳代謝、氨基酸生物合成、內(nèi)質(zhì)網(wǎng)蛋白質(zhì)加工等通路的基因最多(圖6)。重要次生代謝物的生物合成通路包含黃酮生物合成、異黃酮生物合成、花青素生物合成、苯丙素生物合成、黃酮與黃酮醇生物合成等通路，上述通路分別有45、1、5、60、11條基因序列得到注釋(圖7)。

圖6 注釋到KEGG數(shù)據(jù)庫(kù)含有SNP/Indel基因數(shù)目最多的20個(gè)通路

圖7 含有SNP/Indel位點(diǎn)的黃酮類與花青素成分合成相關(guān)基因

3 討論

目前，NCBI數(shù)據(jù)庫(kù)中約有6 000多萬(wàn)條EST序列，尚無(wú)桑葚的EST序列數(shù)據(jù)庫(kù)，這不利于理解桑葚的發(fā)育進(jìn)程及重要營(yíng)養(yǎng)物質(zhì)的積累過(guò)程，制約利用分子標(biāo)記快速開(kāi)展優(yōu)良果桑品種的選育工作。隨著測(cè)序技術(shù)的不斷成熟及成本的降低，利用RNA-seq技術(shù)已建立了多個(gè)物種的EST數(shù)據(jù)庫(kù)[13]。本研究對(duì)桑葚進(jìn)行轉(zhuǎn)錄組測(cè)序，共獲得51 895條Unigene序列，重點(diǎn)分析了Unigene中的SNP/Indel位點(diǎn)特征。3個(gè)時(shí)期桑葚轉(zhuǎn)錄組基因序列中SNP位點(diǎn)的出現(xiàn)頻率依次為1個(gè)/462 bp、1個(gè)/482 bp、1個(gè)/478 bp，低于桉樹基因組(1個(gè)/192 bp)[14]、蘋果基因組(1個(gè)/149 bp)[15]、葡萄基因組(1個(gè)/117 bp)[16]、柿樹轉(zhuǎn)錄組(1個(gè)/253 bp)[17];說(shuō)明不同物種基因組、轉(zhuǎn)錄組數(shù)據(jù)中SNP的出現(xiàn)頻率各不相同，推測(cè)其具有物種特異性。即使同一物種，SNP出現(xiàn)頻率也有所不同。橡膠樹轉(zhuǎn)錄組數(shù)據(jù)中SNP位點(diǎn)出現(xiàn)頻率，有報(bào)道是1/1.5 kb[18]；也有報(bào)道是1/5.2 kb[19]，可能是由測(cè)序材料、深度、檢測(cè)軟件的版本及參數(shù)設(shè)置等多種因素所引起。SNP類型及各類型間比例在不同植物基因組間則相對(duì)穩(wěn)定。桉樹轉(zhuǎn)錄組數(shù)據(jù)中，轉(zhuǎn)換、顛換SNP位點(diǎn)分別占60%、占40%[20]；在太平洋白蝦中，分別為66.8%、33.2%[21]；在桑葚中，分別為61%、38%左右(表2)。在自然選擇過(guò)程中，轉(zhuǎn)換突變?cè)诘鞍拙幋a序列中會(huì)產(chǎn)生同義突變，因此通常情況下SNP的轉(zhuǎn)換類型出現(xiàn)頻率都高于顛換類型[22]。

利用現(xiàn)已建立的多種數(shù)據(jù)庫(kù)可以對(duì)已獲得的基因組、轉(zhuǎn)錄組數(shù)據(jù)進(jìn)行注釋，通過(guò)分析，在黃花魚[23]、樹鼩[24]、褐色砂梨[25]、馬鈴薯[26]、向日葵銹菌[27]、菜豆[28]、玉米[29]等多個(gè)物種中篩選出與重要性狀相關(guān)的SNP/Indel分子標(biāo)記。桑葚轉(zhuǎn)錄組數(shù)據(jù)庫(kù)中含SNP/Indel標(biāo)記基因經(jīng)KOG數(shù)據(jù)庫(kù)注釋，發(fā)現(xiàn)翻譯后修飾、蛋白轉(zhuǎn)換、分子伴侶、信號(hào)轉(zhuǎn)導(dǎo)機(jī)制、翻譯、核糖體結(jié)構(gòu)與生物發(fā)生、碳水化合物運(yùn)輸與代謝、轉(zhuǎn)錄、能量生產(chǎn)與轉(zhuǎn)化、胞內(nèi)運(yùn)輸、分泌、囊泡運(yùn)輸、氨基酸運(yùn)輸與代謝等類別中的基因數(shù)有3 019條，占所有注釋基因數(shù)的47.93%(圖4)，說(shuō)明隨著桑葚的成熟，一系列復(fù)雜的分子機(jī)制發(fā)揮作用，引起桑葚顏色、味道、口感發(fā)生劇烈的變化。從KEGG數(shù)據(jù)庫(kù)注釋結(jié)果分析，含有SNP/Indel標(biāo)記基因數(shù)目最多的通路依次為核糖體、碳代謝、氨基酸生物合成、內(nèi)質(zhì)網(wǎng)蛋白質(zhì)加工、剪接體、RNA轉(zhuǎn)運(yùn)等(圖6)，這些通路與物質(zhì)、能量代謝的分子機(jī)制緊密相關(guān)，進(jìn)一步驗(yàn)證了KOG數(shù)據(jù)庫(kù)的注釋結(jié)果。KEGG注釋結(jié)果中共篩選出122個(gè)與花青素、黃酮類合成有關(guān)的含有SNP/Indel位點(diǎn)的基因(圖7)；這些SNP/Indel位點(diǎn)可能與桑葚品質(zhì)形成的關(guān)鍵基因相連鎖，將成為開(kāi)發(fā)選育優(yōu)良果桑品種的分子標(biāo)記的重要來(lái)源。

4 結(jié)論

本研究通過(guò)RNA-seq技術(shù)建立了桑葚轉(zhuǎn)錄組數(shù)據(jù)庫(kù)，證實(shí)桑葚SNP/Indel位點(diǎn)比較豐富。通過(guò)對(duì)獲得的SNP/Indel標(biāo)記特征、基因功能注釋結(jié)果進(jìn)行詳細(xì)研究，有助于下一步篩選出鑒別力強(qiáng)的分子標(biāo)記及開(kāi)發(fā)出高密度SNP基因分型技術(shù)，加快果用桑樹的品種選育、種質(zhì)資源鑒定等工作，更好地服務(wù)桑葚產(chǎn)業(yè)發(fā)展。