基于線粒體cytb和d-loop基因序列的東南沿海可口革囊星蟲遺傳多樣性分析

許瑞雯，陳興漢，余祥勇

（1.陽江職業技術學院，廣東陽江 529566；2.廣東海洋大學，廣東湛江 524088；3.華南農業大學，廣州 510642）

可口革囊星蟲（Phascolosoma esculenta）為星蟲動物門（Sipunculoidea）革囊星蟲目（Phascolosomaliformes）革囊星蟲科（Phascolosomatidae）革囊星蟲屬動物，俗稱“泥蟲”、“土筍”或“泥丁”，是中國特有的星蟲種類，在我國東南各省區的沿海灘涂均有分布［1］。可口革囊星蟲作為海洋珍品之一，富含蛋白質和多種不飽和脂肪酸等營養物質，被稱為“動物人參”，是當地居民的特色經濟水產生物。近年來我國東南沿海開發建設進程加快，環境污染、灘涂圍墾、超負荷采捕等現象屢禁不絕，導致中國的可口革囊星蟲野生種群種質資源量銳減［2］，市場價格連年攀升。因此相關的資源保護、種苗繁育及增養殖開發引起行業關注［2］，相關學科學者對可口革囊星蟲的人工養殖和苗種繁育做了大量研究［3-6］。然而關于可口革囊星蟲種群遺傳學的研究報道相對較少見，為了該物種資源的合理開發及可持續利用，對其野生群體的遺傳多樣性開展研究十分必要。

自20世紀60年代線粒體DNA（mtDNA）被發現后［7-8］，mtDNA作為核外遺傳物質與核基因組DNA相比，由于其具有母性遺傳、無重組、分子結構簡單、進化速率快、無組織特異性和提取方便等優點［9］，已成為研究近緣種和種內群體間遺傳分化的有力工具［10-12］。其中，線粒體細胞色素b（cytb）和線粒體控制區（d-loop）兩個基因序列是分子標記技術中較常使用的兩個基因序列，其遺傳多樣性在雙齒圍沙蠶（Perinereis aibuhitensis）［13］、裸體方格星蟲（Sipunculus nudus）［14-15］、日本囊對蝦（Penaeus japonicus）［16］、厚殼貽貝（Mytilus coruscus）［17］等海洋無脊椎動物的研究中已被應用。

本研究以線粒體DNA的cytb和d-loop基因序列作為分子標記，對我國東南沿海5個可口革囊星蟲野生地理群體開展遺傳多樣性和遺傳結構分析，以期為這一中國特有的星蟲種類提供人工繁殖和育苗的基礎資料，并為后續研究保護該物種的種質資源及其合理開發利用提供數據。

1 材料與方法

1.1 樣品信息

實驗所用的可口革囊星蟲樣品均采集于2019年，樣品采集地點分別為廣東湛江（ZJ）、廣東陽江（YJ）、廣西欽州（QZ）、福建寧德（ND）、浙江溫州（WZ）5個地方的潮間帶灘涂，具體采集信息見表1。可口革囊星蟲樣品均為野生自然個體，系根據該物種的主要形態學特征“大部收吻肌融合，背腹兩對相距較遠，每對收吻肌的起點較遠；蟲體顏色為淡黃色或棕色，蟲體形狀圓長，類似長頸瓶或燒瓶［1］”進行樣品鑒定。所用樣品均在采樣地進行預處理：每個群體采樣100個以上，然后每組樣品隨機選取30個解剖，取其體壁肌肉以95%乙醇固定、保存。

表1 東南沿海可口革囊星蟲群體采集信息Tab.1 Information of sampling locations of P.esculenta in southeastern coastal areas

1.2 DNA提取、擴增和測序

取可口革囊星蟲體壁肌肉約100 mg剪碎放入1.5 mL的EP（eppendorf）管。采用Universal Genomic DNA Kit（CW2298M，通用型柱式基因組提取試劑盒，康為世紀生物科技股份有限公司）對可口革囊星蟲樣本進行粗提，然后將提取的DNA樣本通過核酸蛋白測定儀（柏點，柏精-BIODROPμLITE，美國）檢測，選取DNA濃度大于50 ng·μL-1的樣品DNA置于-20℃低溫保存，用于后續實驗。

