玉足海參變態附著階段轉錄組分析

吳曉鵬，黃敏偉，陳曉瑛，彭 凱，趙吉臣，鐘 平，劉鳳坤，張業輝，黃 文

1.廣東海洋大學 水產學院，廣東 湛江 524088

2.廣東省農業科學院動物科學研究所/農業農村部華南動物營養與飼料重點實驗室/廣東省畜禽育種與營養研究重點實驗室，廣東 廣州 510640

3.廣東省農業科學院水產協同創新中心，廣東 廣州 510640

4.暨南大學 生命科學學院，廣東 廣州 510632

5.廣東省農業科學院蠶業與農產品加工研究所，廣東 廣州 510640

玉足海參 (Holothurialeucoisplota) 隸屬于楯手目、海參科、海參屬[1]，是一種具有較高經濟價值的熱帶海參。在我國主要分布于海南、廣東等南部沿海海域的巖礁、珊瑚礁、泥沙質淺海。其體形呈圓筒狀，后部比前端粗大，成參可以長到20 cm 以上，口偏于腹面，居維氏器很發達，腹面散生有大量呈乳狀突起的管足且排列不規則[2]。玉足海參具有環境適應能力強、耐應激的特點[3]，在維持近岸海洋生態系統穩定中扮演了關鍵角色，是進行放流移植及島礁生態環境修復工程的理想種類[3-4]。此外，海參自古以來便被視為滋補珍品，在東亞地區被廣泛食用以及用作制藥原料[5-6]。國內外學者對海參進行營養組成測定和營養評價，發現其富含膠原蛋白、多糖、皂苷、維生素、多肽、腦苷酸等多種活性物質。此外，對這些活性因子的研究也在逐步深入，目前已經證實其具有抗腫瘤、抗氧化、抗凝血、抗病毒和提高機體免疫力等多種生理功效[7-9]。

隨著人們對海參營養價值認知的提高，對其需求也日益增大，導致熱帶海參的過度捕撈，使其自然資源急劇減少[10]。目前，全球超過一半的海參資源被過度開發，約30%的海參漁業已經崩潰[11]。市場需求的劇增以及資源的枯竭使得熱帶海參人工繁育和養殖技術的研發顯得尤為迫切[12]。目前，已有幾種熱帶海參成功實現了人工繁殖[3,13-16]，但這些人工繁育技術尚處于探索階段，其繁育規模小、幼苗數量少，與應用于人工養殖及增殖放流仍有一定距離。玉足海參的繁育和養殖研究尚處于起步階段，目前可實現人工催產和孵化[3-4,17]。玉足海參幼蟲在耳狀幼體期間體積逐漸膨大，而在受精約15 d 后 (水溫為29～33 ℃)，體積開始明顯縮小形成樽形幼體，生活狀態也由浮游變為底棲[4]。而這一變態過程的高死亡率是熱帶海參繁育過程中的共性問題，其背后的機制至今不明。本研究基于全長轉錄組測序技術[18-20]，使用Illumina 測序平臺對4 個不同發育時期的海參幼體 (小耳幼體、中耳幼體、大耳幼體和樽形幼體) 進行了高通量測序分析。通過研究其幼體不同時期的基因表達差異，探究與海參生長發育相關的重要信號通路并發掘關聯基因，旨在提高對該物種生長發育機制的認識，為其育苗和養殖提供分子層面的育種研究借鑒，并為熱帶海參增養殖業發展及科學研究提供參考依據。

1 材料與方法

1.1 實驗材料

實驗海參捕撈于深圳大鵬灣海域，并于廣東省茂名市電白區金陽生物技術有限公司養殖基地暫養。玉足海參親本的性腺成熟時期一般在6—8 月，此時是其最佳的繁育時期。因此，本實驗在2022 年7 月進行海參親本的采集、人工繁育及不同發育時期幼體樣品的收集。首先從池中撈出親本進行陰干處理，然后放入桶中待產卵排精。收集精子和卵子，經濕法受精后轉入孵化池進行孵化。其后對其幼體發育過程進行監測，并分別收集小耳幼體、中耳幼體、大耳幼體和樽形幼體4 個生長發育時期的個體 (圖1) ，其中每個時期采集3 份樣品，每份樣品包含3 個同一時期的個體。將收集的樣品立即放入液氮罐中冷凍，運至實驗室于-80 ℃冰箱中保存備用。

