基于重测序技术的塔河马鹿特异性SNP位点筛选

邓伊华，王天娇，王洪亮，董依萌，刘 欣，邢秀梅

（1.东北林业大学野生动物与自然保护地学院，哈尔滨 150006；2.中国农业科学院特产研究所,特种经济动物分子生物重点实验室，长春 130112）

中国养鹿业历史悠久，鹿资源丰富，是世界上养鹿最早、鹿产品加工及应用最广泛的国家之一。中国有8个马鹿亚种，其中新疆分布有3个马鹿亚种，分别为塔里木马鹿（Cervuselaphusyarkandensis）、天山马鹿（Cervuselaphussongaricus）和阿尔泰马鹿（Cervuscanadensisasiaticus）[1]。塔里木马鹿是培育优质高产马鹿的优势资源，在中国近70年的马鹿养殖业中做出了巨大贡献。自20世纪50年代开始，从捕捉野生仔鹿进行驯化和饲养繁殖，到20世纪70年代已经完全依靠自繁自育，最终成功选育的第一个马鹿品种，当地称为塔河马鹿[2]。现代分子生物学研究表明，塔里木马鹿为中国马鹿的祖先，也是世界马鹿的祖先[3]。塔里木马鹿独特的环境适应性决定它们只能生存在塔里木河流域，其耐粗饲、耐受炎热和干燥环境特性一直是学者关注的焦点。由于缺少对塔河马鹿品种遗传资源保护的安全意识，养殖者在追求高产的同时，忽视了对塔河马鹿资源的保护，与天山马鹿和阿尔泰马鹿杂交严重，受胎率低、幼仔存活率低等问题一直未得到解决。系统开展塔河马鹿种质资源资源鉴定工作能合理利用塔河马鹿资源并有效保护塔河马鹿。

近几年来，随着测序技术的不断发展，以第二代测序技术（NGS）为主的高通量自动化测序技术在单核苷酸多态性（SNP）位点的开发和利用方面逐渐凸显出巨大优势。在此基础上，通过全基因组重测序技术对已知基因组序列的不同个体进行基因组重测序，比较参考基因组和全基因组重测序序列，获得大量的SNP位点、拷贝数变异（CNV）、插入缺失位点（InDel）和结构变异位点（SV）等遗传特征[4-6]。全基因组重测序技术在动物变异检测[7]、性状定位[8]、基因挖掘[9]、遗传进化[10]等多个领域均取得了丰硕的成果[11]。Ba等[12]采用ddRAD-seq技术检测到30个圈养个体（7头梅花鹿、6头马鹿和17头杂交鹿）的约320 000个全基因组SNPs，通过筛选并观察杂合度，首次报道了梅花鹿和马鹿的特异性SNP位点。董世武等[13]通过GBS技术对梅花鹿、马鹿及其杂交后代（F1和F2） 4个群体共226个样本进行测序，分别识别出马鹿特异性SNPs位点474个，梅花鹿特异性SNPs位点558个，为梅花鹿、马鹿和杂交鹿的鉴别提供了依据。塔依尔江·麦麦提等[14]利用全基因组重测序技术对塔里木马鹿干旱环境适应相关基因进行筛选和研究，成功筛选出塔里木马鹿过氧化还原酶3（PRDX3）基因。Fan等[15]对500多头梅花鹿、马鹿和杂交鹿进行全基因组重测序，筛选得到1 300多万个SNPs位点，通过生物信息学和芯片设计原则，最终选择1 000个特异性SNPs位点用于芯片合成，随机选取266头鹿进行验证，结果表明，1K梅花鹿SNP芯片能够准确地区分梅花鹿、马鹿和杂交鹿。从简化基因组测序筛选特异性位点到利用重测序技术研发相关基因芯片，测序技术正逐步应用于梅花鹿和马鹿的鉴别研究、基因挖掘等方面。塔河马鹿是所有马鹿种群中规模最大的群，然而目前种群数量不足万只，数量还在持续下降，而且人工饲养的马鹿存在杂交，导致其品质变化，多样性减少，塔河马鹿遗传资源流失严重。因此，需要对塔河马鹿这一宝贵的地方品种资源进行合理的保护与利用，系统地开展塔河马鹿种质资源鉴定工作。

