东北虎豹国家公园梅花鹿的空间遗传格局及其影响因素

段兴汉,吴 峰,张素青,鲍 蕾,王红芳,*

1 生物多样性与生态工程教育部重点实验室,北京 100875 2 东北虎豹国家公园保护生态学国家林草局重点实验室,北京 100875 3 北京师范大学生命科学学院,北京 100875

空间遗传格局是物种在种群内或种群间遗传多样性的空间分布格局,种群遗传学过程是空间遗传格局形成的基础,而景观中斑块、基质的空间配置、物种性状差异可能增加了空间遗传格局的复杂性[1]。然而,在不同物种或者不同生境下,这一期望格局的强度却存在较大变异,主要原因是景观基质可能影响种群间的功能连通度。多项研究表明,地理距离、地形、人造景观、土壤类型、温度、湿度等都对动植物分布有显著的影响[2—5]。因此,探究物种的空间遗传格局,分析格局形成的过程,有效识别可能的基因流障碍,可为制定有效种群保护策略提供依据。

东北梅花鹿(Cervusnipponhortuloru)是我国仅存梅花鹿三个亚种之一,种群数量相对较大。主要分布于我国吉林省、黑龙江以及俄罗斯滨海边疆区连片的广大地区。然而东北梅花鹿的空间分布十分不均匀。在我国,东北梅花鹿种群密度最高的地区为吉林省珲春市境内的中俄边境线附近,而在其他大部分地区东北梅花鹿种群密度低[6]。梅花鹿作为食物链的初级消费者,该物种对植物群落的结构与功能具有很强的塑造作用[7]。此外,梅花鹿也是该地区旗舰物种东北虎(Pantheratigrisaltaica)、远东豹(Pantherapardusorientalis)主要猎物,大型有蹄类动物的可获得性与捕食者密度增加直接相关[8],因此梅花鹿对虎豹种群的分布与密度具有一定的限制作用。

本研究基于非损伤样品对虎豹国家公园的梅花鹿进行遗传多样性和空间遗传格局分析,以期了解该种群现有遗传多样性水平以及影响其空间遗传格局的环境因素和人为因素。

1 材料与方法

1.1 样品采集与处理

本研究于2018年8月至2020年1月多次在东北虎豹国家公园东部的中俄边境附近采集梅花鹿粪便样品110份(图1)。为了满足本研究需要,本研究采取嵌套式采样策略[9],即在大、中、小不同距离等级下确保采集到的一定数量的样品。

图1 梅花鹿粪便样品分布图Fig.1 Map of the sika deer fecal samples

粪便样品采集时佩戴一次性聚乙烯手套(PE手套)或使用一次性筷子,为防止交叉污染,采样工具用后即换。采集的样品原则上仅采集新鲜的粪便。每堆粪便挑取5—8粒,装入50ml离心管中,记录全球定位系统(GPS)坐标,避光处理;为了更好保存粪便样品中的DNA,需要对样品尽快干燥处理:(1)倒入适量95%无水乙醇,液面高度以没过管内样品为宜;(2)浸泡24h后将无水乙醇倒出,加入足量变色硅胶,加入无水乙醇或变色硅胶后,样品可以常温保存,否则应当低温或冷冻保存;(3)每隔一周检查离心管内硅胶的颜色,应在硅胶完全变红之前及时更换,多数样品在更换1—3次即可保持硅胶不再变色。

样品质量信息的现场判断较为重要,因为样品越新鲜,DNA降解越少,质量较差的样品往往导致扩增失败。在野外可根据如下标准判断梅花鹿粪便样品的新鲜程度:

a)分散程度。粪球不成堆分布的,形态大小明显不同的几堆粪便混杂在一起的,应当舍弃。

b)表面物理特征。新鲜粪便表面一般呈现黄色或深褐色,粪球形状饱满,表面光亮、黏滑,没有裂纹,不长霉。粪便在野外暴露一定时间后表面一般会变黑、变硬,表面出现裂纹,裂纹越多说明在野外暴露时间越久。

