Genome-wide Single Nucleotide Polymorphism-based Assay for Phylogenetic Relationship of the Flammulina velutipes

공시균주

공시균주에는 농촌진흥청 버섯과에서 수집 보존 중인 유 전자원 중 팽이버섯속(Flammulina)에 속하면서 종이 다른 계통, 팽이버섯(F. velutipes) 종과 같은 종이지만 variety가 다른 계통, 그리고 동일 종내 수집지역과 특성이 상이한 대표적 계통들을 포함하였다(Table 1).

Table 1. List of Flammulina strains used in this study

ASI, Agricultural Sciences Institute; B, brown; LB, light brown; W, white.

Genomic DNA 추출

시험에 사용된 25균주의 Genomic DNA를 분리하기 위 하여 potato dextrose agar (PDA)에 균사를 접종하여 25ºC 의 incubator에서 7일간 배양한 후 균사를 수확하여 동결 건조시킨 후 실험에 사용하였다.

DNA는 1.5 mL 튜브에 샘플의 양을 tube의 0.05~0.1 mL 정도의 위치로 넣고 뚜껑을 닫고 액체 질소에 담아둔다. 샘 플을 마쇄하기 위해서 급랭한 샘플들이 녹기 전에 액체질 소에서 재빨리 꺼내어 멸균된 막대나 bead를 사용했다. 900 μL 정도 extraction buffer인 cetyl trimethyl ammonium bromide (CTAB) buffer를 첨가한 후 샘플들을 다시 vortexing 작업을 했다. Water bath에서 65ºC로 1시간 정도 처 리한 후 450 μL의 chloroform을 첨가한 후 실온에서 반응 시키기 위해서 약 5분 정도 tapping 해준 다음 12,000 rpm, 20ºC에서 10~15분 원심분리하였다. 상등액을 600 μL 채취 하여 새로운 튜브에 옮긴 후 5 μL RNase-A (10 mg/mL)를 넣고 37ºC에서 30분 처리하였다. 처리 후에 600 μL Isopropanol을 첨가하여, 약 1분 정도 tapping 한 후 12,000 rpm, 20ºC, 10~15분 원심분리한 후 튜브 아래에 pellet은 남기고 상층액만 버렸다. 500 μL의 70% 에탄올을 첨가하여 12,000 rpm, 20ºC, 10분을 원심분리한 후 상층액을 버리고 다시 에 탄올 작업으로 두 번씩 씻은 후, 상온에서 약 30분 건조시 켜 3차 멸균수를 DNA pellet 양에 따라 조절해서 50~100 μL 정도 넣고 pellet이 다 녹도록 한 후 4ºC에 보관했다.

팽이버섯 resequencing data를 이용한 genome-wide SNP 분석

Fig. 1. Flowchart for single nucleotide polymorphism (SNP) detection and selection within Flammulina velutipes strains.

Fig. 1은 팽이버섯 균주간의 계통 분석을 수행하기 위한 분석 절차를 나타낸 것이다. Short read의 정제과정(Solexa QA package)을 통하여 분석 가능한 품질의 데이터를 확보 하고[8], read mapping BWA (Burrows Wheeler Aligner) 을 통해 각 샘플의 consensus sequence를 추출하였다. SAM (Sequence Alignmint/Map) 포맷의 파일을 읽고, 쓰 고, 편집하고, 볼 수 있는 program인 SNP detection (SAM tools)을 사용하였다. 그 과정과 variation을 통해 filter된 SNP을 이용하여 팽이버섯 25개 균주의 시퀀스를 multiple sequence alignment (clustalW: 두 염기서열을 비교하는 EBI의 tool) 수행하고, 이후 phylogenetic tree (MEGA6)을 통해 샘플간의 계통적 유연관계를 분석하였다[9].

Sequence 전처리

