线粒体DNA分析青藏高原东方蜜蜂(Apis cerana)种群遗传结构

刘意秋，周丹银，龚雪阳，冯毅楠，和绍禹*，赵文正*

(1.云南农业大学 动物科学技术学院，云南 昆明 650201；2.云南农业大学 东方蜜蜂研究所，云南 昆明 650201)

种群遗传结构分析的主要内容为：种群内及种群间等位基因/基因型的遗传差异和频率分布。分析结果可为探索种质资源的演化形成过程，定义重要进化单元[1]，规划其保护、发展和利用提供重要基础数据[2]。东方蜜蜂(ApisceranaFabricius，1793)是亚洲的土著蜜蜂种类，其自然分布地域遍及亚洲大陆温带、热带、亚热带地区及其众多岛屿[3]。基于个体数量庞大、适应性和抗逆性较强的特点，东方蜜蜂为其生境内的生态系统提供着十分重要的传粉服务，对于当地生态系统多样性的维持和稳定具有不可替代的作用[4]。此外，在一些植被丰富但耕地稀缺的山区、半山区，饲养东方蜜蜂还为当地居民的经济增收做出贡献。

线粒体DNA(mtDNA)是真核生物的核外遗传物质(之一)，由于具有母系遗传、不发生重组、种内变异丰富[5]等特点而被广泛应用于蜜蜂种群遗传结构分析。利用mtDNAtRNAleu-CO2基因片段对西方蜜蜂开展了大量调查研究[6～11]，这些研究结果与Ruttner(1988)利用形态学方法划分的西方蜜蜂(Apismellifera)种下分类体系基本一致。然而这段mtDNA片段并不十分适宜对东方蜜蜂进行分析，原因在于该片段长度相对较短[12]甚至在某些东方蜜蜂个体中完全缺失[13]。之前一些学者也利用细胞色素C氧化酶亚基1基因(COX1)[14]NADH脱氢酶亚基2基因(ND2)[15]等mtDNA片段对东方蜜蜂开展了遗传多样性分析。本试验中采用的是细胞色素C氧化酶亚基1-亚基2基因片段(CO1-2)。

我国境内的东方蜜蜂又称中华蜜蜂或中蜂。中华蜜蜂的种群遗传结构及遗传多样性状况前人已有研究报道，然而这些研究所得结论不尽相同。杨冠煌等对中华蜜蜂资源进行了详细的形态学和生物学调查，调查结果认为我国境内分布着5个东方蜜蜂亚种，其中滇、川、藏省区分布有西藏亚种(Apisceranaskorikovi)、印度亚种(Apisceranaindica)、阿坝亚种(Apisceranaabansis)和中华亚种(Apisceranacerana)，海南亚种(Apisceranahainana)只分布于海南岛上[16,17]。但Radloff等的形态学研究则认为我国境内的东方蜜蜂分为“阿坝”、“华中华东”和“华南”3个形态类群[18]。在后续的研究中，一些学者利用线粒体DNA[19,20]和微卫星DNA[21]分子标记已经找到了支持海南亚种(形态类群)的分子证据。然而对于青藏高原地区东方蜜蜂种群遗传结构及该区域西藏亚种和阿坝亚种(形态类群)间遗传分化的分子数据目前鲜见报道。鉴于此，我们采集了青藏高原地区(包括滇、川、甘、青、藏等省区)的东方蜜蜂样本并进行线粒体DNA分子标记检测和数据分析，以期调查该区域东方蜜蜂的种群遗传结构，并为研究东方蜜蜂西藏亚种和阿坝亚种(形态类群)的遗传分化积累分子数据。

1 材料与方法

1.1 样本采集与DNA提取

东方蜜蜂样本于2013年5月-2015年6月采集自中国青藏高原及其周边5个省、区的9个采样点，每个点采集7～30群蜜蜂样本，共237群。同一采样点的样本相互间距不超过10 km，样本详细信息见图1和表1。采样时从巢脾上直接抓取50只左右的工蜂，即刻保存于100%的乙醇中并尽快带回实验室贮存于-75 ℃条件下。针对每群蜜蜂样本随机选取1只工蜂，采用经典的酚/氯仿法[13]提取全基因组DNA。

