大竹蟶水管自切后再生過程的轉(zhuǎn)錄組學(xué)分析

孫雪峰，陳愛華，張志東，陳素華，張雨，曹奕，朱艷青，楊家新，蔣建斌，吳楊平*

(1.南京師范大學(xué) 海洋科學(xué)與工程學(xué)院，江蘇 南京 210023；2.江蘇省海洋水產(chǎn)研究所 江蘇省文蛤良種場，江蘇 南通 226007；3.南通市通州區(qū)水產(chǎn)技術(shù)指導(dǎo)站，江蘇 南通 226300)

大竹蟶(Solengrandis)，又稱蟶子王，屬軟體動(dòng)物門(Mollusca)瓣鰓綱(Lamellibranchia)簾蛤目(Veneroida)竹蟶科(Solenida)，在菲律賓、朝鮮、日本及中國南北沿海潮間帶均有分布，是常見的食用貝類之一[1]。大竹蟶體大殼薄，出肉率高，蟶肉中富含蛋白質(zhì)、鈣、鐵、硒和維生素 A 等營養(yǎng)元素，為貝類中的上品，其市場價(jià)格高達(dá)140元/kg。大竹蟶在繁育、養(yǎng)殖及運(yùn)輸過程中頻繁出現(xiàn)水管自切再生現(xiàn)象，水管發(fā)生自切后大竹蟶會(huì)出現(xiàn)伏沙不潛底，有的出現(xiàn)疲軟甚至死亡，導(dǎo)致親貝催產(chǎn)利用率降低，給育苗企業(yè)帶來種質(zhì)浪費(fèi)和經(jīng)濟(jì)損失。在養(yǎng)殖生產(chǎn)中，發(fā)生自切后的大竹蟶生長緩慢，出池時(shí)間延長，養(yǎng)殖風(fēng)險(xiǎn)加大，如果大竹蟶頻繁發(fā)生自切再生會(huì)導(dǎo)致所謂的“長老”現(xiàn)象。在銷售過程中，商品蟶由于受裝卸運(yùn)輸?shù)膽?yīng)激影響，大竹蟶普遍發(fā)生水管自切，因掉落水管的大竹蟶無人購買，只能丟棄，而水管一般占到個(gè)體總質(zhì)量的20%～25%，對(duì)其收益影響較大。因此，大竹蟶水管自切再生問題已成為制約其養(yǎng)殖及銷售的關(guān)鍵因素之一。

再生是指生物體的器官甚至整體受到自然環(huán)境中外力作用后導(dǎo)致創(chuàng)傷或部分丟失，在剩余部分的基礎(chǔ)上長出與損傷或缺失部分相同的形態(tài)和功能組織的修復(fù)過程。組織損傷的修復(fù)和再生是生命科學(xué)研究的重要課題之一，國內(nèi)外對(duì)再生的研究主要集中在一些模式動(dòng)物上，如渦蟲(Schmidteamediterranea)、蠑螈(Ambystomamexicanum)、小鼠(Musmusculus)和斑馬魚(Daniorerio)等[2-5]。再生通常包括已經(jīng)存在組織的重組，成體干細(xì)胞的利用，細(xì)胞的去分化和轉(zhuǎn)分化等，進(jìn)而實(shí)現(xiàn)形態(tài)、結(jié)構(gòu)和功能的重建[6]。目前，國內(nèi)外關(guān)于水生動(dòng)物再生的研究主要有蠑螈的肢體再生[7]、斑馬魚的心臟和尾鰭再生[8-9]、海盤車腕部再生[10]、水螅軀體再生[11]和海參內(nèi)臟團(tuán)再生[12]等。研究中發(fā)現(xiàn)，再生過程對(duì)動(dòng)物的攝食能力、生長和發(fā)育等有較大影響，如刺參出現(xiàn)體質(zhì)量負(fù)生長現(xiàn)象[13]，蟹斷肢后會(huì)出現(xiàn)生長減緩、發(fā)育時(shí)間(蛻皮或蛻皮周期)延長[14]。再生過程同時(shí)還受到包括生長因子、細(xì)胞因子等生物活性分子，以及神經(jīng)、體液和信號(hào)通路等因素調(diào)控，涉及基因、蛋白和細(xì)胞等的作用[15]。在貝類中，僅見關(guān)于菲律賓蛤仔(Ruditapesphilippinarum)、縊蟶(Sinonovaculaconstricta)水管再生的研究[16-18]，而有關(guān)大竹蟶組織再生研究僅見筆者所在課題組發(fā)表的利用代謝組學(xué)方法評(píng)估水管自切對(duì)大竹蟶代謝特性影響的論文[19]。本研究中，以人為切除水管組織的大竹蟶為試驗(yàn)對(duì)象，通過二代測序技術(shù)對(duì)大竹蟶組織再生進(jìn)行轉(zhuǎn)錄組測序研究，并分析水管再生期間轉(zhuǎn)錄組的變化情況，以期為深入研究大竹蟶水管再生機(jī)制提供科學(xué)參考。

