黑龙江省张广才岭北部小飞鼠的遗传多样性与种群历史动态

田新民 廉明栋 宋雅祺 刘小慧 杨孟平 陈红

（牡丹江师范学院生命科学与技术学院，牡丹江 157011）

小飞鼠(Pteromys volans)隶属哺乳纲(Mammalia)啮齿目(Rodentia)鼯鼠科(Pteromyidae)飞鼠属，为树栖夜行滑行类物种。在国外主要分布于欧亚大陆北部，我国主要分布在东北、华北、西北等北方林区(蒋志刚等，2015)。不合理的森林采伐，生境的丧失与破碎化，已导致全球许多地区的小飞鼠种群处于濒危状态(Leeet al.,2008a;Nummertet al.,2020)。小飞鼠作为森林可持续经营的重要指示和伞护物种(Hurmeet al.,2007)，我国也将其列为“易危”与“三有”的重点保护野生动物(王志宝，2000；蒋志刚等，2016)。目前国外在小飞鼠巢址选择、家域、扩散、遗传多样性和分子系统进化等领域已开展系列研究(Hanskiet al.,2000;Selonen and Hanski,2004;Nakama and Yanagawa,2009;Lampilaet al.,2009;Yalkovskayaet al.,2015)，而国内仅在树洞巢穴保温机制方面有相关报道(陈亮，2018)，在生态学其他领域，特别是对这一物种的种群遗传学研究还存在空白。

种群遗传多样性和历史动态是衡量物种生存与利用的重要指标(Haiget al.,1990)，特别是遗传多样性的降低与近交衰退会导致物种灭绝，其遗传信息评价结果可为物种实施科学的保护与管理(Frankhamet al.,2010；张于光等，2019)。我国东北地区曾受森林采伐、盗猎、气候变化等因素的影响，使许多物种面临生存威胁，导致物种的近交、遗传多样性丧失、种群遗传分化等遗传信息发生改变(魏辅文等，2021；田新民，2021；Chenet al.,2022)。黑龙江省张广才岭北部为长白山脉与小兴安岭之间许多物种迁移的重要生态廊道区域，该区域的物种保护与生态修复备受关注(李维平等，2017；田新民，2021)。因此，本研究基于mtDNACyt b、控制区和nDNA的微卫星3种分子标记，分析小飞鼠的遗传特性及系统进化，从分子水平上揭示小飞鼠种群的动态变化及其遗传学机制，为该物种种群遗传资源的保护与管理奠定基础。

1 研究方法

1.1 研究地区概况

研究地区位于黑龙江省张广才岭北部的宾县林区 (北纬 45°30′～ 46°01′，东经 126°55′～ 128°19′)，地处松花江南岸(图1)。该地区属中温带大陆性季风气候，年平均气温3.9℃，年平均降雨量681 mm，无霜期130 d左右，属长白山植物区系。林区经营总面积1 110.81 km2，其中天然林面积675.76 km2，天然林中以阔叶林为主，针阔混交林和针叶林分布极少。人工林面积为272.72 km2，其中以落叶松林为主。分布的树种有紫椴(Tilia amurensis)、胡桃楸 (Juglans mandshurica)、水曲柳(Fraxinus mandshurica)、白桦 (Betula platyphylla)和蒙古栎(Quercus mongolica)等，这些树种为小飞鼠提供了洞巢保障。该区域内的紫貂(Martes zibellina)、黄喉貂(Martes flavigula)、黄鼬(Mustela sibirica)和豹猫(Prionailurus bengalensis)等食肉类为小飞鼠的主要天敌。

图1 本研究的采集样点分布图Fig.1 Sampling sites of this study

1.2 样品采集及DNA提取

2018年与2019年10—12月，对研究地区的小飞鼠毛发样本进行采集。在野外寻找小飞鼠的树洞巢，敲击树干发现个体后，采用网捕的方法在洞口处捕获个体。在个体尾部取20根带有毛囊的毛发样本，收集于封口袋内，标号并GPS定位。在捕获个体尾端剪掉部分毛发标记，避免个体重复采样，最后将个体放归野外。用同样的方法，两年期间，共采集29份毛发样本，-20℃保存备用。此外，在张广才岭中部的黑龙江省牡丹江市三道关林场采集3份毛发样本；大兴安岭北部内蒙古根河市区域的金河林业局采集1份毛发样本(陈亮、张子栋提供)，该4份样本仅用于小飞鼠系统进化关系分析。

