利用可變窗口FST方法檢測不同尾型呼倫貝爾羊尾部脂肪沉積相關(guān)基因

張統(tǒng)雨，樊紅櫻，朱才業(yè)，劉家鑫，鄧天宇，杜立新，王立賢，趙福平

(中國農(nóng)業(yè)科學(xué)院北京畜牧獸醫(yī)研究所，北京 100193)

中國地方綿羊品種根據(jù)尾型分為5類：長脂尾、短脂尾、短瘦尾、長瘦尾和脂臀尾[1]。將脂肪沉積在尾部這是綿羊獨(dú)有的特性。尾部脂肪在逆境條件下(如寒冷等)提供身體能量所需，類似于駱駝的駝峰，屬于適應(yīng)性性狀[2]。但是隨著肥胖、心腦血管以及糖尿病的增加，人們更傾向于高蛋白和低脂肪的肉類[3]。另外，從飼養(yǎng)成本來看，生產(chǎn)1 kg脂肪所消耗的飼料可產(chǎn)生2 kg瘦肉[4]，因此，尾脂沉積過多不利于產(chǎn)生經(jīng)濟(jì)效益。

不同尾型的綿羊品種是通過長期的人工選擇和自然選擇形成的，而選擇也會在基因組上留下烙印，這就是選擇信號(selection signature，SS)[5]。目前，選擇信號檢測方法很多[6]，其中對兩個(gè)及多個(gè)群體進(jìn)行選擇信號檢測時(shí)比較經(jīng)典的方法就是FST[7]，該方法也常應(yīng)用于羊的選擇信號檢測中[3, 8-13]。在現(xiàn)有文獻(xiàn)資料中，對肥瘦尾性狀之間的選擇信號檢測時(shí)[3, 11-13]，都是將不同品種的羊簡單的按照尾型劃分為肥、瘦尾兩類，這樣就很難消除不同遺傳背景對檢測結(jié)果的影響。其實(shí)在中國地方品種中，呼倫貝爾羊(屬于短脂尾綿羊品種)群體內(nèi)存在明顯的尾部大小差異，主要是因?yàn)樵撈贩N內(nèi)存在兩個(gè)品系[14]，這為同一個(gè)品種內(nèi)不同品系開展選擇信號檢測提供了一個(gè)理想的動物資源。

本研究使用 Illumina Ovine 600K BeadChip 高密度SNP芯片，針對呼倫貝爾羊不同品系利用可變窗口FST進(jìn)行選擇信號檢測。旨在挖掘呼倫貝爾羊尾部脂肪沉積的候選基因，揭示影響綿羊尾部脂肪沉積基因，為培育短尾綿羊品種提供理論依據(jù)。

1 材料與方法

1.1 試驗(yàn)材料

在內(nèi)蒙古呼倫貝爾草原的牧場中，從6月齡的呼倫貝爾羊群隨機(jī)挑選2 000只綿羊，然后選出大尾和小尾各150只，并對活體綿羊進(jìn)行尾部長度、尾部最寬處的長度和周長及體尺數(shù)據(jù)進(jìn)行測量，對所選個(gè)體查找系譜信息，使用一次性含有EDTA抗凝的真空采血管進(jìn)行頸靜脈采血，-20 ℃冷凍保存，用于基因組DNA提取。

1.2 試驗(yàn)方法

1.2.1 基因型檢測 首先利用血液基因組DNA提取試劑盒(天根生化科技有限公司)提取血液樣品中的基因組DNA；定制Illumina Ovine SNP600K芯片，獲得601 715個(gè)SNP分型數(shù)據(jù)；通過GenomeStudio軟件進(jìn)行數(shù)據(jù)刪選以及SNP等位基因型頻率的分析，利用PLINK軟件[15]進(jìn)行質(zhì)量控制。質(zhì)控的標(biāo)準(zhǔn)是：平均檢出率大于0.95；平均最小等位基因頻率大于0.02；剔除未定位的和性染色體上的SNPs位點(diǎn)。通過質(zhì)控后，剩余577 401個(gè)SNPs，大尾羊142只，小尾羊146只。再利用Beagle軟件[16]對缺失數(shù)據(jù)進(jìn)行填充，用于后續(xù)的分析。

