서 론
북방종개(Iksookimia pacifica)는 잉어목(Cypriniformes), 미꾸리과(Cobitidae)에 속하는 참종개속(Iksookimia) 어류 로, Kim et al. (1999)에 의해 기름종개속(Cobitis)의 “Cobitis pacifica” 로 신종 보고 되었으나, 이후 반문 및 형태적 특징을 근거로 참종개속으로 전속되었다(Kim, 2009). 북방종개는 분포구계상 참종개속 어류 중 유일하게 동북한 아지역 (Northeast Korea subdistrict)에 서식하는 우리나라 고유 어 종으로, 동해안의 강릉남대천 이북의 독립 하천들에 서식한 다(Kim, 1997; Kim, 2009). 주로 유속이 느리고 하상이 모래 로 이루어진 곳에 서식하여 자갈과 돌바닥에 서식하는 다른 참종개속 어류와 차이를 보인다(Kim, 1997; Choi and Byeon, 2009; Ko, 2015).
북방종개가 서식하는 동해안 하천들은 길이가 짧으며 각 하천들이 오랫동안 격리되어 온 독립하천들로, 이곳에 서식 하는 어류 종들 또한 고립된 환경에 적응하여 제한된 개체군 의 크기를 유지해 오면서 각 하천의 집단들은 서로 다양한 지리적 유전자 분화(geographical genetic differentiation) 양상을 보여왔다(Jang et al., 2006; Bae and Suk, 2015). 특이적으로, 동해안 지역은 2002년부터 2003년까지 연속적 인 대형 태풍으로 인해 심각한 홍수가 일어났으며 이로 인 해 서식지가 교란되고 홍수 복구과정에서 많은 어종들이 자원복원 차원에서 무분별하게 방류된 바 있으며, 최근에 많은 종들이 유입된 것으로 보고되고 있다(Yoon and Kim, 2004; Lee et al., 2010; Ko and Bang, 2013; Byeon and Oh, 2015).
지금까지 보고된 북방종개에 관한 연구는 형태 분류학적 연구(Kim et al., 1999; Kim, 2009)와 핵형분석(Kim and Lee, 1986), 초기생활사(Lee et al., 2011), 생태적 특성 (Choi and Byeon, 2009; Ko, 2015; Ko and Won, 2016), 난막구조(Kim and Park, 1995) 등이 있으며, 한국의 어류적 색목록집에서 하천공사로 서식지가 교란되어 관심대상종 (Least Concern, LC)으로 평가된 바 있다(NIBR, 2011). 유 전자에 기반한 연구는 계통 유전학적 연구(Šlechtováet al., 2008; Kwan, 2015; Perdices et al., 2016)가 있으나, 집단유 전학적 분석에 대한 연구는 아직까지 수행되지 않았다. 집 단유전학적 연구는 집단별 유전자 다양성 및 집단 구조의 정보 등을 제공하기 때문에, 우리나라 멸종위기어종 및 한 국고유어종 등에서 종 복원 및 보전 관리의 중요한 연구 분야로 활용되어왔다(ME, 2009a, 2009b, 2011a, 2011b; MLTM, 2010, 2011).
따라서 본 연구에서는 한국고유종 북방종개의 집단유전 학적 연구를 통해 동해 각 하천들에 서식하는 북방종개 집 단들 간의 유전자 다양성 및 집단 간 분화 양상을 확인하여 생물학적 특징을 밝히고, 더불어 보전 및 관리방안을 제시 하고자 하였다.
연구방법
1.시료확보
북방종개(Iksookimia pacifica)의 표본은 2014년부터 2015년까지 강원도 일대에서 족대(망목 4×4 mm)를 사용하 여 독립하천 10곳에서 북방종개를 채집하였다(Figure 1, Table 1). 집단유전학적 분석을 위한 표본 수는 Pluzhnikov and Donnelly (1996)와 Hudson (2000), Felsenstein (2006) 이 제시한 통계적 검증력을 고려하여 하천 별로 10개체씩을 무작위로 선정하여 20개의 반수체가 포함되도록 하였다. DNA 분석을 위해 가슴지느러미와 배지느러미 1개씩을 절 단하여 99% 알코올에 고정하였으며, 나머지 표본은 10% 포르말린 수용액에 고정하여 표본을 제작하였다.
2.DNA 염기서열 확보
표본의 genomic DNA는 DNeasy Blood and tissue kit (Qiagen, Hilden, Germany)를 사용하여 추출하였으며, 중 합효소 연쇄반응(polymerase chain reaction, PCR)을 실시 하여 DNA 서열을 확보하였다. 1개의 미토콘드리아 유전 자 좌위 mtCOI (Cytochrome-c-oxidase subunit I)와 세 개 의 핵 유전자 좌위 EGR2B (Early growth response 2B), IRBP (Interphotoreceptor retinoid-binding protein), RAG1 (Recombination activating gene 1)에 대해유전자 서열을 확보하였다. 이때 핵유전자의 경우 2n이므로 개체당 2개의 반수체를 각각 확보하여 분석하였다. 본 연구에서는 사용한 미토콘드리아 유전자는 핵유전자에 비해 빠른 유전자 변이 속도로 인해 비교적 짧은 시간 내에 발생하는 종 내 집단 간 분화 및 유전 변이에 대한 정보를 포함하고 있어 일반적 으로 종내 집단유전학 연구에서 많이 사용되는 반면, 주로 모계유전을 통해서만 다음 세대로 전달되어 부계를 통해 전달되는 유전 정보를 설명하지 못하기 때문에 이러한 단점 을 보완하기 위해 핵유전자 좌위들을 같이 분석에 사용하였 다(Avise, 2004; Beebee and Rowe, 2008). 중합효소 연쇄 반응은 nested PCR 방식으로 실시하였고 각 마커의 PCR 조건은 아래에 자세히 기술하였다. 1차와 2차 PCR의 조건 은 mtCOI 서열을 제외하고 모두 동일하였다. mtCOI 마커 의 조건은 initial denaturation 95°C에서 5분; PCR reaction 35 cycles (denaturation 94°C에서 1분, annealing 50°C에서 1분(2차 PCR에서는 55°C로 1분, extension 72°C에서 1분 30초); final extension 72°C에서 5분이었다. IRBP 마커의 조건은 initial denaturation 95°C에서 4분; PCR reaction 35 cycles (denaturation 95°C에서 40초, annealing 53°C에서 40초, extension 72°C에서 1분 5초); final extension 72°C 에서 7분이었다. EGR2B 서열의 경우에는 IRBP 서열 조건 과 동일하였으나, PCR reaction 35 cycles 중 annealing 단 계의 온도만 55°C로 달랐다. RAG1 마커의 조건은 initial denaturation 95°C에서 5분; PCR reaction 35 cycles (denaturation 94°C에서 1분, annealing 단계에서는 touchdown 방식으로 60°C에서 1분 30초로 시작하여 매 cycle마다 0.5°C씩 온도를 낮춰 진행하였으며, 최종적으로 55°C에 도 달한 이후부터는 54°C로 온도를 낮추어 나머지 cycles을 진행, extension 72°C에서 1분 5초); final extension 72°C에 서 7분이었다. PCR 과정 이후, ABI 3730xl sequencer (Applied Biosystems, Foster City, CA)를 이용하여 DNA 서열을 확보하였으며, Geneious v. 6.1.8 프로그램을 사용 하여 확보한 서열을 수정하였다. 최종적으로 분자분석을 위 해 사용된 염기서열의 길이는 mtCOI의 경우 852bp, EGR2B 는 664bp, IRBP는 662bp, 그리고 RAG1은 790bp이었으며, 핵유전자 서열은 모두 exon 부분이었다. 핵유전자 서열들 의 heterozygotic site들은 PHASE 프로그램을 통해 각각의 단상형(haplotype)으로 분리된 후 집단유전학적 분석을 위 해 사용되었으며(Stephens et al., 2001; Stephens and Donnelly, 2003), 분석에 사용된 모든 염기서열은 GenBank database에 등록되었다(accession numbers: KY651384– KY652079).
3.집단유전학적 분석
집단유전학적 분석에 앞서, Genepop v. 4.2 프로그램을 이용하여 확보한 유전자 좌위들과 10개의 집단들에 대해 하디-와인버그 평형(Hardy-Weinberg equilibrium) 및 핵유 전자 좌위들 간의 연관불평형(Linkage disequilibrium)을 시 험하였다(Rousset, 2008). 각 시험은 핵유전자 서열 자료를 genotype data로 변형시킨 후 유전자 빈도에 근거하여 실시되 었으며, 시험의 통계적 유의성(α-value)은 Bonferroni 교정에 의해 계산되었다. 유전자 다양성 및 유전자 별 대립인자 수는 Arlequin v. 3.5 프로그램을 이용하여 계산하였으며 (Excoffier and Lischer, 2010), 대립유전자 풍부도(Allele richness)는 Hp Rare v.1.1 프로그램의 Rarefaction 기법을 이용하여 계산하였다(Kalinowski, 2005). 대립유전자 풍부 도 계산 시, 각 개체군마다의 대립유전자(allele) 수의 차이로 인해 발생할 수 있는 오류를 제거하기 위해 각 집단마다 계산에 사용하는 대립유전자 수를 10개로 통일시켰다.
각 집단 간 유전적 분화도(Genetic differentiation, FST)는 Arlequin v. 3.5 프로그램을 통해 계산하였으며(Excoffier and Lischer, 2010), 3개의 핵유전자 서열 자료들은 PGDSpider v.2.0.5.0를 이용하여 multi-locus genotype 데이터로 전환된 후 분석되었다(Lischer and Excoffier, 2012). 추가적으로 같 은 프로그램을 이용하여 각 하천 별 지리적 거리와 각 집단 간의 유전적 분화도 간의 상관관계를 보기 위해 Mantel test를 실시하였으며, 각 집단 간 지리적 거리는 각 채집지역의 위도 와 경도를 이용하여 채집지역 사이의 직선 거리를 이용하였다 (GeoDatasource; http://www.geodatasource.com/distancecalculator).Table 2
북방종개 집단의 유전학적 구조를 살펴보기 위해 STRUCTURE v. 2.3.4 프로그램의 Bayesian admixture model을 이용하였으며(Pritchard et al., 2000; Hubisz et al., 2009), 핵유전자 좌위들의 multi-locus genotype 데이터를 사용하였다(Lischer and Excoffier, 2012). K 값의 범위를 1부터 10까지 하여 각 K 값마다 1×105번의 burn-in과 5×106MCMC(MarkovchainMonteCarlo)를 각각 5번씩 반 복하였으며, Evanno et al. (2005)의 ΔK 값에 근거하여 가장 적합한 K 값을 도출하였다(Earl, 2012). STRUCTURE 분석 결과에 기반하여 10개의 집단들을 북쪽과 남쪽 그룹으로 나눈 후, Arlequin v. 3.5를 이용하여 AMOVA를 실시하였으 며(Excoffier and Lischer, 2010), Arlequin v. 3.5과 HapStar v. 0.7 (Teacher and Griffiths, 2011)을 이용하여 ‘Minimum spanning tree’ 방식으로 haplotype network를 그렸다.
결 과
본 연구에서는 북방종개(Iksookimia pacifica) 10개 집단에 대한 유전학적 분석을 위해 1개의 미토콘드리아 유전자 좌위 (mtCOI)와 3개의 핵 유전자 좌위(EGR2B, IRBP, RAG1)를 사용하였다. 3개의 핵유전자 좌위들을 시험해 본 결과, 통계 적으로 유의미한 연관불평형 관계를 보이는 좌위는 없었다 (Bonferroni correction; N=10, α<0.005). 본 연구에서 사용된 10개의 집단들은 자산천(JS) 집단만을 제외하고 모두 하디-와 인버그 평형에서 벗어나지 않았으며(Bonferroni correction; N=10, α<0.005), 특히 각 집단에 대해 3개의 핵유전자 좌위를 별도로 고려하여 시험한 결과, 마찬가지로 IRBP 유전자 서열 이 유일하게 자산천 집단에서 유의미한 이형접합자 감소 (heterozygote deficiency, FIS)값을 보여 하디-와인버그 평형 가설을 기각하였다(Bonferroni correction; N=30, α<0.001) (Table 3).
4개의 유전자 좌위들에 대해 10개의 집단들은 거의 집단 내 1개 이상의 유전자 변이를 보였는데, 예외적으로 양양남 대천(YN)과 상운천(SU) 집단은 다른 하천에 서식하는 집 단임에도 불구하고 20개체 모두 동일한 미토콘드리아 유전 자 서열을 지니고 있었다(Figure 1, Table 3). 그 중에서도 양양남대천(YN) 집단은 미토콘드리아 유전자를 포함한 모 든 유전자에서 공통적으로 가장 낮은 유전자 다양성(gene diversity, H)을 보였다. 평균 유전자 다양성(Mean of H) 역시 양양남대천(YN) 집단이 0.23으로 가장 낮았으며, 북 천(BU) 집단이 0.74로 가장 높았다. 유전자 좌위 별 전체 대립유전자 수는 남천(NA) 집단이 22개로 가장 많은 유전 인자 수를 보유하고 있었으며, 배봉천(BB)과 양양남대천 (YN)이 12개로 가장 적은 수를 보유하였다. Rarefaction 방 법에 기반한 대립유전자 풍부도 수치 역시 양양남대천(YN) 집단에서 2.13으로 가장 낮았으며, 반면 남천(NA) 집단이 4.40으로 가장 높았다.
10개의 북방종개 집단들 간의 유전적 분화도(FST)를 살 펴본 결과(Table 4), 대체적으로 유의미한 유전적 분화가 관찰되었다. 먼저 3개의 핵유전자 genotype data는 북방종 개 10개의 집단이 대부분 서로 유전적으로 유의미하게 분화 되어 있음을 보여주었다. 대체적으로 이러한 유의미한 집단 간 유전적 분화가 미토콘드리아 서열의 경우에도 나타났지 만, 몇몇 인접한 하천 그룹들에서는 낮은 분화가 관찰되기 도 하였다. 특히 북천에서 천진천까지의 집단들(북천(BB), 남천(NA), 오호천(OH), 천진천(CJ)) 사이와 양양남대천에 서 강릉남대천까지의 집단들(양양남대천(YN), 상운천 (SU), 연곡천(YG), 강릉남대천(GN)) 사이가 서로 유전적 으로 분화하지 않았음을 보였다.
집단 간의 유전적 분화 정도가 지리적 거리와 상관관계가 있는지 알아보기 위해 Mantel test를 실시하였다. 그 결과, 핵유전자와 미토콘드리아 유전적 분화도 값 모두 집단 간 지리적 거리와 통계적으로 유의미한 상관관계를 보이지 않 았다(핵유전자: Mantel’s r=0.119, P=0.266; 미토콘드리아 유전자: Mantel’s r=0.259, P=0.075).
STRUCTURE 결과에 따라 북방종개 10개의 집단은 천 진천(CJ)과 양양남대천(YN)을 기준으로 북쪽과 남쪽으로 나뉘어지는 2개의 cluster를 형성하였다(K=2) (Figure 2). 북쪽 cluster에는 배봉천부터 천진천 사이에 위치한 집단들 (배봉천(BB), 자산천(JS), 북천(BU), 남천(NA), 오호천 (OH), 천진천(CJ)) 이 포함되며, 남쪽 cluster에는 양양남대 천부터 강릉남대천까지의 집단들(양양남대천(YN), 상운천 (SU), 연곡천(YG), 강릉남대천(GN))이 포함되었다. 특이 적으로 자산천 집단(JS)의 개체들은 지리적으로는 북쪽에 위치하고 있음에도 불구하고 절반 가량의 개체들이 남쪽 cluster에 속하는 것으로 관찰되었다. 이와 같이 남쪽 cluster 의 유전형을 공유하는 개체들을 좀 더 세밀히 알아보기 위 해 더 큰 K값의 막대그래프를 관찰해본 결과, 그 절반의 이질적 개체들은 양양남대천과 강릉남대천 집단이 형성하 는 cluster에 속하였다(Figure 2).
이러한 STRUCTURE 결과에 기반하여 북쪽과 남쪽 cluster를 각각의 그룹으로 설정하여 AMOVA를 실시한 결 과, 미토콘드리아 서열(mtCOI)의 전체 유전 변이의 분산 가운데 64.25%는 그룹 간 차이에 기인하였고, 16.90%는 그룹 내 집단 간의 차이, 18.85%는 집단 내 개체들 간의 차이에 기인하여, 그룹 간의 차이가 북방종개 집단의 미토 콘드리아 유전자 변이에 가장 큰 영향을 미치고 있었다. 반 면, 핵유전자 좌위(EGR2B, IRBP, RAG1)의 전체 유전 변이 분산 가운데 9.74%만이 그룹 간의 차이에 기인하였으며, 17.82%는 그룹 내 집단 간의 차이, 72.44%로 집단 내 개체 들 간의 차이에서 기인하여, 그룹 간의 차이보다는 집단 내 개체들 간의 차이가 북방종개 집단의 핵유전자 변이에 더 큰 영향을 미치고 있음을 보여주었다. 이러한 경향성은 haplotype network에서도 관찰이 되었는데, 미토콘드리아 단상형들은 북쪽과 남쪽 그룹 간의 지리적 구분이 비교적 명확했지만 핵유전자 단상형들은 지리적 구분이 뚜렷하지 않았다(Figure 3).
고 찰
본 연구에서는 한국 고유종 북방종개(Iksookimia pacifica) 가 서식하는 동해의 독립하천들을 확인하고, 집단들의 유전 자 다양성 및 구조를 살펴보았다. 북쪽으로는 배봉천에서 남쪽으로는 강릉남대천에 이르기까지 전체 10개의 하천에서 북방종개가 서식하는 것으로 확인되었으며, 각 하천들은 서 로 지리적으로 단절되어 직접 바다로 유입되기 때문에 북방종 개 집단들 역시 서식지에 정착한 이후 지리적인 유전적 분화 를 겪어왔다고 예상되었다. 본 연구의 결과는 이러한 가설을 지지하였으며, 더불어 동해 하천들의 지리적, 환경적 특징과 관련된 북방종개 집단들의 유전자 다양성 및 구조를 밝혔다.
1.집단의 유전자 다양성
전체 10개의 북방종개 집단들 중 유일하게 자산천(JS) 집단만이 유의미한 하디-와인버그 평형에서 통계적으로 벗 어났는데, 이와 관련하여 집단 내 근친 교배구조(inbreeding), 비무작위 짝짓기(non-random mating), 및 왈런드(Wahlund) 효과 등이 대표적인 원인으로 알려져 있다. STRUCTURE 결과에 기반했을 때, 이 원인들 중에서 왈런드 효과에 의한 영향이 크게 작용했을 것으로 예상된다(Figure 2). 집단 내에 유전적으로 격리 혹은 분화된 다른 아집단들이 존재할 때 집단의 이형접합체(heterozygosity) 빈도가 예상 값보다 감소 하는 현상(heterozygote deficiency)을 왈런드 효과라 일컫는 데, 자산천의 경우 배봉천 다음으로 북쪽에 위치하고 있음에 도 불구하고 집단의 다수의 개체들이 먼 거리에 위치한 남쪽 의 양양남대천(YN)과 강릉남대천(GN)의 유전자형을 공유 하였다. 이는 자산천(JS) 집단이 인위적 방류 등에 의해 교란 되어 왈런드 효과에 의한 유의미한 이형접합자 감소를 초래하 였을 가능성을 제시한다.
북부동해안은 2002년과 2003년에 연속적인 태풍 “루사” 와 “매미"로 인한 대규모 홍수가 발생했었고 이후 하천복원 공사와 더불어 하천복원의 명분으로 한강에 서식하는 어류 뿐만 아니라 동해에 서식하는 많은 어류가 방류되어 급격한 어류상 변화가 보고되어 왔다(Yoon and Kim, 2004; Kim et al., 2006; Lee et al., 2010; Ko and Bang, 2013; Byeon and Oh, 2015). 북방종개의 직접적인 방류기록은 확인할 수 없었으나 많은 기관 및 단체에서 무분별로 어류를 방류 한 것으로 볼 때 충분한 가능성이 있다고 생각되며, 이러한 인위적 요인이 북방종개 집단의 유전자에 영향을 미쳤을 것으로 의심된다.
각각의 집단들에 대한 미토콘드리아 유전자(mtCOI)는 핵유전자 보다 낮은 유전자 다양성을 보였으며, 특히 양양 남대천(YN)과 상운천(SU) 집단의 20개의 개체들은 서로 다른 하천에 서식하고 있음에도 불구하고 오직 하나의 유전 자 서열만을 가지고 있었다(Table 3). 이러한 다양성의 소실 은 최근 이들 집단에서 일어난 병목현상(bottleneck)과 관련 이 있다고 보여진다.
2.유전적 분화도 및 집단유전학적 구조
전체 10개의 북방종개 집단들은 핵유전자에 기반한 Bayesian 분석에 따라 지리적으로 북쪽과 남쪽 그룹으로 나뉘었다. 이러한 핵유전자들의 집단유전학적 구조는 미토 콘드리아 유전자의 AMOVA 분석에서도 일관되게 관찰되 었다. 그럼에도 불구하고 핵유전자들의 AMOVA 결과는 상대적으로 남북간 그룹들 간에 유전적 변이의 분산이 크지 않음을 보여주었다(9.74%). 이러한 미토콘드리아 유전자와 핵유전자들의 차이는 모계유전을 하는 세포 소기관인 미토 콘드리아 유전자보다 핵유전자가 유효집단크기가 4배나 크 다는 일반적인 이론으로부터 설명될 수 있다(Birky et al., 1983). Birky et al. (1983)에 따르면 유효집단크기가 4배나 작은 미토콘드리아 유전자는 핵유전자에 비해 집단간 격리 의 영향을 훨씬 빠르고 민감하게 받는다는 점이 제시되었는 데, 북방종개의 핵유전자와 미토콘드리아 유전자 역시 이러 한 일반적인 민감도 차이를 나타냈다고 해석된다.
북방종개 집단이 한반도에 정착한 역사를 대략적으로 추 론해 본다면, 먼저 조상형 다형성을 포함한 북방종개들이 북쪽과 남쪽으로 크게 나뉘어 정착하였고, 그 후 이 두 그룹 이 지역 내에 있는 개별 하천들로 격리되면서 추가적인 유 전적 분화가 축적되어 갔다고 해석된다. 북방종개와 달리, 동해 북부하천부터 남부하천까지 넓은 분포를 보이는 잔가 시고기(Pungitius kaibarae)의 경우, 조상집단이 북쪽에서 남쪽 방향으로의 점차적인 이동을 통해 서식지에 정착했을 가능성이 제시되었다(Bae and Suk, 2015). 북방종개는 유 전적 분화도와 지리적 거리 사이에 명확한 상관관계가 없는 것으로 나왔다. 따라서 Bae and Suk (2015)의 연구에서 제 시된 북쪽으로부터 남쪽 방향의 이동 양상은 관찰되지 않았 다. 이러한 차이는 각 어류 종의 진화역사, 생활사 및 이동 능력의 차이와 관련이 있을 것으로 추정된다.
미토콘드리아 유전자(mtCOI)의 경우, 핵유전자에 비해 상대적으로 집단 간 유전적 분화도가 컸지만, 일부 이웃하 는 집단들 사이에서는 분화도가 높지 않은 경우도 관찰되었 다(Table 4). 이러한 경우들은 미토콘드리아 유전자 다양성 이 완전히 소실된 양양남대천(YN)과 상운천(SU) 집단들과 나머지 다른 집단들과의 상호비교에서 주로 관찰되었다. 일 반적으로 개체수 병목현상의 유전다양성 감소 효과는 유효 집단크기가 큰 핵유전자에 비해서 유효집단크기가 작은 미 토콘드리아 두드러진다는 사실이 잘 알려져 있다(Fay and Wu, 1999). 양양남대천과 상운천 집단의 미토콘드리아 유 전자 다양성 감소는 태풍과 같은 외부요인에 따른 서식지 훼손과 그 결과 나타난 개체수 병목현상으로 해석된다.
기존 연구들을 보면, 많은 멸종위기 어류들은 같은 종이 라 할지라도 수계별로 분화가 진행되어 독립 수계별로 관리 방안이 필요하다고 보고된 바 있다(MLTM, 2010, 2011; ME, 2011a, 2011b). 이러한 양상은 본 연구에서도 확인되 는데, 동해안 하천에 서식하는 북방종개는 유전적으로 유사 한 몇몇 하천 그룹들로 구별되었다. 따라서 북방종개의 보 전 및 관리를 위해서는 유전적으로 구별되는 그룹 간의 무 분별한 방류는 반드시 지양되어야 하며 그룹 별로 보존대책 이 필요하다고 판단된다.