DNA提取时，从毛发根部毛囊处取约0.5 cm长的一段，用75%的乙醇清洗两次以除菌，蒸馏水再清洗一次，弃废液后于鼓风干燥箱内55℃干燥。将处理后的20根毛囊样本，使用TIANamp Micro DNA Kit试剂盒(天根，北京)，并参考使用说明书对样本进行DNA提取。为进一步验证所用样本来自不同个体，基于后续微卫星分型数据，利用软件Excel mircosatellite tool kit(Park,2001)，寻找数据中相匹配的基因型。不同样本，所有位点的基因型相同或只有1个位点的1个等位基因存在差异，可判定为同一个体(Bellemainet al.,2005)。

1.3 PCR扩增

根据已发表小飞鼠mtDNA全序列(Ryuet al.,2013;Limet al.,2018)，利用Primer Premier v 5.0(Singhet al.,1998)设计引物用于mtDNACyt b和控制区全序列的扩增。Cyt b扩增引物，上游：5′-ACA TGG AAT CTAACC ATG ACC AA-3′，下游：5′-AGA CTT CAT TGT TGG TTT ACA AGA C-3′。控制区扩增引物，上游：5′-GTT CCA CCT TCA ACT CCC AAA-3′，下游：5′-CAT TTT CAG TGC TTT GCT TTG AT-3′。

根据GenBank中报道的微卫星位点，结合Painter等 (2004)、Kiesow 等 (2011)、Zittlau 等(2000)、Jumpa等 (2015)、Hale等 (2001)和 Todd(2000)等文献，选取鼯鼠科(Pteromyidae)小飞鼠指名亚种(P.v.volans)、北美飞鼠(Glaucomys sabrinus)、海南小飞鼠(Hylopetes phayrei)和赤颊林飞鼠 (Hylopetes sagitta)，松鼠科 (Sciuridae)北松鼠(Sciurus vulgaris)5个物种中扩增效果较好的34对引物，以本研究中的小飞鼠DNA为模板进行扩增，筛选微卫星位点。最终获得并使用扩增稳定、具多态性的12个位点，其上游引物5′端进行荧光标记 (Fam、Hex、Tamra或Rox)(表1)。

表1 本研究选取的12对微卫星引物信息Table 1 Details of 12 microsatellite loci used for the study

三种分子标记的PCR反应体系均为20 μL：5 U/μL Ex Taq DNA polymerase(TaKaRa， Japan)0.1 μL， 10 × Buffer 2 μL， 2.5 mmol/L dNTPs 1.6 μL，10 μmol/L 上、下游引物各 0.5 μL，模板 DNA 1.5 μL，PCR grade water(天根，北京)13.8 μL。PCR 反应条件均为：94℃，2 min；(94℃，30 s；Ta，30 s；72℃，30 s)× 35 cycles；72℃，10 min。其中，Cyt b和控制区引物的退火温度(Ta)均为55℃，微卫星引物的退火温度见表1。对于微卫星位点，PCR扩增选择3次重复，其中2次及以上重复均只获得1条PCR产物带的判定为纯合子，3次重复均为2条PCR产物带的判定为杂合子(张于光等，2019)。为了监测污染，在扩增的同时均附加一个不含DNA的阴性对照。所有的引物合成、双向测序和基因分型均由上海生工完成。

1.4 数据分析

mtDNACyt b和控制区序列，首先利用Clustal X 2.1(Larkinet al.,2007)进行序列比对。然后用DnaSP 5.10(Librado and Rozas,2009)计算种群的变异位点数(S)、单倍型数量(H)、单倍型多样性(Hd)和核苷酸多样性(Pi)，并进行错配分布(Mismatch distributions)分析，以通过单倍型差异分布来估测种群历史扩张事件。通过Tajima’s D和Fu’s Fs中性检验来判断种群是否显著偏离中性突变，进一步评价种群的历史扩张事件或瓶颈效应。利用Network 4.6(Bandeltet al.,1999)，基于Median-joining算法，构建小飞鼠的单倍型网络图，分析单倍型间的进化关系。

