鐘 杰，李 利，宋天增，張紅平，郭家中

(1.四川農(nóng)業(yè)大學 動物科技學院，成都 611130；2.西藏自治區(qū)農(nóng)牧科學院 畜牧獸醫(yī)研究所，拉薩 850009)

作為山羊的一個品種標準性狀，耳型主要包括耳朵(耳廓)大小、形狀和直立性。根據(jù)耳型，現(xiàn)有山羊品種可簡單地劃分為直立耳和長垂耳兩類；其中，大多數(shù)山羊品種的耳型為直立耳，例如，薩能山羊、吐根堡山羊、濟寧青山羊等[1]。而著名的長垂耳類型山羊品種則包括波爾山羊和努比亞山羊等[2]。另外，原產(chǎn)于西班牙，現(xiàn)主要飼養(yǎng)于美國的拉曼查山羊，其耳朵極??；但該品種其他性狀正常，尤其是奶用性能突出[3]。總的來說，哺乳動物的耳朵發(fā)育是一個受到多基因協(xié)同調(diào)控的復雜過程。已有研究表明，成纖維細胞生長因子基因家族(Fibroblast growth factor，F(xiàn)GF)的一些成員是調(diào)控哺乳動物耳朵發(fā)育的關(guān)鍵信號[4-6]。例如，F(xiàn)GF10通過調(diào)控Wnt、BMP信號通路，參與到小鼠內(nèi)耳的半規(guī)管、前庭和耳蝸發(fā)育[7]。另外，同源框基因家族(Homeobox，HOX)的一些基因也在哺乳動物耳朵發(fā)育的過程中發(fā)揮重要作用。Kaur等[8]利用轉(zhuǎn)基因小鼠模型發(fā)現(xiàn)，利用肌動蛋白的啟動子特異性干擾HOX-2.2基因表達，會導致小鼠耳朵發(fā)育為畸形的小耳。在畜禽上，Qiao等[9]在二花臉和沙子嶺豬的雜交后代中，發(fā)現(xiàn)HOXA1基因的突變(c.451delinsT>C)與豬的外耳畸形顯著地關(guān)聯(lián)。最近，Wei等[10]通過全基因組分析揭示出MSRB3基因與綿羊耳朵大小相關(guān)；Kumar等[11]利用SNP芯片在7個巴基斯坦山羊品種中，也發(fā)現(xiàn)MSRB3基因與山羊耳朵大小相關(guān)聯(lián)。但總體而言，目前對山羊耳型變異的遺傳機制研究仍較少。

藏山羊養(yǎng)殖是藏區(qū)畜牧業(yè)的重要組成部分；依據(jù)生態(tài)類型，藏山羊可分為高原型與山谷型[12]。傳統(tǒng)上，由于地理和交通原因，不同地區(qū)藏山羊群體間基因交流并不頻繁。Li等[13]利用26個微衛(wèi)星標記對12個藏山羊群體進行遺傳研究，發(fā)現(xiàn)來自不同生態(tài)區(qū)的群體間存在明顯的遺傳分化。近年來，在中國山羊遺傳改良工作中，通過引入其他地區(qū)的優(yōu)良品種與地方品種雜交，取得了較好的效果[14-15]。在藏山羊遺傳改良工作中，其他山羊品種被引入藏區(qū)[16-17]，導致一些地區(qū)的藏山羊群體遺傳結(jié)構(gòu)進一步發(fā)生變化。上述研究表明，青藏高原地區(qū)的藏山羊群體結(jié)構(gòu)值得深入研究。

1 材料與方法

1.1 試驗材料

在西藏自治區(qū)仲巴縣帕江鄉(xiāng)某羊場分別收集來自同一群體的正常耳型(n=30)和小耳型(n=30)藏山羊的耳長、耳寬等表型數(shù)據(jù)。用φ=75%的酒精對羊耳進行消毒處理后，使用耳號鉗采集耳組織樣品，用于基因組DNA的提取與分析。

1.2 DNA提取與高通量測序

分別從兩種耳型(正常耳組和小耳組)藏山羊中，各選擇20個個體的耳朵樣本，利用DNA提取試劑盒(TIANamp Genomic DNA Kit)提取山羊基因組DNA；質(zhì)量檢測合格后，分別將正常耳和小耳組不同個體的DNA樣品等量混合，形成2個DNA池。將混合后的DNA樣品送至北京諾禾致源科技股份有限公司；構(gòu)建高通量測序文庫，并采用Illumina HiSeq X10測序儀進行雙末端 (2×150)測序。

