思南黃牛體重和體尺性狀的GWAS分析

袁 揚，趙 恬，姚 丹，張定紅，王胤晨*，吳鴻友，汪 勇

（1.貴州省農(nóng)科院畜牧獸醫(yī)研究所，貴州貴陽 550005；2.貴州省貴陽市南明區(qū)農(nóng)業(yè)農(nóng)村局，貴州貴陽 550005；3.貴州省羅甸縣農(nóng)業(yè)農(nóng)村局，貴州羅甸 550100）

近年來全基因組關(guān)聯(lián)分析（Genome-Wide Associate Analysis，GWAS）較多應(yīng)用于家畜的遺傳分析和輔助育種。大多數(shù)家畜的GWAS研究都與其經(jīng)濟性狀相關(guān)，如豬的雄烯酮水平、毛被顏色、脂肪性狀等相關(guān)GWAS研究，馬的肌肉表型的研究，綿羊羊角形狀研究，雞的體重影響因素研究和產(chǎn)蛋研究。大多數(shù)農(nóng)作物和家畜主要性狀是多基因共同作用的，遺傳結(jié)構(gòu)較復(fù)雜，難以預(yù)測出對性狀有重要影響的單個基因。GWAS技術(shù)能夠在全基因組水平上搜索突變位點，綜合性狀給出性狀相關(guān)的位點突變和性狀候選基因，極大地促進對家畜復(fù)雜性狀遺傳結(jié)構(gòu)的理解，改善家畜育種計劃。

關(guān)于牛的體重性狀、體尺性狀GWAS研究非常多。Sorbolini等人使用混合線型模型對馬其頓牛生長性狀和胴體性狀進行分析，根據(jù)SNP位置和連鎖不平衡確定了影響性狀的目的區(qū)域，搜索發(fā)現(xiàn)和3個基因上出現(xiàn)了影響脂質(zhì)代謝、骨骼和肌肉的有意義位點。Ede等使用混合線性模型分析SNP芯片發(fā)現(xiàn)與韓國牛體重、胴體性狀相關(guān)的多個SNP位點。Adoligbe等研究了中國本土牛中基因的SNP及其與體尺性狀的關(guān)聯(lián)，研究顯示GDF-10存在SNP位點11471（A/G），該突變在秦川牛、郟縣紅牛和南陽牛中顯著影響體尺性狀，SNP位點G142A對郟縣紅牛臀部寬度、胸寬和胸圍都有顯著影響，認為基因可能是郟縣紅牛繁育篩選過程中體尺性狀的重要遺傳標記。李天科等分析甘南牦?；蚨鄳B(tài)性和生產(chǎn)性狀相關(guān)性，發(fā)現(xiàn)12116（G/A）、12152（C/T）和13041（T/C）3個SNP位點，結(jié)果表明3個位點之間存在弱連鎖平衡，不同突變位點基因型與體斜長、體高、胸圍、體重差異顯著或極顯著，與管圍差異不顯著。郭憲等研究南德溫雜交肉牛發(fā)現(xiàn)1113（G/A）、9569（G/A）和9628（G/A）3個SNP位點與南德溫雜交肉牛的體高、體斜長、胸圍和管圍顯著相關(guān)，對肉牛體尺性狀有影響。

思南黃牛是中國貴州本土優(yōu)良黃牛品種，體質(zhì)緊湊、四肢強健、善于爬山，適于山區(qū)耕作放牧，肉用性能良好。本研究以思南黃牛、西門塔爾牛以及二者雜交改良后代為研究對象，從基因組水平對牛體重和體尺性狀進行關(guān)聯(lián)分析，為肉牛新品種的選育提供技術(shù)參考與借鑒。

1 材料與方法

