基于DNA序列信息的贞琵甲属物种界定分析(鞘翅目:拟步甲科:琵甲亚科)

李秀敏,戎际达,刘一右,范娇娇,任国栋

(河北大学生命科学学院,河北大学生命科学与绿色发展研究院,河北省动物系统学与应用重点实验室,河北保定 071002)

青藏高原是全球生物多样性热点地区之一,以拥有大量的特有物种而闻名(Mittermeieretal., 2005)。贞琵甲属AgnaptoriaReitter是青藏高原的特有属,隶属于鞘翅目Coleoptera拟步甲科Tenebrionidae琵甲亚科Blaptinae琵甲族Blaptini,主要分布在3 000 m以上的高海拔地区,是琵甲族小琵甲亚族中的第三大属,模式种为AgnaptoriarubripesReitter,1887。目前该属已知36 种/亚种,分布于中国西藏、青海、四川和甘肃4省(Reitter, 1887, 1893;Medvedev, 1998, 1999, 2002, 2006, 2008;石爱民等, 2005, 2007; Ren and Liu 2009; Renetal. 2016)。该属的主要形态及生物学特征如下:体中小型,黑至红棕色;第1～2节触角红色/棕色;前胸背板横宽;鞘翅中部两侧平行;足细长,胫节自基部向端部强烈变宽,基部通常为红色或红褐色,跗节黑色,跗节短于胫节;后翅愈合,不能飞行,主要栖息于石块或腐殖质下,夜出性活动,遇到惊扰时,能从腹部特殊的腺体内散发出强烈的刺鼻性气味。

随着全球气候变化等环境因子的改变,物种正在以惊人的速度消失,而越来越多的研究表明生物多样性丧失对生态系统的功能有着巨大影响(Cardinaleetal., 2012)。鉴于高海拔地区的物种分化强烈,以及基于传统的形态特征进行分类的局限性,加快物种发现和认知尤为重要,特别是生物多样性热点地区的物种多样性。近年来,DNA序列信息也被广泛应用于无先验物种信息的物种界定,在越来越多的类群中得到广泛应用,确立了许多物种的分类地位和新种地位(Yeatesetal., 2011; Carstensetal., 2013; Toussaintetal., 2015; Soldatietal., 2017)。尤其是在形态上相近的物种间利用DNA序列信息来划定物种边界,验证形态物种的有效性,已成为快速识别和描述物种、加速发现和认识种类繁多的昆虫世界的一种有效手段。分子物种界定方法结合形态证据进行准确的物种界定,以揭示某一类群的系统发育关系,从新的视角理解生物多样性时空演变格局的形成过程。因此,本研究选取贞琵甲属15个种的30个样品,选择4个基因片段(COI; Cytb; 16S rDNA; 28S rDNA-D2)作为分子标记,构建该属部分物种的系统发育关系,运用分子证据对以形态特征建立的已知物种分类地位进行验证,为该属未知种类的认知以及其物种多样性分布格局的形成和演变奠定基础。

1 材料和方法

1.1 实验材料

本研究选取贞琵甲属15个物种的30个样本作为研究对象,主要来自于近几年野外采集的酒精浸泡标本及部分馆藏干制标本。首先在Nikon SMZ800立体显微镜下观察其形态特征,根据已出版的《中国动物志》拟步甲科(一)及相关分类文献依据其形态特征进行分类鉴定(任国栋等,2016)。所有样本的采集信息见表1。

表1 所用贞琵甲属样本信息

1.2 DNA提取、PCR扩增及测序

