中國肉用西門塔爾牛生長曲線參數(shù)的全基因組關(guān)聯(lián)分析

段星海，安炳星，杜麗麗，常天鵬，梁 忙，楊柏高，高會(huì)江*，俄廣鑫*

(1.西南大學(xué)動(dòng)物科學(xué)技術(shù)學(xué)院，重慶 400715；2.中國農(nóng)業(yè)科學(xué)院北京畜牧獸醫(yī)研究所，北京 100193)

縱向性狀指的是表型值隨著時(shí)間(生命周期、年齡、胎次等)或其他定量因素(生產(chǎn)水平、生理狀態(tài)和環(huán)境條件等)變化的性狀[1]。動(dòng)物的重要經(jīng)濟(jì)性狀如產(chǎn)奶性狀、育肥性狀和產(chǎn)蛋性狀等都屬于縱向性狀[2-4]。與一般的只有單個(gè)記錄的表型性狀相比較而言，縱向性狀更好地描述了畜禽的生長和生產(chǎn)規(guī)律。目前的研究通常擬合生長曲線來反映個(gè)體發(fā)育速度與成熟速率之間的相互關(guān)系[5-7]，這是反映畜禽生長發(fā)育規(guī)律的主要方法之一[8-9]。這些模型一般通過幾個(gè)參數(shù)來描述家畜的體重變化規(guī)律，例如成熟體重(A)、達(dá)到最大生長率的時(shí)間(b)、成熟率(K)等，這些參數(shù)度量個(gè)體整個(gè)發(fā)育時(shí)期的特性[10-11]。對于特定的群體，需通過擬合度來選擇最佳模型擬合生長曲線[12]。已有的研究中，普遍假定參數(shù)(如A和K)作混合線性模型的表型，使用全基因組關(guān)聯(lián)分析(genome wide association study,GWAS) 定位影響生長曲線的數(shù)量性狀基因座(quantitative trait loci，QTL)。全基因組關(guān)聯(lián)分析是掃描全基因組范圍內(nèi)的單核苷酸多態(tài)性(single nucleotide polymorphism,SNPs),以識(shí)別與目標(biāo)性狀變異關(guān)聯(lián)的遺傳變異的分析方法[13]。Das等[14]提出了一種基于隨機(jī)回歸的方法，該模型能夠利用多個(gè)時(shí)間點(diǎn)之間的非加性效應(yīng)和特定的協(xié)方差函數(shù)來描述SNP效應(yīng)標(biāo)記隨時(shí)間的變化。盡管這些基于隨機(jī)回歸的模型是目前評估GWAS縱向性狀最復(fù)雜和最靈活的工具[15-16]，但它不能直接提供動(dòng)物育種項(xiàng)目中通常需要的成熟重量(A)和成熟率(K)等生物學(xué)可解釋的參數(shù)。

自Fitzhugh[5]提出對于動(dòng)物育種中生長曲線的研究以來，肉牛育種中已有大量的研究考慮參數(shù)(主要是A、b和K)和QTL的遺傳相關(guān)，這些遺傳相關(guān)性可能歸因于對多個(gè)參數(shù)具有多效性或緊密連鎖的SNP。Crispim等[17]利用5種常用的生長曲線模型對婆羅門牛的體重?cái)?shù)據(jù)進(jìn)行生長曲線的擬合認(rèn)為，Brody 模型的擬合效果最好，并且發(fā)現(xiàn)，RAB28、BTG1、IL2、APEX2等基因與婆羅門牛的胎兒和肌肉發(fā)育有關(guān)[18-21]；Soares等[22]通過5種 常用的生長曲線模型對婆羅門牛的陰囊周長數(shù)據(jù)進(jìn)行擬合發(fā)現(xiàn)，von Bertalanffy模型擬合效果最好，并且認(rèn)為LYN、CHD7、SOX9、IGF-1等可以作為婆羅門牛繁殖性狀的候選基因[23-26]。中國西門塔爾牛群體適應(yīng)性好，抗病力強(qiáng)，耐粗飼，在我國多種生態(tài)條件下均能表現(xiàn)出良好的生產(chǎn)性能，肉用西門塔爾牛占中國肉牛市場的70%左右[27]，然而尚無對中國西門塔爾牛的多時(shí)間點(diǎn)體重性狀進(jìn)行生長曲線的擬合研究。因此，本研究首先通過3種模型來模擬中國肉用西門塔爾牛的體重生長曲線，評價(jià)擬合度來選擇擬合效果最佳的模型。利用最佳模型估計(jì)的參數(shù)作為GWAS中的表型，鑒定與中國西門塔爾牛生長發(fā)育性狀顯著相關(guān)的SNPs，為中國西門塔爾?；蚪M育種提供可用的突變信息。

