Development and Genetic Diversity Analysis of Microsatellite Markers Using Next-generation Sequencing in Seriola quinqueradiata

차세대 염기서열 분석법을 이용한 방어(Seriola quinqueradiata)의 microsatellite 마커의 개발 및 유전적 특성 분석

Abstract

This study was conducted to develop microsatellite markers in Seriola quinqueradiata using next-generation sequencing. A total of 28,873,374 reads were generated on an Illumina Hiseq2500 system, yielding 7,247,216,874 bp sequences. The de novo assembly resulted in 466,359 contigs. A total of 132 contigs (0.43%), including 60 microsatellite loci, were derived from 30,729 contigs longer than 518 bp. A total of 60 primer sets were designed from the 132 microsatellite loci. A total of 15 polymorphic nuclear microsatellite loci were chosen to evaluate population genetic parameters in the parents and offspring. The mean number of effective alleles was 18.5, ranging from 11 to 30. The observed heterozygosity (H_O) and expected heterozygosity (H_E) ranged between 0.431 and 0.972 with an average of 0.812 and from 0.782 to 0.949 with an average of 0.896, respectively. No significant linkage disequilibrium was observed after Bonferroni revision in any loci. The results show that the 15 polymorphic nuclear microsatellite markers can be used to study the population and conservation genetics of S. quinqueradiata in Korea. To ensure the success of artificial seedling production technology, genetic variations between the parent and offspring populations should be monitored, and inbreeding should be controlled.

본 연구는 차세대 염기서열 분석법(NGS)을 이용하여 방어의 microsatellite 마커를 개발하고, 개발된 마커를 이용하여 방어 집단의 유전적 특성을 분석하기 위해 수행되었다. 차세대 염기서열 분석 장비인 Illumina Hiseq2500를 이용하여 총 28,873,374개의 read들을 얻어 assembly를 수행한 결과, 전체의 약 1.6%에 해당하는 466,359개의 read들이 assembly 되었으며, 이 read들의 총 길이는 7,247,216,874 bp로 확인되었다. 크기가 518 bp 이상이 되는 contig는 30.729개로 나타났으며, 이 중 microsatellite 영역을 포함하는 contig 132개(0.43%)를 1차로 선별하고, PCR 증폭 여부 및 유전자형 분석을 통해 microsatellite 후보 60개를 2차로 선별하였다. 그 중 방어집단의 마커로서 유용한 15개의 microsatellite 마커를 선택하였다. 방어집단을 대상으로 개발된 15개의 microsatellite 마커로 분석한 결과, 관찰된 유효 대립유전자수(N_A)는 평균 18.5(11~30)로 나타났다. 평균 관측치 이형접합도(H_O)와 평균기대치 이형접합도(H_E)는 각각 0.812(0.431~0.972)와 0.896(0.782~0.949)으로 나타났다. 다형성이 관찰된 모든 microsatellite 마커 간의 연관불평형은 나타나지 않았으며, 해산어의 평균 H_E 값인 0.79 이상의 수치를 나타내었다. 따라서 본 연구에서 개발된 15개의 microsatellite 마커는 방어 집단의 유전적 다양성 분석에 유용할 것으로 사료된다.

서론

방어(Seriola quinqueradiata)는 농어목(Perciformes) 전갱이과(Carangidae)에 속하는 어종으로 우리나라를 비롯한 일본, 대만 등 북서태평양에 널리 분포하는 온대성 회유어종이다 [25]. 또한 우리나라와 일본 등에서 인기가 많은 고급어종으로 경제적 가치가 높고, 수요가 지속적으로 증가하고 있어 양식 생산을 위한 종묘생산 및 양성에 관한 연구들이 진행되고 있다[13].

