王 辰,杨娅富,李 立,逯金瑶,李嘉昊,李 波,3*

(1.东北林业大学野生动物与自然保护地学院,黑龙江哈尔滨 150040；2.湖南省生物多样性保护中心,湖南长沙 410000；3.国家林业和草原局野生动植物检测中心,黑龙江哈尔滨 150040)

穿山甲是哺乳纲(Mammalia)鳞甲目(Pholidota)穿山甲科(Manidae)动物的统称。现存3属8种：穿山甲属(Manis)的中华穿山甲(M.pentadactyla)、印度穿山甲(M.crassicaudata)、马来穿山甲(M.javanica)和菲律宾穿山甲(M.culionensis)；地穿山甲属(Smutsia)的南非穿山甲(S.temminckii)和大穿山甲(S.gigantea)；长尾穿山甲属(Phataginus)的树穿山甲(P.tricuspis)和长尾穿山甲(P.tetradactyla),分布于亚洲和非洲的热带和亚热带地区[1-6]。

穿山甲因有食用和药用价值而遭到人类的乱捕滥猎，加之生境被破坏及外来物种入侵等因素,其野生种群数量急剧下降,甚至面临灭绝[7-8]。2016年10月,CIETS缔约国大会通过“所有8种穿山甲物种由附录II 提升至附录 I”的提案[9-10]。由此,为打击穿山甲的非法贸易需要选择合适的分子标记、构建标准的物种鉴定方法。

线粒体基因组作为“扩展条形码”[11],常用于物种鉴定和系统发育关系评估[12-13]。该技术成本较为昂贵,不适合在穿山甲非法贸易频发的东南亚和非洲国家推广使用。线粒体细胞色素b(Cytochrome b,Cytb)基因进化速率适中、含有丰富的遗传信息、结构清晰,被广泛地应用于物种鉴定[14]。目前,基于Cytb基因的穿山甲物种鉴定多是利用脊椎动物或哺乳动物通用引物,所扩增的序列长度长短不一[15-18]。研究表明，利用穿山甲科特异引物可以有效排除非目标物种的干扰[15],Cytb基因序列长度变化可能影响鉴定结果的判读[19]。由此,笔者依据Cytb基因设计了穿山甲科物种鉴定的特异引物,并评估4种长度的Cytb基因序列对穿山甲物种鉴定的影响。

1 材料与方法

1.1 材料收集191份穿山甲样品。其中,174份为已知马来穿山甲样本,包括128份肌肉样,18份皮肤样,12份粪便样,7份毛发样,9份鳞片样,剩余17份样本为未知穿山甲样本,其中包括1份疑似穿山甲鳞片粉末(表1)。利用智人(Homosapiens,n=3)、貉(Nyctereutesprocyonoides,n=3)、朱顶雀(Acanthisflammea,n=3)、东北兔 (Lepusmandshuricus,n=1)、牛(Bostaurus,n=1)、虎(Pantheratigris,n=1)、家猪(Susscrofadomestica,n=1)的DNA作为非穿山甲物种对照样本。上述样本均来自国家林业和草原局野生动植物检测中心样本库。样品采集后于-20 ℃冰箱内保存。同时,从NCBI上下载135条序列穿山甲物种(表1)和部分非穿山甲物种的Cytb基因同源序列用于引物设计和序列综合分析。

1.2 方法

1.2.1DNA提取。穿山甲鳞片用75%乙醇清洗表面,经灭菌去离子水清洗、风干后剪取其表面残留的皮肤组织或毛发,或钻取鳞片粉末。将皮肤、毛发及鳞片粉末分别装入1.5 mL离心管。将肌肉样剪碎后置于1.5 mL离心管。使用AxyPrepTM Multisource Genomic DNA Miniprep Kit试剂盒(爱思进)提取肌肉、毛发、皮肤及鳞片的总DNA。另外,称取180～220 mg的粪便使用QIAamp Fast DNA Stool Mini Kit试剂盒(爱思进)提取总DNA。DNA提取步骤均参考试剂盒说明书,DNA原液存放于-20 ℃冰箱。

