基于重测序数据的昌都黑山羊遗传多样性及群体结构分析

阿旺措吉，仁青措姆，黄舒泓，周先辉，德 吉，索朗达，王小龙，吴玉江，巴 贵*

(1.西藏自治区农牧科学院畜牧兽医研究所，西藏 拉萨 850000;2.西北农林科技大学 动物科技学 院，陕西 杨凌 712100)

藏山羊(Caprahircus)是中国独特的种质资源，在青藏高原的农业、文化、历史、经济乃至宗教等方面都起到了重要的作用，其主要分布在西藏自治区、四川省阿坝和甘孜自治州、青海省玉树自治州、甘肃省甘南、果洛自治州及新疆部分地区等[1]，是中国最古老的羊品种之一，最早记录于公元前3 300-2 000年[2]，外国商人称谓的克什米尔山羊(Cashmere goat)即为藏山羊。长期以来，由于复杂多样的环境、高海拔和极端气候，藏山羊具备了独特的遗传特性[3]，加之其群体大、山羊绒产量高、肉品质及皮革品质优良，现已成为中国宝贵的基因库[4]。

全基因组重测序基于高通量基因测序技术，并依托计算机对研究群体基因组与已发表参考基因组进行比对，可以在基因水平对不同目标间的基因差异进行对比。目前，已有大量使用基因组重测序的方法与研究。如李旭静等[5]利用全基因组重测序对绵羊经济性状基因进行筛选，注释了多个与毛、乳、肉产量相关的选择区域。王统苗等[6]使用全基因组重测序技术对鸭群的遗传资源进行结构分析，将多个不同地区的鸭种资源进行分类与合并。

前人通过应用分子遗传学标记、基因组学、生物信息学等[7]方法，积累了大量关于藏山羊遗传资源多样性的资料。如王杰等[8]采用SSR分子标记对高原型藏山羊、山谷型藏山羊进行遗传多态性研究，发现藏山羊群体遗传多态性丰富且对生存环境的适应能力较强。王永等[9]应用ISSR标记分析了西藏日土藏山羊的遗传多样性，结果显示该品种藏山羊的遗传多态性较丰富，个体间虽有差异但同质性较好。邓娟等[10]通过对细胞色素 b基因(Cytb)全序列进行扩增和测序，研究了西藏地区藏山羊的母系起源及遗传多样性。

本研究收集西藏昌都3个不同地区的藏山羊群体血液DNA并测序。用SNP位点分析、遗传多样性分析、群体遗传结构分析、遗传关系分析及主成分分析等多种分析方式，从基因组水平解释了藏山羊的种质特性，也为后续藏山羊重要生产性状的基因定位、分子育种及遗传资源的保护和利用提供重要的理论依据。

1.1 试验动物

试验样本采自西藏昌都市边坝县边坝镇洛亚码羊场(LYM，n=9)、普玉一村(PY，n=9)及芒康县朱八龙乡(MK，n=10)，共28个样本。采集健康个体的静脉抗凝全血样本，-80 ℃冰箱中保存备用。

1.2 试验方法

1.2.1 DNA提取及建库 血液基因组DNA样品制备：使用购买自华大基因(BGI)的试剂盒，按照说明书将-80 ℃藏山羊静脉血解冻后提取基因组DNA，并冷藏保存。

DNA文库的构建：使用Nanodrop超微量核酸分析仪检测血液基因组DNA样品的浓度，并使用凝胶电泳检测DNA样品的纯度。使用超声波仪打断0.001 mg基因组DNA，选择0.2～0.4 mg大小基因片段。在DNA片段3"端链接DNA接头。使用PCR扩增加入DNA接头的片段并提纯回收，将双链PCR产物解链，破坏未被环化的DNA分子，获得DNA文库。

1.2.2 测序及参考基因组比对 使用华大基因的DNBSEQ-T7基因测序仪器平台进行测序，测序工作由华大基因公司完成。下机后的Raw reads使用trimmomatic软件进行质控分析，测序数据过滤得到Clean reads。使用BWA软件将测序数据比对到参考基因组，结合Samtools软件转换数据格式，Picard软件对比对文件去重，得到最终分析所用文件。

1.3 数据分析

1.3.1 SNP位点分析及质控 处理后比对文件，通过GATK软件的HaplotypeCaller模块进行SNP位点鉴定，之后经过Bcftools及VariantFiltration模块进行质控和过滤，再通过SelectVariants模块以及Bcftools合并各个体的GVCF文件[11]。

