大竹蟶高通量轉錄組測序數(shù)據(jù)組裝和分析

聶鴻濤， 姜力文， 鄭夢鴿， 李東東， 閆喜武

大竹蟶高通量轉錄組測序數(shù)據(jù)組裝和分析

聶鴻濤1，2， 姜力文1， 鄭夢鴿1， 李東東1， 閆喜武1，2

(1.大連海洋大學水產(chǎn)與生命學院，遼寧大連116023；2.遼寧省貝類良種繁育工程技術研究中心，遼寧大連116023)

為進一步開發(fā)大竹蟶Solen grandis的基因資源，采用2代Illumina Hi－seq測序技術對大竹蟶的鰓組織進行了轉錄組測序，構建了轉錄組數(shù)據(jù)庫，獲得338 483 476條Clean Reads數(shù)據(jù)；拼接組裝后獲得190 856條Unigene數(shù)據(jù)，平均長度為1147 bp；與NR、NT、KO、SwissProt、PFAM、GO、KOG等數(shù)據(jù)庫進行Blast信息比對 (E－value為10－5)，共獲得63 337個注釋基因；與NR數(shù)據(jù)庫比對發(fā)現(xiàn)，大竹蟶轉錄組基因序列與長牡蠣Crassostrea gigas具有較高的同源性，為53.3%；將大竹蟶轉錄組的Unigene的功能通過與KOG數(shù)據(jù)庫進行注釋比對劃分為25類；GO數(shù)據(jù)庫注釋可分為三類，即細胞組分、生物過程和分子功能，共包括65個分支；KEGG分析發(fā)現(xiàn)，大竹蟶轉錄組數(shù)據(jù)中按照代謝通路可分為92類，利用Blast蛋白庫比對和Estscan軟件進行ORF預測，獲得長度大于300 nt的ORF共50 681個；通過SSR分析，共獲得73 089個SSR標記。本研究中獲得的轉錄組信息可為今后進行大竹蟶分子標記的開發(fā)和關鍵基因的克隆及功能分析等研究提供基礎數(shù)據(jù)。

大竹蟶；轉錄組；高通量測序；生物信息學分析

大竹蟶Solen grandis在分類上隸屬于軟體動物門、雙殼綱、真瓣鰓目、竹蟶科、竹蟶屬，分布較為廣泛，是中國重要的經(jīng)濟貝類[1]。大竹蟶市場前景較好，經(jīng)濟價值高，因其肉味鮮美，深受消費者喜愛。大竹蟶的出肉率為45%，蛋白質含量高達82.7%[2]，必需氨基酸種類齊全，具有較高的營養(yǎng)價值[3]。目前，有關大竹蟶的研究主要集中在人工育苗技術[4－8]、 繁殖生物學[9－11]、 養(yǎng)殖生態(tài)學[12－14]和分子生物學[15－17]等方面， 而關于大竹蟶組學方面的研究較少。近年來，高通量轉錄組測序技術的廣泛應用[18]，為開展大竹蟶重要經(jīng)濟性狀的分子基礎研究提供了重要的基因資源。

轉錄組測序能夠在沒有基因資源的情況下，通過對某物種的RNA測序分析其生長與代謝規(guī)律，揭示基因與生物學特性之間的關系，同時能夠有效地開發(fā)大量基因資源[18]。該法特別適合基因信息匱乏的物種，并已在長牡蠣、扇貝、菲律賓蛤仔等多個海洋貝類中得到廣泛應用。但關于大竹蟶轉錄組測序的研究還鮮有報道。為了進一步開發(fā)大竹蟶的基因資源，本研究中以大竹蟶野生群體為材料，利用Illumina HiSeq 2500技術對大竹蟶進行轉錄組測序和分析，獲得了大竹蟶轉錄組數(shù)據(jù)庫，得到了與大竹蟶生長、繁殖、代謝、免疫等相關基因資源，并發(fā)現(xiàn)了大量的微衛(wèi)星分子標記 (SSR)，旨在為大竹蟶分子標記開發(fā)、基因克隆和組學研究提供基礎數(shù)據(jù)，同時也為今后大竹蟶功能基因的研究奠定基礎。

1 材料與方法

1.1 大竹蟶的暫養(yǎng)