为探讨小飞鼠系统进化关系，从GenBank下载小飞鼠Cyt b序列 (1 068 bp)共51条 (表2)，合计本研究的18个单倍型，其中宾县样本14个单倍型(代码：CHB1～CHB14)、三道关样本3个单倍型(代码：CHS1～CHS3)和根河样本1个单倍型(代码：CHG)，共计69条序列。利用BEAST 1.6.1(Drummondet al.,2012)，基于贝叶斯法(Bayesian inference,BI)构建系统发育树；并通过jModelTest 2.1.10(Darribaet al.,2012)检验选择最优替换模型及相关参数(HKY+G)。以日本小飞鼠(Pteromys momonga)(单倍型：PM01和PM02，GenBank登录号：AB164675和AB164676)作为外群，以获取发育树的根。构建BI树时，以MCMC运行3 000 000代，每100代抽样1次，舍弃前25%的老化值，并用贝叶斯后验概率评估各节点的支持率。仅显示50%以上的支持率，并对发育树进行美化。

表2 GenBank下载的小飞鼠Cyt b基因信息Table 2 Mitochondrial cytochrome b gene of Siberian flying squirrel from GenBank

基于微卫星分型数据，首先利用Microchecker 2.2.3(van Oosterhoutet al.,2004)检测各位点是否存在无效等位基因(Null allele)，判断位点是否存在严重的分型错误及能够满足遗传学分析要求。然后利用GenAlEx 6.0(Peakall and Smouse,2006)计算种群的等位基因数(Na)、有效等位基因数(Ne)、观测杂合度(Ho)、期望杂合度(He)和近交系数(Fis)。用Excel mircosatellite tool kit(Park,2001)计算多态信息含量(PIC)，PIC＞0.5的位点判定为高多态性位点(Botsteinet al.,1980；田新民，2021)。再利用Genepop 4.0(Raymond and Rousset,1995)检测种群和各位点是否符合Hardy-Weinberg平衡，以及位点间的连锁不平衡，概率检验采用马可夫链法(Markov chain method)，参数为10 000 dememorization、20 batch和5 000 iteration。最后，两个检验以Bonferroni法对显著性进行修正。

此外，基于微卫星数据，利用Bottleneck 1.2(Piryet al.,1999)进一步检测小飞鼠种群是否经历过瓶颈效应，选择两阶变异模型(Two-phase model,TPM)和逐步突变模型(Stepwise mutation model,SMM)进行Wilcoxon检验，重复1 000次，其中TPM被认为是微卫星分析的最适合模型(Ellegren,2004)。TPM选择79%变异遵从SMM，变异系数为 9%(Piryet al.,1999)。利用 Structure 2.3.4(Pritchardet al.,2000)进行种群遗传结构分析，参数为Length of Burnin Period 100 000，Number of MCMC Reps after Burnin 10 000，K=1～6，Number of Iterations 20。运算结果上传Structure Harvester Web v0.6.94(Earl and vonHoldt,2012)进一步分析，以Ln Pr(X|K)曲线最大峰值来判断种群遗传结构的最佳分组数K。

2 结果

2.1 mtDNA序列分析

29个mtDNACyt b全序列(1 140 bp)中共检测到变异位点34个，其中单变异位点11个，简约信息位点23个。碱基含量为A 27.9%、T 31.0%、C 28.6%、G 12.4%，A+T含量(58.9%)高于C+G含量(41.0%)。碱基变异均为转换和颠换，无缺失或插入，转换(颠换)平均值为16.00。在29个样本中得到26个控制区全序列(1 066 bp)，检测到72个变异位点，其中单变异位点24个，简约信息位点47个。碱基含量为A 31.9%、T 30.9%、C 26.0%、G 11.2%，A+T含量 (62.8%)高于C+G含量(37.2%)。碱基变异多为转换和颠换，有1个缺失或插入，转换(颠换)平均值为5.25。两种序列，均存在明显的AT偏倚性，且G含量最低(图2，3)。

2.2 遗传多样性评价

