中國肉用西門塔爾牛生長曲線參數的全基因組關聯分析

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

(1.西南大學動物科學技術學院，重慶 400715；2.中國農業科學院北京畜牧獸醫研究所，北京 100193)

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

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

1 材料與方法

1.1 試驗動物表型的描述性統計

試驗動物來自中國農業科學院北京畜牧獸醫研究所在內蒙古自治區錫林郭勒盟烏拉蓋管理區構建的中國肉用西門塔爾牛資源群體。本研究挑選了同時包含0、6、12、18月齡體重測量數據的個體，共808頭。表1為表型數據的描述性統計信息。

表1 中國肉用西門塔爾牛4個發育時期體重的描述性統計Table 1 The descriptive statistics of body weight at 4 growth stages of Chinese Simmental beef cattle kg

1.2 基因分型和質量控制

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

1.3 生長曲線的擬合

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

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

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

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

1.4 全基因組關聯分析

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

y=Xβ+y*

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

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

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

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

利用錯誤發現率(false discovery rate,FDR)多重檢驗方法對GWAS分析得到的P值進行檢驗，矯正P值的公式為：

P=FDR×n/m

其中，n代表小于FDR值的SNP個數，m代表SNP總個數，FDR設定值為0.05。

1.5 基因功能注釋

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

2 結 果

2.1 生長曲線的擬合

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

表3 群體生長曲線模型各指標的相關估計值Table 3 Estimated values of parameters in growth curve model for population

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

表4 肉用中國西門塔爾牛個體Gomportz模型參數描述性統計Table 4 The descriptive statistics of Gompertz model parameters for Chinese Simmental beef cattle individual

2.2 主成分分析

圖2為基于前兩個主成分將個體聚類的群體結構圖，從圖2中可以明顯看出，試驗動物群體被分在5個 區域，表明個體之間存在比較明顯的群體分層現象，絕大多數個體位于右下角區域，其他4個區域有少量個體分布，所以選取前兩個主成分作為協變量來消除群體分層對關聯分析的影響。

圖2 根據PCA繪制群體結構圖Fig.2 Population structure identified by principal components analysis

2.3 全基因組關聯分析

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

圖3 性狀A、b、K進行GWAS結果的Q-Q圖Fig.3 Q-Q plots of GWAS for A,b,K

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

圖4 性狀A、b、K進行GWAS的Manhattan圖Fig.4 Manhattan plots of GWAS for A,b,K

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

表5 性狀A、b、K的GWAS結果Table 5 The result of GWAS for A,b,K

(續表5 Continued)

3 討 論

3.1 生長曲線的擬合

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

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

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

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

4 結 論

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