試驗用大竹蟶購于大連市新長興市場，大竹蟶個體健康、活動性強，殼長為 (9.4±0.39)cm，殼寬為 (1.58±0.45)cm，殼高為 (1.61±0.56)cm，濕質量為 (43.68±0.73)g。試驗用海水為黑石礁海區(qū)經(jīng)沉淀、沙濾后的海水，儲存?zhèn)溆?。試驗開始前，先將大竹蟶暫養(yǎng)3 d，暫養(yǎng)期間不投喂，水溫為 (12.0±0.3)℃，每天換水1次。用水質分析儀測定海水的鹽度和pH，海水pH值為8.0±0.1， 鹽度為 32.0±0.2。

1.2 方法

1.2.1 大竹蟶總RNA的提取 采用TRIzol方法提取大竹蟶鰓組織的總RNA，測定RNA濃度和質量。如果OD260nm/OD280nm的范圍為1.8～2.2，28S/18S＞2， OD260nm/OD230nm＞1.8， RNA 完 整 值(RIN)＞8.5，則認為提取的RNA合格，符合建庫測序標注。

1.2.2 大竹蟶建庫測序和拼接組裝 大竹蟶轉錄組的建庫測序和拼接組裝委托北京諾禾致源生物信息科技有限公司完成，利用Trinity軟件[19]進行拼接組裝。

1.2.3 基因功能注釋 進行基因功能注釋的數(shù)據(jù)庫包括 NR、 NT、 PFAM、 KOG/COG、 SwissProt、GO、KEGG[20]。拼接得到的Unigene通過Blastx與基因和蛋白數(shù)據(jù)庫進行比對，得到注釋信息，利用WEGO和BlastGO進行GO注釋和功能分類統(tǒng)計，分別用COG和KEGG數(shù)據(jù)庫進行Unigene的功能注釋和代謝通路的定位。

