西藏山溪鯢轉(zhuǎn)錄組組裝與分析

朱文文，郁 川，熊建利，黃 勇

(1.洛陽(yáng)職業(yè)技術(shù)學(xué)院食品與藥品學(xué)院，河南 洛陽(yáng) 471000；2.河南省動(dòng)物疾病與公共衛(wèi)生工程研究中心，河南 洛陽(yáng) 471000；3.綿陽(yáng)師范學(xué)院，生態(tài)安全與保護(hù)四川省重點(diǎn)實(shí)驗(yàn)室，四川 綿陽(yáng) 621000；4.河南科技大學(xué)動(dòng)物科學(xué)院水生動(dòng)物適應(yīng)與進(jìn)化實(shí)驗(yàn)室，河南 洛陽(yáng) 471000)

西藏山溪鯢(Batrachuperustibetanus)隸屬有尾目(Caudata)、小鯢科(Hynobius)、山溪鯢屬(Batrachuperus)，主要分布在青海、甘肅、陜西、四川、西藏等海拔1 500～4 000 m、植被較為豐富的山區(qū)溪流中或泉水石堆下[1-2]。西藏山溪鯢是我國(guó)珍貴、稀有的瀕危水生野生動(dòng)物，已被列入《國(guó)家保護(hù)的有益的或者有重要經(jīng)濟(jì)、科學(xué)研究?jī)r(jià)值的陸生野生動(dòng)物名錄》。由于西藏山溪鯢具有很高的藥用價(jià)值，其也被列入我國(guó)傳統(tǒng)藏藥藥典；同時(shí)具有重要的科研和生態(tài)價(jià)值[3]。近年來(lái)，由于環(huán)境惡化、生態(tài)條件退化和人為過(guò)度開(kāi)發(fā)利用該物種資源，導(dǎo)致其生存空間正在逐漸縮小，種群數(shù)量明顯減少，在IUCN(International Union for Conservation of Nature)紅色名錄中被列為易危的物種[4]。目前，有關(guān)西藏山溪鯢轉(zhuǎn)錄組的研究尚未見(jiàn)報(bào)道。而了解西藏山溪鯢轉(zhuǎn)錄組基因信息，能為后續(xù)科學(xué)合理地利用與保護(hù)該物種基因資源提供理論基礎(chǔ)。

轉(zhuǎn)錄組高通量測(cè)序是近幾年發(fā)展起來(lái)的新技術(shù)，即RNA-Sequencing(RNA-Seq)，具有處理數(shù)據(jù)量大、運(yùn)行成本低、靈敏度高等優(yōu)點(diǎn)，現(xiàn)已成為發(fā)掘物種功能基因的重要研究手段之一[5-6]。該高通量技術(shù)可在沒(méi)有該物種已知全部基因組序列信息的前提下，在較短時(shí)間內(nèi)準(zhǔn)確地得到特定組織或者細(xì)胞在特殊狀態(tài)下全部的轉(zhuǎn)錄組信息，能完整識(shí)別該條件下已知基因的表達(dá)、生理狀態(tài)與特定的分子機(jī)制調(diào)控過(guò)程，并能辨別一些未知的轉(zhuǎn)錄本和遺傳標(biāo)記信息等，在非模式生物轉(zhuǎn)錄組的研究中得到了廣泛應(yīng)用，為進(jìn)一步研究生物學(xué)提供了更全面、便利的平臺(tái)[7-9]。目前，有關(guān)山溪鯢屬物種的研究均集中在屬種形態(tài)分類(lèi)、系統(tǒng)發(fā)育與進(jìn)化和生理生化等方面[10-15]。為挖掘西藏山溪鯢基因數(shù)據(jù)和功能，本研究利用高通量測(cè)序技術(shù)對(duì)西藏山溪鯢進(jìn)行轉(zhuǎn)錄組測(cè)序，并結(jié)合現(xiàn)代生物信息學(xué)分析方法對(duì)測(cè)序得到的序列進(jìn)行拼接、組裝和功能注釋分析。得到的數(shù)據(jù)有利于更全面了解西藏山溪鯢轉(zhuǎn)錄組信息，方便科技工作者實(shí)現(xiàn)數(shù)據(jù)資源共享，為后期開(kāi)展西藏山溪鯢分子遺傳學(xué)及生物多樣性研究提供基礎(chǔ)資料。

