Transcriptome Analysis of Longissimus Tissue in Fetal Growth Stages of Hanwoo (Korean Native Cattle) with Focus on Muscle Growth and Development

한우 태아기 6, 9개월령 등심 조직의 전사체 분석을 통한 근생성 및 지방생성 관여 유전자 발굴

Abstract

The prenatal period in livestock animals is crucial for meat production because net increase in the number of muscle fibers is finished before birth. However, there is no study on the growth and development mechanism of muscles in Hanwoo during this period. Therefore, to find candidate genes involved in muscle growth and development during this period in Hanwoo, mRNA expression data of longissimus in Hanwoo at 6 and 9 months post-conceptional age (MPA) were analyzed. We independently identified differentially expressed genes (DEGs) using DESeq2 and edgeR which are R software packages, and considered the overlaps of the results as final-DEGs to use in downstream analysis. The DEGs were classified into several modules using WGCNA then the modules' functions were analyzed to identify modules which involved in myogenesis and adipogenesis. Finally, the hub genes which had the highest WGCNA module membership among the top 10% genes of the STRING network maximal clique centrality were identified. 913(6 MPA specific DEGs) and 233(9 MPA specific DEGs) DEGs were figured out, and these were classified into five and two modules, respectively. Two of the identified modules'(one was in 6, and another was in 9 MPA specific modules) functions was found to be related to myogenesis and adipogenesis. One of the hub genes belonging to the 6 MPA specific module was axin1 (AXIN1) which is known as an inhibitor of Wnt signaling pathway, another was succinate-CoA ligase ADP-forming beta subunit (SUCLA2) which is known as a crucial component of citrate cycle.

동물의 근섬유는 배아기와 태아기를 거치며 형성하게 되며 출생 후에는 상처 치유를 위한 것 외에 근섬유 수를 늘리는 순수한 근섬유 형성은 없으며, 이미 존재하고 있는 근섬유의 비대로 근육의 성장이 이뤄진다. 따라서 태아기의 근육의 성장과 발달이 성체의 근육량 및 조성에 미치는 영향이 매우 크며 이 시기에 발현되는 유전자 및 기능을 구명하는 것은 최종적으로 육질, 육량에 개선시키기 위한기초 자료로 활용될 수 있을 것이다. 하지만 한우에서의 연구는 전무한 실정이다. 본 연구는한우 태아기 성장 단계별 근육의 성장과 발달에 관여하는 유전자를 찾기 위한 전사체 분석을 수행하였다. 한우 태아기 6, 9개월령 등심 조직 시료에서 생산한 전사체 자료를 대상으로 DESeq2와 edgeR을 활용하여 성장단계별 유전자의 발현량을 분석하여 차등발현유전자군을 추출했으며, 2개 소프트웨어서 공통적으로 추출된 유전자군(6개월령 특이 발현 유전자 913개, 9개월령 특이 발현 유전자 233개)을 차등발현유전자로 구명 하였다. 차등발현유전자군으로 분류하였다. 차등발현유전자군을 활용하여공발현 유전자 네트워크 분석을 구성하였으며, 유사한 발현 양상을 보이는 유전자들을 그룹화하여 6개월령 특이 발현 유전자군 5개, 9개월령 특이 발현 유전자군 2개의 모듈로 분류했다. 각 모듈은 Gene Ontology (GO) 및 KEGG pathway 분석으로 유의한 기능을 확인하였다. 그 결과, 한우 태아기 6, 9개월령 특이 발현 유전자 네트워크 중, 근육과 지방생성 대사회로와 관련된 2개의 모듈에 대해 네트워크 내에 허브 유전자를 선정할 수 있었다. STRING을 활용하여 단백질 상호작용 네트워크를 구성하고, MCC (maximal clique centrality) 점수를 활용하여 상위 10%의 유전자들을 공발현 분석의 모듈내 허브 유전자로 선정하였다. 그 결과 6개월령 특이 발현 유전자군의 모듈에서는 axin1(AXIN1) 유전자, 9개월령 특이 발현 유전자군 모듈에서는 succinate-CoA ligase ADP-forming beta subunit(SUCLA2) 유전자가 허브 유전자로 확인되었다. AXIN1 유전자는 선행 연구를 통해 6개월령에서 9개월령으로 넘어가면서 근섬유 수의 증식이 억제되고 지방생성이 활발히 이뤄지는 것에 핵심적인 역할을 하는 것으로 추정할 수 있었다. 또한, 시트르산 회로의 중요 요소인 SUCLA2 유전자는 소의 태아기 지방 조직 성장단계에 따라 유전자의 발현이 증가된다는 보고에 따라, 지방 대사와 관련된 유전자임을 알 수 있었다. 추후 한우 태아기 6, 9개월령에 특이적으로 발현된 유전자들을 대상으로 근육 및 지방 형성 관련 기능을 검증하는 후속 연구가 필요할 것이다.

서론

