三角帆蚌育珠蚌群體生長性狀和轉錄組差異分析

張愛菊，張根芳，顧志敏，周志明

( 1.浙江省淡水水產研究所，中國水產科學研究院東海水產研究所浙江研究中心，農業農村部淡水漁業健康養殖重點實驗室，浙江省淡水水產遺傳育種重點實驗室，浙江 湖州 313001；2.金華職業技術學院，浙江 金華 321000 )

三角帆蚌(Hyriopsiscumingii)廣泛分布于長江、淮河等流域，是我國特有的河蚌資源，又是淡水珍珠養殖的優良品種。然而, 由于現有的供片蚌和育珠蚌質量參差不齊以及養殖群體種質退化等問題日益突出，使得淡水珍珠產業可持續發展受到制約。因此，開展三角帆蚌的遺傳改良和品種選育對淡水珍珠產業有著極為重要的意義。迄今為止，國內關于三角帆蚌不同地理種群間形態學及其育珠性狀的比較研究等方面已有較多報道[1-4]，均表明三角帆蚌不同地理種群的形態及育珠相關性狀差異明顯。然而，不同三角帆蚌育珠蚌選育群體形態及育珠相關性狀的差異分析尚未有報道。

轉錄組學是從RNA水平上研究基因表達情況的學科[5]，它將基因組學研究帶入了高速發展的嶄新時代[6]。轉錄組學可用于分析不同組織或生理狀態下某些基因表達水平的差異，從而發掘與特定生理功能相關的未知基因，具有信噪比高、分辨率高、應用范圍廣、價格低廉等優勢。目前，用于第二代高通量測序的技術主要有3種：Roche公司的454測序、Illumina公司的Solexa測序以及ABI公司的SOLiD測序[7]。目前，軟體動物中大多數研究主要集中于利用Roche-454測序或Illumina測序技術研究海洋物種的轉錄組學研究[8-10]，而在淡水育珠貝類中，僅見利用Roche-454測序技術研究紫色和白色珍珠層三角帆蚌組織中轉錄組表達差異[11]和利用Illumina測序技術研究三角帆蚌組織總RNA轉錄組[12]的報道。

為探究3個已經過2代選育無顯著形態學差異的三角帆蚌育珠蚌選育群體——鷹蚌、超長蚌和縱紋蚌群體幼蚌之間的差異，筆者運用聚類分析、主成分分析兩種多元統計方法比較分析這3個群體幼蚌的生長性狀指標，并利用Illumina測序技術比較鷹蚌和超長蚌群體肝胰腺組織轉錄組表達的差異情況，以期篩選出具有明顯生長優勢、適合作為插核或插片的育珠蚌群體，為三角帆蚌優良育珠品種的選育和育珠實踐提供科學指導。

1 材料與方法

1.1 試驗材料

以鷹蚌、超長蚌、縱紋蚌3個三角帆蚌幼蚌選育群體為研究對象，于2014年10月在同一養殖池塘中隨機抽取這3個選育群體各40只，暫養于試驗池塘，2 d后進行測量。取樣時，暫養池內水溫23 ℃、pH 7.69、透明度30 cm、亞硝酸氮質量濃度4.98 mg/L、氨氮質量濃度1.38 mg/L。2015年3月再次于同一養殖池塘取樣，隨后暫養，方法同上所述。此時暫養池內水溫10 ℃、pH 7.64、透明度30 cm、亞硝酸氮質量濃度3.05 mg/L、氨氮質量濃度1.16 mg/L。兩次取樣的幼蚌均為2014年春季同批繁殖，且均取自浙江金華威旺養殖新技術有限公司試驗基地。3個群體均從良種場的自然養殖群體中通過家系選育的方法篩選而出，且已經過2代選育。

1.2 生長性狀差異分析

1.2.1 數據測量

