中國北方常見芫菁分子物種界定(鞘翅目: 芫菁科)

程海云, 段家充, 張 超, 潘 昭

(河北大學生命科學學院, 生命科學與綠色發(fā)展研究院, 河北省動物系統(tǒng)學與應用重點實驗室, 河北保定 071002)

芫菁科(Meloidae)隸屬于鞘翅目(Coleoptera)擬步甲總科(Tenebrionoidea)，世界已知3亞科16族133屬近3 000種，廣泛分布于除新西蘭、南極和波利尼西亞群島之外的世界各地；中國記錄2亞科8族27屬近200種(亞種)，各省市自治區(qū)均有記載(Bolognaetal., 2010; Pan and Ren, 2020)。該科昆蟲是一類重要的資源昆蟲，因其復變態(tài)等獨特的生物學特性而被學者廣泛關注，尤其是斑蝥素的藥用價值更是近年芫菁應用研究的熱點之一(楊玉霞和任國棟, 2005; Bolognaetal., 2010)。然而，芫菁科部分物種的外部形態(tài)特征十分近似，在形態(tài)分類研究中經常出現(xiàn)難以區(qū)分、錯誤鑒定的情況，如北方常見物種圓點斑芫菁Mylabrisaulica和蘋斑芫菁Mylabriscalida(Pan and Ren, 2020)、橫紋溝芫菁Hycleussolonicus(Panetal., 2017)等。分子系統(tǒng)學的發(fā)展，特別是DNA條形碼的提出和興起，為解決以上問題提供了可行的方案。Liu等(2016)基于線粒體COI基因片段序列，對西北豆芫菁Epicautasibirica、疑豆芫菁Epicautadubia和中華豆芫菁Epicautachinensis分類研究的有效性進行了探討，提出后兩者是西北豆芫菁的次異名，然而該研究并未使用分子物種界定方法，僅是基于分支單系性得出的結論。

為更加全面地評估分子物種界定方法在芫菁科分類研究中的可行性，本研究以中國北方常見芫菁科昆蟲為研究材料，選取線粒體COI基因片段和核CAD基因片段序列，分別采用自動條形碼間隔探索(automatic barcode gap discovery, ABGD)、廣義混合Yule溯祖模型(generalized mixed Yule coalescent,GMYC)、貝葉斯泊松樹進程(Bayesian Poisson tree processes, bPTP)和貝葉斯系統(tǒng)發(fā)育和系統(tǒng)地理分析(Bayesian phylogenetics and phylogeography, BPP)方法進行分子物種界定。在對比形態(tài)鑒定結果的基礎上，探討以上序列和分析方法界定的準確性和適用性，以期為芫菁科的整合分類研究提供數(shù)據支持。

1 材料與方法

1.1 樣本選取

本研究以中國北方常見芫菁科昆蟲為研究對象，選取58個樣本，經形態(tài)鑒定，共涉及2亞科4族6屬18種；選擇擬步甲科(Tenebrionidae)黃粉蟲Tenebriomolitor為外群。實驗樣本大多來自河北大學博物館館藏，少部分樣本序列下載自GenBank，詳見表1。

1.2 形態(tài)鑒定

選取樣本均使用Nikon SMZ1500立式顯微鏡進行形態(tài)學鑒定。形態(tài)鑒定參考中國芫菁分類相關成果，如：Pan等(2010, 2017)，潘昭等(2010a, 2010b, 2011, 2013), Li等(2020)以及Pan和Bologna(2021)。

1.3 DNA的提取、PCR擴增及測序

從酒精泡制標本中取后足肌肉組織，使用昆蟲基因組DNA提取試劑盒(倍沃醫(yī)學科技)提取總DNA。通過PCR擴增COI和CAD基因片段，COI基因擴增引物為TW-N1284(5′-AAACTAACARCCTTC AAAGYTGT-3′)和COI-R2329(5′-ACHGTAAAYATA TGATGGGCTCA-3′)，反應程序：94℃預變性3 min；94℃變性1 min， 48℃延伸1.5 min， 40個循環(huán)； 72℃復性1.5 min；72℃最終延伸3 min。CAD基因擴增引物為CD439F(5′-TTCAGTGTACARTTYCAYC CHGARCAYAC-3′)和CD688R(5′-TGTATACCTAG AGGATCDACRTTYTCCATRTTRCA-3′)，反應程序： 94℃預變性3 min； 94℃變性30 s， 55℃延伸30 s， 50個循環(huán)； 72℃復性1 min； 72℃最終延伸3 min(Wild and Maddison, 2008)。PCR擴增采用25 μL體系：上下游引物(10 μmol/L)各1 μL，去離子水8.5 μL， PCR Mix 12.5 μL，DNA模板2 μL。 PCR產物經1%瓊脂糖凝膠電泳檢測，送通用生物系統(tǒng)(安徽)有限公司進行雙向測序。