1 材料與方法

1.1 材料

2017年7月，從四川百靈山采集4只外表無(wú)傷、體長(zhǎng)130～150 mm的西藏山溪鯢成體(2♂、2♀)為研究對(duì)象。

1.2 樣品收集

利用MS-222將標(biāo)本麻醉后，分別取4尾成鯢的脾臟、肌肉、腎臟、肝臟、心臟、腸道、皮膚和性腺組織，每個(gè)組織取樣約10 mg，最后將所有組織樣品混為1個(gè)樣品，于-80℃超低溫冰箱保存、備用。

1.3 組織總RNA檢測(cè)

將混合的西藏山溪鯢組織樣品迅速置于裝有液氮的研缽中，研磨成粉狀，根據(jù)Takara公司提供的Trizol Reagent操作說(shuō)明書(shū)完成總RNA的抽提。獲得的總RNA 再經(jīng)過(guò)1.2%瓊脂糖凝膠電泳和Nanodrop-2000核酸蛋白測(cè)定儀檢測(cè)總RNA的完整性和純度。經(jīng)檢驗(yàn)合格的總RNA樣品再進(jìn)行后續(xù)的轉(zhuǎn)錄組測(cè)序。檢測(cè)標(biāo)準(zhǔn)定義為：總RNA量≥10 μg，OD260/280為1.75～2.10，28 S∶18 S≥1.5∶1.0，RIN值≥8.0。

1.4 轉(zhuǎn)錄組測(cè)序與分析