1.2 PCR扩增与DNA测序

根据已发表的东方蜜蜂线粒体基因组全序列(GenBank登录号：NC_014295)设计了细胞色素C氧化酶亚基1-亚基2基因片段(CO1-2)的特异性引物：CO1-2-F(5’→3’) CCTCGACGATACTCAGATTATCC，CO1-2-R(5’→3’) TCCAAGTGATGGGACAGTTCAT，由上海生工生物技术有限公司合成。PCR扩增体系及扩增条件沿用之前的研究[22]。扩增产物交由上海生工生物技术有限公司进行正反向测序，测序引物即PCR扩增引物。

1.3 数据分析

利用Lasergene 7.0软件包中的SeqMan对正、反向DNA序列进行检查、校正、修剪并生成一致序列；从NCBI数据库下载了24条其它地区东方蜜蜂的CO1-2片段作为参考序列(见图3b)，利用MEGA 5.05软件[23]对本次试验测序获得的CO1-2一致序列和加入24条参考序列后的全部序列分别进行Alignment计算(在默认参数下)；使用DnaSP 5.1软件[24]进行单倍型种类、频率、单倍型多样度(Hd)、核苷酸多样度(π)以及种群间遗传距离(Fst)矩阵的统计计算；以GeographicDistanceMatrixGenerator v1.2.3生成种群间地理距离矩阵；使用Arlequin v3.11软件[25]进行Tajima’s D中性检验和分子方差分析(AMOVA)，并对种群间遗传距离Fst矩阵与地理距离矩阵之间实施Mantel检验；使用Network 5.03软件(http://www.fluxus-engineering.com)中的Median-Joining算法对鉴别出的单倍型进行网络关系构建；以苏拉威西蜂(Apiskoschevnikovi)同源片段为外群并加入参考序列后，分别利用邻接法(Neighbour-Joining)、最大简约法(Maximum-Parsimony)[26]和贝叶斯法(Bayesian-Algorithm)[27]对鉴别出的单倍型进行系统发育关系重建。邻接法主要执行参数为：Bootstrap Replications=10000，Substitution Model=Kimura 2-parameter model。最大简约算法主要执行参数为：collapse=maxbrlen，opt=acctran，hsearch start=stepwise，addseq=random，nreps=20，swap=TBR，multrees=yes，steepest=yes，bootstrap nreps=1000。贝叶斯算法主要执行参数为：lset nst=6，rates=invgamma，prset revmat=dirichlet(1,8,1,1,8,1)，statefreq(7,2,2,9)，mcmcp ngen=5000000，nruns=2，nchains=4，temp=0.1，burninfrac=0.3。

2 结果与分析

2.1 DNA序列变异

针对9个地理种群的237只东方蜜蜂样本，本试验测序获得了237条可用的线粒体DNA 序列，序列长度700～702 bp。通过与NCBI数据库中已公布序列进行比对确认为目标片段CO1-2。序列的A+T平均含量为83.4%，与谭宏伟报道的东方蜜蜂线粒体基因组A+T含量(83.96%)[28]相近。序列比对结果共发现30个变异位点，其中单一多态位点4个，简约信息位点26个，发现了两处1-2个碱基的插入/缺失。30个变异位点中包括24个转换位点5个颠换位点和1个既转换又颠换的位点，共形成35种单倍型，依次命名为T1-T35，35种单倍型的序列变异样式见附件。与24条参考序列合并进行分析时同源片段缩短为441 bp，序列比对发现59个变异位点，形成了30种单倍型。其中从样本237条序列中探测到17种单倍型，依次命名为Tib1-Tib17；从参考序列中共发现16种单倍型，3种同前(Tib1、Tib8、Tib10)，13种不同的单倍型依次命名为Ref1-Ref13。

2.2 地理种群内遗传多样性

9个地理种群中(除QHMH外)均发现一定比例(≥25%)的独有单倍型，全部样本的独有单倍型率为80%，而XZLZ地理种群的独有单倍型率为100%(表2)，表明各地理种群遗传结构的独特性较为显著；核苷酸多样度(π)指从群体中随机抽取的两条序列之间平均每个位点上的核苷酸差异数[29]，常用于量度种群内包含的遗传多样性。当一个群体内包含较多的相互间序列差异较大的个体时其π值相对较高，反之则π值相对较低。本次试验中全部样本的π值为0.003 92，比之前研究报道中广东、广西及湖南地区[20,22]东方蜜蜂种群的π值高，推测青藏高原地区可能蕴藏着比华南平原地区更丰富的遗传多样性；Tajima’s D中性检验结果显示没有任何种群呈现显著或极显著的正/负D值，各地理种群近期可能未曾经历明显的瓶颈效应或种群扩张。

表2 东方蜜蜂各地理种群遗传结构指标分析结果Table 2 Results of genetic diversity analysis for Apis cerana geographic populations

注：N样本量；Nh单倍型数；U 独有单倍型数；U % 独有单倍型率；Hd 单倍型多样度；π 核苷酸多样度；D Tajima中性检验D值；*P≤0.05；**P≤ 0.01

Notes: N sample size; Nh number of haplotypes; U number of unique haplotypes; U% percent of unique haplotypes; Hd haploytpe diversity; π nucleotide diversity; D Tajima’s D test; *P≤0.05; **P≤ 0.01

2.3 青藏高原地区种群遗传结构

方差分析结果显示：不论如何分组，群体内变异对总变异的贡献率始终最大(>45.5%)且变异固定指数均达到极显著水平(表3)。其中第5号分组方式的组间变异对总变异的贡献率相对最高(38.42%)且变异固定指数达到显著水平，表明按5号分组方式划分遗传类群相对最合理，可见不同地理种群间存在着较显著的遗传差异；单倍型网络关系图(图2)显示：共探测到7种共享单倍型，占全部单倍型种类(35种)的20%。且共享单倍型最多只出现于3个地理种群，未探测到广泛分布于9个地理种群的共享单倍型。表明青藏高原东方蜜蜂各地理种群间的基因交流十分有限。同一地理种群包含的不同单倍型间呈现相对较大的遗传差异，表现出较显著的种群内变异，这与方差分析的结果一致；Mantel检验结果显示：9个东方蜜蜂地理种群的遗传距离(Fst)矩阵与地理距离矩阵的相关性系数为0.113(P=0.259)，未达到显著水平，据此推测距离隔离并不是形成青藏高原地区东方蜜蜂种群遗传结构的主要因素。

表3 分子方差分析结果Table 3 Results of analysis of molecular variance (AMOVA)

注：分组规则，按行政区划或地理距离远近分组；Fct组间变异固定指数；Fsc组内群体间变异固定指数；Fst群体内变异固定指数；*P≤0.05；**P≤0.01。

Note: grouping rules, grouping according to administrative division or geographic distance;Fctfixation index of variation among groups;Fscfixation index of variation among populations within group;Fstfixation index of variation within populations; *P≤ 0.05;**P≤ 0.01.

图2 35种线粒体DNA单倍型的网络关系图Fig.2 Network relationship of the 35 mtDNA haplotypes注：不同颜色代表不同地理种群；圆圈代表推测存在但未被检测到的单倍型；方点代表突变步骤。Note：different colors represent different populations; circles represent inferred haplotypes that are detected here; square dots represent mutation steps.

与参考序列合并进行分析时获得了30种单倍型，其系统发育关系如图3所示：三种算法推断出的系统发育树拓扑结构完全一致，各大分支的支持率均大于60%。青藏高原地区东方蜜蜂与参考序列中云南、广西、日本、韩国及一条来自印度(KF889329，Ref-6)的东方蜜蜂聚为一支，彼此间遗传差异相对较小；与来自台湾岛、马来西亚、文莱、印度(其它)、俄罗斯、泰国等地区的参考序列分支支持率较高，呈现出相对较大的遗传差异。整体而言，单倍型间系统发育关系呈现显著的地理分布格局。

3 讨论与结论

本文采用的线粒体DNA片段长度约700 bp，从青藏高原地区东方蜜蜂样本中共检测出30个变异位点，占序列总长度的(30/700)4.28%。能够较充分地反映地理种群内/间存在的变异，从而有利于分析种群遗传结构、评价遗传分化程度。利用该片段探测到全部样本的单倍型多样度(Hd)和核苷酸多样度(π)分别为0.939和0.00392，表明青藏高原地区东方蜜蜂种群中蕴藏着较丰富的遗传多样性。

图3 a:30种单倍型的邻接法系统发育树；b:24条参考序列信息Fig.3 a:Neighbor-joining phylogenetic tree of the 30 haplotypes; b:Information about the 24 reference sequences(分支支持率依次为邻接法、最大简约法和贝叶斯法的自检举支持率)注：Tib1, Tib8, Tib10三种单倍型既出现于实验样本又出现于参考序列中；Ref单倍型只出现于参考序列中(Branch support values are bootstrap values from Neighbor-joining、Maximum parsimony and Bayesian inference in turns)Note: Tib1, Tib8 and Tib10 present in both samples of present study and reference sequences; Ref haplotypes present only in reference sequences

西藏地区XZDJ、XZLZ和XZMK采样点检测到的一些单倍型(T23、T25、T27、T29、T33)与四川SCMEK采样点所检测到的单倍型(T14和T15)相互间遗传差异较大，这两个区域恰是东方蜜蜂西藏亚种和阿坝亚种(形态类群)[16,18]的分布区域。这可能是两个亚种间呈现遗传分化的分子证据，对此还有待进一步开展细致的采样与研究。

分子方差分析(表3)结果表明：种群内变异在总变异中所占比重始终最大。Network图(图2)亦显示：地理种群内包含着相互间遗传差异较大的单倍型，GSMX和SCMEK两个种群尤其如此。经与当地工作人员交流了解到：近年来当地曾购买了一些外地的东方蜜蜂蜂群，这些外来蜂群的混入在一定程度上增加了当地东方蜜蜂种群的遗传多样度。这种人为基因交流导致的种群遗传结构“复杂化”现象也不同程度地影响着其它地区的东方蜜蜂。

青藏高原东方蜜蜂不同地理种群间的共享单倍型较少，且这些共享单倍型最多只出现于3个地理种群(图2)，未发现广泛分布于9个地理种群的共享单倍型。表明该地区东方蜜蜂地理种群间的基因交流十分有限，各种群在近期经历了相对独立的繁衍过程。青藏高原是地球上海拔最高的高原，平均海拔4 000～5 000 m。高原上绝大部分区域并不适宜中蜂生存，只在海拔相对较低且有显花植物分布的区域生活着东方蜜蜂[16]。本试验的样本均采集自海拔2 200～2 700 m的区域，这些区域被连绵的高海拔山脉、高原所分隔。这在很大程度上限制了自然分蜂群的飞行范围，阻隔了不同地区间蜂群的自然基因交流。此类地形地貌有利于不同地理种群间遗传分化的加剧和不同亚种(形态类群)的形成，东方蜜蜂西藏亚种和阿坝亚种(形态类群)的形成正是“得益于此”。由此可见：山脉、高原阻隔是塑造青藏高原地区东方种群遗传结构的主要作用因素。

经与参考序列合并分析，单倍型系统发育关系(图3)显示：青藏高原地区与云南、广西及韩国、日本的东方蜜蜂系谱关系相对较近；与台湾岛、马来西亚、文莱、印度、泰国、俄罗斯等地区的系谱关系相对较远。东方蜜蜂的系统发育关系呈现显著的地理分布格局，这与我们之前研究[20]的结论是一致的。

本试验亦存在不足之处：(1) 在与参考序列合并分析时同源片段缩短，图3系统发育树是基于441 bp线粒体DNA片段构建的，一些变异信息可能因此而被“隐藏”；(2) 数据分析中中国大陆其它地区的参考序列相对缺乏，以至于对比分析不够全面。这些有待后续研究进一步加以完善！