소의 근육은 고기로 이용되고, 농가의 수익과 직결되는 생시체중, 이유시체중, 도체중, 근내지방도 등과 같은 경제형질과 밀접하게 연관이 되어있다. 따라서 소의 성장에 따른 근육의 성장 및 발달과정을 살펴보는 것은 개량과 사양 연구 분야에 매우 중요하다. 소의 근생성(myogenesis)은 정자와 난자의 수정 후 약 2개월 간의 배아기를 거치면서 중간엽 줄기세포(mesenchymal stem cells, MSCs)가 증식, 분화, 발달하여 1차 근섬유(Primary myofibers)를 형성하게 되고, 이후 약 7~8개월까지 2차 근섬유(Secondary myofibers)가 형성되어 최종적으로 근육을 이루게 된다[5,44]. 이 과정에서 다양한 유전자들(WNT, PAX3, PAX7, GLI, MRFs, MYF5, MYOD 등)이 발현 및 조절되면서 근육 형성에 직 간접적인 영향을 미친다[13]. 또한, 이 시기에 모든 근섬유가 형성되고 출생 후에는 상처 부위 회복을 위한 것 외에 더 이상의 새로운 근섬유 형성은 없으며 근육의 성장은 근섬유 크기의 증가를 통해 이뤄진다[27, 48, 52]. 배아기에 주로 형성되는 1차 근섬유는 그 수가 상대적으로 많지 않으며, 태아기에 주로 이뤄지고 근형성의 대부분을 차지하는 2차 근섬유가 출생 후 근육의 성장 및 발달에 핵심적인 역할을 하게 된다[13,19].

지방조직은 단순히 동물이 생물로서 에너지를 농축하여 저장하는 역할을 할 뿐만 아니라 주요 내분비 기능을 수행하고[17], 경제동물로서 고기의 풍미와 식감을 증진시키며[40], 경제적 가치를 향상시키는 핵심적인 역할을 한다. 특히 한우의 근내지방도는 경락단가의 90% 이상을 결정[38]할 정도로 소비자 선호도가 크며, 육종 및 사양에 중점이 되는 형질이다. 이러한 근내지방을 형성하는 지방세포(adipocytes)는 중간엽 줄기세포로부터 분화되며 근섬유를 구성하는 근세포(myocytes) 역시 동일하므로 세포 분화과정에서 시기의 차이를 보이며 서로 경합하는 경향을 보인다[14]. 반추 동물의 지방생성(Adipogenesis)은 대략 임신 중기(소의 경우 약 4개월)부터 시작되며, 이 시기는 2차 근생성이 활발히 이뤄지고 있을 시기이므로 어미의 영양상태가 근육과 근내지방 비율에 아주 큰 영향을 미친다[13].

Zhan 등[56]은 태아기 단계 중 수정 후 45, 60, 105일령, 출생 후 3일령 염소의 등심조직에서 성장 단계에 따라 차등 발현되는 유전자를 연구한 결과 골격근섬유(skeletal muscle fiber) 발달에 관여하는 actin alpha 1, skeletal muscle (ACTA1), troponin T1, slow skeletal type (TNNT1) 등의 유전자를 확인하였고, Muráni 등[35]은 돼지 품종(Pietrain, Duroc) 별 태아기 14, 21, 35, 49, 63, 77, 91일령 연구에서 근육 관련 myosin heavy chain 1 (MYH1), myosin heavy chain 2 (MYH2), troponin I2, fast skeletal type (TNNI2) 등의 유전자, Lehnert 등[29]은 소 품종(Wagyu × Hereford, Piedmontese × Hereford)별 수정 후 60, 135, 195일령, 출생 시에 따른 유전 발현 연구를 통해 근육 관련 myostatin (GDF8), 지방 관련 fatty acid binding protein 4 (FABP4), fatty acid binding protein 5 (FABP5) 등의 유전자를 확인하였다. 이처럼 여러 가축에서 태아기 근육의 성장과 발달에 관여하는 유전자들에 관한 연구들이 진행되었지만 한우에 관한 연구는 전무한 실정이다. 본 연구를 통하여 한우 태아기 등심 조직의 유전자 발현 양상을 살펴보고, 근육의 성장과 발달에 관여할 것으로 여겨지는 유전자들을 알아보고자 하였다.

특히, Lehnert 등[29]에 따르면 소는 태아기 3분의 2 지점(약 6개월령)에서 근섬유 수의 고정이 이뤄지고, 출생 준비를 위해 골격근의 극적인 적응을 보인다고 하였다. 따라서 본 연구에서는 한우 태아기 6개월령 및 9개월령의 유전자 발현 양상을 살펴 보고, 근육의 성장과 발달에 관여할 것으로 여겨지는 유전자들을 규명하고자 하였다.

본 연구는 수정 후 6, 9개월령 태아의 등심 조직 추출 후, 전사체 분석을 통해 성장 단계별의 차등발현유전자군(Differentially Expressed Genes, DEGs)을 발굴하였으며, 공발현 네트워크 분석(Co-expression network)을 유전자 상호작용 분석 및 기능 분석을 수행하였다.