방어의 양식산업을 위한 인공종묘생산 연구는 일본의 경우 1960년대부터 시작되었지만[27], 우리나라의 경우에는 자연에서 치어를 포획하여 양성하는 방법으로 종묘생산 및 중간양성 기술개발이 진행되고 있으나 아직까지 인공종묘생산 연구에 있어 미흡한 실정이다[39]. 방어에 관한 국내 연구로는 제주도 주변해역을 중심으로 한 방어의 분포특성[5,21] 및 어획량 변동과 산란을 위한 회유경로에 관한 연구[24] 등 생태학적 특성에 대한 연구가 대부분이며, 유전학적 통계분석을 통해 방어의 어미집단과 이들로부터 태어난 치어집단을 보다 체계적으로 관리할 수 있는 연구는 아직 이루어지지 않고 있다. 방어 어미집단의 유전적 관리가 중요한 이유는 소규모 양식집단에서 어미들의 유전자형을 고려하지 않은 교배로 극단적인 유전 자형의 손실이 일어날 수 있으며[34], 근교약세(inbreeding depression)로 이어져 생산된 자손의 생존에 악영향을 미칠 수있기 때문이다[22]. 또한, 유전적 다양성의 부족으로 유전적 부동과 병목현상으로 인해 질병 저항성 감소, 기형 증가, 성장 저하, 집단 근친도의 증가 등의 부정적인 영향을 미칠 수 있다[2,9].

대표적 양식어종인 넙치와 전복의 경우에도 여러 세대에 걸쳐 한정된 어미로부터 종묘가 생산되어 왔기 때문에 자연집 단에 비해 양식집단이 유전자 빈도 및 대립유전자 수가 크게 감소된다고 보고된 바 있다[7, 14-16, 26, 28, 31, 35, 37, 40].
따라서 우리나라 방어의 인공종묘생산에 있어 어미집단과 치어집단의 체계적인 관리를 위해 유전적 다양성 유지에 관한 연구가 필요한 실정이다.

이를 위해 최근에는 차세대 염기서열(Next-Generation Sequencing, NGS) 분석법을 이용하여 대상어종의 염기서열 정보를 대량으로 확보하고, microsatellite 영역을 직접 탐색함 으로써 microsatellite 마커를 개발하는 방법이 이용되고 있다.[10,39] 참담치[17], 젓새우[18], 해양미생물[23] 등 다양한 해양생물에서 NGS법을 이용한 microsatellite 마커가 개발되고 있다.

본 연구에서는 방어의 유전적 다양성을 분석하기 위해 NGS 법을 사용하여 microsatellite 마커를 개발하고, 개발된 microsatellite 마커를 이용하여 방어의 어미집단과 치어집단을 대상 으로 유전적 특성을 분석하였다.

재료 및 방법

시료 확보 및 DNA 추출

본 연구에 사용한 방어 시료는 국립수산과학원 제주수산연 구소에서 사육 중인 방어의 어미집단(n=72)과 2017년에 종묘로 생산된 치어집단(n=157)의 지느러미 조직을 채집하여 사용 하였다. Genomic DNA의 순수분리와 정제는 자동화 DNA 추출 장비(MagExtractor MFX-6100, Toyobo, Japan)와 DNeasy® 96 Blood & Tissue Kit (Qiagen GmbH, Germany)를 사용하였다[20]. 분석에 사용할 시료를 20 mg씩 E-tube로 옮겨 DNeasy® 96 Blood & Tissue Kit를 사용하여 제조사의 방법에 따라 ATL Buffer 200 μl와 Proteinase K 10 μl를 첨가한 후, 혼합하여 56℃에서 24시간 동안 반응시켰다. 그 후 AL Buffer 200 μl와 99% ethanol 200 μl을 첨가하여 S-Blocks의 DNeasy 96 plates에 옮겨 4,000 rpm (3,000× g)으로 20분간 원심분리하 였다. Column plates를 교체한 후, AW 1 Buffer 500 μl 첨가하여 4,000 rpm (3,000x g)으로 20분간 원심분리하였다. 동일한 방법으로 column plates를 교체한 후, AW 2 Buffer 500 μl 첨가하여 4,000 rpm (3,000× g)으로 40분간 원심분리하였다. Ethanol을 제거한 다음, AE Buffer 100 μl에 녹여 genomic DNA를 분리하였다. 분리한 genomic DNA는 1.8% (w/v) agarose gel을 이용하여 전기영동을 통해 밴드의 유무를 확인한 후, E- Graph Gel Documentation System (ATTO, Korea)을사용하여 농도를 측정하였다.

Microsatellite 탐색