1.4 基因序列處理和系統(tǒng)發(fā)育分析

所得序列使用DNASTAR SeqMan v7.1.0進行檢查和拼接，在NCBI中進行BLAST分析，確定所得序列是否為目的序列， 使用MEGA X (Kumaretal., 2018)中的Clustal W進行序列的比對。所有新得到的序列上傳至GenBank(登錄號見表1)。

使用MEGA X分別對COI和CAD序列數(shù)據集進行序列堿基含量分析及種內和種間遺傳距離分析(Kumaretal., 2018)。使用PhyloSuite v1.2.2 (Zhangetal., 2020)將COI和CAD序列數(shù)據集串聯(lián)，構建3個序列數(shù)據集，即COI,CAD和COI+CAD串聯(lián)序列數(shù)據集。使用PhyloSuite內嵌的MrBayes 3.2(Ronquistetal., 2012)，使用貝葉斯推斷法(Bayesian inference, BI)，分別基于3個數(shù)據集重建系統(tǒng)發(fā)育樹，運行2 000 000世代，采樣頻率為每100代一次，確保分裂頻率分支頻率的標準偏差小于0.01(Ronquist and Huelsenbeck, 2003)，以獲得高支持率的合意樹，生成合意樹之前舍棄25%的老化樣本，最佳堿基替代模型使用ModelFinder (Kalyaanamoorthyetal., 2017)進行評估。 最后，使用Evolview v3 (Subramanianetal., 2019)對系統(tǒng)樹進行可視化編輯。

1.5 物種界定分析

基于3個數(shù)據集，分別使用自動條形碼間隔探索(ABGD),廣義混合Yule溯祖模型(GMYC),貝葉斯泊松樹進程(bPTP)和貝葉斯系統(tǒng)發(fā)育和系統(tǒng)地理分析(BPP)方法進行物種界定分析，具體參數(shù)設置如下：

ABGD：使用ABGD在線分析(https:∥bioinfo.mnhn.fr/abi/public/abgd/abgdweb.html)(金倩等, 2017)，設置參數(shù)X(relative gap width)為0.5，其他參數(shù)為默認值，將得到的序列劃分情況與形態(tài)鑒定的結果進行比較分析。

GMYC(金倩等, 2017)：基于BEAST v1.10.4 (Drummond and Rambaut, 2007)構建超度量二歧分支樹(有根)，選擇嚴格分子鐘。MCMC 鏈為10 000萬代，每1 000代參數(shù)取樣一次，Burn-In默認為鏈長10%，使用Tracer v1.7.2(Rambautetal., 2018)評估系統(tǒng)樹的收斂情況。將含有樹文件的*.trees.txt 文件導入到Tree Annotator v1.10.4生成所需系統(tǒng)樹。Burn In (as trees)根據 Tracer v1.7.2分析結果(Burn-In 鏈長/采樣頻率)設置為1 000，最后生成所需最大分支可信度樹(maximum clade credibility tree)，導入GMYC進行分析。

bPTP(Zhangetal., 2013)：本實驗使用PTP的更新版本bPTP進行界定分析，將GMYC分析中得到的Newick系統(tǒng)樹上傳到PTP在線網站(https:∥species.h-its.org/ptp/)進行分析，設置為有根樹，輸入外群名稱，馬爾科夫鏈代數(shù)為500 000，其他參數(shù)為默認值。

BPP(Yang, 2015)：基于串聯(lián)基因數(shù)據集構建的貝葉斯系統(tǒng)樹，在BPP v4.2.9中利用A11模板輸入指導樹(基于COI序列構建的BI樹)，運行分析。

