大竹蟶水管自切后再生過程的轉錄組學分析

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

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

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

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

1 材料與方法

1.1 材料

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

1.2 方法

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

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

1.2.2 RNA提取、cDNA文庫構建和高通量測序

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

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

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

采用STEM軟件[22]對3個比較組中的并集差異基因進行時間趨勢分析，模塊數量設定為20，獲得顯著的基因表達模式，并選取模式中的基因集進行GO富集分析。

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

表1 引物信息Tab.1 Primer information

2 結果與分析

2.1 大竹蟶水管再生的形態特征及再生長情況

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

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

2.2 轉錄組測序結果和組裝

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

表2 轉錄組測序結果Tab.2 Results of transcriptome sequencing

2.3 Unigenes功能注釋

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

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

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

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

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

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

KEGG通路富集結果顯示，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信號通路(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 差異基因的時間趨勢分析

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

方框左上角數字代表不同基因表達趨勢模式，如Profile5；左下方數字代表差異基因個數；彩色方框表示具有顯著性。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 差異基因的時間趨勢分析Fig.5 Time trend analysis of differentially expressed genes

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

2.6 RT-qPCR驗證

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

圖7 實時熒光定量PCR驗證結果Fig.7 Verification results by real-time quantitative PCR

3 討論

3.1 Unigenes功能注釋信息分析

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

3.2 大竹蟶水管組織受損的早期響應

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

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

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

3.4 大竹蟶水管再生期間的代謝響應

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

4 結論

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

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