Microsatellite 마커 개발에 필요한 방어의 염기서열 정보를 얻기 위해 Illumina Hiseq2500 (Illumina, USA)을 이용하여 방어의 지느러미 조직으로부터 추출한 genomic DNA의 염기 서열을 획득하였다. 생성된 raw sequence를 Quality value 20
으로 cutoff하여 trimmed sequence를 구했으며, 서로 중첩되는 구간을 기반으로 paired-end sequence를 합쳐서 single reads를 생성하였다. 이어 생산된 single reads에 대해 CLC assembly Cell (v4.0)을 이용하여 draft assembly를 수행한 후, 생성된 contig에 single reads를 mapping하여 unmapped reads 서열들을 구하였다. Unmapped reads를 대상으로 NICEM_SSR_FINDER를 이용하여 8개 이상의 microsatellite반복 영역을 포함하며 증폭산물의 크기가 150-400 bp 범위의 SSR 마커를 탐색하였다.

탐색된 microsatellite 영역을 증폭하기 위한 프라이머의 설계는 Primer3 [32] 프로그램을 이용하여 제작하였으며, 최적의 프라이머 길이는 18-22 bp 범위, 증폭 산물의 크기는 120-400 bp, T m 값이 50℃-68℃로 설정하고, 프라이머의 위치와 반복되는 염기서열 등을 고려하여 프라이머를 제작하였다.

Microsatellite 유전자좌의 다형성 분석

염기서열 정보를 바탕으로 탐색한 microsatellite 영역을 증폭할 수 있는 프라이머를 선별하기 위해 방어(n=8)를 대상으로 PCR 분석을 수행하였다. PCR 조성은 10 μl당 10 ng 주형DNA, 10X Taq PCR buffer 1 μl, dNTP mixture (10 mM) 0.2 μl, Hs Taq DNA polymerase (2.5 U/μl; TNTResearch, Korea) 0.1 μl, forward와 reverse 프라이머 각각 0.4 μl (10 mM), 주형 DNA 1 μl로 조성하였다. 각 프라이머의 forward 서열 정방향에 형광물질인 HEX, 6-FAM 및 TAMRA을 합성하여 증폭산물 분획 시 결과물이 형광을 띄게 만듦으로써 Genetic Analyzer 3730XL (Applied Biosystems, USA)를 이용한 PCR 산물의 크기 측정을 가능하게 하였다.

PCR 증폭은 95℃에서 7분간 DNA를 사전변성(preincubation) 처리 후, 95℃에서 50초간 변성(denaturation)시켰다. 이어 58℃에서 50초간 프라이머를 annealing 시킨 후, 72℃에서 50초간 DNA 합성(extension) 과정을 총 35회 반복하였으며, 추가적으로 72℃에서 10분간 신장 과정을 수행하였다. 증폭된 PCR 생산물은 1.8%(w/v) agarose gel로 전기영동하여 밴드의 유무와 증폭된 단편 크기를 확인한 후 희석 배율을 결정하였다.

각각의 PCR 생산물은 형광물질 별로 생산물의 크기가 중첩 되지 않도록 고려하여 GeneScan ^TM 400HD ROX (Applied Biosystems, USA) size standard와 Hi-Di ^TM Formamide (Applied Biosystems, USA)를 혼합하여 95℃에서 3분간 변성시킨 후, ABI PRISM 3730XL를 사용하여 단편의 크기를 분석하였다[3].

Microsatellite 유전통계 분석

방어의 어미집단(n=72)과 치어집단(n=157)을 대상으로 다형성이 확인된 유전자좌에 대한 PCR 분석을 수행하였다. 유전적 다양성 분석을 위해 FSTAT version 2.9.3 [12]를 사용하여 대립유전자 수(number of alleles, NA ), 집단 크기를 보정한 대립유전자 수(allelic richness, AR ), 대립유전자 크기(size range, R)를 구하였으며, 집단 간 유전적 분화 정도, Pairwise F_ST [38] 수치를 측정하여 비교하였다. Arlequin version 3.1
software [8] 및 GENEPOP version 4.0 computer package [33]를 사용하여 집단 및 유전자좌위별 이형접합도(heterozygosity)의 관찰치 이형접합률(observed heterozygosity, H_O )과 Hardy-Weinberg Equilibrium (HWE)에서의 기대치 이형접합률(expected heterozygosity, H_E )을 계산하였다. 또한, GENEPO Pversion 4.0 computer package [33]를 사용하여 집단 내 근친 교배계수(inbreeding coefficients, F IS )의 유의성을 확인하였다[38]. 유전자좌의 다양성 정보량(polymorphic information content, PCI)은 Cervus ver. 3.0.7 [19]의 대립인자빈도 분석법을 적용하여 산출하였다. GeneAlEx 6.4 program [30]을 사용 하여 각각의 마커에 걸쳐 전체 집단과 각 개체군의 private allele 빈도 및 대립인자 빈도의 패턴을 살펴보았다. 또한, 유전적 변이 분포를 이용하여 개체 혹은 집단 간의 기하학적 관계를 2차적으로 보여주는 Principal Coordinates Analysis (PCoA) 를 Codominant Genotypic Distance 바탕으로 분석하였다[30].

