Introduction
한우 산업의 기반이자 생산력의 중심이 되는 암소의 개량은 1970년대부터 기관, 지자체에서 혈통기록, 발육성적, 번식능력 등을 조사하여 관리하였으며, 2010년부터 실시한 축산물이력제 도입으로 축산물품질평가원과 한국종축개량협회에서 개체이력번호 조회를 통해 표현형정보와 정확한 혈통정보 수집이 가능해져 유전능력평가를 위한 기초가 마련되었다(Kim et al., 2020). 전국 암소의 도체중과 근내지방도를 살펴보면, 2017년 352.5 kg, 4.2점에서 2021년 368.5 kg, 4.4점으로 조금씩 개량이 이루어지고 있는 것을 확인할 수 있다(KIAPOE, 2021). 반면, 고능력 송아지 생산 및 효율성에 중요한 영향을 미치는 번식형질 개량에 대한 연구로 번식우 농가에 대한 번식 장애 조사, 환경효과 및 농후사료 급여수준에 따라 번식형질에 미치는 영향 등 분석 등 지속적으로 이루어지고 있으나(Choi et al., 2006; Cho et al., 2016; Jeong et al., 2016; Yang et al., 2017), 번식과 관련된 형질의 평균 유전력은 0.069로 도체형질의 평균 유전력인 0.382에 비해 상대적 낮아 개량에 큰 어려움이 따르는 것으로 판단된다(Utrera and Van Vleck, 2004; Mahnani et al., 2022).
최근 SNP (Single Nucleotide Polymorphism)을 이용한 상용칩이 개발되어 단시간에 대용량 유전체정보 확보가 가능해졌으며, 특히, 식물, 인간 및 동물의 주요형질에 영향을 주는 QTL (Quantitative Trait Loci) 또는 다유전자(polygenic gene) 형질에 대한 후보유전자를 식별하기 위한 GWAS (Genome-Wide Association Study) 연구가 활발히 이루어지고 있다(Kim et al., 2022). 국외의 선행연구를 살펴보면, 홀스타인의 유지방률, 유량, 임신율, 수태율과 관련된 여러 SNP 탐색을 성공하였으며, 현장에 적용하기 위해 추가적인 연구가 지속적으로 이루어지고 있다(Cochran et al., 2013; Penagaricano et al., 2012). 국내의 경우 한우의 번식형질과 관련된 연관 SNP 탐색 및 GWAS 분석에 대한 연구보다는 번식우 농장에 대한 실태 조사 및 형질에 영향을 미치는 환경적 효과에 대한 분석이 주로 이루고 있어 GWAS를 이용한 연구가 부족한 실정이다. 따라서, 본 연구는 암소의 번식능력에서 초산우와 경산우를 대표하는 단일형질과 반복형질을 구분하여 각 형질에 대해 GWAS 분석을 실시한 후 관련된 유의적인 후보유전자를 탐색하였다.
Materials and Methods
공시재료 및 표현형정보
본 연구에 사용된 검정집단은 경상국립대학교 기업 GAST에서 제공받은 경남지역한우 암소 1,077두를 이용하였으며, 검정집단의 모근, 혈액 및 개체이력정보가 수집되었으며, 축산물품질평가원에 개체이력번호 조회를 통해 표현형정보를 수집하였다. 번식형질 중 초산우를 대표하는 형질은 초종부일령(Age at 1st service; FS), 초임일령(Age at 1st conception; FCO), 초산일령(Age at 1st calving; FCA)이고 경산우를 대표하는 형질은 임신기간(Gestation length; GL), 공태일수(Day open; DO), 분만간격(Calving interval; CI), 수태당 종부횟수율(Number of service; NS)이며, Table 1의 방법을 통해 기록하였다.이후, 분석의 용이성을 위해 정규분포를 따르지 않는 이상치를 제거하였고, 최종적으로 초산우 형질에는 874두, 경산우 형질에는 741두가 분석에 사용되었다.
유전체정보 및 QC
GWAS 분석에 사용될 유전체정보를 수집하기 위해 개체의 모근과 혈액에서 추출한 고순도 genomic DNA을 사용하였고, Hanwoo 50K SNP BeadChip (Illumina, USA) 분석을 수행하여 상염색체만 포함된 52,195개의 SNP 마커를 확보하였다. 확보된 유전체정보는 GenomeStudio 2.0 program (Illumina, USA)을 사용하여 Plink1.9 program에 적합한 형식으로 변환시켰고, Quality control을 위해 Plink1.9 program을 사용하여 MAF (minor allele frequency)가 1% 미만, Missing genotype은 10% 이상, HWE (Hardy-Weinberg Equilibrium)은 10-6 미만의 SNP 마커를 선발하였다. 위 옵션 순으로 초산우 형질에서는 9,785개, 982개, 3,526개가 제거되어 최종적으로 37,902개의 SNP 마커가 분석에 사용하였고, 경산우 형질에서는 10,037개, 1,002개, 3,155개가 제거되어 최종적으로 38,001개의 SNP 마커가 분석에 사용하였다.
GWAS model
최종적으로 선발된 SNP 마커를 GEMMA v0.93 (Zhou and stephens, 2012) program에서 제공되는 일변량선형혼합모형 옵션을 이용하여 각 형질에 대하여 GWAS 분석을 수행하였다. GEMMA를 이용하면 집단 및 샘플의 구조를 보정한 GWAS를 수행 할 수 있으며, 통계모형은 다음과 같다.
y = Wα + xβ + u+ ϵ; u ~ MVNn (0, λτ − 1K), ϵ ~ MVNn(0, τ − 1In)
여기서, y는 형질 측정치 벡터이고, W는 고정효과(생년월일, 산차) 공변량 행렬이며, α는 절편을 포함하는 해당 계수 벡터이다. X는 유전자형 벡터이고, β는 각 SNP 마커의 효과의 크기를 나타낸다. u는 임의 효과 벡터로 다변량정규분포를 따르는데, 여기서 τ-1는 잔차분산, λ는 τ-1와 잔차 간의 비율을 나타내며, K는 선발된 SNP 마커 중 상염색체로 계산된 유전체관계행렬을 나타낸다. ϵ는 잔차벡터로서 다변량정규분포를 따르며, In는 항등행렬을 나타낸다. 이후, Genome-wide 임계수준(suggestive threshold level)은 Bonferroni method (1 / number of SNP marker)를 이용하였고, 각 형질에 따라 임계수준을 따르는 SNP 마커를 NCBI database (http://www.ncbi.nih.gov)를 이용하여 후보유전자를 탐색하였다.
Results and Discussion
기초통계량 분석
본 연구에 사용된 번식형질은 초산우를 대표하는 형질인 초종부일령(FS), 초임일령(FCO), 초산일령(FCA)과 경산우를 대표하는 형질인 임신기간(GL), 공태일수(DO), 분만간격(CI), 수태당 종부횟수율(NS)로 나누어 표현형정보를 Table 2에 나타냈다. 각 형질별 평균 및 표준편차를 살펴보면, 초산우 형질에서 초종부일령은 409.7±60.7일, 초임일령은 432.4±71.8일, 초산일령 720.2±71.8일로 나타났고 경산우 형질에서 임신기간은 288.9±4.35일, 공태일수는 98±41.4일, 분만간격은 387.1±41.2일, 수태당 종부횟수율은 0.7±0.2%로 나타났다. 선행연구를 살펴보면, Song. (2020)은 암소검정사업에서 수집한 암소에 대한 연구에서 초임일령과 초산일령이 각각 483.26±94.92일, 770.44±95일로 보고하였고, Song. (2017)은 한우개량사업소에서 제공한 98개의 육종농가에서 수집한 번식자료에서 초종부일령은 463.87±73.914일, 초임일령은 470.17±74.498일, 초산일령은 755.70±75.104일로 보고하였다. 두 선행연구에서 보고한 초산우 형질의 값이 본 연구보다 높게 나타났다. Song. (2017)은 임신기간, 공태일수, 분만간격에 대해 각각 286.23±6.779일, 98.990±86.300일, 385.69±86.472일로 본 연구와 부합하였다. Kwon et al. (2019)는 2010년에서 2017년까지 전북에서 수집한 한우 암소에 대한 연구에서 임신기간과 공태일수를 각각 287.602±4.797일, 88.433±42.194일로 보고하였는데, 공태일수가 본 연구보다 낮게 나타났다. 본 연구결과와 선행연구 간의 기초통계량의 차이는 각 지역마다 사양조건, 농장 위치, 사육환경 차이로 나타난 결과로 판단된다.
GWAS 결과
번식형질에 대한 GWAS 분석을 실시하였고 다중 비교 시 생길 수 있는 오류를 보정하는 방법인 Bonferroni correction (Bland and Altman, 1995)을 Manhattan plot에 적용시켰으며, 사용된 임계점은 초산우 형질은 2.64×10-5 (1/37,903개의 SNP 마커; suggestive threshold level)이고, 경산우 형질은 2.63×10-5(1/38,011개의 SNP 마커; suggestive threshold level)로 나타났다. 초산우 형질에 대한 GWAS 분석 결과를 Manhattan plot과 Q-Q plot으로 Figure 1에 나타냈으며, 각 형질에 따라 유의적인 SNP 마커를 Table 3에 나타냈다. 결과를 살펴보면, 초임일령과 초산일령에서는 유의적인 SNP 마커가 나타나지 않았으나 초종부일령에서 12번 염색체에 ARS-BFGL-BAC-830 (p = 2.44×10-5), 13번 염색체에 BTB-00532606 (p = 1.41×10-5), 24번 염색체에 ARS-BFGL-NGS-33785 (p = 2.08×10-6)가 동정되었다. 각 SNP 마커와 가장 근접해 있는 유전자는 KLHL1 (kelch like family member 1), PCK1 (phosphoenolpyruvate carboxykinase 1), SERPINB10 (serpin family B member 10)로 나타났고 위치상 후보유전자로 선정하였다. KLHL1와 SERPINB10은 유전자의 기능에 대한 정보가 나타나지 않았지만, PCK1는 프로피오네이트를 포도당신생합성 물질로 변환시킴으로 프로게스테론의 활동성을 촉진시키는 기능을 가진다. 증가된 순환 혈장 프로게스테론은 난포발달과 임신율에 영향을 주어 생식 효율성에 도움이 된다고 보고되었다(Bedford et al., 2018; Zhang et al., 2021).
|
Table 3. GWAS results for heifer traits 
|
|
FS: Age at 1st service; Chr: Chromosome; bp: base pair; Allele: Major/Minor; MAF: Minor allele frequency |
경산우 형질에 대해 GWAS 분석 결과를 Manhattan plot과 Q-Q plot으로 Figure 2에 나타냈으며, 각 형질에 따라 유의적 임계수준을 만족하는 SNP 마커를 Table 4에 나타냈다. 결과를 살펴보면, 분만간격과 공태일수에서는 유의적인 SNP 마커가 나타나지 않았지만 임신기간에서 3번 염색체에 BTA-94162-no-rs (p = 2.17×10-5), 20번 염색체에 Hapmap47647-BTA-50780 (p = 2.28×10-5), 수태당 종부횟수율에서 27번 염색체에 ARS-BFGL-NGS-58109 (p = 4.27×10-6), 28번 염색체에 BTB-00973399 (p = 1.42×10-5)로 각각 2개의 SNP 마커가 동정되었다. 초산우 형질과 마찬가지로 동정된 SNP 마커를 NCBI database에서 조회한 결과, 각 SNP 마커와 가장 근접해 있는 유전자는 DEPDC1 (DEP domain containing 1), CDH12 (cadherin 12), FUT10 (fucosyltransferase 10), RAB4A (RAB4A, member RAS oncogene family)로 나타났고 위치상 후보유전자로 선정하였다. DEPDC1, CDH12, FUT10은 유전자 기능에 대한 선행연구가 나타나지 않았지만, RAB4A는 근육 및 지방 세포 생성 과정 중 포도당 수송을 조절하는 유전자이며, 와규에서 근내지방도 생성에 큰 영향을 주는 유전자로 보고되었으나(Hudson et al., 2021), 수태당 종부횟수율과 직접적으로 관련된 선행연구가 없었다.
|
Table 4. GWAS results for cow traits 
|
|
GL: Gestation length; NS: Number of service; Chr: Chromosome; bp: base pair; Allele: Major/Minor; MAF: Minor allele frequency |
초종부일령, 임신기간, 수태당 종부횟수율에서 각각 3, 2, 2개의 유의적인 SNP 마커를 동정되었으며, 가장 근접한 후보유전자의 기능을 확인하였고 초종부일령에서 PCK1, 수태당 종부횟수율에서 RAB4A만이 선행연구가 보고되었다. 두 유전자 중 PCK1는 프로게스테론의 활동성 촉진에 영향을 주는 유전자로 난포발달과 임신율에 영향을 주어 생식 효율성에 도움이 된다고 보고되었지만, RAB4A는 수태당 종부횟수율간의 직접적인 연관을 찾지 못하였다. 앞선 3가지의 형질을 제외한 나머지 4개의 형질에서는 유의적인 SNP 마커가 나타나지 않았는데, 이는 도체형질보다 낮은 유전력과 여러 환경효과에 영향을 받는 번식형질에 대한 원인으로 판단된다(Mahnani et al., 2022). 이후, 환경효과를 고려한 모델의 설정과 검정집단을 지속적으로 늘려 추가적인 연구를 실시한다면, PCK1 유전자와 같은 원인유전자 규명이 가능할 것으로 사료된다.
Conclusion
본 연구에서 경남지역 한우 암소를 대상으로 초종부일령, 초임일령, 초산일령, 임신기간, 분만간격, 공태일수, 수태당 종부횟수율에 대한 GWAS 분석을 실시하였고, 그 결과 초종부일령, 임신기간, 수태당 종부횟수율에서 각각 3, 2, 2개의 유의적인 SNP 마커를 동정하였다. 총 7개의 후보 유전자를 확인하였으며, 그 중 초종부일령에서 PCK1는 프로게스테론의 활동성 촉진에 영향을 주는 유전자로 난포발달과 임신율에 영향을 주어 생식 효율성에 도움이 된다고 보고되었다. 향후 번식형질에 대한 모델 설정 및 검정집단을 지속적으로 늘려 추가 연구를 수행한다면 각 형질에 대한 명확한 원인유전자 규명이 가능할 것으로 판단된다.