表1 穿山甲科试验样本及GenBank来源序列信息Table 1 Experimental samples and sequence information downloaded from GenBank for pangolin

1.2.2引物设计。利用Primer PREMIER 6.0软件基于马来穿山甲Cytb全序列设计穿山甲科特异PCR引物。引物搜寻参数包括：熔解温度(Tm)为60 ℃、引物长度(Length)为18～30 bp、扩增产物长度(amplication length)为70～200 bp,余下参数为默认值。使用MEGA 5.0中的BLAST程序筛选满足设计原则的引物,参考的非穿山甲物种的同源序列包括：兔形目的穴兔(Oryctolaguscuniculus,HQ596486、KT626640)、食肉目的土狼(Canislatrans,KT447695～KT447698)和家犬(CanislupusfamiliarisGU249573、MH714782～MH714784、KJ660981、KJ660982)、翼手目的蝙蝠(Ardopsnichollsi,KJ024748～KJ024751)、偶蹄目的牛(Bostaurus,KT626620)、贫齿目的犰狳(Dasypusnovemcinctus,KU253494、KT626619、KT626620)和灵长目的智人(Homosapiens,JN034134 ～JN034136)。

1.2.3PCR扩增和测序。利用3对引物：L14724/H15149[26]、mcb398/mcb869[27]、DCW-F1/DCW-R1分别扩增穿山甲科自有样本(表1)Cytb基因部分序列。反应体系均为50 μL：5 μL DNA模板,0.5 μL上、下引物,25 μL DNA 聚合酶(Transgen,北京),19 μL灭菌去离子水。PCR在PE9700型DNA扩增仪(PE,美国)进行。L14724/H15149反应条件为：预变性95 ℃/5 min(变性95 ℃/30 s,退火56 ℃/30 s,延伸72 ℃ 1 min),循环40次,延伸72 ℃ 5 min；mcb398/mcb869反应条件为：预变性94 ℃ 5 min(变性94 ℃ 1 min,退火52 ℃ 30 s,延伸72 ℃ 30 s),循环35次,延伸72 ℃ 10 min； DCW-F1/DCW-R1反应条件为：预变性94 ℃ 5 min,(变性94 ℃ 1 min,退火53 ℃ 30 s,延伸72 ℃ 30 s),循环35次,延伸72 ℃ 10 min。扩增产物经1.5%琼脂糖凝胶电泳检测后送至北京擎科生物技术有限公司哈尔滨分公司利用ABI3730XL DNA测序仪进行双向测序。同时将DNA浓度 ≥25 ng/μL的样本送至深圳华大基因科技有限公司,通过鸟枪法构建文库,使用Illu mina HiSeq2500测序仪进行高通量测序。

1.2.4引物有效性测试。引物特异性测试使用NCBI/Primer-BLAST程序验证引物DCW-F1/DCW-R1与GenBank 数据库中非目标物种的结合情况。然后使用7个非穿山甲物种的DNA样品进行实证试验。最后开展混合样品试验,将上述7个物种的DNA样本按浓度均匀混合,再分别与3种不同类型(肌肉、粪便和鳞片)的穿山甲物种DNA抽提液混合,浓度比为1∶1,分别进行PCR扩增及测序。

引物灵敏度测试在肌肉、毛发、组织、鳞片粉及粪便样品中各取2个样品,将DNA抽提液用灭菌去离子水从10 ng/μL依次稀释至1.0、0.1 ng/μL和10、1 pg/μL,使用“1.2.3”中的方法进行PCR扩增及电泳检测。上述引物的有效性测试中单次的PCR中均包含阳性对照(已知马来穿山甲DNA)和阴性对照(灭菌去离子水)。