재료 및 방법

공시 재료

본 시험에 사용된 공시 축은 국립축산과학원 한우연구소의 임신우를 도축하여태아기 6개월령 6두 및 9개월령 8두를 대상으로 등심 조직를 채취하였다. 동물 실험 수행에 관련된 모든 과정은 국립축산과학원 동물실험윤리위원회 승인(승인번호: NIAS-2018-319)을 받아 진행되었다.

전사체자료 분석 및 차등발현유전자군 발굴

한우 태아기 6, 9개월령 등심 조직은 Illumina HiSeq 0000 플랫폼을 활용하여 전사체 자료를 생산하였다. 전사체 자료의 품질 관리는 Trimmomatic (USADELLAB, version 0.36) [8] 소프트웨어를 활용하여 염기서열 양 말단의 베이스 퀄리티 3 이하와 sliding window로 4개 염기 기준 퀄리티 스코어 평균 30(Q30) 이하이면서, 남은 read의 길이가 최소 55 bp를 넘지 않는 염기서열을 제거하였다. 참조 서열은 Ensembl 데이터베이스의 소(Bos taurus) ARS-UCD1.2 버전을 사용하였다. Salmon (Rob Patro, version 0.14.1) [39] 소프트웨어의 index 및 quant 기능을 통해 염기서열 정렬 및 전사체 발현량을 확보하였다. R 패키지인 tximport (Michael Love, version 1.12.3)을 활용하여 유전자 발현량으로 변환했으며 DESeq2 (Michael Love, version 1.24.0) [30], edgeR (Yunshun Chen, version 3.26.5) [42]을 이용하여 차등발현유전자군을 추출하였다. 차등 발현유전자군의 발굴은 DESeq2와 edgeR을 활용하여 공통적으로 존재하는 차등발현유전자군을 추출하여 추후 기능 분석에 활용하였다. DESeq2 소프트웨어는 발현량 정규화 방법으로 median of ratios방식을 사용한다[1]. 또한 샘플 별 참조유전자의 발현량 통계데이터를 도출하는 함수에서 log2 fold change 한계값(lfcThreshold) 파라미터를 설정할 경우 이 값을 넘지 못하는 데이터를 제거하는 독립 필터링(independent filtering) 결과로만 DEGs를 추출한다. edgeR은 trimmed mean of M values (TMM)을 사용하여 유전자 발현량 정규화를 수행한다[43]. filterByExpr 함수로 무의미한 유전자(min count 10, min total count 15)를 필터링한다. 본 연구에서는 DESeq2, edgeR 각각 정규화를 거친 후 DESeq2는 lfcThreshold 파라미터를 1.5로 설정하고 그 결과에서 Benjamini-Hochberg방법으로 조정된 p-value (FDR) < 0.05, log2 fold change 절댓값 1.5를 초과하는 값을 가진 유전자를 차등발현유전자로 추출하였다. edgeR은 filterByExpr 함수로 필터링된 발현량 데이터를 활용해 동일한 옵션으로 차등발현유전자군을 추출하였다. DESeq2, edgeR에서 공통으로 추출된 차등발현유전자군의 발현 양상은 주성분 분석(Principal Component Analysis, PCA) 및 heatmap분석으로 확인하였다. 최종적으로 두 개의 소프트웨어에서 공통으로 존재하는 유전자들을 대상으로 차등발현 유전자군의 공발현 네트워크 및 기능 분석에 활용하였다.

차등발현유전자군의 공발현 네트워크(Co-expression network) 분석

차등발현유전자군의 공발현 네트워크 분석을 위해 R 패키지인 WGCNA (Peter Langfelder, version 1.68) [28]을 사용하였고, 네트워크 분석 소프트웨어인 Cytoscape (Cytoscape Consortium, version 3.7.1) [46]를 활용하여 시각화 하였다. DESeq2의 rlog 함수로 성장단계별 차등발현유전자군 발현량 데이터를 변환하여 WGCNA 분석의 입력 파일로 활용하였다. 인접행렬을 구하기 위한 네트워크 타입을 signed hybrid (두 유전자간의 상관계수가 양수이면 그 값을 그대로 사용하고, 음수이면 0으로 취급)로, 위상 중첩 행렬(topological overlap matrix, TOM) 타입은 unsigned로 설정하였다. 성장단계별 차등발현유전자군 각각에 대해 ‘pickSoftThreshold’ 함수를 이용하여 파워(soft threshold, power)를 결정, 인접행렬, TOM을 구하였다. 1–TOM 행렬을 euclidean 방식으로 거리행렬(distance matrix)화하여 평균 연관(average linkage) 방식으로 계층적 군집화(hierarchical clustering) 후 최소 모듈 크기 30, deepSplit은 2로 설정한 ‘cutreeDynamic’ 함수를 이용 발현되는 양상에 따라 모듈을 나누고, 각각의 모듈을 색깔로 구분하여 Cytoscape로 성장단계별 차등발현유전자군의 네트워크를 시각화하였다.