1 材料與方法

1.1 試驗(yàn)動(dòng)物表型的描述性統(tǒng)計(jì)

試驗(yàn)動(dòng)物來自中國農(nóng)業(yè)科學(xué)院北京畜牧獸醫(yī)研究所在內(nèi)蒙古自治區(qū)錫林郭勒盟烏拉蓋管理區(qū)構(gòu)建的中國肉用西門塔爾牛資源群體。本研究挑選了同時(shí)包含0、6、12、18月齡體重測量數(shù)據(jù)的個(gè)體，共808頭。表1為表型數(shù)據(jù)的描述性統(tǒng)計(jì)信息。

表1 中國肉用西門塔爾牛4個(gè)發(fā)育時(shí)期體重的描述性統(tǒng)計(jì)Table 1 The descriptive statistics of body weight at 4 growth stages of Chinese Simmental beef cattle kg

1.2 基因分型和質(zhì)量控制

采集所有個(gè)體血樣，利用TIANGEN基因組DNA提取試劑盒(天根生化科技有限公司，北京，中國) 從血樣中提取基因組DNA，利用Illumina BovineHD(770K)芯片進(jìn)行基因分型。然后采用PLINK v1.07軟件進(jìn)行質(zhì)量控制，質(zhì)量控制標(biāo)準(zhǔn)為：SNPs檢出率>90%；最小等位基因頻率>10%；哈代-溫伯格平衡檢驗(yàn)(P>10-6)；個(gè)體基因型缺失率<0.1%。質(zhì)控后共有808個(gè)個(gè)體和671 991個(gè) SNPs用于關(guān)聯(lián)分析。

1.3 生長曲線的擬合

利用3種最常用的非線性模型(表2)來描述動(dòng)物的生長曲線，使用SAS 9.4 軟件中的NLIN過程對整理后的數(shù)據(jù)進(jìn)行曲線模型方程參數(shù)的最優(yōu)估計(jì)，采用Gauss-Newton的方法迭代求解[28]，使估計(jì)的參數(shù)殘差平方和最小，收斂標(biāo)準(zhǔn)為10-8。

表2 生長曲線模型Table 2 Growth curve model

R2用來評價(jià)生長曲線模型擬合程度的高低[29]，即擬合度：

選擇擬合效果最好的模型，再次通過SAS 9.4 軟件中的NLIN過程估計(jì)每個(gè)個(gè)體的參數(shù)A、b、K。

1.4 全基因組關(guān)聯(lián)分析

選擇最合適體重?cái)?shù)據(jù)的生長曲線模型，將獲得的生長曲線參數(shù)A、b、K用于全基因組關(guān)聯(lián)分析(GWAS)，分析之前首先使用線性模型對固定因素(場和出生年/月/日)進(jìn)行校正，校正模型如下：

y=Xβ+y*

其中，y是觀測表型(A、b、K估計(jì)值向量)，β代表固定效應(yīng)向量(場和出生年/月/日)，X為對應(yīng)的設(shè)計(jì)矩陣，y*用于后續(xù)的關(guān)聯(lián)分析模型中。

