基于線粒體Cytb和D-Loop序列的大獺蛤5個地理群體遺傳多樣性研究

陳 康，王偉峰，陳秀荔，朱威霖，王煥嶺

(1.華中農業大學水產學院，湖北武漢 430070；2.廣西壯族自治區水產科學研究院，廣西水產遺傳育種與健康養殖重點實驗室，廣西南寧 530021)

大獺蛤(Lutrariamaxima)隸屬雙殼綱(Bivalvia)，異齒亞綱(Heterodonta)，真瓣鰓目(Eulamellibranchia)，蛤蜊科(Mactridae)，獺蛤屬，為營埋棲生活的二倍體貝類，廣泛分布于中國中南部沿海，北至福建省，南至海南省，尤其在北部灣地區較為豐富，目前通過人工養殖，已發展成為廣西、廣東等地淺海水產養殖的主要品種之一[1-3]。近年來，我國沿海地區建設了眾多高能耗、高排放的大型項目，大量生活污水和工業廢水排入海中，導致海水污染加重，赤潮暴發次數增加，水生生物棲息環境被嚴重破壞[4-5]。據相關報道，在北部灣貝類市場檢測的所有貝類生物中均含有麻痹性貝類毒素[6]，以貝類為主的經濟物種面臨巨大的生存挑戰和繁殖壓力。另一方面，人工育種技術不斷成熟[7]，育種范圍相對集中，導致大獺蛤近親繁育加劇，很可能造成種質資源不斷退化。因此，對大獺蛤野生群體的資源評估顯得十分必要。

為了促進大獺蛤資源的可持續利用，本研究采用分析遺傳多樣性的方式為其遺傳育種和種質資源的評估提供重要的依據[8]。由于多態性信息豐富、遺傳穩定性強和檢測手段相對簡便等特點，DNA分子標記在群體遺傳結構及其分化研究方面具有無可比擬的優勢[9]。其中線粒體DNA(mitochondrial DNA，mtDNA)是閉合、環狀的雙鏈DNA分子，其作為遺傳信息的重要載體，具有母系遺傳、結構簡單、進化速度快、幾乎無重組及不同的區域進化速度存在差異等特點，使得mtDNA成為一種重要的分子標記[10]。在mtDNA中，細胞色素b(Cytochrome b，Cytb)和mtDNA控制區(又稱D-環區，control region displacement loop，D-Loop) 的進化速度存在差異且變異速率快[11-12]。因此，利用Cytb和D-loop序列分析和對比不同群體mtDNA多態性，能夠揭示種內或種間的親緣關系[13]。近年來，通過分析Cytb和D-Loop 序列的堿基信息檢測水生生物遺傳多樣性已被廣泛應用，如基于Cytb基因序列發現中國養殖區和日本原產地的蝦夷扇貝(Patinopectenyessoensis)出現明顯的遺傳分化[14]；王劍平等[15]基于Cytb基因發現洞庭湖河蜆(Corbiculafluminea)遺傳分化弱，變異主要來自群體內；GUO等[16]利用Cytb和D-Loop序列得出西藏雅魯藏布江6個地理群體的裸裂尻魚(Schizopygopsisyounghusbandi)的遺傳多樣性水平較低，并且遺傳變異來自群體內部，這些研究為水生生物種質資源的開發與利用奠定了理論依據。

目前，關于大獺蛤的研究主要集中在營養成分分析與評價[17]、繁殖特性[18]及人工育苗[7,19]等方面。而關于大獺蛤種質資源評估的研究鮮有報道，相關研究僅見李斌等[20]采用形態學特征和RAPD技術，對3個地理群體的施氏獺蛤(L.siebaldii)與越南大獺蛤(L.maxima)遺傳差異進行了比較分析。因此為了評估我國不同地理群體大獺蛤資源，本研究利用Cytb和D-Loop序列，對我國大獺蛤5個地理群體的遺傳多樣性、群體結構和群體歷史動態進行分析，以期揭示大獺蛤種質資源的遺傳背景，為大獺蛤的人工繁育和資源保護提供理論依據。

