종 보존 및 자원관리를 위해 집단유전학적 연구가 크게 활기를 띠고 있다. 집단의 유전자 분석은 한 종의 유전적 다양성을 평가하는데 좋은 방법으로(Crandall et al., 1999) 특히 미토콘드리아 DNA 염기서열 분석이 효과적인 방법 중 하나이다(Buonaccorsi et al., 2001). 미토콘드리아 DNA는 보존영역과 변이영역을 가지는데 이중 특히 치환속도가 빠른 control region 영역이 집단의 유전 변이 분석에 많이 사용된다. 또한, cytochrome b (cytb) 영역도 종간 변이 뿐 아니라 종내 변이에 관한 다양한 연구 보고가 있어(Kim et al., 2006a; López et al., 2006; Murray et al., 2008), 종내 지역 집단을 구분하는데 유용한 마커로 널리 활용되고 있다(Grant and Bowen, 1998; Colborn et al., 2001; Muss et al., 2001; Kim et al., 2006a; Liu et al., 2006a; Liu et al., 2007). 그외 계수 및 계측과 같은 형태 형질에 의한 다변량분석도 오래 전부터 지역 집단을 분리하는데 사용된 전통적인 기법 중 하나이다(Turan et al., 2006; Kim et al., 2008; Murta et al., 2008; Myoung and Kim, 2014).

