基于RNA-Seq 技術(shù)的鯪轉(zhuǎn)錄組分析

許建，趙建，徐禮鳴，崔軍，李強(qiáng)，朱新平，徐鵬

(1.中國水產(chǎn)科學(xué)研究院 生物技術(shù)研究中心，北京100041;2.中國水產(chǎn)科學(xué)研究院 珠江水產(chǎn)研究所，廣東 廣州510380)

鯪Cirrhina molitorella 俗稱土鯪、鯪公、花鯪，是中國珠江流域地區(qū)的特有種，華南地區(qū)“四大家魚”(鰱Hypophthalmichthy smolitrix、鳙H.snobilis、草Ctenopharyngodon idella、鯪)之一，僅在廣東省年養(yǎng)殖產(chǎn)量就在20 萬t 左右。鯪肉質(zhì)細(xì)嫩、味鮮美、產(chǎn)量大、價格適中，是市場的暢銷水產(chǎn)品。鯪也可入藥，具有健筋骨、活血行氣、逐水利溫之功效。目前，生長速度慢和不耐寒等問題是鯪產(chǎn)業(yè)發(fā)展的瓶頸，快速生長的鯪品系不僅會提高傳統(tǒng)養(yǎng)殖地區(qū)漁民的養(yǎng)殖積極性，也可以推廣到生長期較短的北方地區(qū)，能極大地提高鯪在中國淡水漁業(yè)中的地位。所以，選育快速生長的鯪新品系是解決鯪產(chǎn)業(yè)發(fā)展中關(guān)鍵問題的有效途徑。

近年來，珠江水產(chǎn)研究所通過對野生鯪資源及其生長參數(shù)進(jìn)行調(diào)查，在西江群體中篩選出一個早期生長較快的群體，有效地進(jìn)行了鯪的保種和擴(kuò)繁，可作為進(jìn)一步選育的基礎(chǔ)群。然而，傳統(tǒng)家系選育方法歷時長、工作量大，而分子標(biāo)記輔助育種則可以大大節(jié)省時間和勞動力，是一種快速有效的育種措施，但受限于鯪基因組信息較少，迄今為止，僅有少量微衛(wèi)星和零星單核苷酸多態(tài)性(SNP)標(biāo)記在種群遺傳分析中進(jìn)行了初步應(yīng)用［1-3］，遠(yuǎn)遠(yuǎn)不能滿足標(biāo)記輔助育種的要求，亟須找出一種快速、大量地獲得鯪遺傳信息的途徑。此外，對鯪營養(yǎng)、生理生化等方面的研究，也亟須獲取鯪相關(guān)的功能基因，而目前僅僅依賴從近緣模式魚類斑馬魚Danio rerio 基因組獲取相關(guān)序列信息，然后進(jìn)行繁瑣的分子克隆實(shí)驗(yàn)才能獲取相關(guān)基因和序列。因此，盡快建立鯪轉(zhuǎn)錄組數(shù)據(jù)庫具有重要的意義。鑒于此，本研究中開展了鯪的首個高通量轉(zhuǎn)錄組研究，采集代表性鯪種群中多個個體的組織樣本，采用第二代基因組測序技術(shù)進(jìn)行深度的轉(zhuǎn)錄組測序，利用生物信息學(xué)分析流程，系統(tǒng)地進(jìn)行了序列清洗、基因拼接、全長序列獲取、基因功能注釋、重復(fù)序列和元件分析評估、微衛(wèi)星和SNP 位點(diǎn)挖掘等分析研究，并建立了鯪轉(zhuǎn)錄組數(shù)據(jù)庫和網(wǎng)站，旨在方便同行隨時調(diào)取數(shù)據(jù)，實(shí)現(xiàn)數(shù)據(jù)共享。

1 材料與方法

1.1 材料

試驗(yàn)用鯪采自珠江水產(chǎn)研究所實(shí)驗(yàn)基地，取生長良好的鯪幼魚10 尾，體長約為3 cm，體質(zhì)量約為3 g。

1.2 方法

1.2.1 樣品的采集和總RNA 的抽提 將10 尾魚迅速置于液氮預(yù)冷的研缽中，邊加液氮邊研磨至粉末狀，混合后使用Invitrogen 公司的TriZol 試劑進(jìn)行總RNA 抽提。使用安捷倫生物分析儀2100和紫外分光光度儀檢測總RNA 的質(zhì)量和數(shù)量。

