戴瑞奶綿羊產(chǎn)羔數(shù)全基因組關(guān)聯(lián)分析

米布農(nóng)，張立果，烏日漢，郭玉林，王春偉，徐全忠，馮爽，李光鵬，蘇小虎*，張立,*

（1.內(nèi)蒙古大學(xué)生命科學(xué)學(xué)院，內(nèi)蒙古呼和浩特 010070；2.烏蘭察布市家畜改良工作站，內(nèi)蒙古烏蘭察布 012200）

戴瑞羊（Dairy Meade，DM）是新西蘭以本地柯泊華茲羊?yàn)槟副?、東弗里生羊?yàn)楦副?，?jīng)過20 多年培育形成的奶綿羊新品種，具有泌乳性能優(yōu)良、繁殖力較強(qiáng)（220%）等特點(diǎn)。我國(guó)歷史上沒有專用的奶綿羊品種，為填補(bǔ)這一空白，內(nèi)蒙古大學(xué)于2016 年引入了戴瑞奶綿羊胚胎，通過胚胎移植獲得了國(guó)內(nèi)純種戴瑞奶綿羊，隨后又以本地小尾寒羊?yàn)槟副荆魅鹧驗(yàn)楦副?，進(jìn)行新品種培育。

奶綿羊新品種培育的核心是建立高效的選育技術(shù)，由于綿羊產(chǎn)羔性狀受到多基因調(diào)控，遺傳力低且遺傳進(jìn)展緩慢，傳統(tǒng)育種方式不能在短時(shí)期內(nèi)取得理想的遺傳進(jìn)展，這在很大程度上阻礙了綿羊產(chǎn)業(yè)的快速發(fā)展。隨著測(cè)序技術(shù)和生物信息學(xué)的高速發(fā)展，全基因組關(guān)聯(lián)分析（Genome Wide Association Study，GWAS）已成為準(zhǔn)確鑒定產(chǎn)羔性狀遺傳標(biāo)記和候選基因的重要手段。目前在綿羊產(chǎn)羔性狀方面有很多研究報(bào)道。Smolucha等對(duì)波蘭山綿羊群體進(jìn)行GWAS 研究，發(fā)現(xiàn)了1 個(gè)重要的SNP位于中。Hernández-Montiel 等通過Illumina Ovine SNP50 分型進(jìn)行GWAS 分析，發(fā)現(xiàn)了與Pelibuey 羊產(chǎn)羔數(shù)相關(guān)的染色體水平上顯著關(guān)聯(lián)的3 個(gè)SNPs、54 個(gè)潛在SNPs 以及10 個(gè)潛在候選基因。Abdoli 等對(duì)伊朗Lori-Bakhtiari 綿羊品種進(jìn)行全基因組關(guān)聯(lián)分析，在3 號(hào)染色體上鑒定出2 個(gè)SNPs 和1 個(gè)新的候選基因。本研究基于前期已形成的戴瑞奶綿羊群體，通過對(duì)其產(chǎn)羔性狀進(jìn)行測(cè)定和基因組重測(cè)序，利用GWAS 方法挖掘與產(chǎn)羔性狀相關(guān)的候選基因及關(guān)鍵位點(diǎn)，為奶綿羊新品種的培育提供理論技術(shù)支持。

1 材料與方法

1.1 實(shí)驗(yàn)材料 本實(shí)驗(yàn)所用的135 只戴瑞奶綿羊繁育群體（15～18 月齡）均來源于內(nèi)蒙古自治區(qū)烏蘭察布市蒙天然牧業(yè)科技發(fā)展有限公司。表型均為頭胎產(chǎn)羔數(shù)。取頸靜脈EDTA 抗凝血2 mL，-20℃保存?zhèn)溆谩?/p>

1.2 基因組DNA 提取及檢測(cè) 將采集的血樣送至北京康普森公司，提取基因組DNA 及全基因組重測(cè)序。

1.3 基因型數(shù)據(jù)獲取及質(zhì)量控制 采用Illumina Novaseq 6000 平臺(tái)對(duì)奶綿羊基因組進(jìn)行重測(cè)序。有效測(cè)序數(shù)據(jù)通過BWA 0.7.17軟件比對(duì)到參考基因組，比對(duì)結(jié)果經(jīng)SAMTOOLS 1.7軟件去除重復(fù)；并將其映射到綿羊Oar_rambouillet_v1.0 參考序列；映射的讀段通過SAMTOOLS 分析獲取變異位點(diǎn)的分型；使用Plink 1.9進(jìn)行質(zhì)控，質(zhì)控標(biāo)準(zhǔn)為：剔除SNP 缺失大于10%的位點(diǎn)，剔除最小等位基因頻率小于1%的位點(diǎn)和哈迪-溫伯格平衡檢驗(yàn)值小于10的位點(diǎn)。