1 材料與方法

1.1 實驗材料與DNA提取

大獺蛤分別采自廣西北海(BH)、廣西潿洲島(WZD)、廣東湛江(ZJ)、福建廈門(XM)和福建福州(FZ)5個地區(圖1)，取前后閉殼肌組織保存于95%乙醇中，用于基因組DNA的提取，每個群體的個體數及遺傳多樣性參數見表1。采用醋酸銨法[21]從肌肉組織中提取大獺蛤基因組DNA，1%瓊脂糖凝膠電泳和Nanodrop測定DNA的質量和濃度。

1.2 PCR擴增

根據NCBI數據庫中大獺蛤mtDNA全序列(Genbank：NC_036766.1)利用Primer Premier 5.0設計引物用于擴增Cytb和D-Loop序列。引物序列分別為，CytbF：5′ -TGCGGCTGTCTGGTATTGA-3′、CytbR：5′-AACCCTTTCATCGTCCCACTA -3′和D-Loop F：5′-ATTAGAATACGCCGTTGAAG-3′、D-Loop R：5′-GAGTAGTTACATCCTGCTTAC-3′，由武漢天一輝遠生物科技有限公司合成。

PCR總反應體積為10 μL，反應體系為：10×PCR buffer 1 μL，dNTPs 0.4 μL，雙蒸水6.6 μL，DNA模板1 μL，ESTaqDNA聚合酶0.2 μL，左右引物各0.4 μL。擴增程序為：預變性 94℃ 3 min，變性94℃ 30 s，退火30 s(D-Loop為52℃，Cytb為54℃)，延伸72 ℃(D-Loop為45 s，Cytb為105 s)，共32個循環，終延伸72℃ 7 min。擴增產物經1%瓊脂糖凝膠電泳檢測后，選擇目的條帶明顯且特異性好的PCR產物送武漢天一輝遠生物科技有限公司測序。

1.3 數據分析

基于測序峰圖使用DNAstar軟件包[22]中的Seqman軟件對測序結果進行分析校對，以確認序列質量，并與NCBI數據庫中公布的大獺蛤mtDNA進行對比，獲得不同大獺蛤群體的Cytb和D-Loop序列。MEGA7.0[23]計算Cytb和D-Loop序列的堿基組成并使用鄰接法(neighbor-joining，NJ)構建系統發育樹，采用Bootstrap(重復次數1 000)檢驗聚類樹各分支置信度。DNAsp5[24]獲得Cytb和D-Loop序列的核苷酸多樣性指數(nucleotide diversity，Pi)、平均核苷酸差異數(average number of nucleotide differences，K)、單倍型數(number of haplotypes，H)、單倍型多樣性指數(haplotype diversity，Hd)及多態位點數(number of polymorphic sites，S)。運用Arlequin3[25]計算群體間的遺傳分化系數和擴張參數τ，并結合中性檢驗和核苷酸不配對分布曲線來推測大獺蛤種群的歷史動態情況；根據公式τ= 2ut和T=t×(代時)，估算種群擴張時間T[26]。其中u=μk，μ為研究序列變異速率，k為研究的序列長度(bp)；大獺蛤性成熟約為1年，則代時為1[18]。本研究中，Cytb和D-loop的變異速度分別為(1.0%～2.5%)/百萬年[27]、(3%～10%)/百萬年[28]。

圖1 大獺蛤采樣分布圖 Fig.1 Sampling sites for L. maxima注：BH、WZD、ZJ、XM、FZ分別代表廣西北海、廣西潿洲島、廣東湛江、福建廈門、福建福州；●代表采樣點的位置Note：BH, WZD, ZJ, XM， FZ represent Beihai of Guangxi, Weizhou Island of Guangxi, Zhanjiang of Guangdong, Xiamen and Fuzhou of Fujian respectively；●indicates sampling site