1.3.2 遗传多样性分析 利用PLINK软件，设定参数滑动窗口大小为50 kb，步长为20 kb，计算3个群体的核苷酸多样性PI值[12]。通过PopLDdecay软件计算3个群体的LD衰减，得到群体内的LD连锁不平衡信息[13]，最后利用R语言脚本以及perl语言脚本处理结果生成可视化图形。

1.3.3 群体遗传结构分析 将包含SNP信息的VCF文件通过PLINK软件转换为BED及PED文件[12]，输入到Admixture软件中推断群体遗传结构，设置群体亚群数目K=2到12计算交叉验证误差(Cross-validation error)，通过绘制Cross-validation error图确定最佳分群数[14]，利用R语言脚本及ggplot2软件包对结果生成可视化图形[15]。

1.3.4 遗传关系分析及主成分分析 利用PLINK软件生成PED、MAP及BED文件，根据LD进行过滤，通过编写的perl脚本计算遗传距离，生成Neighborhood-Join法构建的进化树文件，得到的meg文件输入到MEGA软件生成可视化的NJ树[16]，利用在线网站Interactive Tree Of Life(https://itol.embl.de/)进行进化树可视化的美化[17]。

使用EIGENSOFT的smartpca程序进行主成分分析[18]，得到的主成分信息利用R语言脚本及ggplot2软件包生成可视化图形[15]。

2.1 藏山羊重测序数据质量分析及参考基因组比对

本次研究利用华大基因DNBSEQ-T7平台对28个藏山羊个体血液样本进行全基因组重测序(表1)。所得下机数据初始Raw reads数量范围为155 118 313～272 533 434，经过软件质控过滤，Clean reads数量范围在151 935 132～267 865 582，过滤比率在97.83%～98.35%之间，过滤比率是指数据清洗后reads数量与原始数据(Raw reads)之间的比值，比值越高则测序质量越高。Q30碱基比率均高于95%，测序质量Q值指的是测序过程碱基识别(Base Calling)过程中，对所识别的碱基给出的错误概率，Q30即为碱基错误率0.1%，正确率99.9%，因此Q30是评估测序质量标准的重要一环，且GC含量也位于40%～50%之间，说明样本的建库良好，测序质量达到重测序分析标准。

表1 28个藏山羊个体重测序数据质量控制信息

使用bwa程序mem模块将28个藏山羊重测序数据比对到山羊参考基因组ARS1(NCBI accession: GCF_001704415.1)，比对情况如表2所示。使用测序数据质控后reads数与参考基因组reads数的匹配比对率(Mapping rate)以及平均覆盖度(Mean coverage)作为评估数据可信度、以及分析正确度的指标，通常比对率与覆盖度越高，则分析数据的可信度越高。比对上reads数在303 831 381～535 839 901之间，比对率在99.61%～99.93%之间，平均覆盖度在15.2262×～26.5768×之间，说明比对情况良好，可用于后续SNP位点分析及群体结构分析。

表2 28个藏山羊个体重测序数据与参考基因组比对情况

2.2 藏山羊全基因组SNP位点鉴定

经前期质控比较及SNP筛选，在3个采样群体28个个体中共检测到13 726 331个SNP位点，不同的染色体上SNP位点数目见表3。由于SNP位点是DNA序列上发生碱基改变的位点，这些位点由于突变的不确定性而在不同群体与个体中具有鲜明特点。总体来说，鉴定出的SNP位点数量与染色体的长度呈正相关。

表3 28个藏山羊相对于参考基因组的SNP位点信息

2.3 藏山羊遗传多样性分析

对3个群体进行核苷酸多样性分析，最终得到LYM、MK、PY3个采样群体核苷酸多样性(PI)分别为0.00182、0.00178、0.00183，除个别位点存在差异外，整体分布情况接近，差异极小，表明3个采样群体SNP情况接近，如图1所示。同核苷酸多样性分析结果类似，3个采样群体的连锁不平横衰减速度也极为相近，整体趋势及进化速度较为一致，证明采样群体间关系紧密。

图1 昌都黑山羊采样群体核苷酸多样性

2.4 藏山羊群体遗传结构

使用admixture分析对3个采样群体，以及阿里、班戈、措勤、林芝、尼玛、那曲、日土、山南共8个西藏其它地区藏山羊品种或群体进行分析。由于群体遗传结构是由群体的表型种类、频率、基因种类、基因型种类及频率等组成的。而地理环境因素之间的差异，导致群体间产生不同的变异与分化，根据变异的数量与方式不同，可以将群体按照亲缘关系划分为不同亚群。对最优亚群数K值的评估，结果如图2所示，当交叉验证误差到达最小值时，K=5为假定的最佳分群数目。K=5时群体遗传结构如图3所示，3个采样群体的群体遗传结构类似，但LYM采样群体及PY采样群体存在一定的杂交特征，这可能是由于羊群组群时混合其他群体导致的。

图2 不同K值下的交叉验证误差值分布

图3 K=5时藏山羊种群遗传结构图

2.5 藏山羊遗传关系及主成分分析

通过加入尼玛(NM)、措勤(CQ)、林芝(LZ)、日土(RT)、山南(SN)、阿里(AL)、班戈(BG)、那曲(NQ)共8个西藏其它地区藏山羊品种或群体，同本研究的3个采样群体进行Neighborhood-Join法构建进化树以及主成分分析。NJ树的拓扑结构如图4所示，本研究中3个采样群体所在分支离西藏其他地区藏山羊群体遗传分支距离较远，且3个采样群体共同聚为一支，群体间亲缘关系更近，可能为同一群体。

图4 藏山羊NJ进化树分析结果

使用主成分分析，3个采样群体间相互距离较近，并与其他藏山羊群体明显分离，如图5所示。进化树分析以及主成分分析结果与遗传结构分析结果一致，说明本研究3个采样群体亲缘关系以及遗传距离较近，为共同品种藏山羊。

图5 藏山羊主成分分析结果

在公元前3 000-公元前2 000年便有了关于藏山羊的记载[2]，由于西藏地区海拔高、温度低的环境特点，西藏地区山羊的适应性与低海拔地区山羊品种差异明显，因此藏山羊对于中国其他山羊品种的生产性能以及相关环境适应性研究具有重要作用。同时，西藏高原地区的地理差异导致不同地区之间存在长期的地理隔离[19]，可能在藏山羊种群内产生不同亚种。王杰等[8]、王永等[9]和秦国庆等[1]采用分子标记法对西藏地区山羊的种群多样性进行分析，发现不同地区藏山羊群体差异较大，具有较高的种群内多样性。

遗传多样性是生物多样性的一个层次，是某地区内所有生物遗传信息的集合，是物种在长期进化过程中不断发生变异所产生的特征。遗传多样性越丰富，自然环境的抗逆性也越强。本研究利用了全基因组SNP筛选获得了高质量的SNP集合，并利用此集合进行了核苷酸多样性分析、连锁不平横衰减计算。并加入了西藏地区其他品种山羊SNP数据集进行群体遗传结构分析以及进化树的构建，研究结果对于中国的遗传资源保护具有重要意义。

昌都黑山羊是西藏地区东部地区特有的高原型山羊群体，具有抗逆性强、板皮厚满、皮毛质软细小、瘦肉率较高、肉质鲜嫩、膻味小等优良性状，在当地羊产业发展过程中发挥着重要作用。本研究使用昌都黑山羊的3个群体(LYM、PY和MK)共28个藏山羊个体的重测序数据，经遗传多样性分析，3个群体的SNP分布与数量均相似，且连锁不平衡的整体趋势与衰减速度也较为一致，可以认定昌都地区3个群体为同一品种。后续通过加入8个西藏地区其他品种的SNP数据进行群体遗传结构分析、NJ法进化树构建以及主成分分析，发现昌都黑山羊的3个群体与西藏地区其他品种遗传距离较远且分化明显。

综上所述，3地采样的昌都黑山羊群体为同一品种，且品种内存在一定的杂交。遗传进化关系表明，该品种与西藏其他地区的山羊品种遗传分化明显，可被认为是独立的西藏山羊群体，为后续新遗传资源的鉴定提供了重要的理论支撑。本研究通过28个昌都黑山羊重测序数据获得了高精度的SNP数据集，可用于后续昌都黑山羊分子育种及生物信息学分析，如BLUP育种值计算、种间基因渗入研究、环境适应性研究等，为发展昌都地区优势特色畜牧业、提高昌都黑山羊的良种化程度提供理论依据。

本研究发现昌都黑山羊具有丰富的遗传多样性；
昌都地区3个黑山羊群体内亲缘关系与遗传结构较近，且相较西藏其他8个地区藏山羊群体遗传分化较大。研究结果对于中国家畜遗传资源保护具有重要意义。