全基因組選擇信號揭示綿羊毛囊發(fā)育及脫毛性狀相關(guān)的候選基因

雷志惠，張利平*，趙洪昌，朱韶華，盧曾奎，郭婷婷，孫渭博，赫 雪，岳耀敬*

(1.甘肅農(nóng)業(yè)大學(xué)動物科學(xué)技術(shù)學(xué)院，蘭州 730070； 2. 中國農(nóng)業(yè)科學(xué)院蘭州畜牧與獸藥研究所，蘭州 730050)

羊毛在養(yǎng)羊業(yè)中具有重要經(jīng)濟(jì)價值，提高羊毛的產(chǎn)量和品質(zhì)是增加毛用綿羊經(jīng)濟(jì)價值的有效途徑。隨著羊毛產(chǎn)量的增加和剪毛工藝的發(fā)展，羊毛價值與剪毛成本之間呈現(xiàn)出收支失衡的現(xiàn)象[1]，因此，人們嘗試在保留產(chǎn)毛性能的基礎(chǔ)上，通過挖掘與脫毛性狀相關(guān)的候選基因來培育無需剪毛的綿羊品種[2-3]，以降低剪毛成本并有效控制毛叢寄生蟲滋生，維持動物機(jī)體健康。但目前綿羊脫毛性狀的相關(guān)研究鮮有報(bào)道，其遺傳機(jī)制尚不清楚。

綿羊經(jīng)過長期的品種選育和品種質(zhì)量提高等過程，使其被毛分化出粗毛、細(xì)毛和半細(xì)毛等多種類型，其中，非洲杜泊羊等無絨毛型綿羊，具有明顯的脫毛性狀，能夠適應(yīng)多種氣候條件和生態(tài)環(huán)境，羊毛可隨氣候變化自然脫落，無需剪毛[4]；而高山美利奴羊等細(xì)毛型綿羊，不具有脫毛性狀，被毛基本由同質(zhì)細(xì)毛組成，長度7～10 cm，整齊均勻，剪毛量4.2～9.8 kg，是毛紡工業(yè)的優(yōu)質(zhì)原料[5]；薩?？搜虻戎忻途d羊，羊毛細(xì)度50～58支，長度7～8 cm，剪毛量2.3～6.0 kg，主要用于毛線、毛毯和工業(yè)用氈等[6]。多種羊毛類型主要受基因調(diào)控、毛囊結(jié)構(gòu)及其生長發(fā)育之間的差異等因素影響，諸多研究表明，F(xiàn)GF家族成員與毛囊的生長周期和分化有密切關(guān)系[7]：在毛囊生長時期，F(xiàn)GF2和FGF5均有調(diào)控毛囊發(fā)育的作用[8]；FGF5可抑制毛發(fā)生長，當(dāng)過表達(dá)時會縮短毛囊生長期，促使毛囊進(jìn)入退行期[9]；FGF7可以抑制毛囊細(xì)胞的分化[10]。毛囊主要分布在皮膚的表皮和真皮中，真皮中的毛乳頭細(xì)胞可誘導(dǎo)調(diào)控毛囊生長發(fā)育，在毛囊形態(tài)學(xué)發(fā)生和周期性生長調(diào)控中起主導(dǎo)作用[11]。

隨著生物技術(shù)快速發(fā)展，選擇信號分析可鑒定出選擇過程中基因組上留下的痕跡。在諸多選擇信號分析方法中，F(xiàn)st可利用多態(tài)性數(shù)據(jù)，檢測群體之間等位基因頻率差異[12]，并有效檢測群體中受正向選擇的位點(diǎn)；θπ Ratio是基于基因雜合度的選擇性消除分析的方法[13]，可通過兩個群體之間θπ的比值來確定受選擇位點(diǎn)所在的群體。Guo等[14]對6個表型不同的山羊品種進(jìn)行選擇信號分析，發(fā)現(xiàn)了影響被毛顏色(IRF4、EXOC2、RALY、EIF2S2、KITLG)、生長性狀(LDB2)和繁殖性狀(KHDRBS2)的相關(guān)基因。Li等[15]利用Fst和θπ Ratio方法對絨山羊重測序數(shù)據(jù)進(jìn)行選擇信號分析，篩選出FGF5、SGK3、IGFBP7、ROCK1等可能參與羊絨纖維形成的基因。Islam等[16]利用Fst和XP-EHH對6種不同山羊群體進(jìn)行選擇信號分析，篩選出MARF1、SYCP2、TMEM200C、SF1、ADCY1和BMP5等與山羊繁殖性狀相關(guān)的候選基因。因此，利用選擇信號分析方法挖掘綿羊毛囊發(fā)育及脫毛性狀的關(guān)聯(lián)基因并了解其功能[17]，對探究其遺傳機(jī)制具有重要意義。