表1 大獺蛤5個群體的遺傳多樣性參數Tab.1 Genetic diversity parameters of five populations of L. maxima

注：S：多態位點數；H：單倍型數；Hd：單倍型多樣性指數；Pi：核苷酸多樣性指數；K：平均核苷酸差異數

Note:S: number of polymorphic sites;H: number of haplotypes;Hd: haplotype diversity;Pi: nucleotide diversity;K: average number of nucleotide differences

2 結果與分析

2.1 大獺蛤線粒體Cytb和D-Loop的序列特征和遺傳多樣性

測序結果經軟件處理后分別獲得了大獺蛤Cytb基因序列(1 155 bp)和D-Loop序列(412 bp)。利用MEGA7軟件計算序列的堿基組成，結果顯示5個群體間Cytb和D-Loop序列的堿基含量基本一致，Cytb基因和D-Loop序列4種堿基T、C、A、G的平均含量分別占其全長的43.60%、13.53%、20.61%、22.26%和23.31%、25.25%、39.29%、12.16%，且兩序列的A+T含量明顯高于G+C含量(表2)，表現出了AT偏移，這與大多數水生動物的已知模式一致。使用DNAsp5對大獺蛤5個地理群體的遺傳參數進行計算，結果顯示Cytb基因在5個群體的161個個體的核苷酸多樣性指數為0.002 3±0.001 4，單倍型多樣性為0.834 5±0.030 3，平均核苷酸差異數為2.631 3，共檢索到89個多態位點和73種單倍型。D-Loop序列在5個群體中159個個體的核苷酸多樣性指數為0.001 4±0.001 3，單倍型多樣性為0.445 6±0.049 8，平均核苷酸差異為0.562 7，共發現25個多態位點和27種單倍型。其中，北海群體的單倍型多樣性指數Cytb序列為0.903 2，D-Loop序列為0.584 9，該群體的單倍型也是比較豐富的(表1)。

2.2 大獺蛤種群的遺傳分化分析

在種群遺傳學中，遺傳分化指數(Fst)可反映群體間的遺傳分化程度，有學者將遺傳分化程度分為弱分化(Fst<0.05)、中度分化(0.05≤Fst≤0.25)、高度分化(Fst>0.25)3個等級[29]。本研究基于Cytb和D-Loop序列分析，北海、福州、潿洲島、廈門、湛江5個群體Fst分別在-0.008 0～0.010 1(P>0.05)和-0.014 3～-0.007 2(P>0.05)之間變動，同時5個地理群體間的遺傳距離分別在0.001 8～0.002 9和0.000 8～0.002 0之間(表3,表4)。結果表明，每個群體之間遺傳分化程度都較弱(Fst<0.05)，遺傳距離較近，說明5個群體之間具有較高的遺傳同質性。

基于Cytb對大獺蛤群體AMOVA分析表明：99.68%的差異屬于群體內差異，0.32%為群體間差異(表5)；基于D-Loop對大獺蛤群體AMOVA分析表明：101.16%的差異屬于群體內差異，-1.16%為群體間差異(表5)，表明群體內遺傳變異遠遠高于群體間的變異。

2.3 系統進化樹

基于Cytb基因(圖2-A)和D-Loop序列(圖2-B)的兩個系統發育樹都分為兩大支系，未形成明顯的單系群，5個地理群體間相互交叉，無明顯的地理群體聚類，與上述群體間遺傳距離弱、群體間分化指數Fst值低的結果相符。

表2 大獺蛤5個群體Cytb和D-Loop序列的堿基組成Tab.2 Base composition of Cytb and D-Loop sequences in five populations of L. maxima

注：BH、WZD、ZJ、XM、FZ分別代表廣西北海、廣西潿洲島、廣東湛江、福建廈門、福建福州