1.2.4 ORF預測與SSR位點搜索 利用Blast比對到NR和SwissProt蛋白庫提取ORF序列，沒比對到的用Estscan軟件預測已注釋Unigene的ORF(Open Reading Frame)。微衛(wèi)星標記的搜索使用MISA軟件，對Unigenes進行SSR重復序列搜索(http://pgrc.ipk－gatersleben.de/misa/misa.htmh/)。

2 結果與分析

2.1 大竹蟶測序數(shù)據(jù)質量評價及序列組裝分析

經(jīng)HiSeq2500高通量測序獲得338 483 476條讀序 (clean data)，數(shù)據(jù)量為50.77 G，測序數(shù)據(jù)已上傳至SRA數(shù)據(jù)庫 (PRJNA385244)，GC含量為39.60%，Q30值為91.89%。利用Trinity軟件對上述獲得的clean reads進行組裝，得到305 485條transcripts，平均長823 bp，N50為1536 bp；進一步去冗組裝獲得190 856條Unigenes，平均長度為1147 bp，N50為1875 bp。拼接得到的Unigene的N50為1875 bp。Unigene長度分布見圖1。

2.2 Unigene的功能注釋

將組裝得到的190 856條Unigene與已知的NR、NT、 KO、 SwissProt、 PFAM、 GO、 KOG 這 7個數(shù)據(jù)庫進行功能注釋。如表1所示，全部Unigene在7個不同數(shù)據(jù)庫中的同源比對數(shù)目不同，分別為 46 657、 7857、 19 044、 34 627、 50 377、50 481和25 168條。至少在一個數(shù)據(jù)庫中獲得注釋信息的Unigene數(shù)量為63 337條，占總Unigene的33.18%；在7個數(shù)據(jù)庫中都得到注釋的基因數(shù)量為4744條 (表1)。

圖1 大竹蟶Unigene長度分布Fig.1 Sequence length distribution of unigenes assembled from Illumina reads of the gill transcriptome in grand jacknife clam Solen grandis

表1 大竹蟶基因注釋成功率統(tǒng)計Tab.1 Number of genes annotated in different database in grand jacknife clam Solen grandis

從圖2－A可以看出，比對到NR數(shù)據(jù)庫中的46 657條 Unigenes，有 16.4%的 E值分布在 le－15～1e－5， 14.4%分布在 le－30～1e－15， 9.7%分布在 le－45～1e－30， 7.9%分布在 le－60～1e－45，15.4%分布在 le－100～1e－60， 17.8%分布在0～1e－100。從圖2－B可以看出，從匹配的物種來源分析，大竹蟶Unigene注釋到長牡蠣Crassostrea gigas的占 53.3%，注釋到帽貝Lottia gigantea的占11.7%，注釋到加州海兔Aplysia californica的占8.0%，在紫海膽Strongylocentrotus purpura和文昌魚Branchiostoma floridae中注釋到的Unigene均占1.9%，其余23.2%注釋到其他物種。

圖2 NR注釋的E－value分布和物種分布Fig.2 E－value distribution and species classification of sequences matched to the NR database

2.3 Unigene的功能分類

在基因的GO注釋中，有50 481條Unigene(26.44%)獲得GO數(shù)據(jù)庫注釋 (圖3)。注釋到的3大類包括:細胞組分分為20個亞類，生物過程分為25個亞類，分子功能分為10個亞類。從圖3可見:生物過程獲得的注釋信息最多，為146 264條Unigene，占全部注釋信息的76.64%，生物過程里涉及代謝過程、細胞過程和對刺激的應激效應的Unigene比較多，分別有23 801、29 427、9526條；其次是細胞組分，分為20個亞類，共有87 560條Unigene，占全部注釋信息的45.88%，涉及到細胞部分、膜、細胞和細胞器的Unigene比較多，分別為15 763、11 013、15 763和10 280條；注釋信息相對最少的是分子功能，共有58 699條Unigene(30.76%)，排在前 2位的是結合與催化活性，Unigene數(shù)量分別為 27 847、19 094條。另外，KOG數(shù)據(jù)庫功能注釋結果顯示，KOG分類中得到25種，包括25 168條Unigene，其中，信號轉導機制包含Unigene最多 (5884條)，其次是一般功能預測 (5569條)，轉錄后修飾、蛋白折疊和分子伴侶 (2502條)，胞內運輸、分泌和囊泡運輸 (1252條)，轉錄 (1232條)，細胞骨架 (1232條)。

2.4 Unigene代謝途徑分析

在KEGG數(shù)據(jù)庫中分析了190 856條Unigene，得到注釋的Unigene有30 914條 (16.2%)，KEGG代謝途徑有230種。其中，細胞黏附的Unigene數(shù)量最多 (568條)，其次是 Rap1信號通路 (512條)、涉及內吞作用 (482條)、PI3K－Akt信號通路 (472條)、催產(chǎn)素信號通路 (417條)、Ras信號通路 (403條)、cAMP信號通路 (402條)、嘌呤代謝 (293條)等代謝通路。

圖3 大竹蟶轉錄組Unigene GO功能分類Fig.3 GO(Gene Ontology)categorization(biological process， cellular component and molecular function)of the unigenes in the gill transcriptome of grand jacknife clam Solen grandis

2.5 ORF預測與SSR位點分析

利用Blast比對到NR和SwissProt蛋白庫提取ORF序列，沒比對到的用Estscan預測Unigene，得到ORF共86 152個，其中，長度在300 nt以上的Unigene數(shù)量為50 681條，占58.83%(圖4)。同時，利用MISA軟件搜索到微衛(wèi)星標記73 089個，其中單堿基微衛(wèi)星最多，單堿基型 (50 676個，占69.33%)，其他類型依次為二堿基型重復(15 449個，占21.14%)、三堿基型 (5692個，占7.79%)、四堿基型 (1250個，占1.71%)，五堿基型和六堿基型最少，分別為20、2個，各占0.027%和0.003%(表2)。對大竹蟶微衛(wèi)星位點的分析，可為大竹蟶分子標記的開發(fā)與應用研究提供重要資源。

圖4 大竹蟶ORF長度分布Fig.4 ORF length distribution in the gill transcriptome of grand jacknife clam Solen grandis

表2 大竹蟶SSR分析統(tǒng)計結果Tab.2 Number of SSR marker in the gill transcriptome of grand jacknife clam Solen grandis

3 討論

3.1 大竹蟶鰓組織轉錄組組裝與注釋

近年來，高通量測序技術已在許多水產(chǎn)動物的轉錄組研究中得到廣泛應用，包括魚類[21－23]、貝類[24－26]和蝦蟹類等[27－29]。 本研究中， 采用高通量測序技術對大竹蟶的鰓組織進行轉錄組測序，拼接組裝后獲得190 856條Unigene，平均長度為1147 bp，其中有65 683條Unigenes的長度在1 kb以上。所有的Unigene有63 337條 (33.18%)成功獲得功能注釋。但仍有127 519條 (66.82%)序列定位不清楚，原因可能是基因數(shù)據(jù)庫中大竹蟶基因資源少，導致得到功能注釋基因偏少。對拼接組裝獲得的190 856條Unigenes進行了代謝途徑分析和功能分類。KOG功能注釋包括25 168條Unigenes，分為25個KOG亞類。本文中獲得了大竹蟶轉錄組數(shù)據(jù)庫，得到了與大竹蟶生長繁殖、合成代謝與免疫抗逆等相關基因資源，完善了大竹蟶的基因數(shù)據(jù)庫信息，為今后大竹蟶功能基因的研究奠定了基礎。

3.2 大竹蟶鰓組織轉錄組基因資源

近年來，有越來越多的水產(chǎn)經(jīng)濟動物轉錄組和基因組測序結果的發(fā)表，基因數(shù)據(jù)庫中的資源得到不斷地豐富和完善。在基因資源非常充足的情況下，很多未知功能的基因將被注釋到具體的功能，因此，可以利用這個途徑來挖掘一些新的功能基因。鰓是貝類重要的呼吸和濾食器官，研究表明，鰓在貝類的免疫方面也發(fā)揮重要作用[30－31]。轉錄組研究可以進行不同組織測序數(shù)據(jù)的比較分析，但本研究中選取了大竹蟶的鰓組織進行建庫測序，通過基因功能注釋發(fā)現(xiàn)了一些重要的免疫基因，如凋亡抑制因子 (Inhibitors of apoptosis，IAP)、腫瘤壞死因子 (Tumor necrosis factors，TNF)和補體 C3(complement component C3)等。本研究結果不僅為今后的基因克隆和功能研究提供了重要的基礎數(shù)據(jù)，同時為研究大竹蟶生長、免疫、繁殖、代謝等途徑中功能基因的分子機制奠定基礎。

3.3 大竹蟶轉錄組微衛(wèi)星分子標記

與傳統(tǒng)測序方法相比，高通量測序技術具有操作簡單、通量大、效率高等優(yōu)點。能夠快速地獲得成千上萬條基因序列，并開發(fā)出大量的微衛(wèi)星分子標記資源[32]。

本研究中，利用MISA微衛(wèi)星搜索軟件對大竹蟶轉錄組基因序列進行搜索，共發(fā)現(xiàn)微衛(wèi)星分子標記73 089個。在得到的6種微衛(wèi)星重復類型中，單堿基型重復數(shù)量最多，為50 676個，占所有微衛(wèi)星標記的69.33%，這與對櫛孔扇貝的研究結果類似[33]。目前，關于大竹蟶微衛(wèi)星標記開發(fā)的報道較少，在大竹蟶中僅有少量的微衛(wèi)星標記可以利用。因此，有必要進一步豐富大竹蟶的微衛(wèi)星標記資源。

本研究中進行了大竹蟶轉錄組的初步探究，大竹蟶的基因數(shù)據(jù)庫資源得到了補充和完善，為今后開展大竹蟶分子遺傳學和組學方面的相關研究奠定了重要基礎。本研究中獲得了大量的大竹蟶微衛(wèi)星候選標記，可為今后大竹蟶微衛(wèi)星分子標記開發(fā)、遺傳連鎖圖譜構建、家系分析和系譜鑒定、物種鑒定以及群體遺傳學等研究提供重要資源。

[1] 齊鐘彥.中國經(jīng)濟軟體動物[M].北京:中國農(nóng)業(yè)出版社，1998:233－234.

[2] 戴聰杰.大竹蟶軟體部分營養(yǎng)成分分析及其評價[J].集美大學學報:自然科學版，2002，7(4):304－308.

[3] 戴聰杰.大竹蟶軟體部分的氨基酸組成分析[J].莆田學院學報，2002，9(3):32－35.

[4] 陳愛華，姚國興，張志偉.大竹蟶生產(chǎn)性人工繁育試驗[J].海洋漁業(yè)，2009，31(1):66－72.

[5] 宋賢亭，于瑞海，馬培振，等.大竹蟶室內人工育苗技術研究[J].海洋湖沼通報，2015(4):56－60.

[6] 王雪梅，路宜華，豐愛秀，等.大竹蟶健康苗種培育新模式的研究[J].水產(chǎn)養(yǎng)殖，2012，33(8):14－16.

[7] 侯和要，牟乃海，宋全山，等.大竹蟶人工繁殖技術研究[J].齊魯漁業(yè)，2004，21(6):32－35.

[8] 楊輝.大竹蟶人工育苗技術研究[J].河北漁業(yè)，2015(5):41－69.

[9] 肖國強，柴雪良，邵艷卿，等.大竹蟶的繁殖生物學[J].海洋科學，2009，33(10):21－25.

[10] 吳楊平，陳愛華，姚國興，等.大竹蟶胚胎發(fā)生及稚貝發(fā)育基本特征[J].動物學雜志，2012，47(4):74－81.

[11] 李碧全，陳東輝，鄭杰民，等.大竹蟶性腺觀測與人工催產(chǎn)技術的研究[J].集美大學學報:自然科學版，2011，16(2):97－102.

[12] 侯和要，王君霞，彭作波，等.不同鹽度對大竹蟶存活的影響[J].齊魯漁業(yè)，2004，21(5):5－6，4.

[13] 陳愛華，姚國興，張志偉，等.溫度、鹽度和底質對大竹蟶稚貝生長及存活的影響[J].熱帶海洋學報，2010，29(5):94－97.

[14] 陳愛華，吳楊平，姚國興，等.底質和溫度對大竹蟶苗種生長存活的復合影響[J].江蘇農(nóng)業(yè)科學，2012，40(5):208－210.

[15] 張?zhí)?，劉相全，孫國華，等.大竹蟶AFLP分子標記反應體系的建立與優(yōu)化[J].安徽農(nóng)業(yè)科學，2011，39(29):17923－17926.