1 材料與方法

1.1 材料

試驗(yàn)在江蘇省海洋水產(chǎn)研究所呂四試驗(yàn)基地進(jìn)行。以項(xiàng)目組自繁養(yǎng)殖8～10 cm殼長的大竹蟶為試驗(yàn)材料，選取無破損、活力好的大竹蟶于室內(nèi)1 000 L圓缸中暫養(yǎng)7 d，圓缸內(nèi)底鋪20 cm厚的細(xì)沙，加注沙濾海水600 L。暫養(yǎng)期間，海水鹽度為28，pH為8.3，水體溶氧量為7.6～8.0 mg/L，水溫為26～30 ℃。暫養(yǎng)期間缸內(nèi)海水保持連續(xù)充氣，每日上午換水一次，換水量為總水量的50%；選擇密度為(15～20)×105ind./mL 的球等鞭金藻(Isochrysisgalbana)，早、中、晚定時(shí)往暫養(yǎng)缸內(nèi)加入120 L藻水，保證大竹蟶處于正常養(yǎng)殖狀態(tài)。

1.2 方法

1.2.1 試驗(yàn)設(shè)計(jì)及樣品采集 以預(yù)試驗(yàn)中大竹蟶水管自切后的生長情況為參考，分別設(shè)置自切前對(duì)照組，以及自切后7 h、7 d和20 d 3個(gè)選擇組。由于在pH、水壓等外部因素刺激下大竹蟶會(huì)發(fā)生水管自切的隨意性，水管斷裂位置的不同造成斷裂水管長短不一，因此，在正式試驗(yàn)開始時(shí)，采用人工切除水管的方法統(tǒng)一將水管切除至水管基部(圖1)。從暫養(yǎng)缸中逐一取出大竹蟶，切除水管后的大竹蟶共計(jì)60條，置于重新準(zhǔn)備的1 000 L圓缸中飼養(yǎng)，飼養(yǎng)條件與暫養(yǎng)期間一致；飼養(yǎng)期間觀察缸中大竹蟶存活情況，發(fā)現(xiàn)不鉆沙個(gè)體及時(shí)清除。分別取自切前，以及自切后7 h、7 d和20 d 4個(gè)時(shí)間點(diǎn)的大竹蟶水管基部組織樣品，每組設(shè)置3個(gè)重復(fù)，樣品對(duì)應(yīng)標(biāo)號(hào)記為BC、ATa、ATb和ATc，所取組織樣本用液氮速凍后置于-80 ℃超低溫冰箱中保存?zhèn)溆谩?/p>

圖1 大竹蟶水管組織的切除Fig.1 Removed siphon of Solen grandis

1.2.2 RNA提取、cDNA文庫構(gòu)建和高通量測序

采用Trizol試劑盒(Invitrogen公司)提取總RNA，采用10 g/L瓊脂糖凝膠電泳檢測RNA的完整性，采用NanoDrop微量分光光度計(jì)測定RNA純度。RNA檢測合格后使用帶有Oligo(dT)的磁珠富集mRNA，再使用片段緩沖劑將得到的mRNA裂解為片段，在ProtoScript Ⅱ逆轉(zhuǎn)錄酶體系下用隨機(jī)引物反轉(zhuǎn)錄合成cDNA的第一條鏈，隨后用RNase H降解mRNA，在DNA聚合酶Ⅰ體系下以dNTPs為原料合成cDNA的第二條鏈。采用QiaQuick PCR 提取試劑盒純化cDNA，再加上polyA末端并連接測序接頭。利用瓊脂糖凝膠電泳篩選200 bp左右的連接產(chǎn)物并進(jìn)行PCR擴(kuò)增，經(jīng)AMPure XP beads純化后的PCR產(chǎn)物用以構(gòu)建cDNA文庫，并通過IlluminaNova Seq 6000測序平臺(tái)進(jìn)行高通量測序。