圖1 玉足海參幼體 4 個發育時期的形態注: a.小耳幼體 (A 期)； b.中耳幼體 (B 期)；c.大耳幼體 (C 期)；d.樽形幼體 (D 期)。Fig.1 Morphology of H.leucospilota larva at four developmental stagesNote: a.Early auricularia (Stage A); b.Mid auricularia (Stage B); c.Late auricularia (Stage C); d.Doliolaria (Stage D).

1.2 RNA 提取與轉錄組測序

按常規Trizol 法提取玉足海參幼體樣品的總RNA，并標記如下：小耳幼體 (A 期)、中耳狀幼體(B 期)、大耳狀幼體 (C 期)、樽形幼體 (D 期)，每個時期的3 個重復樣品用1、2 和3 區別?？俁NA 的質量濃度大于100 ng·μL-1，OD260與OD280比值介于 1.8～2.2，OD260與OD230比值大于 2.0，確保實驗室中所提取的RNA 無污染、無降解。將提取的RNA 樣品送至上海美吉生物有限公司，經文庫構建后采用Illumina HiSeqTM 6000 高通量測序平臺進行轉錄組測序。

1.3 轉錄組從頭組裝

測序所獲得的初始圖像數據經Base Calling 處理轉成序列數據，使用fastp 軟件去掉reads 中的接頭序列和低質量的序列，從而得到高質量的測序數據(Clean reads)。使用Trinity 軟件在默認程序設置下對測序數據進行組裝，首先將reads 拼接成長的片段contigs，然后去除掉冗余contigs，將contigs 進行分組，連接成兩端不能再延長的轉錄片段(Unigenes)。

1.4 轉錄組生物信息學分析