Note: BH, WZD, ZJ, XM， FZ represent Beihai of Guangxi, Weizhou Island of Guangxi, Zhanjiang of Guangdong, Xiamen and Fuzhou of Fujian respectively

表3 基于Cytb基因的大獺蛤群體間遺傳距離及遺傳分化指數FstTab.3 Genetic distance and fixation index (Fst) between every two populations of L. maxima based on Cytb gene

注：對角線下數據為群體間遺傳距離，對角線上數據為遺傳分化指數Fst；BH、WZD、ZJ、XM、FZ分別代表廣西北海、廣西潿洲島、廣東湛江、福建廈門、福建福州；Fst值均無顯著性(P>0.05)

Note：Data below the diagonal mean genetic distance, above the diagonal mean fixation indexFst; BH, WZD, ZJ, XM， FZ represent Beihai of Guangxi, Weizhou Island of Guangxi, Zhanjiang of Guangdong, Xiamen and Fuzhou of Fujian respectively; allFstvalues were not significant (P>0.05)

表4 基于D-Loop序列的大獺蛤群體間遺傳距離及遺傳分化指數FstTab.4 Genetic distance and fixation index (Fst) between every two populations of L. maxima based on D-Loop sequence

注：對角線下數據為群體間遺傳距離，對角線上數據為遺傳分化指數Fst；BH、WZD、ZJ、XM、FZ分別代表廣西北海、廣西潿洲島、廣東湛江、福建廈門、福建福州；Fst值均無顯著性(P>0.05)

Note：Data below the diagonal mean genetic distance, above the diagonal mean fixation indexFst; BH, WZD, ZJ, XM, FZ represent Beihai of Guangxi, Weizhou Island of Guangxi, Zhanjiang of Guangdong, Xiamen and Fuzhou of Fujian respectively; allFstvalues were not significant (P>0.05)

表5 大獺蛤5個地理群體Cytb和D-Loop序列的AMOVA結果Tab.5 AMOVA result of five populations based on Cytb and D-Loop sequences

注： Va、Vb 分別表示群體間變異和群體內變異

Note: Va, Vb mean variations among populations and within populations

圖2 基于Cytb(A)基因和D-Loop(B)序列構建大獺蛤單倍型NJ樹Fig.2 The NJ phylogenetic tree of L. maxima based on Cytb (A) and D-Loop(B) haplotypes注：每個群體的私有單倍型都用圓圈標記，而群體之間的其他共享單倍型則不標記。BH、WZD、ZJ、XM、FZ分別代表廣西北海、廣西潿洲島、廣東湛江、福建廈門、福建福州Note：The private haplotypes of each population are marked with circles, and the other shared haplotypes among populations are not marked. BH, WZD, ZJ, XM, FZ represent Beihai of Guangxi, Weizhou Island of Guangxi, Zhanjiang of Guangdong, Xiamen and Fuzhou of Fujian respectively

2.4 種群動態

由于大獺蛤5個地理群體間的遺傳分化程度較弱，故分析其歷史動態信息時將其看作一個整體，并依據Tajima’sD[30]和 Fu’sFS[31]值中性檢驗和核苷酸不配對分布曲線來判定大獺蛤在過去是否發生了種群歷史擴張。大獺蛤Cytb和D-Loop序列的中性檢驗分析結果如表6所示。基于Cytb序列各群體的Tajima’D和Fu’sFS，均為負值，統計學分析都存在顯著差異(P<0.05)，將5個群體作為一個整體進行中性檢驗，其Tajima’sD和 Fu’sFS值分別為-2.229 0(P<0.01)和-12.658 8(P<0.01)，Cytb序列核苷酸錯配分布圖為單峰(圖3-A)。基于D-Loop序列各群體的Tajima’D和Fu’sFS，均為負值，統計學分析都存在顯著差異(P<0.05)，將5個群體作為一個整體進行中性檢驗，其Tajima’sD和 Fu’sFS值分別為-1.943 2(P<0.05)和-6.248 6(P<0.05)，D-Loop序列核苷酸錯配分布圖為單峰(圖3-B)。這些結果表明大獺蛤在歷史上發生過明顯的種群擴張。在擬合度檢驗中(表6)，基于Cytb和D-Loop序列分析，5個群體的偏離方差(sum of squared deviation)和Raggedness統計量(R)的值都較小，且大部分不存在顯著差異(P>0.05)。從整體來看，這兩個參數的值都較小，均未達到顯著水平(P>0.05)，同樣表明大獺蛤歷史上經歷了種群擴張。基于Cytb和D-Loop序列的平均擴張參數τ值分別為1.523 8、1.038 7，推測出種群擴張時間距今約為1.3×104～6.6×104年，表明大獺蛤種群擴張時間約在更新世晚期。