樣品檢測(cè)合格后，利用 Oligo (dT)磁珠法純化出mRNA，然后將mRNA截為片段，經(jīng)過(guò)PCR方法得到西藏山溪鯢cDNA文庫(kù)，最后建好的cDNA文庫(kù)利用Illumina HiSeq 4000平臺(tái)技術(shù)進(jìn)行上機(jī)測(cè)序，由杭州聯(lián)川生物技術(shù)股份有限公司完成測(cè)序。原始 reads 經(jīng)過(guò)去接頭并且過(guò)濾掉低質(zhì)量及長(zhǎng)度過(guò)短序列后，得到高質(zhì)量測(cè)序數(shù)據(jù)(clean reads)。利用 Trinity 軟件的Paired-end拼接方法對(duì)clean reads進(jìn)行Denovo組裝，得到unigenes。由于西藏山溪鯢目前沒(méi)有基因組數(shù)據(jù)，本研究以報(bào)道的墨西哥蠑螈(Ambystomamexicanum)(https://www.ncbi.nlm.nih.gov/data-hub/genome/GCA_002915635.3/)基因組序列作參考基因組，使用STAR 軟件進(jìn)行比對(duì)[16]，再用StringTie 軟件基于參考基因組注釋文件對(duì)所有轉(zhuǎn)錄本進(jìn)行整合組裝[17]。使用MiscroSAtellite(MISA)軟件進(jìn)行SSR鑒定和分析。應(yīng)用Rsem軟件進(jìn)行基因表達(dá)定量分析，基因表達(dá)量采用FPKM值(Fragments per kilobase of transcript per million fragments mapped)表示。最后，利用Swiss-prot (Swiss prot protein database)、Nr (Non-redundant protein sequences)、KEGG (Kyoto encyclopedia of genes and genomes)、KOG (Eukaryotic ortholog groups)和 Pfam (Protein families database)和GO (Gene ontology)6個(gè)公共數(shù)據(jù)庫(kù)與組裝得到的西藏山溪鯢unigenes進(jìn)行Blast序列比對(duì)[18]。選擇閥值條件為E value <1e-10，進(jìn)行功能注釋。

2 結(jié)果與分析

2.1 RNA檢測(cè)

提取的混合組織總RNA 樣品呈現(xiàn)完整清晰的28 S、18 S和5 S帶型，OD260/280值為2.03，28 S∶18 S值接近1.80，且RIN值為8.0，說(shuō)明提取到的總RNA質(zhì)量較高，符合后續(xù)轉(zhuǎn)錄組測(cè)序建庫(kù)要求。

2.2 測(cè)序質(zhì)量控制

采用Illumina Hiseq 4000測(cè)序得到西藏山溪鯢轉(zhuǎn)錄組的數(shù)據(jù)。對(duì)獲得的原始數(shù)據(jù)經(jīng)過(guò)測(cè)序質(zhì)量控制，得到49 924 038 bp的clean reads，包含7.49 G的堿基數(shù)據(jù)。測(cè)序質(zhì)量顯示：堿基Q20(序列質(zhì)量不低于20的堿基所占百分比)占96.99%，Q30占93.61%(>85%)，GC含量平均值為47.18%，表明測(cè)序堿基組成的結(jié)果較好，組裝質(zhì)量完整性較高，能用于下一步分析。

2.3 轉(zhuǎn)錄組數(shù)據(jù)組裝與分析

在去除低質(zhì)量數(shù)據(jù)和進(jìn)行質(zhì)控后得到的clean reads，采用Denovo方法進(jìn)行序列拼接后，總共獲得43 626條轉(zhuǎn)錄本，序列長(zhǎng)度為46 293 115 bp，其中N50片段序列的長(zhǎng)度為1 822 bp，平均長(zhǎng)度為1 061 bp。在獲得轉(zhuǎn)錄本序列的基礎(chǔ)上，經(jīng)Trinity軟件進(jìn)行組裝，參數(shù)選用 Trintity 的省缺參數(shù) Kmer=25[19]，然后拼接好的片段進(jìn)一步合成。最終得到了36 252條unigenes，共33 976 485 bp，序列大小范圍為201～20 766 bp，得到的平均長(zhǎng)度為937 bp，其中N50為1 549 bp(表1)。對(duì)每條unigenes長(zhǎng)度統(tǒng)計(jì)相應(yīng)的unigenes數(shù)量，其中長(zhǎng)度范圍在200～500 bp的最多，有16 775條，占總數(shù)的46.27%；其次長(zhǎng)度在500～1 000 bp之間的有9 091條，占總數(shù)的25.1%；長(zhǎng)度在1 000～2 000 bp之間的有6 263條，占總數(shù)的17.28%；長(zhǎng)度大于2 000 bp以上的最少，有4 123條，占總數(shù)11.37%，隨長(zhǎng)度增加基因數(shù)量逐漸減少，說(shuō)明測(cè)序的序列、數(shù)據(jù)拼接和組裝質(zhì)量較高。

表1 測(cè)序組裝的序列結(jié)果統(tǒng)計(jì)

2.4 SNP檢測(cè)分析

單核苷酸多態(tài)性定義為基因組上單個(gè)核苷酸產(chǎn)生的變異，是進(jìn)行分子標(biāo)記鑒定、輔助育種和遺傳圖譜構(gòu)建等非常重要的一種遺傳標(biāo)記方法。本研究利用轉(zhuǎn)錄組作為參考序列，使用BWA和Samtools軟件對(duì)西藏山溪鯢的外顯子區(qū)域進(jìn)行SNP發(fā)掘，結(jié)果顯示總共得到 3 100個(gè)SNP 位點(diǎn)，包括A-G、C-T、A-C、A-T、C-G和G-T六種顛換類(lèi)型的SNP。在所有顛換類(lèi)型的SNP中，A-G和C-T兩種顛換類(lèi)型的比例最高，占所有SNP位點(diǎn)的58.94%；C-G顛換類(lèi)型占的比例最低，僅為6.87%。A-C、A-T和G-T這3種 SNP顛換類(lèi)型有相似的比例，占總量的34.19%(表2)。

表2 SNP位點(diǎn)的分類(lèi)

2.5 CDS預(yù)測(cè)與微衛(wèi)星SSR位點(diǎn)分析

通過(guò)Genscan軟件預(yù)測(cè)其CDS序列，預(yù)測(cè)結(jié)果顯示：總共檢測(cè)出17 472條unigenes可被編碼，占全部36 252條unigenes的48.2%，未檢測(cè)到CDS的unigenes有18 780條。長(zhǎng)度在200～500 bp所編碼的氨基酸占比最高，預(yù)測(cè)到超過(guò)2 000 bp的編碼氨基酸序列有1 045條，能編碼氨基酸平均長(zhǎng)度為367.6 bp。使用MISA軟件對(duì)36 252條unigenes進(jìn)行搜索SSR位點(diǎn)，總共有2 779條序列被檢測(cè)到具有SSR位點(diǎn)(表3)。其中，單核苷酸至六核苷酸重復(fù)類(lèi)型均有被檢測(cè)到，SSR類(lèi)型單核苷酸出現(xiàn)頻率最多，有1 473個(gè)；其次為二核苷酸，有715個(gè)；490個(gè)SSR具有三核苷酸位點(diǎn)；四和五核苷酸SSR位點(diǎn)的數(shù)量分別為56個(gè)和38個(gè)；SSR位點(diǎn)為六核苷酸的數(shù)量最少，僅為7個(gè)。這些SSR位點(diǎn)可為后續(xù)進(jìn)行分子標(biāo)記鑒定提供引物設(shè)計(jì)基礎(chǔ)。

表3 SSR位點(diǎn)統(tǒng)計(jì)

2.6 基因表達(dá)量分析

在轉(zhuǎn)錄組測(cè)序中，用RPKM方法計(jì)算基因的表達(dá)水平(表4)，表示在每百萬(wàn)reads中來(lái)自某一基因每千堿基長(zhǎng)度的read讀數(shù)。本研究以RPKM≥0.1作為基因表達(dá)標(biāo)準(zhǔn)，在已獲得的36 252個(gè)unigenes中，有基因表達(dá)量的序列為36 226條。對(duì)不同RPKM 區(qū)間的基因數(shù)量進(jìn)行統(tǒng)計(jì)發(fā)現(xiàn)，其中RPKM值在3.57～15.00的基因最多，為17 099條，占到了47.2%；其次是RPKM值位于0.30～3.57的基因，數(shù)量為13 868條；RPKM值位于15～60的基因有4 000條，高表達(dá)的基因RPKM值>60為1 243條；而RPKM值在0.1～0.3的低表達(dá)基因最少，僅16條，占所有表達(dá)基因的0.04%。上述結(jié)果表明Illumina HiSeq測(cè)序技術(shù)能夠檢測(cè)到極低水平基因的表達(dá)。

表4 FPKM值密度分布

2.7 Unigenes功能注釋

在6個(gè)數(shù)據(jù)庫(kù)中對(duì)組裝得到的36 252條unigenes進(jìn)行Blast比對(duì)，并進(jìn)行功能注釋。注釋成功的unigene基因數(shù)目在不同數(shù)據(jù)庫(kù)中所占比例有所差別。如表5所示，在Swiss-port數(shù)據(jù)庫(kù)獲得注釋的unigenes有16 465條，占總數(shù)的45.42%；在Nr數(shù)據(jù)庫(kù)中獲得注釋的unigenes有18 749條，注釋比例最大，達(dá)到了51.72%；在Pfam數(shù)據(jù)庫(kù)中獲得注釋的unigenes有13 983條，占38.57%；在KEGG數(shù)據(jù)庫(kù)注釋的unigenes有10 607條，占29.26%；在KOG數(shù)據(jù)庫(kù)中獲得注釋的unigenes有15 704條，占43.32%；在GO數(shù)據(jù)庫(kù)中獲得注釋的unigenes有14 242條，占39.29%。在Nr數(shù)據(jù)庫(kù)中，按物種分布統(tǒng)計(jì)，與西部錦龜(Chrysemyspicta)匹配度最多，為13.6%；其次是熱帶爪蟾(Xenopustropicalis)和綠海龜(Cheloniamydas)，分別為9.5%和7.8%，最低的是非洲爪蟾(Xenopuslaevis)，為4.2%，而與其他物種蛋白質(zhì)無(wú)匹配的unigenes占55.8%(圖1)。

表5 序列功能注釋

2.8 KOG功能分類(lèi)

將獲得的西藏山溪鯢unigenes 在COG數(shù)據(jù)庫(kù)中進(jìn)行功能注釋?zhuān)煞譃?5類(lèi)(圖2)，在COG中注釋到的unigenes涉及功能類(lèi)別較廣，與生命活動(dòng)相關(guān)的占大部分。其中，基因數(shù)注釋最多的是一般功能預(yù)測(cè)類(lèi)，有2 490條。其次，信號(hào)轉(zhuǎn)導(dǎo)機(jī)制與翻譯后修飾、蛋白質(zhì)周轉(zhuǎn)、伴侶類(lèi)基因，分別有2 210條和1 140條；表明遺傳信息的傳遞在西藏山溪鯢生理活動(dòng)中極為活躍。值得注意的是，注釋到細(xì)胞運(yùn)動(dòng)類(lèi)基因最少，僅有40條，這說(shuō)明該物種與它們高海拔獨(dú)特的生活環(huán)境和遷徙能力弱有關(guān)。