1.3 原始讀段質(zhì)量控制與比對

首先，過濾掉原始數(shù)據(jù)中包含接頭的讀段，然后使用Trimmomatic[18](v0.36)過濾低質(zhì)量讀段(參數(shù)：LEADING:20,TRAILING:20,SLIDINGWINDOW:4:20,MINLEN:50)。使用BWA[19](v0.7.12)將高質(zhì)量讀段比對到山羊參考基因組[20](assembly ARS1，https://asia.ensembl.org/index.htm)。使用Picard(v2.10.6，http://broadinstitute.github.io/picard/)去除重復讀段，然后使用GATK[21](v3.8-0)進行局部重比對和堿基質(zhì)量值校正。

1.4 SNPs的檢測與質(zhì)量控制

為了更準確地分析藏山羊群體特征，本研究還使用前期獲得的5個山羊群體的測序數(shù)據(jù)[22]。應(yīng)用Varscan[23](v2.4.3)檢測波爾山羊(BE)、美姑山羊(MG)、金堂山羊(JT)、南江黃羊(NJ)、藏山羊(措勤，ZZ)以及正常耳組的藏山羊(NE)、小耳組藏山羊(SE)共7個混池樣品的SNPs(參數(shù)：mpileup2snp-min-coverage 15-min-reads2 2-min-avg-qual 30-min-var-freq 0.01)；然后對原始SNPs進行過濾，只保留最小等位基因頻率(MAF)> 0.05和缺失基因型個體數(shù)量≤ 1的SNPs。另外，從7個群體的SNPs數(shù)據(jù)集里提取并合并正常耳和小耳組的SNPs數(shù)據(jù)，重新計算等位基因頻率，過濾掉MAF<0.05的位點，得到仲巴縣藏山羊群體的SNPs，最后使用SnpEff[24]對SNPs進行注釋。

1.5 全基因組選擇信號的檢測與基因富集分析

參考已有研究[22]，采用滑動窗口的方法(窗口長度為100 kb，步長為25 kb)，在全基因組范圍內(nèi)，計算正常耳組和小耳組藏山羊之間的群體分化指數(shù)(FST)以及期望雜合度比值(HPratio)。將分布在前5%窗口的FST值和HP比值對應(yīng)窗口的交集區(qū)域定義為受選擇區(qū)域。根據(jù)注釋結(jié)果，與受選擇區(qū)域重合的基因被確定為受選擇 基因。

利用Enrichr[25]對篩選出的受選擇基因進行MGI哺乳動物表型富集分析，篩選與哺乳動物耳朵發(fā)育相關(guān)的條目(P<0.05)。最后，計算相關(guān)選擇區(qū)域內(nèi)SNPs位點的等位基因頻率差異(△AF=|AFSE-AFNE|，AFSE和AFNE分別是小耳組和正常耳組的突變型等位基因頻率)；參考已有研究[26]，選擇△AF>0.8的SNPs位點進一步鑒定影響藏山羊耳型的遺傳變異。

1.6 主成分分析

使用GCTA[27](v1.26.0)進行主成分分析，探討本研究中7個山羊品種的遺傳組成特點。

2 結(jié)果與分析

2.1 耳長和耳寬表型值特征

由圖1-B所示，正常耳組藏山羊的耳朵長度、寬度均值分別為12.3 cm和6.4 cm；小耳組藏山羊耳朵長度、寬度均值分別為5.1 cm和4.0 cm。t檢驗結(jié)果表明，正常耳組的耳長、耳寬均極顯著高于小耳組(P<0.01)。

A. 藏山羊小耳與正常耳 Tibetan goats with normal and small ear size；B. 兩種耳型的耳長與耳寬的比較 Comparison of the ear length and width between the two types of ear；紅點和藍點分別代表耳長和耳寬 Red and blue points represent the ear length and width,respectively;SE和NE分別為小耳組和正常耳組 SE and NE represent the small ear group and normal ear group,respectively.“**”表示t檢驗差異極顯著(P<0.01) The “**” indicate significant difference underttest(P<0.01)

圖1 兩種耳型藏山羊耳長與耳寬的比較
Fig.1 Comparison of ear length and width between Tibetan goats andother two different ear types

2.2 讀段比對和SNPs檢測

基于正常耳組藏山羊構(gòu)建的測序文庫共包含421 120 994條高質(zhì)量讀段，測序深度約為22.1×；基于小耳組藏山羊構(gòu)建的測序文庫共包含345 842 498條高質(zhì)量讀段，測序深度約為21.2×。序列比對結(jié)果顯示，正常耳組和小耳組分別有414 280 705和340 355 762條高質(zhì)量讀段比對到山羊參考基因組，比對率分別為98.38%和 98.41%。

但是不同于祖母所代表的墨西哥傳統(tǒng)女性，賽利亞意識到依靠男性并不能真正使自己擺脫困境，并鼓勵女性實現(xiàn)自身價值，這是她不同于墨西哥傳統(tǒng)女性的女性覺醒。

基于混池測序數(shù)據(jù)，首先在7個山羊群體中共檢測到6 959 360個高質(zhì)量SNPs。進一步在仲巴藏山羊群體中檢測到3 662 933個SNPs(表1)，SNPs的平均密度為2.82個/kb。

表1 仲巴縣藏山羊29條常染色體的SNPs數(shù)量分布及密度Table 1 Distribution and density of SNPs across 29 autosomes in Tibetan goats from Zhongbacounty

2.3 選擇信號的鑒定與受選擇基因的富集分析

如圖 2 所示，依據(jù)5%的閾值，共有4 899 個窗口的FST值大于0.12，其中2號染色體的 56.475～56.575 Mb區(qū)域在兩種耳型藏山羊群體間分化程度最高(FST=0.36)?；蜃⑨尳Y(jié)果顯示，與該窗口重合的基因包括STAT1和GLS。HP比值的計算結(jié)果顯示，共有4 918 個窗口的HP比值大于1.21(前5%的閾值)，最大值為 5.26。結(jié)合FST與HP比值，最終共有894 個窗口被確定為受選擇區(qū)域。根據(jù)基因注釋結(jié)果，在受選擇區(qū)域內(nèi)總共鑒定到 372 個受選擇基因；其中，25% 的基因(93 個)位于7號染色體，并且在這93 個基因中，56%的基因(52 個)位于45.05～59.76 Mb區(qū)域內(nèi)。

圖2 兩種耳型藏山羊全基因組的選擇信號分布Fig.2 Genome-wide distributions of selection signals betweenTibetan goats andother two different ear types

如圖3所示，通過對 372 個受選擇基因進行富集分析，共鑒定到10 個可能與哺乳動物耳朵發(fā)育相關(guān)的表型條目(P<0.05)，具體包括：破骨細胞數(shù)量減少(MC4R、SPARC、CAMK4、RGS10、CCR2)；普通髓系祖細胞形態(tài)異常 (TSTA3、KMT2A、STAT1、TNFRSF9、TYK2)；背根神經(jīng)節(jié)形態(tài)異常(UCHL1、KMT2A、TCOF1、NDN、GDAP1)；皮膚易脫落(ERRFI1、SPINK5、CYP7A1)；小型骨形態(tài)異常(FAM210A、ENAM、SLC4A4、PPARD)；中線腭架融合(TCOF1、PDGFC、GLI2)；骨強度降低(FAM210A、BBX、PPARD)；成骨細胞分化受損(SATB2、PPARD)；顳骨形態(tài)異常(TCOF1、SATB2)；額骨形態(tài)異常(TCOF1、SATB2、GLI2)。

圖3 與哺乳動物耳朵發(fā)育相關(guān)的表型條目Fig.3 MGI-MP terms related to ear development in mammals

2.4 選擇信號區(qū)域內(nèi)SNPs位點的等位基因頻率差異

如圖4所示，在全基因組FST值最大的滑動窗口(chr2：56.475～56.575 Mb)內(nèi)共有77個SNPs位點，其中30 個的等位基因頻率在兩種耳型藏山羊之間差異較大(△AF>0.5)，其中1個位于GLS基因內(nèi)含子的SNP位點(chr2:56559507,C>T)的△AF高達0.87。

對于7 號染色體，共有1 884 個SNPs 位點的△AF>0.5，10 個SNPs位點的△AF>0.8，其中7 個均在45.05～59.76 Mb區(qū)域內(nèi)；該區(qū)域內(nèi)△AF的峰值位于SPINK5基因上游約5 kb處，對應(yīng)SNPs位點(c.-4283A>G、c.-4295G>A)的△AF分別為0.9和0.83(圖5)。

圖4 2號染色體56.475～56.575 Mb區(qū)域內(nèi)SNPs位點的等位基因頻率差異Fig.4 Allele frequency differences of SNPs in the region of 56.475-56.575 Mb on chromosome 2

圖5 7號染色體SNPs位點的等位基因頻率差異Fig.5 Allele frequency differences of SNPs onchromosome 7

2.5 主成分分析

主成分分析的結(jié)果顯示(圖6)，第一和第二特征向量分別解釋25.95%和19.04%的遺傳變異，7個山羊群體聚集成4個簇。其中，來自仲巴縣的正常耳(NE)組和小耳(SE)組藏山羊聚在一起；來源于四川盆地的金堂黑山羊(JT)、南江黃羊(NJ)、美姑山羊(MG)聚類為一簇；而波爾山羊(BE)和來自西藏措勤縣的藏山羊(ZZ)均未與其他群體聚類到一起。

圖6 7個山羊群體的主成分分析Fig.6 Principal component analysis of seven goat populations

3 討 論

耳型是現(xiàn)代家畜品種的一個重要表型性狀；家畜耳型的特征主要表現(xiàn)在耳朵大小、形狀和直立性三方面。盡管有關(guān)山羊耳型變異的研究已有報道[3,11]；但總體來說，人們對山羊耳型變異的分子機制仍知之甚少。本研究基于全基因組SNPs信息，利用比較基因組學的分析方法，首先發(fā)現(xiàn)，2號染色體的56.475～56.575 Mb基因組區(qū)域在兩種耳型藏山羊之間遺傳分化程度最高，并且多個位點表現(xiàn)出較高的等位基因頻率差異。雖然GLS基因與該區(qū)域重合部分最多，但目前未見GLS基因與耳朵發(fā)育有直接或間接關(guān)聯(lián)的報道。

值得注意的是，在染色體水平上，本研究鑒定到的受選擇基因在7號染色體上分布最多，且主要集中在45.05～59.76 Mb區(qū)域。而Brito等[3]利用努比亞山羊、波爾山羊、薩能山羊等長耳、中等耳長的山羊品種作為對照，鑒定到7 號染色體的48.5～57.3 Mb區(qū)域是影響拉曼查山羊小耳的一個基因組區(qū)域，這與本研究結(jié)果較一致。已有研究表明，該區(qū)域內(nèi)存在多個與哺乳動物耳朵發(fā)育相關(guān)的基因；例如，SPINK5基因通過調(diào)控皮膚表皮蛋白分解速度參與皮膚形成過程[28]。而ANXA6基因通過調(diào)控ERK-MAPK信號通路的活性來調(diào)節(jié)終末軟骨細胞分化[29]。Xiang等[30]利用小鼠模型發(fā)現(xiàn)，敲除POU4F3基因?qū)е聝?nèi)耳毛細胞的退化。已有研究表明，F(xiàn)GF10基因調(diào)控小鼠內(nèi)耳發(fā)育[7]，而NDST1基因在FGF信號轉(zhuǎn)導過程中發(fā)揮作用[31]。因此，推測NDST1基因間接參與小鼠內(nèi)耳的形成。簡言之，基于拉曼查山羊與本研究中小耳型藏山羊相似的小耳表型，該基因組區(qū)域可以被認為是影響山羊耳型變異的遺傳位點，但是具體的因果突變?nèi)孕柽M一步研究。另外，在豬和綿羊上報道的與耳型變異相關(guān)的基因HOXA1[9]、MSRB3[10]在本研究中并未被鑒定到，表明在哺乳動物中存在多種耳型變異機制。

雖然生活在青藏高原及其毗鄰地區(qū)的固有山羊品種均被稱為藏山羊[32]，但青藏高原幅員遼闊，不同地區(qū)的藏山羊群體在體型、毛色、角型等性狀上表現(xiàn)出顯著的區(qū)別?；谖⑿l(wèi)星標記的研究表明，藏山羊具有豐富的遺傳多樣性，但來自不同生態(tài)區(qū)的藏山羊群體間存在明顯的遺傳分化[13,33]。利用外顯子組測序技術(shù)，Song等[34]比較了2個西藏絨山羊品種和7個低海拔絨山羊品種的群體分化指數(shù)，結(jié)果表明，西藏絨山羊與低海拔絨山羊之間的遺傳分化程度更高。最近，Deng等[35]通過對線粒體D-loop區(qū)的分析，發(fā)現(xiàn)拉薩、曲水、薩嘎等10個地區(qū)的藏山羊群體間遺傳分化明顯。在本研究中，仲巴藏山羊與金堂黑山羊、南江黃羊等低海拔品種的遺傳距離更近，與措勤地區(qū)的藏山羊較遠，表明上述兩個地區(qū)藏山羊群體的遺傳組成差異較大。