공발현네트워크(Co-expression Network) 내 차등발현유전자군의 기능 구명

비슷한 발현 양상으로 그룹화된 차등발현유전자군은 Cytoscape의 플러그인 ClueGO (Bernhard Mlecnik, version 2.5.4)[7]를 사용하여 Gene Ontology (GO) [2]의 Biological process(BP) 기능 분석을 수행하였다. P-value는 0.05, GO Tree Interval 4~10 Level, 최소 유전자 수 2, 전체 유전자 중 4%, Kappa Score 0.4, 나머지는 기본 옵션으로 설정하였다. 또한, R 패키지인 clusterProfiler (Guangchuang Yu, version 3.12.0)[55]를 이용하여 pvalueCutoff = 0.1로 옵션을 설정하여 Kyoto Encyclopedia of Genes and Genomes (KEGG) [26] 경로 분석을 수행하였다.

모듈별 허브 유전자 선정

공발현 네트워크 분석 결과, 근내지방 및 근육 형성 관련 유의한 기능을 보유한 모듈들을 대상으로 각 모듈별 허브 유전자를 찾기 위해 STRING (ELIXIR Core Data Resources, version 11.0) [49]을 활용한 단백질-단백질 상호작용 네트워크로 허브 유전자 후보군을 추출한 후 이 중 공발현 네트워크의 모듈내 지위(module membership)값이 가장 높은 유전자를 허브 유전자로 분류하였다. 모듈내 지위가 높을수록 그 유전자의 모듈에 대한 대표성이 높다는 것을 의미한다. STRING은 Cytoscape에서 플러그인인 stringApp (John Morris, version 1.4.2)으로 Confidence cutoff 0.4, Maximum additional interactors 0의 옵션을 설정하여 각 모듈 유전자들에 대한 네트워크 데이터를 활용하였으며, cytoHubba (Chia-Hao Chin, version 0.1) [12] 플러그인의 MCC (maximal clique centrality)방식으로 유전자별 MCC점수를 뽑아 상위 10% 이내인 유전자를 후보군으로 선정하였다. Chin 등은 효모 단백질-단백질 상호작용 네트워크에서 MCC방식이 허브 유전자를 찾는데 높은 정확도를 보였다고 보고하였다[12].

결과

전사체 분석 및 차등발현유전자군 발굴

한우 태아 6개월령 및 9개월령의 등심 조직 샘플을 품질관리 및 참조서열 정렬 수행 결과, 총 402,325,883개의 리드 중 약 97.54%가 품질관리과정을 통과하였다. Salmon을 이용한 참조 서열 매핑 과정에서 최종적으로 272,605,173개의 염기서열이 정렬되어 원시 데이터 기준 약 67.76%의 매핑 비율을 보였다(Supplementary 1). 샘플별 발현량 데이터의 PCA 결과를 Fig. 1A와 같이 확인할 수 있었다. DESeq2 소프트웨어 분석 결과, 태아기 6개월령에서 특이적으로 발현되는 유전자 1,233개, 9개월령에서 특이적으로 발현되는 유전자 333개로 나타났다. edgeR 소프트웨어를 활용한 분석결과는 6개월령에서 특이적으로 발현되는 유전자 1,598개, 9개월령에서 특이적으로 발현되는 유전자 1,125개를 확인하였다. 두 소프트웨어에서 공통적으로 나타난 차등발현유전자군은 6개월령 특이 발현유전자 913개, 9개월령 특이 발현 유전자 233개가 확인되었다(Table 1). 공통적으로 나타난 1,146개의 차등발현유전자군에 대한 발현 양상은 Fig. 1B과 같이 확인되었다.

Fig. 1. (A) PCA distinguishes the gene expression pattern of samples. (B) Heatmap shows differentially expressed genes expression pattern.

Table 1. Differentially expressed genes (DEGs) between 6 and 9 months in Hanwoo longissimus tissue identified using DESeq2 and edgeR software

6개월령 특이 발현 유전자 913개 중 DESeq2 기준 log2 fold change 차이가 큰 상위 10개의 유전자, 9개월령에 특이적으로 발현되는 상위 10개 유전자는 각각 Table 2, Table 3으로 정리하였다.

Table 2. Top 10 of 6 months post-conceptional age specific differentially expressed genes in Hanwoo longissimus tissue identified using DESeq2

Table 3. Top 10 of 9 months post-conceptional age specific differentially expressed genes in Hanwoo longissimus tissue identified using DESeq2

차등발현유전자군의 공발현네트워크 분석