결과 및 고찰

염기서열 탐색 및 microsatellite 다형성 확인

국립수산과학원 제주수산연구소에서 사육 중인 방어 1개체를 표준시료로 선정하여 차세대 염기서열 분석법인 Illumina Hiseq2500를 이용하여 다량의 염기서열을 획득하였다. 방어의 염기서열을 분석한 결과, 총 28,873,374개의 염기서열 read를 얻었으며, 총 길이의 합은 7,247,216,874 bp이었다. 총 28,873,374개 read의 assembly를 수행한 결과, 전체의 약 1.6%에 해당하는 466,359개의 read가 assembly에 사용되었 으며, 이들 read의 총 길이는 602,570,785 bp로 확인되었다(Table 1).

Table 1. Summary of NGS sequencing results

NICEM_SSR_FINDER를 이용하여 SSR (simple sequence repeat)을 포함하여 프라이머를 디자인할 수 있는 서열을 탐색한 결과, Di-nucleotide가 6번 이상 반복되는 [CA] _n , [GA]_n , [AT]_n , [GC]_n 구간을 포함하는 서열의 총 수는 18,414개이었으며, 비율은 각각 71.3%, 18.7%, 9.8%, 0.2%로 [CA]n 의 반복이 가장 많이 확인되었다. Tri-nucleotide가 5번 이상 반복되는 구간은 총 10,246개이며, 이중 [ATA]_n 이 12.3%로 가장 많았다. Tetra-nucleotide가 5번 이상 반복되는 서열은 총 2,069개로 [AGAT]_n 의 반복이 40.5%로 가장 많이 확인되었다. Di-, Tri-, Tetra-nucleotide의 전체 비율은 각각 59.9%, 33.3%, 6.7%로 Di-nucleotide가 가장 높았으며, Tetra-nucleotide가 가장 낮게 나타났다(Table 2). 염기서열을 분석한 결과, 총 30,729개의 후보 서열 중 반복 염기의 종류, 반복 수, 프라이머의 위치, 증폭 산물의 크기 등을 고려하여 1차 선별한 후보서열 132개에 대해 PCR 증폭 여부 등을 확인하여 microsatellite 마커 60개의 후보를 2차 선별한 후, 형광표지 프라이머를 제작하였다. 유전 자형 분석을 통해 대립유전자 수, 크기 등을 고려하여 최종적 으로 microsatellite 마커 15개를 개발하였다(Table 3).

Table 2. Summary of characteristics of SSR

Table 3. Characteristics of 15 microsarellite loci developed in S. quinqueradiata

Annealing temperature was 58℃ for all markers.

유전적 다양성 및 특성 분석

개발된 15개의 microsatellite 마커를 이용하여 방어 어미집단(n=72)과 치어집단(n=157)을 대상으로 유전적 다양성을 분석하였다. 마커의 다형성 정도를 나타내는 PIC값이 0.5 이상이면 충분한 개체식별력을 갖는 것으로 보고되었는데[4], 본 연구에서 개발된 microsatellite 마커 15개의 PIC값은 0.755-0.939범위로 평균 0.881의 수치로 나타나, 개발된 마커가 상당히 높은 개체식별력을 나타내는 것으로 판단된다(Table 3). 또한 Ohara et al. (2003)가 보고한 microsatellite probe 기반의 microsatellite 마커를 개발하는 연구에 따르면 집단 내의 대립유 전자 수(N_A )는 4-13개이며, H _E 값은 0.85-0.91으로 보고된 것과 비교해 볼 때 방어 어미집단의 대립유전자 수(11-30개)와 H_E값(0.782-0.949)이 비슷하거나 높게 나타났으며, 평균 H_o와 H _E 값이 각각 0.812, 0.896로 크게 차이가 나타나지 않고, 해산어의 평균 H _E 값인 0.79 이상의 수치를 나타내었다[6]. 따라서 15개의 microsatellite 마커는 집단의 유전적 다양성을 확인하는데 유용한 마커로 판단된다.