选取中、后足肌肉组织,使用EZNA©昆虫DNA提取试剂盒(Omega Bio-tek, USA),参照说明书中的操作步骤提取DNA,获得的DNA存放于-20℃冰箱中保存用于PCR扩增。本文选取了4对线粒体和核基因进行扩增,引物由通用生物系统(安徽)有限公司合成(Simon, 1994; Folmeretal., 1994; Belshawetal., 1997; Simmonsetal., 2001)。COI基因的引物序列为2183: 5′-CAACA TTTATTTTGATTTTTTGG-3′;3014: 5′-TCCAA TGCACTAATCTGCCATATTA-3′。Cytb基因的引物序列为revcb2: 5′-TGAGGACAAATATCATTTTGA GGW-3′; rebcbj: 5′-TCAGGTCGAGCTCCAATTC ATGT-3′。16S基因的引物序列为13398: 5′-CG CCTGTTTATCAAAAACAT-3′;12887: 5′-CCG GTCTGAACTCAGATCAT-3′。28S基因的引物序列为3665: 5′-AGAGAGAGTTCAAGAGTACGTG-3′; 4068: 5′-TTGGTCCGTGTTTCAAGACGGG-3′。PCR采用TaKaRa Ex Taq进行PCR扩增(TaKaRa, 中国大连),反应体系均为25 μL,PCR Mix 12.5 μL,10 μmol/L上下游引物各1 μL,DNA模板1.5～5 μL,ddH2O补足至25 μL。PCR反应条件为:94℃预变性4 min;30个循环包含94℃变性1 min 30 s,50～58℃复性30 s,72℃延伸1 min,最后72℃延伸8 min。PCR扩增产物经1%琼脂糖凝胶电泳检测,并用凝胶成像系统拍照。确定为目标条带的PCR产物由通用生物系统(安徽)有限公司完成双向测序,获得DNA序列信息。

1.3 序列数据处理