1.2.2 選擇信號檢測方法 本研究采用群體分化系數(shù)FST法進(jìn)行呼倫貝爾羊品種內(nèi)不同品系間選擇信號檢測。群體分化指數(shù)FST可用來度量群體分化程度，進(jìn)而估計(jì)品系間的分化程度。本研究利用PLINK軟件[15]計(jì)算每個(gè)SNP位點(diǎn)的FST值，其計(jì)算基本原理是按照Weir和Cockerham[17]描述的無偏估計(jì)方法，其計(jì)算公式：

(1)

其中，MSG為檢測的群體內(nèi)部位點(diǎn)的誤差均方，其計(jì)算公式：

(2)

MSP是檢測的群體之間位點(diǎn)的均方差，其計(jì)算公式：

(3)

nC指校正后群體間的平均樣本大小，其計(jì)算公式：

(4)

通過公式(1)計(jì)算出每個(gè)位點(diǎn)的FST值。再根據(jù)每條染色體上每個(gè)SNP的FST值大小利用三次光滑樣條估計(jì)方法確定窗口大小和數(shù)量[18]，主要利用R語言中GenWin包完成。然后通過整個(gè)基因組范圍內(nèi)的每個(gè)窗口平均FST值和SNP個(gè)數(shù)，構(gòu)建全基因組范圍的W統(tǒng)計(jì)量：

(5)