根據可口革囊星蟲已發布的線粒體基因（GenBank登錄號為NC_012618），采用Primer Premier5.0軟件［18］進行線粒體cytb序列和d-loop序列的引物設計。其中cytb序列的正、反向擴增引物分別是：PE-F4645（TGGATGACTACTACGATTC）和 PE-R5475（TGGCCTAGTCAGTCATATGG）；d-loop序列的正、反向擴增引物分別是：PE-F10225（ACTGAATATTCTAAAGTGGCAG）和PE-R11025（GTGGAAGTCCTATTAGGCTAAC）。以上引物均由廣州瑞科基因科技有限公司合成。

根據2×TaqMaster MixPCR擴增試劑盒（code：CW0682L，康為世紀生物科技股份有限公司）的說明書要求和步驟，配置反應體系40μL，其中2×TaqMasterMix（Dye）20μL，正反引物（濃度10μmol·L-1）各2μL，模板DNA 2μL，加滅菌水使反應體系至40μL。PCR反應程序為：94℃預變性3 min；98℃變性15 s，50℃退火30 s，68℃延伸1 min，共35個循環；最后68℃延伸5 min。PCR擴增產物經1%瓊脂糖凝膠電泳檢測，選出顯示亮帶的可口革囊星蟲cytb和d-loop擴增片段進行純化測序。

1.3 實驗數據分析

測序后的序列使用DNAStar軟件包中的EditSeq程序讀取，以SeqMan程序對序列進行拼接、編輯，然后用MegAlign程序對序列進行比對并輔以人工校正。修正后的cytb和d-loop序列采用DnaSP 6.0軟件［19］對5個地理群體的可口革囊星蟲樣本分別進行4種堿基構成及數量、多態位點、單倍體數量等運算分析，并檢測單倍型在各群體中的分布。

利用MEGA 7.0軟件［20］基于K2-P雙參數模型計算可口革囊星蟲5個群體兩兩之間的群體間遺傳距離，并構建單倍型鄰接樹（1 000次自展重復）。

利用Arlequin 3.5.2.2軟件［21］計算5個群體的基因多樣性指數（H）與核苷酸多樣性指數（π），并進行分子生物學方差分析（AMOVA分析），估算種群內和種群間遺傳變異的分布情況；利用Tajima’s D和Fu’s Fs檢驗來分析可口革囊星蟲群體的歷史動態。

2 結果與分析

2.1 可口革囊星蟲群體的序列差異及堿基組成

可口革囊星蟲的cytb和d-loop基因序列通過PCR擴增，并去除低質量堿基及非目標基因片段后，用于測序的cytb序列堿基長度主要為576 bp，而d-loop序列堿基長度主要為385 bp和384 bp兩種，其中大部分是堿基序列，長度為385 bp，占序列總數的75.70% 。通過對序列進行抽樣分析可知，序列堿基長度存在多態性的原因是空白位點的插入。經統計，4種堿基含量在可口革囊星蟲5個地理群體的cytb和d-loop序列中差異小：cytb的A含量為28.11%～28.17%，G含量為12.70%～12.73%，C含量為25.53%～25.61%，T含量為33.50%～33.62%；d-loop的A含量為38.22%～39.30%，G含量為11.44%～11.58%，C含量為18.76%～18.87%，T含量為31.37%～31.45%。其中兩個基因序列均表現出明顯的A+T含量高于C+G含量，有明顯的堿基偏倚性，與可口革囊星蟲的線粒體堿基特征相符［22］。

2.2 可口革囊星蟲群體cytb和d-loop序列結構及差異