对于mtDNA数据，在29条Cyt b序列中共检测到14个单倍型，其中单倍型多样性指数(Hd)为0.909±0.031，核苷酸多样性指数 (Pi)为(0.616±0.065)%。26条控制区序列中存在15个单倍型，其中Hd为0.945±0.024，Pi为(1.698±0.068)%(表3)。两种序列中仅有1只个体的稀有单倍型均占多数，Cyt b和控制区都为9个，分别占各自所有单倍型的64.29%和60.00%(图2～4)。

图2 小飞鼠种群14个Cyt b单倍型变异位点.N：具有相同单倍型的个体数量Fig.2 The variable sites of 14 cytochrome b haplotypes identified in Siberian flying squirrel population.N:Number of individuals with each haplotype

图3 小飞鼠种群15个控制区单倍型变异位点.N：具有相同单倍型的个体数量Fig.3 The variable sites of 15 control region haplotypes identified in Siberian flying squirrel population. N: Number of individuals with each haplotype

图4 基于Cytb(A)与控制区(B)序列构建的单倍型网络图.圆形面积代表单倍型频次，空心点代表缺失的单倍型Fig.4 Median-joining network of haplotypes based on Cytb(A)and control region(B).The size of each circle represents the frequency of haplotypes,and small hollow dot represents missing haplotype

表3 基于Cyt b与控制区序列的种群遗传多样性分析Table 3 Analysis of population genetic diversity based on cytochrome b(Cyt b)and control region(CR)

基于微卫星分型数据，证实29份毛发样本DNA均来自不同个体。12个微卫星位点经过Microchecker检测，未发现有无效等位基因，位点间均不存在连锁不平衡。29只个体均获得12个位点的准确基因型，其中种群平均等位基因数(Na)为(13.167±2.725)个，有效等位基因数 (Ne)为(7.714±1.887)个，各位点的有效等位基因数显著少于实际等位基因数(P＜0.01)。12个微卫星位点的PIC为0.806～0.908，皆为高多态性位点(PIC＞0.5)。遗传多样性显示，种群的观测杂合度 (Ho)为0.727±0.175，期望杂合度 (He)为0.864 ± 0.033(表4)。

表4 基于12个微卫星位点的种群遗传多样性分析Table 4 Analysis of population genetic diversity based on 12 microsatellite loci

2.3 种群历史动态分析

mtDNACyt b基因和控制区序列构建的错配分布图呈明显的多峰曲线状态(图5)。基于中性检验，Tajima’s D和Fu’s Fs的结果也均未显著偏离零(P＞0.05)，种群没有显著偏离中性突变(表3)。两种分析方法结果表明，小飞鼠种群经历了比较复杂的种群历史，没有发生过历史扩张事件。微卫星的种群瓶颈效应检测显示，TPM和SMM模型中均未发现显著的杂合度过剩(P＞0.05)，同时等位基因频率没有显著偏离正态L-shaped分布。再结合Cyt b和控制区的中性检验结果，均未检测到小飞鼠种群近期的遗传瓶颈效应。

图5 基于Cyt b(A)和控制区(B)序列的小飞鼠种群错配分布图Fig.5 Mismatch distribution for Siberian flying squirrel population based on Cyt b(A)and control region(B)

微卫星数据显示，种群观测杂合度显著小于期望杂合度(P＜0.05)，呈现一定程度的杂合度不足，种群的近交系数0.159，为极显著偏离零的正值(P＜0.01)(表4)，表明小飞鼠种群存在近亲繁殖情况。根据STRUCTURE软件的种群遗传结构表明，当K=1时，Ln Pr(X|K)具有最大平均值，且重复计算的波动性最小(图6)，说明小飞鼠种群内不存在显著的遗传分化。

图6 基于微卫星数据STRUCTURE聚类结果的Ln Pr(X|K)变化趋势Fig.6 Changing trends of Ln Pr(X|K)from STRUCTURE clustering results for the microsatellite

2.4 系统发育关系分析

