Introduction
지구의 평균기온은 지난 133년간 0.85°C (0.65~1.06°C) 상승하였으며, 최근 지구 온난화가 지속 되고 있다. 세계 온실가스 배출량의 지속적인 증가로 기온 상승은 가속화될 예정이다(기상청, 2015). 또한 국내 기후도 온난화가 진행되고 있으며, 지난 30년(1981~2010년)간 연평균 기온은 1.2°C, 연평균 폭염일수(일 최고기온이 33°C 이상인 날의 연중일수)는 지역별로 차이가 있지만 수도권은 6.6일, 대구 23.2일로 증가하였고, 21세기 후반 연평균 기온은 1981~2010년보다 5.7°C 상승할 전망이다(환경부, 2015). 국내에서 사육 비중이 높은 젖소인 홀스타인 품종은 고온에 약하여 주변온도가 임계점(27°C)을 초과하였을 시 고온스트레스를 받는다고 보고되어 있다(Armstrong, 1994; Nguyen et al., 2016). 고온스트레스 환경에서 발생되는 문제점은 사료섭취량, 산유량, 유질 및 수태율이 감소하여 낙농 생산비 증가와 수입 감소를 초래한다. 기후변화에 따라 낙농산업에서 고온스트레스의 문제점은 점점 커지고 있다(Beede and Collier, 1986; Hansen, 2007; Hammami et al., 2013; Nguyen et al., 2016). 온도와 습도를 이용한 온습도지수(Temperature-Humidity index, THI)는 열부하(heat load)의 지표형질로 이용하면 고온스트레스 상태의 진단이 가능하다. 또한 THI는 검정 기록과 함께 이용하면 산유능력의 변화를 고온스트레스의 판단지표로 할 때, 고온 저항성을 증진을 위한 선발이 가능하다고 보고되었다(Ravagnolo and Misztal, 2000; Nguyen et al., 2016). 지구온난화의 진행으로 젖소의 고온스트레스 환경 기후변화에 맞는 육종계획 수립이 필요하다.
본 연구는 고온 스트레스 저항성증진을 위한 평가모형 선정의 기초자료로 활용하기 위해 국내의 기후변화 추세를 파악하고, 비유곡선과 불완전 Gamma모형(Incomplete gamma distribution model) 을 이용하여 추정된 비유곡선 모수 값들과 THI와의 상관관계 분석을 실시하였다.
Materials and Methods
기상자료
기상 자료는 기압, 기온, 풍속, 상대습도 등 14개 요소를 측정할 수 있는 종관기상관측 장비 (Automated Surface Observing System, ASOS)가 설치 된 79개 지역에서 2000년 1월 1일부터 2016년 12월 31일까지 측정된 일 최고온도, 일 평균온도, 일 최저온도, 일 평균상대습도자료 (https://data.kma.go.kr)를 이용하였다.
온습도지수 THI는 NRC(1971)가 개발한 가축의 온습도지수 예측모형을 이용하여 다음과 같이 계산하였다.
THI=(1.8×T+32)-[0.55-0.0055×RH)×(1.8×T-26.8)
여기서, THI는 온습도지수, T는 해당 지역의 검정일 최고온도(°C), RH는 해당지역의 검정일 평균상대습도(%)이다.
THI 계산 시 일 최고기온과 일평균상대습도를 이용하여 THI를 계산하고 이를 일 최고 THI라고 정의하였다.
젖소 생산 형질자료
농협 경제지주 젖소 개량사업소에서 진행하는 유우군 능력검정사업에 참여하고 있는 홀스타인 암소 중 2000년 1월 1일부터 2016년 12월 31일까지 분만기록을 가지며 월 1회 검정일 검정기록이 3회 이상인 초산 암소 11,236두의 이상치를 제거한 116,611건의 검정일 기록을 분석에 이용하였다. 분석에 이용된 자료의 기초통계량은 Table 1에 제시하였다. 1산차 검정일 기록을 가지는 개체들에 대해 개체별로 비유곡선(불완전 Gamma 모형)을 추정하고 그 중 표준곡선모양을 가지며 기상자료를 이용 가능한 4,809두에 대하여 비유곡선 모수들과 고온스트레스 요인들과의 상관관계를 분석하였다.
비유곡선 추정
비유곡선의 모수를 추정하기 위하여 Wood(Wood, 1967)의 불완전 gamma 모형을 이용하였으며, 모형은 다음과 같다.
여기서, 
는 분만후 t번째 착유일의 산유량, t는 분만 후 일수, A는 분만직후 첫 번째 검정일 유량, b는 비유곡선의 상향경사도, c는 비유곡선의 하향경사도를 표시하는 상수이며, e는 자연대수이다.
그리고 최고유량 도달시기(
), 최고 추정유량(
) 및 비유지속성(s)은 각각 , 
, 
및 
로 구할 수 있으며, Wood모형은 비유곡선의 모수는 log 선형식 
의 미분값으로 추정이 용이하다 (Rekaya et al., 2000; Won et al., 2014; Won et al,. 2016).
비유곡선함수의 추정은 SAS 9.2의 NLIN Procedure (Nonlinear procedure) 를 사용하였고, 편도함수의 지정이 필요하지 않은 DUD 방법(Doesn’t Use Derivative)을 이용하였다.
고온스트레스 정도
고온스트레스와 비유곡선의 상관관계를 알아보기 위하여 고온스트레스 노출일수(HSday)는 최고 유량 도달 기간 중 THI가 72이상인 일 수로 정의하였고, 고온스트레스 노출비율은 다음과 같이 계산하였다.
HSraito = HSday÷
여기서, HSraito는 고온스트레스 노출비율, HSday는 고온스트레스 노출일수, 
는 최고유량 도달시기이다.
Results and Discussion
국내 기후변화 추세
Figure 1에 본 연구에서 이용한 기상자료로부터 지난 15년간 국내의 기후변화 추세를 제시하였다. 환경부(2015)는 지난 과거 1981~2010년 까지 연평균 기온 1.2°C(0.41°C/10년) 상승하였다고 보고하였다. 2000년도와 2016년의 온도차이는 최고 온도는 0.9°C, 평균 온도는 1.1°C, 최저 온도는 1.2°C 증가하였으며, 온습도지수(Temperature-Humidity index, THI)도 1.5unit 상승하여 기온과 유사한 상승 추세를 보였다.
Figure 2에는 홀스타인 암소의 임계점 72(THI) 이상인 일수를 제시하였다. 홀스타인 암소는 THI가 72 이상이면 생산능력이 감소한다고 보고되어 있으며(NRC, 1971; St-Pierre et al., 2003; Akyuz et al., 2010), 국내에서는 365일 중 최소 120일 이상이 최고 THI가 72 이상이였으며, THI가 72이상인 일수가 매년 0.88일 증가하는 추세를 보였다. 또한 국내 기후는 온대기후에서 아열대 기후로 변화하고 있는 중이라고 보고되었다(IPCC, 2012). 따라서 국내의 기후는 홀스타인 암소의 생산성에 악영향을 미치는 기후로 점차 변하고 있다고 사료된다.
비유곡선
Figure 3에는 THI 그룹에 따른 단순 비유곡선을 나타냈다. 산유량에 대한 사전 분석결과 THI 그룹, 분만계절, 분만월령은 모두 유의적인 영향을 미치는 것으로 나타났다 (p<0.001). THI 그룹 중 에서는 스트레스를 가장 크게 받을 것으로 예상한 집단 (THI > 89)에서 일착유량이 가장 낮게 나타났다. 지역별 기상자료는 젖소 농가의 직접 측정한 기상 자료는 아니다. 더욱더 정확한 평가를 위해서는 젖소 농가별 기상 자료를 사용해야겠지만 농가별 기상자료 수집에는 한계가 있다. 농가의 위치와 가장 인접한 관측소의 기상자료를 적용하기 위하여 농가가 속한 시군 관측소의 자료를 찾고, 해당 시군 관측소가 없는 경우에는 해발 고도와 산과 바다의 위치를 감안하여 지리적으로 가장 인접한 인근 시군의 관측자료를 이용하였다.
Figure 4에는 분만계절에 따른 평균 비유곡선을 나타냈다. 분만계절에 따른 착유량을 보면 여름(4~8월)에 분만한 암소들의 일 착유량이 다른 계절에 분만한 암소들보다 낮게 나타났다. 이는 여름철의 높은 기온과 사료섭취량의 감소에 의한 것으로 타 연구 결과와 유사하였다(Tekerli et al., 2000; 원정일 2014).
Figure 5에는 분만월령에 따른 평균 비유곡선을 제시하였다. 분만시 월령이 낮은 어린 개체(18~24 Month) 들의 일착유량이 24개월령 이상에서 분만한 개체들에 비해 모두 낮게 추정이 되었다.
비유곡선 모수 추정
Table 2 에는 Wood 모형으로 추정한 비유곡선 모수 값 A, b, c과 이들 곡선 형태 모수로부터 계산한 비유지속정도(s), 최고유량 도달일 및 최고 유량 추정량에 대한 기초 통계량을 제시하였다. 비유곡선을 추정한 모든 개체의 추정 모수 중에서 정규분포를 벗어난 이상치를 제거한 개체 중 기상자료와 결합이 가능한 개체들의 비유곡선 만을 분석에 이용하였다. 분만직후 유량을 나타내는 모수 값(A)은 18.07 ± 5.40, 비유곡선에서 상향경사도를 나타내는 모수 값(b)는 0.175(±0.081) 및 하향경사도를 나타내는 모수 값(c)는 0.00218 (±0.00113)으로 나타났다. 최고 추정유량은 31.95kg/일, 최고 유량 도달 시기는 83.69일 및 비유지속성(s)는 7.32으로 추정 되었다.
Table 3에는 고온스트레스 노출일수(Heatday(day)), 고온스트레스 노출 비율(Heatday(%)), 최고 유량 도달할 때까지의 가장 높은 THI(Max_THI)와 최고 유량 도달 할 때까지의 THI의 평균(Mean_THI)의 기초통계량을 제시하였다.
최고 유량 도달까지 평균 적으로 고온스트레스에 노출 일수 평균 30.83일, 노출 비율은 36.12% 그리고 최고유량 도달까지의 평균 THI는 63.45로 고온스트레스의 기준점보다는 낮았다. 겨울에 분만을 한 개체는 최고유량 도달 시기까지 고온스트레스 노출일수가 0이다. 고온스트레스에 노출되는 분만계절이 여름인 개체로만 추정하여 분석을 하면 고온스트레스에 노출 일수, 노출 비율, 최고유량 도살시기까지의 평균 THI가 상승할 것이다.
Table 4에는 고온스트레스 정도와 wood모형에서 추정된 최고 추정유량, 최고 유량 도달시기 및 비유지속성과의 상관관계를 제시하였다. 최고 추정유량은 고온스트레스의 지표들과 상관관계가 없었지만 최고유량 도달시기에서는 고온스트레스 노출일수와 최고 유량 도달할 때까지의 일 최고 THI와는 양의 상관관계를 보였다. 따라서 고온스트레스 노출 일수가 길어지거나 최고유량 도달시기 중 한번이라도 고온스트레스에 노출이 되면 최고유량 도달 시기는 길어지는 것으로 사료된다. 비유지속성에서도 최고유량도달시기와 유사한 추세를 보였지만 최고유량도달시기와 비유지속성의 상관관계는 0.84 (p<0.001) 로 매우 높았다. 그러나 비유지속성과 총 유량과의 관계는 알 수가 없어 비유지속성증가가 생산성 증대와 관련이 있다고 판단할 수는 없다.
Conclusions
지구 온난화의 영향으로 국내의 기후는 해마다 온도와 THI는 증가하고 있으며, 고온스트레스 노출 일수 또한 증가 추세였다. 따라서 젖소에 악영향을 미치는 기후로 점점 변하고 있어 국내에서는 고온스트레스를 대비를 해야 한다.
지역별 기상자료는 젖소 농가의 직접적인 기상 자료는 아니다. 더욱더 정확한 평가를 위해서는 젖소 농가 별 기상 자료를 사용해야 한다. 하지만 농가 별 기상자료 수집에는 한계가 있다. 따라서 지역별 기상자료를 사용과 동시에 농장의 효과를 고려하면 THI는 고온스트레스 저항성 관한 지표로 이용 가능 할 것으로 사료된다.
분만계절 별 비유곡선은 차이가 있었으며, 분만월령 별 비유곡선도 차이가 있었으며, 고온스트레스 노출 일수가 길어지거나 최고유량 도달시기 중 한번이라도 고온스트레스에 노출이 되면 최고유량 도달 시기는 길어졌다. 따라서 추후 고온저항성증진을 위한 유전능력 평가 모형에는 분만월령, 분만계절 및 착유일(Days in milk)을 고려해야 할 것이다.
또한 본 연구에서 고온 스트레스 정도와 총 생산량에 대한 추가적인 연구가 필요하다고 생각된다.
Summary
본 연구는 고온 스트레스 저항성증진을 위한 평가모형 선정의 기초자료로 활용하기 위해 실시하였다. 기상자료는 종관기상관측 장비(ASOS)가 설치 된 79개 지역에서 2000년 1월 1일부터 2016년 12월 31일까지 측정된 일 최고온도, 일 평균온도, 일 최저온도, 일 평균상대습도자료를 이용하였으며, 젖소의 검정자료는 2000년 1월 1일부터 2016년 12월 31일까지 수집된 암소 11,236두의 이상치를 제거한 116,611건의 검정 기록을 분석에 이용하였다. 2000년도 이후로 지속적으로 최고 온도, 평균 온도, 최저 온도, 온습도지수(Temperature-Humidity index, THI) 및 고온스트레스 노출일수는 증가하는 추세를 보였다. 분만계절, 분만 월령 및 THI 그룹별 비유곡선은 차이가 있었다. 고온스트레스 노출 일수가 길어지거나 최고유량 도달시기 중 한번이라도 고온스트레스에 노출이 되면 불완전 감마모형으로 추정된 최고유량 도달 시기는 늦어졌다. 따라서 추후 고온저항성증진을 위한 유전능력 평가 모형에는 분만월령, 분만계절 및 착유일 (Days in milk)을 고려해야 하며, 고온 스트레스 정도와 총 생산량에 대한 추가적인 연구가 필요하다고 생각된다.
