山东小毛驴全基因组选择信号检测

陈建兴，童家兴，张孝忠，张向阳，王宇鑫，孙玉江,5

(1.赤峰学院 化学与生命科学学院，内蒙古 赤峰 024000； 2.青岛农业大学 动物科技学院，山东 青岛 266109；3.内蒙古东阿黑毛驴牧业有限公司，内蒙古 赤峰 024328； 4.东阿阿胶股份有限公司，山东 聊城 252200；5.东营职业学院，山东 东营 257091)

山东小毛驴(SDL)过去主要产地在胶东半岛、沂蒙山区和鲁中平原[1]，旧称胶东小毛驴，因现在中心产区在海阳等地，又称海阳小毛驴[2]。山东小毛驴体质外形与华北驴类同，属小型驴，具有体型小、挽力大、耐粗饲、抗逆性强等特点[2]。现代农业机械化和交通业的快速发展，驴的役用地位迅速降低，使得山东小毛驴种质资源迅速衰减，遗传多样性下降，育种潜力逐渐丧失。家畜遗传资源问题是全球性生物资源问题的组成部分，与人类未来的生存与发展紧密相关[3]。驴产业又是我国的特色产业、民生产业和创新产业[4]，随着我国社会经济发展和中国特色社会主义建设的不断推进，独特的驴遗传资源将会越来越珍贵，对山东小毛驴的种质资源的保护开发将会日趋重要。

不同品种驴在体型外貌、生长性能、抗病力和适应性等方面存在较大差异，这些差异是自然选择与人工选择共同作用于一些目标基因定向选择出来的结果，而选择信号是与选择相对应的基因组信息，是选择在基因组上留下的印迹，通常表现为受选择的DNA片段或位点多态的降低或者基因的纯合[5]。随着单核苷酸多态性(SNP)芯片和第二代测序技术成本的不断降低，利用选择信号分析来揭示引起畜禽性状表型差异的遗传机制的相关研究日益普遍。吕世杰等[6]利用选择性清除方法筛选郏县红牛和中国荷斯坦奶牛2个品种间差异的基因组区域，通过与动物数量性状座位(QTL)数据库中牛繁殖性状相关QTLs进行比对，认为CFDP1、CFDP2和FAM204A基因可优先作为牛繁殖性状相关候选基因。PETERSEN等[7]使用固定指数分析寻找了马的基因组的选择信号，研究结果表明，MSTN、ECA11和DMRT3基因受到了选择，另外通过关联分析发现，MSTN内含子的1个SNP和启动子的1个InDel与肌纤维比例有关。樊英智[8]通过对我国6个不同类型的家驴品种群体(共57头)进行全基因组混池重测序，分析了不同类群基因组间的选择信号，找到与驴品种表型如体尺、毛色等性状有关联的11个候选基因。

通过对不同驴品种进行基因组重测序，检测不同驴品种群体间基因组的选择信号的差异，有助于揭示品种群体的进化历史，了解重要表型性状形成的遗传基础。目前，对山东小毛驴基因组选择信号的检测还未见报道。利用群体遗传分化系数(Fst)和核苷酸多样性比值(π ratio)方法，对山东小毛驴与德州驴的三粉类群(DZS)、德州驴的乌头类群(DZW)、广灵驴(GL)、华北驴(NC)等4个驴群体间的选择信号差异进行分析，以期筛选出山东小毛驴的选择信号区域，探讨山东小毛驴特性的形成原因，旨在为山东小毛驴的保护和利用提供参考。

1 材料和方法

1.1 供试动物及样品

本研究采用4个驴品种(5个群体)共60头驴作为供试动物(表1)。对所有样本均采集颈静脉血液10 mL于EDTA抗凝管中，置于-20 ℃冰箱短暂保存后于干冰保鲜盒中快递至美吉生物公司(上海)进行后续研究。

1.2 测序数据

利用超声波将检测合格的样品基因组DNA片段化形成随机片段，对片段化的DNA进行末端修复、3′端加A、连接测序接头后，再利用磁珠吸附富集400 bp左右的随机片段，经PCR扩增形成测序文库。构建好的测序文库通过Illumina HiSeqTM平台进行测序，测序策略为Illumina PE150。测序数据已上传到NCBI的SRA数据库，SRA号为SAMN14484743—SAMN1448 4802。

1.3 SNP、InDel检测与注释

