利用兩種統(tǒng)計模型對中國肉用西門塔爾牛屠宰性狀的全基因組關聯(lián)分析

常天鵬，夏江威，寶金山，金生云，朱 波，徐凌洋，陳 燕，張路培，高 雪，李俊雅，高會江*

(1.中國農(nóng)業(yè)科學院北京畜牧獸醫(yī)研究所，北京 100193； 2.內蒙古錫林郭勒盟烏拉蓋管理區(qū)獸醫(yī)局，烏拉蓋 026321)

屠宰性狀是肉牛產(chǎn)業(yè)中最為重要的一類性狀，直接影響著牛肉產(chǎn)量和經(jīng)濟效益，因此，對于屠宰性狀的選育是肉牛育種者的重點研究方向。借助于現(xiàn)代分子生物學技術，分子標記輔助選擇(MAS)顯著縮短了選育周期，節(jié)省了時間和成本投入，成為當前育種工作的有效方法[1]。而實現(xiàn)肉牛屠宰性狀MAS的關鍵步驟就是目標性狀的數(shù)量性狀基因座(QTLs)及候選基因定位[2]。目前，關于世界上主要肉牛品種屠宰性狀的候選基因定位工作已經(jīng)有了許多的研究報道，其中對胴體重的研究較多。A.Takasuga等[3]通過對日本黑牛進行QTLs定位，在6號染色體38 Mb和14號染色體41.7和49.8 Mb均檢測到了影響胴體重的QTLs；A. K.Lindholm-Perry等[4]通過對一個包括14個品種血統(tǒng)的雜交牛群體進行基因定位，在6號染色體38 Mb鑒定到影響胴體重的LAP3、NCAPG、LCORL等候選基因；A.G.Doran等[5]借助于GWAS方法對荷斯坦牛的胴體性狀進行分析，找到了一系列影響胴體性狀的QTLs、候選基因及信號通路。

在GWAS時代到來之前，研究者主要利用QTL定位策略定位候選基因。QTL定位常用的試驗群體類型有BC、F2及RIL群體等，這些群體是精心設計的，交配類型也是明確的，其目的在于使其基因型分離以增加表型的豐富性，而且還方便推測偽標記的基因型[6]。QTL定位比較經(jīng)典的方法有區(qū)間定位法(IM)、復合區(qū)間定位法(CIM)[7-8]等。但是由于標記密度低，達不到精確定位的層次。隨著高密度芯片技術的發(fā)展，GWAS以其較高的精確性和效率，成為當前探索復雜性狀遺傳機理的最有效的方法之一。GWAS選用的群體一般是自然群體，基于歷史重組造成標記與QTL的連鎖不平衡來定位真實QTL的位置。然而由于樣本量要求大，自然群體經(jīng)常會有群體分層現(xiàn)象，從而導致假陽性過高。目前，常用的GWAS方法有線性混合模型(LMM)和貝葉斯方法[9-10]。通過QTL定位與GWAS方法的比較，不難發(fā)現(xiàn)，由于QTL定位精心的試驗設計，能夠很好地避免GWAS存在的問題，因此先前在QTL定位積累的統(tǒng)計方法對GWAS研究具有很好的借鑒意義。

本試驗擬嘗試兩種不同的統(tǒng)計模型：一種是線性混合模型(LMM)，另一種是結合復合區(qū)間作圖策略的線性混合模型(CIM-LMM)，基于Illumina BovineHD (770K)芯片分型獲得基因型數(shù)據(jù)，對中國肉用西門塔爾牛的胴體重和骨重性狀進行關聯(lián)分析，旨在篩選影響目標性狀的顯著SNPs并挖掘相關的功能基因；同時，比較兩個方法獲得的結果，為GWAS的研究方法提供新策略。

1 材料與方法

1.1 試驗材料

本研究中試驗動物來自中國農(nóng)業(yè)科學院北京畜牧獸醫(yī)研究所在內蒙古錫林郭勒盟烏拉蓋管理區(qū)構建的中國肉用西門塔爾牛資源群體。該群體自2008年開始建立并記錄各個階段表型數(shù)據(jù)。為了便于數(shù)據(jù)的測量和采集，將5～9月齡的牛運送至北京金維福仁清真食品有限公司肉牛養(yǎng)殖場，在該場期間按照統(tǒng)一的飼養(yǎng)管理方法集中育肥，每3個月測量一次體尺和體重數(shù)據(jù)，當飼養(yǎng)至18～24月齡時進行分批屠宰分割。截止2016年12月底，共收集到1 301頭肉牛的表型數(shù)據(jù)和基于Illuminate BovineHD (770K)芯片分型的基因型數(shù)據(jù)。