1.1 實驗材料 實驗動物群體選自貴州省石阡縣聚鳳鄉(xiāng)彭寨村肉牛養(yǎng)殖合作社，在相同條件下飼養(yǎng)1年，對其中112頭思南黃牛（12月齡～18月齡）、17頭西門塔爾牛（12月齡～18月齡）以及19頭本地雜交牛（18月齡～24月齡）共計148頭牛，測量樣本的體重、體高、體長、胸寬、胸圍、管圍等體重和體尺性狀數(shù)據(jù)（測量方法參考DB52/T 1257.2-2017），采集148頭牛耳靜脈血樣本并提取DNA。

1.2 實驗儀器與試劑 磁珠法血液基因組DNA提取試劑盒（北京天根生化科技有限公司 DP329），通用臺式冷凍離心機（epServices 5805）。樣品送至北京諾禾致源生物信息科技有限公司完成全基因組關(guān)聯(lián)分析。

1.3 血液基因組提取及GBS文庫構(gòu)建及測序

1.3.1 血液基因組DNA提取 使用磁珠法血液基因組DNA提取試劑盒（DP348-03）對血液樣品的基因組DNA進行提取，操作步驟按照試劑盒說明書進行。

1.3.2 血液基因組DNA質(zhì)量檢測 DNA 3.5 μL（2 μL DNA 樣品，1.5 μL 6×Loading Buffer，混勻）上樣于1%瓊脂糖凝膠，120 V電泳18 min，凝膠成像系統(tǒng)觀察結(jié)果。

1.3.3 GBS文庫構(gòu)建及測序 檢測合格的DNA樣品經(jīng)限制性內(nèi)切酶對基因組進行酶切，加上帶有barcode的接頭后，對每個樣品進行擴增，然后對樣品進行混合，選擇需要的片段進行文庫構(gòu)建，利用Illumina HiSeq測序平臺進行雙端150測序（北京諾禾致源生物信息科技有限公司）。步驟如下：

1）酶切：0.1～1 μg基因組DNA用限制性內(nèi)切酶進行酶切，以得到適合的marker密度；

西北農(nóng)林科技大學從美國引育的早熟新品種，為皇家嘎拉芽變。2006年通過了陜西省果樹品種審定委員會審定。果實圓錐形，平均單果重180克；蓋色濃紅，著色有條紋；果皮光滑，果點大；果肉黃白色，風味酸甜適度；可溶性固形物含量14.6%，可滴定酸含量0.22%，果肉硬度11.3公斤/平方厘米；肉質(zhì)硬脆，汁液多，有香氣，為鮮食品種。果實發(fā)育期120天左右。豐產(chǎn)性好。紅蓋露樹體適應(yīng)性強，耐瘠薄，對早期落葉病、白粉病有一定的抗性。

2）加P1和P2接頭：酶切后的片段兩端加P1和P2 Adapter（可與酶切DNA缺口互補）；

3）片段選擇：PCR擴增兩端分別含有P1和P2接頭的tag序列，DNA片段pooling，電泳回收需要區(qū)間的DNA；

4）高通量測序：Cluster制備，上機測序。

5）文庫構(gòu)建完成后，使用Qubit 2.0進行初步定量，稀釋文庫至1 ng/μL，隨后使用Agilent 2100對文庫的insert size進行檢測，insert size符合預(yù)期后，使用Q-PCR方法對文庫的有效濃度進行準確定量（文庫有效濃度＞2 nM），以保證文庫質(zhì)量。

6）庫檢合格后，把不同文庫按照有效濃度及目標下機數(shù)據(jù)量的需求pooling后進行Illumina Hiseq PE150測序。

1.4 基因組測序、比對和SNP檢測 對Illumina HiSeq測序平臺雙端測序下機原始數(shù)據(jù)預(yù)處理，過濾掉含有接頭序列的測序序列、低質(zhì)量測序序列、N（未測出的堿基）含量超過該條測序序列長度10%的成對測序序列，最后得到高質(zhì)量測序序列。以歐洲普通牛（）為參考基因組（ftp://ftp.ensembl.org/pub/release-86/fasta/bos_taurus/），基因組大小為2 649 682 029堿基對（不含N），將所有樣本作為一個群體與參考基因組比對，使用bwa軟件（參數(shù)：mem-t4 -k32 -M）將雙端高質(zhì)量測序序列比對到參考基因組，使用samtools軟件得到變異的vcf文件，共獲得了3 601 237 個SNP位點，經(jīng)過測序深度3x、Miss 0.2和次要等位基因頻率（MAF）＞0.01的條件過濾后，最后共獲得441 548個SNP位點用于后續(xù)分析。