本研究采用BWA软件[9]将高质量测序数据比对到参考基因组序列 (https://www.ncbi.nlm.nih.gov/genome/7038)上，利用GATK软件[10]进行比对后校正，并进行SNP和Small InDel标记的检测；利用SnpEff软件[11]和参考基因组的基因预测信息进行变异功能注释，得到SNP、InDel的功能注释信息。

表1 供试驴信息

1.4 Fst与π ratio计算

Fst计算公式：Fst=(MSP-MSG)/[MSP+(n-1)MSG][12]。其中，MSP为群体间均方差，MSG为群体内均方差，n为校正后平均样本大小。Fst可用来评价群体间的分化程度，该值越接近1说明两群体间分化程度越高，越接近于0说明两群体间分化程度非常有限。对于高质量SNP (次等位基因频率maf不低于0.05，缺失率miss为0，样本测序深度不低于5)，利用VCFtools软件[13]计算了两两群体间的Fst值(2 Mb窗口，10 kb步长滑窗)。Π=∑ijxixjπij，其中，xixj分别代表第i个和第j个序列的对应频率，πij则为2个序列之间不同位点所占百分比。Π是遗传变异的1个量化值，用于表征某一种群多态性的强弱，通常用于衡量种群内或种群间的多样性，该值不依赖于样本大小[14]。每个群体的π ratio值也使用VCFtools软件[13]进行计算，同样是2 Mb窗口，10 kb步长滑窗。

1.5 群体间选择信号检测

本研究采用基于Fst和π ratio的方法对5个驴群体间的选择信号进行检测。Fst和π ratio分别选取阈值0.95和0.05(分位数)，关联Fst和π ratio提取相应候选区域(取重叠区域)，并提取相应区域内的变异位点信息。采用选择性清除方法，分别对SDL与其余4个驴群体进行比较分析，检测到相应的选择信号区域。选择信号强弱的判定，根据VCFtools软件筛选出有选择性消除位点相应的注释，选择注释级别为HIGH而排除MODERATE和LOW的结果。对选择信号区域中存在的基因，通过在GeneCards数据库 (www.genecards.org) 或NCBI Gene (www.ncbi.nlm.nih.gov/gene/?term=)中查询基因注释来确定基因功能。

2 结果与分析

2.1 5个驴群体变异检测结果

5个驴群体的全基因组重测序，共获得740.57 G高质量数据，平均每个样本获得了12.34 G的数据。将高质量测序数据比对到参考基因组之后，总共获得了10 096 033个 SNP(SDL、DZS、DZW、NC、GL分别占72.03%、77.01%、71.62%、74.99% 和70.34%) 和1 311 358个InDel(SDL、DZS、DZW、NC、GL分别占75.50%、79.71%、74.46%、76.86%和73.48%)。整体来看，2种变异(SNP和InDel)的最大比例均出现在德州驴的三粉类群中，而最小比例出现在广灵驴群体中(表2)。然而，杂合SNP和InDel数目以及π ratio的最大值都出现在山东小毛驴群体中，而最低值都出现在华北驴群体中。观测杂合度和杂合SNP数、杂合InDel数一样，最大值出现在山东小毛驴群体中，最低值出现在华北驴群体中。转换颠换比(Ts/Tv)的最大值出现在山东小毛驴群体中，最低值出现在德州驴的三粉类群中(表2)。

表2 5个驴群体的变异和遗传多样性指数

2.2 群体间遗传分化分析

通过计算Fst，群体间的遗传分化程度能够得到估量，计算结果见表3。Fst值从DZS和NC群体间的0.007 80到DZW和SDL群体间的0.011 460。显然，各Fst值都非常低，接近于0(WRIGHT认为[15]，该值处于0～0.05，表明群体间遗传分化很小，可以不用考虑)，表明这5个驴群体间的分化水平极低。

表3 5个驴群体间遗传分化系数

2.3 山东小毛驴与其他4个驴群体的选择性清除结果

选择性清除指新的有利突变会增加其频率并固定下来，导致其相邻核苷酸序列的差异下降或消除的过程[16]。为了检测到可能的选择性消除位点，采用基于Fst和π ratio的方法搜索寻找了驴的基因组中的高度固定的区域，筛选过程见图1，选择Fst值大于0.95分位阈值并且π ratio值小于0.05分位阈值的区域，筛选的具体结果见表4。

蓝色点是筛选的候选区域，Fst值大于0.95分位阈值并且π ratio值小于0.05分位阈值

由表4可知，山东小毛驴与其他4个驴群体共检测到正向选择区域95个，德州驴的乌头类群最多，共30个，华北驴最少，只有17个。为进一步了解这些信号选择区域的功能，对落入选择信号区域的39个基因及基因功能进行了生物信息学分析。山东小毛驴与德州驴的乌头类群比较，群体间信号选择区落入的基因数最多，有13个，与华北驴群体比较信号选择区落入的基因数最少，只有5个(表4)。经统计，与德州驴的三粉类群相比，群体间检测出强选择信号候选基因5个，分别是Fsip1、AHNAK2、CTAGE2、CYP3A12、LOC106830441；与德州驴的乌头类群相比，群体间检测出强选择信号候选基因5个，分别是NKG2DL1、KLK1E2、CTAGE2、FAM170A、LOC106823932；与广灵驴相比，群体间检测出强选择信号候选基因6个，分别是NKG2DL1、AHNAK2、CTAGE2、FAM170A、LOC106823932、LOC106848008；与华北驴相比，群体间检测出强选择信号候选基因2个，分别是CTAGE2、FAM170A。落入强选择信号的10个基因(表4中强选择信号基因列中显示的所有基因)，大多与免疫、生殖以及细胞作用和代谢相关。

表4 山东小毛驴与其他驴群体信号选择区域

3 结论与讨论

本研究采用基于Fst和π ratio的方法对山东小毛驴和德州驴的三粉类群、德州驴的乌头类群、广灵驴、华北驴等5个驴群体进行了选择信号分析。SNP和InDel的最大比例出现在德州驴的三粉类群中，这可能与德州驴自古以来都是优秀大型驴种有关[17]，一直与各地驴种之间存在交流，不断有优秀的个体引入该群体，而SNP和InDel的最小比例出现在广灵驴群体中，这与广灵驴目前处于保种状态是相一致的。因为处于保种状态，很少有其他驴品种与之交流，必然导致群体近交系数增大，纯合度提高，而群体多样性以及SNP和InDel的比例下降。然而，杂合SNP和InDel数目以及π ratio的最大值都出现在山东小毛驴群体中。尽管山东小毛驴各驴体高明显较德州驴小，而且毛色多为灰色，但很可能与其地理位置相近的德州驴有杂交，也可能与地缘位置较近的河北省的一些驴种有基因交流，而杂合SNP和InDel数目以及π ratio的最低值都出现在华北驴群体中，可能与华北驴距离各驴种都较远，而且也可能与采集的样本来自偏远的蒙古族村子有关，与其余驴种很少有杂交，几乎没有外源基因导入该群体。

遗传分化通常是生物长期进化的产物，受交配系统、生物历史和基因流等生物特征的影响。WRIGHT[15]认为，如果Fst<0.050 0，各群体间就几乎没有分化。本研究中，各驴群体间的Fst值介于 0.007 80～0.011 460，明显低于具有混合交配系统或几年寿命的物种的Fst值，表明各驴群体在遗传上很相似。这也凸显了对这些驴品种进行保护的必要性，不及时进行保护，现今犹存的一些群体特征会随着群体间杂交而逐渐消失。

Fsip1是一种生精细胞特异性表达蛋白，也是精子鞭毛纤维鞘的成分，在成年动物的睾丸组织有超高水平的表达，其次为中枢神经系统也有一定表达，已知参与组装AKAP4辅助调节蛋白激酶A (PKA)。LIU等[18]发现，Fsip1在HER2过表达型乳腺癌中能够与HER2直接结合调控乳腺癌细胞的生长和侵袭。随后，LIU等[19]发现，Fsip1能够通过诱导自噬，减少线粒体生成及增强和激活AMPK途径来调节三阴性乳腺癌细胞的耐药性。AHNAK2是AHNAK2基因编码的1个较大的核蛋白，该蛋白质可能通过与钙通道蛋白相结合在钙信号调节中发挥重要作用[20]。NKG2DL1是NKG2D这种免疫受体蛋白的配体，该受体因为具有抗病毒和抗肿瘤功能近些年引起了广泛关注[21]。KLK1E2基因编码的是丝氨酸蛋白酶的1个亚基，该酶具有广泛的生理功能，越来越多的证据表明其参与癌症发生，有些分子具有作为癌症和其他疾病新的生物标志物的潜力。有研究表明，KLK1在流感病毒感染早期就干预抗病毒防御，调节流感感染的严重程度，慢性阻塞性肺病患者KLK1表达降低可能导致流感恶化[22]。或许受到选择的这些基因对于山东小毛驴适应严酷环境的能力至关重要，这些基因才会在山东小毛驴群体内逐渐固定下来。这暗示山东小毛驴可能在免疫相关过程中经历了一定的选择作用，对揭示山东小毛驴免疫性状的遗传机制具有一定的意义。这也说明相对于其他驴种，山东小毛驴在抗病力强、繁殖力强这些优良性状上经历了较强的人工选择。

本研究利用5个驴群体共计60个个体的全基因组重测序数据，对山东小毛驴群体选择信号进行了检测分析，检测到选择信号的区域共计95个，落入这些区域并且选择信号强的候选基因有10个。相比于其他4个驴群体，山东小毛驴在免疫、生殖等性状上经历了人工选择。

致谢：感谢内蒙古阿鲁科尔沁旗太极天驴集团有限公司的钟勇先生和舒蕾先生在采集华北驴样本时的大力支持；感谢国家级广灵驴保种场姜正广先生和许增孝先生在采集广灵驴样本时的大力支持；感谢东阿阿胶股份有限公司嵇传良先生在采集德州驴样本时提供的大力支持；感谢山东海阳小毛驴保种繁育基地由松利先生在采集山东小毛驴样本时提供的大力支持。