選擇的目標性狀為胴體重和骨重，胴體重是指牛經(jīng)過宰殺放血后，除去皮、頭、蹄、尾、內臟以及生殖器(母牛去除乳房)后的軀體部分重量；骨重是指胴體經(jīng)過標準流程分割后剩余所有骨頭的總重量。屠宰以及胴體分割流程完全按照國家標準GB/T 27643-2011進行。表1顯示了胴體重和骨重兩個性狀的描述性統(tǒng)計量，包括表型數(shù)據(jù)的記錄數(shù)、最小值、最大值、平均值、標準差及遺傳力。

表1胴體重、骨重性狀的描述性統(tǒng)計量

Table1Thedescriptivestatisticsofcarcassweightandboneweight

性狀Trait記錄數(shù)Number最小值Min.最大值Max.平均值Mean標準差SD遺傳力*h2胴體重/kgCarcassweight1217162.6486.0271.345.630.38骨重/kgBoneweight121720.2080.0040.676.5150.41

*性狀遺傳力由R 環(huán)境下asreml包估計

*Heritabilities were estimated by asreml package in R environment

1.2 基因分型和質量控制

采集1 301個個體血樣，利用TIANGEN基因組DNA提取試劑盒(天根生化科技有限公司，北京,中國)從血樣中提取出基因組DNA，并利用Illumina BovineHD (770K)芯片進行基因分型。然后采用PLINK v1.07軟件進行質量控制，質量控制標準為：(1)SNPs檢出率>90%；(2)最小等位基因頻率>10%；(3)哈代-溫伯格平衡檢驗(P>10-6)；(4)個體基因型缺失率<0.1。經(jīng)過質控之后，共有1 217個個體和671 990個SNPs用于關聯(lián)分析。

1.3 統(tǒng)計模型

利用兩種統(tǒng)計模型：線性混合模型(LMM)和復合區(qū)間定位-線性混合模型(CIM-LMM)進行目標性狀的關聯(lián)分析。

LMM為GWAS最常用的統(tǒng)計方法：

y=Xβ+Zkγk+ξ+ε

(1)

其中,X表示固定效應系數(shù)矩陣，β是固定效應向量；Zk表示第k個SNP的指示變量，γk表示第k個SNP的效應值(作固定效應處理)；ξ～N(0,Kφ2)稱作微效多基因效應，K表示基于標記推斷的親緣關系矩陣，φ2表示多基因效應的方差；ε表示隨機殘差并且服從多元正態(tài)分布N(0,Iσ2)。該方法最早是由J.M.Yu等[10]提出，也稱作常用的Q+K模型，本試驗中該方法由Tassel 5.0軟件中的LMM實現(xiàn)。

CIM-LMM的統(tǒng)計模型：

y=Xβ+Zk-1γk-1+Zkγk+Zk+1γk+1+ξ+ε

(2)

在統(tǒng)計檢驗方面，這兩種方法都采用P<0.05作為篩選顯著SNPs的標準。為了解決多重假設檢驗引起的假陽性問題，本研究采用Bonferroni法進行P值校正[12]。

1.4 群體結構分析

在進行GWAS分析時，群體結構的影響不可忽略。本研究采用主成分分析法(Principal component analysis，PCA)檢測群體結構[13]。通過PCA獲得基因型數(shù)據(jù)特征值和特征向量，然后選擇適當?shù)奶卣飨蛄孔鳛閰f(xié)變量加入關聯(lián)分析統(tǒng)計模型中，從而實現(xiàn)對群體結構的校正。本研究中PCA通過R程序中的princomp函數(shù)實現(xiàn)。

1.5 固定效應顯著性檢驗