1.3 数据分析利用MegAlign将获得的二代测序结果转成Fasta格式,参考NCBI上的同源序列使用MEGA5.2[28]手工校对；利用SeqMan软件拼接一代测序获得的Cytb基因序列,并结合测序峰图手工校对。利用DNASP V.6.0.70计算该研究所得序列的单倍型。使用MEGA5.2软件中的clustalW对序列排序并校正,计算基于Kimura-2-Parameter双参数(K2P)模型[29]的各穿山甲的种内和种间遗传距离,用boot-strap test 作1 000次置信度分析。使用GenBank数据库中的BLAST工具,对所有序列进行比对,并进行同源性分析比对,以确定种属。通过PhyloSuite V.1.2.2[30]基于最大似然法(mL)和贝叶斯法(BI)重建系统发育树。利用ModelFinder(选取BIC标准,模型分别选择IQ-TREE和Mrbayes)得到最佳分区方案和核苷酸替换模型。应用IQ-TREE和Mrbayes分别重建 mL和BI系统发育树(兔形目的穴兔同源序列为外群)。

2 结果与分析

2.1 引物有效性设计得到的引物 DCW-F1：5′- CCTATTCGCATACGCAATC-3′；DCW-R1：5′- AAGGATAGAGGCTACTTGTC-3′。该引物经NCBI/Primer-BLAST程序检索,与中华穿山甲同源序列的结合率为90.91% ～100%,其余7种穿山甲同源序列的结合率 ≥95%。实证试验表明，该引物未扩增出任何非目标物种条带。所有混合样品的扩增子测序图谱均为单一的可识别峰图,所得序列与混合前试验所得序列一致。1份疑似穿山甲鳞片粉末经引物mcb398/mcb869扩增,测序图谱为难以识别的混合峰图；而经引物DCW-F1/DCW-R1扩增,得到单一、可识别测序峰图。这些结果都说明该引物具有较高的特异性。灵敏度测试电泳结果：肌肉样DNA稀释至10 pg/μL时有可见的电泳条带；毛发和鳞片样DNA可稀释至0.1 ng/μL,皮肤和粪便样DNA可稀释至1 ng/μL(图1)。这表明该PCR系统的灵敏度足以检测这5种不同类型的样品。

2.2 序列分析该研究自测4组Cytb基因序列的单倍型数量随着序列长度变短而减少。具体是：106条MJ-L(1 140 bp)中存在70个单倍型,160条MJ-M(762 bp)中存在69个单倍型,160条MJ-S(421 bp)中存在28个单倍型,60条MJ-DCW(211 bp)中存在19个单倍型。

4组序列中所有已知穿山甲样本序列与马来穿山甲同源序列的相似度均大于99%。除MJ-M中样品PA-M-1、MJ-DCW中样品MP-DCW-1与中华穿山甲同源序列的相似度略小于95%(94.88%、93.27%)；MJ-DCW中部分样本与树穿山甲同源序列的相似度略大于95%(95.73%),其余样品与对应的穿山甲同源序列的相似度均大于98%(表2)。其中，前3组序列(1 140、762、421 bp)与已知4种穿山甲同源序列的相似度值变化很小,但最短的MJ-DCW中除马来穿山甲外,其余3种穿山甲同源序列的相似度有变小的趋势。

MJ-L、MJ-M、MJ-DCW中种内遗传距离最大值和最小值分别来自树穿山甲(0.113 6、0.117 4、0.149 4)和马来穿山甲(0、0、0.004 8)； MJ-S中种内遗传距离最大值和最小值分别来自树穿山甲(0.123 9)和中华穿山甲、马来穿山甲(0)(表3)。南非穿山甲和大穿山甲的前3组序列中,种间遗传距离最小值是种内遗传距离最大值的10倍以上；印度穿山甲和长尾穿山甲的前3组序列中,该值为5～10倍；马来穿山甲的MJ-DCW和树穿山甲的MJ-S、MJ-DCW中种间遗传距离最小值小于种内遗传距离最大值；其余种穿山甲的种间遗传距离最小值是种内遗传距离最大值的1～5倍(图2)。其中马来穿山甲、大穿山甲和树穿山甲的前3组序列中种内遗传距离无明显变化,MJ-DCW中明显变大,其余种穿山甲的4组序列变化很小；印度穿山甲/大穿山甲、印度穿山甲/长尾穿山甲、中华穿山甲/马来穿山甲、中华穿山甲/长尾穿山甲、马来穿山甲/长尾穿山甲、马来穿山甲/大穿山甲前3组序列的种间遗传距离变化很小,MJ-DCW中明显变小；南非穿山甲和长尾穿山甲中,MJ-S显著大于MJ-L、MJ-M；其余种穿山甲4组序列的种间遗传距离变化很小。

