全基因组选择信号揭示绵羊毛囊发育及脱毛性状相关的候选基因

雷志惠，张利平*，赵洪昌，朱韶华，卢曾奎，郭婷婷，孙渭博，赫 雪，岳耀敬*

(1.甘肃农业大学动物科学技术学院，兰州 730070； 2. 中国农业科学院兰州畜牧与兽药研究所，兰州 730050)

羊毛在养羊业中具有重要经济价值，提高羊毛的产量和品质是增加毛用绵羊经济价值的有效途径。随着羊毛产量的增加和剪毛工艺的发展，羊毛价值与剪毛成本之间呈现出收支失衡的现象[1]，因此，人们尝试在保留产毛性能的基础上，通过挖掘与脱毛性状相关的候选基因来培育无需剪毛的绵羊品种[2-3]，以降低剪毛成本并有效控制毛丛寄生虫滋生，维持动物机体健康。但目前绵羊脱毛性状的相关研究鲜有报道，其遗传机制尚不清楚。

绵羊经过长期的品种选育和品种质量提高等过程，使其被毛分化出粗毛、细毛和半细毛等多种类型，其中，非洲杜泊羊等无绒毛型绵羊，具有明显的脱毛性状，能够适应多种气候条件和生态环境，羊毛可随气候变化自然脱落，无需剪毛[4]；而高山美利奴羊等细毛型绵羊，不具有脱毛性状，被毛基本由同质细毛组成，长度7～10 cm，整齐均匀，剪毛量4.2～9.8 kg，是毛纺工业的优质原料[5]；萨福克羊等中毛型绵羊，羊毛细度50～58支，长度7～8 cm，剪毛量2.3～6.0 kg，主要用于毛线、毛毯和工业用毡等[6]。多种羊毛类型主要受基因调控、毛囊结构及其生长发育之间的差异等因素影响，诸多研究表明，FGF家族成员与毛囊的生长周期和分化有密切关系[7]：在毛囊生长时期，FGF2和FGF5均有调控毛囊发育的作用[8]；FGF5可抑制毛发生长，当过表达时会缩短毛囊生长期，促使毛囊进入退行期[9]；FGF7可以抑制毛囊细胞的分化[10]。毛囊主要分布在皮肤的表皮和真皮中，真皮中的毛乳头细胞可诱导调控毛囊生长发育，在毛囊形态学发生和周期性生长调控中起主导作用[11]。

随着生物技术快速发展，选择信号分析可鉴定出选择过程中基因组上留下的痕迹。在诸多选择信号分析方法中，Fst可利用多态性数据，检测群体之间等位基因频率差异[12]，并有效检测群体中受正向选择的位点；θπ Ratio是基于基因杂合度的选择性消除分析的方法[13]，可通过两个群体之间θπ的比值来确定受选择位点所在的群体。Guo等[14]对6个表型不同的山羊品种进行选择信号分析，发现了影响被毛颜色(IRF4、EXOC2、RALY、EIF2S2、KITLG)、生长性状(LDB2)和繁殖性状(KHDRBS2)的相关基因。Li等[15]利用Fst和θπ Ratio方法对绒山羊重测序数据进行选择信号分析，筛选出FGF5、SGK3、IGFBP7、ROCK1等可能参与羊绒纤维形成的基因。Islam等[16]利用Fst和XP-EHH对6种不同山羊群体进行选择信号分析，筛选出MARF1、SYCP2、TMEM200C、SF1、ADCY1和BMP5等与山羊繁殖性状相关的候选基因。因此，利用选择信号分析方法挖掘绵羊毛囊发育及脱毛性状的关联基因并了解其功能[17]，对探究其遗传机制具有重要意义。

本研究以3种羊毛类型(无绒毛型、细毛型和中毛型)共21个品种的绵羊群体作为研究对象，利用Illumina Ovine SNP 50K芯片分型数据，基于Fst和θπ Ratio方法进行选择信号分析，挖掘绵羊毛囊发育及脱毛性状的候选基因，以期为绵羊的分子育种提供重要参考。

1 材料与方法

1.1 试验样本与数据处理