1.6 全基因組關(guān)聯(lián)分析 以歐洲普通牛全基因組序列為參考基因組（ftp://ftp.ensembl.org/pub/release-86/fasta/bos_taurus/，Ensembl UMD3.1），使 用ANNOVAR軟件對牛體尺性狀相關(guān)的SNPs進行功能注釋，過濾后的高質(zhì)量SNP使用GEMMA軟件混合線性模型進行關(guān)聯(lián)分析。分析模型為：y=X+Z+W+e，其中，y為所要研究的表型性狀；X為固定效應(yīng)（Fixed Effect），影響y的其他因素，包括群體結(jié)構(gòu)、性別、年齡等因素；Z，標記效應(yīng)（Marker Effect）；W，隨機效應(yīng)（Random Effect），一般指個體的親緣關(guān)系；e，隨機誤差。采用Bonferoni correction方法來確定閾值。最后將關(guān)聯(lián)到的結(jié)果用R軟件繪制曼哈頓圖和Q-Q圖進行展示。

2 結(jié)果

2.1 樣本表型測量結(jié)果 本實驗測量了牛的體重（Body Weight，BW）、體高（Body High，BH）、體長（Body Length，BL）、胸寬（Chest Breadth，CB）、胸圍（Chest Circumference，CC）、管圍（Circumference of Cannon Bone，CCB）性狀，這6個體尺性狀的描述性統(tǒng)計見表1。

表1 樣本體尺性狀的描述性統(tǒng)計

2.2 GBS文庫測序和SNP檢測結(jié)果 本次選取了148頭牛的樣本數(shù)據(jù)，測序共產(chǎn)生98.396 G Raw data，平均每個樣品664.835M Raw data，過濾后的Clean data共98.394G，平均每個樣品664.823M。以歐洲普通?；蚪M為參考序列，將每個樣本測序得到的clean data和基因組進行比對，統(tǒng)計得到的平均比對率、平均測序深度和平均覆蓋度見表2。

表2 測序質(zhì)量

質(zhì)控后數(shù)據(jù)與參考基因組比較，在思南黃牛群體中鑒定出SNP位點3 601 237個，平均SNP密度為100.627個/Mb。對獲得的SNP質(zhì)控后，獲得高質(zhì)量SNP位點441 548個，平均SNP密度為12.338個/Mb。

2.3 群體結(jié)構(gòu)分析結(jié)果 主成分分析將牛群分為3個類群，樣本牛群具有明顯分層（圖1A）。通過顯著性檢驗糾正不同遺傳因素和繁殖條件導致的種群分層。SNP位點構(gòu)建的進化樹在整體上聚為一個大分支和兩個小分支，結(jié)果和主成分分析以及種群結(jié)構(gòu)分析相吻合（圖1B）。Structure軟件分析種群結(jié)構(gòu)，K值為3時，148個樣本聚類為3簇，此時最低交叉驗證值CV error（K=3）=0.464 84，與主成分分析結(jié)果相近（圖1C）。

圖1 基于SNPs的群體結(jié)構(gòu)分析

2.4 連鎖不平衡分析結(jié)果 連鎖不平衡分析是全基因組關(guān)聯(lián)分析的基礎(chǔ)。本研究對148個牛外周血樣本進行連鎖不平衡分析，計算等位基因相關(guān)系數(shù)（），分析148個樣本441 548個SNP位點在?；蚪M的LD共線性。根據(jù)位點之間的物理距離，得到連鎖不平衡度（LD度）（圖2）。結(jié)果顯示，LD度隨距離增加迅速下降。經(jīng)過注釋發(fā)現(xiàn)，SNP在基因之間分布最多，其次是內(nèi)含子，最后是外顯子。