DESeq2 및 edgeR에서 공통적으로 확인된 6개월령 특이 발현 유전자 913개의 WGCNA 계층적 군집화 후 공발현 네트워크 분석 결과, Fig. 2A와 같은 네트워크를 확인할 수 있었고, Fig. 2B와 같은 구조를 이루며 5개의 모듈(blue, brown, green, turquoise, yellow)로 나뉘었다. 이들 모듈은 각 blue 129, brown 97, green 73, turquoise 308, yellow 93개의 유전자로 구성되었다(Supplementary 2).

Fig. 2. 6 months post-conceptional age specific differentially expressed genes network analysis by WGCNA. (A) The color of each node shows a module to which the node belongs. The size of label and node represents the connectivity of the node. (B) Gene cluster dendrogram. Each color of the dynamic tree cut graph is considered as module except grey. There are five modules in the graph (blue, turquoise, brown, green, yellow).

233개의 9개월령 특이 발현 유전자의 공발현 네트워크는 총 2개의 모듈(blue, turquoise)로 군집화 한 것을 확인할 수 있었으며(Fig. 3B), 네트워크를 Fig. 3A와 같이 시각화 하였다. blue 모듈은 60개, turquoise 모듈은 132개의 유전자로 구성되었다(Supplementary 2).

Fig. 3. 9 months post-conceptional age specific differentially expressed genes network analysis by WGCNA. (A) The color of each node shows a module to which the node belongs. The size of label and node represents the connectivity of the node. (B) Gene cluster dendrogram. Each color of the dynamic tree cut graph is considered as module except grey. There are two modules in the graph (blue, turquoise).

공발현네트워크로 분류된 차등발현 유전자군의기능 분석

clueGO 소프트웨어를 활용한 BP분석 결과 6개월령 특이 발현 유전자군의 blue모듈은 prostaglandin biosynthetic process 등의 19개의 유의한 기능을 보유하고 있음을 확인했으며(Fig. 4A), KEGG pathway 분석에서는 유의한 기능이 없었다. Brown모듈은 transcription from RNA polymerase I promoter 등 8개의 유의한 기능을 보였으며(Fig. 4B), KEGG는 유의한 기능이 나타나지 않았다. green모듈은 35개의 기능을 가지고 있었으며(Fig. 4C), Notch signaling pathway, Axon guidance, VEGF signaling pathway 등의 경로가 유의한 것으로 나타났다(Fig. 5A). turquoise 모듈은 5개의 유의한 기능이 나타났으며(Fig. 4D), Wnt signaling pathway, Axon guidance, Hippo signaling pathway 등의 KEGG pathway 경로를 보였다(Fig. 5B). Yellow모듈은 17개 기능이 유의하였으나(Fig. 4E), KEGG pathway 분석에서는 유의한 경로가 없는 것으로 확인 되었다.

Fig. 4. Gene ontology biological process term enrichment of the 6 months post-conceptional age specific differentially expressed genes. (A) Blue, (B) brown, (C) green, (D) turquoise, (E) yellow module.

Fig. 5. Kyoto encyclopedia of genes and genomes enrichment of the 6 months post-conceptional age specific differentially expressed genes. (A) Green, (B) turquoise module.

9개월령 특이 발현 유전자군 중 blue모듈은 15개의 유전자가 관여하고 있었으며, acetyl-CoA biosynthetic process from pyruvate, oxaloacetate metabolic process 과 같은 에너지 대사 관련 기능들이 나타났다(Fig. 6A). 또한, clusterProfiler 소프트웨어를 활용하여 KEGG pathway 분석한 결과 Pyruvate metabolism, Carbon metabolism, Propanoate metabolism 등의 대사회로가 유의함을 확인하였다(Fig. 6B). turquoise 모듈은 8개의 유의한 clueGO BP term을 확인했으며(Fig. 6C), Calcium signaling 경로 등이 유의한 것을 확인하였다(Fig. 6D).

Fig. 6. Functional analysis of the 9 months post-conceptional age specific differentially expressed genes. (A)(B) Blue module functional enrichment. (C)(D) Turquoise module functional enrichment. (A)(C) Gene ontology biological process term enrichment. Each label shows low level term. Each colors of the bars in the graph represents group of the terms. (B)(D) Kyoto encyclopedia of genes and genomes enrichment.

한우 태아기 6, 9개월령 차등발현유전자군의 허브 유전자 선정

공발현 네트워크 분석을 통해 비슷한 발현 양상을 보이는 유전자군들이 그룹화되어 모듈로 구성 되었고, 각 모듈은 단백질 상호작용 네트워크 분석을 통해 그룹 내에 존재하는 유전자들과 연결 정도를 분석하여 모듈별 허브 유전자를 선정하였다. 6개월령에서 특이적으로 발현되는 유전자 중 turquoise모듈의 308개 유전자를 대상으로 STRING 단백질 상호작용 분석 결과 Fig. 7A와 같은 단백질 상호작용 네트워크가 구축 되었고, 이 중 MCC 점수 상위 10%인 28개의 유전자(Table 4)를 후보로 선정하였다. 이 중 공발현 네트워크 모듈내지위값이 가장 높은 axin 1 (AXIN1)을 허브 유전자로 선정하였다. 9개월령에서 특이적으로 발현되는 유전자 중 blue모듈의 60개 유전자는 단백질 상호작용 분석 결과, Fig. 7B와 같은 결과를 보였으며, 이 중 MCC 점수 상위 10%인 5개의 유전자를 허브 유전자 후보군(Table 5)으로 선정하였다. 이 중 공발현 네트워크에서 모듈내 지위값이 가장 높은 succinate-CoA ligase ADP-forming beta subunit (SUCLA2) 유전자를 허브 유전자로 선정하였다.

