東海帶魚(Trichiurus japanicus)肝臟轉(zhuǎn)錄組SSR和SNP特征分析

章 霞，柳敏海，李凌剛，徐志進，李偉業(yè)，殷小龍，傅榮兵

(浙江省舟山市水產(chǎn)研究所，浙江 舟山 316000)

東海帶魚(Trichiurusjapanicus)，俗稱刀魚、白帶魚和鱗刀魚，屬于亞熱帶性魚類，主要分布于我國的黃渤海、東海以及南海海區(qū)[1-2]，曾為我國重要的經(jīng)濟捕撈魚類，年產(chǎn)量最高曾達50余萬噸，約占世界同種魚漁獲量的70%[3-4]。帶魚是目前東?！八拇蠛．a(chǎn)”中唯一能形成較大漁汛的傳統(tǒng)捕撈對象，在東海漁業(yè)中具有舉足輕重的位置。近年來，由于過度的捕撈，東海帶魚漁業(yè)資源已經(jīng)呈衰退趨勢[5]，捕撈群體出現(xiàn)了小型化、低齡化、產(chǎn)量下降等現(xiàn)象，可見東海帶魚資源保護刻不容緩。目前，國內(nèi)外對帶魚的研究主要集中于其生物學特征[6]、年齡與生長[7]、攝食習性[1，8]、漁業(yè)資源[4]、種群結(jié)構(gòu)分析[9]等方面，關(guān)于帶魚的分子標記開發(fā)、基因克隆和組學研究等現(xiàn)代分子生物實驗技術(shù)應(yīng)用于資源保護方面的研究寥寥。

表達序列標簽SSR(Expressed sequence tag SSR，EST-SSR)和單核苷酸多態(tài)性(Single nucleotide polymorphism，SNP)是開展物種遺傳多樣性、分子標記開發(fā)、目標性狀定位的重要技術(shù)手段。EST-SSR是指從cDNA文庫或轉(zhuǎn)錄組數(shù)據(jù)中獲得的與功能基因直接相關(guān)的SSR位點[10-11]，SNP是指個體在基因組水平上發(fā)生單個核苷酸的變異，從而引起的DNA序列多態(tài)性。目前已有許多從水生動物轉(zhuǎn)錄組數(shù)據(jù)庫中篩選SSR、SNP標記的成功案例，如基于轉(zhuǎn)錄組測序曼氏無針烏賊(Sepiellamaindroni)[12]和凡納濱對蝦(Litopenaeusvannamei)[13]SNP的開發(fā)，大口黑鱸(Micropterussalmonides)馴食相關(guān)[14]、三疣梭子蟹(Portunustrituberculatus)生長相關(guān)SNP位點的鑒定[15]等。但關(guān)于東海帶魚EST-SSR和SNP的研究未見報道。

本研究從已獲得轉(zhuǎn)錄組數(shù)據(jù)(85 710個Unigene序列)中，利用軟件篩選獲得SSR和SNP，并對其進行統(tǒng)計分析，為后續(xù)開發(fā)東海帶魚的SSR和SNP標記提供科學數(shù)據(jù)。

1 材料與方法

1.1 試驗材料與數(shù)據(jù)來源

東海帶魚樣本于2017年在舟山市六橫海域獲得，樣本數(shù)為6尾[體質(zhì)量為(42.34±3.16)g]，肛長為[(15.37±1.46)cm]，標記為T1、T2、T3、T4、T5、T6。取其肝臟送上海凌恩生物科技有限公司，通過Illumina HiSeq(儀器型號Hiseq X-10)高通量測序技術(shù)(Illumina Xten PE150)獲得轉(zhuǎn)錄組數(shù)據(jù)(RNA seq，結(jié)果未發(fā)表)，數(shù)據(jù)量統(tǒng)計見下表1。

表1 東海帶魚轉(zhuǎn)錄組數(shù)據(jù)量統(tǒng)計

1.2 轉(zhuǎn)錄組SSR的篩選及分析

利用MISA軟件(MIcroSAtellite identification tool)從東海帶魚的獨立基因序列中進行SSR搜索，篩選標準：二核苷酸重復(fù)數(shù)大于等于6，三核苷酸、四核苷酸、五核苷酸、六核苷酸大于等于5；運用Excel軟件對轉(zhuǎn)錄組的SSR各類型比例、序列分布和特征進行分析和統(tǒng)計。