꼼치(Liparis tanakae)는 쏨뱅이목(Order Scorpaeniformes) 꼼치과(Family Liparidae)에 속하는 어류로 우리나라, 일본, 러시아, 중국 등지에 분포하며(Chyung, 1977; Yamada et al., 1995), 우리나라에서는 “물메기탕”으로 이용되는 대표적인 겨울 어류이다(Kim et al., 1986; Ustadi et al., 2005). 꼼치는 1년생으로 침성부착난을 낳고 죽는데(Kosaka et al., 1971; Ustadi et al., 2005), 이때 통발이나 자망 등에 부착된 꼼치난을 식용으로 이용하기도 한다. 해양수산연구정보포털(http://www.nifs.go.kr)에서 제공되는 어획량 자료에 의하면 꼼치는 2005년 1,184톤이었으나 2011년 2,804톤으로 2배 이상 급격히 증가하여 남획의 가능성이 우려되며, 이를 방지하기 위하여 향후 효율적인 자원관리 방안이 마련될 필요가 있다.

꼼치에 관한 국내 연구로는 난발생과 부화자어의 형태(Kim et al., 1986), 성어의 식성(Huh, 1997), 계절에 따른 출현 특성(Huh and An, 2000; An and Huh, 2002), 유어의 식성(Kwak and Huh, 2003), 자치어의 형태(Kim et al., 2011), 자치어의 식성(Choi et al., 2015) 등이 있다. 국외에서는 선택적인 생활사(Kawasaki et al., 1983), 자치어의 형태(Okiyama, 1988), 형태와 생태(Yamada et al., 2007) 등이 있을 뿐 합리적인 자원관리를 위한 집단유전학적 연구는 전무한 실정이다. 꼼치는 침성부착난을 낳고, 치어기(14 mm 전후) 이후에 배지느러미가 흡반으로 변형되어(Okiyama, 1988) 해류에 의한 확산 정도가 다소 미약할 것으로 추정되므로 지역 집단을 형성할 가능성이 높다고 판단된다. 따라서 본 연구는 우리나라 3개 해역(동해, 서해, 남해)에 서식하는 꼼치를 대상으로 분자 및 형태 분석을 통해 지역 집단간 차이 유무를 파악하고자 한다.

2011년 1월부터 2012년 12월까지 동해(경상북도 포항), 서해(충청남도 보령), 남해(경상남도 진해)에서 총 90개체를 채집하여 형태 및 분자 분석을 실시하였다(Fig. 1). DNA분석용 근육은 99% 에탄올로 고정시킨 후 −20℃에서 냉동 보관하였으며, 어체는 5% 포르말린에 7일간 고정시켰다가 3일간 흐르는 물로 세척한 후 70% 에탄올에 넣어 보관하였다. 본 연구에 사용된 표본은 부경대학교(Pukyong National University, PKU)에 등록, 보관하였다.

계수형질 4개, 계측형질 20개는 Nakabo (2002)의 방법에 의거 측정하였다. 계측형질은 버니어캘리퍼스로 0.1 mm 단위까지 측정하였고, 계측치는 로그 변환 후 주성분분석(PCA) 및 정준판별분석(CDA)을 통해 지역집단의 분리 여부를 파악하였다. 각 지느러미 기조수 및 척추골수는 radiograph (SOFTEX HA-100, Japan)를 이용하여 측정하였으며 비모수검정법인 Kruskal-Wallis test를 통해 지역 집단간 차이 유무를 파악하였다. 모든 통계 분석은 SPSS ver 12.01을 사용하였다.

Genomic DNA 추출

Genomic DNA는 99% 에탄올에 보관된 꼼치의 근육조직을 Chelex 100 resin (Bio-rad, USA) 150 μL에 넣은후, Thermal cycler (Bio-rad MJ mini PTC-1148, USA)를 이용하여 60℃와 99℃에 각각 20분, 25분간 두었다. 추출된 Genomic DNA는 4℃에서 냉장보관하였다.

중합효소연쇄반응

미토콘드리아 DNA cytochrome b (cytb) 영역을 증폭시키기 위해 Palumbi et al. (1996) 의 L-GLUDG (5′-TGA CTT GAA RAA CCA YCG TTG-3′)와 H-CB3 (5′-GGC AAA TAG GAA RTA TCA TTC-3′) primer를 이용하였다. 10X PCR buffer 3 μL, 2.5 mM dNTP 2.4 μL, GluDGD primer 1 μL, CB3 primer 1 μL, Ex-Taq DNA polymerase 0.1 μL (GnP, Korea) 를 섞은 혼합물에 total DNA 3 μL를 첨가한 후, 총 30 μL가 될 때까지 3차 증류수를 넣고 Thermal cycler (Bio-rad MJ mini PTC-1148, USA)를 이용하여 다음과 같은 조건으로 PCR을 수행하였다. Initial denaturation 94℃에서 5분; PCR reaction 35 cycles (denaturation 94℃에서 1분, annealing 59℃에서 1분, extension 72℃에서 1분); final extension 72℃에서 5분 실시하였다. 정제는 ExoSAP-IT (United States Biochemical Corporation, USA)를 이용하였다. 염기서열은 ABI 3730XL sequencer (Applied biosystems Inc., USA)에서 ABI Bigdye terminator cycle sequencing ready reaction kit v3.1 (Applied Biosysterms Inc., USA)를 이용하여 얻었다.

자료분석

미토콘드리아 DNA cytb 영역의 염기서열은 BioEdit ver 7의 ClustalW (Thompson et al., 1994)를 이용하여 정렬하였다. 유전거리는 Mega 5 (Tamura et al., 2011)에서 Kimura-2-Parameter모델(Kimura, 1980)을 이용하여 구하였다. 분지도는 최소근접네트워크(Minimum spanning network)를 이용하여 작성하였고 bootstrap은 1,000번 수행하였다. 외집단은 Liparis liparis (EU492094)로 설정하였다. 지역 집단의 전이율 및 전환율은 Arlequin ver 3.5.1.2 (Excoffier et al., 2005)를 사용하여 계산하였고, 염기다양성(π, Nei and Li, 1979) 및 유전자형다양성(h, Nei, 1987)도 Arlequin 프로그램으로 계산하였다. 지역 집단간 유전적 분화 정도를 나타내는 고정지수(F_ST)는 Rice (1989)의 방법으로 구하였고, 이때 유의성 검증을 위해 1,000번의 bootstrap을 실시하였다. 지역 집단의 유전적 다양성이 중립성과 개체군 평형을 따르는지 알기 위해 Tajima’s D (Tajima, 1989) 및 Fu’s Fs (Fu, 1997) 테스트를 실시하였다. 또한, 집단 크기의 변천 역사는 Rogers and Harpending (1992)의 Mismatch distribution 분석을 통해 추정하였고, 집단의 팽창 시기는 Li (1997), 팽창 전후의 집단크기는 Rogers and Harpending (1992)의 방법을 따라 계산하였다.

꼼치의 지역간 계수 및 계측 결과는 Table 1과 같고, 각 지역 집단의 계수형질의 빈도 분포는 Table 2에 나타내었다. 지역 집단간 계수형질 차이 유무를 파악하기 위한 Kruskal-Wallis test 결과, 뒷지느러미 수에서만 유의한 차이를 보였다(P<0.05) (Table 3).

주성분분석(PCA)을 실시한 결과, 1축 고유치는 14.227 (기여율: 71.137%), 2축은 1.344 (기여율: 6.722%), 3축은 0.851 (기여율: 4.253%)이었다(Table 4). 총 20개의 계측형질 중 전장(고유치: 0.936), 체장(0.933) 및 등지느러미 앞길이(0.9358)가 1축에서 가장 큰 기여율을 나타내었고, 흡반과 항문사이 길이(0.772)가 2축에서 가장 큰 기여율을 나타내었다(Table 4). PCA 결과 1축으로는 3개 지역이 대부분 중첩되어 보였지만 2축으로는 보령이 나머지 두지역과 약간 구분되었다(Fig. 2).

정준판별분석(CDA)을 실시한 결과, 각 지역 집단은 매우 높은 판별적중율(보령과 포항은 96.6%, 진해는 100%)을 보였다. 1축의 고유치는 3.454 (기여율: 59.1%), 2축은 2.388 (기여율: 40.9%)이었고(Table 5), 1축으로는 주둥이에서 흡반까지 길이(0.159), 2축으로는 흡반에서 항문까지 길이(0.243)가 3개 지역을 분리하는데 가장 기여한 형질이었다(Table 5). CDA 결과 1, 2축으로 3개 지역이 명확하게 구분되었고, 1축으로는 진해가 다른 2개 지역과 잘 분리되었으며, 2축으로는 보령과 포항이 잘 분리되었다(Fig. 3).

유전다양성

미토콘드리아 DNA cytb 762 bp가 보령(25개체), 진해(25) 및 포항(23)의 3개 지역 집단에서 확보되었으며, 이들 염기서열은 모두 NCBI에 등록하였다(KU362804-KU362876). Cyt b 분석 결과 6개의 위치에서 염기치환이 확인되었으나, 삽입이나 결실은 없었다. 각 지역 집단내 유전거리는 0.0-0.3%, 지역 집단간 유전거리는 0.0-0.4%로 집단간 차이는 집단내 차이와 같거나 약간 더 컸다.

한국산 꼼치의 유전자형은 7개로, 그 중 LTcytb1이 55개체에서 나타나 가장 우점하였고(75.34%), 이 유전자형은 3지역(보령, 진해, 포항)에서 모두 발견되었다. 그 다음으로 LTcytb2가 9개체에서 나타나 우점하였는데(25.3%), 이 유전자형은 진해와 포항개체에서만 발견되었다(Table 6). 지역별 유전자형다양성(h) 및 염기다양성(π)은 Table 7에 나타내었다. 유전자형다양성은 진해가 0.5767로 가장 높았고, 포항이 0.5099로 그 다음으로 높았으며, 보령이 0.0800으로 가장 낮았다. 염기다양성은 세 지역 모두 낮았는데, 진해와 포항이 0.0076의 값을 나타냈고, 보령이 0.000229로 낮았다. 지역 집단간 유전적 분화 정도를 나타내는 고정지수(F_ST)는 보령과 진해가 0.17 (P<0.05), 보령과 포항이 0.07 (P<0.05)로 유의미한 차이를 나타냈다. 반면, 진해와 포항은 0.03 (P>0.05)으로 유의미한 차이를 보이지 않았다(Table 8).

집단구조

최소근접네트워크(MSN)는 1개의 거대 유전자형(LTcytb1) 주위로 5개의 작은 유전자형이 1개의 염기 차이로 연결된 별모양을 나타내었다(Fig. 4). 세 지역(보령, 진해, 포항)은 모두 Tcytb1이란 거대 유전자형을 공유하고 있으며, 진해와 포항만 LTcytb2의 단일 유전자형을 가졌다. 또한, 진해는 포항과 LTcytb3의 유전자형도 공유하였다. 나머지 유전자형(LTcytb4-Ltcytb7)은 각 지역마다 독특하게 가지는 것으로 나타났다. 꼼치 73개체의 염기서열과 외집단인 Liparis liparis (EU492094)와의 염기서열 분석 결과, 꼼치의 종내 유전거리는 0.0-0.3%였고, L. liparis와의 유전거리는 12.3-12.7%의 범위를 나타내었다.

집단변천사

집단의 팽창 연대를 추정할 수 있는 매개변수(τ)는 보령이 3.000으로 최대를 보였고, 다음으로 진해가 높은 값(0.828)을 보였다(Table 9). 포항은 가장 낮은 값(0.713)을 보였다. 팽창 전 돌연변이 매개변수(θ₀)는 세 지역 모두 0 이었지만, 팽창 후 돌연변이 매개변수(θ₁)는 진해와 포항이 99999.0로 최대를 나타내었고, 보령만 0.091로 매우 낮았다(Table 9). Tajima’s D test는 세 지역 모두 음의 값을 보였으나, 집단의 급격한 팽창 모델과 달랐다(P>0.05). Fu’s Fs test는 세 지역 모두 음의 값을 보였으며, 그 중 포항만 유의미한 값을 나타냈다(P<0.05). Mismatch distribution 분석 결과에서는 3개의 지역 모두 급격한 팽창 모델을 잘 반영해 주는 것으로 나타났으며(보령: P=0.25, 진해: P=0.11, 포항: P=0.68), 세 집단은 mismatch distribution에서 1개의 peak를 보였다(Fig. 5).

꼼치의 계측형질 20개로 주성분분석을 실시한 결과, 세 지역(보령, 진해, 포항)이 약간중첩되는 경향을 보였으나, 제 2축을 기준으로 보령이 나머지 두 지역(진해, 포항)보다 약간 높은 값을 보이면서 구분되었다(Fig. 2). 반면, 정준판별분석은 세 지역을 명확히 구분지어 주었다(Fig. 3). 이처럼 종내 형태 특성에 따라 지역 집단이 구분되며, 이를 자원관리의 정책 수단으로 제시되기도 한다. 예를 들면, Turan et al. (2006)은 터키 주변해역에 서식하는 Pomatomus saltatrix의 다변량 형태 분석을 통해 세 지역 집단으로 구분할 것을 주장한 바 있다. 또한, 한국 및 일본산 국수뱅어(Salanx ariakensis)를 대상으로 정준판별분석 결과 두 지역집단이 잘 구분되어(Kim et al., 2006b), 서로 다른 관리단위(management unit)로 간주되어야 한다고 주장되었다(Kim et al., 2007). 비슷한 시기에 Murta et al. (2008)은 대서양과 지중해산 전갱이류(Trachurus trachurus)의 형태를 분석한 결과 6개의 관리단위로 구분해야 한다고 주장하였다. 최근 한국산 전어(Konosirus punctatus) 여러 지역 개체를 대상으로 정준판별분석 결과 남해 중부(여수)를 중심으로 2개의 morphotype으로 구분되나, 남해 동부(부산)에서 2개 morphotype 사이 교잡 가능성이 제기되어 관리단위를 명확히 하기 위해 추가 연구가 필요하다고 보고한 바 있다(Myoung and Kim, 2016).

지역 집단이 어느 정도 잘 구분되는가는 판별적중률을 통해 알 수 있다. 한국산 꼼치의 경우, 판별적중률은 진해가 100%로 최대였고, 보령과 포항은 96.6%로 1개체를 제외하면 채집된 지역의 형태 특성을 잘 반영해 주었다. Silva (2003)는 판별분석에서 집단별 판별적중률이 85% 이상이 되면 형태적으로 구분이 가능하다고 주장하여, Silva (2003)의 기준에 의하면 한국산 꼼치는 세 지역 집단으로 구분 가능하며, 주둥이에서 흡반까지의 길이가 중요한 기여 형질이다.

형태 형질은 유전 형질의 영향을 가장 많이 받지만 그 외 성장 과정에서 환경 요인의 영향도 많이 받는다. 특히 난 또는 자치어 시기에 환경 요인의 영향을 지배적으로 받는다(Lindsey, 1988; Turan, 2000). 예로, Scott (1972)는 까나리속(Ammodytes) 어류의 경우 수온과 같은 환경 요인이 계수 형질이나 성장율에 큰 영향을 미칠 수 있다고 주장하였다. 한편 다른 학자들은 계수 형질과 달리 계측 형질이 다양한 환경 변화에 민감하게 반응한다고 주장하였다(Lindsey, 1988; Swain et al., 2001). 한국산 꼼치의 경우 정준판별분석 결과 세 지역 집단으로 명확히 구분되었는데, 이는 세 지역의 해양 환경 특성이 다른 데서 비롯된 것이다. 더욱이 꼼치는 점착성 알을 낳고, 흡반을 가지는 저서 어종이므로(Huh, 1997), 소규모 이동을 하거나 정착성 어종으로 판단된다. 따라서, 꼼치의 형태 특성이 지역별 차이를 보인 것은 해양 환경 차이와 정착성이란 생물 고유의 특성에서 비롯된 것으로 사료된다.

계측형질을 이용한 다변량분석은 중요한 수산자원의 관리단위와 관리범위를 정하는데 매우 중요한 자료를 제공해 준다. 본 연구에서 꼼치는 정준판별분석으로 3개 지역집단으로 명확하게 구분되었으므로 3개의 관리단위를 설정할 필요가 있다고 생각되지만, 분자분석 결과에서는 2개 집단(보령 vs. 진해와 포항)으로 구분되었기 때문에 향후 보다 민감도가 높은 microsatellite DNA에 의한 추가 분석이 필요하다.

세 지역 집단은 대부분 542번째 염기서열에서 G를 가지고 있었지만, 진해와 포항 9개체에서는 A를 가져 보령과 의미 있는 차이를 보였다. 뿐만 아니라 362번째 염기서열에서도 대부분 T를 가졌으나 진해와 포항 3개체에서만 C를 가져 역시 보령과 의미 있는 차이를 보였다. 이러한 관계는 F_ST를 통해서도 알 수 있는데, 보령과 진해는 0.17, 보령과 포항은 0.07로 유의미한 차이 (P<0.05)를 보인 반면, 진해와 포항은 차이가 없었다(Table 8). Wright (1965)는 두 지역 집단간 F_ST가 0.05-015 범위에 있으면 약간 차이가 있고, 0.15-0.25 범위에 있으면 상당한 차이가 있다고 보아, 본 연구에서 보령과 진해(0.17)가 보령과 포항(0.07)에 비해 더 오랜 진화 역사를 가진 것으로 추정된다. 그러나 공유 유전자형은 보령과 포항이 16개, 보령과 진해가 15개로 1개 밖에 차이가 나지 않는 점을 고려하면 F_ST가 보여주는 것과 달리 보령-포항이 보령-진해와 유사한 시기에 분화가 이루어졌을 것으로 사료된다. 이처럼 유전적 분화가 진행된 경우 진화적으로 의미 있는 단위(Evolutionary significant unit: ESU) 또는 관리단위(Management unit: MU)의 두 관점에서 논의가 가능하다(Crandall et al., 2000). 예들 들어, Kim et al. (2006a)은 한국 및 일본산 국수뱅어의 cytb 염기서열에 의한 F_ST는 0.43을 보였으나 두 지역집단이 단계통을 형성하지 않는 때문에 결론을 내리지 못했다. 그러나 이후 Kim et al. (2007)은 동일개체를 대상으로 AFLP 추가 연구를 통해 두 지역집단이 잘 분리됨을 밝혔으나 지역집단을 정의내릴 수 있는 고유 형질이 없다는 점에서 서로 다른 관리단위(MU)로 보아야 한다고 주장하였다. 이와 유사하게, 꼼치의 경우도 보령을 진해 및 포항과는 별개의 독립된 관리 단위로 보는 것이 타당해 보인다.

집단의 진화 역사는 유전자형 다양도와 염기 다양도를 이용하여 추정할 수 있다. 유전자형 다양도와 염기 다양도가 낮은 집단은 최근 병목현상이 일어난 집단으로 볼 수 있다(Grant and Bowen, 1998). 본 연구에서 보령이 유전자형과 염기 다양도에서 낮은 값을 나타냈으며, 진해와 포항은 높은 유전자형 다양도와 낮은 염기 다양도를 나타내어 병목현상이 일어난 뒤 급격한 팽창 사건이 있었음을 알 수 있다(Grant and Bowen, 1998).

일반적으로 미토콘드리아 DNA cytb 영역의 돌연변이율은 100만년당 2%로 알려져 있다(Liu et al., 2006b). 이에 따라 꼼치의 개체군 팽창 연대를 계산해 보면 18,713-78,740년전에 팽창이 일어난 것으로 추정된다. 한반도에 마지막 빙하기는 플라이스토세에 있었고, 이 시기에 서해지역의 해수면은 현재보다 130 m 낮아 육지 형태를 보였다(Lambeck et al., 2002; Kitamura and Kimoto, 2006; Khim et al., 2007). 수심이 얕은 한반도와 일본 사이는 해수면의 하강으로 육지로 연결되었고 그 결과 동해와 남해는 거의 분리되다시피 하였다. 까나리와 전어를 대상으로 한 집단유전 연구에서 빙하기의 지형적 격리가 어류와 같은 해양생물의 지리적 분화에 큰 영향을 줄 수 있음을 알 수 있다(Kim et al., 2006a; Myoung and Kim, 2014). 따라서, 꼼치는 과거 빙하기에 바다가 분리됨에 따라 환경적 차이로 인한 이소적 분화를 겪었을 것이다. 또한, 꼼치가 침성부착난을 낳고, 거의 이동이 없는 어종임을 고려한다면 이러한 지역간 차이는 시간이 흐를수록 더욱 커질 것으로 전망된다.