[16] 楊頂瓏，韋秀梅，楊建敏，等.大竹蟶(Solen grandis)鐵蛋白基因的克隆及其在轉錄水平上對微生物多糖的應答[J].海洋與湖沼，2013，44(3):664－669.

[17] Yuan Y，Li Q，Kong L F，et al.The complete mitochondrial genome of the grand jackknife clam，Solen grandis(Bivalvia:Solenidae):a novel gene order and unusual non－coding region[J].Molecular Biology Report，2012，39(2):1287－1292.

[18] Sánchez C C，Weber G M，Gao G T，et al.Generation of a reference transcriptome for evaluating rainbow trout responses to various stressors[J].BMC Genomics，2011，12:626.

[19] Grabherr M G，Haas B J，Yassour M，et al.Full－length transcriptome assembly from RNA－Seq data without a reference genome[J].Nature Biotechnology，2011，29(7):644－652.

[20] Kanehisa M，Goto S.KEGG:Kyoto encyclopedia of genes and genomes[J].Nucleic Acids Research，2000，28(1):27－30.

[21] Ju Z，Dunham R A，Liu Z.Differential gene expression in the brain of channel catfish(Ictalurus punctatus)in response to cold acclimation[J].Molecular Genetics and Genomics， 2002， 268(1):87－95.