1.2.3 測序數(shù)據(jù)組裝與Unigenes注釋 測序得到的原始序列(raw reads)通過Fastp 0.18.0軟件過濾低質(zhì)量序列后獲得有效序列(clean reads)，采用Trinity軟件[20]對(duì)clean reads進(jìn)行組裝，并挑出每個(gè)基因中最長的轉(zhuǎn)錄本作為Unigene，將得到的Unigenes與COG、GO、KEGG、KOG、Pfam、Swiss-Prot和NR等數(shù)據(jù)庫進(jìn)行比對(duì)，獲得基因功能注釋。

1.2.4 差異基因(differentially expressed genes，DEGs)的表達(dá)及時(shí)間趨勢分析 以FPKM(fragments per kilobase of transcript per million mapped reads)值作為基因表達(dá)量高低的指標(biāo)，采用DESeq軟件進(jìn)行組間差異分析，篩選FDR(false discovery rate)<0.05且|log2(fold change)|>1的基因作為差異基因[21]，再將得到的差異基因通過GO、KEGG數(shù)據(jù)庫進(jìn)行富集分析。

采用STEM軟件[22]對(duì)3個(gè)比較組中的并集差異基因進(jìn)行時(shí)間趨勢分析，模塊數(shù)量設(shè)定為20，獲得顯著的基因表達(dá)模式，并選取模式中的基因集進(jìn)行GO富集分析。

1.2.5 RT-qPCR驗(yàn)證 根據(jù)測序結(jié)果，以β-actin為內(nèi)參基因[23]，隨機(jī)挑選6個(gè)差異基因進(jìn)行熒光定量(表1)。具體步驟：基于選擇基因的CDS序列，采用Primer 5軟件設(shè)計(jì)引物；以反轉(zhuǎn)的cDNA為模板，采用熒光定量試劑盒(SYBR Green)，在ABI公司的7 300 plus熒光定量PCR儀上進(jìn)行操作。PCR反應(yīng)體系為(20 μL)：2×SuperReal PreMix Plus 10 μL，正、反向引物各0.6 μL，cDNA模板1 μL，50×ROX Reference Dye 1 μL，用RNase-free ddH2O補(bǔ)足至20 μL。擴(kuò)增程序：95 ℃下預(yù)變性15 min；95 ℃下變性10 s，56 ℃下退火30 s，72 ℃下延伸31 s，共進(jìn)行40個(gè)循環(huán)。使用2-ΔΔCt法計(jì)算基因相對(duì)表達(dá)量。

表1 引物信息Tab.1 Primer information

2 結(jié)果與分析

2.1 大竹蟶水管再生的形態(tài)特征及再生長情況

從圖2可見：試驗(yàn)開始后5 h時(shí)，除5條大竹蟶未鉆沙外，其余大竹蟶均鉆沙；7 h后觀察到沙面有氣孔，部分大竹蟶水管伸出沙面濾食，發(fā)現(xiàn)被切除水管的基部增厚，外緣處新生多個(gè)凸起；7 d后清晰可見大竹蟶水管增長，平均增長(6.52±1.02)mm；20 d后大竹蟶水管增長明顯，平均增長(38.67±2.11)mm。在此期間，除10、15 d時(shí)分別有2、4條大竹蟶死亡外，最終大竹蟶存活率為81.67%。

S—水管；SR—水管基部。S—siphon；SR—siphon root.圖2 大竹蟶水管恢復(fù)生長Fig.2 Siphon regrowth of Solen grandis

2.2 轉(zhuǎn)錄組測序結(jié)果和組裝

從表2可見，原始序列經(jīng)質(zhì)控后共得到241 840 528條有效序列，Q20、Q30分別高于97%和92%，說明測序結(jié)果良好。有效序列經(jīng)Trinity組裝后共獲得235 719條Transcripts和103 909條Unigenes，Transcripts與Unigenes的N50分別為1 523、1 377 bp。

表2 轉(zhuǎn)錄組測序結(jié)果Tab.2 Results of transcriptome sequencing