SSR發(fā)生頻率=含SSR的Unigene數(shù)/Unigene總數(shù)；

SSR分布的平均距離=總Unigene長度/搜索到的SSR數(shù)量。

1.3 轉(zhuǎn)錄組SNP的篩選及分析

以組裝好的轉(zhuǎn)錄本為模板序列，將質(zhì)控后的高質(zhì)量序列使用BWA(BWA mem-k 32)與其進行比對。使用Samtools對比對的sam文件進行格式轉(zhuǎn)換(Samtools view-bS-t ref_genome.fa.fai xx.sam>xx.bam)和排序(samtools sort-m 5920000000 xx.bam xx.sort)。對高質(zhì)量的bam文件使用GATK“UnifiedGenotyper”功能(java-Xmx20G-jar GenomeAnalysisTK-2.7-2/GenomeAnalysisTK.jar-T UnifiedGenotyper-R ref_genome.fa-I xx.sort.bam-metrics xx.SNP.metrics-o xx.SNP.vcf-stand_call_conf 50-stand_emit_conf 50-dcov 2000-glm SNP)進行SNP和indel鑒定，并運用Excel軟件進行相應(yīng)的數(shù)據(jù)統(tǒng)計。

2 結(jié)果與分析

2.1 東海帶魚肝臟轉(zhuǎn)錄組序列SSR的特征分析

2.1.1 SSR重復(fù)基元的分布

經(jīng)數(shù)據(jù)組裝、去除冗余和進一步拼接后，最終獲得85 710個Unigene序列，總長度68 387 598 bp?；诖私Y(jié)果上進行SSR檢測分析，結(jié)果發(fā)現(xiàn)，能識別SSR總數(shù)為49 311個，分布在35 098條Unigene序列中。其中10 162條Unigene含有一個以上SSR位點，復(fù)合型SSR數(shù)目為4 687個。SSR的發(fā)生頻率為40.95%，平均1.33 kb出現(xiàn)1個SSR位點。

2.1.2 東海帶魚轉(zhuǎn)錄組序列中SSR的數(shù)量和分布特點

各堿基類型中的不同序列組成的微衛(wèi)星序列占比差異較大。由數(shù)據(jù)得，東海帶魚中共有68種重復(fù)基元。二、三、四、五、六堿基重復(fù)基元分別有4、10、26、17、9種。單核苷酸重復(fù)類型含量約占所有堿基重復(fù)類型的50.79%，比例最高，其中(A/T)n出現(xiàn)頻率最高，占單核苷酸類型的95.53%；二堿基重復(fù)SSR含量約占總數(shù)的30.63%，其中出現(xiàn)的頻率最多的是(AC/GT)n，占比約為72.16%；三堿基重復(fù)SSR占所有堿基重復(fù)類型的14.34%，四堿基重復(fù)SSR占比1.12%，五堿基、六堿基重復(fù)SSR所占比例較少。(AGG/CCT)n、(AAAT/ATTT)n、(AAAAT/ATTTT)n、(ACCAGG/CCTGGT)n在三、四、五、六堿基重復(fù)基元類型出現(xiàn)頻率最多，在各自類型中占比分別是31.39%、16.64%、20.00%、31.25%(表2)。在所有的堿基重復(fù)模式中，各種重復(fù)基元中總SSR比例前五的依次為(AC/GT)n(22.10%)、(AG/CT)n(5.56%)、(AGG/CCT)n(4.5%)、(AGC/CTG)n(3.53%)、(AT/AT)n(2.87%)，其他不同類型重復(fù)基元SSR占總SSR的比例分布見圖1。

注：“其他”表示頻率小于 0．50%的重復(fù)基元類型。

Note:“Other motifs” denoted the repeat motifs with frequency below 0.50%.

2.1.3 SSR長度分析