圖2 連鎖不平衡圖

2.5 全基因組關(guān)聯(lián)分析 主成分分析結(jié)果作為協(xié)變量，利用GEMMA軟件混合線性模型對牛的6個體尺性狀表現(xiàn)型進行全基因組關(guān)聯(lián)分析（GWAS），得到牛體尺性狀的GWAS正態(tài)化結(jié)果曼哈頓圖和Q-Q圖（圖3）。圖3左側(cè)曼哈頓圖顯示了GWAS的優(yōu)勢位點，信號與原假設(shè)之間的偏差在右側(cè)以Q-Q圖的形式展示。通過Q-Q圖觀察發(fā)現(xiàn)，觀測值與期望值基本在一條斜線上，每個性狀Q-Q圖結(jié)果僅在尾部區(qū)域稍微偏離值分布，表明混合線性模型不存在群體分層現(xiàn)象（λ接近1），SNP關(guān)聯(lián)分析的假設(shè)與統(tǒng)計分布沒有偏離，模型適合該資源群體。曼哈頓圖中水平線是Bonferroni顯著性閾值，本次關(guān)聯(lián)分析共發(fā)現(xiàn)18個基因中的51個SNP與體重、體尺性狀變化相關(guān)，篩選出體重、體高、體長、胸寬、胸圍、管圍等身體特征數(shù)據(jù)相關(guān)的基因分別有7、4、2、1、0個和4個（表3）。

表3 牛體重和體尺性狀相關(guān)候選基因列表

圖3 混合線性模型對體重（A）、體高（B）、體長（C）、胸寬（D）、胸圍（E）、管圍（F）進行GWAS分析

3 討 論

肉牛是全球主要肉類來源之一，體尺性狀是衡量肉牛品種質(zhì)量的指標。體尺性狀和牛的健康、生產(chǎn)率、壽命息息相關(guān)。通過GWAS分析篩選牛體尺性狀相關(guān)的候選基因，分別得到23、5、3、6、1、13個和體重、體高、體長、胸寬、胸圍、管圍等體重、體尺性狀相關(guān)的SNP，這些SNP分別在體重的7個、體高的4個、體長的2個、胸寬的1個和管圍的4個基因內(nèi)或基因附近。

3.1 體重性狀關(guān)聯(lián)基因 牛的體重性狀相關(guān)基因有（）、、（）、（）、（）、（）和（）。是牛肉品質(zhì)的差異表達基因，影響牛的營養(yǎng)攝入，與牛體重性狀相關(guān)。小鼠中IFT57蛋白的無效突變表現(xiàn)出纖毛病的特征，如上肢和下肢多指、上頜突縮短、神經(jīng)管缺陷等。人的純合突變會導致未分類的口面指綜合征（Oral-facial-digital Syndrome，OFDS），表現(xiàn)為口腔、面部和四肢畸形，同時也會導致身材矮小。腺苷脫氨酶（）基因可能通過脂類代謝影響牛皮下脂肪積累?；蚺c牛的性格發(fā)展有關(guān)，這影響了牛的活動頻率、情緒控制和順從性等性格，可能導致牛體重的改變。降鈣素基因相關(guān)肽（）能夠影響哺乳動物骨密度，基因敲除小鼠的骨密度會增加，還導致小鼠去脂體重增加，說明信號參與小鼠骨骼分解代謝，同時可能與腸代謝有關(guān)。主要功能是促進消炎和免疫抑制性環(huán)境，多態(tài)性與人類炎癥性腸病有關(guān)。此外，關(guān)于眼皮膚白化病2型基因（）之前并無關(guān)于體重性狀的報告，是一個新的體重性狀關(guān)聯(lián)候選基因。