形態學測量參數為殼長、殼寬、殼高等11個可量性狀，采用游標卡尺等測量工具(精確至0.01 mm)，測量方法參照文獻[1]的方法，測量部位與參數見圖1。測量性狀有：殼長(AB)、殼高(OH)、殼寬(WI)、全高(FG)，以及殼頂至前端長(OA)、殼頂至后端長(OB)、殼頂至鉸合部前緣長(OC)、殼頂至鉸合部后緣長(OD)、殼頂至殼頂上方突起長(OE)、帆狀頂點至前端長(FA)、帆狀頂點至后端長(FB)。質量參數采用電子天平等工具測量(精確至0.1 g)，為濕質量、殼質量、內臟團質量等3個參數。

圖1 形態學測量位點

1.2.2 數據計算

為消除蚌體規格大小對參數值的影響，將每只蚌的所有形態學參數分別除以殼長(AB)予以校正，得出10個形態學比例性狀。同時，參考文獻[4，13]的方法計算體質量指數(IBW)、殼質量指數(ISW)、肥滿度(CF)和殼寬指數(IST)。按下式計算每個選育群體的絕對生長率和相對生長率：

ρ1=(L2-L1)/t

ρ2=(L2-L1)/L1/t×100%

式中，ρ1為絕對生長率(cm/月或g/月)，ρ2為相對生長率(%/月)，L1為第一次采樣時殼長、殼寬、殼高或體質量(cm或g)，L2為第二次采樣時殼長、殼寬、殼高或體質量(cm或g)，t為采樣間隔月數(月)。

1.2.3 數據處理

應用Excel 2007與SPSS 16.0軟件進行數據處理。數據表示為平均值±標準差，以P<0.05為差異顯著，P<0.01為差異極顯著。對10個形態學比例性狀參數與4個指數參數進行統計學分析。采用的多元統計方法主要為聚類分析、主成分分析。其中，聚類方法為歐式距離的最短系統聚類法[1,14]。

1.3 轉錄組差異表達分析

1.3.1 組織樣本獲取

形態學測量參數和質量參數測量完成后，挑選2015年選取的批次中殼體完整、活力好的鷹蚌和超長蚌2個群體幼蚌各9只，清除其體表附著生物及污物，活體解剖，迅速取肝胰腺組織，用液氮速凍后分別放在對應的群體凍存管中，存于-80 ℃冰箱中備用。

1.3.2 總RNA的提取

試驗用槍頭在DEPC-H2O中過夜后高壓滅菌，置于100 ℃烘箱中烘干備用。采用TRIzol總RNA提取試劑(Invitrogen, USA) 分別提取2個選育群體肝胰腺樣本的總RNA。

1.3.3 轉錄組測序

參照文獻[12]的方法，由上海美吉生物醫藥科技有限公司提供測序服務。

1.3.4 測序數據處理及分析

Illumina HiseqTM2000測序得到的原始圖像數據處理參照文獻[12]的方法進行質控后，使用TopHat軟件分別將鷹蚌組和超長蚌組質控后得到的clean reads與文獻[12]的de novo組裝結果比對和組裝。組裝好的結果按照文獻[12]的方法進行功能注釋與功能分類。

根據每個基因/轉錄本在樣本中的FPKM值，分析差異表達分析。對篩選出的差異表達基因，使用EdgeR軟件對所有基因在各組樣本中的表達進行差異顯著性分析，篩選標準為：FDR(False Discovery Rate)<0.05，并且︱logFC︱≥1。使用Goatools軟件進行GO富集分析，使用KOBAS軟件進行KEGG富集分析,滿足條件FDR<0.05所對應的GO條目和KEGG代謝通路即為顯著富集的GO條目代謝通路和KEGG代謝通路。

2 結 果

2.1 3個選育群體間生長性狀比較

2.1.1 性狀參數比較

方差分析結果顯示，殼頂至鉸合部后緣長/殼長、殼頂至殼頂上方突起長/殼長在3個群體間均無顯著性差異(P>0.05)；鷹蚌群體的殼寬/殼長為(0.23±0.003)，顯著高于超長蚌的(0.20±0.018)和縱紋蚌的(0.21±0.003)(P<0.05)，超長蚌與縱紋蚌殼寬/殼長之間差異不顯著(P>0.05)；殼高/殼長、全高/殼長、殼頂至前端長/殼長、殼頂至后端長/殼長、殼頂至鉸合部前緣長/殼長、帆狀頂點至前端長/殼長在3個群體間差異顯著(P<0.05)，均呈現縱紋蚌>鷹蚌>超長蚌的趨勢；鷹蚌群體帆狀頂點至后端長/殼長顯著高于縱紋蚌和超長蚌(P<0.05)。綜合分析可知， 鷹蚌群體偏寬型，縱紋蚌次之；縱紋蚌群體偏高型，鷹蚌次之；超長蚌偏長型，鷹蚌次之。