固定效應的校正可以降低數(shù)據(jù)資料的不平衡性和關聯(lián)分析的假陽性率。對于固定效應的分析，要把盡可能的影響因素加入到校正模型中，然后對各因素做顯著性分析，將對表型值有顯著影響的因素放入統(tǒng)計模型，不顯著的影響因素去掉。本研究根據(jù)資源群體的實際情況，考慮的固定因素為出生年、季節(jié)、出生場效應，并將進場重、育肥天數(shù)、群體結構作為協(xié)變量。通過R程序中的GLM函數(shù)實現(xiàn)固定效應的顯著性檢驗，其模型為y=u+year+season+farm+weight+fatten_day+pca+e，其中y為性狀表型值，u為群體均值，year為出生年效應，season為出生季節(jié)效應，farm出生場效應，weight為進場重協(xié)變量，fatten_day為育肥天數(shù)協(xié)變量，pca為群體結構協(xié)變量，e為隨機殘差。

1.6 基因注釋

利用生物信息學網(wǎng)站Ensemble中的BioMart模塊將檢測到的顯著SNPs比對到牛的基因組(BostaurusUMD 3.1)中，依據(jù)SNPs的物理位置在其上下游尋找候選基因。然后通過NCBI網(wǎng)站的Gene數(shù)據(jù)庫查找相關基因的生物學功能，并結合前人報道對候選基因進行功能注釋。

2 結 果

2.1 群體結構及固定效應檢驗

圖1為基于前兩個主成分將個體聚類的群體結構圖。從圖中明顯可以看出，個體間存在很大程度的群體分層。絕大多數(shù)個體集中于右下角區(qū)域，另外4個區(qū)域有少量個體分布。所以選取前兩個主成分作為協(xié)變量來消除群體分層對關聯(lián)分析的影響。在固定效應的顯著性檢驗中，年效應、進場重及育肥天數(shù)協(xié)變量對表型值有顯著影響。因此，選擇將這些顯著的固定效應及協(xié)變量一起放入統(tǒng)計模型進行關聯(lián)分析。

圖1 根據(jù)PCA繪制群體結構圖Fig.1 Population structure identified by principal components analysis

2.2 兩種模型的QQ圖分析

利用LMM和CIM-LMM對目標性狀的GWAS結果的P值作QQ圖，橫坐標表示期望P值的負對數(shù)，縱坐標表示實際觀察P值的負對數(shù)，見圖2?？梢钥闯觯褂肔MM分析得到的P值，觀察值與期望值吻合度很好，大多數(shù)的點都位于對角線上，這表明此模型較好的校正了系統(tǒng)偏差。由CIM-LMM分析結果的得到的QQ圖前端顯示了一定的拖尾現(xiàn)象，其原因在于將目標標記k的鄰近標記做隨機效應處理，導致許多微效或無效標記的效應值壓縮為零，而統(tǒng)計檢驗時P值接近或等于1，而后端觀察值與期望值基本吻合。在統(tǒng)計顯著性方面，CIM-LMM方法得到的-lg(P)值比LMM方法高2～5個數(shù)量級，而且識別出較多偏離期望值的SNPs。

2.3 兩種模型鑒定的SNPs和候選基因

圖3展示了通過LMM和CIM-LMM對目標性狀GWAS分析的Manhattan圖。對于這兩個性狀，以Bonferroni法矯正后的P值作為閾值，分別在全基因組范圍內共找到8和7個顯著SNPs，其中2個SNPs與胴體重及骨重都有關聯(lián)。通過Ensembl網(wǎng)站BioMart模塊鑒定了與這些SNPs密切相關的11個基因。對于胴體重，LMM方法共檢測到3個顯著相關的SNPs，分別位于第5、14、17號染色體上。這3個SNPs分別鄰近或坐落于C12ORF74、RIMS2、BT.88981基因。CIM-LMM方法共檢測到了8個顯著相關的SNPs，其中包含了LMM方法識別的3個SNPs，剩余的5個SNPs分別位于第3、10、14、16號染色體上，鄰近或坐落于PBX1、GCNT4、ALDH1A2、FAM84B、DUSP10基因。對于骨重，LMM方法共檢測到4個顯著相關的SNPs，分別位于第5、6、14號染色體上，鄰近或坐落于C12ORF74、LCORL、RIMS2基因。CIM-LMM方法共找到了7個顯著相關的SNPs，其中包含了LMM方法檢測的4個SNPs，剩余的3個SNPs分別位于第6、14號染色體上，鄰近于WDFY3、FER1L6基因。具體如表2和表3所示。

3 討 論

3.1 群體結構的分層現(xiàn)象