2.9 GO功能注釋

根據(jù)得到的注釋信息進(jìn)行分類(lèi)，共有14 242條unigenes被注釋?zhuān)?9.29%。按照GO功能分類(lèi)方式將注釋到的unigenes主要分為3大類(lèi)(生物過(guò)程、細(xì)胞組分和分子功能)，如圖3所示。這3個(gè)大類(lèi)別又被詳細(xì)地劃分為50個(gè)功能亞類(lèi)小組。其中大類(lèi)生物過(guò)程包含25個(gè)不同的亞類(lèi)功能組，這也是三大類(lèi)中所含功能類(lèi)別最多的一類(lèi)，注釋到轉(zhuǎn)錄、DNA依賴性和轉(zhuǎn)錄調(diào)控、DNA依賴性的unigenes占最多，分別有1 183條和932條，占比分別為72.27%和60.13%；而與RNA剪接相關(guān)unigenes所占數(shù)量最少，僅140條，占12.02%；在細(xì)胞組分中，有15個(gè)亞類(lèi)，注釋到細(xì)胞核的unigenes最多(2 798條)，占90.12%；注釋到細(xì)胞膜的unigenes最少，僅為138條。分子功能類(lèi)別中又劃分為10個(gè)亞類(lèi)，注釋到ATP結(jié)合相關(guān)功能的unigenes數(shù)量最多，有1 855條，其次是注釋到鋅離子結(jié)合的unigenes，有1 668條。而注釋到蛋白絲氨酸/蘇氨酸激酶活性相關(guān)的unigenes數(shù)量最少，有377條序列。