Fig. 7. STRING protein-protein interaction network. Yellow nodes with black border are candidates of the hub gene. (A) Turquoise module of the 6 months post-conceptional age specific differentially expressed genes. (B) Blue module of the 9 months post-conceptional age specific differentially expressed genes.

Table 4. Top 10% of STRING maximal clique centrality of 6 months post-conceptional age specific DEGs in turquoise module

Table 5. Top 10% of STRING maximal clique centrality of 9 months post-conceptional age specific DEGs in blue module

고찰

한우 태아기 6, 9개월령 등심 조직의 차등발현유전자군의 발현 양상

9개월령 대비 6개월령 등심 조직에서 발현량이 높게 나타난 유전자들 중 가장 큰 차이를 보인 reticulon 4 receptor like 2 (RTN4RL2)은 중추신경계의 뉴런에서 주로 발현되며 중추신경계의 손상 시 축삭의 성장을 제한하는 것으로 알려져 있다[23]. 그 다음으로 큰 차이를 보인 myomaker, myoblast fusion factor (MYMK)은 근세포 표면에서 발현되어 근세포들이 융합하는데 역할을 하는 것으로 보고 되었으며[33], 세번째로 zinc finger protein, FOG family member 1 (ZFPM1)는 줄기세포의 적혈구 분화를 억제하는 전사 조절 인자로 알려져 있다[53]. 이들 중 MYMK은 근세포들이 서로 융합하여 다핵 세포(multinucleated cells)인 근관(myotube)을 형성하도록 유도한다. 이는 줄기세포에서 근세포로 분화 후 근관, 근섬유, 근육 형성까지의 전체 과정에서 볼 때 실질적인 근생성이 일어나는 단계라고 볼 수 있으며, 한우 태아기 6개월령에 근육 형성이 활발히 일어나고 이후 9개월령에는 줄어들었음을 추측할 수 있다.

9개월령 등심조직에서 특이적으로 발현된 유전자들 중, 6개월령 대비 가장 큰 차이를 보인 유전자는 ENSBTAG00000001858로 novel gene이었고, Basic Local Alignment Search Tool(BLAST) [24] 검색 결과 voltage dependent anion channel 2(VDAC2)와 높은 유사성(query cover 59% percent identity 98.65%)을 보였다. VDAC2는 전압의존적(voltage dependent)이온 채널 단백질인 VDAC2를 암호화하고 있는 유전자이다[32]. VDAC2 단백질은 Bcl-2 homologous antagonist killer(BAK) 단백질을 음성 조절하여 Mitochondrial Apoptosis을조절하는 것으로 알려져 있다[11]. 다음으로 큰 차이를 보인 myosin heavy chain 2 (MYH2)는 myosin의 구성원 중 하나인 myosin heavy chain 2a 단백질을 코딩하고 있는 유전자로 근수축에 관여한다[45]. 세번째로 큰 차이를 보인 ADP-ribosyltransferase 3 (ART3)는 NAD⁺에서 표적 단백질로 ADP-ribose 잔기를 전이시키는 효소를 암호화하고 있는 유전자로 세포분열, DNA 수선, 면역 반응 조절 등에 관여하는 것으로 알려져 있다[47]. 이들 중 MYH2는 빠른 근육으로 분류되는 Type 2A 근섬유를 형성하는 유전자로 태아 단계에서 주로 발현되어 배아 동형체(isoform)를 대체하는 것으로 보고 되었다[45]. 또한, MYH2 외에도 같은 태아 동형체인 myosin heavy chain 1 (MYH1), myosin heavy chain 4 (MYH4) 등이 9개월령 특이유전자로 나타나고, 배아 동형체인 myosin heavy chain 3(MYH3), myosin light chain 족의 배아 동형체인 myosin light chain 4 (MYL4) 등이 6개월령 특이 유전자(Supplementary 2)로 나타난 것은 선행 연구와 유사하게 한우에서도 이 시기에 근육 타입이 변화되고 있음을 시사한다.

차등발현유전자군의 공발현네트워크 분석

일반적으로 연결 정도에 가중치를 적용하기 위한 soft threshold power가 150, 450으로, 보통 1에서 30 사이의 수가 사용되는데 비해 상대적으로 매우 높게 설정되었다. 이는 차등발현유전자군의 한계치를 조금 높게 하여 추출한 결과 매우 밀접한 상관관계에 있는 유전자들로만 집단이 구성된 것으로 사료된다. 즉, 유전자간의 상관 계수가 상향 평준화되어 분포범위가 매우 좁아 모듈화할 수 있는 변별력이 매우 떨어지므로 이를 보정하기 위해 높은 수준의 파워를 설정하였다.

