西藏山溪鯢轉錄組組裝與分析

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

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

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

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

1 材料與方法

1.1 材料

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

1.2 樣品收集

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

1.3 組織總RNA檢測

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

1.4 轉錄組測序與分析

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

2 結果與分析

2.1 RNA檢測

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

2.2 測序質量控制

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

2.3 轉錄組數據組裝與分析

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

表1 測序組裝的序列結果統計

2.4 SNP檢測分析

單核苷酸多態性定義為基因組上單個核苷酸產生的變異，是進行分子標記鑒定、輔助育種和遺傳圖譜構建等非常重要的一種遺傳標記方法。本研究利用轉錄組作為參考序列，使用BWA和Samtools軟件對西藏山溪鯢的外顯子區域進行SNP發掘，結果顯示總共得到 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預測與微衛星SSR位點分析

通過Genscan軟件預測其CDS序列，預測結果顯示：總共檢測出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類型單核苷酸出現頻率最多，有1 473個；其次為二核苷酸，有715個；490個SSR具有三核苷酸位點；四和五核苷酸SSR位點的數量分別為56個和38個；SSR位點為六核苷酸的數量最少，僅為7個。這些SSR位點可為后續進行分子標記鑒定提供引物設計基礎。

表3 SSR位點統計

2.6 基因表達量分析

在轉錄組測序中，用RPKM方法計算基因的表達水平(表4)，表示在每百萬reads中來自某一基因每千堿基長度的read讀數。本研究以RPKM≥0.1作為基因表達標準，在已獲得的36 252個unigenes中，有基因表達量的序列為36 226條。對不同RPKM 區間的基因數量進行統計發現，其中RPKM值在3.57～15.00的基因最多，為17 099條，占到了47.2%；其次是RPKM值位于0.30～3.57的基因，數量為13 868條；RPKM值位于15～60的基因有4 000條，高表達的基因RPKM值>60為1 243條；而RPKM值在0.1～0.3的低表達基因最少，僅16條，占所有表達基因的0.04%。上述結果表明Illumina HiSeq測序技術能夠檢測到極低水平基因的表達。

表4 FPKM值密度分布

2.7 Unigenes功能注釋

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

表5 序列功能注釋

2.8 KOG功能分類

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

2.9 GO功能注釋

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

注：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.細胞骨架。

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

2.10 KEGG功能分類

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

3 討論

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

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

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