注：A.RNA 加工和修飾;B.染色質(zhì)結(jié)構(gòu)與動(dòng)力學(xué);C.能量產(chǎn)生和轉(zhuǎn)換;D.細(xì)胞周期調(diào)控、細(xì)胞分裂、染色體;E.氨基酸運(yùn)輸和代謝;F.核苷酸運(yùn)輸和代謝;G.碳水化合物的運(yùn)輸和代謝;H.輔酶運(yùn)輸和代謝;I.脂質(zhì)運(yùn)輸和代謝;J.翻譯、核糖體結(jié)構(gòu)和生物合成;K.轉(zhuǎn)錄;L.復(fù)制、重建和修復(fù);M.細(xì)胞壁/細(xì)胞膜/膜結(jié)構(gòu)的生物合成;N.細(xì)胞運(yùn)動(dòng);O.翻譯后修飾、蛋白質(zhì)周轉(zhuǎn)、伴侶;P.無(wú)機(jī)離子轉(zhuǎn)運(yùn)與代謝;Q.次生代謝產(chǎn)物的合成、轉(zhuǎn)運(yùn)和代謝;R.普通功能預(yù)測(cè);S.未知功能;T.信號(hào)轉(zhuǎn)導(dǎo)機(jī)制;U.胞內(nèi)運(yùn)輸、分泌和囊泡運(yùn)輸;V.防御機(jī)制;W.胞外結(jié)構(gòu);Y.核結(jié)構(gòu);Z.細(xì)胞骨架。