將該次轉錄組測序獲得的所有轉錄本與6 大數據庫分別進行比對注釋分析。著重對GO 和KEGG 進行分析。GO 注釋包含3 個大分支，即細胞組分 (Cellular component,CC)、生物過程 ((Biological process,BP) 和分子功能 (Molecular function,MF)；所有已知基因序列都會被注釋到其中，隨著分支的展開，基因所具有的功能也越來越詳細的注釋出來。KEGG 通路可分為細胞過程、藥物開發、代謝、遺傳信息處理、環境信息處理、生物體系統和人類疾病7 大類別。隨著各類別的展開，對類別所包含的注釋基因描述就越清晰。采用FPKM 方法來測量基因的表達量水平。在比較不同時期的樣本后，使用DEGseq2 軟件包篩選一定條件下的差異表達基因，所設定的閾值為 |log2(FoldChange)|>1。按照通用默認值設定Padjust<0.05，則該基因為顯著差異表達基因。采用 Goatools 軟件對差異表達的基因進行GO 功能富集分析，當經過校正的Padjust<0.05 時，則認為此GO 功能類別存在顯著富集情況。KEGG 通路富集分析計算原理同GO 功能富集分析一致，即經過校正的Padjust<0.05 時，則認為此KEGG 通路存在顯著富集情況。

1.5 實時熒光定量PCR

隨機選擇6 個在玉足海參發育過程中表達量顯著變化的基因進行實時定量PCR 檢測，對轉錄組數據的可靠性進行驗證。用Trizol 試劑提取玉足海參幼體總RNA (上海百賽生物技術有限公司)，使用HiScript II Q RT SuperMix for qPCR (+gDNA wiper)(南京諾唯贊生物科技股份有限公司) 合成cDNA 第一鏈。按照 TB Green?Premix Ex Taq? II(上海百賽生物技術股份有限公司)說明書進行實時定量 PCR，并使用 SPSS 24.0 軟件 進行獨立樣本t檢驗，P<0.05 為顯著性差異，熒光定量檢測數據結果使用2—ΔΔCt方法計算分析。以β-Actin作為內參基因，所用引物見表1，所設計引物的擴增效率和特異性使用cDNA 樣品的標準曲線確定和依據熔解曲線分析確定。玉足海參的12 個樣品均進行檢測，每個樣品技術重復3 次。

表1 引物信息Table 1 Primer information

2 結果

2.1 測序數據統計與組裝結果

采集4 個不同發育時期的海參幼體樣品進行全長轉錄組測序分析。通過高通量測序后獲得了87.38 Gb Raw data。在經過原始數據過濾處理后，每個樣品的清潔讀數從42 974 554 到59 419 634 不等，有效讀數比率為98.9%。詳細統計結果見表2，包括原始堿基、有效堿基、Q20、Q30、GC 含量等。隨后通過Trinity 軟件進行denovo 組裝，共得到93 528 個Unigenes。

表2 轉錄本組裝結果統計分析Table 2 Statistics analysis of transcript assembly results

2.2 轉錄組注釋

將獲得的Unigenes 分別與GO、KEGG、COG、NR、Pfam、Swissprot 數據庫進行比對，獲得玉足海參幼體Unigenes 在各個數據庫的注釋信息。結果有48 087 個Unigenes 得到注釋，占整個組裝轉錄本的51.41%，GO、KEGG、COG、NR、Pfam 和Swissprot 數據庫中的同源序列注釋條數分別為36 422 (38.94%)、31 013 (33.16%)、36 608(39.14%)、47 228 (50.49%)、35 775 (38.25%)和35 567 (38.04%)。NR 物種注釋表明(圖2)，玉足海參與仿刺參 (Apostichopusjaponicus) 同源性最高(73%)，其次是綠海膽 (Lytechinusvariegatus)、紫色球海膽 (Strongylocentrotuspurpuratus)、長棘海星(Acanthasterplanci)和蝠海星(Patiriaminiata)。GO 功能注釋顯示基因功能聚類排名前三位的分別是細胞組分類別的細胞部分、分子功能類別的組合和催化活性 (圖3)。KEGG 通路注釋顯示功能基因參與排名前三位的生物學通路分別是環境信息處理類別的信號轉導、人類疾病類別的癌癥概述和神經退行性疾病(圖4)。

圖2 NR 數據庫比對統計圖Fig.2 Comparison chart of NR database

圖3 Unigene GO 功能注釋結果Fig.3 GO annotation results of unigene

圖4 Unigene KEGG 功能注釋結果Fig.4 KEGG annotation results of unigene

2.3 差異表達轉錄組分析

為探究玉足海參幼體變態發育過程中的基因表達變化情況，本實驗對海參的4 個不同發育時期的轉錄組測序文庫進行了比較分析，共檢測到43 720個差異表達基因，其中A 期與B、C、D 期比較，表達量差異顯著的基因數目分別為17 732 個 (上調9 863 個，下調 7 869 個)、17 189 個 (上調10 025個，下調 7 164 個)、24 081 個 (上調13 700 個，下調10 381 個)；B 期與C、D 期比較，表達量差異顯著的基因數目分別為11 757 個 (上調 6 229 個，下調 5 528 個)、20 990 個 (上調10 877 個，下調10 113個)； C 期與D 期比較，表達量差異顯著的基因數目為11 319 個 (上調6 727 個，下調4 592 個)。A期與D 期差異基因數最懸殊，表明玉足海參發育至樽形幼體后，基因表達數量發生顯著變化，使機體得以適應變態發育過程。相鄰組間 (A_vs_B,B_vs_C,C_vs_D) 上調和下調基因數量雖然基本均呈下降趨勢 (圖5)，但上調基因遠多于下調基因，說明在變態發育過程中，更多的基因被釋放表達。同時也表明一些基因只在特定時期的表達發生顯著變化，這可能與特定的生長發育表達調控相關。

圖5 相鄰組間差異表達基因火山圖注：A_vs_B.B 期以A 期為參照基準比較；B_vs_C.C 期以B 期為參照基準比較；C_vs_D.D 期以C 期為參照基準比較。后圖同此。Fig.5 Volcano map of differentially expressed genes between adjacent groupsNote: A_vs_B.Stage B was compared with Stage A; B_vs_C.Stage C was compared with Stage B; C_vs_D.Stage D was compared with Stage C.The same case in following figures.

2.4 GO 和KEGG 富集分析

為探究玉足海參響應變態發育的分子機制，對差異表達的基因進行GO 功能富集和KEGG 通路分析。相鄰組間 (A_vs_B,B_vs_C,C_vs_D) 差異基因分別富集于579、371、441 個GO 功能類別。其中，富集頻率排名靠前的主要是分子功能、細胞成分、催化活性等與細胞成長相關的生物學功能。KEGG 通路富集分析顯示A_vs_B、B_vs_C 和C_vs_D 的差異基因分別涉及342、340 和339 個通路，而且絕大部分被富集的信號通路與生長發育密切相關，包括細胞增殖分化相關通路 (細胞周期、基因復制等)、細胞凋亡通路 (癌癥途徑、Wnt等)、糖類代謝及脂類相關通路等 (圖6)。說明玉足海參變態發育是一個復雜的過程，受到多個基因的調控。

圖6 相鄰組間差異表達基因 KEGG 通路分析Fig.6 KEGG pathway analysis of differentially expressed genes between adjacent groups

2.5 與生長相關的基因分析

隨著玉足海參幼體的生長發育，其體型大小在變態附著階段會發生轉折。首先在耳狀幼體時期幼體體積會不斷膨大成樹突狀幼蟲 (圖1-c)，而后在樽形幼體時期幼體體積逐漸收縮成橢圓狀幼蟲(圖1-d)。對相鄰組間差異基因分析的結果顯示，與生長性狀調控相關的基因在特定發育時期顯著差異表達。如隨著幼體體積的變化，cox1和tll1基因表達量在A 至C 期逐漸上調，然后在D 期顯著下調。ctsl2基因在B 期顯著上調而在D 期顯著下調，c-lectin2基因僅在D 期顯著上調，survivin基因只在D 期顯著下調。這說明機體依據變態階段不同生長發育時期的反饋任務去執行相應的表達程序，控制有效調節因子去執行細胞增殖、分化和凋亡功能，使其能夠正常渡過變態附著階段[21-22]。

2.6 實時熒光定量PCR 分析

隨機挑選6 個差異表達顯著的基因，對目標轉錄本進行定量分析以驗證玉足海參轉錄組測序結果數據和分析的準確性。結果如圖7 所示，實時熒光定量PCR 分析得到的表達趨勢與Illumina 測序得到的表達趨勢相似，但幾個基因的表達差異幅度存在一定差異。這可能是兩種方法的靈敏度不同所致。據報道，RNA-Seq 對目的評估基因更為敏感，說明使用測序分析所驗證的差異表達轉錄本較為可靠[23]。

圖7 qRT-PCR 驗證 RNA-seq 結果Fig.7 RNA-Seq results verified by qRT-PCR

3 討論

3.1 注釋信息分析

玉足海參幼蟲發育需經歷浮游生活到附著生活的轉變，這個過程受到環境和遺傳因素的共同影響[24-25]。本研究對組裝得到的Unigenes 進行生物學功能注釋，共檢測到 43 720 個差異表達基因。通過NR 物種相似度功能注釋分析，有三分之一的Unigenes 與庫中已知基因同源，相似度排名靠前的均為海洋無脊椎動物，玉足海參與這些物種具有較高的親緣性，尤其是以仿刺參、海星等為代表的棘皮動物相似度最高；它們的許多基因是由共同祖先基因分化演變而成。本實驗所獲得的NR 序列比對呈現出注釋率數值較小，這可能是由于國際公共數據庫收錄的相近物種基因信息較為匱乏，并且也受到實驗條件和測序水平等因素的限制。另外有一小部分序列未注釋匹配成功，可能原因是這些序列與現有數據庫中的基因序列匹配度較低，可能為無翻譯能力的RNA 或者無表達功能的結構域序列[26]；它們有可能參與重要的生命活動，并且在生物體內發育功能和過程中起重要作用。對3 個相鄰組間GO 注釋富集分析 (表3)，其中，生物過程最主要的GO 類別聚集是細胞過程、代謝過程，細胞組分最主要的GO 類別聚集是細胞部分、膜部分，分子功能最主要的GO 類別聚集是組合、催化活性。不同比較組間轉錄本GO 類別與功能富集多數與細胞的完整性及其功能相關，這意味著隨著生長周期的延長，細胞內的碳水化合物、脂質等活性物質含量逐漸積累，細胞數量和體積也逐漸增長，相關轉錄本的功能多與細胞周期調控、細胞組合及催化活性緊密相關。在大耳到樽形幼體生長模式的轉換過程中，細胞黏附類別富集頻率顯著增加，其生物學功能是通過細胞的識別與黏附使具有相同表面特性的細胞聚集在一起形成組織和器官，這很可能是幼體從浮游期到附著期特定器官生長變化的關鍵功能。同樣的，對相鄰比較組間轉錄本參與的KEGG 注釋通路進行富集分析，發現其注釋多與神經和肌肉發育代謝過程相關。這說明在變態附著發育期間，很多控制神經、肌肉和皮膚等組織表達通路里的轉錄因子很可能是影響生物體正常發育的關鍵因素，如肌動蛋白細胞骨架調節、ECM-受體相互作用、鈣等信號通路，這有待更為深入的研究。

表3 相鄰組間差異表達基因 GO 富集分析 (前 5)Table 3 GO enrichment analysis of differentially expressed genes between adjacent groups (Top five)

3.2 變態附著階段差異基因分析

玉足海參在浮游幼體期主要是以海洋藻類等浮游有機質為食，大耳幼體后期變態為樽形幼體附著后轉而以淺層沉積物質為主要攝食對象。因此，這兩種截然不同的生活方式可能涉及到兩套非常不同的基因調控網絡[27]。對4 個時期幼體觀察研究發現，樽形幼體時期的死亡率達到54%，這在玉足海參變態發育以及整個生長過程中為最高[3]。然而，目前針對這一階段的基因調控基礎仍不清楚，基因層面的機制受控很可能使變態延遲或變態失敗，最終導致幼蟲可能因能量消耗過多和環境變化而死亡[28]。分析4 個時期的差別，A 和D 期之間的差異基因數目最多 (24 081 個)，而C 和D 期之間的差異基因數目最少 (11 319 個)。表明有較多的調控基因參與了幼體發育過程，多種機制的協同與拮抗作用反饋影響幼體發育的相關基因的激活與抑制。此外，相鄰時期表達的差異基因數量逐漸減少且放緩。這表明變態發育引起的玉足海參體內的反應從最開始不穩態到過渡新穩態，至樽形幼體時期達到一個新的平衡。分析比較組間基因表達量差異結果，有高FPKM 值的基因與能量代謝、蛋白質代謝、催化活性、細胞骨架等基本過程相關。它們在體內干細胞正常發育、組織維持和血管生成中可能發揮了重要作用。如在B 期顯著上調和D 期顯著下調的ctsl2基因，它是一種重要的半胱氨酸蛋白酶，與溶酶體蛋白酶類似，能夠清除細胞中多余的蛋白質[29]；研究發現，ctsl2可以使仿刺參三倍體的生長速度快于二倍體[30]。此外，也有一些基因雖然并無顯著差異的表達，但在變態附著階段中可能起到至關重要的作用。如作為執行凋亡程序的caspase-3基因，在很多生物體內的細胞凋亡過程中起重要作用，在玻璃海鞘 (Cionaintestinalis)幼體變態發育的研究中，capsase-3依賴性細胞凋亡是通過磷酸化 ERK 執行了尾部的退化過程[31-33]，其在玉足海參D 期的表達水平較高。還有作為負調節動物骨骼肌發育和生長因子的肌肉生長抑制素(myostatin) 基因，能通過抑制成肌細胞增殖和多核肌管分化來調節肌肉質量[34-36]，其在玉足海參A 和D 期的表達水平較高。因此，差異顯著基因可能在玉足海參變態附著的過程中，對機體適應新的生長模式起非常重要的作用。它們可以作為玉足海參選育潛在的候選基因，推動整個幼體繁育的生物學過程研究。

3.3 關鍵發育信號通路分析

在本次轉錄組測序中，KEGG 注釋顯示玉足海參4 個不同發育時期幼體的功能基因涉及細胞周期、Wnt、PI3K-Akt、磷酸戊糖、胰島素信號轉導、蛋白質消化吸收、膽固醇代謝、MAPK 信號通路等8 個與變態發育密切相關的信號通路。如細胞周期信號通路中的cdk、cyclin基因家族網絡調控生長發育在其他水生動物的研究中已有報道，cdk基因通過與cyclin基因結合形成復合物，從而調節細胞增殖和分化。對文蛤 (Meretrixmeretrix)的研究發現，cdk基因在各實驗組的表達情況與其殼長生長正相關，推測cdk基因在文蛤早期苗種生長中發揮了重要作用[37]。與胚胎發育相關的Wnt 信號通路也有相關報道，在對半球美螅水母(Clytiahemisphaerica) 研究中，wnt3基因通過影響受體基因的表達量來維持其胚胎發育的體軸模式系統中口端和反口端的穩定[38]。在對貝螅(Hydractinia) 的研究中，通過分析frizzled和β-catenin這兩個Wnt 信號通路關鍵基因，以及gsk-3基因的阻斷，證實了經典Wnt 信號通路參與了貝螅功能干細胞的調控過程[39]。而PI3K-Akt、磷酸戊糖等通路的基因對生長發育的調控在其他生物中也有報道[40-42]，它們在細胞的存活、分化、生長、運動和凋亡等多種生理和病理過程中起到重要作用。這些通路中的基因在玉足海參變態發育過程中如何調控關鍵生物學過程有待進一步研究。

此外，根據實驗觀察 (圖1)，玉足海參在大耳幼體變態成樽形幼體的過程中，顏色逐漸加深、體型逐漸變小，這個過程中的通路調控網絡主要是控制幼體發生了細胞凋亡、增殖分化等生命活動，最終大耳幼體在這些程序正常執行之后變成了樽形幼體。此龐大的運作機制是一個最為復雜且極為重要的步驟，這可能是導致其高死亡率的主要原因。對這一時期富集通路的分析發現，癌癥途徑發揮了關鍵作用 (圖6)，其富集注釋基因數量顯著上升，富集頻率在這一時期也居首位。它涉及到很多生長調控和凋亡基因，正常的細胞凋亡會抑制致癌轉化，并控制器官形成和免疫反應。如此途經的ap-1基因是信號傳導的第三信使，在D 期顯著上調表達。它的基因產物能夠誘導D N A 與保守的bZIP 結構結合，其活性在信號傳導過程中會被MAPK 信號通路的磷酸化調控[43]。通過影響相關免疫、生長基因的轉錄表達過程，進而參與生物體內的免疫系統防御、細胞增殖分化、編程性死亡等過程[44]。該基因在青蛤 (Cyclinasinensis) 受到免疫刺激后顯著上調，其不僅影響免疫應答過程，同時也可能參與了細胞的編程性死亡[43]。該基因也參與了青海湖裸鯉 (Gymnocyprisprzewalskii) 肌肉生長和分化過程，其在肌肉中的表達量顯著上調[45]。此外，survivin基因作為凋亡蛋白抑制因子也在這時期顯著下調，survivin基因能抑制fas、bax及caspase3/7基因誘導的細胞凋亡[46]。這說明癌癥途徑參與了細胞凋亡響應的重要過程，具體在玉足海參變態附著階段如何調控生物學過程也有待進一步研究。

4 結論

本實驗通過轉錄組分析研究了玉足海參變態附著階段的分子調控機理，結果顯示玉足海參幼體生長發育和細胞凋亡發生的信號轉導途徑具有多樣性，其差異基因和通路彼此間相互促進或制約形成了錯綜復雜的調控網絡。針對玉足海參附著階段死亡率過高的問題，本研究推測主要與細胞是否正常執行凋亡相關，對這些相關差異基因和通路進行內在機理的探索，有助于揭示玉足海參的生長發育規律，為后續開展玉足海參的規?；庇夹g研究奠定基礎，同時也可為其他熱帶海參的人工繁育提供參考。