[22] Long Y，Li L C，Li Q，et al.Transcriptomic characterization of temperature stress responses in larval zebrafish[J].PLoS One，2012，7(5):e37209.

[23] Mininni A N，Milan M，F(xiàn)erraresso S，et al.Liver transcriptome a-nalysis in gilthead sea bream upon exposure to low temperature[J].BMC Genomics，2014，15:765.

[24] Zhang L L，Li L，Zhu Y B，et al.Transcriptome analysis reveals a rich gene set related to innate immunity in the eastern oyster(Crassostrea virginica)[J].Marine Biotechnology，2014，16(1):17－33.

[25] Milan M，Coppe A，Reinhardt R，et al.Transcriptome sequencing and microarray development for the Manila clam，Ruditapes philippinarum:genomic tools for environmental monitoring[J].BMC Genomics，2011，12:234.

[26] Shi M J，Lin Y，Xu G R，et al.Characterization of the Zhikong scallop(Chlamys farreri)mantle transcriptome and identification of biomineralization－ related genes[J].Marine Biotechnology，2013，15(6):706－715.

[27] Zhao Q，Pan L Q，Ren Q，et al.Identification of genes differentially expressed in swimming crab Portunus trituberculatus response to low temperature[J].Aquaculture，2015，442:21－28.