從圖1可以看出，本研究所用的中國肉用西門塔爾牛群體呈現(xiàn)一定的群體分層。所謂群體分層，主要由于自然選擇、人工選擇、遺傳漂變等諸多因素導致的群體內出現(xiàn)了等位基因頻率不同的亞群體。它會導致一些非原因等位基因同真實QTLs形成連鎖不平衡而表現(xiàn)出與目標性狀關聯(lián),從而導致偽關聯(lián)或假陽性[14]。若不對群體分層加以校正，那么群體分層效應會被誤認為是QTLs效應，降低關聯(lián)分析準確性。因此，本研究選取PCA結果中前兩個主成分作為協(xié)變量來校正群體分層，從而避免由于群體分層引起的假陽性問題。

a、c.利用LMM分別對胴體重、骨重GWAS結果的QQ圖;b、d.利用CIM-LMM對胴體重、骨重GWAS結果的QQ圖a，c. QQ plots for the GWAS results of carcass weight and bone weight by LMM, respectively; b，d. QQ plots for the GWAS results of carcass weight and bone weight by CIM-LMM, respectively圖2 兩種模型對目標性狀分析結果的QQ圖Fig.2 QQ plots for the analyzing results of target traits by two statistic models

a、c. 利用LMM分別對胴體重、骨重GWAS結果作的曼哈頓圖;b、d. 利用CIM-LMM對胴體重、骨重GWAS結果作的曼哈頓圖a，c. Manhattan plots for the GWAS result of carcass weight and bone weight by LMM, respectively; b，d. Manhattan plots for the GWAS result of carcass weight and bone weight by CIM-LMM, respectively圖3 兩種模型對目標性狀分析結果的曼哈頓圖Fig.3 Manhattan plots for the analyzing results of target traits by two statistic models

表2兩種模型對胴體重的GWAS分析結果

Table2TheGWASresultofcarcassweightbytwostatisticmodels

模型Model多態(tài)位點SNP染色體Chromosome物理位置Position距離*/bpDistance候選基因Gene校正P值CorrectedP-valueLMMBovineHD0500006528522558100161123C12ORF743.07E-19BovineHD14000174551462769117withinRIMS22.58E-17BovineHD170002134017730075223479BT.889811.06E-10CIM-LMMBovineHD030000147634776413249130PBX11.90E-08BovineHD0500006528522558100161123C12ORF744.64E-22BovineHD1000001969106318056114132GCNT47.38E-08BovineHD10000156321052349382withinALDH1A23.40E-08BovineHD14000043711415588175121944FAM84B2.39E-08BovineHD14000174551462769117withinRIMS21.85E-21BovineHD16000071661625759434137313DUSP101.76E-08BovineHD170002134017730075223479BT.889811.51E-11

*表示SNP與鄰近基因的物理距離，單位為堿基對(bp),“within”表示落入基因內，下表同

*.Mean the physical distance between SNP and the nearest gene and the units are bp, “within” indicate that the SNP located in the nearest genes. The same as below

表3兩種模型對骨重的GWAS分析結果

Table3TheGWASresultofboneweightbytwostatisticmodels

模型Model多態(tài)位點SNP染色體Chromosome物理位置Position距離/bpDistance基因Gene校正P值CorrectedP-valueLMMBovineHD0500006528522558100161123C12ORF741.03E-11BovineHD0600010952639990876998873LCORL2.95E-08BovineHD06000109566399978801005877LCORL3.91E-08BovineHD14000174551462769117withinRIMS28.00E-11CIM-LMMBovineHD0500006528522558100161123C12ORF743.31E-13BovineHD0600010952639990876998873LCORL4.14E-08BovineHD06000109566399978801005877LCORL2.72E-08BovineHD0600028304610151470230400WDFY36.63E-09BovineHD1400005015141760061617032FER1L61.88E-12BovineHD1400005019141761244428860FER1L62.05E-12BovineHD14000174551462769117withinRIMS23.71E-13

3.2 兩種模型分析結果比較

兩種模型GWAS結果的QQ圖顯示，LMM方法結果的理論值與觀察值基本吻合，表明該模型較好的校正了群體分層；然而顯著SNPs檢測效力一般，對于兩個性狀，在全基因組范圍內分別僅發(fā)現(xiàn)3、4個SNPs與之關聯(lián)。CIM-LMM方法利用基于符合區(qū)間定位(CIM)的策略，在目標標記的兩側加入側翼標記并做隨機效應處理，考慮了標記之間的連鎖關系同時壓縮無效標記的效應；結果除了檢測到LMM發(fā)現(xiàn)的所有SNPs，還檢測到更多的SNPs，同時各SNPs的檢驗P值更加顯著，較之于LMM方法更有優(yōu)勢。