本研究选用国内外21个绵羊品种的290只个体作为研究对象，依据羊毛类型分为无绒毛型(hair)、细毛型(fine wool)和中毛型(medium wool)，各品种名称、缩写和样本量信息见表1。本试验所使用的21个绵羊品种Illumina Ovine SNP 50K芯片的基因分型数据来源于ISGC(https://www.sheephapmap.org/)、NRSP(https://www.animalgenome.org/sheep/community/)、WIDDE(http://widde.toulouse.inra.fr/)和OSF(https://www.ontariosheep.org/)4个数据库。为了将芯片基因分型数据与现有绵羊基因组注释文件相对应，将SNP位置信息按照Ovis_aries_v4.0参考基因组的位置重新排序。使用PLINK(v1.90b)软件[18]进行质控，标准如下：1)SNP检出率大于90%；2)次等位基因频率MAF大于0.01；3)哈代温伯格平衡P值大于10-6；4)仅对常染色体上的选择信号进行检测[17,19]。经过以上质控后共得到290个绵羊个体的34 436个SNPs用于后续选择信号分析。

表1 本研究中绵羊品种信息

1.2 试验方法

1.2.1 群体背景分析 为了解所选试验样本的聚类情况和遗传关系，本研究通过PLINK(v1.90b)软件(参数：plink--file bestqc.comsnp--pca 3 header tabs)进行主成分分析，利用R软件包(ggplot)[20]对所得结果进行可视化展示。

为从群体层面上研究所选群体的分层情况，本研究基于贝叶斯模型进行STRUCTURE分析，利用ADMIXTURE软件[21](参数：for k in 1 2 3 4; do admixture--cv admixture_bestqc.comsnp.bed $k| tee admixture_bestqc.comsnp/log${k}.out; done;)对各个位点的基因型数据进行聚类模型的构建。

1.2.2 选择信号检测 本研究利用Fst和θπ Ratio两种方法开展选择信号检测，利用VCFtools (0.1.15)软件[22]进行数据分析，运用滑动窗口的方法，以1 000 kb作为一个滑动窗口，100 kb作为步长，计算染色体窗口内的Fst值[23-24]。提取Fst和θπ Ratio均为top 5%[12]的窗口内位点作为显著SNP位点，即为选择信号的候选位点。

1.2.3 候选基因检测与注释 将Fst和θπ Ratio筛选出top 5%的SNPs作为本次试验的“离群位点”，离群位点上、下游各50 kb视为选择信号作用区域，参照绵羊Ovis_aries_v4.0基因组信息对选择信号作用区域进行基因注释，注释所得到的基因称作“候选基因”。

2 结 果

2.1 群体背景分析

对21个绵羊群体的基因分型数据进行主成分分析和群体结构分析，结果如图1所示。PC1能够解释12.63%的遗传变异，很明显地将中毛型绵羊从所有群体中区分开；在PC1正值方向，细毛羊基本聚集在一起，部分无绒毛型绵羊与细毛型绵羊混在一起，PC2能够解释9.15%的遗传变异，在PC2正值方向，部分无绒毛型绵羊单独聚在一起。

图1 不同绵羊群体的主成分分析Fig.1 Principal component analysis of different sheep populations

为证实PCA所得结果的准确性，另采用ADMIXTURE软件进行分析。通过比较不同K值的交叉验证错误率发现，当K=10时交叉验证的错误率最小(图2)，因此CV error (K=10)是最佳的建模选择，ADMIXTURE结果如图3所示。绵羊各群体的亲缘关系较远，起源相对复杂，不同羊毛类型的绵羊之间存在明显的基因交流。

2.2 选择信号分析

2.2.1Fst分析 用PLINK(v1.90b)软件对质控后的34 436个SNPs位点进行计算，得到3个不同毛型群体之间成对的遗传分化指数Fst值，并绘制Fine-Hair和Medium-Hair两组中每个位点的Fst值经验分布图，如图4所示。在全基因组水平取Fine-Hair组和Medium-Hair组top 5%的SNPs作为受选择位点，分别有1 469个和1 049个SNPs位点高于阈值线。其中，在Fine-Hair组间存在14个位点位于Fst分布的尾端(Fst>0.25)，在Medium-Hair组间存在77个位点位于Fst分布的尾端(Fst>0.25)，表明在Fine-Hair组和Medium-Hair组中均存在遗传分化。

2.2.2 θπ Ratio分析 将细毛型绵羊与中毛型绵羊作为参考群体，无绒毛型绵羊作为试验群体进行θπ Ratio分析(表2)。分别取Fine-Hair组和Medium-Hair组top 5%的SNPs作为受选择位点，结果均有1 221个SNPs位点高于阈值线。当核苷酸多态性比率θπ Ratio越偏离1时，基因组受选择程度越高，其中在Fine-Hair组间存在46个位点偏离较大(θπ Ratio>1.4)，在Medium-Hair组间存在69个位点偏离较大(θπ Ratio>1.4)。

图2 绵羊群体的交叉验证Fig.2 Cross validation of sheep populations

2.2.3Fst和θπ Ratio的基因定位和功能注释 利用Fst和θπ Ratio在绵羊基因组中筛选显著基因，其中在Fine-Hair组中筛选到608个显著基因(图5A)；在Medium-Hair组中筛选到473个显著基因(图5B)。根据两组中获得的显著基因绘制Venn diagram，在交叉区域内获得177个基因(图6)，经功能注释，最终筛选出SOX18、EXT2、ACP2、ALX4、FGF1和LRP4等与毛囊发育及脱毛性状相关的基因(表3)。

图3 绵羊群体的STRUCTURE分析Fig.3 STRUCTURE analysis of sheep populations

图4 试验群体间Fst值分布图Fig.4 The distribution of Fst values in experimental groups

3 讨 论

本研究采用群体间等位基因频率变化的方法计算Fst，以衡量群体分化程度和遗传距离，寻找受选择的潜在位点。但基于单位点估计的方法会导致假阳性结果，为提高检测的准确性，选用滑动窗口的方法降低假阳性出现的概率[25]；θπ Ratio根据基因杂合度进行分析，计算两个群体之间的比值，基因组受选择程度越高，θπ Ratio就越偏离1[15,26]；因此，本研究通过上述两种方法筛选显著SNP位点，选取top 5%的SNPs作为受选择位点[27]，可有效提高选择信号筛选的准确性。对3种羊毛类型的绵羊群体进行遗传背景和主成分分析发现，不同品种所组成的群体，根据其产地和亲缘关系分布在不同位置，结果与预期一致；根据群体结构分析，研究涉及的群体祖先来源广泛，血缘组成复杂，存在明显的基因交流，主成分分析和群体结构分析的结果一致，为选择信号分析提供准确依据。采用Fst和θπ Ratio两种方法进行相互验证，可提高基因定位可靠性，两种方法的交叉区域内共获得177个强烈受选择基因。

表2 试验群体间θπ Ratio 值的分布

A. Fine-Hair的选择信号分布；B. Medium-Hair的选择信号分布A. The selection signal distribution of Fine-Hair; B. The selection signal distribution of Medium-Hair图5 绵羊常染色体选择信号分布分析Fig.5 Analysis of selection signal distribution on sheep autosomes

图6 Fine-Hair与Medium-Hair两组中的重叠基因Fig.6 Overlapping genes of Fine-Hair and Medium-Hair

表3 无绒毛羊中与毛囊发育及脱毛性状相关的基因

以交叉区域内的基因为基础，结合前人研究结果进行基因功能注释，根据注释结果可知，成纤维细胞生长因子(FGF)家族成员参与机体许多不同类型细胞的增值和分化，在调节毛囊发育中起重要作用[28-29]，FGF家族由许多肝素结合蛋白组成，能够影响细胞的有丝分裂，促进迁移和分化以及诱导血管生成[30]。此外，成纤维细胞生长因子1(FGF1)在毛囊中表达并可能调节其生长[31]。本研究中，在绵羊Chr5中检测出SNP位点5-50965861(Fst1=0.418 6，Fst2=0.348 0，θπ Ratio1=1.051 7，θπ Ratio2=1.537 8)，经功能注释，发现了与皮肤及毛囊发育相关的FGF1基因，该位点在无绒毛型绵羊群体中受到强烈选择，据此推断，该基因对无绒毛型绵羊群体羊毛脱落具有重要作用；在Chr13中检测出SNP位点13-53089456(Fst1=0.362 2，Fst2=0.157 0，θπ Ratio1=1.014 6，θπ Ratio2=1.015 0)，位于SOX18(SRY (sex determining region Y) -box 18)内，该位点在无绒毛型绵羊群体受到强烈选择。SOX18是SOX转录因子F(SOX F)亚家族的重要成员之一，在人和小鼠上分别位于Chr20和Chr2，在血管、淋巴管发育过程中发挥重要作用[32]。SOX蛋白属于SRY相关基因家族编码的转录因子，在机体中参与多种组织器官发育过程与调控机制[33]。此外，该蛋白还参与血管生成、毛囊发育和淋巴管生成等生物学过程[34]。稀毛症-淋巴水肿-毛细血管扩张综合征通常表现为淋巴水肿、眉毛稀疏、发育迟缓和脱发等症状[35]，Irrthum等[36]对患有稀毛症-淋巴水肿-毛细血管扩张疾病的临床患者及后代进行相关研究，发现转录因子SOX18基因在这些患者中存在隐性突变。Pennisi等[37]发现SOX18在小鼠胚胎发育中的血管内皮细胞和毛囊中均有表达，此外，在两个不同的突变小鼠等位基因中均发现SOX18存在点突变，与野生型对比的结果表明，含有上述突变的融合蛋白缺乏激活转录能力，综上，SOX18的异常表达会引发毛发稀少-淋巴水肿-毛细血管扩张综合征，并确定SOX18是心血管和毛囊形成的关键基因[32]。因此推测，SOX18在绵羊毛囊形成中具有关键作用，SOX18的异常表达可能与绵羊毛囊的发育及脱毛性状有关，但仍需进一步验证。在绵羊Chr15的SNP位点15-74875572(Fst1=0.139 7，Fst2=0.198 9，θπ Ratio1=1.305 0，θπ Ratio2=1.954 7)注释到LRP4，发现该位点在无绒毛绵羊中受到强烈选择，孙向东[38]研究发现，LRP4为低密度脂蛋白受体家族成员之一，编码的蛋白质可能是Wnt信号的调节剂，在皮肤和大脑等组织中广泛表达[39]。Ahn等[40]的研究结果表明，在细胞培养过程中，LRP4过表达可抑制调控毛囊周期性发育的Wnt信号激活[41-42]；此外，BMP家族基因是毛囊发育及周期性生长关键因子，调控毛囊生长发育、胚胎发育及骨骼生长等[43]，研究表明LRP4可参与BMP信号调节[44]，综上，LRP4对Wnt和BMP通路及信号转导存在直接或间接抑制作用，其机理目前尚不清楚。Zhang等[45]发现,LRP4在乳腺和毛囊等皮肤附属物中表达，与机体毛囊发育关系密切。在绵羊Chr15中检测出SNP位点15-75217088(Fst1=0.101 7，Fst2=0.179 4，θπ Ratio1=1.155 4，θπ Ratio2=1.138 6)，位于ACP2内，该基因突变会导致小鼠的小脑和皮肤畸形，引起生长迟缓和毛发延迟出现等现象[46]；在Chr15中检测出SNP位点15-72418538(Fst1=0.150 4，Fst2=0.167 2，θπ Ratio1=0.022 0，θπ Ratio2=1.035 6)注释到ALX4，发现该基因突变可能与Wnt/β-catenin信号相互作用影响表皮分化，导致毛囊形成和分化异常，使机体出现毛发脱落的现象[47]。

研究发现，上述6个与绵羊毛囊发育相关的候选基因均处于Chr5、Chr13和Chr15，因此，后续研究可将以上染色体作为重点区域进行深入挖掘。另外，毛囊具有自我更新和周期性生长等特点，是研究器官大小、组织形态、干细胞再生和内环境稳态等的理想模型[48]。可为深入挖掘绵羊的分子遗传机制提供重要的参考依据，也为研究绵羊毛囊发育及生理性脱毛提供新思路。

4 结 论

本研究利用Fst和θπ Ratio方法对绵羊SNP芯片的分型数据进行选择信号检测，对强烈受到选择的SNPs位点进行注释得到177个受强烈选择的基因，其中FGF1、SOX18、ALX4和LRP4与毛囊的发育周期、毛发形成以及毛囊和皮脂腺的部分细胞具有密切联系，且主要位于Chr5、Chr13和Chr15，这些发现可为绵羊毛囊发育机制的研究提供参考，也可将其作为后期功能验证的候选基因。