董瑞蘭 孫國強 于光輝

(青島農(nóng)業(yè)大學(xué)動物科技學(xué)院，青島 266109)

畜牧養(yǎng)殖對環(huán)境造成的污染問題已引起全球的廣泛關(guān)注[1]。奶牛氮排泄量與飼糧的粗蛋白質(zhì)(CP)含量和消化率密切相關(guān)[2]。在奶業(yè)生產(chǎn)中，為了提高產(chǎn)奶量和乳品質(zhì)而提高飼糧CP的含量，導(dǎo)致70%～80%的攝入氮通過糞便和尿液排泄到水體、陸地和大氣中[3-4]，進(jìn)一步造成環(huán)境的惡化。我國農(nóng)業(yè)農(nóng)村部開展的全國農(nóng)業(yè)污染源普查中，畜禽養(yǎng)殖污染物的排放是基于動物的產(chǎn)排污系數(shù)核算的[5]。這種方法操作起來方便，但是糞污養(yǎng)分含量隨飼料的來源和組成不同而發(fā)生變化，用固定的系數(shù)估計變化的排放數(shù)值是不準(zhǔn)確的。因此，通過影響氮排泄量的指標(biāo)開發(fā)預(yù)測模型是準(zhǔn)確評估養(yǎng)分排放的有效方法[6]，有助于提高糞污養(yǎng)分預(yù)測的準(zhǔn)確性。而且，鑒于尿液中含有的氮比糞中的氮更易揮發(fā)[7]，通過排泄路徑(尿液和糞便)預(yù)測氮排放量更有意義。

關(guān)于奶牛的氮排泄量已有較多預(yù)測方程報道，有根據(jù)乳中尿素氮(milk urea nitrogen，MUN)[8-9]、氮攝入量(nitrogen intake，NI)[10-12]、瘤胃降解蛋白(rumen degradable protein，RDP)[13]預(yù)測尿氮(urinary nitrogen，UN)的一元線性方程；也有根據(jù)NI、干物質(zhì)采食量(dry matter intake，DMI)和產(chǎn)奶量(milk yield，MY)組合[14]，或MUN、CP和DMI組合[15]預(yù)測UN的多元線性方程。有根據(jù)NI[10-11]或DMI[12]預(yù)測糞氮(fecal nitrogen，F(xiàn)N)的一元線性方程；也有根據(jù)NI和DMI[14]或有機(jī)物采食量(organic matter intake，OMI)組合[16-17]預(yù)測FN的二元方程?？偟?total nitrogen，TN)可根據(jù)MY[13]或NI[11-12,18]單獨預(yù)測或根據(jù)NI和DMI組合[14]，體重(body weight，BW)、CP與MY組合[18]以及BW、CP與DMI組合[19]預(yù)測。Dong等[20]用荷斯坦奶牛數(shù)據(jù)庫評估了上述預(yù)測模型，發(fā)現(xiàn)預(yù)測的準(zhǔn)確性普遍偏低。因此，本研究旨在用荷斯坦奶牛的氮排泄數(shù)據(jù)庫，構(gòu)建準(zhǔn)確的UN和FN預(yù)測模型，并驗證準(zhǔn)確性，為奶牛生產(chǎn)提供依據(jù)。

1 材料與方法

1.1 建立數(shù)據(jù)庫