然后使用R v3.4.2中的基因組關(guān)聯(lián)與預(yù)測集成工具(GAPIT)包(http://www.maizegenetics.net/gapit)進(jìn)行主成分分析(PCA)和親緣關(guān)系矩陣的計(jì)算[30]，GWAS模型如下：

y*=Xβ+Ws+Zμ+е

其中，β代表固定效應(yīng)向量(PCA前3個(gè)特征向量)，X為固定效應(yīng)向量的關(guān)聯(lián)矩陣；W是SNP基因型指示載體，3種基因型AA、AB、BB被編碼為0、1、2；s代表SNP效應(yīng)向量；Z為多基因效應(yīng)向量的關(guān)聯(lián)矩陣；μ代表微效多基因效應(yīng)向量，服從正態(tài)分布(0,Kσ2)；e表示隨機(jī)殘差，并且服從正態(tài)分布N(0,Iσ2)。

利用錯(cuò)誤發(fā)現(xiàn)率(false discovery rate,FDR)多重檢驗(yàn)方法對GWAS分析得到的P值進(jìn)行檢驗(yàn)，矯正P值的公式為：

P=FDR×n/m

其中，n代表小于FDR值的SNP個(gè)數(shù)，m代表SNP總個(gè)數(shù)，F(xiàn)DR設(shè)定值為0.05。

1.5 基因功能注釋

依據(jù)生物信息學(xué)網(wǎng)站Ensembl中的BioMart模塊將檢測到的顯著SNPs比對到牛的基因組(BostaurusUMD 3.1) (http://www.animalgenome.org)中，依據(jù)SNPs的物理位置在其上、下游尋找相關(guān)候選基因。然后通過NCBI網(wǎng)站(https://www.ncbi.nlm.nih.gov/)的Gene數(shù)據(jù)庫查找相關(guān)基因的生物學(xué)功能，并結(jié)合相關(guān)文獻(xiàn)報(bào)道對候選基因進(jìn)行功能注釋。

2 結(jié) 果

2.1 生長曲線的擬合

表3給出了3種生長曲線模型的擬合參數(shù)結(jié)果和R2，R2最高的為Gompertz模型，達(dá)到了0.954，Logistic模型和Brody模型的R2均為0.951；Gompertz、Logistic和Brody模型參數(shù)A的值依次為617.900、551.000和1 458.500，參數(shù)b的值依次為2.740、9.304和0.976，參數(shù)K的值依次為0.153、0.273和0.024；圖1顯示了3種模型曲線與體重均值曲線的吻合程度，可以看出Gompertz模型和體重均值的曲線基本吻合，而其他兩種模型曲線有所偏差。表4給出了Gompertz模型參數(shù)的描述性統(tǒng)計(jì)，結(jié)果顯示，成熟體重(A) 的最大值為985.10，最小值為388.94，平均成熟體重為629.11，標(biāo)準(zhǔn)差為88.71；達(dá)到最大生長率的時(shí)間(b)最大值與最小值分別為9.17和2.24，平均值為2.82，標(biāo)準(zhǔn)差為0.43；成熟率(K) 的最大值和最小值分別為0.40和0.08，平均值為0.16，標(biāo)準(zhǔn)差為0.03。

表3 群體生長曲線模型各指標(biāo)的相關(guān)估計(jì)值Table 3 Estimated values of parameters in growth curve model for population

圖1 生長曲線模型趨勢圖Fig.1 Pattern trend chart of growth curve model

表4 肉用中國西門塔爾牛個(gè)體Gomportz模型參數(shù)描述性統(tǒng)計(jì)Table 4 The descriptive statistics of Gompertz model parameters for Chinese Simmental beef cattle individual

2.2 主成分分析

圖2為基于前兩個(gè)主成分將個(gè)體聚類的群體結(jié)構(gòu)圖，從圖2中可以明顯看出，試驗(yàn)動(dòng)物群體被分在5個(gè) 區(qū)域，表明個(gè)體之間存在比較明顯的群體分層現(xiàn)象，絕大多數(shù)個(gè)體位于右下角區(qū)域，其他4個(gè)區(qū)域有少量個(gè)體分布，所以選取前兩個(gè)主成分作為協(xié)變量來消除群體分層對關(guān)聯(lián)分析的影響。

圖2 根據(jù)PCA繪制群體結(jié)構(gòu)圖Fig.2 Population structure identified by principal components analysis

2.3 全基因組關(guān)聯(lián)分析

圖3給出了利用線性混合模型對生長曲線參數(shù)性狀進(jìn)行GWAS結(jié)果的Q-Q圖，其中，橫坐標(biāo)表示期望P值的負(fù)對數(shù)，縱坐標(biāo)表示實(shí)際觀察P值的負(fù)對數(shù)?？梢钥闯?，成熟體重性狀(A)的Q-Q圖效果最好，大多數(shù)點(diǎn)都位于對角線上，表明觀察值與期望值的吻合度很好，末尾SNPs的顯著性也比較好；達(dá)到最大生長率的時(shí)間(b)的Q-Q圖末端與對角線的吻合度有部分偏差；成熟率(K)的Q-Q圖SNPs的顯著性略差。

圖3 性狀A(yù)、b、K進(jìn)行GWAS結(jié)果的Q-Q圖Fig.3 Q-Q plots of GWAS for A,b,K

圖4為利用線性混合模型對生長曲線參數(shù)性狀進(jìn)行GWAS結(jié)果的Manhattan圖。對于性狀成熟體重(A)，共檢測到了9個(gè)顯著的SNPs，主要分布在4、7、10、11、15和22號(hào)染色體上，臨近或坐落于PLIN3、BSN、KCNS等基因；其中P值最小的SNP位點(diǎn)是ARS-BFGL-NGS-14531，位于7號(hào)染色體，P值為9.55×10-7。對于達(dá)到最大生長率的時(shí)間(b)來說，共檢測到了49個(gè)顯著的SNPs，主要分布在1、3、5、9、12、14和23號(hào)染色體上，臨近或坐落于TMCO1、ANGPTL2、IGF-1等基因；其中P值最小的位點(diǎn)是BovineHD0900028514，位于9號(hào)染色體上，P值為4.43×10-8。對于成熟率(K)，共

圖4 性狀A(yù)、b、K進(jìn)行GWAS的Manhattan圖Fig.4 Manhattan plots of GWAS for A,b,K

檢測到了7個(gè)顯著的SNPs，集中位于22和25號(hào)染色體上，臨近或坐落于GRM7和SHISA9基因；P值最小的位點(diǎn)為BovineHD2200005378，位于22號(hào)染色體上，P值為3.24×10-6。具體信息如表5所示。

表5 性狀A(yù)、b、K的GWAS結(jié)果Table 5 The result of GWAS for A,b,K

(續(xù)表5 Continued)

3 討 論

3.1 生長曲線的擬合

表3給出了3種曲線模型的擬合結(jié)果，3種模型的擬合度R2都在0.95以上，模擬效果最好的是Gompertz模型，R2達(dá)到了0.954，而其他兩種模型R2雖然也有0.951，但是根據(jù)參數(shù)發(fā)現(xiàn)，兩種模型的成熟體重分別為551.0和1 458.5 kg，而中國肉用西門塔爾牛的成熟體重在600～800 kg[27]，兩種模型的成熟體重偏差較大，說明Logistic模型和Brody模型可能不適合本次肉用中國西門塔爾牛體重?cái)?shù)據(jù)的擬合。從圖1中也可以看出，擬合效果最好的是Gompertz模型，其曲線與原始數(shù)據(jù)曲線基本重合，這與梁永虎等[29]對于體重?cái)?shù)據(jù)所選的生長曲線模型相一致。表4給出了Gompertz模型參數(shù)的描述性統(tǒng)計(jì)，模型中參數(shù)A代表成熟體重，然而，A也出現(xiàn)了一些極值個(gè)體，如388.94 kg，對比原始記錄中18月齡體重僅為346 kg，分析原因可能是個(gè)體發(fā)育異常，生長速度較慢等；參數(shù)b代表達(dá)到最大生長率的時(shí)間，說明中國肉用西門塔爾牛在3月齡左右生長速度最快，最大值為9月齡，最小值在2月齡，而烏拉蓋地區(qū)的肉用西門塔爾牛在出生到1周歲這一時(shí)期內(nèi)生長速度最快，符合本研究結(jié)果，因此選擇Gompertz模型獲得的參數(shù)A、b、K作為表型來進(jìn)行GWAS。