2.3 Unigenes功能注釋

分別在COG、GO、KEGG、KOG、Pfam、Swiss-Prot和NR數(shù)據(jù)庫中對(duì)Unigenes進(jìn)行注釋分析。從表3可見，在COG數(shù)據(jù)庫中注釋到7 586條Unigenes，在GO、KEGG、KOG、Pfam、Swiss-Prot和NR數(shù)據(jù)庫中分別注釋到15 586、21 669、16 332、18 540、16 910和29 326條Unigenes，累計(jì)被注釋的Unigenes為30 027條，注釋率為28.90%，表明依舊有大量的未知基因有待研究。

表3 Unigenes功能注釋分析Tab.3 Analysis of Unigenes functional annotation

2.4 差異基因的篩選及GO、KEGG分析

差異基因篩選結(jié)果顯示，BC vs ATa組共篩選出7 098個(gè)差異基因，其中，3 343個(gè)基因上調(diào)，3 755個(gè)基因下調(diào)；BC vs ATb組共篩選出5 075個(gè)差異基因，其中，3 172個(gè)基因上調(diào)，1 903個(gè)基因下調(diào)；BC vs ATc組共篩選出4 245個(gè)差異基因，其中，2 561個(gè)基因上調(diào)，1 684個(gè)基因下調(diào)；862個(gè)基因在3組中共表達(dá)(圖3)。

GO功能注釋結(jié)果顯示，差異基因被分到生物過程(biological process)、細(xì)胞組分(cellular component)和分子功能(molecular function)3個(gè)類別中。Padj值為多重假設(shè)檢驗(yàn)校正后的P值，Padj值越小，表明富集程度越高，為了進(jìn)一步了解差異基因所行使的分子功能，本研究中篩選了分子功能類別中顯著富集(Padj<0.05)的GO條目。從表4可見，BC vs ATa組顯著富集到4個(gè)GO條目，其中，鈣離子結(jié)合(calcium ion binding)條目的富集程度最高且富集的基因數(shù)目也最多，171個(gè)差異基因富集到該條目；BC vs ATb、BC vsATc兩組中富集程度最高的條目均為幾丁質(zhì)結(jié)合(chitin binding)。

表4 3組差異基因顯著富集的GO條目Tab.4 GO Terms with significant enrichment of three groups of DEGs

KEGG通路富集結(jié)果顯示，BC vs ATa、BC vs ATc兩組無通路被顯著富集(Padj>0.05)，BC vs ATb組共有5條通路被顯著富集(Padj<0.05)，依次為氨基糖與核苷酸糖代謝(amino sugar and nucleotide sugar metabolism)、氮素代謝(nitrogen metabolism)、黏著斑(focal adhesion)、溶酶體(lysosome)和MAPK信號(hào)通路(MAPK signaling pathway)。圖4展示了富集程度前20的KEGG通路。

圖3 差異基因的篩選Fig.3 Screening of differentially expressed genes

圖4 BC vs ATb組差異基因的KEGG富集通路(前20)Fig.4 KEGG enrichment pathway of differentially expressed genes (top 20)

2.5 差異基因的時(shí)間趨勢分析

為了解差異基因隨時(shí)間的變化情況，取3個(gè)比較組中差異基因的并集，采用STEM軟件進(jìn)行時(shí)間趨勢分析，模塊數(shù)量設(shè)定為20。從圖5可見，富集基因數(shù)前3的模式依次為Profile5、Profile17、Profile12，分別富集到1 966、1 093、981條基因。Profile5的基因表達(dá)量在7 h時(shí)呈下降趨勢，在7、20 d時(shí)表達(dá)量恢復(fù)并保持自切前水平，Profile17的基因表達(dá)量自切后高于自切前，Profile12的基因表達(dá)量在7 h時(shí)保持不變，隨后上升并保持不變。

方框左上角數(shù)字代表不同基因表達(dá)趨勢模式，如Profile5；左下方數(shù)字代表差異基因個(gè)數(shù)；彩色方框表示具有顯著性。The number on the upper left of the box denotes different trends，as Profile5.The number on the lower left represents the number of differential genes. Colored box representation is significant.圖5 差異基因的時(shí)間趨勢分析Fig.5 Time trend analysis of differentially expressed genes