荷斯坦奶牛氮排泄數(shù)據(jù)庫(包括32個發(fā)表研究的99個觀察值)來自Dong等[20]。采用文獻(xiàn)發(fā)表的荷斯坦奶牛氮排泄數(shù)據(jù)構(gòu)建數(shù)據(jù)庫。在線檢索的數(shù)據(jù)庫包括中國知網(wǎng)(CNKI，http://www.cnki.net/)、Science Direct(http://www.sciencedirect.com/)和Journal of Dairy Science(http://www.journalofdairyscience.org/)。確立相關(guān)文獻(xiàn)研究的標(biāo)準(zhǔn)是：1)研究對象是荷斯坦奶牛；2)從已出版文獻(xiàn)的數(shù)據(jù)庫檢索出在中國進(jìn)行的研究，代表現(xiàn)代奶牛養(yǎng)殖的實際條件，特點是使用全混合日糧(TMR)飼養(yǎng)和采用自動擠奶設(shè)備；3)需提供荷斯坦奶牛氮排泄的輸出變量信息，包括UN、FN、MUN、TN和尿氮占總氮的比例(UN/TN)；4)需提供關(guān)于飼糧的CP含量、奶牛的NI、DMI、BW和MY及乳成分等特征輸入變量信息；5)若研究中某些變量未報道，但可通過其他已知變量計算的也包括在數(shù)據(jù)庫中。例如，飼糧的CP含量可通過DMI和NI計算獲得，TN為UN與FN的總和，UN/TN為UN與TN的比值。剔除觀察值變量缺失的5個研究[21-25]，剩余27個研究的89個觀察值對應(yīng)的所有變量均已知，即如果文獻(xiàn)的觀察值提供了關(guān)于NI、DMI、CP、BW、粗飼料比例(forage proportion，F(xiàn)P)、MY、乳氮(milk nitrogen，MN)、MUN、氮的全消化道表觀消化率(apparent total tract nitrogen digestibility，TTND)、乳氮占氮攝入量的比例(MN/NI)、氮攝入量占干物質(zhì)采食量的比例(NI/DMI)、飼料利用效率(feed efficiency，F(xiàn)E，MY與DMI的比值)、氮沉積量(nitrogen retention，NR)、UN、FN、TN、UN/TN等所有變量信息，則被包含在數(shù)據(jù)庫中。如果有觀察值的某個變量未知則不被包含在數(shù)據(jù)庫中。

數(shù)據(jù)庫經(jīng)JMP統(tǒng)計軟件馬氏離群分析(Mahalanobis outlier analysis，MOA)認(rèn)定為極值的數(shù)據(jù)，其學(xué)生化殘差在-2.5～2.5的范圍之外，被鑒定為離群異常值。通過分析，共計8個觀察值被鑒定為離群異常值，剔除離群異常值后，最后剩余27個研究的81個觀察值(隨機(jī)選擇2/3數(shù)據(jù)用于建模，1/3數(shù)據(jù)用于分析模型預(yù)測值與實際觀察值之間的回歸關(guān)系)進(jìn)行最終的分析。

1.2 構(gòu)建預(yù)測模型

1.2.1 篩選最佳預(yù)測變量

模型的最佳預(yù)測變量采用SAS 8.1軟件Corr過程程序的皮爾森相關(guān)系數(shù)(r)確定[1]。依次分析輸入變量(NI、DMI、CP、MUN、BW、TTND、NI/DMI、MY、MN、MN/NI、FP、NR、FE)與輸出變量(UN、FN、TN、UN/TN)之間的相關(guān)關(guān)系。

1.2.2 自變量和因變量之間的回歸關(guān)系

經(jīng)1.2.1中篩選出的最佳預(yù)測變量作為自變量，將奶牛的氮排泄量(UN、FN、TN、UN/TN)作為因變量，分析自變量和因變量之間的一元和多元回歸關(guān)系。

進(jìn)行多元回歸關(guān)系統(tǒng)計分析時，采用逐步回歸分析法(step wise selection，SWS)剔除不顯著(P>0.05)的變量，篩選出最終保留在模型中的自變量。多個變量之間的共線性用方差膨脹因子(variance inflation factors，VIF)檢驗，它是衡量每個自變量與模型中其他變量之間緊密相關(guān)的指標(biāo)。預(yù)測變量若具有較大的VIF(>2.5)，可導(dǎo)致模型發(fā)生較大的膨脹誤差。當(dāng)確定存在共線性時，預(yù)測氮排放量時皮爾森相關(guān)系數(shù)最小的變量即被剔除。

1.2.2.1 一元線性模型

自變量和因變量之間的一元回歸關(guān)系用以下方程分析：

y=bx+a。

式中：x指因變量；y指UN(g/d)、FN(g/d)、TN(g/d)和UN/TN；a是固定常數(shù)；b指對應(yīng)變量的系數(shù)。

1.2.2.2 多元線性模型

自變量和因變量之間的多元回歸關(guān)系用以下方程分析：

y=b1(x1)+b2(x2)+b3(x3)+…+bn(xn)+a。

式中：x指因變量；y指UN(g/d)、FN(g/d)、TN(g/d)和UN/TN；a是固定常數(shù)；b1，b2，b3，…，bn指各對應(yīng)變量的系數(shù)。

1.2.3 模型的擬合效果

用SAS 8.1軟件分析均方根誤差(root mean square error，RMSE)和阿凱克信息論準(zhǔn)則(Akaike’s information criterion，AIC)值。RMSE和AIC值是模型擬合好壞的指標(biāo)，其值越小代表模型擬合效果越好[1]。

1.3 模型預(yù)測值與實際觀察值之間的回歸關(guān)系分析

將模型輸出的預(yù)測值作為自變量，實際觀察值作為因變量，繪圖分析模型預(yù)測值與實際觀察值之間線性回歸關(guān)系的顯著性，獲得決定系數(shù)(R2)，同時確定回歸線和等線(y=x)的偏離程度。

1.4 數(shù)據(jù)的統(tǒng)計與分析

用SAS 8.1軟件的混合回歸模型(mixed regression model，MRM)過程程序分析自變量和因變量之間的一元和多元回歸關(guān)系，P<0.05為回歸關(guān)系顯著，P<0.01為回歸關(guān)系極顯著。通常認(rèn)為，模型預(yù)測值與實際觀察值之間R2越大，則獲得的預(yù)測模型越準(zhǔn)確。在用SAS 8.1軟件程序分析時，不同的研究作為隨機(jī)變量也包括在模型中，用于消除不同研究之間的試驗誤差。模型預(yù)測值和實際觀察值之間的線性關(guān)系用以下方程分析：

y=bx+a。

式中：x指UN(g/d)、FN(g/d)、TN(g/d)和UN/TN的模型預(yù)測值；y指UN(g/d)、FN(g/d)、TN(g/d)和UN/TN的實際觀察值；a是固定常數(shù)；b指對應(yīng)變量的系數(shù)。

2 結(jié)果與分析

2.1 建模數(shù)據(jù)庫

由表1可知，奶牛的BW、CP、DMI、TTND、NI平均值分別為565 kg、15.4%、19.20 kg/d、69.8%和470 g/d，變動范圍分別為407～755 kg、11.0%～19.2%、11.57～28.10 kg/d、58.9%～78.0%和214～720 g/d。奶牛的UN和FN的變動范圍較大，介于58.2～329.0 g/d、81.0～193.0 g/d，平均值分別為168.0和141.0 g/d。奶牛的TN變動范圍為139.2～481.3 g/d。UN/TN平均值為0.54，變動范圍為0.39～0.69。

表1 基于動物和飼糧因子以及氮輸入和輸出變量的建模數(shù)據(jù)庫描述性統(tǒng)計

2.2 奶牛氮排泄量與動物和飼糧因子之間的皮爾森相關(guān)系數(shù)

由表2可知，UN分別與NI、DMI、CP、MUN、TTND、NI/DMI、MY和MN呈顯著或極顯著正相關(guān)關(guān)系(P<0.05或P<0.01)。FN分別與NI、DMI、CP、NI/DMI、MY、MN和FE呈極顯著正相關(guān)關(guān)系(P<0.01)，而與TTND和FP呈顯著或極顯著負(fù)相關(guān)關(guān)系(P<0.05或P<0.01)。TN分別與NI、DMI、CP、MUN、NI/DMI、MY和MN呈極顯著正相關(guān)關(guān)系(P<0.01)。UN/TN分別與DMI、MUN、TTND和FP呈顯著或極顯著正相關(guān)關(guān)系(P<0.05或P<0.01)，而與NR呈極顯著負(fù)相關(guān)關(guān)系(P<0.01)。

表2 奶牛氮排泄量與動物和飼糧因子之間的皮爾森相關(guān)系數(shù)

2.3 一元線性預(yù)測模型

2.3.1 UN預(yù)測模型

預(yù)測UN的模型見表3(模型1～6)。輸入變量NI、CP、NI/DMI、MY、MN/NI、MN對回歸關(guān)系都有顯著影響(P<0.05)。根據(jù)RMSE和AIC值可知，基于單一變量的UN模型擬合效果為：MN/NI

表3 奶牛尿氮、糞氮和總氮的一元線性預(yù)測模型

2.3.2 FN預(yù)測模型

預(yù)測FN的模型見表3(模型7～12)。輸入變量NI、MY、MN、DMI、NI/DMI、CP對回歸關(guān)系都有顯著影響(P<0.05)。當(dāng)單獨用NI(模型7)預(yù)測FN時，其RMSE和AIC值最低。但單獨用CP(模型12)預(yù)測FN時，其RMSE和AIC值最高。用NI預(yù)測的FN與實際觀察值之間的比較見圖1-b，回歸分析顯示，基于NI的回歸模型預(yù)測值和實際觀察值之間的R2為0.84(P<0.01)。

2.3.3 TN預(yù)測模型

預(yù)測TN的模型見表3(模型13～18)。輸入變量NI、MY、MN、NI/DMI、CP、FP對回歸關(guān)系都有顯著影響(P<0.05)。當(dāng)單獨用NI(模型13)預(yù)測TN時，其RMSE和AIC值最低。但單獨用FP(模型18)預(yù)測FN時，其RMSE和AIC值最高。用NI預(yù)測的TN與實際觀察值之間的比較見圖1-c，回歸分析顯示，基于NI的回歸模型預(yù)測值和實際觀察值之間的R2為0.93(P<0.01)。

圖1 尿氮(a)、糞氮(b)、總氮(c)和尿氮占總氮的比例(d)的實際觀察值和模型預(yù)測值之間的關(guān)系

2.3.4 UN/TN預(yù)測模型

UN/TN預(yù)測模型見表3(模型19～22)。輸入變量TTND、MN/NI、CP、FE對回歸關(guān)系都有顯著影響(P<0.05)。當(dāng)單獨用TTND(模型19)預(yù)測UN/TN時，其RMSE和AIC值最低。單獨用TTND預(yù)測的UN/TN與實際觀察值之間的比較見圖1-d，回歸分析顯示，基于TTND的回歸模型預(yù)測值和實際觀察值之間的R2為0.90(P<0.01)。

2.4 多元線性模型

2.4.1 UN預(yù)測模型

預(yù)測UN的多元線性回歸模型見表4(模型23和24)。在多個輸入變量中，CP、DMI與NI/DMI組合以及DMI與NI/DMI組合對UN多元線性模型有極顯著影響(P<0.01)。模型23的RMSE和AIC值與模型24相比較低。用模型23預(yù)測的UN與實際觀察值之間的比較見圖2-a，回歸分析表明，模型23的預(yù)測值和實際觀察值之間的R2為0.90(P<0.01)。

表4 奶牛尿氮、糞氮和總氮的多元線性預(yù)測模型

2.4.2 FN預(yù)測模型

預(yù)測FN的多元線性回歸模型見表4(模型25和26)，在多個輸入變量中，MY與NI組合以及MY與DMI、FP組合對FN多元線性模型有極顯著影響(P<0.01)。模型25的RMSE和AIC值較模型26低。用模型25預(yù)測的FN與其實際觀察值之間的比較見圖2-b，回歸分析表明，模型25的預(yù)測值和實際觀察值之間的R2為0.87(P<0.01)。

圖2 尿氮(a)、糞氮(b)、總氮(c)和尿氮占總氮的比例(d)的實際觀察值和模型預(yù)測值之間的關(guān)系

2.4.3 TN預(yù)測模型

預(yù)測TN的多元線性回歸模型見表4(模型27～29)，在多個輸入變量中，NI與NI/DMI組合、DMI與NI/DMI組合以及DMI與MY組合對TN多元線性模型有極顯著影響(P<0.01)。模型27的AIC和RMSE值是3個TN預(yù)測模型中最低的。模型29的AIC和RMSE值是3個TN預(yù)測模型中最高的。用模型27預(yù)測的TN與其實際觀察值之間的比較見圖2-c，回歸分析顯示，模型27的預(yù)測值和實際觀察值之間的R2為0.92(P<0.01)。

2.4.4 UN/TN預(yù)測模型

預(yù)測UN/TN的多元線性回歸模型見表4(模型30和31)。NI、TTND、FP與NR組合以及NI、TTND與NR組合對UN/TN多元線性模型有極顯著影響(P<0.01)。模型31的RMSE和AIC值最低。用模型31預(yù)測的UN/TN與實際觀察值之間的比較見圖2-d，回歸分析顯示，模型31的預(yù)測值和實際觀察值之間的R2為0.89(P<0.01)。

3 討 論

3.1 建模數(shù)據(jù)庫

數(shù)據(jù)庫涵蓋了奶牛泌乳初期、中期和后期各產(chǎn)乳階段。但需指出的是，由于研究使用了發(fā)表文獻(xiàn)數(shù)據(jù)，收集的數(shù)據(jù)是各試驗處理的觀察值平均值，而不是動物個體的試驗數(shù)據(jù)，所以本研究未考慮動物之間的個體差異。與美國1990—1995年泌乳奶牛數(shù)據(jù)庫[12]DMI范圍10.3～28.2 kg/d相比，荷斯坦奶牛的最低和最高DMI數(shù)值與之接近。將荷斯坦奶牛的氮排泄數(shù)據(jù)庫與Spek等[15]建立的北歐和北美數(shù)據(jù)庫相比，中國NI的平均值低于北歐(470 g/d vs. 485 g/d)和北美(470 g/d vs. 637 g/d)。同樣，UN的平均值也低于北歐(168 g/d vs. 185 g/d)和北美(168 g/d vs. 212 g/d)；FN的平均值也低于北歐(141 g/d vs. 159 g/d)和北美(141 g/d vs. 223 g/d)。

3.2 UN預(yù)測模型

用單一變量預(yù)測UN時，基于NI(模型1)預(yù)測的準(zhǔn)確性(AIC值=510.2，R2=0.91)最高。關(guān)于NI和UN之間存在顯著的相關(guān)關(guān)系已有較多報道[10-12]。Castillo等[10]和Kebreab等[11]研究發(fā)現(xiàn)，NI和UN之間呈指數(shù)關(guān)系，原因可能為在氮平衡為正的情況下，UN不是直接測定，而是用物料守恒法計算差值獲得，由此導(dǎo)致了UN的偏高估計。荷斯坦奶牛排泄的UN與NI之間呈現(xiàn)正相關(guān)線性關(guān)系，而且用NI可準(zhǔn)確預(yù)測奶牛的UN，該結(jié)果與Kebreab等[11]和Reed等[12]的結(jié)論一致。本研究中UN和NI之間的回歸關(guān)系斜率與Reed等[12](0.34 vs. 0.33)和Kebreab等[11](0.34 vs. 0.38)報道的范圍接近，而且本研究的回歸關(guān)系常數(shù)項值更低(4.8 vs. 12.0和4.8 vs. 20.0)。

當(dāng)用多個變量預(yù)測UN時，模型23提供了較好的預(yù)測方法(R2=0.90)。Spek等[15]用DMI和CP對歐洲和北美奶牛的UN進(jìn)行預(yù)測，發(fā)現(xiàn)在輸入變量DMI和CP的基礎(chǔ)上增加第3個變量MUN，比單獨用MUN或CP以及同時用MUN和CP的組合預(yù)測得更準(zhǔn)確。本研究中，模型23的回歸關(guān)系常數(shù)項與Spek等[15]報道的一致，而且AIC值更低(488.3 vs. 988.0)。

3.3 FN預(yù)測模型

用單一變量NI、MY、MN、DMI、NI/DMI、CP預(yù)測FN時，模型7的擬合效果和準(zhǔn)確性最高。Castillo等[10]基于英國的奶牛數(shù)據(jù)發(fā)現(xiàn)，F(xiàn)N與NI存在線性函數(shù)關(guān)系(FN=0.21NI+52.3)。此方程的斜率與本研究的模型7一致，但是Castillo等[10]建立的模型的常數(shù)項值約為本研究中模型7的1.2倍。Vasconcelos等[26]也報道了奶牛FN與NI之間存在線性關(guān)系。Kebreab等[11]報道了奶牛FN與NI之間也存在類似的線性關(guān)系，但FN的斜率為0.28，截距數(shù)值為10。Schuba等[27]研究發(fā)現(xiàn)NI和FN之間呈二次關(guān)系，即FN=117.3+3.24×NI2-0.004 3×NI，但其AIC值高達(dá)1 509，說明FN和NI之間二次關(guān)系擬合效果低。

用多個變量預(yù)測荷斯坦奶牛排泄的FN時，基于NI和MY的FN預(yù)測模型(模型25)擬合性略高，預(yù)測效果較好(R2=0.87)。將多元線性模型(模型25)與一元線性模型(模型7)相比，發(fā)現(xiàn)多元線性模型的預(yù)測值和實際觀察值之間的R2較高，而常數(shù)項值也低于一元線性模型。

3.4 TN預(yù)測模型

將NI、CP、MY、MN、FP與NI/DMI作為單一輸入變量時，基于NI(模型13)預(yù)測TN的準(zhǔn)確性最高(AIC值=505.6，R2=0.93)。用多個變量預(yù)測TN時，模型27預(yù)測的準(zhǔn)確性較高。將基于NI的TN預(yù)測模型(模型13)分別與基于NI的UN預(yù)測模型(模型1)和基于NI的FN預(yù)測模型(模型7)相比，模型13的擬合效果優(yōu)于模型1和模型7。

3.5 UN/TN預(yù)測模型

考慮到UN在TN中占較大比例，建立了UN/TN預(yù)測模型。與MN/NI(AIC值=-157.2，模型20)、CP(AIC值=-139.8，模型21)或FE(AIC值=-148.3，模型22)相比，將TTND(AIC值=-158.0，模型19)作為單一輸入變量的模型預(yù)測效果最好。通過UN/TN與動物或飼糧因子的多變量回歸關(guān)系及逐步回歸分析發(fā)現(xiàn)，變量TTND與UN/TN的相關(guān)系數(shù)最高。因此，TTND被證明是預(yù)測奶牛UN/TN的最佳預(yù)測變量，該發(fā)現(xiàn)與在肉牛中[1]的研究結(jié)果一致。Huhtanen等[14]將CP作為奶牛UN/TN的預(yù)測變量，發(fā)現(xiàn)二者存在顯著的線性關(guān)系。但本研究結(jié)果發(fā)現(xiàn)用TTND預(yù)測UN/TN比CP效果更好，原因是TTND與CP密切相關(guān)。飼糧CP和氮的消化率越高，就會導(dǎo)致越多的氮以瘤胃氨或小腸氨基酸的形式被吸收。那么，超過奶牛營養(yǎng)需要的剩余氮量就會進(jìn)入尿液，從而增加了UN的比例。需要注意的是，本研究僅基于發(fā)表文獻(xiàn)進(jìn)行建模分析，最好通過試驗獲得數(shù)據(jù)進(jìn)一步驗證其預(yù)測性。

4 結(jié) 論

① 基于NI的一元UN、FN和TN回歸模型以及基于TNND的一元UN/TN回歸模型的RMSE和AIC值最低，一元回歸模型擬合效果好。

② 基于DMI、CP與NI/DMI組合的多元UN回歸模型，基于NI與MY組合的多元FN回歸模型，基于NI與NI/DMI組合的多元TN回歸模型的擬合效果好。