注：A1、A2为毛发样；B1、B2为肌肉样；C1、C2为皮肤样；D1、D2为鳞片样；E1、E2为粪便样。DNA稀释浓度为10 ng/μL至1 pg/μL；M.DL1200 Maker(1 200、900、700、500、300、100 bp)。Note:A1 and A2 are hair samples, B1 and B2 are muscle samples, C1 and C2 are skin samples, D1 and D2 are scale samples, and E1 and E2 are faecal samples.The diluted concentration of DNA is from 10 ng/μL to 1 pg/μL.M is the DL 1 200 marker (1 200，900，700，500，300，100 bp).图1 灵敏度测试的琼脂糖凝胶电泳Fig.1 Agarose gel electrophoresis of sensitivity test

表2 BLAST比对结果Table 2 BLAST results

4组序列的ML和BI系统发育树均显示出了支持良好的拓扑结构(图3)。穿山甲物种的分化节点具有较高的置信度(60% ～100%)(图3～4)。MJ-L、MJ-M、MJ-S中8种穿山甲可被分为2个大支：中华穿山甲、菲律宾穿山甲、马来穿山甲和印度穿山甲为一支,南非穿山甲、大穿山甲、长尾穿山甲及树穿山甲为一支。而MJ-DCW中8种穿山甲的拓扑结构发生改变,南非穿山甲和大穿山甲脱离非洲穿山甲分支并入亚洲穿山甲分支(图4)。这表明随着序列长度变短,8种穿山甲的分支会逐渐合并。菲律宾穿山甲和MJ个体与已知马来穿山甲聚为一支,PA-1 和MP-DCW-1与已知中华穿山甲聚为一支,PA-3 和SG-DCW-1、2与已知大穿山甲聚为一支,PA-6、10、11和 MT-DCW-1-7、SJ -2018-72-2、SJ-2018-72-3与已知树穿山甲聚为一支。另外，有3条序列被证明分类错误：原标注为中华穿山甲的JN411577应为马来穿山甲，原标注为中华穿山甲的NC_016008应为马来穿山甲，原标注为长尾穿山甲的NC_004027应为树穿山甲。

表3 4组序列的种间遗传距离和种内遗传距离Table 3 The intra-d and inter-d based on 4 groups of sequences

注：CR.印度穿山甲,PE.中华穿山甲,JA.马来穿山甲,GI.大穿山甲,TE.南非穿山甲,PH.长尾穿山甲,TR.树穿山甲。Note:CR.M.crassicaudata, PE.M.pentadactyla, JA.M.javanica, GI.S.gigantea, TE.S.tem minckii, PH.P.tetradactyla, TR.P.tricuspis.图2 4组序列种内遗传距离和种间遗传距离差异分布Fig.2 Distributions of the intra-d and inter-d on 4 groups sequences

3 讨论