東海帶魚肝臟轉(zhuǎn)錄組中共發(fā)現(xiàn)有44 625個完整的SSR，bp片段長度由11至374不等，平均長度24.03 bp。在SSR序列中，主要是重復(fù)長度小于20 bp的序列，有30 997條，占總數(shù)的69.46%，其中在11～15 bp的數(shù)量最多，16～20 bp次之，具體見圖2。SSR數(shù)量與長度(重復(fù)次數(shù))的關(guān)系如圖3所示，各種SSR的數(shù)量隨著重復(fù)次數(shù)的增加而減少。單堿基重復(fù)次數(shù)達到23次時，SSR的數(shù)量突然增加，之后繼續(xù)下降，二堿基重復(fù)達到11次時也出現(xiàn)類似的情況。三堿基重復(fù)數(shù)量隨重復(fù)次數(shù)的增加而減少，其他SSR(四、五、六堿基重復(fù))的數(shù)量呈此變化規(guī)律。在各個堿基類型重復(fù)次數(shù)達到一定數(shù)值時，各自的下降速率變緩慢，其中單堿基、二堿基、三堿基、其他堿基的重復(fù)次數(shù)點分別為25次、13次、9次和7次，具體見圖3。

注：“其他”表示四、五、六堿基重復(fù)SSR的總和。

Note:“Others” denoted the total of trinucleotide，tetranucleotide，hexnucleotide repeats.

2.2 東海帶魚肝臟轉(zhuǎn)錄組序列SNP的特征分析

2.2.1 SNP顛換轉(zhuǎn)、換信息

在所獲得的轉(zhuǎn)錄組測序結(jié)果中，通過利用Samtools和GATK軟件共篩選含有SNP位點的Unigenes為67 376條，占總轉(zhuǎn)錄組測序獲得的Unigenes數(shù)的78.61%。這些序列上獲得SNP位點數(shù)為876 199個，平均每76.8 bp含有1個SNP。其中顛換335 218個，轉(zhuǎn)換540 981個，轉(zhuǎn)換約為顛換的1.61倍。轉(zhuǎn)換類型中C：G→T：A的發(fā)生頻率高于T：A→C：G；顛換類型中T：A→G：C比例最高，占比10.50%(表3)。

表3 東海帶魚SNP顛換、轉(zhuǎn)換信息

2.2.2 SNP堿基轉(zhuǎn)換類型分析

對東海帶魚肝臟轉(zhuǎn)錄組序列堿基轉(zhuǎn)換類型進行分析，在非編碼區(qū)的SNP數(shù)量為578 020；在編碼區(qū)的SNP總數(shù)298 179，其中同義SNP 為203 863個，非同義SNP為93 024個，stopgain SNP為1 189， stoploss SNP為103。東海帶魚在編碼區(qū)的堿基轉(zhuǎn)換類型中C→T占比最高，為22.09%，其次為G→A類型，占比21.91%，T→A占比最低，為2.92%，具體見表4。

表4 東海帶魚SNP在編碼區(qū)的堿基轉(zhuǎn)換類型數(shù)據(jù)分析

2.3 東海帶魚肝臟轉(zhuǎn)錄組序列SSR和SNP多態(tài)性分析

根據(jù)東海帶魚肝臟轉(zhuǎn)錄組序列結(jié)果，本研究提供各原始引物各10對，可用于今后的帶魚多態(tài)性分析和驗證。引物信息見表5。

表5 東海帶魚肝臟轉(zhuǎn)錄組SSR和SSR部分引物信息

3 討論

新一代的測序技術(shù)出現(xiàn)后，通過轉(zhuǎn)錄組測序篩選獲得SSR和SNP已成為開發(fā)應(yīng)用于遺傳多樣性分析、分子育種等分子標記的重要技術(shù)手段，此技術(shù)相較于基因組篩序獲得分子標記，不僅能大大降低成本，且能快速精準地獲得具有較高通用性的有用序列，使得開發(fā)利用分子標記更加有效便捷[16-17]。