1.4 樣本群體結(jié)構(gòu)分析 為避免群體分層導(dǎo)致的假陽性結(jié)果，采用Plink 1.9 對(duì)SNP 分型數(shù)據(jù)計(jì)算IBS（Identity-By-State）遺傳距離矩陣和主成分分析（Principal Components Analysis）。通過Plink 1.9 軟件構(gòu)建親緣關(guān)系矩陣（G 陣），應(yīng)用R 語言的gg plots包的heatmap 函數(shù)繪制熱圖，采用解釋方差最大的前2個(gè)主成分進(jìn)行散點(diǎn)圖繪制以分析群體結(jié)構(gòu)。

1.5 全基因組關(guān)聯(lián)分析及統(tǒng)計(jì)推斷 應(yīng)用GEMMA 0.98.1軟件進(jìn)行混合線性模型分析，模型如下：

其中，y 為表型向量，X為群體結(jié)構(gòu)效應(yīng)和場(chǎng)、年、季等固定效應(yīng)，Ζγ為待檢驗(yàn)標(biāo)記效應(yīng)，～（0，）為多基因效應(yīng)，～（0，）為殘差效應(yīng)。多基因效應(yīng)中的K 為標(biāo)記推斷的親緣關(guān)系矩陣。最終，每個(gè)SNP 位點(diǎn)對(duì)應(yīng)一個(gè)關(guān)聯(lián)值。根據(jù)Bonferroni 校正法，確定產(chǎn)羔性狀全基因組關(guān)聯(lián)分析顯著性閾值（0.05/21 020 909）、潛在顯著閾值（1/21 020 909）。曼哈頓圖和分位數(shù)圖由GEMMA 繪制。

1.6 基因注釋及候選基因篩選 利用ENSEMBL 網(wǎng)站下載綿羊基因組序列Oar_rambouillet_v1.0（ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/002/742/125/GCF_00274 2125.1_Oar_rambouillet_v1.0/GCF_002742125.1Oar_ra mbouillet_v1.0_genomic.fna.g），使用ANNOVAR 軟件將顯著SNP 注釋到其對(duì)應(yīng)的基因上。使用cluster Profiler 包將注釋到的候選基因根據(jù)GO 和KEGG 數(shù)據(jù)庫(kù)進(jìn)行基因功能富集分析。通過文獻(xiàn)查閱基因生物學(xué)功能，綜合分析確定影響奶綿羊產(chǎn)羔性狀的候選基因。

2 結(jié)果與分析

2.1 產(chǎn)羔數(shù)表型描述性統(tǒng)計(jì) 對(duì)100 只奶綿羊的產(chǎn)羔數(shù)進(jìn)行描述性統(tǒng)計(jì)，平均值是1.88，標(biāo)準(zhǔn)差是1.02，最小值是1，最大值是4，變異系數(shù)是54%。

2.2 質(zhì)控結(jié)果 對(duì)22 683 429 個(gè)SNP 進(jìn)行了基因分型，在進(jìn)行質(zhì)量控制時(shí)，SNP 缺失大于10% 的位點(diǎn)，剔除696 899 個(gè)；最小等位基因頻率小于1% 的位點(diǎn)，剔除0 個(gè)；不符合哈迪-溫伯格平衡檢驗(yàn)、值小于10的位點(diǎn)，剔除69 658 個(gè)。剩余135 個(gè)樣本和21 020 909個(gè)SNP 用于后續(xù)關(guān)聯(lián)分析。

2.3 樣本群體結(jié)構(gòu)分析 由圖1 可以看出，實(shí)驗(yàn)群體分成了3 個(gè)部分，說明該群體存在人工選育所引起的群體結(jié)構(gòu)差異。因此，在后續(xù)分析中，將主成分分析的前3個(gè)組分作為協(xié)變量加入到全基因組關(guān)聯(lián)分析模型中以消除群體分層的影響。

圖1 PCA 群體結(jié)構(gòu)圖

