Introduction
한우는 국가단위 검정체계를 기반으로 씨수소를 선발하여 집단을 개량, 유지하고있는 한국의 고유품종이다. 한우의 씨수소 선발은 당대검정에서 후보씨수소 선발을 거쳐서 후대검정 종료 시 자손의 도체성적을 반영한 유전능력평가결과를 바탕으로 선발하고 있다. 이와 같은 유전능력평가의 핵심요인은 개량기관 간 자료공유와 협업을 통해 구축된 국가단위 유전능력평가체계와 한우 개량이 시작된 80년대 이래 꾸준히 쌓여온 검정자료를 포함한 기준집단(reference population)이 존재하기 때문이다. 기준집단은 표현형자료, 혈통자료와 유전체자료로 구성되며 한우 기준집단의 경우, 매년 약 1,300두의 검정자료, 약 2,000두의 유전체자료가 축적되고 있다. 그 중에서도 유전체자료는 유전체 선발이 핵심적인 평가방법이 된 현재로써 기존에 사용해오던 혈통을 대신하여 유전관계행렬(genetic relationship matrix, GRM)을 제공함으로써 평가결과의 정확도에 큰 비중을 차지하고 있다(Jang et al., 2019).
유전체자료는 세부적으로 Beadchip의 제조사에 따라 그 구성이 달라질 수 있고 제조사에서 배포하는 플랫폼에 따라서도 차이가 발생한다. 따라서 고유품종의 구분이나 연구를 목적으로 제조사의 상용 Beadchip을 사용하지 않고 해당 집단에서 발견될 수 있는 변이를 포함한 Custom Beadchip을 사용할 수 있으며 한우 또한 국가단위 유전능력평가를 위해 Custom Beadchip을 사용하여 유전체자료를 확보하고 있다. 한우 기준집단의 유전체자료 platform은 Illumina bovine v2, Illumina bovine v3 그리고 Hanwoo custom v1 세 가지가 있고 평가를 위해 기존에 개체수가 가장 많았던 Illumina bovine v2 칩을 기준으로 imputation 후 평가에 사용해오고 있다(Teo et al., 2007; Melka et al., 2011). 하지만 2019년도를 기준으로 Hanwoo custom v1으로 구성된 개체수가 더 많아지면서 현재는 Beadchip platform 별로 imputation된 유전체자료가 유전능력평가결과에 어떤 영향을 미치는지 검증이 필요한 시점에 도달했다.
현재 국가단위 유전능력평가에서 사용되는 조화행렬 기반의 ssGBLUP (single-step genomic best linear unbiased prediction)은 Henderson’s mixed model equation을 기반으로 하는 유전능력평가 프로그램에서 유전체기반 관계행렬에 95%의 가중치를, 혈통기반 개체행렬에 5%의 가중치를 두고 조화행렬로 혼합하여 사용하고 있다(Legarra et al., 2014; Fragomeni et al., 2017). 그렇기에 유전체자료의 비중이 더 클 수밖에 없는 실정이며 유전관계행렬의 구성 방법에 따라 결과가 다소 달라질 여지가 있다. 유전관계행렬의 구성방법으로는 VanRaden과 Yang의 방법(VanRaden, 2008; Yang et al., 2010)이 대표적으로 알려져 있으나 편향된 관계행렬 값을 조정하기 위해 사용되는 보정계수(scaling factor, SF)는 일반적으로 0 (Falconer and Mackay, 1996; Momin et al., 2021)을 기본값으로 취하여 사용하고 있다. 이와 같은 보정계수는 집단의 연관불평형에 따라 다르게 조정될 필요가 있고 그 차이에 대한 결과가 조명되어야 할 필요가 있다. 따라서 본 연구에서는 한우 국가단위 유전능력평가에 있어서 유전체자료의 칩 플랫폼과 유전관계행렬을 구성할 때, 그리고 집단의 편향을 조정하는 보정계수를 달리하였을 때 그 결과를 비교하기 위해 수행되었다.
Materials and Methods
유전능력평가용 기준집단 및 평가모형
Beadchip platform과 유전관계행렬 보정계수에 대한 육종가 추정 및 정확도 비교를 위해 기준집단에는 동일한 표현형 21,881두와 혈통 88,202두의 기록이 사용되었다. 유전체자료 또한 활용된 두 수는 15,777두로 동일하나, imputation을 위해 적용된 chip platform에 따라 단일염기다형성(single nucleotide polymorphism) 마커의 구성 및 밀도 그리고 공통된 대립유전자에 대해 imputation을 거쳐 빈도에 미세한 차이가 있다. 마커의 개수는 Hanwoo custom v1, Illumina bovine v2, Illumina bovine v3 각각 52,007, 52,791개 그리고 51,173개가 있고 육종가 추정 등을 위해 1종오류를 유발할 가능성이 있는 5% 수준 미만의 희소대립유전자(minor allele frequency)를 제외시킨 42,509, 39,745개 그리고 39,184개의 마커를 추출하고 각 플랫폼을 대표하는 마커로 활용하였다(Table 1).
Table 1. Number of markers by platform and after quality control.
|
|
SNP, single nucleotide polymorphism. |
평가모형은 현재 국가단위 평가체계에서 사용되는 성장형질인 12개월령 보정체중(adjusted 12 months weight, WT12)과 도체형질인 도체중(carcass weight, CWT), 등심단면적(eyemuscle area, EMA), 등지방두께(backfat thickness, BFT), 근내지방도(marbling score, MS) 4개의 다형질 평가모형이며 성장형질과 도체형질에 따라 아래와 같은 수식들로 표현될 수 있다.
(1)
여기에서 Y 는 성장형질과 4개의 도체형질 표현형을 나타내고, X, Z는 각각 고정효과와 임의효과에 대한 계수행렬을 나타내며, β는 검정차수, 출생지, 우사 및 도축일로 구성된 동기군의 고정효과 및 동기군에 대한 개체의 측정일령과 도축일령을 나타낸다. u는 개체의 육종가, e는 개체의 잔차효과를 나타낸다. 개체 및 잔차 분산은 다음과 같이 표현될 수 있다.
(2)
(3)
여기에서 G0, R0, A는 각각 상가적 유전분산-공분산행렬, 잔차분산-공분산 행렬, 그리고 상가적 혈연계수 행렬을 나타낸다. 평가를 위한 유전모수 또한 국가단위 평가체계에서 활용되는 고정된 모수를 기반으로 유전능력이 평가되었다(Supplementary Table 1). 보정계수에 따른 유전관계행렬 구성, 육종가 추정과 정확도의 산출은 MTG2.0 (Lee and Werf, 2006) 프로그램이 사용되었다.
Imputation 및 보정계수 등
육종가 추정 등에 앞서 유전체자료의 임퓨테이션은 Fimpute (Sargolzaei et al., 2014) 프로그램을 활용하여 진행되었으며 기준 칩을 구성하는 두 수는 Hanwoo custom v1, Illumina bovine v2, Illumina bovine v3 각각 7,911, 3,720두 그리고 4,146두로 총 15,777두이다(Table 1). 앞서 언급된 보정계수는 아래의 일반화된 유전관계행렬 구성 수식으로 표현된다.
(4)
위의 수식에서 Aij는 개체간의 유전관계행렬을, d는 대각성분 예측값을, L은 전체 사용된 마커의 수를, Xij는 개체의 마커에 대한 유전자형가를, Pl는 마커에 대한 기준집단에서의 대립유전자 빈도를, 마지막으로 α는 보정계수를 나타낸다. Speed 등(Speed et al., 2012; 2017)에 의해 발표된 보정계수는 연관불평형에 따라 편향될 수 있는 유전적 구조를 조율하는 하나의 변수로 작용된다. 본 연구에서 보정계수는 기본값인 0에서 ± 0.5만큼 입력하여 유전관계행렬을 구성하였으며 Beadchip platform 3가지마다 각각 유전관계행렬을 구성하여 비교하여 총 9세트의 분석 과정이 수행되었다.
우리나라 국가단위 한우 평가체계의 특성상 부계 기준으로 정확도를 비교하기 위해 비교 대상우는 유전체와 혈통 전체 89,303두 중 자손을 6두 이상 보유한 씨수소(Korean proven-bull, KPN) 850두(가축검정기준 제1장 11조 3항)로 한정하였다. 유전체자료의 일반적인 전처리 작업은 Plink v1.9 (Purcell et al., 2007)를 사용하였고 유전능력평가결과의 정렬, 기초통계 및 자료시각화는 R v3.6.3을 사용하였다.
Results and Discussion
Table 2에서는 imputation Beadchip reference platform에 따라 각각의 보정계수(-0.5, 0, 0.5)로 구성된 GRM을 기반으로 유전능력을 추정하였을 때 추정된 육종가 정확도의 차이를 표현하였다. Beadchip platform 간의 정확도 평균의 차이를 비교했을 때 보정계수가 -0.5로 구성된 GRM에서는 모든 형질에 대한 전체 정확도 평균에서 Hanwoo custom v1 platform이 Illumina bovine v2, Illumina bovine v3 platform보다 각각 0.3, 0.2%p 수준으로 높은 것을 확인하였다. 보정계수가 0인 경우에도 마찬가지로 Hanwoo custom v1 platform이 0.2, 0.1%p 정도 높았고 보정계수가 0.5인 경우, Hanwoo custom v1과 Illumina bovine v3 platform 간 정확도의 차이가 없으나 Illumina bovine v2에서 0.1%p 수준으로 정확도가 높았다. 결과적으로 platform 간의 비교에서 보정계수가 -1인 평가결과를 제외하고 모두 Hanwoo custom v1 platform의 정확도가 더 높았음을 확인할 수 있다. 형질을 기준으로 보았을 때에는 등지방두께의 정확도 차이가 Hanwoo custom v1 플랫폼의 모든 보정계수에 대해서 가장 높게 나타났다. 형질 간 정확도의 차이는 형질에 따라 다양한 변수가 작용할 수 있겠으나 본 연구에서는 각 형질 모두 동일한 기록 수를 기반으로 수행되었기에 표현형 기록 수에 따른 영향은 적을 것으로 보인다.
Table 3 에서는 같은 reference Beadchip platform 내에서 설정된 보정계수에 따라 추정 육종가의 정확도 차이를 표현하였다. 모든 platform에서 가장 추정 육종가의 정확도 평균이 높게 나타난 것은 보정계수가 -0.5로 설정된 GRM을 기준으로 유전능력이 추정된 경우였으며 -0.5, 0, 0.5 순으로 정확도가 높은 것을 확인하였다. Hanwoo custom v1 platform에서는 보정계수 -0.5 대비 0, 0.5에서 전체 정확도 평균의 차이가 0.6, 1.6%p 수준으로 차이가 발생하였고 Illumina bovine v2에서는 0.5, 1.4%p, Illumina bovine v3 platform에서는 0.5, 1.4%p로 정확도의 차이를 확인하였다. 형질을 기준으로 보았을 때에는 근내지방도의 정확도가 모든 경우에서 보정계수 -0.5인 경우 가장 높게 나타나는 것을 확인할 수 있다. 이는 육종가 상관(Fig. 1)과도 유사한 경향을 보인다.
Fig. 1. Correlation of marbling score genomic breeding value by Beadchip platform and genomic relationship matrix (GRM) using each scaling factor (SF). v1, Hanwoo custom v1 platform; v2, Illumina bovine v2 platform; v3, Illumina bovine v3 platform; (1), scaling factor 0.5; (0), scaling factor 0; (-1), scaling factor 0.5.
Fig. 1부터 5는 reference Beadchip platform 및 보정계수에 따른 비교 그룹 간의 육종가 상관을 나타내었다. 모든 형질에 대해서 Beadchip platform과 보정계수를 불문하고 0.992 이상의 상관계수를 나타내고 있으며 같은 platform 내에서는 모든 형질에 대해 0.995 이상의 상관 정도를 보인다. 그 외에도 거의 상관 정도가 유사하지만 모든 형질들에서 보정계수 -0.5과 0.5에서는 히트맵에서 나타내는 바와 같이 상관 정도가 떨어지는 경향을 보였다.
Fig. 2. Correlation of adjusted 12 months weight genomic breeding value by Beadchip platform and genomic relationship matrix (GRM) using each scaling factor (SF). v1, Hanwoo custom v1 platform; v2, Illumina bovine v2 platform; v3, Illumina bovine v3 platform; (1), scaling factor 0.5; (0), scaling factor 0; (-1), scaling factor 0.5.
Fig. 3. Correlation of eye muscle area genomic breeding value by Beadchip platform and genomic relationship matrix (GRM) using each scaling factor (SF). v1, Hanwoo custom v1 platform; v2, Illumina bovine v2 platform; v3, Illumina bovine v3 platform; (1), scaling factor 0.5; (0), scaling factor 0; (-1), scaling factor 0.5.
Fig. 4. Correlation of eye muscle area genomic breeding value by Beadchip platform and genomic relationship matrix (GRM) using each scaling factor (SF). v1, Hanwoo custom v1 platform; v2, Illumina bovine v2 platform; v3, Illumina bovine v3 platform; (1), scaling factor 0.5; (0), scaling factor 0; (-1), scaling factor 0.5.
Fig. 5. Correlation of backfat thickness genomic breeding value by Beadchip platform and genomic relationship matrix (GRM) using each scaling factor (SF). v1, Hanwoo custom v1 platform; v2, Illumina bovine v2 platform; v3, Illumina bovine v3 platform; (1), scaling factor 0.5; (0), scaling factor 0; (-1), scaling factor 0.5.
위의 정확도를 비교한 결과들로 볼 때 Hanwoo custom v1과 보정계수 -0.5의 GRM 기반 유전능력평가 결과가 타 platform 및 보정계수보다 정확도를 약소하게나마 더 개선시킬 수 있다는 점을 시사하고 있다. 그리고 한우 국가단위 유전능력평가 체계에서 사용하는 reference imputation platform인 Illumina bovine v2와 보정계수 0을 사용한 정확도와 본 연구에서 가장 높은 평균 정확도를 보인 Hanwoo custom v1에 보정계수 -0.5을 사용한 결과와 비교해 봤을 때, 기존의 platform 및 보정계수보다 새로운 platform 기반 유전체자료와 보정계수를 쓸 때 0.8%p 수준으로 정확도를 향상시킬 수 있다는 점을 보여준다. 위와 같은 칩 플랫폼 간의 정확도 차이는 reference Beadchip의 개체 수의 차이에 의함으로 사료된다(Table 1). 그리고 이러한 차이의 발생은 유효한 마커로 갱신된 Beadchip platform 또한 하나의 요인이 될 수 있으며 이는 Hanwoo custom v1의 경우에 기존 상용 Illumina bovine v2 및 v3보다 주요변이정보를 더 많이 포함함과 동시에 동형접합체(monomorphic marker) 개수의 비중을 줄임으로 실제 GRM을 구성할 때 반영되는 마커의 정보가 확연히 달라질 수 있다는 점에서 추정 정확도에 영향을 끼칠 수 있다고 보여진다. 동형접합체의 개수와 희소대립유전자, 실제 평가에 사용된 마커의 개수에서 볼 때 전체 마커의 구성 개수는 비슷한 수준이지만 Beadchip platform 간 동형접합체의 수와 실제 평가에 사용된 마커의 개수에서 차이가 있음을 확인할 수 있다. 그 외에 평가에 동일한 유전체자료 15,777두가 사용되었더라도 각각의 Beadchip platform기반으로 imputation을 거치면서 사용된 platform 구성개체 수에 따른 차이도 고려해보아야 한다. 추가로, imputation을 진행할 때 개체 간의 연관된 haplotype을 추적함에 있어서 Hanwoo custom v1 platform 이 다른 두 platform보다 동형접합체가 아닌, 한우의 질병 및 연관형질 등 한우 가계와 직접적으로 연관된 변이 정보를 더 많이 내포할 가능성이 있어 영향을 끼쳤을 것으로 보인다.
보정계수와 관련하여 선행된 연구(Momin et al., 2021)에서 보정계수 수치에 따른 집단 간의 유전상관 등 그 집단에 적합한 보정계수의 설정에 대한 중요성을 언급하고있다. 이는 평가할 집단에 맞는 보정계수를 찾고 그에 따른 GRM을 구성하는 것이 평가의 정확도를 향상시키는데 충분히 기여할 수 있는 요인임을 시사하고 있다. 본 연구에서는 다른 품종 등의 집단을 구성 및 비교하지는 않았으나 실제 유전능력평가결과의 정확도를 보정계수에 따라 비교함으로써 기존에 사용되던 보정계수 대비 현재의 평가집단을 좀 더 잘 설명할 수 있는 보정계수에 가까워졌다고 판단하였다. 이러한 연구결과는 향후 한우 국가단위 유전능력평가체계의 정확도를 개선하는데 있어 도움이 될 수 있을 것으로 사료된다.
Conclusion
본 연구에서는 한우 국가단위 유전능력평가를 수행함에 있어 사용되는 유전체자료의 imputation reference Beadchip platform과 보정계수에 따른 유전능력평가의 정확도를 비교하고 그 결과를 제시하였다. Beadchip platform에 대한 비교 결과로는 Hanwoo custom v1이 다른 플랫폼보다 모든 형질들에 대한 정확도 평균을 0.1 - 0.3%p 수준으로 증가시킬 수 있음을 밝혔으며 보정계수의 경우, 사용된 평가집단에 한정하여 보정계수를 -0.5로 설정하여 GRM을 구성하고 그에 따른 유전능력평가 결과를 비교했을 때 0.5 - 1.6%p 정도 향상시킬 수 있음을 확인하였다. 이러한 연구 결과는 평가집단의 정보를 더 많이 담을 수 있는 customized Beadchip platform을 사용할 때 더 유리할 수 있다는 점을 시사하며 마찬가지로 평가집단에서 개체간의 유전관계행렬을 더 정확하게 만들 수 있는 보정계수를 찾아 사용해야 한다는 점 또한 내포하고 있다. 국가단위 한우 평가체계를 개선함에 있어 본 연구결과가 정확도 향상을 위한 배경이 될 수 있을 것으로 보인다.
Acknowledgements
본 연구는 농촌진흥청 바이오그린 연계 농생명혁신기술개발사업 “국가단위 가축 유전능력평가용 범용 통계분석 패키지 개발(PJ01609901)” 연구과제의 지원 하에 수행되었음.
Authors Information
Soo Hyun Lee, https://orcid.org/0000-0001-5257-2068
Chang Gwon Dang, https://orcid.org/0000-0003-1026-0167
Mina Park, https://orcid.org/0000-0001-7078-9463
Seung Soo Lee, https://orcid.org/0000-0001-7662-9508
Young Chang Lee, https://orcid.org/0000-0002-7442-4186
Jae Gu Lee, https://orcid.org/0000-0001-5970-8714
Hyuk Kee Chang, https://orcid.org/0000-0001-9763-1842
Ho Baek Yoon, https://orcid.org/0000-0002-5300-6081
Chung-il Cho, https://orcid.org/0000-0003-3810-0632
Sang Hong Lee, https://orcid.org/0000-0001-9701-2718
Tae Jeong Choi, https://orcid.org/0000-0002-8299-9342