3.3 候選基因的功能分析

本研究共發(fā)現(xiàn)11個不同的候選基因，分布于3、5、6、10、14、16、17號染色體上，且多集中于6號和14號染色體上，一些候選基因已在其他牛種或物種上被相繼研究。其中，PBX1基因位于BTA3上，據(jù)報道，它與人類的骨密度有關，并且被認為是影響印第安人2型糖尿病的1個候選基因[15]。LCORL、WDFY3基因位于BTA6上，LCORL已經(jīng)被證實與人類的骨骼骨架尺寸和馬的肩胛骨高度有關[16]，而且也被發(fā)現(xiàn)它是影響牛飼料攝入量和生長的調控基因[4]；WDFY3已經(jīng)被確定與仔豬的出生重有關，可能參與豬的胎兒發(fā)育過程[17]。GCNT4、ALDH1A2基因位于BTA10上，研究表明，牛的GCNT4基因參與脂肪酸從頭合成的途徑和黏液生物合成[18]；ALDH1A2基因在合成視黃酸的過程中起著關鍵作用，視黃酸是維生素A代謝的中間產(chǎn)物，而維生素A又主要影響動物的身體生長和骨骼發(fā)育[19]。FAM84B、RIMS2、FER1L6基因都位于BTA14上，研究發(fā)現(xiàn)，F(xiàn)AM84B可能是影響人類食道鱗狀細胞癌的候選基因[20]，在牛上該基因被發(fā)現(xiàn)與荷斯坦公牛精子活力有關系[21]；RIMS2基因是人類中樞神經(jīng)系統(tǒng)突觸膜胞外分泌的調控基因，并且與人類退化性腰脊柱側凸(DLS)的患病風險有關[22-23]；而韓牛的一項研究表明，F(xiàn)ER1L6基因與韓牛牛肉的大理石花紋形成密切關聯(lián)，并認為該基因是影響牛肉質量的候選基因[24]。

4 結 論

本研究分別使用LMM和CIM-LMM兩種統(tǒng)計模型，對中國肉用西門塔爾牛的胴體重和骨重性狀進行了全基因組關聯(lián)分析。分別鑒定到8、7個顯著相關的SNPs,其中有2個SNPs與這兩個性狀都有關聯(lián)，并利用NCBI及Ensemble網(wǎng)站查找到影響目標性狀的11個候選基因；同時還探討了GCNT4、ALDH1A2、LCORL和WDFY3等基因的功能。此外，我們也對2種統(tǒng)計模型的分析結果進行比較，結果顯示，對于本試驗群體CIM-LMM方法較之于LMM方法具有更高的統(tǒng)計效力。本研究為解析中國肉用西門塔爾牛屠宰性狀的遺傳機理做了探索，為后期的驗證研究提供了方向，也為GWAS方法研究提供了新思路。

參考文獻(References)：

[1] DEKKERS J C M. Commercial application of marker-and gene-assisted selection in livestock: strategies and lessons[J].JAnimSci, 2004, 82 (E-Suppl): E313-E328.

[2] DUNHAM R A. Gene mapping, quantitative trait locus mapping and marker-assisted selection[M]//DUNHAM R A. Aquaculture and Fisheries Biotechnology: Genetic Approaches. Wallingford, UK: CABI Publishing, 2004.

[3] TAKASUGA A, WATANABE T, MIZOGUCHI Y, et al. Identification of bovine QTL for growth and carcass traits in Japanese Black cattle by replication and identical-by-descent mapping[J].MammGenome, 2007, 18(2): 125-136.

[4] LINDHOLM-PERRY A K, SEXTEN A K, KUEHN L A, et al. Association, effects and validation of polymorphisms within the NCAPG - LCORL locus located on BTA6 with feed intake, gain, meat and carcass traits in beef cattle[J].BMCGenet, 2011, 12(1): 103.

[5] DORAN A G, BERRY D P, CREEVEY C J. Whole genome association study identifies regions of the bovine genome and biological pathways involved in carcass trait performance in Holstein-Friesian cattle[J].BMCGenomics,2014, 15: 837.

[6] 魏巨龍. 混合線性模型解析數(shù)量性狀遺傳基礎的研究[D]. 北京: 中國農(nóng)業(yè)大學, 2016.