3.1 引物有效性穿山甲物种鉴定中利用脊椎动物或哺乳动物通用引物包括F-/R-primer[15]、GVL14724/H15149[17-18]、L14724/H15149[20]、mcb398/mcb869[21]、 L15601/H15748[31],它们具有可识别物种范围广、适用性好的特点 。但上述通用引物难以鉴别混有其他非穿山甲物种的疑似检材，部分引物因扩增片段较长，对检材的质量要求偏高。该研究所设计的穿山甲科特异引物DCW-F1/DCW-R1,通过实证试验、混合试验和盲测试验均得到与目标片段一致的DNA序列,未产生非目标物种条带。这表明该引物能有效、特异性识别穿山甲科物种。同时,穿山甲的非法贸易经常涉及肌肉、皮肤、鳞片样本,野外调查时可以收集非损伤性样本(如粪便和毛发)。上述涉案和野外调查收集样本可能会存在不同程度降解,导致抽提的DNA质量不佳[32]。新设计引物的扩增产物长度仅为211 bp,灵敏度试验结果也证明其良好的适用性,可以满足不同来源的5种类型样本。

3.2 序列分析与判别标准Cytb基因是研究种间和种内遗传分化程度的良好指标,常用于物种鉴定[33]。Cytb基因序列长度变化可能影响鉴定结果的判读[19],因此通常使用2种以上的序列分析方法(如BLAST比对法、遗传距离法、系统发育树法)开展物种鉴定[34]。通常,BLAST比对法中当待测样本与已知单一物种同源序列相似度高于95%时可以判定为同一物种[35]。该研究只有2份样本与已知中华穿山甲同源序列的相似度低于该鉴定标准,这可能与已知中华穿山甲Cytb基因序列单倍型数量偏低有关。其实脊椎动物物种鉴定中上述判别标准并不是唯一的[36]。鉴于多数样本与相应穿山甲物种同源序列相似度均高于98%,建议使用相似度98%作为穿山甲物种判别标准。

穿山甲科物种中除马来穿山甲的MJ-DCW和树穿山甲的MJ-S、MJ-DCW外，剩余的种间最小遗传距离均大于种内最大遗传距离,符合遗传距离法物种鉴定判别标准[37]。其中,不同种穿山甲、同种穿山甲的不同组序列间的样本量差异较大,且随着序列长度变短,序列相似性变大,单倍型数量逐渐减少,物种间样本量的不平衡会导致遗传距离法结果不稳定[38],使得穿山甲种内遗传距离和种间遗传距离数值差异较大。此外,树穿山甲的种内遗传距离均远大于其他种穿山甲,高于0.05,且在MJ-S、MJ-DCW中种间最小遗传距离大于种内最大遗传距离,这表明在树穿山甲中极有可能存在亚种[39-41]或隐蔽物种[18]。

系统发育树(BI树和ML树)法鉴别物种的判别标准：同一物种的不同个体聚为单系分支[42]。该研究中4组序列构建的系统发育树,待测样品均与已知物种的同源序列分别聚为单系(除马来穿山甲)。过去菲律宾穿山甲常被认为是马来穿山甲的亚种,分布区重叠因海平面上升导致其被地理隔离[3],马来穿山甲种群具有较高的遗传多样性。且该研究中菲律宾穿山甲的样本量较少,需要进一步增加其Cytb基因单倍型,以确定它们之间的关系[26]。序列长度发生改变时系统发育树部分分支的置信度值，甚至整个拓扑结构都会产生一定差异[43]。前3组序列构建的系统发育树将亚洲穿山甲和非洲穿山甲分为独立分支,这与线粒体基因组构建的系统发育树[11]结构相同；而MJ-DCW构建的系统发育树则与之不同；且4组序列随着序列长度变小,其各物种分支的置信度有变小的趋势。 因此,利用系统发育树法鉴别物种时要考虑序列长度对结果判读的影响。

注：节点后验概率/自举值从上至下依次为MJ-L、MJ-M、MJ-S；不包含外群。Note:The posterior probability/bootstrap value are MJ-L，MJ-M，MJ-S from top to bottom;out groups are not included in the phylogenetic tree.图3 基于MJ-L构建系统发育树Fig.3 Phylogenetic tree based on MJ-L

图4 基于MJ-DCW构建系统发育树(不包含外群)Fig.4 Phylogenetic tree based on MJ-DCW (exculded out groups)