c)软硬程度。一般来说粪便越软越新鲜。

d)被枯枝落叶、积雪等覆盖的一般不新鲜,可根据降雪时间判断粪便质量。

e)梅花鹿卧迹内或1m范围内、拨雪取食痕迹附近一般可寻获粪便,因此可根据卧痕、咬痕、足迹等新鲜程度判断粪便新鲜程度。

f)冬季采集的粪便一般都能够满足样品DNA提取和分型的需要,但是冬季采集时要注意不要采集表面开裂、位于落叶层以下或与土壤黏连的粪便。

1.2 粪便基因组提取、个体鉴定

本研究使用QIAamp Fast DNA Stool Mini Kit (QIAamp)进行粪便基因组提取。为确认粪便样品所属物种,选用Riaz等[10]开发的鹿科动物通用引物12srRNA对粪便DNA进行扩增。引物序列如下,F: TAG AAC AGG CTC CTC TAG;R: TTA GAT ACC CCA CTA TGC。PCR反应体系为40μL:20μL 2×EasyTaq PCR SuperMix,0.8μL BSA,引物各0.8μL(10μMol),4μLDNA模板,加13.6μL超纯水补足。聚合酶链式反应(PCR)程序为:94℃预变性5min;94℃变性30s,51℃退火30s,72℃延伸30s,35个循环;72℃延伸10min。本研究测序和引物合成由北京擎科新业生物技术有限公司完成。测序结果质量使用Codoncode Aligner 6.0.2进行检测。利用美国国立生物技术信息中心(NCBI)数据库进行数据比对(https://www.ncbi.nlm.nih.gov/),设置匹配度100%,确定粪便样品所属物种。

1.3 微卫星的筛选和分型

本研究根据过去使用率较高的有蹄类动物微卫星序列,筛选出9对稳定扩增且多态性较高的引物(表2)。扩增体系:Takara Taq HS 0.5μL,10X PCR Buffer (Mg2+Plus) 2μL,dNTP Mixture 2μL,引物各0.5μL,BSA 1μL,DNA模板 3μL,加超纯水至20μL。扩增程序:95℃预变性5min;95℃变性30s,退火30s,72℃延伸1min,进行35个循环;72℃终延伸10min。PCR扩增产物使用ABI 3730测序仪测序。微卫星数据的读取和校正在Genemapper 4.0中完成。

表1 各种变量对梅花鹿种群遗传格局影响的假设Table 1 Hypothesis of the influence of different variables on genetic pattern of sika deer population

表2 野生东北梅花鹿种群9个微卫星位点分析Table 2 The analysis of 9 microsatellites locus in wild sika deer population

由于粪便DNA质量差,等位基因容易丢失。为了避免这种现象,本研究采用多管扩增法(multi-tube PCR)进行个体基因型确定[11],基本方法如下:杂合子由至少两次相同的分型结果确定,而纯合子则由3次相同分型结果确定。每份样品每个位点先进行4次重复扩增,判断基因型结果;对仍然无法确定分型情况的样品再进行3次独立扩增,判断基因型。此后仍无法确定基因型的样品,则将该个体的该位点处理为缺失数据。

1.4 数据分析

1.4.1遗传多样性

利用GIMLET 1.3.1[12]计算微卫星位点的鉴别能力:同一性概率(PID)和同怀群同一性概率(PIDsibling);使用MICRO-CHECKER[13]检查各微卫星位点是否存在无效等位基因、等位基因丢失的现象。使用CERVUS和FSTAT 2.9.3.2[14]计算位点等位基因数量(Na)、有效等位基因数量(Ne)、观测杂合度(HE)、期望杂合度(HO)、多态信息含量(PIC)等种群遗传学指标;使用Genepop 3.4[15]检查每个位点的哈温平衡和连锁平衡。

使用CERVUS[16]进行个体鉴别。

1.4.2聚类分析

使用STRUCTURE 2.3.4[17]软件对鉴定的所有梅花鹿个体遗传结构进行分析。计算过程采用混合祖先模型和等位基因关联模型,开始时的预迭代次数为105,马尔可夫链蒙特卡罗迭代(MCMC)重复次数为104次,K值范围定为K=1—6,对每个K进行20次重复运算。利用Structure Harvester最大似然值以及ΔK随K值的变化,综合判断梅花鹿种群的最优分组。

使用GenAlEx 6.5[18]对鉴定的所有梅花鹿个体做主坐标分析。

1.4.3空间自相关分析

有限的扩散能力常常导致种群在遗传距离上具有显著的空间自相关模式[19]。本研究首先使用MSA 4.05[20]软件计算个体两两之间的遗传距离,距离指标为共有等位基因比例。个体间的地理距离,以及空间自相关检验都在GenAIEx中进行。空间自相关检验的置换数设置为999次。

1.4.4空间变量对遗传分化影响分析

本研究主要分析了海拔、坡度、坡向、地表起伏率(Surface Relief Ratio, SRR)、人类干扰5个变量对梅花鹿种群遗传结构的影响。这5个变量多被认为与大中型哺乳动物扩散阻碍相关[21]。前四个变量来自于分辨率30m×30m的DEM 数据(http://www.usgs.gov/)。人类干扰使用道路、居民点图层进行模拟。

各种环境变量的阻力模型通过ArcGIS 10.1中的“空间分析工具”下的“重分类工具”建立。

(1)海拔:随海拔升高温度降低,因此海拔可能影响梅花鹿的食物可获得性,对梅花鹿来说可能存在一个最适海拔。本研究利用逆高斯方程模拟[22]海拔高度对扩散的影响:

海拔高度分为9组,每组间距100m。Rmax为最大阻力值,分别设为2、10、100、500、1000。Eopt为梅花鹿扩散最适海拔高度,分别设为200、300、400、500m。ESD分别设为50、100、200m。基于以上3个参数,建立60个海拔高度阻力模型。

(2)坡度:从运动时的能量消耗角度推测,梅花鹿可能更倾向于利用较缓坡度。本研究用逆高斯方程模拟坡度对梅花鹿扩散的阻碍作用:

坡度分为9组,每组间距7°。Rmax为最大阻力值,分别设为2、10、100、500、1000,θopt为梅花鹿扩散最适坡度,分别设为5°、10°、15°、20°、30°,θSD分别设为2°、5°、10°。基于以上3个参数,建立75个坡度阻力模型。

(3)坡向:坡向不同,温度、湿度等微生境不同,从而可能影响梅花鹿的生境选择和扩散通道的选择。用热负荷指标函数模拟坡向对梅花鹿的扩散阻碍[22]:

坡向分为9组,8组为坡向组,一组为平面组(-1°)。坡向组每组范围为45°,平面组赋值为Rmax/2,其他按上述函数进行赋值。Rmax分别设为2、10、100、500、1000,θopt分别设为0°、45°、90°、135°、180°、225°、270°、315°。x设为0.5、2、10。基于以上3个参数,建立坡向阻力模型120个。

(4)地表起伏率:地表复杂程度将影响梅花鹿运动时的能量消耗,地形复杂程度越高,梅花鹿扩散时的能量消耗越大。SRR对动物扩散的阻力由以下函数进行模拟[22]:

R=(SRRx)×Rmax+1

建立SRR阻力模型之前,首先下载Geomorphometry and Gradient Metrics 2.0工具(https://github.com/jeffreyevans/GradientMetrics),并使用ArcCatalog加载至工具箱中,然后打开ArcMap,使用该工具中的Surface Relief Ratio工具对研究区域的DEM进行重分类:分析窗口设置为“Circle”,半径分别设置为2、10、25、50,分析单元为“Cell”。再对4个不同半径的SRR栅格图层分别进行重分类,SRR值共分为10组,每组组距为0.1。Rmax分别设为2、10、100、500、1000,x分别为0.5、2、10。建立SRR阻力模型60个。

(5)建立人类干扰阻力模型6个。根据居民点、道路分布情况,首先在ArcMap中使用BufferWizard工具,对道路、村、镇矢量分别设置1km、1.5km、3km缓冲区;然后使用Union工具,将3个缓冲区图层合并,在属性表中新建一列,使用“Field calculator”将该列赋值为10,同时用同样的方法将研究区矢量图的该列赋值为1;再次使用Union工具将两个图层合并。最后使用“Circuitscape for Arcgis”工具,将矢量转换为栅格格式。使用重分类功能,将不包含居民点和道路的栅格的阻力值设置为1,其他栅格的阻力值R分别设为2、10、100、250、500、1000。

(6)建立空间距离阻力模型1个。对DEM栅格进行重分类,将所有栅格赋值为1,以便后续分析中排除空间自相关对遗传距离的影响。

(7)利用生境适宜性模型建立15个阻力模型。首先使用MaxEnt 3.4.1[23]制作生境适宜性图层,选取的变量为海拔、坡度、坡向、地表起伏率等4个环境变量以及与居民点距离、与省道距离2个人类干扰变量。这些变量数据来源及处理方法如上所述。考虑到变量间可能存在相关性影响模型精度,先对所有变量因子做相关性分析,去除相关系数绝对值大于0.75的变量。共计379个梅花鹿出现点数据来源于梅花鹿粪便样品坐标以及本实验室近几年红外相机数据。本研究选择线性、二次方、乘积和多拐点阈值关系等4种运算方法的综合运用;由于测试样本比例对模型预测精度有一定影响,所以本研究首先以不同的测试样本比例进行模拟,每个比例重复10次,筛选出平均AUC最高(25%,0.861)的一个生境适宜性图层。

随后依据该模型建立15个阻力模型。由于该图层取值范围为H=0—1,值越大代表生境适宜性越高,所以根据以下方式进行重分类:

R=(1-H)x×Rmax+1

生境适宜性值共分为10组,每组组距为0.1。Rmax分别设为2、10、100、500、1000,x分别为0.5、2、10。

阻力模型建立完成后,使用Circuitscape for ArcGIS工具生成各变量的不同模型的阻力矩阵。计算模式为成对的模式、连接附近8个栅格。

使用R包ecodist[24]进行偏曼特尔检验,在排除空间自相关(IBD)的影响后,分析遗传距离(GD)和景观阻力(LR)的相关性。可表示为GD-LR|IBD,其中GD为遗传距离,LR为景观变量阻力,IBD为空间自相关。对任一景观变量模型进行检验时,当GD—LR|IBD的结果显著,而GD—IBD|LR的结果不显著时,才认定该模型能够解释个体间的遗传变异模式[21]。

2 结果与分析

2.1 物种鉴定和微卫星分型

使用12srRNA对110份粪便样品进行扩增,其中101份样品测序成功。通过基于局部比对算法的搜索工具(BLAST)分析,除了7份样品为狍(Capreoluscapreolus)粪便外,其他94份样品均为梅花鹿粪便样品。利用Gimlet软件计算发现,这9对微卫星位点的联合识别分辨率较高。PIDunbiased=7.892×10-10,PIDbiased=1.468×10-9,PIDsibling=3.663×10-4;通过Cervus软件进行个体鉴别,在成功扩增的87份样品中,鉴别出86个个体。在所使用的9个微卫星位点中,所有位点都符合遗传连锁平衡,RT1、RM188、BM6506显著的偏离哈温平衡,且可能存在无效等位基因。由于这三个位点在梅花鹿种群遗传多样性研究中广泛使用,且RT1和RM188具有较高的等位基因数量和遗传多样性,为个体鉴定提供了更精准的信息,本研究保留这三个位点。

2.2 种群遗传多样性

本研究所使用的9对微卫星的等位基因数量在4—11之间,共扩增出68个等位基因(表2)。平均每个位点有7.56个等位基因,其中RT1等位基因最多,有11个等位基因,BM6506位点等位基因最少,有4个等位基因。这9个梅花鹿微卫星位点平均多态信息含量(PIC)为0.683,除BM6506(PIC=0.410)之外,其余位点的多态信息含量都在0.631—0.798之间,由于PIC>0.5为高度多态性,0.25≤PIC≤0.5为中度多态性,因此,BM6506为中度多态性位点,其余均为高度多态性位点。观测杂合度位于0.346—0.886之间,平均观测杂合度为0.689;期望杂合度位于0.472—0.842之间,平均期望杂合度为0.721。期望杂合度普遍高于观测杂合度。

2.3 种群遗传结构

STRUCTURE得出聚类数量估计值。L(K)峰值出现在K=1(图2),因此该地区梅花鹿最有可能的分群数量为K=1,当K>1时,种群中的个体被分到不同群的比例大致相当,说明本研究的目标种群栖息的景观中不存在显著的屏障效应。

图2 采用贝叶斯聚类分析方法得到的聚类数量估计值 Fig.2 Estimated number of clusters inferred using Bayesian clustering analysisL(K)的峰值出现在K=1处

主坐标分析(PCoA)表明(图3),所有样品点未形成明显的聚集分布。结合STRUCTURE和PCoA结果,推测该地区没有明显的遗传分组,该种群处于哈温平衡。

图3 梅花鹿种群的主坐标分析Fig.3 Principal Coordinates (PCoA) Analysis of the sika deer population图中各点为本研究的全部样品

Ta:退火温度 Annealing temperature;N:样本量 Sample Size;Na:等位基因数量 No. Alleles;Ne:有效等位基因数量 No. Effective Alleles;PIC:多态信息含量 Polymorphic Information Content;HO:观测杂合度 Observed Heterozygosity;HE为期望杂合度 Expected Heterozygosity;*代表该位点偏离哈温平衡并且存在无效等位基因

2.4 空间自相关分析

对该地区东北梅花鹿种群85km范围内的空间自相关分析表明,该地区的梅花鹿出现一定的空间结构。在第一个(0—5km,r=0.017,P=0.003)、第二个(5—10km,r=-0.009,P=0.043)空间等级表现为显著的空间自相关系数。在0—85km内的所有距离等级内,r值有正有负,呈现明显的波动格局(图4)。

图4 梅花鹿种群空间自相关分析Fig.4 Spatial autocorrelation analysis of sika deer population

对该地区东北梅花鹿种群10km范围内的空间自相关分析表明,该地区的梅花鹿出现一定的空间结构。在第1个(0—1km,r=0.043,P=0.002)、第9个(8—9km,r=-0.030,P=0.039)、第10个(9—10km,r=-0.025,P=0.016)空间等级表现为显著的空间自相关系数。在0—10km内的所有距离等级内,r值有正有负,呈现明显的波动格局(图4)。

2.5 空间变量对遗传分化影响分析

利用ArcGIS建立共计336个阻力模型,其中海拔阻力模型60个,坡度阻力模型75个,坡向阻力模型120个,地表起伏率阻力模型60个,人类干扰阻力模型6个,生境适宜性阻力模型15个。

在所有336个模型中,依据海拔、坡向、地表起伏率、人类干扰假设建立的246个阻力模型与遗传距离之间的关系并不显著(P>0.05)。在依据坡度假设建构的75个阻力模型中,只有1个模型与遗传距离有显著的正相关关系(P<0.05,表3),该模型同时也是在控制空间自相关影响后,在所有模型中与遗传距离相关性最高(r=0.088)的模型。综合所有变量的15个生境适宜性模型阻力模型与遗传距离的关系也都不显著。

表3 排除空间自相关影响后最佳阻力模型Table 3 Best resistance models when partialling out the effect of IBD

3 讨论

3.1 遗传多样性

遗传多样性是衡量生物种群健康程度的一项重要指标。研究发现,虎豹国家公园梅花鹿种群每个位点平均有7.56个等位基因,平均观测杂合度为0.689,平均期望杂合度为0.721。梅花鹿各地区亚种的遗传多样性已有许多报道(表4),刘鑫鑫[30]、王美楠[31]、朱明月[32]的目标种群均为野生东北亚种,且采样点与本研究采样区域临近,均位于东北虎豹国家公园内或临近区域。本研究期望杂合度数值略低于或等于以上三项研究。东北梅花鹿在珲春保护区、汪清保护区、黑龙江绥阳老爷岭保护区等地呈现连续分布的格局,该地区没有阻碍梅花鹿扩散的高速公路、大型山脉或大型河流。同时,各项研究采样时间也在5年之内,可认为该地区梅花鹿种群没有较大波动。

表4 不同地区野生梅花鹿基于微卫星的遗传多样性比较Table 4 Comparison of genetic diversity based on microsatellite of wild sika deer in different regions

通过与其他地区梅花鹿亚种对比发现,东北梅花鹿期望杂合度普遍高于其他亚种(表4)。东北梅花鹿可能比其他地区梅花鹿亚种保存有更多的遗传多样性。

栖息地范围以及生境质量可能是影响梅花鹿遗传多样性的一个重要因素。我国境内的野生东北梅花鹿种群靠近中俄边境线,主要分布在东北虎豹国家公园东部地区。而在俄罗斯一侧的滨海边疆区(Primorsky Krai)存在着西南、中部、东南等3个相互联系着的大型梅花鹿种群,其中,西南种群分布范围位于中俄边境附近[35]。经粗略估计,俄罗斯梅花鹿种群的分布面积约为中国境内野生东北梅花鹿分布区的17倍。中国的梅花鹿种群很可能是俄罗斯一侧梅花鹿的扩散种群。同时,俄罗斯滨海边疆区人烟稀少,人类对生态环境的干扰较小,有利于梅花鹿种群的繁衍和扩散。更大的分布面积和更小的人类干扰则意味着可以供养更大的梅花鹿种群,就可能会产生较高的遗传多样性。

在中国的其他梅花鹿分布区,梅花鹿原有生境逐渐被人类挤占,栖息地破碎化不断加剧。例如,梅花鹿四川亚种种群被限制在四川若尔盖县冻列、崇尔、热尔、占哇、降扎5个地区[36];梅花鹿华南亚种种群分布区被分割为位于江西、安徽和浙江的桃红岭、天目山、黄衮山、大会山等多个相互隔离的区域内[37]。

3.2 遗传结构及影响因素分析

通过分析该种群的遗传结构可知,本研究区梅花鹿种群没有明显的遗传结构,这可以由研究区范围内的地理、水文、人类干扰情况所解释。该地区不存在大型山脉、大型河流和高等级的道路,因此也没有能够影响梅花鹿扩散的客观障碍。在梅花鹿分布区内,虽然存在一些硬化道路,但在野外可经常观察到有鹿群穿越。此外,保护区内的居民点规模较小,也不太可能影响梅花鹿扩散。

对该地区东北梅花鹿种群10km范围内的空间自相关分析表明,在第1个(0—1km)、第9个(8—9km)、第10个(9—10km)空间等级表现为显著的空间自相关系数,而其他距离等级的个体之间没有显著的相关关系。这说明遗传相关性仅存在于临近的个体之间。据此可推测,该地区梅花鹿扩散距离为1km左右。

利用GPS项圈方法,Kim等[38]得出韩国俗离山国家公园的梅花鹿的年度家域面积为(2.24±1.50)km2;Yu[39]在日本日光国立公园测得雄鹿夏季家域面积(1.93±0.51)km2,雌鹿夏季家域面积(1.12±0.57)km2;而Takafumi[40]发现日本钏路湿地的梅花鹿平均家域面积为(6.8±1.8)km2。有研究表明,梅花鹿家域面积随着单位面积内个体数量的增加而减小。虽然梅花鹿存在季节性迁徙现象,但是梅花鹿一般在生长季又会迁回并利用原来的栖息地,具有较强的家域固定性[41],因此,梅花鹿在交配时可能也会选择空间临近的个体。本研究结果在一定程度上与前人研究相符。

在梅花鹿的栖息地,除空间距离以外,还有很多其他环境因素可潜在影响梅花鹿种群的遗传模式。过去大量研究集中于梅花鹿的生境选择偏好,而基本没有研究探讨环境因素对梅花鹿遗传格局的影响。本研究中,只有坡度被证明与该地区梅花鹿种群遗传多样性的空间格局相关。对梅花鹿生境选择和占域研究表明,梅花鹿倾向于利用坡度小于25°的区域[42],坡度平缓的区域,梅花鹿的占域率也相对较高[43]。这与梅花鹿生存策略相一致。梅花鹿倾向于在缓坡进行求偶、交配等活动,为了节约能量,梅花鹿也倾向于从平缓地带扩散。

纳入研究的其他几个变量被认为是影响动物种群遗传模式的潜在因素,如Ma等[21]发现地表起伏率和坡向影响了秦岭大熊猫的遗传格局;美洲旱獭的遗传格局受到海拔、道路、树木胸高直径、植被覆盖等因素的影响[44]。同时,河流、坡向、水系等因素对美洲鼠兔种群的基因流造成显著影响[22]。然而,当栖息地高度连接或栖息地与非栖息地之间的阻力差异较小,即景观结构无法对基因流造成限制时,往往无法检测到景观遗传效应[45]。

本研究虽然为每个变量构建一系列不同的阻力模型,但仍然未能探明其余变量与梅花鹿遗传格局的相关性。其他几个环境变量对梅花鹿遗传格局的影响可能受到遗传取样策略、景观要素分布格局等因素的共同影响。由于国家公园西部梅花鹿分布较为稀少,粪便采集难度大,因此本研究未能将研究区向西拓展,未能将国家公园西部更大范围、更为复杂的各种环境特征纳入景观遗传学分析中,客观限制了本研究的深度和广度。Manel等[46]研究表明,样本量的多少以及样本的空间布局将极大地影响后续空间统计方法检测模式的能力,非最佳采样方案需要更高的样本量。此外,有关模拟研究证明,一种景观变量在栖息地和非栖息地间若存在较高的对比度,则对该变量对种群遗传的影响的检测将更容易[47],而本研究所在区域未见有明显的边界将梅花鹿种群相互分隔开来。

遗传位点取样的偏差也可能引起检验偏差。当微卫星位点存在无效等位基因时,即在该位点存在不能成功扩增的等位基因,一方面可能低估群体的遗传多样性,另一方面也容易使亲缘关系分析结果产生偏差。尽管如此,本文的遗传多样性估计结果基本和文献报道一致(表4)。因此,尽管少部分使用的微卫星位点存在无效等位基因对结果存在影响,但这种影响是有限的。

若想进一步验证梅花鹿的扩散是否倾向于避开陡坡或与其他变量相关,研究可结合无线电项圈技术以及运动生理学相关研究,如将由项圈得到实际运动路径与模型模拟路径进行比较分析,或从生理学角度探讨穿越不同景观要素时给动物带来的能量损耗。为进一步研究东北虎豹国家公园各种环境变量对动物遗传格局的影响,可选择旗舰物种东北虎、东北豹,或全域分布的狍、野猪等作为研究对象,可进一步为野生动物扩散廊道设计提供指导。高通量测序技术的发展为检测大量遗传标记与环境变量之间的关联提供技术支撑,景观基因组学相关方法将弥补传统景观遗传学方法的不足,可以更加高效和准确地检测动物种群的微进化格局。