1.2.2 cDNA 文庫的構(gòu)建和測序 取100 ng RNA樣品(10 條魚混合樣)，利用羅氏公司的Ovation RNA - Seq 試劑盒(NuGEN Technologies，SanCarlos，CA)合成cDNA，通過末端修復(fù)、連接接頭和純化，獲得鯪樣品的cDNA 文庫。將該cDNA 文庫用454 平臺測序，運(yùn)行通量為1 個run。

1.2.3 de novo 拼接及數(shù)據(jù)分析 采用454 測序技術(shù)，原始的測序數(shù)據(jù)被保存為454 特有的SFF 文件格式，由于沒有可用的參考基因組數(shù)據(jù)，使用Newbler 2.8 的-cdna 模式進(jìn)行de novo 拼接，并使用-vt 參數(shù)去除載體序列，用-vs 參數(shù)去除核糖體RNA 序列，最終獲取初步拼接結(jié)果，并保留拼接后大于100 bp 的contig。

1.2.4 SNP 分析 為了進(jìn)行SNP 鑒定和分析，將Newbler 清洗后的SFF 格式文件通過PERL 腳本轉(zhuǎn)換成fastq 文件，以Newbler 軟件拼接的轉(zhuǎn)錄組結(jié)果作為參考序列，采用BWA和SAMtools 軟件對轉(zhuǎn)錄組SNP 進(jìn)行鑒定，測序深度大于10 且測序質(zhì)量值大于20 作為過濾閾值，其余參數(shù)為默認(rèn)值。

1.2.5 重復(fù)序列分析及微衛(wèi)星鑒定 對側(cè)翼長度大于50 bp 的微衛(wèi)星序列采用Msatfinder 2.0.9程序進(jìn)行微衛(wèi)星鑒定。其中鑒定二、三、四、五、六核苷酸的重復(fù)，閾值分別被設(shè)為8、5、5、5、5。

1.2.6 功能注釋 使用在線生物信息學(xué)分析軟件BlastX 將組裝得到的轉(zhuǎn)錄組數(shù)據(jù)與NCBI 的非冗余蛋白數(shù)據(jù)庫、斑馬魚蛋白數(shù)據(jù)庫和UniProt 蛋白數(shù)據(jù)庫進(jìn)行比對，獲得同源蛋白匹配結(jié)果，e -value閾值設(shè)定為1E - 5。為進(jìn)一步通過NCBI Entrez Gene 數(shù)據(jù)庫和Ensembl 斑馬魚基因組數(shù)據(jù)的注釋，通過同源比對進(jìn)行GO 注釋，并使用WEGO 進(jìn)行生物過程、分子功能和細(xì)胞成分的富集分析。

2 結(jié)果與分析

2.1 轉(zhuǎn)錄組組裝

采用454 測序技術(shù)獲得了鯪轉(zhuǎn)錄組的數(shù)據(jù)。對原始數(shù)據(jù)進(jìn)行統(tǒng)計(jì)，共獲得1 297 479 條reads，總堿基數(shù)為486 586 191 bp。最大的reads 長度為1195 bp，最小的reads 長度為40 bp，平均reads 長度為375 bp。

鯪的轉(zhuǎn)錄組拼接結(jié)果顯示，共得到19 962 條contigs，平均contig 長度為1269 bp，N50 為1509 bp，最大的contig 長度為20 289 bp，最小的contig長度為101 bp。對每條contig 長度統(tǒng)計(jì)相應(yīng)的contig 數(shù)量，結(jié)果如圖1所示，contig 長度主要分布在600 ～800 bp。

圖1 鯪轉(zhuǎn)錄組contig 長度的分布圖Fig.1 Length freguency of assembled transcriptome contigs in mud carp

2.2 功能注釋及GO 分析

將組裝的contig 與NCBI 的非冗余蛋白庫(NR)、UniProt 蛋白庫和Ensembl 斑馬魚蛋白庫進(jìn)行BlastX 比對。結(jié)果顯示，總共有13 657 條contigs在NR 蛋白數(shù)據(jù)庫中被比對出，其中包括10 577個特異蛋白(表1)。

表1 NR、Uniprot和Zebrafish 數(shù)據(jù)庫的BlastX 比對結(jié)果Tab.1 BlastX search results against NR，Uniprot and Zebrafish database