Gene Ontology 활용 공발현네트워크 내 유전자군의 기능분석

6개월령 특이 발현 유전자군 blue 모듈은 peptidyl-lysine hydroxylation, peptidyl-proline modification 등이 유의적으로 나타났다. Peptidyl-lysine hydroxylation과 peptidyl-proline modification 중, peptidyl-proline hydroxylation는 콜라겐(collagen) 생합성 과정과 관련되어 있다고 알려져 있다[36]. 콜라겐은 근육과 뼈를 이어주는 연결 조직인 건(tendon)의 핵심적인 구성원이며, 근육의 발달 과정에 영향을 미치는 것으로 알려져 있다[4,21]. 이 대사회로의 모듈내 참여 유전자는 matrix metallopeptidase 14 (MMP14), prolyl 3-hydroxylase 3 (P3H3), prolyl 3-hydroxylase family member 4 (inactive)(P3H4), procollagen-lysine,2-oxoglutarate 5-dioxygenase 3(PLOD3), FKBP prolyl isomerase 10 (FKBP10)이었으며, 이 모듈내 연결 정도가 가장 높은 유전자는 BAF chromatin remodeling complex subunit BCL7C (BCL7C)로(Supplementary 2), 염색질 재구성에 관여하는 복합체의 구성원에 해당된다. 이 복합체는 전사를 조절할 뿐만 아니라 구성원의 종류에 따라 기능적 다양성이 크며, 쥐의 배아 줄기세포의 유지와 다능성 지속 여부에 관여한다고 알려져 있다[22]. 하지만 BCL7C이 골격근이나 콜라겐 형성에 미치는 직접적인 영향에 관한 연구는 현재까지 없는 것으로 보고 되었다.

turquoise 모듈은 establishment or maintenance of apical/ basal cell polarity기능이 유의한 것으로 나타났다. 세포의 극성은 세포 분열시 방향을 결정하게 되고 골격근에서 정점(apical)과 기저부(basal)의 비대칭 분할은 정점 자식 세포의 근섬유 접촉을 끊어 myogenic regulatory factors (MRFs)의 상향 조절에 저항하도록 만든다[18,51]. 이 대사회로와 관련된 유전자는 discs large MAGUK scaffold protein 4 (DLG4), fascin actin-bundling protein 1 (FSCN1), leucine rich repeats and calponin homology domain containing 4 (LRCH4), protein tyrosine kinase 7 (Inactive) (PTK7), SH3 domain binding protein 1 (SH3BP1), transcription factor 15 (TCF15)였으며, 모듈내 연결 정도가 가장 높은 유전자는 zinc finger protein 444 (ZNF444)였다(Supplementary 2). ZNF444는 돌연변이에 의한 암 발생에 관한 연구가 주를 이뤘으며, 근 발달 과정에 관한 연구는 찾을 수 없었다. 이 유전자의 주요 paralog인 zinc finger protein 202 (ZNF202)는 인간의 지방 대사에 관여하며, 고밀도지단백질 조절자인 ATP binding cassette transporter A1 (ABCA1)의 억제자로 저알파지질단백혈증(hypoalphalipoproteinemia)과 관련이 있는 것으로 보고되었다[25].

9개월령 특이 발현 유전자군의 blue 모듈은 acetyl-CoA biosynthetic process from pyruvate, cellular respiration, tricarboxylic acid cycle 등의 세포 호흡, 에너지 대사 관련 기능이 유의한 것으로 나타났다. Remels 등[41]은 근세포의 분화 과정에서 분화 전 대비 세포 호흡 및 에너지 대사가 증가한다고 보고하였다. 이러한 기능에 dihydrolipoamide S-acetyltransferase (DLAT), dihydrolipoamide dehydrogenase (DLD), phosphoglycerate kinase 1 (PGK1), malate dehydrogenase 1(MDH1), NADH dehydrogenase subunit 4L (ND4L), NADH: ubiquinone oxidoreductase core subunit S1 (NDUFS1), nipsnap homolog 2 (NIPSNAP2), solute carrier family 25 member 12 (SLC25A12), succinate-CoA ligase ADP-forming beta subunit (SUCLA2) 유전자가 관여 되었다. 모듈내 연결 정도가 가장 높은 유전자는 S-phase kinase associated protein 1(SKP1)로(Supplementary 2), SCF ubiquitin ligase 단백질 복합체의 구성원 중 하나다. 유비퀴틴 의존적 단백질분해(ubiquitin dependent proteolysis)에 관여하며 β-catenin 레벨을 조절하는 것으로 알려져 있다[3].

한우 태아기 6, 9개월령 등심 조직 형성 관련 허브 유전자발굴