分析結果顯示，殼質量指數、殼寬指數、體質量指數、肥滿度在3個群體間差異顯著(P<0.05)，鷹蚌群體的殼質量指數(21.26±0.402)和體質量指數(68.09±0.583)均低于超長蚌的(28.59±0.895，74.76±0.991)和縱紋蚌的(22.19±0.383，68.27±0.638)群體的(P>0.05)，而其殼寬指數(0.13±0.001)顯著高于超長蚌的(0.12±0.001)和縱紋蚌的(0.12±0.001)(P<0.05)，且其肥滿度(36.08±0.481)高于超長蚌的(33.76±0.577)和縱紋蚌的(32.15±0.527)群體(P>0.05)。故推測鷹蚌群體具有相對較肥滿的內臟團，且在殼寬生長上具有明顯優勢，可能為更適合內臟囊插核的三角帆蚌養殖群體。

2.1.2 聚類分析

對所有樣本14個參數進行聚類分析得出(圖2)，鷹蚌群體與縱紋蚌群體首先聚在一起，表示鷹蚌與縱紋蚌的距離較短，在形態上可能更為接近。

圖2 3個三角帆蚌選育群體聚類分析

2.1.3 主成分分析

主成分分析結果顯示，只有前4個特征根大于1，因此SPSS只提取了前4個主成分。第一主成分的總方差貢獻率為28.2%，第二主成分為16.74%，第三主成分為16.68%，第四主成分為8.61%，貢獻率累積達到70.24%(表1)。4個主成分的累計貢獻率較高，基本可以反映不同選育群體之間的生長差異。由表1可見，主成分1在體質量指數、殼質量指數和殼高/殼長上有較大的載荷，主要從體質量和殼高上反映形態和育珠性能；主成分2在殼寬指數、殼寬/殼長上有較大載荷，從殼寬方面反映形態和育珠性能。主成分3在帆狀頂點至前端長/殼長上有較大載荷，表現為三角帆蚌縱向的性狀方面；主成分4在全高/殼長、殼頂至后端長/殼長上有較大載荷，也主要表現為三角帆蚌縱向的性狀方面。

表1 4個主成分的載荷值和方差貢獻率

注：*表示載荷值>0.700.

3個育珠蚌選育群體第1、2主成分的散布情況見圖3。由圖3可見，3個群體均有較多的重疊區域，尤其是鷹蚌群體和縱紋蚌群體，進一步驗證其形態更為相近。

圖3 3個三角帆蚌選育群體第1、2主成分散布

2.1.4 生長速率比較

超長蚌群體和鷹蚌群體的殼長絕對生長率顯著高于縱紋蚌群體(P<0.05)，但鷹蚌群體和超長蚌群體之間差異不顯著(P>0.05)；鷹蚌群體的殼寬絕對生長率顯著高于其他2個群體(P<0.05)，而其殼高絕對生長率顯著高于縱紋蚌群體(P<0.05)，但與超長蚌群體差異不顯著(P>0.05)；鷹蚌群體和超長蚌群體的體質量絕對生長率極顯著高于縱紋蚌群體(P<0.01)(圖4)。同時，3個選育群體的殼長、殼寬、殼高和體質量的相對生長率呈現與對應參數的絕對生長率相同的趨勢(圖5)。因此，綜合分析認為，鷹蚌群體在殼寬和殼高生長上具有明顯優勢，超長蚌群體在殼長生長方面具有明顯優勢。

2.2 轉錄組差異分析

2.2.1 測序產量統計及質量評價

測序共計獲得61.92×106條原始序列，堿基數達到62.5×108bp，近94.20%的序列測序質量值均在Q20以上。質控后獲得58.09×106條高質量序列，堿基數56.79×108bp，近98.48%的序列測序質量值均在Q20以上。