表6 大獺蛤5個群體的中性檢驗、擬合度檢驗Tab.6 Neutrality test and test of goodness of fit in five populations of L. maxima

注：BH、WZD、ZJ、XM、FZ分別代表廣西北海、廣西潿洲島、廣東湛江、福建廈門、福建福州

Note：BH, WZD, ZJ, XM, FZ represent Beihai of Guangxi, Weizhou Island of Guangxi, Zhanjiang of Guangdong, Xiamen and Fuzhou of Fujian respectively

圖3 基于Cytb (A)和D-Loop (B)序列的大獺蛤5個群體核苷酸不配對分布圖Fig.3 Mismatch distribution of five populations based on Cytb (A) and D-Loop (B) sequences

3 討論

遺傳多樣性又稱基因多樣性，生物種群的遺傳多樣性是其物種進化和環境適應的基礎，遺傳多樣性缺乏可能會導致物種的資源衰退甚至瀕臨滅絕，豐富的遺傳多樣性能夠確保物種的延續，同時也為物種進化提供充足的潛力。

3.1 遺傳多樣性

大獺蛤5個群體的線粒體Cytb基因序列中，檢測到89個多態位點和73個單倍型，單倍型多樣性為0.834 5±0.030 3；D-Loop基因序列中，檢測到25個多態位點和27個單倍型，單倍型多樣性為0.445 6±0.049 8，整體呈現出中等偏高的遺傳多樣性，但Cytb和D-Loop的核苷酸多樣性僅分別為0.002 3±0.001 4和0.001 4±0.001 3，處于一個較低水平。一般來說，當群體數量足夠大并且在很長一段時間內保持穩定，該群體中的高水平遺傳變異是可以維持的，而當群體經歷嚴重的瓶頸效應時，隨機遺傳漂變將導致遺傳變異的大量喪失[11]。根據GRANT和BOWEN[32]提出的分類標準，以Hd=0.5、Pi=0.005為界，本研究結果與高單倍型多樣性和低核苷酸多樣性的模式相符。該遺傳多樣性模式的種群在歷史上有經歷擴張的可能，即某種群在數量急劇減少后，存留的優良個體歷經瓶頸效應或建群效應，突然快速繁衍發展為一個大種群的現象，在擴張過程中，種群數量的增加導致了單倍型多樣性的提高，而無充足的時間去積累核苷酸變異，所以會產生單倍型多樣性較高而核苷酸多樣性較低的遺傳多樣性模式[33]。

基于Cytb基因獲得的5個大獺蛤地理群體的總體遺傳多樣性指數大于D-Loop序列的結果，通常情況下D-Loop序列獲得的種群遺傳變異大于Cytb基因[10]。目前有人發現mtDNA控制區D-Loop變異速率低于Cytb的現象，趙亮等[34]對太湖新銀魚(Neosalanxtaihuensis)Cytb和D-Loop序列采用貝葉斯(Bayes)MCMC模擬估測出Cytb相對速率為1.000±0.131，而D-Loop相對速率為0.859±0.261，Cytb的進化速率相對比D-Loop高。

3.2 種群遺傳結構