本研究利用全基因组重测序技术，通过对32份塔河马鹿进行全基因组测序，结合实验室前期获得的塔河马鹿、天山马鹿和阿尔泰马鹿共88份重测序数据（未发表），以组装到染色体级别的梅花鹿基因组序列作为参考，深入挖掘塔河马鹿与天山马鹿、阿尔泰马鹿之间全基因组水平上的差异，筛选出特异性强，稳定性高的塔河马鹿特异性SNP位点，构建塔河马鹿分子遗传标记，为塔河马鹿种质资源收集、整理、保护和评价工作奠定基础。

1 材料与方法

1.1 材料

国家特种经济动物资源共享平台收集并保存的1999年在新疆地区采集的32份塔河马鹿血液样本。

1.2 方法

1.2.1 血液基因组DNA的提取 采用酚-氯仿法提取32份塔河马鹿血液基因组DNA，提取的DNA用1.0%琼脂糖凝胶电泳和核酸定量蛋白仪（Agilent 5400，安捷伦）进行质量检测，经检测合格后-20 ℃保存备用。

1.2.2 DNA文库的构建和测序 检测合格的DNA样品用于DNA文库构建，32份塔河马鹿DNA分别建库和测序。建库和测序由北京诺禾致源科技股份有限公司完成。将DNA样品通过Covaris超声波破碎仪随机打断，经末端修复、加A尾、加测序接头、纯化、PCR扩增等步骤完成整个文库制备。文库检测合格后，把不同文库按照有效浓度及目标下机数据量的需求，利用Illumina HiSeq/MiSeq对32个样本进行双末端低深度测序，随后对原始下机数据进行预处理，测序得到原始测序序列（sequenced reads）或者原始数据（raw reads：含有带接头的、低质量的reads）。将得到的序列进行GC含量分布检测和测序质量分析。为了保证信息分析质量，对raw reads进行过滤，得到clean reads，基于clean reads进行后续分析。

1.2.3 测序数据比对 使用BWA-MEM软件[16]将过滤后的clean reads和实验室前期获得的重测序数据一同与组装到染色体水平的梅花鹿参考基因组序列（MHL_V1）进行比对，并对本次32份样品的比对率（mapping rate）、覆盖度（coverage）和测序深度（depth）进行统计。

1.2.4 基于SNP位点的生物信息学分析

1.2.4.1 SNP位点检测和注释 使用SAMTOOLS软件[17]对BAM文件进行排序，删除重复读取的数据。利用GATK 4.0.2软件进行变异检测，过滤结果输出到VCF文件保存。使用软件SNPEff对得到的SNP位点进行注释，统计SNP位点在33条染色体上的分布。

1.2.4.2 分子进化树构建 使用SNPhylo软件对所有样品进行关系推断，构建分子进化树，统计候选特异性SNP位点。

1.2.4.3 候选特异性SNP位点筛选 将1999年收集整理的塔河马鹿样品作为一个群体，其余马鹿作为另一个群体，使用VCFTOOLS[18]计算Fst值，将计算得到的Fst值按照降序排序，设定阈值Fst≥0.25,淘汰不符合条件位点，得到每条染色体剩余SNP数目用于构建塔河马鹿遗传标记。

1.2.4.4 塔河马鹿特异性遗传标记构建 统计筛选后每条染色体上的SNP数目并将所有SNP位点进行编码，基因型AA编码为0；基因型AB编码为1；基因型BB编码为2，基因型缺失编码为-9；编码后的SNP位点数据构成矩阵，使用R语言4.0.4版本prcomp函数[19]进行主成分分析（PCA），计算每条染色体上SNP标记的得分，根据排名选取每条染色体排名前1 500名的SNP，统计33条染色体筛选出来的SNP位点数，计算PC1和PC2两个主成分右奇异向量最大的列。选取前100个SNPs构成特征SNP子集。提取前100个SNPs集合的特征值，计算协方差矩阵。分别对33条染色体筛选得到的SNP位点和前100个SNPs位点绘制主成分分析聚类图。 最后统计前100个SNPs在染色体上的分布情况。

2 结 果

2.1 血液基因组DNA质检结果

1.0%琼脂糖凝胶电泳检测结果表明，得到的塔河马鹿DNA目的条带单一，无明显降解和杂质污染。核酸定量蛋白仪检测其浓度和纯度都在要求范围之内，符合后续建库要求。其中DNA样品总量在0.09182～4.14943 μg之间。

2.2 测序数据统计与质量评估