2.2.2 de novo 序列組裝結果

利用Trinity軟件進行de novo拼接，將58.09×106條高質量序列拼接成92 347條非冗余unigene。unigene長度為351～28 861 bp，平均長度1150.61 bp。

圖4 3個三角帆蚌選育群體的絕對生長率變化Tukey Test分析，柱狀圖上同個參數中不同小寫字母的表示群體組間差異顯著(P<0.05)，不同大寫字母的表示群體組間差異極顯著(P<0.01)，下同.

圖5 3個三角帆蚌選育群體的相對生長率變化

2.2.3 unigene功能注釋

利用GO數據庫，可以將基因通過GO功能分類，共有11 174條unigenes獲得了GO信息，歸為生物學過程，分子功能和細胞組件3個分支。細胞組件分支中，歸類到細胞組分和細胞條目的unigenes最多，均為5181條，分別占22.47%。其次為細胞器和細胞器組分，分別為3670條(15.92 %)和2168條(9.40%)；分子功能分支中，結合功能的unigenes最多，5772條(45.27%)，其次為催化活性(5262條，41.27%)和轉運活性(494條，3.87%)；生物學過程分支中，歸為細胞生理過程的unigenes最多(17.78%)，代謝過程次之，占14.82%。

2.2.4 差異表達基因的篩選

鑒于鷹蚌群體和超長蚌群體明顯的殼體生長優勢，主要對這兩個群體的轉錄組差異表達進行闡述。以超長蚌每個基因在樣本中的FPKM值為對照，共篩選出差異表達unigenes 63 098條，其中顯著性差異表達基因388個。與超長蚌群體相比，鷹蚌群體中顯著上調的基因有244個，顯著下調的基因有144個，分析顯示，大部分基因與代謝、免疫功能等生物學過程相關。

2.2.5 GO功能富集分析

GO功能富集分析可以給出兩樣本間差異表達基因所參與的主要生物學過程。通過GO功能富集分析，共有3571個差異表達基因富集上了140個GO號，其中顯著富集的GO號80個，歸類到生物學過程類別的GO號41個，分子功能類別的GO號20個，歸類到細胞組分類別的GO號19個。分析結果顯示，顯著富集的GO條目主要參與各類物質代謝過程的調節、生物合成過程的調節以及基因表達的調節等生物學過程。

2.2.6 KEGG通路富集分析

KEGG代謝通路富集分析可以給出兩樣本間差異表達基因所參與的最主要代謝通路。本次研究結果顯示，2個選育群體中的差異表達基因共富集上了241條通路，其中注釋到癌癥通路中差異表達基因最多，共有46個；其次為泛蛋白介導的蛋白酶解通路、HTLV-I感染通路和弓形蟲(Toxoplasmagondii)病通路，差異表達基因分別為45、40、38個。KEGG 通路富集分析結果顯示，共有2條顯著富集的通路，為半胱氨酸和蛋氨酸代謝通路和檸檬酸循環通路，注釋到前者通路中差異表達基因為 12個，而注釋到后者通路中差異表達基因為13個。進一步分析后得出，這2條顯著富集的通路分別屬于代謝通路中的氨基酸代謝和糖代謝通路。

3 討 論

3.1 3個群體幼蚌生長性狀差異比較分析

錢榮華等[1]利用多元分析方法對我國5大淡水湖的三角帆蚌進行分析，從形態學上探討了地理種群間的差異；董志國等[2]運用3種多元分析方法對3個地理種群雜交與自交后的9個群體的10個形態比例性狀進行了比較研究；聞海波等[3-4]先后對長江和淮河兩個水系的三角帆蚌親本形態特征及3個地理種群三角帆蚌形態等性狀進行了比較。這些研究均表明，三角帆蚌不同地理種群的形態及育珠相關性狀差異明顯。