本研究中在東海帶魚肝臟轉(zhuǎn)錄組Unigene序列中可識別SSR總數(shù)為49 311個，發(fā)生頻率為40.95%，平均1.33 kb出現(xiàn)1個SSR位點，相較于其他水生動物SSR發(fā)生頻率較高。據(jù)報道，大竹蟶(Solenstrictus)[18]的SSR發(fā)生頻率為38.30%，黃姑魚(Nibeaalbiflora)[19]為39.30%，凡納濱對蝦[14]為16.76%，馬氏珠母貝(Pinctadamartensi)[20]為13.34%，縊蟶(Sinonovaculaconstricta)[21]為8.89%。SSR發(fā)生頻率的巨大差異可能與不同物種、組織、測序平臺的差異性有關(guān)[22]。另本研究結(jié)果表明帶魚肝臟轉(zhuǎn)錄組中的SSR數(shù)量較為豐富，總數(shù)為49 311個，共有68種重復(fù)基元。除單核苷酸重復(fù)類型外，二堿基重復(fù)類型數(shù)量最多，約占總數(shù)30.63%；三堿基重復(fù)類型次之，占比為14.34%；且呈現(xiàn)隨著堿基數(shù)量增加，微衛(wèi)星逐步減少的現(xiàn)象，符合眾多真核生物的轉(zhuǎn)錄組微衛(wèi)星組成情況[23-24]。同時，二堿基重復(fù)SSR中的GC/CG核心基元含量僅為二堿基重復(fù)SSR總數(shù)的0.1%，推測GC/CG重復(fù)SSR在所有物種的基因組中都極為罕見[11]，但不排除是由于物種不同和微衛(wèi)星搜索條件設(shè)定沒有統(tǒng)一的標準，造成以上結(jié)果。

具有多態(tài)性的SSR位點才可能開發(fā)為微衛(wèi)星標記，而多態(tài)性高低則決定SSR的可用性。一般動植物的多態(tài)性與其長度呈正相關(guān)[25-26]。Tenmykh等[27]認為，SSR長度≥20 bp，微衛(wèi)星多態(tài)性較高，12 bp

本研究中篩選獲得SNP位點的Unigenes為67 376條，占總轉(zhuǎn)錄組測序獲得的Unigenes數(shù)的78.61%。這些序列上獲得SNP 位點數(shù)為876 199個，平均每76.8 bp含有1個SNP。其中顛換335 218個，轉(zhuǎn)換540 981個，理論上發(fā)生轉(zhuǎn)換的概率與發(fā)生顛換概率的比值應(yīng)該等于0.5 (1∶2)，但有些生物的比值常常>0.5，這種差異被稱為“轉(zhuǎn)換偏差”[28]。在本文中發(fā)現(xiàn)轉(zhuǎn)換概率約為顛換的1.61倍。大于理論值，說明存在轉(zhuǎn)換偏差，這可能與堿基組成及進化過程的選擇機制有關(guān)，表明堿基轉(zhuǎn)換突變并不是隨機產(chǎn)生[28]。另在研究中發(fā)現(xiàn)C/T發(fā)生頻率在轉(zhuǎn)換類型中最高，這可能與在CG序列上出現(xiàn)的高頻率，胞嘧啶極易甲基化脫去氨基而形成胸腺嘧啶等因素有關(guān)[29]。

對東海帶魚轉(zhuǎn)錄組序列堿基轉(zhuǎn)換類型進行分析，在非編碼區(qū)的SNP數(shù)量為578 020；在編碼區(qū)的SNP總數(shù)298 179，其中同義SNP 為203 863個，非同義SNP為93 024個。編碼區(qū)內(nèi)的同義SNP造成的編碼序列的變化不會引起氨基酸序列變化；而非同義SNP則會影響蛋白質(zhì)序列，導(dǎo)致生物性狀改變[30]。因此，開發(fā)東海帶魚編碼區(qū)非同義SNP(93 024個)具有重要的生物學意義。

綜上所述，本研究基于東海帶魚幼魚的肝臟轉(zhuǎn)錄組測序的基礎(chǔ)上，進行了SSR和SNP的檢索分析，并對SSR重復(fù)序列特征和SNP多態(tài)性特征進行了研究，補充和完善了帶魚基因數(shù)據(jù)庫信息，今后可開展東海帶魚轉(zhuǎn)錄組序列GO分類、COG分類注釋和代謝通路注釋(KEGG pathway)，EST-SSR、EST-SNP相關(guān)性狀驗證分析等工作，為帶魚的高密度遺傳連鎖圖譜構(gòu)建、物種鑒定以及種質(zhì)資源保護等提供重要科學支撐。