2.4 全基因組關(guān)聯(lián)分析 圖2 為產(chǎn)羔性狀的分位數(shù)-分位數(shù)（Quantile-Quantile）圖，結(jié)果顯示，大多數(shù)點(diǎn)都位于對(duì)角線上，表明觀察值與期望值的吻合度很好，末端觀測(cè)值與預(yù)期值的強(qiáng)烈偏差表明，與產(chǎn)羔性狀顯著相關(guān)的SNP 比預(yù)期的要多。圖3 為產(chǎn)羔性狀的曼哈頓圖，虛線代表全基因組潛在顯著閾值（4.76×10），實(shí)線代表全基因組顯著水平閾值（2.38×10）。由圖3 可知，沒有位點(diǎn)達(dá)到基因組顯著水平，有23 個(gè)SNPs 達(dá)到基因組潛在顯著水平。表1 為潛在關(guān)聯(lián)的SNP 位點(diǎn)，這些位點(diǎn)上下游基因有7 個(gè)。8 號(hào)染色體上1 個(gè)潛在顯著SNP（19 826 095）位于（未表征的基因）下游542 161 bp 處，（甘露糖苷酶Alpha 1A 級(jí)成員1）上游1 252 314 bp 處；15 號(hào)染色體上20個(gè)潛在顯著SNP 位于（甲基轉(zhuǎn)移酶15）下游約1.1 Mb 左右，（鉀電壓門控通道亞家族4）上游約0.6 Mb 左右；19 號(hào)染色體上1 個(gè)潛在顯著性SNP（18 205 638）位于（未表征的基因）下游343 407 bp 處，（CAMP 調(diào)節(jié)的磷酸蛋白21）上游983 325 bp 處；22 號(hào)染色體的1 個(gè)潛在顯著性SNP（19 684 838）位（細(xì)胞質(zhì)聚腺苷酸化元素結(jié)合蛋白3）內(nèi)含子區(qū)域。

表1 與產(chǎn)羔性狀潛在顯著關(guān)聯(lián)的SNPs

圖2 產(chǎn)羔性狀Quantile-Quantile 圖

圖3 產(chǎn)羔性狀全基因組關(guān)聯(lián)分析曼哈頓圖

3 討 論

GWAS 分析已成為研究禽畜重要經(jīng)濟(jì)性狀遺傳變異的有效方法。而有關(guān)戴瑞奶綿羊產(chǎn)羔性狀的GWAS研究未曾有報(bào)道。本研究對(duì)奶綿羊群體產(chǎn)羔性狀進(jìn)行全基因組關(guān)聯(lián)分析，共得到23 個(gè)潛在顯著SNP，分別位于綿羊8 號(hào)、15 號(hào)、19 號(hào)和22 號(hào)染色體上，通過候選基因注釋，共得到7 個(gè)候選基因，如下：

3.1基因 8 號(hào)染色體上SNP 19 826 095 鄰近和基因未曾有研究報(bào)道。又稱甘露糖苷酶Alpha 1A 級(jí)成員1。基因本體分析表明，主要在大分子糖基化的生物過程（GO:0043413）、鈣離子結(jié)合的分子功能（GO:0005509）和高爾基體膜的細(xì)胞成分（GO:0000139）中起作用，并參與N-聚糖生物合成（KEGG: ssc00510）和內(nèi)質(zhì)網(wǎng)中的蛋白質(zhì)合成（KEGG:ssc04141）。內(nèi)質(zhì)網(wǎng)是哺乳動(dòng)物真核細(xì)胞內(nèi)的重要膜性細(xì)胞器，分布廣，在蛋白質(zhì)合成、加工、轉(zhuǎn)運(yùn)、調(diào)節(jié)鈣離子動(dòng)態(tài)平衡和細(xì)胞信號(hào)轉(zhuǎn)導(dǎo)等方面發(fā)揮重要作用。母源mRNA 的翻譯對(duì)卵母細(xì)胞成熟和胚胎發(fā)育非常重要。翻譯的蛋白質(zhì)必須在內(nèi)質(zhì)網(wǎng)中正確折疊才能維持卵母細(xì)胞成熟和胚胎發(fā)育。而內(nèi)質(zhì)網(wǎng)腔大量未折疊或錯(cuò)誤折疊蛋白質(zhì)的積累會(huì)破壞內(nèi)質(zhì)網(wǎng)穩(wěn)態(tài)，引起內(nèi)質(zhì)網(wǎng)應(yīng)激。Lin 等研究表明蛋白質(zhì)糖基化或二硫鍵受到破壞也可能會(huì)引起內(nèi)質(zhì)網(wǎng)應(yīng)激，嚴(yán)重的應(yīng)激或長(zhǎng)時(shí)間的應(yīng)激會(huì)引發(fā)mRNA 剪接，促進(jìn)細(xì)胞凋亡。在卵母細(xì)胞中，這將導(dǎo)致發(fā)育停滯和表觀遺傳學(xué)變化。因此，可能通過影響蛋白質(zhì)糖基化或者內(nèi)質(zhì)網(wǎng)中蛋白質(zhì)合成而影響卵母細(xì)胞成熟發(fā)育。卵母細(xì)胞的發(fā)育狀況直接影響產(chǎn)羔表型，故可能是影響產(chǎn)羔性狀的候選基因。