WEI J L. Application of Linear Mixed Model (LMM) to dissect genetic basis of quantitative traits[D]. Beijing: China Agricultural University, 2016. (in Chinese)

[7] ZENG Z B. Precision mapping of quantitative trait loci[J].Genetics, 1994, 136(4): 1457-1468.

[8] LANDER E S, BOTSTEIN D. Mapping Mendelian factors underlying quantitative traits using RFLP linkage maps[J].Genetics, 1989, 121(1): 185-199.

[9] FERNANDO R L, GARRICK D. Bayesian methods applied to GWAS[M]//GONDRO C, VAN DER WERF J, HAYES B. Genome-Wide Association Studies and Genomic Prediction. Totowa, NJ: Humana Press, 2013: 237-274.

[10] YU J M, PRESSOIR G, BRIGGS W H, et al. A unified mixed-model method for association mapping that accounts for multiple levels of relatedness[J].NatGenet, 2006, 38(2): 203-208.

[11] WEI J L, XU S Z. A random-model approach to QTL mapping in multiparent advanced generation intercross (MAGIC) populations[J].Genetics, 2016, 202(2): 471-486.

[12] JOHNSON R C, NELSON G W, TROYER J L, et al. Accounting for multiple comparisons in a genome-wide association study (GWAS)[J].BMCGenomics, 2010, 11: 724.

[13] PRICE A L, PATTERSON N J, PLENGE R M, et al. Principal components analysis corrects for stratification in genome-wide association studies[J].NatGenet, 2006, 38(8): 904-909.

[14] WANG D, SUN Y, STANG P, et al. Comparison of methods for correcting population stratification in a genome-wide association study of rheumatoid arthritis: principal-component analysis versus multidimensional scaling[J].BMCProc, 2009, 3(Suppl 7): S109.

[15] THAMEEM F, WOLFORD J K, BOGARDUS C, et al. Analysis ofPBX1 as a candidate gene for type 2 diabetes mellitus in Pima Indians[J].BiochimBiophysActa, 2001, 1518(1-2): 215-220.

[16] SORANZO N, RIVADENEIRA F, CHINAPPEN-HORSLEY U, et al. Meta-analysis of genome-wide scans for human adult stature identifies novel loci and associations with measures of skeletal frame size[J].PLoSGenet, 2009, 5(4): e1000445.

[17] ZHANG L F, ZHOU X, MICHAL J J, et al. Genome wide screening of candidate genes for improving piglet birth weight using high and low estimated breeding value populations[J].IntJBiolSci, 2014, 10(3): 236-244.

[18] WU J H Y, LEMAITRE R N, MANICHAIKUL A, et al. Genome-wide association study identifies novel loci associated with concentrations of four plasma phospholipid fatty acids in the de novo lipogenesis pathway: Results from the cohorts for heart and aging research in genomic epidemiology (CHARGE) consortium[J].CircCardiovascGenet, 2013, 6(2): 171-183.

[19] OHOKA Y, YOKOTA-NAKATSUMA A, MAEDA N, et al. Retinoic acid and GM-CSF coordinately induce retinal dehydrogenase 2 (RALDH2) expression through cooperation between the RAR/RXR complex and Sp1 in dendritic cells[J].PLoSOne, 2014, 9(5): e96512.

[20] CHENG C X, CUI H Y, ZHANG L, et al. Genomic analyses reveal FAM84B and the NOTCH pathway are associated with the progression of esophageal squamous cell carcinoma[J].GigaScience, 2016, 5: 1.

[21] HERING D M, OLENSKI K, KAMINSKI S. Genome-wide association study for poor sperm motility in Holstein-Friesian bulls[J].AnimReprodSci, 2014, 146(3-4): 89-97.

[22] KIM K T, LEE J S, LEE B W, et al. Association between regulating synaptic membrane exocytosis 2 gene polymorphisms and degenerative lumbar scoliosis[J].BiomedRep, 2013, 1(4): 619-623.

[23] KAESER P S,DENG L B, FAN M M, et al. RIM genes differentially contribute to organizing presynaptic release sites[J].ProcNatlAcadSciUSA, 2012, 109(29): 11830-11835.

[24] RYU J, LEE C. Identification of contemporary selection signatures using composite log likelihood and their associations with marbling score in Korean cattle[J].AnimGenet, 2015, 45(6): 765-770.