32份塔河马鹿血液基因组DNA样品经全基因组重测序，共产生raw reads 869 841 762 900 bp，去除低质量序列后得到clean reads 868 791 354 600 bp，平均每个样品27 149 729 831.25 bp。测序质量Q20≥96.96%、Q30≥92.09%，GC含量在43.29%～46.30%之间，没有发生碱基分离的现象。此次建库测序成功，测序质量较高，获得的数据量、质量和序列长度均符合后续数据分析要求。

2.3 测序数据比对情况

以组装到染色体级别的梅花鹿全基因组为参考序列（MHL_V1），对每个样品的比对率（mapping rate）、覆盖度（coverage）和测序深度（depth）进行统计，32份样品的测序数据的平均比对率为98.06%，平均覆盖度为97.66%，平均深度为6.787。

2.4 基于SNP位点的生物信息学分析

2.4.1 SNP位点的特征分析 将32份样品测序得到的有效数据与试验前期得到的测序数据与参考基因组（MHL_V1）进行比对，过滤后得到20 139 122个高质量的SNPs位点。对每条染色体上的SNP位点分布情况进行统计（图1）。SNP位点基本随机均匀分布在每条染色体上，其中4号染色体上分布最多，有1 132 360个，26号染色体上分布最少，有317 705个，SNP在染色体上的分布可能与染色体的长度有关，染色体越长，SNP位点数量越多。用SNPEff统计所有染色体上变异类型的分布情况表明，1 156 947（6.76%）个SNPs分布在上、下游区域；基因间的变异最多，有11 118 374个（65.02%）；外显子区域有128 027个（0.75%）；内含子区域有4 701 559个（27.47%）（图2）。

图1 每条染色体上的SNP数目Fig.1 The number of SNP on each chromosome

图2 变异类型在染色体上的分布Fig.2 Distribution of variant types on chromosomes

2.4.2 基于SNP位点构建分子进化树 基于SNP位点构建分子进化树，结果显示， 塔河马鹿、天山马鹿和阿尔泰马鹿区分明显，其中塔河马鹿单独汇聚成一支，天山马鹿和阿尔泰马鹿聚类位置相近，根据进化树聚类位置选择样本进行分析，过滤后位点数为12 050 781（图3）。

图3 基于SNPs的分子进化树Fig.3 Molecular evolutionary tree based on SNPs

2.4.3 候选特异性SNP位点的筛选 将1999年收集整理的塔河马鹿样品作为一个群体，其余马鹿作为另一个群体，使用VCFTOOLS软件淘汰不符合条件位点，统计每条染色体的候选特异性SNP位点，筛选后的SNP在染色体上的分布结果显示，剩余位点数544 717个。其中SNPs位点在4号和5号染色体上分布最多，其余染色体上位点分布较为均匀（图4）。

图4 Fst筛选后的SNP在染色体上的分布Fig.4 Distribution of SNP on chromosomes after Fst screening

2.4.4 塔河马鹿特异性遗传标记的构建 从每条染色体选取排名前1 500的SNPs位点，33条染色体共筛选得到49 500个SNPs位点。49 500个SNP位点及其前100个SNPs位点主成分分析结果表明，当位点数下降至100时，1999年的塔河马鹿依然能准确与天山马鹿和阿尔泰马鹿区分开，区分效力基本没有下降。大多数阿尔泰马鹿分布PC1的右边，少量分布与天山马鹿相近，群体分布较为分散，个体血源可能不纯，此外，1999年的塔河马鹿和近年的塔河马鹿在聚类图上区分明显，大部分现在的塔河马鹿与天山马鹿分布相近，天山马鹿分布在PC1的左上角，个体联系较为紧密（图5、6）。

图5 49 500个SNPs的主成分分析Fig.5 PCA of 49 500 SNPs

图6 前100个SNPs的主成分分析Fig.6 PCA of top 100 SNPs

通过计算筛选后选取前100个SNPs位点作为塔河马鹿特异性分子遗传标记，前100个SNPs位点在染色体上的分布可以看到，15和30号染色体的分布最多，其中15号染色体上分别有3个C/C、2个T/T,G/G和A/A基因型各1个，30号染色体有4个T/T，C/C、G/G、A/A基因型各1个。1、4、5和23号染色体上无分布。得到的塔河马鹿优势基因型在大部分染色体上均有分布，分布情况良好（图7）。

图7 前100个SNPs在染色体上的分布Fig.7 Distribution of top 100 SNPs on chromosomes