2 結果

2.1 基因序列特征

COI基因片段序列經拼接、比對、剪切處理后的長度為561 bp。使用MEGA X軟件中的Kimura雙參數(shù)模型(Kimura 2-parameter model)進行序列堿基含量分析及種內、種間遺傳距離分析(表2和3)，結果顯示：序列的T, C, A和G平均含量分別為33.5%, 22.3%, 26.8%和17.5%，平均A+T含量顯著高于G+C含量，變異多發(fā)生在堿基第3位；Hycleussolonicus,Mylabrisaulica,Epicautamegalocephala,Mylabrissibirica,Lyttacaraganae,Meloeproscarabaeus,Meloeviolaceus和Epicautasibirica的種內平均遺傳距離分別為0.92%, 0.71%, 0.83%, 0.30%, 1.66%, 1.34%, 0.11%和0.24%，其余10個物種樣本量不足，無法進行種內遺傳距離分析；18種芫菁的種間遺傳距離在9.49%～22.10%之間，平均值為18.36%。

表2 基于COI基因片段序列計算得出芫菁科種內平均遺傳距離和種內遺傳距離范圍Table 2 Average intraspecific genetic distance and range of intraspecific genetic distance of Meloidaecalculated based on COI fragment sequences

CAD基因片段序列處理后長度為730 bp。使用MEGA X 軟件中的Kimura 雙參數(shù)模型進行序列堿基含量分析及種內種間遺傳距離分析(表4和5)，結果顯示：T, C, A和G平均含量分別為33.5%, 22.3%, 26.8%和17.5%，平均A+T含量顯著高于G+C含量，沒有堿基變異偏向性，密碼子的1, 2和3位均有發(fā)生變異；Hycleussolonicus,Mylabrisaulica,Epicautamegalocephala,Epicautasibirica,Lyttacaraganae,Meloeproscarabaeus和Meloeviolaceus的種內平均遺傳距離為1.71%, 0.65%, 0.91%, 0.73%, 0.46%, 0.27%和0.06%，Mylabrissibirica沒有變異，其余10個物種樣品量不足，無法進行種內遺傳距離分析; 18種芫菁的種間遺傳距離是處于0.43%～22.33%之間，平均值為13.85%。

表4 基于CAD基因片段序列計算得出芫菁科種內平均遺傳距離和種內遺傳距離范圍Table 4 Average intraspecific genetic distance and range of intraspecific genetic distance ofMeloidae calculated based on CAD fragment sequences

2.2 系統(tǒng)發(fā)育

利用PhyloSuite內嵌的ModelFinder(Kalyaanamoorthyetal., 2017)進行最佳堿基替代模型分析， 實現(xiàn)Akaike’sInformation Criterion (AIC)選擇最適模型，3個數(shù)據集的最適模型均為GTR+F+I+G4。

對58頭芫菁樣本進行BI分析，基于COI,CAD和COI+CAD串聯(lián)序列數(shù)據集分別建樹，顯示出高度相似的拓撲結構(圖1)，聯(lián)合數(shù)據集分支支持度高于單基因建樹，高度支持這18個種的單系性。除Hycleussolonicus,Meloevariegatus,Meloebrevicollis和Epicautaobscurocephala這幾個分支的后驗概率(posterior probability, PP)<1，其余種劃分PP均等于1。

2.3 物種界定

本實驗對58個芫菁樣本(不含外群)的分子物種界定結果如圖1所示，具體如下：

ABGD劃分結果包括初始劃分(initial partition)和遞歸劃分(recursive partition)兩種情況。 基于COI基因片段和COI+CAD串聯(lián)基因片段序列的劃分結果顯示，初始劃分較為穩(wěn)定，遺傳距離在0.001～0.1之間，被劃分為18個分子操作分類單元(molecular operational taxonomic unit, MOTU)，支持18個形態(tài)種，且與BI樹拓撲結構一致。遞歸劃分先驗值浮動較大，不夠穩(wěn)定?；贑AD基因片段序列的ABGD劃分結果顯示，屬級分類單元的劃分和形態(tài)劃分吻合，而物種階元上則將Meloeproscarabaeus和Meloeviolaceus劃分為一個MOTU，物種劃分紊亂，如圖2(A, B, C)。