3.2 GWAS和候選基因的功能分析

參數(shù)A、b、K的Q-Q圖整體效果都較好，大多數(shù)點(diǎn)都位于對角線上，表明觀察值和期望值的吻合度很高；Manhattan圖顯示出成熟體重(A)性狀檢測到了9個(gè)顯著的SNPs，多數(shù)位于7和11號(hào)染色體上，對于性狀達(dá)到最大生長率的時(shí)間(b)來說，共檢測到了49個(gè)顯著的SNPs，對于成熟率(K)來說，檢測到了7個(gè)顯著的SNPs，集中分布于22和25號(hào)染色體上。

通過基因功能注釋篩選到了與顯著SNPs關(guān)聯(lián)的基因，并且其中有部分基因已經(jīng)被鑒定出與生長發(fā)育性狀有關(guān)系。對于成熟體重(A)來說，ARS-BFGL-NGS-14531位點(diǎn)的顯著性最好，它臨近PLIN3基因，PLIN3基因是脂肪組織中脂肪分解和甘油三酯儲(chǔ)存的重要調(diào)節(jié)因子[31]，與皮脂腺脂肪生成相關(guān)的脂肪生成途徑交織在一起，如去飽和甘油三酯的合成[32]；3個(gè)顯著的SNPs都臨近KCNS3基因，而KCNS3基因被證明與脂肪百分比(%BF)性狀顯著關(guān)聯(lián)[33]，可作為發(fā)育性狀的候選基因；在性狀達(dá)到最大生長率的時(shí)間(b)中，有兩個(gè)顯著的SNPs坐落于TMCO1基因內(nèi)，其與PRKAG3基因有顯著的相互作用，這種作用可能會(huì)影響肌肉發(fā)育[34]；2個(gè)顯著SNPs臨近ANGPTL2基因，有研究表明，其可能作為新型的脂肪細(xì)胞因子[35]，與動(dòng)物的脂肪發(fā)育有顯著關(guān)聯(lián)；1個(gè)顯著SNP坐落于CFB基因內(nèi)，CFB基因被鑒定出與豬的出生仔豬數(shù)有關(guān)[36]；4個(gè)顯著的SNPs坐落于或臨近IGF-1基因，而IGF-1基因被認(rèn)為在生長發(fā)育中起著中心作用[37]；4個(gè)顯著的SNPs坐落于或臨近ALPL基因，有研究表明，肥胖者的白細(xì)胞中ALPL的基因表達(dá)水平明顯高于瘦者，表明ALPL基因可能與肥胖有關(guān)[38]；有1個(gè)顯著的SNP臨近ASPH基因，而ASPH基因參與調(diào)控肉牛胴體和肉質(zhì)性狀[39]；2個(gè) 顯著的SNPs坐落于EPHA4基因內(nèi)，其是豬生殖性狀的潛在候選基因之一[40]；成熟率(K)的7個(gè) 顯著SNPs集中分布于22和25號(hào)染色體上，并且發(fā)現(xiàn)了兩個(gè)關(guān)聯(lián)基因，分別為GRM7和SHISA9基因，其中SHISA9基因被發(fā)現(xiàn)與生長發(fā)育有關(guān)[41]。

4 結(jié) 論

本研究首先分別利用Gompertz、Logistic、Brody 3種 生長曲線模型對中國肉用西門塔爾牛的體重進(jìn)行生長曲線的擬合，選擇R2最高的Gompertz模型所獲得的參數(shù)A、b、K作為表型進(jìn)行GWAS，定位到了一些與生長發(fā)育性狀相關(guān)的候選基因，為其他縱向數(shù)據(jù)的研究提供了參考，也為調(diào)控肉牛生長發(fā)育進(jìn)而提高產(chǎn)肉量的育種工作提供了新的候選分子標(biāo)記。