遺傳距離是反映不同物種或不同種群之間變異的常用參數。5個群體之間的遺傳距離幾乎處于同一水平，說明各群體之間的基因交流不受影響。ARNAUD等[35]認為在幼蟲期具有浮游特性的海洋雙殼類生物具有很強的分散能力，大獺蛤具有大約12 d的浮游期[18]，在此期間一定程度上受到我國東南海暖流和東南部沿岸流的影響[36]，從而被動地隨著海流擴散和基因交流；此外我國大獺蛤苗種培育[19]的成功和工廠化育苗技術[7]的日漸成熟，導致各地大獺蛤群體之間交流加劇。

遺傳分化指數Fst是作為對種群遺傳分化評估的另一個重要指標。基于Cytb的Fst值在-0.008 0～0.010 1之間，基于D-Loop的Fst值在-0.014 3 ～ -0.007 2之間，Fst值結果也顯示5個群體間存在負值現象，表明群體內差異大于群體間差異，同時根據WRIGHT[29]定義的標準和AMOVA分析結果表明，群體內的變異遠大于群體間的變異，說明各地理群體間不存在顯著的遺傳分化。基于Cytb基因和D-Loop序列構建的NJ樹中，5個群體的單倍型相互混雜，未形成特定的地理聚群。這些結果表明各區域個體間的交配是隨機的，沒有顯著的遺傳分化。黃榮蓮等[37]在對貽貝(Pernaviridis)的研究中發現環境脅迫會對群體結構產生影響。本研究采樣點福州(26°08′ N、119°30′ E)、廈門(24°48′ N、118°09′ E)和北海(21°48′ N、109°12′ E)為亞熱帶季風性氣候；湛江(21°27′ N、110°37′ E)和潿洲島(21°04′ N、109°12′ E)為熱帶季風性氣候，5個采樣地點的氣候、溫度和鹽度等環境因素都大致相同，這可能是遺傳分化弱的原因之一。

3.3 種群歷史動態

通過分析種群的核苷酸不配對分布曲線的單峰或多峰類型，以及是否符合中性檢驗，來推斷大獺蛤種群近期是否出現過種群擴張事件[38]。當種群偏離中性檢驗，Tajima’D檢驗會給出一個較低值。Fu’sFS值的負數越大表明由于種群擴張或種群選擇導致中性偏離。整體來說，基于大獺蛤Cytb和D-Loop的Tajima’D和Fu’sFS值均為負值，且統計檢驗顯著(P<0.05)，核苷酸錯配分布曲線為單峰，說明大獺蛤種群歷史上可能經歷過瓶頸和種群擴張，上述的高單倍型多樣性和低核苷酸多樣性的分布模式也證明了這一點。根據平均τ值估算的大獺蛤種群擴張時間距今約為1.3×104～ 6.6×104年，這個時期屬于更新世晚期，中國沿海海洋魚類的種群擴張大多發生在該時期[39]。該時期由于冰期和間冰期的交替，導致海平面劇烈變化，使得海平面下降120 ～ 140 m，現生海洋物種在近期歷史上的擴張大多受到這種冰期和間冰期交替的影響[40-41]。因此，大獺蛤種群可能經歷了同樣的種群擴張。

綜上所述，本研究基于線粒體Cytb和D-Loop序列對大獺蛤5個地理群體的遺傳多樣性進行了分析。整體來說，各群體間遺傳距離較小，遺傳分化程度較弱，遺傳變異水平較低，可作為一個管理保護單位，且近期發生過種群擴張事件；結合群體內遺傳距離和平均遺傳多樣性指數，可以得出大獺蛤的遺傳多樣性處于一個中等的水平。但近年來東南沿海地區生態環境破壞的加劇和人工繁殖技術的日臻成熟，勢必會對大獺蛤的種質資源造成潛在的威脅。因此在合理的開發利用下，有必要對野生資源進行保護，以保證其行業的可持續性發展。