基于3個數(shù)據集的GMYC分析所劃分的MOTUs結果皆與形態(tài)種相符。單閾值分析，基于COI基因片段序列劃分的似然實體數(shù)量為19，置信區(qū)間為19～22，基于CAD基因片段序列劃分的似然實體數(shù)量為19，置信區(qū)間為1～49，基于COI+CAD串聯(lián)序列劃分的似然實體數(shù)量為19，置信區(qū)間為15～21，如圖3-5所示。

bPTP分析結果表明，基于COI基因序列劃分得到21個MOTUs，將Mylabrisaulica,Meloeproscarabaeus和Lyttacaraganae各自分為了2個MOTUs；基于CAD基因序列劃分為17個MOTUs，將Meloeproscarabaeus和Meloeviolaceus歸為一個MOTU；基于COI+CAD基因串聯(lián)序列得到18個MOTUs，和形態(tài)學劃分結果一致，如表6-8所示。

表6 基于COI序列通過bPTP方法得到的芫菁科物種劃分Table 6 Species division of Meloidae based on COI sequence by bPTP method

表7 基于CAD序列通過bPTP方法得到的芫菁科物種劃分Table 7 Species division of Meloidae based on CAD sequence by bPTP method

表8 基于COI+CAD串聯(lián)序列通過bPTP方法得到的芫菁科物種劃分Table 8 Species division of Meloidae based on concatenated COI+CAD sequence by bPTP method

BPP是基于給定指導樹進行分析，得到18個MOTUs的后驗概率(PP)是0.98925，遠大于其他劃分方式的PP，即支持將58個樣本劃分為18個獨立物種，支持指導樹給定的18個物種，與形態(tài)學劃分一致。

3 討論和結論

Liu等 (2016)首次應用分子系統(tǒng)學方法對芫菁科昆蟲的物種劃分問題進行了探討并提出了2個新異名，然而此工作依據的是外部形態(tài)、遺傳距離和系統(tǒng)樹拓撲結構的綜合結果，并未使用其他分子物種界定方法進行分析。

本研究對中國北方常見的芫菁科6屬18個形態(tài)種進行分子物種界定，總體而言所得結果基本與形態(tài)劃分一致，僅基于CAD序列的ABGD分析結果和基于單基因的bPTP分析結果有所差異(圖1)。該結果一方面驗證了形態(tài)種劃分結果，例如目前使用觸角末節(jié)形狀對具有相似鞘翅斑紋的圓點斑芫菁和蘋斑芫菁進行區(qū)分(Pan and Ren, 2020)，該形態(tài)界定結果得到了分子物種界定結果的支持；另一方面也顯現(xiàn)出對芫菁科昆蟲相對可靠的可用于分子物種界定的DNA序列和分析方法，即多基因片段串聯(lián)序列優(yōu)于COI序列優(yōu)于CAD序列，ABGD方法最簡便、準確性相對較高且穩(wěn)定，在一定程度上可減少人為界定的誤差。在進行ABGD和bPTP分析時，發(fā)現(xiàn)基于CAD的劃分對物種劃分不清晰，因此推測是由于CAD基因進化速率較慢，種間差異小，所以導致種的劃分紊亂。

此外，參考分子物種界定結果，基于COI基因種內和種間的遺傳距離分析結果表明，在芫菁科中，當遺傳距離在1.66%以下可認為是同一物種的種內變異，而大于9.49%則可認為是不同的物種，而在1.66%～9.49%之間則需要綜合更多證據謹慎處理(表2和3)?；贑AD基因種內和種間的遺傳距離和物種界定綜合分析結果表明，在芫菁科中，當遺傳距離在1.71%以下可認為是同一物種的種內變異，而大于1.71%則有成為不同物種的可能，在1.71%～2.77%之間則需要綜合更多證據謹慎處理(表4和5)。

綜上，基于COI和CAD的獨立和串聯(lián)數(shù)據集和ABGD, GMYC, PTP和BPP 4種分析方法均適用于芫菁科昆蟲的分子物種界定，但有一定的誤差，需多基因、多方法綜合考慮。此外，形態(tài)和分子物種界定方法的結合能夠有效提高芫菁科昆蟲物種劃分的準確性。