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

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

(河北大學(xué)生命科學(xué)學(xué)院, 生命科學(xué)與綠色發(fā)展研究院, 河北省動物系統(tǒng)學(xué)與應(yīng)用重點(diǎn)實(shí)驗(yàn)室, 河北保定 071002)

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

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

1 材料與方法

1.1 樣本選取

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

1.2 形態(tài)鑒定

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

1.3 DNA的提取、PCR擴(kuò)增及測序

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

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

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

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

1.5 物種界定分析

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

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

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

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

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

2 結(jié)果

2.1 基因序列特征

COI基因片段序列經(jīng)拼接、比對、剪切處理后的長度為561 bp。使用MEGA X軟件中的Kimura雙參數(shù)模型(Kimura 2-parameter model)進(jìn)行序列堿基含量分析及種內(nèi)、種間遺傳距離分析(表2和3)，結(jié)果顯示：序列的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的種內(nèi)平均遺傳距離分別為0.92%, 0.71%, 0.83%, 0.30%, 1.66%, 1.34%, 0.11%和0.24%，其余10個物種樣本量不足，無法進(jìn)行種內(nèi)遺傳距離分析；18種芫菁的種間遺傳距離在9.49%～22.10%之間，平均值為18.36%。

表2 基于COI基因片段序列計算得出芫菁科種內(nèi)平均遺傳距離和種內(nèi)遺傳距離范圍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ù)模型進(jìn)行序列堿基含量分析及種內(nèi)種間遺傳距離分析(表4和5)，結(jié)果顯示：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的種內(nèi)平均遺傳距離為1.71%, 0.65%, 0.91%, 0.73%, 0.46%, 0.27%和0.06%，Mylabrissibirica沒有變異，其余10個物種樣品量不足，無法進(jìn)行種內(nèi)遺傳距離分析; 18種芫菁的種間遺傳距離是處于0.43%～22.33%之間，平均值為13.85%。

表4 基于CAD基因片段序列計算得出芫菁科種內(nèi)平均遺傳距離和種內(nèi)遺傳距離范圍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內(nèi)嵌的ModelFinder(Kalyaanamoorthyetal., 2017)進(jìn)行最佳堿基替代模型分析， 實(shí)現(xiàn)Akaike’sInformation Criterion (AIC)選擇最適模型，3個數(shù)據(jù)集的最適模型均為GTR+F+I+G4。

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

2.3 物種界定

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

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

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

bPTP分析結(jié)果表明，基于COI基因序列劃分得到21個MOTUs，將Mylabrisaulica,Meloeproscarabaeus和Lyttacaraganae各自分為了2個MOTUs；基于CAD基因序列劃分為17個MOTUs，將Meloeproscarabaeus和Meloeviolaceus歸為一個MOTU；基于COI+CAD基因串聯(lián)序列得到18個MOTUs，和形態(tài)學(xué)劃分結(jié)果一致，如表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是基于給定指導(dǎo)樹進(jìn)行分析，得到18個MOTUs的后驗(yàn)概率(PP)是0.98925，遠(yuǎn)大于其他劃分方式的PP，即支持將58個樣本劃分為18個獨(dú)立物種，支持指導(dǎo)樹給定的18個物種，與形態(tài)學(xué)劃分一致。

3 討論和結(jié)論

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

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

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

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