6개월령 특이 발현 유전자군의 공발현 네트워크 중 turquoise 모듈은 KEGG pathway 분석(FDR adjusted p-value 0.011) 결과, dishevelled segment polarity protein 2 (DVL2), SRY-box transcription factor 17 (SOX17), nuclear factor of activated T cells 4 (NFATC4), C-terminal binding protein 1(CTBP1), transcription factor 7 like 1 (TCF7L1), NKD inhibitor of WNT signaling pathway 1 (NKD1), AXIN1, NKD inhibitor of WNT signaling pathway 2 (NKD2), Wnt family member 4 (WNT4), frizzled class receptor 8 (FZD8) 유전자 등으로 구성된 Wnt signaling pathway가 유의하게 나타났다. Wnt signaling pathway는 동물의 근육과 지방 형성 조절에 매우 중요한 경로라고 보고되었다[15]. 특히 Wnt signaling pathway의 다양한 경로 중 canonical Wnt signaling pathway 즉, β-catenin signaling pathway는 골격근에서 paired box 3(PAX3)와 glioma-associated oncogene (GLI) 계열 인자들을 발현시켜 근생성을 조절하는 것으로 보고되었다[9, 10, 20]. 또한, 말단 분화(terminal differentiation)를 유도하여 지방생성을 일으키는 것으로 알려진 peroxisome proliferator-activated receptor gamma (PPARγ)을 β-catenin이 억제하여 지방생성을 조절하는 것으로 보고되었다[34,37]. 이러한 β-catenin을 분해하는 복합체의 구성원으로써 AXIN1은 β-catenin signaling pathway를 억제하는 것으로 알려져 있다[6,31]. 또한, Figeac과 Zammit (2015)[16]은 C2C12 근세포를 이용한 실험에서 근세포의 증식 과정 대비 분화 과정에서 유의하게 AXIN1의 발현량이 줄었으며, AXIN1을 siRNA 매개 넉다운 시킨 쥐의 위성세포는 증식이 억제되고 조숙한 분화가 유도됨을 발표하였다. 한우 태아기 6개월에서 9개월령으로 성장 하면서 등심 조직의 AXIN1 발현량이 감소한 현상은 태아기 중 후반 기간에 근조직 줄기세포들의 증식이 점차 감소하는 대신 근육, 지방 등의 분화로 이어지고 있을 것으로 추정 된다. Du등(2015)[14]이 소에서 상처 부위 재생을 위한 것 외에 순수한 근섬유 생성은 수정 후 약 7개월령에서 중단 되고, 비슷한 시기에 지방 분화가 점차 증가한다는 결과를 얻어, 본 연구결과와 비슷한 양상을 보이는 것을 확인할 수 있었다. 하지만 근세포의 증식, 분화 과정에서 이 모듈내 유전자들에 미치는 AXIN1의 작용 기작은 명확하지 않으며 AXIN1과 높은 연관성을 보이는 Wnt signaling pathway 대사회로와 근육 세포의 성장, 발달과의 연관성 연구가 추가적으로 진행되어야 할 것이다.

9개월령에서 특이적으로 발현된 유전자군의 공발현 네트워크 중 blue 모듈에서acetyl-CoA biosynthetic process from pyruvate 기능과 관련하여 DLAT, DLD, PGK1, tricarboxylic acid cycle에는 DLAT, DLD, MDH1, ND4L, NDUFS1, SLC 25A12, SUCLA2 등의 유전자가 존재 하였다. Zhou 등(2012)[57]은 지방이 많은 Rongchang돼지 품종에서 지방이 적은 품종인 Landrace보다 MDH1 유전자가 높은 발현량을 보인다고 발표했다. Yin 등(2014)[54]은 사람의 지방전구세포(cell line SW872) 분화 실험에서 분화 유도 후 14일(분화 후기)된 세포의 DLD 유전자가 대조구(분화 전 지방전구세포)보다 10배 이상 차등 발현되었다고 보고하였다. Taga 등(2012)[50]은 수정 후 110, 180, 210, 260일령 소(Charolais, Blond d'Aquitaine) 태아의 신장 주위 지방 조직을 채취하여 단백질 발현 분석 결과 DLAT, DLD, MDH1, NDUFS1, SUCLA2 등이 4개 성장단계에 따라 발현이 증가함을 보고하였다. 이러한 연구 결과들을 기반으로 blue 모듈에 구성되어 있는 유전자군이 한우 태아기 6, 9개월령 지방 형성에 관여한다는 가능성을 제시할 수 있다. 특히, 모듈 내 허브 유전자로 선정된 SUCLA2 유전자가 핵심적인 역할을 할 것으로 추정된다.

감사의 글

본 연구는 농촌진흥청 연구사업(PJ01361502) 연구비를 지원받아 수행하였으며, 연구비 지원에 감사 드립니다.

The Conflict of Interest Statement

The authors declare that they have no conflicts of interest with the contents of this article.