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

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

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

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

1 材料與方法

1.1 材料

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

1.2 樣品收集

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

1.3 組織總RNA檢測

將混合的西藏山溪鯢組織樣品迅速置于裝有液氮的研缽中，研磨成粉狀，根據(jù)Takara公司提供的Trizol Reagent操作說明書完成總RNA的抽提。獲得的總RNA 再經(jīng)過1.2%瓊脂糖凝膠電泳和Nanodrop-2000核酸蛋白測定儀檢測總RNA的完整性和純度。經(jīng)檢驗合格的總RNA樣品再進行后續(xù)的轉(zhuǎ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)錄組測序與分析

樣品檢測合格后，利用 Oligo (dT)磁珠法純化出mRNA，然后將mRNA截為片段，經(jīng)過PCR方法得到西藏山溪鯢cDNA文庫，最后建好的cDNA文庫利用Illumina HiSeq 4000平臺技術(shù)進行上機測序，由杭州聯(lián)川生物技術(shù)股份有限公司完成測序。原始 reads 經(jīng)過去接頭并且過濾掉低質(zhì)量及長度過短序列后，得到高質(zhì)量測序數(shù)據(jù)(clean reads)。利用 Trinity 軟件的Paired-end拼接方法對clean reads進行Denovo組裝，得到unigenes。由于西藏山溪鯢目前沒有基因組數(shù)據(jù)，本研究以報道的墨西哥蠑螈(Ambystomamexicanum)(https://www.ncbi.nlm.nih.gov/data-hub/genome/GCA_002915635.3/)基因組序列作參考基因組，使用STAR 軟件進行比對[16]，再用StringTie 軟件基于參考基因組注釋文件對所有轉(zhuǎn)錄本進行整合組裝[17]。使用MiscroSAtellite(MISA)軟件進行SSR鑒定和分析。應用Rsem軟件進行基因表達定量分析，基因表達量采用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個公共數(shù)據(jù)庫與組裝得到的西藏山溪鯢unigenes進行Blast序列比對[18]。選擇閥值條件為E value <1e-10，進行功能注釋。

2 結(jié)果與分析

2.1 RNA檢測

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

2.2 測序質(zhì)量控制

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

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

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

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

2.4 SNP檢測分析

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

表2 SNP位點的分類

2.5 CDS預測與微衛(wèi)星SSR位點分析

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

表3 SSR位點統(tǒng)計

2.6 基因表達量分析

在轉(zhuǎn)錄組測序中，用RPKM方法計算基因的表達水平(表4)，表示在每百萬reads中來自某一基因每千堿基長度的read讀數(shù)。本研究以RPKM≥0.1作為基因表達標準，在已獲得的36 252個unigenes中，有基因表達量的序列為36 226條。對不同RPKM 區(qū)間的基因數(shù)量進行統(tǒng)計發(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條，高表達的基因RPKM值>60為1 243條；而RPKM值在0.1～0.3的低表達基因最少，僅16條，占所有表達基因的0.04%。上述結(jié)果表明Illumina HiSeq測序技術(shù)能夠檢測到極低水平基因的表達。

表4 FPKM值密度分布

2.7 Unigenes功能注釋

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

表5 序列功能注釋

2.8 KOG功能分類

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

2.9 GO功能注釋

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

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

注：1.轉(zhuǎn)錄、DNA依賴性;2.轉(zhuǎn)錄調(diào)控、DNA依賴性;3.蛋白質(zhì)轉(zhuǎn)運;4.多細胞器官發(fā)育;5.細胞分化;6.凋亡;7.細胞黏附;8.蛋白質(zhì)水解;9.信號轉(zhuǎn)導;10.細胞分化;11.細胞周期;12.有絲分裂;13.小G蛋白介導信號轉(zhuǎn)導;14.DNA 修復;15.轉(zhuǎn)運; 16.RNA加工;17.翻譯;18.轉(zhuǎn)錄正調(diào)控;19.胞內(nèi)信號轉(zhuǎn)導;20.染色質(zhì)修飾;21.轉(zhuǎn)錄負調(diào)控;22.細胞內(nèi)蛋白質(zhì)轉(zhuǎn)運;23.跨膜輸送;24.精子形成;25.RNA剪切;26.細胞核;27.必需膜;28.細胞漿;29.胞液;30.細胞質(zhì)膜;31.內(nèi)質(zhì)網(wǎng)膜;32.線粒體;33 胞外區(qū);34.胞核;35.核質(zhì);36.細胞骨架;37.高爾基體膜;38.微管;39.高爾基氏復合體;40.細胞膜;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功能分類

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

3 討論

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

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

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