注：1.轉(zhuǎn)錄、DNA依賴性;2.轉(zhuǎn)錄調(diào)控、DNA依賴性;3.蛋白質(zhì)轉(zhuǎn)運(yùn);4.多細(xì)胞器官發(fā)育;5.細(xì)胞分化;6.凋亡;7.細(xì)胞黏附;8.蛋白質(zhì)水解;9.信號(hào)轉(zhuǎn)導(dǎo);10.細(xì)胞分化;11.細(xì)胞周期;12.有絲分裂;13.小G蛋白介導(dǎo)信號(hào)轉(zhuǎn)導(dǎo);14.DNA 修復(fù);15.轉(zhuǎn)運(yùn); 16.RNA加工;17.翻譯;18.轉(zhuǎn)錄正調(diào)控;19.胞內(nèi)信號(hào)轉(zhuǎn)導(dǎo);20.染色質(zhì)修飾;21.轉(zhuǎn)錄負(fù)調(diào)控;22.細(xì)胞內(nèi)蛋白質(zhì)轉(zhuǎn)運(yùn);23.跨膜輸送;24.精子形成;25.RNA剪切;26.細(xì)胞核;27.必需膜;28.細(xì)胞漿;29.胞液;30.細(xì)胞質(zhì)膜;31.內(nèi)質(zhì)網(wǎng)膜;32.線粒體;33 胞外區(qū);34.胞核;35.核質(zhì);36.細(xì)胞骨架;37.高爾基體膜;38.微管;39.高爾基氏復(fù)合體;40.細(xì)胞膜;41.ATP結(jié)合;42.鋅離子結(jié)合;43.蛋白質(zhì)結(jié)合;44.DNA 結(jié)合;45.金屬離子結(jié)合;46.RNA結(jié)合;47.鈣離子結(jié)合;48.結(jié)合;49.特異序列DNA結(jié)合;50.蛋白絲氨酸/蘇氨酸激酶活性。

2.10 KEGG功能分類(lèi)

對(duì)西藏山溪鯢所有的 unigenes基因進(jìn)行KEGG通路注釋?zhuān)灿?0 607條 unigenes 得到注釋?zhuān)@些注釋到的unigenes涉及到生物系統(tǒng)、代謝、遺傳信息處理、環(huán)境信息處理和細(xì)胞過(guò)程5個(gè)大類(lèi)30亞類(lèi)的通路信息。如圖4 所示，這些注釋到unigenes 分布于257個(gè)已知功能的代謝通路中，其中有較多的unigenes涉及信號(hào)轉(zhuǎn)導(dǎo)通路，共有1 058條；其他的幾個(gè)分別為MAPK 信號(hào)通路(320條)、Wnt信號(hào)通路(191條)、Calcium信號(hào)通路(190條)、ErbB信號(hào)通路(130條)和TGF-beta信號(hào)通路(118條)，這些代謝通路都與環(huán)境信息處理大類(lèi)中信號(hào)轉(zhuǎn)導(dǎo)相關(guān)。注釋到細(xì)胞通訊通路的unigenes有862條與細(xì)胞過(guò)程有關(guān)，占第二位；其次是注釋到與免疫系統(tǒng)通路相關(guān)的unigenes有757條；推測(cè)這與西藏山溪鯢在遭受到病原微生物入侵中產(chǎn)生的免疫應(yīng)答過(guò)程有著重要功能，同時(shí)也富集到與環(huán)境適應(yīng)相關(guān)的unigenes，為28條，這可能與西藏山溪鯢適應(yīng)高山生活環(huán)境的特殊性有一定關(guān)系。

3 討論

目前高通量測(cè)序技術(shù)在缺少基因組信息的非模式生物研究中已被廣泛應(yīng)用[20-23]。當(dāng)前尚無(wú)西藏山溪鯢轉(zhuǎn)錄組研究工作的報(bào)道。本研究中，利用這種測(cè)序技術(shù)對(duì)西藏山溪鯢組織進(jìn)行了轉(zhuǎn)錄組測(cè)序，在沒(méi)有參考基因組的情況下對(duì)其進(jìn)行了拼裝，共產(chǎn)生了36 252條unigenes，檢測(cè)到的unigenes平均長(zhǎng)度為937 bp，其中長(zhǎng)度≥2 000 bp以上的有4 123條。表明本次測(cè)序質(zhì)量較高，數(shù)據(jù)組裝效果很好；也說(shuō)明高通量測(cè)序是一種可靠性較好、能高效獲取非模式生物基因序列的方法。對(duì)所有獲得的unigenes在6大數(shù)據(jù)庫(kù)中進(jìn)行比對(duì)，結(jié)果都得到了注釋。但在物種注釋上，注釋結(jié)果相似度最高的是西部錦龜，也僅為13.6%；仍有10 562條(55.8%)unigenes注釋的物種不明確，一方面原因可能是基因數(shù)據(jù)庫(kù)中山溪鯢屬物種基因資源偏少，影響功能注釋的基因；另一方面可能有些unigenes是西藏山溪鯢特有的基因，后續(xù)需要進(jìn)一步的研究。