3.2 管圍性狀關(guān)聯(lián)基因 牛的管圍性狀關(guān)聯(lián)基因有（）、（）、（）和（）。也被稱為，與原因不明的腸道疾病克羅恩病有關(guān)，該基因可能與消化功能相關(guān)。RAB7B是一種GTPase，位于溶酶體、高爾基體和反式高爾基體網(wǎng)絡(luò)，調(diào)節(jié)內(nèi)體到高爾基體的運輸，基因能夠影響牛采食量與體重，與日本黑牛飼料利用和生長特性相關(guān)。在牛飲食受限時，該基因表達量也發(fā)生變化，對牛飲食限制及恢復(fù)，牛瘤胃上皮轉(zhuǎn)錄組分析顯示基因出現(xiàn)差異表達。基因還與牛背最長肌的生長相關(guān)，也是小牛成骨不全癥可能的候選基因?；虻腟NP與豬的管圍性狀相關(guān)聯(lián)，處于豬剩余飼料攝入量相關(guān)的重要候選QTL區(qū)域。

3.3 體長性狀關(guān)聯(lián)基因 牛的體長相關(guān)的基因有（）和（）。位于牛的腿型、蹄型關(guān)聯(lián)區(qū)域，有報道附近存在脂肪酸性狀相關(guān)的SNP位點，與牛的骨骼生長和脂類代謝相關(guān)，一定程度上能夠解釋其與體長性狀的關(guān)聯(lián)?；虼嬖谌私?、遠端結(jié)腸的差異甲基化位點，與消化系統(tǒng)的甲基化相關(guān)，由此可能會影響牛的體長性狀。

3.4 身高性狀關(guān)聯(lián)基因 牛的身高性狀關(guān)聯(lián)基因有（）、（）、（）和（）。和都位于8號染色體上，都是牛的拷貝數(shù)變異基因。基因和基因都位于內(nèi)羅爾牛背最長肌組織中鋅含量有關(guān)QTL區(qū)域內(nèi)。COL12A1是一種膠原蛋白，表達量隨著牛背最長肌的發(fā)育而逐漸下降，該基因位于青海高原牛和金川牛胸椎骨長度選擇壓力區(qū)域，牛胸椎骨的節(jié)數(shù)和長度直接影響牛的骨架，從而和身高產(chǎn)生關(guān)聯(lián)。的突變?nèi)笔е录∪夤δ苄院土α總鲗芰ο陆?，此外，基因還是牛背最長肌組織中鐵和銅含量重要的差異表達基因。SNBT1是一種外周膜蛋白，參與調(diào)控肌肉收縮過程，是快紅肌和慢紅肌的表達差異最大基因之一，同時處于牛背最長肌組織中的單不飽和脂肪酸相關(guān)區(qū)域，可能也在牛的脂質(zhì)代謝中起到一定作用。

3.5 胸圍性狀關(guān)聯(lián)基因 牛的胸圍性狀關(guān)聯(lián)基因有（）。位于影響魚骨長的QTL區(qū)域，能夠影響魚頭骨大小。還是豬體內(nèi)共表達網(wǎng)絡(luò)中肌肉肉質(zhì)性狀關(guān)聯(lián)基因。

綜上所述，在以上篩選的基因中：能夠影響動物體型；和能夠影響脂質(zhì)代謝；影響牛的性格；影響動物骨骼生長發(fā)育；和影響動物消化系統(tǒng)；和影響動物的采食量；影響牛背最長肌組織中礦物質(zhì)含量。這些基因共同作用，綜合影響了牛的體重和體尺性狀。

4 結(jié) 論

本研究通過對112頭思南黃牛、17頭西門塔爾牛以及19頭本地雜交牛共計148個個體生產(chǎn)性狀進行GWAS分析，篩選出體重、體尺性狀相關(guān)的24個SNPs位于18個基因附近，發(fā)現(xiàn)和分別是體重、體高、體長、胸寬和管圍相關(guān)性狀最有相關(guān)性的候選基因，這些候選基因?qū)⒂兄谶M一步探索思南黃牛的遺傳改良性狀，為后續(xù)思南黃牛遺傳信息精細作圖提供數(shù)據(jù)，為思南黃牛育種研究提供參考。