選取Profile5、Profile17和Profile12 3個(gè)模式的基因集進(jìn)行GO功能富集分析。結(jié)果顯示，Profile5中共有529條基因注釋到離子結(jié)合(ion binding)條目，占總注釋基因數(shù)的53.77%，其次注釋到基因數(shù)較多的條目為蛋白質(zhì)結(jié)合(protein binding)、鈣離子結(jié)合(calcicium ion binding)等(圖6)；Profile17富集基因數(shù)較多的條目包括受體活性(receptor activity)和分子傳導(dǎo)活性(molecular transducer activity)等，富集程度最高的條目為G蛋白偶聯(lián)受體活性(G-protein coupled receptor activity)；Profile12中基因主要顯著富集到氨基糖代謝(amino sugar metabolic process)和幾丁質(zhì)代謝(chitin metabolic process)等條目。

2.6 RT-qPCR驗(yàn)證

為了驗(yàn)證測序結(jié)果的可靠性，隨機(jī)挑選了NR數(shù)據(jù)庫注釋下的6條顯著差異基因進(jìn)行熒光定量PCR(RT-qPCR)驗(yàn)證，結(jié)果表明(圖7)，相對(duì)定量趨勢與RNA-Seq的趨勢總體一致，表明測序結(jié)果可靠。

圖7 實(shí)時(shí)熒光定量PCR驗(yàn)證結(jié)果Fig.7 Verification results by real-time quantitative PCR

3 討論

3.1 Unigenes功能注釋信息分析

相關(guān)研究中，楊旻珉[24]獲得了羅氏沼蝦(Macrobrachiumrosenbergii)肝胰腺轉(zhuǎn)錄組142.38 Gb clean data，Q30在94.49%以上，組裝得到252 123條轉(zhuǎn)錄本和135 239條Unigenes，N50為1 737 bp，注釋率為29.75%。楊軼博[25]基于高通量測序?qū)θ毡救菧u蟲(Dugesiajaponica)再生過程進(jìn)行轉(zhuǎn)錄組分析，組裝獲得163 329條轉(zhuǎn)錄本和101 183條Unigenes，Q30在90.5%以上，轉(zhuǎn)錄本的N50為1 346 bp，注釋率為30.78%。聶鴻濤等[26]開展了大竹蟶高通量轉(zhuǎn)錄組測序數(shù)據(jù)組裝和分析，得到305 485條轉(zhuǎn)錄本，去冗組裝獲得190 856條Unigenes，N50為1 875 bp，Q30為91.89%，注釋率為33.18%。殷麗坤等[27]對(duì)低氧脅迫下中華圓田螺(Cipangopaludinacathayensis)肝臟組織進(jìn)行高通量測序，發(fā)現(xiàn)只有11.46% Unigenes得到注釋。本研究中，基于高通量測序技術(shù)對(duì)大竹蟶水管自切前后的水管基部組織進(jìn)行轉(zhuǎn)錄組測序，獲得235 719個(gè)轉(zhuǎn)錄本和103 909條Unigenes，Q30達(dá)92%以上，將其與多個(gè)數(shù)據(jù)庫比對(duì)共有30 027個(gè)Unigenes得到注釋，注釋率為28.90%，部分基因未得到注釋可能是數(shù)據(jù)庫中近源物種基因注釋信息較少，還有大量的未知基因有待研究。因此，加強(qiáng)大竹蟶全基因的分析研究，將有助于提高基因的注釋率。

3.2 大竹蟶水管組織受損的早期響應(yīng)

生物體內(nèi)具有不依賴于轉(zhuǎn)錄水平的快速響應(yīng)損傷的保守分子機(jī)制，研究發(fā)現(xiàn)，多種小分子物質(zhì)可以作為信號(hào)分子對(duì)組織損傷做出響應(yīng)(如Ca2+、ATP、ROS等)[6]，這些信號(hào)分子將損傷信號(hào)傳遞給細(xì)胞或組織，從而引起一系列響應(yīng)損傷的生物學(xué)過程[28]。Ca2+作為生物體內(nèi)含量最高的金屬離子，在組織受損的快速應(yīng)答中發(fā)揮重要作用，其通過受損細(xì)胞膜內(nèi)流進(jìn)入細(xì)胞與多種蛋白結(jié)合[29]，提高多種酶活性繼而引導(dǎo)組織修復(fù)和再生[30]，同時(shí)Ca2+借助ATP還可以將受損信號(hào)向受損組織更遠(yuǎn)的方向傳播[31]。此外，ROS也對(duì)鈣離子傳播受損信號(hào)具有重要影響[32]。本研究中，Profile5基因集的GO功能注釋發(fā)現(xiàn)，自切后7 h大量差異基因顯著富集于離子結(jié)合條目(圖6)，其中，鈣離子結(jié)合條目有98條基因富集到該分子功能條目且富集程度最高。這表明，在大竹蟶水管再生早期或大竹蟶對(duì)受損組織的早期響應(yīng)過程中，鈣離子很可能充當(dāng)信號(hào)分子，在傳遞受損信號(hào)中具有重要作用。