我國東南沿海5個可口革囊星蟲地理群體150個樣本個體中：cytb基因序列的多態位點118個，其中轉換位點31個，顛換位點87個（顛換和轉換比2.81）（表2）；d-loop基因序列的多態位點153個，其中轉換位點66個，顛換位點87個（顛換和轉換比1.32）（表3）。在5個可口革囊星蟲地理群體中：廣東陽江群體的cytb序列的多態位點最多（31個），福建寧德群體的多態位點最少（16個）；廣東湛江群體的d-loop序列的多態位點最多（36個），福建寧德群體的多態位點最少（22個）。

表2 可口革囊星蟲群體的cytb序列差異Tab.2 Base variation of cytb sequence of P.esculenta populations

表3 可口革囊星蟲群體的d-loop序列差異Tab.3 Base variation of d-loop sequence of P.esculenta populations

2.3 可口革囊星蟲群體cytb和d-loop序列單倍型組成及分布

兩個基因序列定義的單倍型在各個地理群體間的分布如表4和表5所示。我國東南沿海5個可口革囊星蟲地理群體150個樣本個體通過cytb序列定義了62種單倍型，通過d-loop序列定義了79種單倍型。

表4 可口革囊星蟲群體cytb序列單倍型分布Tab.4 Haplotype distribution of cytb sequence of P.esculenta populations

續表4

通過表4中5個地理群體的cytb序列單倍型數量來看，陽江群體（22個）＞湛江群體（21個）＞欽州群體（20個）＞寧德群體（13個）＞溫州群體量（12個）。從單倍型的分布情況來看，Hap21在5個地理群體均有分布，也是在5個地理群體中分布最廣的單倍型，樣本數為23，占總樣本數15.33%，其余單倍體大部分分散分布在各個地理群體中。

通過表5中5個地理群體的d-loop序列單倍型數量來看，湛江群體（25個）＞陽江群體（23個）＞欽州群體（22個）＞溫州群體（20個）＞寧德群體（16個）。從單倍型的分布情況來看，Hap15和Hap67在5個地理群體均有分布，其中：Hap15單倍型的樣本數為11，占總樣本數7.33%；Hap67單倍型的樣本數為17，占總樣本數11.33%。其余單倍體大部分分散分布在各個地理群體中。

表5 可口革囊星蟲群體d-loop序列單倍型分布Tab.5 Haplotype distribution of d-loop sequence of P.esculenta populations

續表5

2.4 可口革囊星蟲群體cytb和d-loop序列遺傳多樣性

如表6所示，可口革囊星蟲cytb序列的基因多樣性指數（H）為0.871 3～0.963 2，核苷酸多樣性指數（π）為0.032 487～0.040 975。其中基因多樣性指數（H）和核苷酸多樣性指數（π）最高的均為欽州（QZ）群體（H＝0.963 2，π＝0.040 975）。而基因多樣性指數（H）最低的是溫州（WZ）群體（H＝0.871 3），核苷酸多樣性指數（π）最低的則是寧德（ND）群體（π＝0.032 487）。

表6 可口革囊星蟲cytb序列遺傳多樣性參數Tab.6 Genetic diversity parameter of cytb of P.esculenta

如表7所示，可口革囊星蟲d-loop序列的基因多樣性指數（H）為0.942 5～0.981 6，核苷酸多樣性指數（π）為0.033 247～0.046 151。其中基因多樣性指數（H）和核苷酸多樣性指數（π）最低的均是寧德（ND）群體（H＝0.942 5，π＝0.033 247），而基因多樣性指數（H）最高的是湛江（ZJ）群體（H＝0.981 6），核苷酸多樣性指數（π）最高的則是陽江（YJ）群體（π＝0.046 151）。

表7 可口革囊星蟲d-loop序列遺傳多樣性參數Tab.7 Genetic diversity parameter of d-loop of P.esculenta

2.5 可口革囊星蟲的群體遺傳結構

基于K2-P雙參數模型，構建單倍型鄰接樹（圖1，圖2）。各地理群體的單倍型均廣泛分布在單倍型鄰接樹上，大部分節點分支支持率較低（＜50%），單倍型聚類沒有展示出地域性特征。單倍型鄰接樹的拓撲結構簡單，未呈現明顯的地理譜系結構。