基于69个mtDNACyt b单倍型构建的贝叶斯系统发育树显示，小飞鼠存在3个明显的遗传谱系，包括远东、欧亚大陆北部和日本北海道3个分支，本研究张广才岭和大兴安岭的样本单倍型归属为远东谱系。俄罗斯东部、韩国与中国样本单倍型共同形成远东谱系，3个地区之间无共享单倍型，并且俄罗斯、韩国与中国单倍型分别存在两个独立亚分支。俄罗斯中西伯利亚(欧亚大陆东北部)与俄罗斯西部(欧亚大陆西北部)的两个亚分支构成了欧亚大陆北部谱系，其中蒙古单倍型归属于欧亚大陆东北部亚分支。日本北海道地区的单倍型独立形成北海道谱系(表2，图7)。

图7 基于Cyt b构建的小飞鼠贝叶斯系统发育树Fig.7 Phylogenetic tree of Siberian flying squirrel obtained for cytochrome b gene in Bayesian inference

3 讨论

遗传多样性的丧失会导致动物对环境变化的适应能力下降，甚至物种的灭绝，评价种群的遗传多样性，是物种保护的重要内容(Frankhamet al.,2010)。单倍型多样性(Hd)和核苷酸多样性(Pi)是衡量种群mtDNA遗传变异程度的重要指标，其中Pi值考虑了单倍型在种群中所占的比例，在评价遗传多样性时更为精确(Neigel and Avise,1993；马英等，2019)。本研究控制区的Pi值(1.695%)明显高于Cyt b(0.616%)，Nummert等(2020)对爱沙尼亚、芬兰和俄罗斯的小飞鼠种群研究中也呈现此规律。mtDNA非编码控制区的进化速度较快、承受的选择压力较小 (Krojerová-Prokešováet al.,2013；佐日古丽·伊斯马伊力等，2019)，在许多动植物和人类中也普遍存在基因非编码区多态性高于编码区的现象(Crochet and Desmarais,2000)。微卫星数据显示，种群平均观测杂合度 (Ho)和期望杂合度 (He)分别为0.727和0.864。本研究中Pi值与Lee等(2008a)关于黑龙江省小飞鼠种群(Pi=0.614%)和韩国种群(Pi=0.616%)的研究基本一致，仅低于俄罗斯滨海边疆区(Pi=0.950%)；远高于爱沙尼亚、芬兰、俄罗斯大部分地区和日本北海道地区。与其他地区相比，本研究的遗传多样性处于最高水平(表5)，同时参考Grant和Bowen(1998)、Yuasa等 (2007)种群高单倍型多样性(Hd≥0.5)和高核苷酸多样性(Pi≥0.5%)的标准，认为张广才岭北部小飞鼠种群的遗传多样性丰富。但是研究发现，种群中Cyt b和控制区的稀有单倍型均在60%以上，微卫星各位点的有效等位基因数显著少于实际等位基因数(P＜0.01)，说明小飞鼠的基因频率分布处于一个不稳定的状态，未来种群发展中存在等位基因丢失的风险。因此，应加强张广才岭北部小飞鼠种群的保护与管理，防止种群遗传多样性的急剧下降。

表5 基于mtDNA和微卫星位点的鼯鼠科主要物种种群遗传多样性参数比较Table 5 Comparison of population genetic diversity at mtDNA and microsatellite loci in the Pteromyidae

景观格局和物种的迁移能力是影响种群遗传结构的最重要因素(Yuasaet al.,2007;Pérez-Esponaet al.,2008)。本研究的宾县地区有高速G1011、国道G221和哈佳铁路东西方向穿过，被分成南、北两块隔离区域，采样的区域主要集中在南部，该林区无明显的景观阻隔。同时小飞鼠作为滑行类物种，具有一定的迁移能力(Selonen and Hanski,2004,2012)，保证了个体间的基因交流，因此种群内没有出现显著的遗传分化。本研究的张广才岭北部区域、尚志和方正林区受多条高速、高铁和国道，以及延寿农田区等因素的阻隔，其森林生境呈现破碎化分布，未来待加大取样范围，进一步评价北部区域小飞鼠种群的遗传结构。基于mtDNA和微卫星，本研究均未检测到小飞鼠种群的近期遗传瓶颈效应。由于小飞鼠种群历史资料的匮乏，缺少数量变化的准确信息。鉴于之前的大量森林采伐，据主管部门报道的多次滥捕行为，本研究地区小飞鼠种群数量可能较历史时期下降，而数量下降的程度还未能表现遗传上的种群瓶颈效应。受森林采伐导致的生境破碎化，爱沙尼亚和芬兰地区的小飞鼠种群数量一直呈下降趋势，分别成为当地的濒危与易危物种，也未检测到遗传瓶颈效应(Nummertet al.,2020)。