根據(jù)特異蛋白注釋結(jié)果進(jìn)行GO 分析，結(jié)果有7314 條contigs 有GO 注釋，包含5381 個特異蛋白。將注釋信息整理成WEGO 所需的輸入文件，共分為分子功能(molecular function)、生物途徑(biological process)和細(xì)胞成分(cellular component)3 個大類(圖2)。在生物過程中，與細(xì)胞過程(cellular process)(GO:0009987，GO 注釋條目2694，GO 注釋條目占總條目的比例為50.1%)和代謝過程(metabolic processes)(GO:0008152，GO注釋條目2133，39.6%)相關(guān)的基因產(chǎn)生了顯著的富集。對于分子功能，連接(binding)(GO:0005488，GO 注釋條目2716，50.5%)是最為主要的成分，其次是催化活性(catalytic activity)(GO:0003824，GO 注釋條目1793，33.3%);而細(xì)胞(cell)(GO:0005623，GO 注釋條目2675，49.7%)和細(xì)胞組分(cell part)(GO:0043226，GO 注釋條目2675，49.7%)是最具代表性的類別的細(xì)胞成分。

圖2 鯪轉(zhuǎn)錄組的GO 功能二級分類Fig.2 Gene ontology(GO)(level 2)for transcriptome in mud carp Cirrhina molitorella

2.3 微衛(wèi)星預(yù)測

從2058 條contigs 中，總共鑒定出2438 個微衛(wèi)星。這些微衛(wèi)星包括二堿基、三堿基、四堿基、五堿基和六堿基重復(fù)(閾值分別被設(shè)定為8、5、5、5、5)。以側(cè)翼序列50 bp 為閾值篩選，獲得側(cè)翼序列大于50 bp 的微衛(wèi)星共1379 個，并為后續(xù)的PCR 驗(yàn)證設(shè)計(jì)了相關(guān)的引物(表2)。

表2 鯪轉(zhuǎn)錄組序列中微衛(wèi)星的分布情況Tab.2 Statistics of microsatellites of transcriptome from mud carp Cirrhina molitorella

2.4 SNP 分析

利用轉(zhuǎn)錄組作為參考序列，使用BWA和SAMtools 程序?qū)N的外顯子區(qū)域進(jìn)行SNP 發(fā)掘，結(jié)果顯示，總共得到5014 個SNP 位點(diǎn)，包括A -G、C-T、A-C、G-T、A -T和G -C 六種類型的SNP。在所有類型的SNP 中，A -G和C -T 類型的比例最高，占所有SNP 位點(diǎn)的72.6%，A -C、G-T、A-T和G -C 這4種SNP 類型有相似的比例，占總量的27.4%(表3)。

表3 鯪轉(zhuǎn)錄組SNP 位點(diǎn)的分類Tab.3 Classification of SNPs identified transcriptome from mud carp Cirrhina molitorella

2.5 全長cDNA 的鑒定

根據(jù)Ensembl 斑馬魚蛋白數(shù)據(jù)庫，通過Target-Identifier在線工具對所有的contigs 進(jìn)行完整ORF查找。結(jié)果共鑒定出5931 個全長cDNA 序列(e -value 為1E-5)(表4)，全長cDNA 的長度分布如圖3所示。

2.6 鯪和鯉轉(zhuǎn)錄組的比較

在鯉科魚類中，鯉的全魚轉(zhuǎn)錄組已經(jīng)被測定［4］。通過比較，鯪的總reads 數(shù)要少于鯉，但平均reads 長度比鯉長。鯉的contig 條數(shù)也要比鯪多，鯪的N50 為1509 bp，這表明Newbler 軟件用于轉(zhuǎn)錄組組裝的分析結(jié)果比較理想。利用NCBI 的非冗余蛋白庫進(jìn)行比對，結(jié)果顯示，13 657 條contigs得到注釋，占總contigs 數(shù)的68.4%，鯉為76.2%;此外，兩者具有GO 注釋的特異蛋白的比例也基本一致，分別為50.9%和49.8%(表5)。

表4 Target Identifier 分析結(jié)果Tab.4 Summary of the Target Identifier report

圖3 鯪轉(zhuǎn)錄組全長cDNA 的長度分布圖Fig.3 Length distributions of putative full - length cDNAs in mud carp transcriptome

