서론
유전능력평가를 위한 전통적인 방법은 표현형 자료(Phenotype), 혈통정보(Pedigree) 및 BLUP을 활용하여 육종가를 추정하였으며, 이를 기반으로 선발(Selection) 및 도태(Culling)를 실시하여 유전적 개량을 진행해왔다. 최근 다양한 SNP 패널(Affymetrix, Illumina, Neogen-GeneSeek 등)플랫폼의 개발 및 판매로 개량을 위한 유전체 정보가 적용될 수 있었다. 유전체 선발(Genomic selection)은 유전체 정보를 활용하여 유전체육종가를 추정 후 선발에 적용하는 방법으로 전통적인 BLUP방법보다 세대간격 단축, 신뢰도 향상 등의 장점이 있다. 유전체 선발을 위한 유전체 육종가 추정에서 가장 큰 문제는 제한된 개체 정보를 이용하여 수많은 SNP 마커의 효과를 추정하는데 있으며, 이를 해결하기 위해 Meuwissen et al (2001)은 SNP 마커 효과를 균일하게 낮추는 Ridge Regression BLUP (RR-BLUP)과 Bayes A, Bayes B 등 Bayesian 방법이 유전체 선발에 적용이 되었다. 이후 효율적인 유전체 선발을 위한 Least Absolute Angle and Selection Operator (LASSO) (Usai et al., 2009), Elastic Net (Zou and Hastie 2005), Support Vector Regression with a linear kernel (Moser et al. 2009; Xu et al. 2018) 및 Bayesian LASSO, BayesCπ 및 BayesDπ (de los Campos et al. 2009; Habier et al. 2011) 등과 같은 다양한 방법들이 제안되었다. 또한, Bayes 접근법과 같은 SNP 마커 효과 추정 방법은 BLUP을 이용해 추정 된 육종가(또는 유전체육종가)와 SNP 정보를 이용하여 추정 된 MBV (Molecular Breeding Value)를 혼합 공식(Kachman, 2008; VanRaden, 2008)을 이용하여 유전체 정보가 더해진 육종가(GE-EBV)를 추정하는 방법과 혈연관계인 NRM (Numerator Relationship Matrix)과 GRM (Genomic Relationship Matrix)를 이용한 H-matrix (IBS, Identitiy by status)를 이용하는 single-step GBLUP (Misztal et al. 2009) 등이 유전체 선발을 위해 이용되고 있다.
따라서, 본 연구에서는 한우 도체형질의 유전능력평가 시 유전체 선발 적용을 위해 다양한 사전 분포 적용이 가능한 Bayes 접근법을 활용하여 유전체 정확도 및 유전체 육종가 신뢰도 추정을 통하여 한우 도체형질의 유전체 선발 기초자료로 활용하고자 수행하였다.
재료 및 방법
1. 공시자료
1) 종속 변수(Response Variable)
유전체 분석을 위한 종속변수는 Kim et al (2022)에서 추정된 유전모수와 육종가를 활용하였다. 종속변수는 추정 육종가 및 신뢰도를 이용하여 부모의 효과가 제외된 DEBVexcPA, 부모의 효과가 포함된 DEBVincPA 및 육종가(Estimated Breeding Value)를 추정하여 사용하였다. 또한, Deregressed EBV는 개체마다 다른 신뢰도를 가지기 때문에, 이러한 이질적 분산(Heterogeneous variance)을 설명하기 위해 아래와 같은 가중치(Weighting Factor) 공식(Garrick et al., 2009)을 이용하여 모형에 적용하였다.
여기서, ri2은 종속 변수의 신뢰도, h2는 유전력, c는 SNP 정보에 의해 설명되지 않는 유전적 분산의 비율(c = 0.4, Saatchi et al., 2012), 역회귀육종가(DEBVincPA, DEBVexcPA)로 전환한 후 신뢰도가 0.01 이하인 개체를 제거하였다.
2) 유전체 자료
본 연구를 위하여 전국의 한우농가에서 다양한 방법으로 채취된 8,413두의 DNA를 추출하였다. 유전체 정보를 생성하기 위해 Axiom Bovine 60k version3 (Affymetrix Inc., 2006) SNP패널을 이용하였다. 유전체 분석의 정확성을 높이기 위한 Quality Control은 SNP가 성염색체에 상에 있거나 염색체상의 위치(Position)가 확인되지 않은 SNP를 제거하였으며, SNP Call rate가 95% 이하, Minor Allele Frequency (MAF)가 0.01 이하, Hardy-Weinberg 불평형(HWE)을 측정한 Chi-square (χ2) 값이 95%를 초과한 마커를 제거한 총 64,973개의 SNP 마커를 이용하였다. 유전체 정보를 가진 8,413두 중 중복개체, 혈통오류, 종속변수와 매칭이 되지 않는 개체는 제거하여 유전체 정보가 있는 참조 축군 6,616두를 분석에 이용하였다. 분석에 사용된 보증씨수소는 정액을 이용하여 DNA추출을 진행하였으며, 거세우는 축산물품질평가원에서 제공받은 조직 샘플을 이용하여 DNA추출을 실시하였고, 암소의 경우, 모근을 채취하여 모근에 있는 DNA를 추출하여 분석에 이용하였다.
3) 혈통 및 도체자료
혈통자료의 경우 6,616두와 관련된 5계대의 자료를 한국종축개량협회에서 수집하였다. 분석에 이용된 자료는 총 5,153,168두에 대한 혈통자료를 이용하였다. 도체자료는 5,153,168두 중 도체
자료를 가지고 있는 개체 중 결측치 및 이상치를 사전제거 한 후 남은 2,376,865두에 대해서 축산물품질평가원에서 자료를 수집하여 분석에 이용하였다.
이용된 형질은 도체중(Carcass Weight, CW), 등심단면적(Eye Muscle Area, EMA), 등지방두께(Backfat Thickness, BF) 및 근내지방도(Marbling Score, MS)의 4개 형질을 고려하였다.2. 통계 분석
1) Bayesian Method
유의미한 SNP 마커의 효과를 추정하기 위하여 Bayes 접근법 중 BayesB, BayesC 방법을 이용하였으며, 다음의 통계모형이 SNP 마커 효과를 추정하기 위해 이용되었다.
여기서, yi는 종속 변수(Deregressed EBV, EBV), μ는 전체 집단의 평균, k는 SNP 마커의 수, Zij는 i번째 개체의 j번째 마커의 대립유전자 상태(0, 1, 2), uj는 SNP 마커 효과이며, δj는 SNP 마커의 모델 상 존재 유무(0 또는 1)이다.
BayesB 방법은 각 SNP 마커 마다 다른 분산을 가지게 되며, π값의 설정에 따라 SNP 마커 효과를 0으로 가정하는 혼합 모델(Mixture model)을 이용하며, SNP 마커 효과에 대한 사전 분포(Prior distribution)는 t-분포(t-distribution)를 사용한다. 또한, BayesC 방법 역시 BayesB 방법과 같이 π값의 설정에 따라 SNP 마커 효과를 0으로 가정하는 혼합 모델을 이용하며 SNP 마커 효과에 대한 사전 분포는 정규분포(Normal-distribution)를 사용하며, SNP 마커는 등분산(Common variance)을 가지는 가정으로 샘플링 되어진다(Haiber et al., 2011). 모든 절차는 Gensel4R(Garrick 및 Fernando., 2013) 프로그램을 이용하여 추정하였다.
2) 유전체 정확도 (Accuracy of Molecular Breeding Value )
유전체 정확도를 추정하기 위하여 5개의 집단(5-fold) 교차 검증을 이용하였으며, 교차 검증은 K-means, 임의 군집 방법을 이용하였다(Figure 1).
여기서, K-means 군집 방법은 혈통 정보를 이용하여 혈연 행렬을 추정한 후 이를 거리 행렬(Distance Matrix, D matrix)로 변환하여 집단 내 혈연관계를 최대화, 집단 간 혈연관계를 최소화하였다. D matrix에 대한 공식은 다음과 같다.
여기서, dij는 i번째 개체와 j번째 개체 사이의 거리, aij는 i번째 개체와 j번째 개체의 상가적 유전 관계(Additive genetic relationship), aii(ajj)는 개체의 NRM의 대각 성분(개체의 혈연계수)이며, 이는 근친 효과를 제거하고 D matrix의 대각 성분을 0으로 만든다.
군집화를 통해 참조 집단(Reference population)은 훈련 집단(Training set) 및 검증 집단 (Validation set)으로 나뉘며, 검증 집단 내 개체의 MBV는 훈련 집단에서 추정된 모든 SNP 효과의 사후 평균(posterior mean) 합으로 계산되어지며, 공식은 다음과 같다.
여기서, MBVi는 검증 집단에서 i번째 개체의 MBV, Zij는 i번째 개체의 j번째 SNP의 유전자형(Genotype, i.e., AA=0, AB=1, BB=2), ûj는 훈련 집단의 McMC 샘플링에서 5,000번의 번인구간 이후의 SNP 마커 효과의 사후 평균(Posterior mean)이다.
각 검증 집단에서 추정된 유전체 육종가(MBV)와 종속변수를 이변량 개체 모형(Bivariate animal model)을 이용한 유전상관을 유전체 정확도(Genomic accuracy)로 정의하였으며, 종속변수(T)와 MBV의 유전상관(rg(T, MBV))은 다음과 모형과 같다.
여기서, β1 및 β2는 고정효과의 벡터(β1은 형질의 평균, β2는 군집화 그룹 수에 대한 효과), α1 은 형질, α2는 MBV의 임의 상가적 유전 효과의 벡터이며, Var(α1) = Gσα12, Var(α2) = Gσα22 및 Cov(α1, α2) = Gσα1α2이며, σα12는 형질의 상가적 유전 분산, σα22는 MBV의 상가적 유전분산(σα22/σα12는 마커에 의해 설명되는 형질의 상가적 유전분산)이며, σα1α2는 형질과 MBV의 유전상관이다. G는 같은 그룹에 있는 개체사이의 혈통 기반의 NRM에 기초한 A의 0이 아닌 요소로 이루어져 있으나, 다른 그룹에 있는 개체들 사이의 공분산은 다음과 같이 0이 된다.
이 접근법은 집단 개수 모두에 걸쳐 정보를 효과적으로 모아 유전모수를 추정한다. ε1 및 ε2는 두 형질에 대해 상호 관련이 없는 임의 잔차 효과의 벡터이며, Var(ε1) = Iσε12 및 Var(ε2) = Wσε22, cov(ε1, ε2) = 0 이며, 여기서 I는 항등 행렬(identity matrix)이며, W는 종속변수의 신뢰도에 기초한 가중치를 포함하는 대각행렬이다(SNP 마커 효과의 추정에 사용된 것과 동일한 가중치). 이변량 개체 모형을 활용하기 위해 ASREML 4.1(Gilmour et al., 2015)을 이용하였다.
3) Genomic Enhanced-EBV 추정
Kachman(2008)에 의해 제시된 “Correlated traits” 방법을 이용하여 BLUP 방법에 의해 추정된 육종가 및 SNP 마커 효과의 합으로 추정된 MBV간 혼합 공식을 통하여 GE-EBV를 추정하였으며, 공식은 다음과 같다.
여기서, Ri2는 i번째 개체의 신뢰도, rg2는 유전체 정확도, MBVi는 i번째 개체의 유전체 육종가이며, EBVi는 BLUP에 의해 추정된 i번째 개체의 육종가이다.
결과 및 고찰
1. 군집화 결과
군집화 방법으로 K-means, 임의 방법 등 여러 가지 방법들이 있는데 이전 연구에서 K-means와 임의 군집화를 이용한 연구와 일치한 결과를 나타냈다는 자료가(Camposet al (2018); Boddhireddy et al (2014)) 있어 K-means 방법을 이용한 5-군집 방법을 이용하여 MBV의 정확도와 유전체 정보가 더해진 육종가(GE-EBV)를 추정하였다. 결과는 각각 Table 1에 나타냈다.
K-means clustering을 이용한 방법은 Random clustering을 이용한 그룹에 비하여 그룹 내 혈연관계가 높게 나타났으며, 이는 K-means 클러스터링은 NRM(Numeric Relationship Matrix)을 이용함으로 혈연관계가 높은 개체가 한 그룹에 편성되기 때문인 것으로 사료된다.(Camposet al (2018); Boddhireddy et al (2014)). K-means 클러스터링은 훈련 집단과 검증 집단 사이의 비율 (그룹 내/그룹 간)이 타 논문보다 높게 나타났는데(Lee et al, 2017), 이는 한 그룹에 편성되는 개체들이 더 높은 혈연관계를 가지는 개체로 편성되기 때문인 것으로 사료되며, 개체 수와 비교한다면 정액의 치우침이 발생했다고 사료된다. 또한, 3그룹에서 낮은 개체 수를 보였는데 이는 기존 집단과 혈연관계가 적은 외부 유전자원이 한 그룹에 편성되었기 때문인 것으로 사료된다.
2. 가중치
Figure 2는 6,619두의 유전체 자료에 대한 가중치 분포를 나타냈다. 0.4부터 2.3까지 다양한 가중치 분포를 나타내고 있다. 신뢰도가 높을수록 가중치 값이 상승하는데 이를 증명하듯, 1.0부터 2,3까지의 가중치를 받은 개체는 전부 한우 씨수소였으며, 0.6에서 0.8 사이를 받은 씨수소가 있는 것도 확인을 했지만 이들은 최근에 보증 종모우 정액으로 사용된 씨수소로 후대가 많지 않아 신뢰도가 낮은 것으로 나타났다. 0.4에서 0.6으로 많이 나타나는 개체들은 주로 암소와 거세우들 이였으며, 거세우들은 0.5~0.6의 값을 나타냈고, 암소는 0.4~1.5로 거세우보다 다양한 경향을 나타냈다. 따라서, 정확도를 이용한 가중치는 암소 및 거세우 자료가 후대검정을 진행한 씨수소들 보다 정확도가 떨어지므로 이를 보정하기 위하여 가중치를 이용하여 역회귀육종가를 구하였다.
3. MBV 정확도
Table 2에서는 Axiom 플랫폼을 이용하여 종속변수인 역회귀육종가와 육종가들을 가지고 Bayesian B와 C 방법으로 분포가 0.99일 경우에 대하여 정확도를 나타냈다. 역회귀육종가와 육종가 사이의 정확도 차이를 보면 Bayes B 방법은 근내지방도에서 12%의 차이를 보였으며 도체중에서 26%의 큰 차이를 보였다. Bayes C방법에서도 근내지방도에서 12%의 차이를 보였으며 도체중과 등심단면적에서 26%의 정확도 차이를 보였다. Lee at al. (2019)을 보면 일반적으로 Bayesian방법에서 역회귀육종가와 육종가 사이의 정확도 차이는 약 6%이하정도의 차이를 보였는데 결과가 많이 차이나는 이유는 첫 번째로 유전체 자료의 참조집단이 부족했으며, 씨수소의 유전체 자료가 암소와 거세우보다 월등히 적었고 무엇보다도 Lee at al. (2019)의 논문에서 사용한 자료는 씨수소의 자료만 이용하였지만 이번 연구에서는 암소와 씨수소의 후대축인 거세우의 자료가 더 많이 사용됨으로 인해서 정확도가 기존 연구보다 많이 떨어진다는 것이라고 사료된다. 추후 유전체 자료가 늘어나고, 씨수소 유전체 자료가 추가된다면 정확도가 상승하고, 농가 또는 지역단위 암소평가에 더 정확한 자료가 될 것이라고 사료된다. Table 3은 앞서 육종가가 유전체 육종가보다 정확도가 떨어진다는 결과가 나왔으므로 역회귀육종가와 육종가인 종속변수와 Bayesian 방법에서 편의가 어느 정도 나왔는지 확인하기 위하여 회귀값을 산출하였다. 1의 값으로부터 어느 정도 떨어져 있는지를 나타낸 것으로 Bayes B방법의 역회귀육종가값이 편의가 가장 많이 발생한 것으로 나타났으며, Bayes C방법의 육종가 값이 편의가 가장 적게 발생하였다. Song et al. (2020)에서는 정확도가 더 높게 추정이 되었지만 편의가 많이 발생한 반면, 이번 연구에서는 정확도가 낮게 추정되었지만 편의가 상대적으로 적게 발생하였다. 반면, Bayes B방법보다 Bayes C방법이, 육종가값이 역회귀육종가보다 편의가 적게 발생하였다는 점은 동일한 것으로 나타났다. 앞선 결과와 마찬가지로 유전체 자료의 참조집단의 크기가 작으며, 암소 및 거세우와 혈통으로 연결되어 있는 씨수소의 자료가 추가된다면 농가 또는 지역단위 암소평가에 이용할 수 있는 자료로 이용할 수 있을 것이라 사료된다.
4. 육종가 순위변화
Table 4는 ssGBLUP 방법으로 구한 유전체육종가와 유전체 정보가 포함된 육종가에 대하여 육종가의 순위를 나타냈다. 개체의 상태를 나타내는 지표(Identitiy by status, IBS)인 MBV와 종속변수의 Blending에 의해 GE-EBV를 구하였다. GE-EBV의 순위상관은 Bayes C (π=0.99, Response-variable=EBV, Cross Validation=5-fold Clustering)와 Single-Step 방법으로 추정한 유전체육종가의 상관을 나타냈다. 육종가 순위의 변화를 통해 더 높은 신뢰도를 바탕으로 우수한 능력의 한우 선발이 가능할 것으로 사료된다. EBV와 GEBV간의 상관관계를 형질별로 분석하여 Figure 3-6에 나타내었다. Single Step 방법에서는 도체중, 등심단면적, 등지방두께, 근내지방도순으로 0.89, 0.92, 0.91, 0.94의 상관을 보였으며, Two Step 방법에서는 0.93, 0.94, 0.93, 0.93으로 나타났다. 도체중, 등심단면적, 등지방두께에서는 Two Step 방식이 EBV와 비슷하게 나타났으나 근내지방도에서는 Single Step 방식이 상관이 더 높았다. 하지만 이 차이는 극히 미미하므로 큰 차이를 보이지 않는다고 사료된다.