3.2和基因 本研究發(fā)現(xiàn)，在15 號(hào)染色體上有10 個(gè)SNPs 位點(diǎn)與產(chǎn)羔數(shù)存在關(guān)聯(lián)，分布在0.6～1.1 Mb，這個(gè)區(qū)域內(nèi)有和基因。是甲基轉(zhuǎn)移酶15。Haute 等研究表明對(duì)于線粒體中蛋白質(zhì)的合成和氧化磷酸化過程十分重要。大多數(shù)哺乳動(dòng)物的線粒體都是母性遺傳，在紡錘體形成和卵母細(xì)胞受精等過程中卵母細(xì)胞的線粒體產(chǎn)生的ATP 是最直接的能量來源?？赡芡ㄟ^影響線粒體的生物發(fā)生而影響卵母細(xì)胞成熟或胚胎發(fā)育，進(jìn)而影響產(chǎn)羔性狀。因此，該基因可能是與產(chǎn)羔性狀相關(guān)的候選基因。是鉀離子電壓門控通道。此基因主要富集在皮質(zhì)醇的合成與分泌的信號(hào)通路上（KEGG：ssc04927）。皮質(zhì)醇通常稱為應(yīng)激激素，由腎上腺產(chǎn)生，主要在新陳代謝、免疫反應(yīng)以及生殖等方面發(fā)揮作用。Breen 等研究發(fā)現(xiàn)血漿皮質(zhì)醇濃度升高會(huì)抑制雌二醇的分泌，降低對(duì)該信號(hào)的神經(jīng)內(nèi)分泌反應(yīng)，影響卵泡發(fā)育和排卵等生物過程。有研究報(bào)道皮質(zhì)醇可能直接通過鉀通道調(diào)節(jié)來發(fā)揮它的某些神經(jīng)元作用。因此，可能通過調(diào)節(jié)皮質(zhì)醇的分泌進(jìn)而影響排卵。該基因可能是影響產(chǎn)羔性狀的候選基因。

3.3基因是cAMP 調(diào)節(jié)的磷蛋白，也是cAMP 依賴性蛋白激酶的磷蛋白底物，其主要分布在大腦特定的區(qū)域，該區(qū)域有豐富的D1 多巴胺受體。Ivkovic 等研究發(fā)現(xiàn)可能是D1 多巴胺受體刺激腺苷酸環(huán)化酶后引起細(xì)胞內(nèi)級(jí)聯(lián)反應(yīng)的第3 個(gè)信使。有研究表明多巴胺D1 受體的激活會(huì)增加的磷酸化水平。Tsou 等研究表明D1 多巴胺受體位于紋狀體神經(jīng)通路的神經(jīng)末梢，ARPP-21 的磷酸化可能介導(dǎo)了多巴胺對(duì)這些末梢的細(xì)胞內(nèi)作用。James 等研究認(rèn)為多巴胺神經(jīng)末梢不透明帶中的D1 多巴胺受體參與了促性腺激素釋放的控制，并且在生殖過程中可能具有生理功能。由此猜測(cè)可能通過影響多巴胺受體而影響與羊發(fā)情的相關(guān)通路，因此該基因可能是影響產(chǎn)羔性狀的重要候選基因。

3.4基因 SNP（19 684 838）位于（細(xì)胞質(zhì)聚腺苷酸化元素結(jié)合蛋白3）內(nèi)含子區(qū)域。屬于細(xì)胞質(zhì)聚腺苷酸化元素結(jié)合家族的成員。細(xì)胞質(zhì)聚腺苷酸化元素結(jié)合家族通過控制減數(shù)分裂細(xì)胞周期進(jìn)程來影響母體mRNA 的翻譯?；蜃⑨尫治霰砻鳎饕诤颂求w多肽合成的啟動(dòng)、激活、持久、抑制或終止的分子功能（GO:0045182）和卵母細(xì)胞減數(shù)分裂的信號(hào)通路中（KEGG：ssc04114）發(fā)揮作用。Wang 等的研究表明受到抑制后，牛卵丘細(xì)胞具有異常的細(xì)胞周期，S 期和G2/M 期的細(xì)胞數(shù)量增加，細(xì)胞發(fā)生明顯增殖。因此可能通過影響卵丘細(xì)胞減數(shù)分裂的細(xì)胞周期來影響產(chǎn)羔，該基因也可能是影響產(chǎn)羔性狀的重要候選基因。

通過對(duì)以上候選基因的功能分析，和最可能是影響產(chǎn)羔性狀的重要候選基因。后期仍需要進(jìn)一步對(duì)挖掘到的候選基因進(jìn)行功能驗(yàn)證。

4 結(jié) 論

本研究發(fā)現(xiàn)與產(chǎn)羔性狀潛在關(guān)聯(lián)的23 個(gè)SNP，注釋出7 個(gè)候選基因。通過基因功能分析可知，和是最有可能影響產(chǎn)羔性狀的候選基因。本研究結(jié)果為深入揭示奶綿羊產(chǎn)羔性狀的遺傳效應(yīng)提供了參考，為奶綿羊新品種培育提供了新的研究線索。