圖1 基于可口革囊星蟲cytb基因序列構建的單倍型鄰接樹Fig.1 Neighbor-joining tree of P.esculenta based on cytb gene sequence

圖2 基于可口革囊星蟲d-loop基因序列構建的單倍型鄰接樹Fig.2 Neighbor-joining tree of P.esculenta based on d-loop gene sequence

可口革囊星蟲5個群體兩兩之間的群體間遺傳距離如表8所示：基于cytb序列的遺傳距離計算結果顯示，遺傳距離最大的有兩對群體（YJ群體和QZ群體，WZ群體和QZ群體）均為0.006 2，最小的也有兩對群體（ZJ群體和ND群體，WZ群體和ND群體）均為0.005 4；基于dloop序列的遺傳距離計算結果顯示遺傳距離最大的是YJ群體和QZ群體為0.012 7，最小的是WZ群體和ND群體為0.010 3。

表8 不同群體可口革囊星蟲的群體間遺傳距離Tab.8 Genetic distance between population of P.esculenta in different populations

如表9所示，不同群體可口革囊星蟲cytb和d-loop序列的分子方差分析（AMOVA）結果一致表明，群體內部的遺傳變異（cytb：98.56%；d-loop：98.61%）比例顯著高于群體間的遺傳變異（cytb：1.44% ；d-loop：1.39%）。這表明了可口革囊星蟲群體的遺傳變異在不同地理群體間并不明顯，不存在明顯的地理種群分化。

表9 不同群體可口革囊星蟲cytb和d-loop序列的分子方差分析Tab.9 AMOVA analysis of cytb and d-loop in different populations of P.esculenta

2.6 可口革囊星蟲的群體歷史動態

從上文的分子方差分析的結果得出，5個地理群體間的遺傳分化不顯著。如表10所示，基于cytb和d-loop序列的Tajima’s D和Fu’s Fs檢驗結果都是負值，并且統計檢驗均有顯著性差異（P＜0.05）。綜合分析Tajima’s D和Fu’s Fs檢驗結果顯著且偏離中性這一現象，可推測出可口革囊星蟲群體近期經歷過群體快速擴張事件。

表10 可口革囊星蟲cytb和d-loop序列的中性檢驗Tab.10 Neutral test of cytb and d-loop in P.esculenta

3 討論

廣義定義中的遺傳多樣性包括了生物所攜帶遺傳信息的總和，不過一般提到的遺傳多樣性是指單一物種的種內或不同地理群體的遺傳變異總和。物種的遺傳多樣性水平密切影響著該物種適應環境的能力、生存能力和進化潛能等。如多樣性水平越高，物種就越能適應棲息地環境的變化。

本次研究通過基因多樣性指數（H）和核苷酸多樣性指數（π）兩個評價遺傳多樣性的常用參數，對東南沿海的廣西欽州（QZ）、廣東湛江（ZJ）、廣東陽江（YJ）、福建寧德（ND）、浙江溫州（WZ）5個地區的可口革囊星蟲進行遺傳多樣性分析。分析結果顯示，5個地理群體均表現出基因多樣性指數（H）處于較高水平、而核苷酸多樣性指數（π）處于較低水平的遺傳多樣性特征。該特征在金丹璐等［2］基于線粒體COI基因、王帥［23］基于線粒體COI和16SrRNA基因序列研究我國東南沿海其他地理群體的可口革囊星蟲的遺傳多樣性中也有體現：不同地理群體，兩個不同序列的分析結果都表現出較高水平的基因多樣性和較低水平的核苷酸多樣性的分布模式，其原因是基因多樣性的累積可在較短時間內提升，但核苷酸變異的累積需要的時間較長［24］。宋素霞等［14］、王劍平等［25］、陳康等［26］研究發現，cytb和d-loop基因序列在雙圍齒沙蠶（Perinereis aibuhitensis）［13］、光裸方格星蟲（Sipunculus nudus）［14］等與可口革囊星蟲同為海洋無脊椎類物種的遺傳多樣性中，都表現出這種高基因多樣性水平和低核苷酸多樣性水平的特征。該遺傳多樣性模式在這些海洋無脊柱動物的表現，反映了這些生物群體可能在較近的歷史時期經歷過快速擴張，即在較短的時間內由一個小規模的繁殖群體衍生為一個大群體［27］。