3 讨 论

中国家养马鹿遗传背景单一，塔河马鹿由于高强度的人工选择和近交导致品种退化现象日益严重。近年来,随着杂交改良的推进，产茸重量上的绝对优势致使杂交鹿在养殖上占据的市场空间越来越大，养殖者盲目追求鹿茸产量而无序杂交，造成了塔河马鹿纯种资源的匮乏。杂交鹿从肉眼上难以与塔河马鹿区分，给塔河马鹿的保种、育种工作带来了极大困难。目前二代测序技术凭借其高通量和低成本的优势在动物基因组测序中得到了广泛的应用，在全基因组水平上筛选特异性SNP位点[20]，结合相关检测技术和生物信息学分析，推断动物品种构成，预测杂种或纯种[21]。朱涛等[22]选取7个群体的56只鸭进行重测序，筛选出686个存在于特定群体中的SNPs和InDel并随机选取7个进行了验证，结果完全符合预期。岳阳等[23]对小梅山猪和其他11个猪品种进行简化基因组测序，筛选出5个小梅山猪的特异性SNPs位点并进行验证，建立小梅山猪的SNP分子标记鉴别方法。本研究以大量马鹿的重测序数据为基础，采用本团队组装到染色体的梅花鹿基因组作为参考，更准确地检测塔河马鹿SNP的变异情况，利用生物信息学软件分析SNP位点，提高计算效率，降低研发成本，为后续塔河马鹿鉴别的相关技术开发及应用提供多样化的选择。

本试验对32份塔河马鹿DNA样本进行全基因组重测序，结合试验前期的56份马鹿重测序数据，与组装到染色体级别的梅花鹿参考基因组进行比对，通过生物信息学分析并得出结论。结果表明，32个塔河马鹿样品测序和比对情况良好，数据量符合预期。筛选的SNP位点基本分布在基因间和内含子区域，与Wang等[24]研究广西地方鸭中筛选的SNP位点变异分布结果一致，SNP的分布差异可能与DNA在不同基因区域的功能差异有关[25]。通过构建分子进化树能进一步分析马鹿的遗传多样性并挖掘个体间的差异[26]，从分子进化树上可以看到天山马鹿和阿尔泰马鹿分支相近，塔河马鹿单独聚成一支，与天山马鹿和阿尔泰马鹿相距较远，与苏莹等[27]、Shao等[28]研究结果一致。王小鹏[29]利用60K基因芯片结合主成分分析和系统进化树筛选特异位点构建猪特异性遗传标签,刘月丽等[30]利用主成分分析和随机森林算法结合筛选高质量SNP位点用于羊的品种分类。从聚类图中可以看到1999年的塔河马鹿与其他马鹿区分明显，与现在未知血源状况的塔河马鹿也能准确区分，进一步说明现在的塔河马鹿样本不纯，已经进行了杂交，混有天山马鹿或阿尔泰马鹿的血源。通过筛选得到的塔河马鹿特异性SNP位点，不仅将早年的塔河马鹿与天山马鹿和阿尔泰马鹿区分开，还能与现在杂交的塔河马鹿区分。SNP标记作为第3代分子标记被认为是目前应用前景最好的遗传标记，其在动植物基因组中均分布广泛，稳定性高且富有代表性。SNP标记不再以DNA片段的长度变化作为检测手段，而是直接以序列变异作为标记，完全摒弃了凝胶电泳而采用测序技术或最新的DNA芯片技术。本研究采用全基因组重测序技术，测序深度和覆盖度高，得到的塔河马鹿特异性位点范围广、准确性高，另一方面选取的参考群体范围广，32份塔河马鹿重测序数据结合之前56份马鹿的重测序数据（11头塔河马鹿、11头天山马鹿、34头阿尔泰马鹿）使特异性位点的筛选结果更加准确。

4 结 论

本研究结果表明，32份塔河马鹿重测序数据过滤后筛选得到了20 139 122个高质量的SNPs位点，建树过滤后SNP位点数为12 050 781个，过滤后共得到544 717个SNPs位点用于塔河马鹿分子遗传标记构建，最终筛选得到了100个具有高度多态性、稳定性良好的SNPs位点，结果为塔河马鹿的精准鉴别和核心种质的筛选提供理论依据，为塔河马鹿种质资源的收集、整理、保存和评价工作的顺利开展奠定坚实基础。