本研究以3種羊毛類型(無絨毛型、細(xì)毛型和中毛型)共21個品種的綿羊群體作為研究對象，利用Illumina Ovine SNP 50K芯片分型數(shù)據(jù)，基于Fst和θπ Ratio方法進(jìn)行選擇信號分析，挖掘綿羊毛囊發(fā)育及脫毛性狀的候選基因，以期為綿羊的分子育種提供重要參考。

1 材料與方法

1.1 試驗(yàn)樣本與數(shù)據(jù)處理

本研究選用國內(nèi)外21個綿羊品種的290只個體作為研究對象，依據(jù)羊毛類型分為無絨毛型(hair)、細(xì)毛型(fine wool)和中毛型(medium wool)，各品種名稱、縮寫和樣本量信息見表1。本試驗(yàn)所使用的21個綿羊品種Illumina Ovine SNP 50K芯片的基因分型數(shù)據(jù)來源于ISGC(https://www.sheephapmap.org/)、NRSP(https://www.animalgenome.org/sheep/community/)、WIDDE(http://widde.toulouse.inra.fr/)和OSF(https://www.ontariosheep.org/)4個數(shù)據(jù)庫。為了將芯片基因分型數(shù)據(jù)與現(xiàn)有綿羊基因組注釋文件相對應(yīng)，將SNP位置信息按照Ovis_aries_v4.0參考基因組的位置重新排序。使用PLINK(v1.90b)軟件[18]進(jìn)行質(zhì)控，標(biāo)準(zhǔn)如下：1)SNP檢出率大于90%；2)次等位基因頻率MAF大于0.01；3)哈代溫伯格平衡P值大于10-6；4)僅對常染色體上的選擇信號進(jìn)行檢測[17,19]。經(jīng)過以上質(zhì)控后共得到290個綿羊個體的34 436個SNPs用于后續(xù)選擇信號分析。

表1 本研究中綿羊品種信息

1.2 試驗(yàn)方法

1.2.1 群體背景分析 為了解所選試驗(yàn)樣本的聚類情況和遺傳關(guān)系，本研究通過PLINK(v1.90b)軟件(參數(shù)：plink--file bestqc.comsnp--pca 3 header tabs)進(jìn)行主成分分析，利用R軟件包(ggplot)[20]對所得結(jié)果進(jìn)行可視化展示。

為從群體層面上研究所選群體的分層情況，本研究基于貝葉斯模型進(jìn)行STRUCTURE分析，利用ADMIXTURE軟件[21](參數(shù)：for k in 1 2 3 4; do admixture--cv admixture_bestqc.comsnp.bed $k| tee admixture_bestqc.comsnp/log${k}.out; done;)對各個位點(diǎn)的基因型數(shù)據(jù)進(jìn)行聚類模型的構(gòu)建。

1.2.2 選擇信號檢測 本研究利用Fst和θπ Ratio兩種方法開展選擇信號檢測，利用VCFtools (0.1.15)軟件[22]進(jìn)行數(shù)據(jù)分析，運(yùn)用滑動窗口的方法，以1 000 kb作為一個滑動窗口，100 kb作為步長，計(jì)算染色體窗口內(nèi)的Fst值[23-24]。提取Fst和θπ Ratio均為top 5%[12]的窗口內(nèi)位點(diǎn)作為顯著SNP位點(diǎn)，即為選擇信號的候選位點(diǎn)。

1.2.3 候選基因檢測與注釋 將Fst和θπ Ratio篩選出top 5%的SNPs作為本次試驗(yàn)的“離群位點(diǎn)”，離群位點(diǎn)上、下游各50 kb視為選擇信號作用區(qū)域，參照綿羊Ovis_aries_v4.0基因組信息對選擇信號作用區(qū)域進(jìn)行基因注釋，注釋所得到的基因稱作“候選基因”。

2 結(jié) 果

2.1 群體背景分析

對21個綿羊群體的基因分型數(shù)據(jù)進(jìn)行主成分分析和群體結(jié)構(gòu)分析，結(jié)果如圖1所示。PC1能夠解釋12.63%的遺傳變異，很明顯地將中毛型綿羊從所有群體中區(qū)分開；在PC1正值方向，細(xì)毛羊基本聚集在一起，部分無絨毛型綿羊與細(xì)毛型綿羊混在一起，PC2能夠解釋9.15%的遺傳變異，在PC2正值方向，部分無絨毛型綿羊單獨(dú)聚在一起。

圖1 不同綿羊群體的主成分分析Fig.1 Principal component analysis of different sheep populations

為證實(shí)PCA所得結(jié)果的準(zhǔn)確性，另采用ADMIXTURE軟件進(jìn)行分析。通過比較不同K值的交叉驗(yàn)證錯誤率發(fā)現(xiàn)，當(dāng)K=10時交叉驗(yàn)證的錯誤率最小(圖2)，因此CV error (K=10)是最佳的建模選擇，ADMIXTURE結(jié)果如圖3所示。綿羊各群體的親緣關(guān)系較遠(yuǎn)，起源相對復(fù)雜，不同羊毛類型的綿羊之間存在明顯的基因交流。

2.2 選擇信號分析

2.2.1Fst分析 用PLINK(v1.90b)軟件對質(zhì)控后的34 436個SNPs位點(diǎn)進(jìn)行計(jì)算，得到3個不同毛型群體之間成對的遺傳分化指數(shù)Fst值，并繪制Fine-Hair和Medium-Hair兩組中每個位點(diǎn)的Fst值經(jīng)驗(yàn)分布圖，如圖4所示。在全基因組水平取Fine-Hair組和Medium-Hair組top 5%的SNPs作為受選擇位點(diǎn)，分別有1 469個和1 049個SNPs位點(diǎn)高于閾值線。其中，在Fine-Hair組間存在14個位點(diǎn)位于Fst分布的尾端(Fst>0.25)，在Medium-Hair組間存在77個位點(diǎn)位于Fst分布的尾端(Fst>0.25)，表明在Fine-Hair組和Medium-Hair組中均存在遺傳分化。

2.2.2 θπ Ratio分析 將細(xì)毛型綿羊與中毛型綿羊作為參考群體，無絨毛型綿羊作為試驗(yàn)群體進(jìn)行θπ Ratio分析(表2)。分別取Fine-Hair組和Medium-Hair組top 5%的SNPs作為受選擇位點(diǎn)，結(jié)果均有1 221個SNPs位點(diǎn)高于閾值線。當(dāng)核苷酸多態(tài)性比率θπ Ratio越偏離1時，基因組受選擇程度越高，其中在Fine-Hair組間存在46個位點(diǎn)偏離較大(θπ Ratio>1.4)，在Medium-Hair組間存在69個位點(diǎn)偏離較大(θπ Ratio>1.4)。

圖2 綿羊群體的交叉驗(yàn)證Fig.2 Cross validation of sheep populations

2.2.3Fst和θπ Ratio的基因定位和功能注釋 利用Fst和θπ Ratio在綿羊基因組中篩選顯著基因，其中在Fine-Hair組中篩選到608個顯著基因(圖5A)；在Medium-Hair組中篩選到473個顯著基因(圖5B)。根據(jù)兩組中獲得的顯著基因繪制Venn diagram，在交叉區(qū)域內(nèi)獲得177個基因(圖6)，經(jīng)功能注釋，最終篩選出SOX18、EXT2、ACP2、ALX4、FGF1和LRP4等與毛囊發(fā)育及脫毛性狀相關(guān)的基因(表3)。

圖3 綿羊群體的STRUCTURE分析Fig.3 STRUCTURE analysis of sheep populations

圖4 試驗(yàn)群體間Fst值分布圖Fig.4 The distribution of Fst values in experimental groups

3 討 論

本研究采用群體間等位基因頻率變化的方法計(jì)算Fst，以衡量群體分化程度和遺傳距離，尋找受選擇的潛在位點(diǎn)。但基于單位點(diǎn)估計(jì)的方法會導(dǎo)致假陽性結(jié)果，為提高檢測的準(zhǔn)確性，選用滑動窗口的方法降低假陽性出現(xiàn)的概率[25]；θπ Ratio根據(jù)基因雜合度進(jìn)行分析，計(jì)算兩個群體之間的比值，基因組受選擇程度越高，θπ Ratio就越偏離1[15,26]；因此，本研究通過上述兩種方法篩選顯著SNP位點(diǎn)，選取top 5%的SNPs作為受選擇位點(diǎn)[27]，可有效提高選擇信號篩選的準(zhǔn)確性。對3種羊毛類型的綿羊群體進(jìn)行遺傳背景和主成分分析發(fā)現(xiàn)，不同品種所組成的群體，根據(jù)其產(chǎn)地和親緣關(guān)系分布在不同位置，結(jié)果與預(yù)期一致；根據(jù)群體結(jié)構(gòu)分析，研究涉及的群體祖先來源廣泛，血緣組成復(fù)雜，存在明顯的基因交流，主成分分析和群體結(jié)構(gòu)分析的結(jié)果一致，為選擇信號分析提供準(zhǔn)確依據(jù)。采用Fst和θπ Ratio兩種方法進(jìn)行相互驗(yàn)證，可提高基因定位可靠性，兩種方法的交叉區(qū)域內(nèi)共獲得177個強(qiáng)烈受選擇基因。

表2 試驗(yàn)群體間θπ 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 無絨毛羊中與毛囊發(fā)育及脫毛性狀相關(guān)的基因

以交叉區(qū)域內(nèi)的基因?yàn)榛A(chǔ)，結(jié)合前人研究結(jié)果進(jìn)行基因功能注釋，根據(jù)注釋結(jié)果可知，成纖維細(xì)胞生長因子(FGF)家族成員參與機(jī)體許多不同類型細(xì)胞的增值和分化，在調(diào)節(jié)毛囊發(fā)育中起重要作用[28-29]，F(xiàn)GF家族由許多肝素結(jié)合蛋白組成，能夠影響細(xì)胞的有絲分裂，促進(jìn)遷移和分化以及誘導(dǎo)血管生成[30]。此外，成纖維細(xì)胞生長因子1(FGF1)在毛囊中表達(dá)并可能調(diào)節(jié)其生長[31]。本研究中，在綿羊Chr5中檢測出SNP位點(diǎn)5-50965861(Fst1=0.418 6，F(xiàn)st2=0.348 0，θπ Ratio1=1.051 7，θπ Ratio2=1.537 8)，經(jīng)功能注釋，發(fā)現(xiàn)了與皮膚及毛囊發(fā)育相關(guān)的FGF1基因，該位點(diǎn)在無絨毛型綿羊群體中受到強(qiáng)烈選擇，據(jù)此推斷，該基因?qū)o絨毛型綿羊群體羊毛脫落具有重要作用；在Chr13中檢測出SNP位點(diǎn)13-53089456(Fst1=0.362 2，F(xiàn)st2=0.157 0，θπ Ratio1=1.014 6，θπ Ratio2=1.015 0)，位于SOX18(SRY (sex determining region Y) -box 18)內(nèi)，該位點(diǎn)在無絨毛型綿羊群體受到強(qiáng)烈選擇。SOX18是SOX轉(zhuǎn)錄因子F(SOX F)亞家族的重要成員之一，在人和小鼠上分別位于Chr20和Chr2，在血管、淋巴管發(fā)育過程中發(fā)揮重要作用[32]。SOX蛋白屬于SRY相關(guān)基因家族編碼的轉(zhuǎn)錄因子，在機(jī)體中參與多種組織器官發(fā)育過程與調(diào)控機(jī)制[33]。此外，該蛋白還參與血管生成、毛囊發(fā)育和淋巴管生成等生物學(xué)過程[34]。稀毛癥-淋巴水腫-毛細(xì)血管擴(kuò)張綜合征通常表現(xiàn)為淋巴水腫、眉毛稀疏、發(fā)育遲緩和脫發(fā)等癥狀[35]，Irrthum等[36]對患有稀毛癥-淋巴水腫-毛細(xì)血管擴(kuò)張疾病的臨床患者及后代進(jìn)行相關(guān)研究，發(fā)現(xiàn)轉(zhuǎn)錄因子SOX18基因在這些患者中存在隱性突變。Pennisi等[37]發(fā)現(xiàn)SOX18在小鼠胚胎發(fā)育中的血管內(nèi)皮細(xì)胞和毛囊中均有表達(dá)，此外，在兩個不同的突變小鼠等位基因中均發(fā)現(xiàn)SOX18存在點(diǎn)突變，與野生型對比的結(jié)果表明，含有上述突變的融合蛋白缺乏激活轉(zhuǎn)錄能力，綜上，SOX18的異常表達(dá)會引發(fā)毛發(fā)稀少-淋巴水腫-毛細(xì)血管擴(kuò)張綜合征，并確定SOX18是心血管和毛囊形成的關(guān)鍵基因[32]。因此推測，SOX18在綿羊毛囊形成中具有關(guān)鍵作用，SOX18的異常表達(dá)可能與綿羊毛囊的發(fā)育及脫毛性狀有關(guān)，但仍需進(jìn)一步驗(yàn)證。在綿羊Chr15的SNP位點(diǎn)15-74875572(Fst1=0.139 7，F(xiàn)st2=0.198 9，θπ Ratio1=1.305 0，θπ Ratio2=1.954 7)注釋到LRP4，發(fā)現(xiàn)該位點(diǎn)在無絨毛綿羊中受到強(qiáng)烈選擇，孫向東[38]研究發(fā)現(xiàn)，LRP4為低密度脂蛋白受體家族成員之一，編碼的蛋白質(zhì)可能是Wnt信號的調(diào)節(jié)劑，在皮膚和大腦等組織中廣泛表達(dá)[39]。Ahn等[40]的研究結(jié)果表明，在細(xì)胞培養(yǎng)過程中，LRP4過表達(dá)可抑制調(diào)控毛囊周期性發(fā)育的Wnt信號激活[41-42]；此外，BMP家族基因是毛囊發(fā)育及周期性生長關(guān)鍵因子，調(diào)控毛囊生長發(fā)育、胚胎發(fā)育及骨骼生長等[43]，研究表明LRP4可參與BMP信號調(diào)節(jié)[44]，綜上，LRP4對Wnt和BMP通路及信號轉(zhuǎn)導(dǎo)存在直接或間接抑制作用，其機(jī)理目前尚不清楚。Zhang等[45]發(fā)現(xiàn),LRP4在乳腺和毛囊等皮膚附屬物中表達(dá)，與機(jī)體毛囊發(fā)育關(guān)系密切。在綿羊Chr15中檢測出SNP位點(diǎn)15-75217088(Fst1=0.101 7，F(xiàn)st2=0.179 4，θπ Ratio1=1.155 4，θπ Ratio2=1.138 6)，位于ACP2內(nèi)，該基因突變會導(dǎo)致小鼠的小腦和皮膚畸形，引起生長遲緩和毛發(fā)延遲出現(xiàn)等現(xiàn)象[46]；在Chr15中檢測出SNP位點(diǎn)15-72418538(Fst1=0.150 4，F(xiàn)st2=0.167 2，θπ Ratio1=0.022 0，θπ Ratio2=1.035 6)注釋到ALX4，發(fā)現(xiàn)該基因突變可能與Wnt/β-catenin信號相互作用影響表皮分化，導(dǎo)致毛囊形成和分化異常，使機(jī)體出現(xiàn)毛發(fā)脫落的現(xiàn)象[47]。