在GO注釋的unigenes中，與轉(zhuǎn)錄、DNA依賴性、細(xì)胞核和ATP結(jié)合相關(guān)的基因最多，這可能與西藏山溪鯢組織的生長(zhǎng)、細(xì)胞的增殖與分化和能量代謝密切相關(guān)。在KOG注釋的unigenes中，獲得了西藏山溪鯢轉(zhuǎn)錄組數(shù)據(jù)庫(kù)，得到了與西藏山溪鯢生長(zhǎng)發(fā)育、生物合成與代謝相關(guān)基因資源。此外，從KEGG通路分類(lèi)結(jié)果看：共有10 607條西藏山溪鯢unigenes參與到257個(gè)已知功能的代謝通路中。其中參與西藏山溪鯢信號(hào)通路的unigenes最多，為1 058條，其中最多的是MAPK 信號(hào)通路和Wnt信號(hào)通路的unigenes數(shù)目很多，分別為320條和191條，說(shuō)明這些信號(hào)分子在西藏山溪鯢生命活動(dòng)與代謝活動(dòng)中起著重要的生理作用。其次是參與細(xì)胞通訊通路的unigenes有862條，這可能與西藏山溪鯢在生長(zhǎng)過(guò)程中要不斷適應(yīng)自身環(huán)境有關(guān)。其中unigenes被注釋到免疫系統(tǒng)通路中占第三，有757條。這些涉及免疫通路相關(guān)的基因主要包括Toll樣受體、T細(xì)胞抗原活化分子、干擾素刺激基因模式識(shí)別受體、補(bǔ)體成分和抗菌肽基因等，表明在西藏山溪鯢生長(zhǎng)中可能形成了其特有的天然免疫機(jī)制。本研究還發(fā)現(xiàn)有48條unigenes注釋到環(huán)境適應(yīng)，推測(cè)可能與西藏山溪鯢特定的低溫棲息環(huán)境和適應(yīng)高海拔生活特點(diǎn)有關(guān)。

此外，高通量測(cè)序技術(shù)的另一個(gè)優(yōu)勢(shì)是能快速地從大量基因序列中獲得SSR分子標(biāo)記資源，能被廣泛用于動(dòng)植物的進(jìn)化論和遺傳學(xué)研究[20]。本研究中利用MISA軟件查找測(cè)序的數(shù)據(jù)，檢測(cè)到SSR總數(shù)為2 779條，其中單堿基型的數(shù)量最多，為1 473條，所占比例超過(guò)50%。這與在其他水生物種轉(zhuǎn)錄組報(bào)道有相似的結(jié)果。例如，岳華梅等[24]利用該測(cè)序技術(shù)對(duì)興國(guó)紅鯉 (Cyprinuscarpiovar.singuonensis)進(jìn)行了SSR標(biāo)記篩選，發(fā)現(xiàn)單堿基型占的比例最大，為47.86%。Zhou X X等[25]對(duì)刺參(Apostichopusjaponicus)進(jìn)行轉(zhuǎn)錄組序列分析，發(fā)現(xiàn)單堿基型有9 154條，所占比例為75.56%。Huang Y等[26]對(duì)大鯢(Andriasdavidianus)進(jìn)行轉(zhuǎn)錄組SSR研究，也發(fā)現(xiàn)單堿基型有25 100條，所占比例達(dá)到了84.3%。表明單堿基型的SSR分子標(biāo)記類(lèi)型可能普遍存在于水生動(dòng)物中。這些數(shù)據(jù)的獲得，極大地豐富了西藏山溪鯢轉(zhuǎn)錄本信息和基因資源，可為西藏山溪鯢相關(guān)性狀的基因進(jìn)行深入定位與克隆，并為后續(xù)同屬物種群體遺傳多樣性分析與連鎖圖譜構(gòu)建、分子標(biāo)記開(kāi)發(fā)、評(píng)估和保護(hù)其遺傳資源、適應(yīng)性進(jìn)化機(jī)制等研究提供分子基礎(chǔ)支撐。