3.3 大竹蟶水管再生對(duì)受體活性的影響

G蛋白偶聯(lián)受體(GPCR)是真核生物中最大的一類膜蛋白家族，具有接受并傳導(dǎo)多種細(xì)胞外信號(hào)分子的功能，其在受損組織再生過程中發(fā)揮重要作用。周曉春[15]從轉(zhuǎn)錄組水平發(fā)現(xiàn)，蠑螈肢體再生過程中GPCR介導(dǎo)的信號(hào)傳導(dǎo)活動(dòng)可促進(jìn)細(xì)胞分化及組織再生。郭會(huì)萍[33]在斑馬魚性腺和心臟再生研究中發(fā)現(xiàn)，GPCR可以調(diào)控p-PKC/p-Akt活性，影響受損組織芽基生成及可塑性。本研究中，Profile17基因集GO功能注釋結(jié)果顯示，以GPCR活性為代表的受體活性條目被顯著富集，其自切后表達(dá)量均高于自切前，RT-qPCR驗(yàn)證了GPR84在自切后表達(dá)量顯著高于自切前。這表明，在大竹蟶水管再生過程中，GPR84很可能充當(dāng)接收和傳導(dǎo)各種再生相關(guān)信號(hào)的信號(hào)分子，并能激活下游信號(hào)通路。

3.4 大竹蟶水管再生期間的代謝響應(yīng)

組織再生是一個(gè)耗能過程，DNA合成、細(xì)胞有絲分裂和細(xì)胞遷移均會(huì)消耗大量的能量[34]，糖類代謝是生物體內(nèi)能量的主要來源，途徑包括糖酵解和三羧酸循環(huán)等。本研究中，BC vs ATb組GO功能注釋和KEGG通路富集分析顯示，多個(gè)與糖代謝相關(guān)的條目與通路被顯著富集，如幾丁質(zhì)代謝、氨基糖代謝、核苷酸糖代謝和氮素代謝等。Profile12中基因集GO富集結(jié)果同樣顯示，幾丁質(zhì)代謝、氨基糖代謝等代謝過程被顯著富集，并且其基因表達(dá)量在水管再生早期(自切后7 h)無明顯變化，但隨著再生的繼續(xù)，7 d和20 d時(shí)基因表達(dá)明顯上調(diào)。這與Wang等[17]在菲律賓蛤仔水管再生過程的研究中發(fā)現(xiàn)，再生中后期大量基因被富集到能量代謝和生物合成的結(jié)果類似。這表明，大竹蟶水管在響應(yīng)早期損傷后，物質(zhì)代謝及生物合成過程可能大量增加以滿足組織再生的需求。

4 結(jié)論

1)通過高通量測序技術(shù)分析大竹蟶水管再生不同時(shí)間點(diǎn)的轉(zhuǎn)錄組變化情況發(fā)現(xiàn)，再生不同時(shí)間點(diǎn)基因表達(dá)有所差異，差異基因主要富集在離子結(jié)合、受體活性和能量代謝等過程中。

2)初步推測大竹蟶水管再生早期，以鈣離子為代表的信號(hào)分子在傳遞受損信號(hào)過程中發(fā)揮重要作用，同時(shí)多種受體活性增強(qiáng)響應(yīng)水管再生過程，隨著再生的進(jìn)行，能量代謝和生物合成過程大量增加以滿足組織再生的需求。本研究中初步分析了大竹蟶水管再生過程中的轉(zhuǎn)錄組變化情況，為深入研究水管再生機(jī)制提供了有益參考。