2.7 構(gòu)建鯪轉(zhuǎn)錄組數(shù)據(jù)庫

為了更便捷地應(yīng)用鯪轉(zhuǎn)錄組的分析結(jié)果，根據(jù)轉(zhuǎn)錄組裝配序列、基因注釋結(jié)果、蛋白家族、SNP、微衛(wèi)星等信息構(gòu)建了鯪轉(zhuǎn)錄組數(shù)據(jù)庫和網(wǎng)站(http://genomics.cafs.ac.cn/atd_www/)。通 過 該數(shù)據(jù)庫，可以進(jìn)行鯪基因的序列調(diào)取、序列比對、斑馬魚同源序列下載等，為鯪的基因克隆等基礎(chǔ)遺傳研究和分子育種提供了有力的支持。

3 討論

目前，基于RNA -Seq 的二代測序技術(shù)已經(jīng)改變了轉(zhuǎn)錄組的研究方式［5-6］，RNA 測序在揭示轉(zhuǎn)錄組的復(fù)雜性、基因的鑒定、相關(guān)標(biāo)記的分析、非編碼RNA 分析和可變剪切分析等方面起著重要的作用［7-9］。RNA 測序技術(shù)在魚類的研究中也發(fā)揮了重要作用。利用該技術(shù)在斑馬魚早期發(fā)育的不同階段，鑒定出差異基因，為斑馬魚的發(fā)育提供了理論基礎(chǔ)［10］。在進(jìn)化研究方面，Zheng等［11］利用該技術(shù)，通過斑馬魚鰾的轉(zhuǎn)錄組同人肺的比較，發(fā)現(xiàn)了魚類鰾和人肺是同源器官的分子生物學(xué)證據(jù);Shin等［12］通過對3種南極魚轉(zhuǎn)錄組的研究鑒定出一些抗寒基因;Liu等［7］對斑點(diǎn)叉尾鮰Ictalurus furcatus升溫至其失去平衡的試驗(yàn)中鑒定出魚適應(yīng)高溫脅迫的相關(guān)基因及生理過程;草魚經(jīng)GCRV 感染后，在頭腎中鑒定出很多與免疫相關(guān)的基因和代謝途徑，為培育出抗GCRV 的品種提供了理論基礎(chǔ)［13］。在性別差異方面，Sun等［14］利用該技術(shù)，通過比較斑點(diǎn)叉尾鮰雌雄性腺的基因表達(dá)差異，鑒定出影響性別的決定基因，為研究魚的性別決定及分化提供了可靠的資料。另外，將RNA-Seq 技術(shù)用于大量的SNP 鑒定，為遺傳育種、種群遺傳學(xué)和資源評估等研究提供了豐富的標(biāo)記資源，如應(yīng)用于鯉C.carpio L.［15］、斑點(diǎn)叉尾鮰［9］、虹鱒Oncorhynchus mykiss［16］、紅鰭東方鲀Takifugu rubripes［17-18］等魚類中。

表5 鯪和鯉轉(zhuǎn)錄組的比較Tab.5 Comparison of transcriptome between mud carp Cirrhina molitorella and common carp Cyprinus carpio L.

本研究中，采用454 測序技術(shù)對鯪進(jìn)行轉(zhuǎn)錄組測序，分析鯪轉(zhuǎn)錄組的特征，共獲得1 297 479 條reads，組裝后得到19 962 條contigs，其中7314 條contigs 獲得了GO 的功能注釋。另外，鑒定出5931個具有完整的ORF 的全長cDNA 序列，并且共鑒定出2438 個微衛(wèi)星和5014 個SNP 位點(diǎn)。以上結(jié)果均整合到了鯪轉(zhuǎn)錄組數(shù)據(jù)庫中。近年來，對鯪的遺傳和分子生物學(xué)方面的研究，初步闡述了鯪的遺傳多樣性［2］，開展了少量的基因克隆與分子標(biāo)記試驗(yàn)［1］，并探討了不同地理分布的鯪的進(jìn)化關(guān)系以及與肌肉發(fā)育相關(guān)的基因。然而，在分子育種層面，由于基礎(chǔ)數(shù)據(jù)資源的缺乏，尚無相關(guān)的研究進(jìn)展。本研究中獲得了鯪的轉(zhuǎn)錄組序列、基因注釋結(jié)果和數(shù)千個分子標(biāo)記，使該領(lǐng)域的研究人員可以大大加快研究的步伐，利用已有的家系樣本，可對目前關(guān)心鯪的生長、抗寒等性狀進(jìn)行深入地分析。本研究中獲得了鯪的轉(zhuǎn)錄組數(shù)據(jù)，并在此基礎(chǔ)上對序列進(jìn)行了注釋，獲得了大量微衛(wèi)星及SNP 信息，構(gòu)建了鯪的數(shù)據(jù)庫和網(wǎng)站，為深入開展鯪的生理生化機(jī)制研究、分子標(biāo)記輔助的遺傳育種、種群遺傳學(xué)和資源評估等研究提供了豐富的標(biāo)記資源。

［1］鐘茂春，鄭光明，趙建，等.鯪Myf5 基因克隆及其SNPs 分析［J］.中國水產(chǎn)科學(xué)，2010，17(4):681 -688.

［2］張丹丹，鄭光明，朱新平，等.西江野生鯪與養(yǎng)殖群體的遺傳分析［J］.華南農(nóng)業(yè)大學(xué)學(xué)報，2009，30(3):81 -85.

［3］劉佳瑤，趙建，鄭光明，等.鯪微衛(wèi)星DNA 分子標(biāo)記的篩選與遺傳多樣性分析［J］.基因組學(xué)與應(yīng)用生物學(xué)，2012，31(4):374 -380.

［4］Ji P，Liu G，Xu J，et al.Characterization of common carp transcriptome:sequencing，de novo assembly，annotation and comparative genomics［J］.PLoS One，2012，7(4):e35152.

［5］Wang Z，Gerstein M，Snyder M.RNA-Seq:a revolutionary tool for transcriptomics［J］.Nat Rev Genet，2009，10:57 -63.

［6］Anisimov S V.Serial analysis of gene expression(SAGE):13 years of application in research［J］.Curr Pharm Biotechnol，2008，9:338 -350.

［7］Liu S，Wang X，Sun F，et al.RNA -Seq reveals expression signatures of genes involved in oxygen transport，protein synthesis，folding and degradation in response to heat stress in catfish［J］.Physiol Genomics，2013，45(12):462 -476.

［8］Liu S，Zhang Y，Zhou Z，et al.Efficient assembly and annotation of the transcriptome of catfish by RNA - Seq analysis of a doubled haploid homozygote［J］.BMC Genomics，2012，13:595.

［9］Liu S，Zhou Z，Lu J，et al.Generation of genome-scale gene-associated SNPs in catfish for the construction of a high - density SNP array［J］.BMC Genomics，2011，12:53.

［10］Vesterlund L，Jiao H，Unneberg P，et al.The zebrafish transcriptome during early development［J］.BMC Dev Biol，2011，11:30.

［11］Zheng W，Wang Z，Collins J E，et al.Comparative transcriptome analyses indicate molecular homology of zebrafish swim bladder and mammalian lung［J］.PLoS One，2011，6(8):e24019.

［12］Shin S C，Kim S J，Lee J K，et al.Transcriptomics and comparative analysis of three Antarctic notothenioid fishes［J］.PLoS One，2012，7(8):e43762.

［13］Chen J，Li C，Huang R，et al.Transcriptome analysis of head kidney in grass carp and discovery of immune - related genes［J］.BMC Vet Res，2012，8:108.

［14］Sun F，Liu S，Gao X，et al.Male -biased genes in catfish as revealed by RNA - Seq analysis of the testis transcriptome［J］.PLoS One，2013，8(7):e68452.

［15］Xu J，Ji P，Zhao Z，et al.Genome-wide SNP discovery from transcriptome of four common carp strains［J］.PLoS One，2012，7(10):e48140.

［16］Salem M，Vallejo R L，Leeds T D，et al.RNA - Seq identifies SNP markers for growth traits in rainbow trout［J］.PLoS One，2012，7(5):e36264.

［17］Cui J，Liu S，Zhang B，et al.Transciptome analysis of the gill and swim bladder of Takifugu rubripes by RNA-seq［J］.PLoS One，2013，9(1):e85505.

［18］Cui J，Wang H，Liu S，et al.SNP discovery from transcriptome of the swim bladder of Takifugu rubripes［J］.PLoS One，2013，9(3):e92502.