1.2.3 候選基因的檢測和注釋 參照NCBI數(shù)據(jù)庫(https://www.ncbi.nlm.nih.gov/)和CSIRO數(shù)據(jù)庫(https://www.livestockgenomics.csiro.au/sheep/oar3.1.php)的Ovis aries 3.1 基因組信息，對篩選出來的受選擇位點(diǎn)進(jìn)行基因注釋。以選擇信號發(fā)生區(qū)域核心SNP為中心，上下游各擴(kuò)展1 000 kb 為選擇區(qū)段。將落在這個(gè)選擇區(qū)段內(nèi)的基因定義為選擇信號的“候選基因”。

1.2.4 生物信息學(xué)分析 利用Gorilla[19](http://cbl-gorilla.cs.technion.ac.il/)進(jìn)行分析，包括細(xì)胞組分(cellular component)、分子功能(molecular function)和生物學(xué)過程(biological process)分析。

本研究采用R語言[20](https://www.r-project.org/)進(jìn)行數(shù)據(jù)處理、分析以及試驗(yàn)結(jié)果的整理。

2 結(jié) 果

2.1 呼倫貝爾羊全基因組范圍內(nèi)可變窗口分布

通過三次光滑樣條估計(jì)方法確定基因組范圍內(nèi)23 144 個(gè)可變窗口。窗口包含的SNPs個(gè)數(shù)從1到775個(gè)，其中包含最多數(shù)目的SNP窗口位于chr12: 35 241 750～43 798 950 bp處，該窗口大小為8.56 Mb。圖1展示了窗口內(nèi)SNP頻數(shù)直方，其中包括≤30個(gè)SNPs的窗口數(shù)所占總窗口數(shù)的97.7%。

圖1 窗口內(nèi)SNP頻數(shù)直方圖Fig.1 Histogram of the SNPs frequency within windows

2.2 候選基因的確定

呼倫貝爾羊全基因組范圍內(nèi)的W統(tǒng)計(jì)量分布情況見圖2，以W統(tǒng)計(jì)量為11.02作為閾值，即P<0.001，共檢測出27個(gè)顯著的基因區(qū)段。按照W統(tǒng)計(jì)量的大小列在表1中，主要分布在1、2、3、6、7、11、14、15、20、21、22和24號染色體上，最顯著的區(qū)段位于26號染色體上。

通過基因注釋，在27個(gè)顯著區(qū)段內(nèi)共找到了337個(gè)候選基因(表1)。其中最顯著的區(qū)段Chr26: 0.00～0.11，找到了兩個(gè)候選基因：ERICHI、DLGAP2。ERICHI基因功能目前報(bào)道的比較少，而DLGAP2報(bào)道與神經(jīng)元突觸功能有關(guān)，是精神分裂癥的易感基因[21]。通過數(shù)據(jù)庫查找發(fā)現(xiàn)與脂代謝相關(guān)的基因只有22個(gè)(表2)。

2.3 基因富集分析

本研究通過在線工具Gorilla(http://cbl-gorilla.cs.technion.ac.il/)，將FDR閾值設(shè)為0.001，對獲得的337個(gè)候選基因進(jìn)行基因的GO富集分析。分析結(jié)果見表3，從表3可以看出受到選擇的基因主要集中在小分子代謝過程、有機(jī)氮化合物代謝過程等生物學(xué)過程。在細(xì)胞組分上主要集中在細(xì)胞內(nèi)組成成分上，總共發(fā)現(xiàn)了197個(gè)基因，占總基因的58.46%，在所有富集條目中數(shù)量占比最多。而分子功能富集分析沒有發(fā)現(xiàn)受選擇基因。

圖2 呼倫貝爾羊全基因組W統(tǒng)計(jì)量值分布Fig.2 Genome-wide distribution of W-statistic values in Hulun Buir sheep

表1W統(tǒng)計(jì)量篩選的選擇區(qū)段及候選基因

Table1SelectedregionsidentifiedbyW-statisticandannotatedcandidategenes

染色體Chromosome物理位置/Mb Physical distanceSNP數(shù)目No. of SNPW統(tǒng)計(jì)量W-statistic基因名稱Gene name260.00～0.11222.14ERICHI, DLGAP22111.71～112.14521.99PALLD, CBR4, ANXA10, LGSN, OCA2, HERC2, NIPA1, NIPA2, CYFIP1, TUBGCP5394.05～94.212119.06ZNF638, DYSF, CYP26B1, EXOC6B, SPR, SFXN52033.55～33.891916.52PRP2, UBE2N, PRL, HDGFL11106.54～106.64817.60ETV3, FCRL5, FCRL4, FCRL3, FCRL2, FCRL1, KIRREL1, CD1D, CD1E, OR1OT2619.95～20.05916.53TBCK, NPNT, GSTCD, INTS12, ARHGEF38, PPA2, TET21516.04～16.13614.88PIWIL4, FUT4, C11ORF97, CWF19L2, GUCY1A2, ALKBHB, ELMOD1, SLN, SLC35F2, RAB39A539.54～39.781215.22OR6F1, OR13G1, OR2W3, TRIM58, OR2AK2, OR2AJ1, OR2L13, OR2M4, OR2T27, OR2G6, LYPD8, ZNF672, ZNF692, SH3BP5L, PGBD2, PLPP2, MIER2, CDC34, HCN2, BSG, TPGS1, GZMM, MADCAM1

(轉(zhuǎn)下頁 Carried forward)

(轉(zhuǎn)下頁 Carried forward)

表2篩選出已報(bào)道與脂肪合成代謝相關(guān)的候選基因

Table2Thecandidategenesrelatedtolipidmetabolisminpreviousstudies

染色體Chromosome物理位置/MbPhysical distance基因Gene基因名稱或功能Gene name or function1140.59～140.69NRIP1/RIP140核酸受體交互蛋白1 nuclear receptor interacting protein 1/受體交互蛋白140 receptor interacting protein 1401105.54～107.64CD1DCD1d分子 CD1d molecule1260.41～262.55ABCG1三磷酸結(jié)合盒轉(zhuǎn)運(yùn)體G1 ATP binding cassette subfamily G member 12192.77～194.92TMEFF2表皮生長因子(EGF)類似結(jié)構(gòu)域和兩個(gè)卵泡抑素類似結(jié)構(gòu)域組成的跨膜蛋白 transmembrane protein with EGF like and two follistatin like domains 22215.96-218.12IGFBP2胰島素樣生長因子結(jié)合蛋白2insulin like growth factor binding protein 2393.05～95.21ZNF638鋅指蛋白638 zinc finger protein 638

(轉(zhuǎn)下頁 Carried forward)

表3選擇區(qū)段內(nèi)基因GO富集條目匯總

Table3SummaryofenrichmentanalysisofGOtermsforgenesinselectedregions

GO 功能分類GO function term基因數(shù)目Gene numberP值P-valueFDRGO生物學(xué)過程 GO biological processGO:0007608-感覺嗅覺 GO:0007608-sensory perception of smell75.26×10-61.73×10-2GO:0007606-化學(xué)刺激的感覺知覺 GO:0007606-sensory perception of chemical stimulus75.26×10-68.65×10-3GO:0007600-感覺知覺 GO:0007600-sensory perception72.19×10-42.40×10-1GO:0019748-次級代謝過程 GO:0019748-secondary metabolic process23.21×10-42.64×10-1GO:0044281-小分子代謝過程-GO:0044281-small molecule metabolicprocess365.09×10-43.35×10-1GO:1901564-有機(jī)氮化合物代謝過程GO:1901564-organonitrogen compound metabolic process755.46×10-42.99×10-1GO細(xì)胞組分 GO cellular componentGO:0044424-細(xì)胞內(nèi)組成成分 GO:0044424-intracellular part1973.02×10-41.56×10-1GO:0030027-板狀偽足 GO:0030027-lamellipodium35.61×10-41.45×10-1GO:0005769-早期內(nèi)體 GO:0005769-early endosome48.91×10-41.53×10-1

3 討 論

3.1 利用可變窗口FST值進(jìn)行選擇信號檢測的優(yōu)勢

因?yàn)橄嗤旧w上的SNPs之間受到連鎖不平衡和重組率等因素的影響，在選擇信號檢測時(shí)，經(jīng)常將基因組分成等距離區(qū)間，然后計(jì)算每個(gè)區(qū)間內(nèi)所有單點(diǎn)FST的平均值，將其作為整個(gè)區(qū)間的值。一般區(qū)間劃分方法有兩種：一種就是固定窗口，另外一種就是滑動窗口。固定窗口就是將整個(gè)基因組分成相同片段的區(qū)間，而滑動窗口就是相鄰窗口之間有一定的重合。但是，第一種方法對窗口大小確定具有一定的隨機(jī)性和主觀性，而第二種方法使得相鄰區(qū)間統(tǒng)計(jì)值之間又具有很強(qiáng)的相關(guān)性，因此，這兩種方法都會降低檢測功效。于是Beissinger 等[18]提出三次光滑樣條估計(jì)方法，為確定窗口大小和數(shù)量的新方法。該方法不僅可以消除噪音的干擾，又不會減低檢測信號的功效。這也是本研究采用可變窗口FST值的主要原因。

3.2 呼倫貝爾羊品系間差異分析

呼倫貝爾羊由巴爾虎品系(半橢圓狀尾)和短尾品系(小桃狀尾)組成，主要分布在內(nèi)蒙古的呼倫貝爾牧區(qū)草原上[14]。短尾品系的形成可能與引入了前蘇聯(lián)的后貝加爾羊血液有關(guān)[22]。目前兩個(gè)品系除了在尾型上存在很大差異外，其他外貌特征基本上一致[14]。為了能更好的區(qū)分這兩個(gè)品系，本研究對尾部脂肪重量測量后選擇兩端重量的個(gè)體進(jìn)行分析，使得檢測結(jié)果更加準(zhǔn)確。

3.3 與前期研究結(jié)果的比較

與前期研究[3, 11-13]相比，本研究采用了不同的群體、試驗(yàn)設(shè)計(jì)、檢測方法和極顯著的判斷標(biāo)準(zhǔn)(P<0.001)，所以并沒有找到與前期研究完全重合的基因區(qū)段，但也找到了與前期研究臨近的顯著區(qū)段。Chr3: 94.05～94.21 Mb與Yuan等[11]研究所檢測的區(qū)段Chr3: 93.33～93.98 Mb相距0.07 Mb；Chr22: 40.92～41.01 Mb與劉真等[3]在顯著區(qū)段Chr22: 40.23～40.73 Mb相距0.19 Mb。此外，在Chr1: 140.59～140.69 Mb區(qū)段也驗(yàn)證了Xu等[23]通過GWAS找到的與綿羊尾部脂肪沉積相關(guān)的候選基因NRIP1(又叫RIP140)。因此，這些區(qū)段所對應(yīng)的基因或者鑒定的基因可作為控制綿羊尾脂重要的候選基因。

3.4 脂肪沉積相關(guān)基因功能討論

本研究檢測到了337個(gè)候選基因，有22個(gè)基因在已報(bào)道文獻(xiàn)中與脂肪代謝相關(guān)。其中與豬脂肪代謝相關(guān)的基因有IGFBP2、AKT2、LIPE、MDH2和NPNT。IGFBP2可在脂肪組織中表達(dá)，與豬眼肌面積、肥肉率、瘦肥比、花油和色值等胴體和肉質(zhì)性狀顯著相關(guān)[24-26]。AKT2可以調(diào)節(jié)肌肉與脂肪細(xì)胞內(nèi)的葡萄糖轉(zhuǎn)運(yùn)[27]，也與豬脂肪沉積相關(guān)[28]。LIPE是一種脂肪酶，可影響金華豬胴體脂肪的合成[29]。MDH2在三羧酸循環(huán)中可催化蘋果酸生成丙酮酸，并產(chǎn)生還原型輔酶Ⅱ(NADPH)，而NADPH是動物體內(nèi)脂肪酸合成及碳鏈延長過程中的重要輔酶和脂肪酸合成的供氫體，體內(nèi)NADPH的供應(yīng)情況直接影響脂肪及類脂的合成。另外，Zhou等[30]研究發(fā)現(xiàn)，MDH2基因在榮昌豬不同部位脂肪組織中存在差異化表達(dá)。NPNT基因可以有效降低肉牛肌內(nèi)脂肪含量和大理石花紋的密度[31]。由PRLR介導(dǎo)的PRL可直接調(diào)節(jié)脂肪組織中的代謝過程[32-33]。解偶聯(lián)蛋白-2(UCP2)是線粒體轉(zhuǎn)運(yùn)蛋白超家族的成員，在脂肪組織中高度表達(dá)，而脂肪酸進(jìn)入脂肪組織的總量又可調(diào)節(jié)UCP2基因的表達(dá)[34-36]。另外，褐色脂肪組織UCP1具有明顯的解偶聯(lián)活性，對于維持小型哺乳動物的體溫至關(guān)重要[37]。NRIP1基因可作為調(diào)控脂粒和葡萄糖代謝的重要因子[38-40]。ZNF638是脂肪生成的新型和早期調(diào)節(jié)劑[41]。TMEFF2基因在白色脂肪組織中高度表達(dá)[42]。ACOX1是脂肪細(xì)胞內(nèi)參與脂肪酸氧化相關(guān)的酶[43]。DDB1直接與SIRT7結(jié)合，從而抑制E3泛素連接酶復(fù)合物的形成，最終導(dǎo)致參與脂肪代謝的核受體TR4含量增加[44]。DDB1還參與組成WDTC1的一個(gè)結(jié)合元件，而WDTC1可調(diào)節(jié)脂肪相關(guān)基因的轉(zhuǎn)錄，其功能的喪失會導(dǎo)致脂肪細(xì)胞增加，從而引起脂肪沉積[45]。RYR3基因是影響中國人群循環(huán)脂聯(lián)素的QTL，而脂聯(lián)素是脂肪細(xì)胞分泌的細(xì)胞因子，可影響脂代謝和糖代謝[46]。

與奶牛乳脂相關(guān)的基因有FADS1、FADS2和ACADSB。脂肪酸脫氫酶1(FADS1)和脂肪酸脫氫酶2(FADS2)是多不飽和脂肪酸合成途徑的關(guān)鍵酶。FADS1基因與奶牛乳脂率[47]和雞肉脂肪含量[48]密切相關(guān)。FADS2中的rs174575基因多態(tài)性與乳汁ALA水平相關(guān)聯(lián)[49]。而且FADS1/2基因的多態(tài)性與脂肪酸代謝和脂肪組織炎癥相關(guān)，最終會導(dǎo)致體重減輕[50-51]。短支鏈?；o酶A脫氫酶(ACADSB)在乳脂代謝中起到重要作用，對細(xì)胞中脂肪含量和甘油三酯代謝起到調(diào)控作用[52]。

另外還有一些基因與脂肪代謝相關(guān)的疾病有關(guān)。CD1D、ABCG1、SAP30BP、RAB37和ARAP1與2型糖尿病相關(guān)。2型糖尿病患者因?yàn)橐葝u素功能不足，使糖代謝功能減弱，脂肪沉積增加，最終導(dǎo)致肥胖[53]。脂肪細(xì)胞CD1D在脂肪iNKT細(xì)胞的刺激中起關(guān)鍵作用，可導(dǎo)致高脂飲食小鼠抗炎反應(yīng)[54]，也可破壞肥胖小鼠內(nèi)臟脂肪組織中的免疫平衡，最終引起2型糖尿病的發(fā)生[55]。三磷酸結(jié)合盒轉(zhuǎn)運(yùn)體G1(ATP-binding cassette G1，ABCG1)在膽固醇逆向轉(zhuǎn)運(yùn)過程中發(fā)揮著重要作用，并參與脂肪代謝[56]，ABCG1與動脈粥樣硬化[57]、肥胖[58-59]、糖尿病[60]等代謝疾病密切相關(guān)。Chen等[61]通過單個(gè)組織中表達(dá)數(shù)量性狀位點(diǎn)(eQTL)鑒定出SAP30BP是肥胖和2型糖尿病的候選基因。RAB37基因甲基化之后容易造成成糖化血紅蛋白(HbA1c)生成，而糖化血紅蛋白與2型糖尿病、肥胖等疾病有關(guān)，RAB37基因間接影響著脂肪組織的代謝[62]。ARAP1表達(dá)量增加會導(dǎo)致胰島素功能降低，從而導(dǎo)致患2型糖尿病的風(fēng)險(xiǎn)增加[63]。

4 結(jié) 論

本研究利用高密度SNP綿羊芯片，針對呼倫貝爾羊品種內(nèi)不同品系開展全基因組選擇信號檢測，并利用最新提出的三次光滑樣條估計(jì)方法確定窗口大小和數(shù)量，共檢測出27個(gè)受選擇區(qū)段，其中3個(gè)區(qū)段驗(yàn)證了已有的報(bào)道。這些區(qū)段注釋了337個(gè)候選基因，其中有22個(gè)候選基因在相關(guān)文獻(xiàn)中已報(bào)道與脂肪代謝相關(guān)。這些候選基因主要富集在細(xì)胞內(nèi)組成成分、有機(jī)氮化合物代謝過程以及小分子代謝過程等條目中。這為綿羊尾型選育、品種改良以及后續(xù)功能基因驗(yàn)證提供了重要依據(jù)。