[28] Shekhar M S，Kiruthika J，Ponniah A G.Identification and expression analysis of differentially expressed genes from shrimp(Penaeus monodon)in response to low salinity stress[J].Fish＆Shellfish Immunology，2013，35(6):1957－1968.

[29] Chen K，Li E C，Xu T Y，et al.Transcriptome and molecular pathway analysis of the hepatopancreas in the Pacific white shrimp Litopenaeus vannamei under chronic low－salinity stress[J].PLoS One，2015，10(7):e0131503.

[30] Wei X M，Yang J M，Liu X Q，et al.Identification and transcriptional analysis of two types of lectins(SgCTL－1 and SgGal－1)from mollusk Solen grandis[J].Fish ＆ Shellfish Immunology，2012，33(2):204－212.

[31] Nie H T，Liu L H，Huo Z M，et al.The HSP70 gene expression responses to thermal and salinity stress in wild and cultivated Manila clam Ruditapes philippinarum[J].Aquaculture，2017，470:149－156.

[32] Li H J，Liu M，Ye S，et al.De novo assembly，gene annotation，and molecular marker development using Illumina paired－end transcriptome sequencing in the clam Saxidomus purpuratus[J].Genes ＆ Genomics，2017，39(6):675－685.

[33] 趙柏淞.櫛孔扇貝(Chlamys farreri)BAC文庫的構建及其基因組特征分析[D].青島:中國海洋大學，2013.

Transcriptome data assembly and analysis of grand jackknife clam Solen grandis through high－throughput sequencing

NIE Hong－tao1，2， JIANG Li－wen1， ZHENG Meng－ge1， LI Dong－dong1， YAN Xi－wu1，2
(1.College of Fisheries and Life Science， Dalian Ocean University， Dalian 116023， China； 2.Engineering Research Center of Shellfish Culture and Breeding in Liaoning Province， Dalian 116023， China)

The Illumina Hi－seq sequencing technology was used to sequence the transcriptome in gills of grand jackknife clam Solen grandis and the data were analyzed by bioinformatics method to study the transcriptome of the grand jackknife clam.The transcriptome library of grand jackknife clam was found to contain 338 483 476 reads，and 190 856 Unigenes were obtained by assembling the Scaffolds in the transcriptome library，with an average length of 1147 bp.A total of 63 337 annotated genes were obtained compared with NR， KO， SwissProt， PFAM，GO， KOG and other databases for Blast(E－value 10－5).Comparing with the NR database， we found that the gene sequence in grand jackknife clam had a high homology(53.3%)with the Pacific oyster(Crassostrea gigas).Unigenes in the transcriptome of grand jackknife clam was divided into 25 classes according to the function by comparing Unigene with the KOG database.GO annotations was divided into 3 major categories:biological processes，cellular components and molecular functions including a total of 65 branches.The KEGG database revealed that the transcriptome data of grand jackknife clam was divided into 92 types according to the metabolic pathway.The ORF were predicted by Blast with protein database and estscan， 50 681 ORF with the length of more than 300 nt， and total of 73 089 SSR markers was found through repeat sequence motifs search and analysis.The transcriptome information in this study can provide the gene resources for cloning and functional analysis and molecular markers development of grand jackknife clam in the future.

Solen grandis； transcriptome； high－throughput sequencing； bioinformatics analysis

Q142；S968.3

A

10.16535/j.cnki.dlhyxb.2017.06.004

2095－1388(2017)06－0658－06

2017－03－30

現(xiàn)代農(nóng)業(yè)產(chǎn)業(yè)技術體系建設專項 (CARS－48)；遼寧省農(nóng)業(yè)領域青年科技創(chuàng)新人才項目 (2014004)；大連市科技計劃項目(2014B11NC092，2016RQ065)

聶鴻濤 (1984—)，男，博士，副研究員。E－mail:htnie＠dlou.edu.cn

閆喜武 (1962—)，男，博士，教授。E－mail:yanxiwu＠dlou.edu.cn