在Hardy-Weinberg平衡检验中发现，本研究共有6个微卫星位点以及种群整体上偏离平衡。Hardy-Weinberg平衡是建立在自然种群随机交配假设基础上的，认为在濒危物种种群中出现不符合平衡的现象，主要原因有无效等位基因、种群亚结构和近亲繁殖导致的杂合度不足(Ladeet al.,1996；任鹏等，2017)。本研究未检测到无效等位基因和种群亚结构，因此排除这两个因素导致其偏离Hardy-Weinberg平衡的可能。Fis是评价种群是否存在近亲繁殖的一个重要参数，当Fis为正值时表示群体内存在近交(Weir and Cockerham,1984；任鹏等，2017)。本研究6个位点Fis均为正值，种群的Fis为0.159，也为显著偏离零的正值，因此认为种群内的近亲繁殖是导致微卫星位点偏离Hardy-Weinberg平衡的主要原因。小飞鼠种群数量的下降，与尚志和方正等林区的隔离，可能是影响种群内个体近亲繁殖的重要因素。同样在爱沙尼亚，因小飞鼠种群的生境破碎化和数量急剧下降也表现出了近亲繁殖(Nummertet al.,2020)。本研究较高的种群遗传多样性表明，小飞鼠在其近期历史中可能具有较大的有效种群规模(Balakrishnanet al.,2003)，而种群内高比例的稀有单倍型和等位基因为个体近亲繁殖的表现。近亲繁殖可导致物种遗传多样性的降低与丧失，对环境的适应性下降，最终可能面临濒危或灭绝的威胁(Lampilaet al.,2009;Nummertet al.,2020)。基于上述分析认为，亟待加强张广才岭北部小飞鼠种群的保护，连通破碎化生境增加种群间的基因交流，防止遗传多样性的下降与丧失。

本研究证实小飞鼠存在3个明显的进化分支，远东、欧亚大陆北部和日本北海道，结果与其他研究基本一致(Oshidaet al.,2005;Leeet al.,2008a;Nummertet al.,2020)。其中，Lee等(2008a)和Nummert等(2020)研究中仅包括中国黑龙江省5个样本得到的5个单倍型，其中3个单倍型(CHJ3、CHJ4和CHJ5)分别与本研究宾县单倍型(CHB12、CHB13和CHB4)为共享单倍型。在远东谱系中，韩国与俄罗斯东部单倍型分处两个亚分支内，两个地区之间存在明显的遗传分化(Leeet al.,2008a)，而中国样本单倍型在两个亚分支中都有分布，说明本研究地区的小飞鼠与这两个区域之间存在一定的基因交流。与相似分布的花鼠(Tamias sibiricus)和北松鼠不同，在mtDNA研究中，韩国、俄罗斯东部和中国东北3个地区的花鼠存在明显的遗传分化(Leeet al.,2008b)，而北松鼠却不存在(Leeet al.,2008a)，说明3个物种存在不同的进化历史。俄罗斯中西伯利亚(欧亚大陆东北部)与俄罗斯西部(欧亚大陆西北部)的两个亚分支共同形成了欧亚大陆北部谱系，两个区域没有显示出更明显的系统地理结构，说明乌拉尔山脉(Ural Mountains)似乎不是小飞鼠明显的遗传屏障。北海道种群形成最为明显的一个独立分支，研究认为北海道种群的系统地理分化可能是荷斯廷(Holsteinian)间冰期的结果，在欧亚种群内部分化之前，该种群已经与欧亚种群分离(Oshidaet al.,2005)。本研究中大兴安岭地区仅有1个根河样本，并且单倍型与其中宾县单倍型共享，未来应增加取样量和研究区域，进一步探讨我国东北地区及全国范围内种群的遗传结构和亚种分类问题，为小飞鼠的科学保护与管理提供参考。