本次試驗的三角帆蚌是自洞庭湖、鄱陽湖等天然水域中引進的親本，經過兩代選育形成的3個不同家系的群體，在形態上既相似又有一定程度的差異，成熟鷹蚌個體腹緣具有鷹嘴狀凹刻，成熟超長蚌個體的殼長/殼高明顯大于鷹蚌和縱紋蚌，而縱紋蚌的成熟個體貝殼上則具有顯著的縱紋。但這3個群體的幼蚌則無類似明顯的形態學特征。為此，采用多元統計分析方法，對其生長性狀參數進行了比較分析。聚類分析結果中，鷹蚌群體與縱紋蚌群體距離最短，形態更為接近，而超長蚌群體相對這2個群體形態差異較大。主成分散布圖進一步驗證了這一結果。通常，形態特征是由遺傳因子與環境因子共同作用的結果，本試驗的3個三角帆蚌均采自同一地區的同一池塘，基本可以排除環境因子因素。究竟是什么遺傳因素導致這些群體之間形態的變異，以及這些群體之間的形態變異究竟多大程度地反應在遺傳物質與等位基因的變異上，還有待于進一步的深入研究。

3.2 貝殼形態、生長性狀與選育

珍珠質量除受養殖環境因子影響外[15]，還與多種因素相關，而其中育珠蚌親本選擇是一個重要因素。Wada[16]認為，育珠蚌親本的選擇應主要以珍珠層顏色、貝殼形態及生長等指標為主。白志毅等[17]認為，三角帆蚌選育應以殼質量(或體質量)和殼寬為主要選育目標，李夢軍等[18]群體選育結果也證明了這一觀點。由于本研究的目的是篩選出具有明顯生長優勢、適合作為插核或插片的育珠蚌群體，不考慮珍珠顏色的選育，故認為三角帆蚌親本選育應以貝殼形態及生長等指標為主。鑒于鷹蚌和超長蚌群體幼蚌較快的生長速率，主要考慮貝殼形態方面的因素。本研究對3個三角帆蚌選育群體幼蚌的貝殼形態測定與分析結果表明，鷹蚌幼蚌的殼寬指數和肥滿度均高于超長蚌和縱紋蚌，因此鷹蚌群體的育珠性能優于其他2個群體，推斷其作為育珠蚌更適合內臟囊插核，將有利于大顆粒優質有核珍珠的培育。同時，從貝殼形態上看，超長蚌為長扁型，在殼長生長方面具有明顯優勢，作為育珠蚌更適合外套膜插片，可以增加插入的細胞小片數量，有利于無核珍珠的培育。

3.3 高通量測序技術與選育

近年來，運用高通量測序技術輔助水產經濟物種的遺傳選育工作在我國越來越受到重視。目前，對于三角帆蚌的轉錄水平研究主要采用cDNA文庫/EST測序技術[19]、cDNA文庫/抑制消減雜交技術[20]、cDNA文庫/Roche測序技術[11]、Illumina/Solexa測序技術[12]，但這些研究主要集中于三角帆蚌整體水平或者珍珠顏色選育等方面，而關于有核珍珠/無核珍珠育珠蚌選育方面的研究略有匱乏。

目前，軟體動物外套膜轉錄組高通量測序的研究較多，且主要用于生物礦化分析[8,11]。然而，鑒于本研究的目的以及鷹蚌和超長蚌各自的殼體生長優勢，同時考慮肝胰腺在軟體動物眾多生命大分子的合成及分解代謝等方面的重要功能[21]，很大程度上影響其生長發育等因素，本研究進一步應用高通量測序技術對同齡鷹蚌和超長蚌2個選育群體的幼蚌肝胰腺轉錄組進行深度測序。結果發現，篩選出差異表達unigenes 63 098條，其中顯著性差異表達基因388個，顯著富集的GO條目主要參與各類物質代謝過程的調節、生物合成過程的調節以及基因表達的調節等生物學過程，顯著富集的KEGG通路為半胱氨酸和蛋氨酸代謝通路和檸檬酸循環通路，均屬于代謝通路。故推測這2條代謝通路及其中的差異表達基因可能是造成鷹蚌和超長蚌2個選育群體生長性能差距的主要遺傳因素。

以上測序結果和分析為三角帆蚌的遺傳育種提供了豐富的分子生物學信息，有利于今后的家系選育工作。然而，這些結果僅是在轉錄組水平生物信息學分析上的一種推測，今后仍需進一步進行驗證。