차세대 염기서열 분석을 통해 생산된 팽이버섯의 sequencing data는 농촌진흥청 버섯과에서 운영하는 공개 데이 터베이스(http://112.220.192.2/mushroom/)에서 수집하였 다.

팽이버섯 25균주는 차세대 염기서열 분석 (Hi-seq 2000; Theragen Etex, Suwon, Korea) 데이터로 구성되었으며, 표 준유전체인 KACC42780 (4019-20)과 비교 분석하여 genome- wide SNP를 추출하였다[7]. 서열의 품질(sequence quality)을 측정하고, 기준 이상 품 질의 서열로 filtering 하는 Illumina Hiseq platform/Solexa QA (v.1.13) package의 DynamicTrim과 LengthSort 프로그 램을 이용하여 전처리 과정을 진행하였다[8]. Dynamic Trim은 phred score에 따라 short read의 양쪽 끝의 bad quality base를 잘라내고 양질의 cleaned read로 정제하는 과정을 수행하였다. LengthSort는 DynamicTrim에서 너무 많은 base가 잘린 read를 제거하는 과정을 수행하였다. Dynamic Trim의 phred score ≥ 20을, LengthSort 과정은 short read length ≥ 25bp 작업을 수행하였다[10, 11].

표준유전체의 alignment

BWA는 short read들로 표준유전체(reference genome)에 mapping하는 것이 게놈을 위한 가장 기본적인 Bioinformatics의 분석이라고 할 수 있다. 그 전처리 과정을 통과한 cleaned reads를 BWA (0.6.0-r104)을 이용하여 표준유전체 로의 alignment를 수행하였고, mapping은 표준유전체와 sequence한 샘플간의 raw SNP(In/Del)를 detection하기 위 한 선행 과정으로서 Binary Alignment/Map(BAM) format 의 파일을 생성하며, 적용한 옵션은 다음과 같다. Maximum number of gap extensions (-e) = 50, seed length (-l) = 30, maximum differences in the seed (-k) = 1, number of threads (-t) = 16, mismatch penalty (-M) = 6, gap open penalty (-O) = 15, gap extension penalty (-E) = 8로 설정 하였다. 유전체 재분석 데이터는 표준유전체와의 alignment 후에, 각각 1개의 파일로 합친 후, SNP 분석을 진행 하였다[12].

Raw SNP detection 및 consensus sequence 추출

BAM과 SAM은 모두 sequence를 저장한 정보를 가지고 있으며, BAM 파일은 SAM 파일과 동일하지만, reference sequence names, length들이 헤더에 포함되어 있다. 또한 프로그램으로 원하는 sequence를 빨리 찾기 위한 indexing 기법을 활용하였다. Clean reads를 표준유전체에 mapping 하여 생성된 BAM format의 파일을 SAMtools (v0.0.16)를 이용하여 raw SNP (IN/Del)을 detection하고, consensus sequence를 추출하였다. SAM tools 사용 옵션 이외에는 기 본값을 사용하였다. Minimum mapping quality for SNPs (-Q) = 30, minimum mapping quality for SNPs (-q) = 15, minimum read depth (-d) = 3, minimum read depth (-D) = 100, min indel score for nearby SNP filtering (-G) = 30, SNP within INT bp around a gap to be filtered (-W) = 15, window size for filtering dense SNPs (-W) = 15로 설정하였다. 추출된 SNP 결과에서 정확도 높은 SNP를 선 발하기 위해 mapping quality 값을 기본값 25보다 높은 30 으로 적용하고, raw reads를 갖는 SNP로 filtering 하였다 [13].

SNP matrix의 작성

팽이버섯 분석 간의 SNP 비교분석을 수행하기 위해 샘 플간 통합 SNP matrix를 작성하였다. 각 샘플을 표준유전 체와 비교하여 얻은 raw SNP position을 후보로 하여 합집 합의 list를 구축하고, 이 때 빈 영역(non-SNP loci)은 sample의 consensus sequence로부터 채워 넣는 filling 과정을 거쳐 matrix를 작성하였다[14]. 이후 샘플간의 SNP 비교를 통해 mis-calling된 SNP를 필터하고, SNP 유형 구분 기준 은 SNP read depth 정확도 높은 것이 90% 이상일 때 homozygous이며, SNP read depth 50%일 때 heterozygous라 고 하며, homozygous와 heterozygous 유형으로 구분할 수 없는 경우는 기타로 표현하였다(Table 3).

유연관계 분석

공시균주의 유전체 재분석 데이터를 확보한 후, National Center for Biotechnology Information (NCBI)의 Basic Local Alignment Search Tool (BLAST)을 이용하여 DNA 데 이터베이스와 유사한 염기서열을 비교하였다. Phylogenetic analysis (계통발생분석)은 MEGA vesion 6 program의 neighbor- joining, maximum likelihood, minimum evolution 방법에 의해 분석하였고 sequence distance는 1,000 replication의 bootstrap value에 의해 계산하였다[9, 15, 16]. Neighbor-joining는 근접군 연결 방법으로 가장 속도가 빠 르며 격리된 근접군 쌍 사이와 나머지 가지를 연속해서 삽 입하여 분석 이후 가장 가까운 근접군 쌍을 통합하였다.

Genome-wide SNPs 분석

팽이버섯 25계통의 유전체 재분석 데이터를 표준유전체 (KACC42780)와 비교하여 genome-wide SNP를 선발하였 다. Total reads은 ASI 4188의 51,018,766개로 가장 많았으 며, 단핵균주인 mono19에서 3,331,674개로 가장 낮은 것을 볼 수 있었다. Total reads중 mapped reads는 ASI 4178에 서 18,503,265개 (68.62%)로 보였으며, ASI 4188에서는 9,416,665개 (18.46%)로 낮은 것을 볼 수 있다. Readsmapping region은 ASI 4175에서 가장 많은 32,294,244개 (90.61%)로 보였고, ASI 4190에서 22,653,709개 (63.56%) 로 가장 낮은 것을 볼 수 있었다. 각 균주의 read mapping 율을 통해 팽이버섯 자원내에서도 변이가 아주 크다는 것 을 알 수 있었다. 즉 표준유전체와 변이가 큰 균주일수록 mapping율이 낮고 N count (bp)가 높아지는 것을 알 수 있었다(Table 2).

이는 품종에 따른 표준유전체의 서열 유사 정도에 따라 나타난 polymorphism의 차이로 보인다[17, 18]. 각 균주의 SNP를 homozygous SNP와 heterozygous SNP로 추출하였 다. ASI 4188의 homozygous SNP는 473,518개이고 heterozygous SNP는 3,124개, 그리고 구분이 안 되는 other가 10,513개로 총 476,642개의 SNP가 선발되었다. 확보된 genome-wide SNPs의 분포를 보면 homozygous SNP가 가 장 많은 균주는 ASI 4188로 473,518개로 검출되었으며, 반 면 ASI 4175는 가장 적은 89,009개를 보였다. 또한 heterozygous SNP에서는 가장 많은 균주는 ASI 4178로 86,287개 이었으며, 가장 적은 균주는 단핵균주인 mono19로 1,945 개를 볼 수 있었다(Table 3).

Table 2. Summary of alignment against the Flammulina velutipes reference genome

ASI, Agricultural Sciences Institute.

aReads-mapping region: consensus sequence of base = N (un-mapping region) deprived area. (Reads-mapping region / total length of reference genome) × 100

Table 3. Summary of SNPs detection obtained from Flammulina strains

SNP, single nucleotide polymorphism; ASI, Agricultural Sciences Institute.

공시군주의 대부분이 이핵균주이므로 이형 핵의 존재로 hetero SNP가 검출되는 것은 당연하겠으나 이핵균주들 간 의 편차가 가장 많은 것이 86,287개로 가장 적은 것이 2,096개로 보였으며 그 이유를 이형핵만으로는 설명하기 어렵다. 또한 단핵균주인 momo19는 핵이 하나이므로 hetero type이 나오지 않을 것으로 기대했으나 1,945개가 검 출되었다. 그 이유는 추후에 더 깊은 연구가 필요하겠으나 균류에서 흔히 나타나는 염색체의 polymorphism과 관련이 있을 것으로 생각한다[19].

Genome-wide SNP에 의한 유연관계 분석

그룹 간의 근연적 거리를 살펴보기 위해, maximum likelihood 방법으로 보았을 때 Group I은 F. velutipes var. 계 통인 ASI 4062, 4148, 4195이며, Group II는 ASI 4188 (F.elastica), ASI 4190 (F. fennae), ASI 4194 (F. rossica) 균주 가 포함되었으며, 서로 다른 학명을 가졌으므로 상당히 떨 어진 모습을 나타냈다. 그 외 F. velutipes 19개 계통은 같은 그룹으로 나타났다. 유연관계 분석은 radiation에 정확성을 알기 위해 3번 반복한 결과이며, 학명에 따라 그룹 간에는 같은 결과가 확인되었다(Fig. 2). 이와 같은 유연관계 분석 은 ITS 분석에 의해서도 손쉽게 얻어질 수 있으나 [6, 20], genome wide의 SNP정보를 이용해도 동일한 결과를 보였 다. 이는 ITS 염기서열의 변이가 전체 게놈 내에서 일어난 진화적 변화 과정을 잘 반영하는 것으로 판단할 수 있다.

백색집단과 갈색집단의 유연관계 비교

백색집단과 갈색집단의 유연관계 비교 Neighbor-joining, maximum likelihood, minimum evolution을 2번 반복으로 유연관계 분석을 하였다. 각 그룹별 로 백색팽이군 공통 SNP와 갈색팽이군 공통 SNP를 비교 하여 자실체 색깔을 반영할 수 있는 SNP를 선발하고자 하 였다.

팽이버섯의 계통간 비교에서 neighbor-joining 결과는 백 색 group에는 한국 팽이버섯과 일본 팽이버섯이 함께 포함 되었으며, 그 중 ASI 4175, 4103, 4028는 갈색임에도 불구 하고 백색 팽이 group에 포함되었다. 갈색 group에서는 한 국과 일본, 중국의 팽이가 한 그룹을 형성하였다. Neighborjoining, maximum likelihood, minimum evolution의 방법 을 이용하여 백색 팽이버섯과 갈색 팽이버섯 계통간 색깔 별 구분을 유연관계로 분석하고자 시도하였으나 색깔에 따 른 group은 이루어지지 않았으며, 방법에 따라 grouping은 비슷한 경향을 보였다. 이와 같은 결과는 분석에 사용된 공시균주의 선발에서 제한된 수에 기인했을 수도 있으며, genome-wide의 SNP를 사용함으로써 색깔 이외의 형질에 대한 변인 통제가 제대로 되지 않아 색깔 이외의 유전요인 에 관련된 SNP가 분석에 동등하게 포함된 것이 원인으로 생각해 볼 수 있었다. 그러므로 phylogenetic tree를 만드는 방법에 따라 백색 팽이버섯과 갈색 팽이버섯 group간 구분 할 수는 없으나, Group I과 Group II로 분석되었던 Fig. 2 의 결과와는 일치하였다(Fig. 3).

Fig. 2. Phylogenetic tree based on SNPs of Flammulina strains representing their phylogenetic relationship. SNP, single nucleotide polymorphism.

Fig. 3. Phylogenetic trees of SNPs among Flammulina velutipes strains inferred by neighbor-joining, maximum likelihood, and minimum evolution methods. SNP, single nucleotide polymorphism.

Fig. 4에서는 백색 팽이버섯 중 한국 품종과 일본 품종간 을 neighbor-joining으로 유연관계를 분석하였다. 한국 백 색 균주인 ASI 4210, 4166, 4178과 일본 백색 균주인 ASI 4209, 4167이 하나의 그룹으로 형성되었다. Neighbor-Joining는 2번 반복으로 같은 결과가 나타났으며, 갈색 균주들 과 백색균주가 3번째 노드에서 구분되었다. Phylogenetic tree 상에서 한국 백색 균주와 일본 백색 균주의 유전적 상 동성이 매우 높음을 확인할 수 있었다. 비록 사용된 균주 수는 많지 않으나 이와 같은 결과를 이용하여 버섯을 발생 시키지 않고 버섯 색깔을 예측할 수 있는 분자 마커를 개발 할 수 있을 것으로 기대되며 더 나아가 상업적으로 이용되 고 있는 백색 품종들간에 지리적 기원을 밝힘으로써 종자 주권 확립에 기여할 수 있을 것으로 기대한다.

Fig. 4. Phylogenetic relationship of Korean and Japanese strains producing white fruit body.