研究發(fā)現(xiàn)，上述6個與綿羊毛囊發(fā)育相關(guān)的候選基因均處于Chr5、Chr13和Chr15，因此，后續(xù)研究可將以上染色體作為重點(diǎn)區(qū)域進(jìn)行深入挖掘。另外，毛囊具有自我更新和周期性生長等特點(diǎn)，是研究器官大小、組織形態(tài)、干細(xì)胞再生和內(nèi)環(huán)境穩(wěn)態(tài)等的理想模型[48]?？蔀樯钊胪诰蚓d羊的分子遺傳機(jī)制提供重要的參考依據(jù)，也為研究綿羊毛囊發(fā)育及生理性脫毛提供新思路。

4 結(jié) 論

本研究利用Fst和θπ Ratio方法對綿羊SNP芯片的分型數(shù)據(jù)進(jìn)行選擇信號檢測，對強(qiáng)烈受到選擇的SNPs位點(diǎn)進(jìn)行注釋得到177個受強(qiáng)烈選擇的基因，其中FGF1、SOX18、ALX4和LRP4與毛囊的發(fā)育周期、毛發(fā)形成以及毛囊和皮脂腺的部分細(xì)胞具有密切聯(lián)系，且主要位于Chr5、Chr13和Chr15，這些發(fā)現(xiàn)可為綿羊毛囊發(fā)育機(jī)制的研究提供參考，也可將其作為后期功能驗(yàn)證的候選基因。