测序得到的DNA序列首先使用DNASTAR SeqMan v7.1.0软件查看测序峰图并进行拼接,为了确保序列的准确性,对照双向测序峰图,逐一对所有碱基进行人工校对剪切,确保读取正确,然后将得到的序列在NCBI(https://www.ncbi. nlm.nih.gov/)中BLAST搜索进行同源性分析,确定所得到的序列是否为目标序列,并保证序列方向一致。同一基因的目标序列按照分析需要通过BioEdit v7.0.9.0(Hall, 1999)软件导入到一个.fas文件中,然后将.fas文件导入软件PhyloSuite中,用MAFFT进行序列比对;然后用Concatenate Sequence将4个基因片段串联起来;用PartitionFinder2(unlingked)进行分区,选择最适合的分区方案;用ModelFinder对优化的数据进行最优模型选择用于系统发育分析研究(Lanfearetal, 2017; Kalyaanamoorthyetal, 2017)。

1.4 系统发育树构建

系统发育树的构建采用目前公认的最大似然法(ML),选用与贞琵甲属关系较近的小琵甲属的9个物种作为外群。ML分析使用了IQ-TREE v1.6.6(Nguyenetal., 2015)在线软件,4个基因片段的联合数据集在Auto+R6替代模型下进行运算,bootstrap分析选用Ultrafast,重复抽样次数设为默认值为1 000,启用SH-aLRT,设置为默认值1 000,快速自举值(uBV)大于或等于95%被认为是强有力的支持。

1.5 分子物种界定分析

对于分子物种界定分析,最近研究显示有些方法可能会低估或高估物种的数量(Dellicouretal., 2018; Luoetal., 2018)。因此,综合使用多种分子物种界定方法,以更准确地评估物种的分类地位。本研究使用如下3种方法的8种方案分析物种界限:

(1)ASAP(Assemble Species by Automatic Partitioning)在线分析软件(Puillandreetal., 2012;https://bioinfo.mnhn.fr/abi/public/asap/asapweb. html)。为了避免与缺失数据相关的距离估计偏差,本文仅对COI和Cytb的单基因进行了ASAP物种界定分析。将比对好的.fas文件导入程序,基于K2P模型计算遗传距离,其余参数使用缺省值,即种内差异先验值P(prior intraspecific divergence)为0.001～0.1,将最小相对间隙宽度(relative gap width)X值设置为0.5,其余参数为系统默认设置。

(2)GMYC(Generalized Mixed Yule Coalescent)物种界定方法(Ponsetal., 2006;Reidetal., 2012)。先利用BEAST v2.4.7(http://www. beast2.org/)获得超度量(ultrametric)二叉分支树,使用jModelTest选择最佳核苷酸替代模型。首先,运行了1个联合基因数据集的BEAST分析(去掉所有外群),然后又对COI和Cytb数据集进行了单基因的BEAST分析(设置相同)。利用BEAUti打开FAS格式文件,依次约束内群为单系群、设置核苷酸替代模型、分子钟、系统树、MCMC链等,并生成xml文件;BEAST下运行xml文件;Tracer v1.6对BEAST运算结果进行分析,Burn-In默认为链长10%,全部ESS>200表示数据运行很好;再利用TreeAnnotator v1.8.2生成所需系统树;转换的系统发育树在R语言包“splits”读取并计算得到物种界定结果,以与其它分子界定方法、形态鉴定结果对照分析。

(3)PTP(Poisson Tree Processes)物种界定方法(Zhangetal., 2013)。本研究使用PTP在线程序运行(https://species.h-its.org/ptp/),共分为3组数据进行PTP分析:去掉所有外群后的联合基因的系统发育树、COI基因系统发育树和Cytb基因系统发育树。运行时将MCMC代数设置为100 000代,其它参数使用默认值。

2 结果与分析

2.1 系统发育关系

基于贞琵甲15个物种的30个样本的4个基因片段:COI,Cytb,16S rDNA和28S rDNA的DNA序列数据集在IQ-TREE软件下构建了ML系统发育树(图1)。总的来说,贞琵甲属15个物种中的13个物种的种间关系得到了较高的支持(80%的节点uBV ≥ 95%),与外群小琵甲属部分物种构成了2个分支,内群分支的单系性也得到了很好的支持(uBV=100)。

图1 基于30-taxa的3个线粒体基因和1个核基因序列(COI, Cytb, 16S rDNA, 28S rDNA),2 327 bp的串联数据集获得的贞琵甲属Agnaptoria的系统发育树和物种界定结果Fig.1 Phylogenetic relationships within the genusAgnaptoria and species boundaries. Phylogenetic tree based on mitochondrial and nuclear DNA sequences. A similar topology was recovered under a maximum likelihood approach.注:树的每个节点的支持度由大似然法推断得到的超快bootstrap值(uBV)表示。彩色柱形条带表示形态物种数(黑色),以及8种物种划分方案的分子物种界定分析结果(红色、橙色、紫色、蓝紫色、粉色、绿色、橘黄色和蓝色);对于单基因COI/Cytb的界定结果中,灰色柱形条带表示缺失的序列信息;21个MOTUs为ASAP、GMYC和PTP分子物种界定分析得出的结果。Note: Support for each node was represented by ultrafast bootstrap values (uBV) obtained for the maximum likelihood tree. Vertical coloured bars delineate extant morphospecies (black colour), and the results of the eight separated molecular analyses delimiting species (red, orange, purple, bluish-purple, pink, green, yellow, blue). For the analyses using COI or Cytb genes only, we used grey colour to delineate specimens for which sequencing failed. Numbers on the right represented codes for 21 putative MOTUs delineated using ASAP, GMYC and PTP analyses.

2.2 物种界定结果

利用3种分子物种界定分析方法(ASAP、GMYC和PTP)对单基因片段(COI和Cytb)和联合基因基因(COI,Cytb,16S rDNA和28S rDNA)片段的数据集对30个贞琵甲样品进行了物种界定,其中GMYC和PTP两种物种界定方法既使用了联合基因数据集,也使用了单基因数据集,而ASAP方法只使用了单基因数据集。3种分子物种界定结果(共8种方案)推断出的MOTUs的数基本一致,PTP方法的3种方案界定出的MOTU数以及与其它两种方法在部分MOTUs上显示出不同的结果,且界定的物种数略多。通过形态学鉴定,30个样本为13个已知物种和2个待定种,其中13个物种的形态鉴定结果与分子界定结果完全吻合(图1)。

ASAP方法对单基因COI和Cytb的界定结果一致,为确定最终物种的数量,本研究倾向于先验值在0.01～0.02之间,30个样品界定为17个物种;GMYC方法使用单基因与联合基因界定的结果基本一致,为17～18个物种(COI基因界定出17个。Cytb基因界定出17个,一个物种基因缺失;联合基因界定出18个),不同的是COI基因将MOTU19与MOTU20界定为1个物种,而Cytb基因和联合基因数据集界定为2个物种。

PTP方法的3种方案界定出的MOTU数以及与上述两种方法的在部分MOTU上显示出不同的结果,且界定的物种数略多于上述两种方法(COI基因界定出19个物种;Cytb基因界定出19个物种,一个物种基因缺失;联合基因界定出18个物种)。Cytb基因将MOTU1、MOTU2和MOTU3界定为3个物种,而Cytb基因和联合基因数据集将其界定为1个物种;Cytb基因将MOTU17、MOTU18、MOTU19和MOTU20界定为一个物种,COI基因与其它7种方案结果不同,将MOTU17和MOTU18界定为2个物种,而将MOTU19和MOTU20界定为1个物种。

3 结论与讨论

不同的分子物种界定方法界定物种的结果并不完全一致,但是无论使用单基因还是联合基因,分子种界定结果与形态分类结果是大致吻合的,对已知物种分类地位的验证以及提高物种鉴定效率有很大帮助。29个样品的COⅠ基因序列通过ASAP物种界定方法得到了18个MOTUs,通过GMYC方法得到了17个MOTUs,通过PTP方法得到了17个MOTUs;28个样品的Cytb基因序列通过ASAP和GMYC分子物种界定方法均得到了18个 MOTUs,通过PTP方法得到了17个MOTUs。除PTP方法中的Cytb基因外,其它7种方案都将A.markana的5个样品界定为一个MOTU;这一结果与已有研究认为的GMYC和PTP方法会高估MOTUs的数量较为一致(Lecocqetal., 2015; Dellicour &Flot, 2018; Luoetal., 2018)。除PTP方法中的Cytb基因的界定结果与形态鉴定结果一致外,其它7种方案均将A.nigriceps的4个样品界定为2～3个MOTUs。通过检视不同采样地点的黑头贞琵甲A.nigriceps,确认其形态特征并无差异。综合考虑,本研究将来自不同采样地点的样品暂时作为一个物种处理。总的来说,ASAP、GMYC和PTP方法的8种方案界定结果与形态鉴定结果均呈现较高的吻合度,证明了这3种方法综合运用对贞琵甲属进行分子物种界定的有效性。

青藏高原是印度板块与欧亚板块两大陆地板块的剧烈碰撞形成的地球上海拔最高、最年轻的高原。随着青藏高原隆升过程中的地质变化,受到山脉、河流的阻隔,形成“岛屿”化的地理格局,某一种群被限制于特定高度分布并在特定的小尺度生境内繁衍,并以小种群形式存在,从而限制了种群间的基因交流,为遗传物质的分化——新物种的快速形成创造了机会(Lietal., 2021)。青藏高原这种极端环境下物种形成时间较短,种间识别特征差异不显著,仅依据传统的形态特征很难将其区分开,分子界定结果通常也会出现多于形态物种的现象(Lietal., 2021)。物种的形成是一个十分复杂的过程,其形成模式一直受到广泛关注(Futuyma and Mayer, 1980; Taylor and Friesen, 2017)。对分化期较短的物种,个体或种群间因生殖隔离而产生的遗传距离上的差异,因其极小的形态差异,从而可能导致基于进化树和遗传距离的分子物种界定结果与传统的形态鉴定不完全吻合。青藏高原特异性环境下产生的一些“疑难种”在很多类群中都普遍存在,对于这些物种的界定边界在哪里,还需继续积累更多不同采样地点样本的数量,检视大量标本找出其形态变化规律,并通过增加基因数量或者通过基因组数据来进行其遗传距离的分析,使研究结果尽可能符合客观事实,反应物种真实的分类地位。

运用DNA序列与传统的形态学方法相结合鉴定物种,是准确确定物种分类地位的科学方法,脱离传统分类学而单纯开展分子物种经鉴定得出的结果并不完全可靠,因此,传统的形态分类学方法与分子鉴定必须结合,才能准确确定物种的分类地位,提高物种鉴定的效率和准确性。本文初步构建了贞琵甲属部分物种的DNA序列信息数据库,为贞琵甲的快速、准确鉴定提出了可行方法,对贞琵甲的物种演化和地理分布格局研究也具有重要意义。