通過對5個地理群體共150個可口革囊星蟲樣品的分析，共檢測出cytb單倍型類型62種，dloop單倍型類型79種，其中在廣東陽江（YJ）和廣東湛江（ZJ）兩個地理群體中cytb和d-loop兩個不同基因的單倍型數量均比另外3個地理群體多，兩個不同基因序列的單倍型分布中只有1～2個單倍型在5個群體中均有分布。這表明可口革囊星蟲5個地理群體的兩種基因序列單倍型隨機分散在不同群體中，某些單倍型是多個群體共有，部分單倍型是個別群體特有。這樣的單倍型分布現象，從生物學角度可解釋為地理群體因一起近期的外來基因入侵事件及隨后的遺傳分化共同作用產生的結果。

群體遺傳結構主要受時間和環境變化的影響，是基因序列的變異在物種或種群中的一種非隨機分布，是生物遺傳資源應用及物種多樣性保護的基礎。不同群體可口革囊星蟲cytb和d-loop序列的分子方差分析結果以及它們的單倍型鄰接樹均一致，表明可口革囊星蟲群體的遺傳變異在不同地理群體間并不明顯，說明該物種的5個地理群體沒有顯著的遺傳分化，遺傳差異較小，也印證了近年來發表的基于COI基因或cytb基因分析其他地理群體可口革囊星蟲遺傳變異主要存在于群體間的結論［2，23，28］。導致該現象的原因可能是：1）可口革囊星蟲棲息在沿海的潮間帶，地理障礙少，而其早期生活史中的擔輪幼蟲和海球體幼蟲時期會隨著海水的運動發生遠距離擴散，導致不同地理群體之間有了基因交流［29］；2）因人工增養殖技術日漸成熟，養殖規模擴大，許多養殖戶從多個地方采購苗種投放到當地養殖場，而由于養殖場排水等導致的養殖種群逃逸現象時有發生，使各可口革囊星蟲地理群體間的基因交流頻繁［2］；3）因可口革囊星蟲廣泛分布在我國東南沿海地區的潮間帶，生境條件相似，間接影響了各個不同地理群體的遺傳變異方向。另外，通過綜合遺傳距離的兩種算法結果，福建寧德（ND）群體與其他群體兩兩之間的遺傳距離是最近的，可能與當地養殖場人工投放的苗種逃逸導致基因混淆有關［2］。因為采樣點牛姆嶼所在的內海灣是閩東地區最大的可口革囊星蟲養殖基地——霞浦縣長春鎮土筍養殖基地，因當地采挖的星蟲苗種滿足不了養殖需求，養殖戶從廣西、浙江等地大量采購苗種投放養殖，由于潮汐的影響人工苗種逃逸現象必然存在，這導致當地野生地理群體的基因混淆現象嚴重。

綜上分析，廣東湛江（ZJ）、廣東陽江（YJ）、廣西欽州（QZ）、福建寧德（ND）、浙江溫州（WZ）5個可口革囊星蟲地理群體兩個不同序列的分析結果都表現出較高水平的基因多樣性和較低水平的核苷酸多樣，總體遺傳多樣性水平較高；遺傳變異在不同地理群體間并不明顯，說明該物種的5個地理群體沒有顯著的遺傳分化，遺傳差異較小；另外，從兩個不同序列的Tajima’s D和Fu’s Fs檢驗結果可推測出，可口革囊星蟲群體近期經歷過群體快速擴張事件。本研究結果可以為可口革囊星蟲種質資源的保護與可持續開發利用提供參考數據。