이에 반해 이들로부터 태어난 치어집단의 유전적 다양성은 크게 차이를 보이는데 치어집단의 평균 유효 대립유전자 수 (N_A )는 6.267개이며, 관측치 이형접합도(H_o)와 기대치 이형접 합도(H_E )는 각각 평균 0.805와 0.651로 어미집단의 유전적 다양성보다 낮은 경향을 나타내었다. 이러한 결과는 치어집단의 모든 유전자좌에서 Hardy-Weinberg Equilibrium (HWE) 이탈이 관찰(p<0.0045, Bonferroni correction) 되는 것으로 확인 되었는데(Table 4), 이는 가계 생산 시 어미집단 중 일부 개체의 관여도가 높다는 것을 의미하며, 이에 따라 유전적 다양성이 감소하였음을 알 수 있다[10]. 양식집단과 같이 한정된 집단 에서 유전적 다양성의 감소는 산란양의 감소와 생존율 및 성장률 저하의 원인으로 알려져 있어[1], 한정된 수의 어미로부터 세대를 거듭할수록 치어집단의 유전적 다양성 감소는 더욱 더 심화될 것으로 사료된다[36]. 유전적 변이의 분포를 통해 집단 간의 기하학적 관계를 보여주는 주성분 분석 결과에서도두 집단의 유전학적 구조는 명확한 차이를 나타내었으며, 각개체들의 특성을 2차원 그래프로 표현했을 때, coord 1축은 전체 데이터 특징의 17.38%를 나타내고 coord 2축은 23.81%를나타내었다. 이 또한 대부분이 한정된 어미 개체 간의 교배를 통해 출생된 개체로서 유전적 차이가 낮은 것으로 분석되었다 (Fig. 1).

Table 4. Genetic variabilities at 15 micrasatellite markers in the parents and offspring population of S. quinqueradiata

N_A=number of alleles per locus; A_R=allelic richness; R=size range of allelesa; H_E=expected heterozygosity; H_O,=observed heterozygosity; P_IC=polymorphic information content; F_IS=in breeding coefficient.

*Not in conformity with Hardy-Weinberg Equilibrium (p<0.003, Bonferroni-corrected value.)

Fig. 1. Principal coordinate analyses (PC OA) of parents and offspring opulations, S. quinqueradiata based on Nei’s genetic distances.

본 연구에서는 방어의 유전적 다양성 분석을 위해 microsatellite 마커를 개발하고자 차세대 염기서열 분석(NGS)을 수행하여 염기서열을 확보하였다. 확보된 염기서열에서 직접 icrosatellite를 추출하는 방법을 통하여 132개의 후보 microsatellite 마커를 선정한 후, 유전자형의 분석을 통해 최종적으로 15개의 다형성 마커를 개발하였다. 개발된 마커로 방어의 어미집단과 치어집단을 대상으로 적용한 결과, 충분한 다형성이 나타났으며, 또한 연관불평형성이 나타나지 않아 방어집단의 유전적 다양성의 분석에 유용한 마커로 확인되었다. 방어에 대한 상업적 관심이 높아짐에도 불구하고, 유전적 정보의 분석이나 다양한 microsatellite 마커가 개발되지 않은 상황을 고려할 때 방어의 유전자원 보존을 위해 보다 많은 microsatellite 마커의 개발이 필요하며, 더 나아가 지속적인 방어의 어미집단과 치어집단의 유전학적 통계분석을 통해 한정된 수의 어미로부터 치어집단의 유전적 다양성의 결핍을 예방할수 있을 것으로 사료된다.

감사의 글

이 논문은 2019년도 국립수산과학원 수산과학연구사업 수산생명자원의 확보, 분석 및 통합관리(R2019033)와 수산유전 자원의 탐색 및 활용(R2019030)의 지원으로 수행된 연구입니다.

The Conflict of Interest Statement

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