밀(Triticum aestivum)은 식량 작물 중 두 번째로 많이 생산되며 옥수수, 벼와 함께 세계 3대 식량 작물이다. 보통밀, 빵밀로도 불리며 밀 속(Triticum)에서도 가장 많이 재배된다(Kang et al. 2010). 밀은 추파 또는 춘파, 적립계 또는 백립계 등 유전적 특성을 가진 다양한 품종들이 가진 장점과 높은 수량성 및 품질을 보유하고 있다(Lee et al. 2002). 이처럼 오늘날 요구되는 식량작물로서의 광범위한 기능적 요소를 모두 충족하는 장점으로 인해 오래 전부터 인류 생존과 관련하여 필수적인 작물 중 하나로 널리 소비되고 있다(Choi et al. 2018).

밀은 타 작물과 다르게 대부분의 경우 가공의 단계를 거쳐야 식품으로서의 가치가 높아지기 때문에 가공적성(end use quality)은 밀에서의 주요한 요인으로 여겨지며 특히 생물적, 비생물적 스트레스 피해는 수확량 뿐만 아니라 식품으로서 밀의 가치를 떨어뜨리는 주요 방해 요소가 될 수 있다(Kang et al. 2008). 최근 이상 기후로 인한 재배 환경 변화와 강수량의 변화로 밀의 피해는 더 심화되고 있는 추세이며, 태풍이나 기습적 폭우 등의 환경적 요인에 의해 작물의 도복이나 과습 피해의 발생 조건이 빈번해지면서 식량 작물의 대표적 수분장해인 수발아의 발생과 그로 인한 피해는 점차 확대되고 있다(Shin et al. 2013). 수발아(pre-harvest sprouting; PHS)는 종실이 출수기 이후 수확을 앞둔 등숙 기간에 잦은 강수로 인한 과수분 조건과 발아억제물질의 용해 및 분해되어 종자의 휴면타파와 함께 이삭에 달린 채로 발아가 진행되는 것을 말하며(Schopfer et al. 1979, Nielsen et al. 1984, Biddulph et al. 2008, Kang et al. 2014) 수발아가 시작되면 배유의 α-amylase의 활성 증가와 함께 탄수화물이 당으로 분해되면서 밀 종실의 품질을 극도로 저하 시키는 재해로 전세계적으로 광범위한 피해를 일으키는 대표적인 비생물적 스트레스라고 할 수 있다(Andreoli et al. 2006, Farashah et al. 2011, Kang et al. 2014). 특히 국내 밀의 육종 방향은 균일하고 빠른 발아가 주된 목적 형질이었으며 주요 장려품종들은 휴면(seed dormancy)에 약하여 수발아에 취약한 점이 있어 최근의 기후 문제와 더불어 보완이 필요한 것으로 판단된다.

최근 작물 분자 육종 연구에서 생물정보학의 활용은 트렌드를 넘어 필수적인 연구 요소 중 하나이다. 특히 RNA 분석 연구에서 마이크로어레이(microarray) 기술과 더불어 NGS 기술을 통해 정량화된 mRNA를 시퀀싱하는 RNA-seq 기술은 후보 유전자의 주석화(annotation) 및 각 유전자의 정확한 발현량 분석이 가능하다. 특히, 알려져 있는 유전자의 범위를 벗어난 분석도 가능하기때문에 밀과 같이 유전체 연구가 완벽하지 않은 작물에서도 유전자 선발 및 육종 소재 개발에 많은 도움이 되고 있다(Duan et al. 2012, Pearce et al. 2015).

본 연구는 수발아에 감수성인 국내 장려 품종 ‘금강’과 금강을 모본으로 돌연변이 육종을 통해 수발아 내성을 가진 ‘전주 377호’를 대상으로 수발아 조건 실험을 실시하였으며 이를 대상으로 대량전사체 분석을 하였다. 수발아 조건에서의 발현 유전자 기작(mechanism) 및 기능 분석을 실시하였으며 수발아 관련 후보 DEGs (differentially expressed genes) 발굴을 통한 발현 양상을 파악하여 밀의 휴면과 발아에 대한 이해를 통해 향후 저항성 밀 품종 육성에 이용하고자 위함이다.

식물 재료 및 수발아 처리

본 실험에서는 국립 식량과학원에서 분양 받은 국내 장려 품종 ‘금강’과 수발아 내성 품종인 ‘우리’ 그리고 돌연변이 육종법으로 개발된 금강 돌연변이 계통 ‘전주 377호’를 대상으로 실시하였다. 각 계통은 4°C 에서 42일 간 춘화처리(vernalization) 후 포트에 1개체 씩 파종하였고 각 개체 당 첫 이삭(main tiller)을 기준으로 출수 후 40일(Days after fertilization)까지 생육하였다.

‘전주 377호’는 국립식량과학원에서 2007년 국내 장려 품종인 ‘금강’밀을 방사선 처리하여 개발하였다. 돌연변이원은 Co⁶⁰이며 100 Gy의 선량을 이용하여 M₁ 변이체를 만든 후 bulk 재배로 M₄까지 세대를 진전시켰으며 이후 수발아 저항성 선발을 통하여 계통 진행으로 F₈ 세대까지 진전하여 개발한 돌연변이 밀 계통이다(RDA 2017).

수발아 실험은 2021년 6월 16일부터 7일간 공주대 온실에서 실시되었다. 800×1100 cm 규격의 배드에 원예용 상토를 살포하여 평균 2 cm 높이의 인공적인 노지 환경을 조성하였고 직경 5.5 cm의 분수 호스를 활용하여 일평균 12시간 급수하여 강우 조건을 유지하였다. 샘플은 수발아 처리가 1일된 ‘금강’(K1), ‘전주377호’(KM1)의 식물체를 포함한 이삭을 대조군으로 사용하였으며, 7일간 수발아 처리한 ‘금강’(K7), ‘전주377호’(KM7) 이삭을 실험군으로 정하였다. 시료는 수발아 처리 전 후에 품종 당 독립적인 3개체씩 수확하였으며 액화질소를 활용해 급속 냉각하여 다음 실험 전까지 -70°C 초저온냉동고에서 보관하였다.

RNA-seq 라이브러리 구축 및 시퀀싱

Total RNA 추출은 GeneAll 사의 RibospinTM Seed/Fruit (GeneAll® Biotechnology, Seoul, Korea)를 이용해 수행하였으며 제조사의 프로토콜에 따라 추출하였다. cDNA 합성은 Power cDNA Synthesis Kit (iNtRON Biotechnology, Seoul, Korea)를 사용하여 수행하였다. RNA-seq 분석을 위한 라이브러리는 TruSeq™ RNA library prep kit (Illumina, CA, USA)를 사용하여 구축되었고 100-bp paired end sequencing은 Illumina HiSeqTM- 2500 platform (Illumina, CA, USA)에서 수행되었다.

RNA-seq 데이터 전처리 (read mapping and filtering)

De novo assembly는 velvet, oases 어셈블러(Zerbino et al. 2008, Schulz et al. 2012)를 이용하여 진행되었으며, 이후 Short reads의 전처리 과정을 위하여 SolexaQA (v.1.13) package (Cox et al. 2010)의 DynamicTrim과 LengthSort 프로그램을 사용하였다. Assembly 서열에 bowtie2 (v2.-1.0) software (Langmead et al. 2012)를 활용하여 read mapping을 하였으며, DESeq library (Anders et al. 2010)를 이용하여 Normalization을 실시하였다.

RNA-seq 데이터 및 DEGs 분석

Assembled transcript의 기능적 분류 및 파악을 위해 유전자 서열을 활용하여 KOG (eukaryotic orthologous groups) 및 GO (gene ontology) 분석을 수행하였고(Ashburner et al. 2000), 근연종인 벼(Oryza sativa)의 orthlog 유전자 서열을 통해 KEGG (Kyoto encyclopedia of genes and genomes) 분석을 수행하였다.

샘플 간 유의하게 발현되는 유전자(DEGs; differentially expressed genes)선발은 DESeq library에 포함된 방법론(Anders et al. 2010)을 사용하여 분석하였다. 각 유전자에 mapping된 발현 값이 샘플 간에 비교에서 2배 이상 발현의 차이를 구분하는 2 fold change방법과 adjust E-value이 0.01 이하를 만족하는 binomial test 방법을 사용하였다. log2 Fold Change의 값이 1보다 크면 up-regulation, -1보다 작으면 down-regulation 되었다고 판단하였다. 유의하게 발현된 유전자 정보를 이용하여 유전자 발현 패턴을 확인하기 위하여 R의 amap (Lucas 2018)과 gplot (Warners et al. 2015) library를 이용하여 hierarchical clustering 분석을 실시하였으며 피어슨계수(pearson’s correlation)를 이용하여 발현 유사 정도를 분석하였다.

DEGs 발현 검정

Quantitative Real-time PCR은 QIAGEN의 Rotor-Gene Q (QIAGEN, Hilden, Germany) Real Time-PCR machine을 이용하였으며 Rotor-Gene SYBR Green PCR kit (QIAGEN, Hilden, Germany)를 사용하여 실시하였다. 선발된 DEGs 특이 프라이머는 Primer3 소프트웨어를 사용하여 설계하였으며 qRT-PCR 조건은 제조사에서 제시한 thermal cycling 조건을 따랐다. Housekeeping gene으로는 TaActin을 활용하여 유전자 발현을 정규화 시켰으며 각 sample에 대하여 생물학적 반복(biological repeat), 기술적 반복(technical repeat)을 각각 3회 실시하고 상대 발현 수준은 2^-ΔΔCt 방법(Livak et al. 2001)을 통해 정규화(normalization)하여 분석하였다.

수발아 처리 및 발아에 따른 표현형 측정

RNA-seq 진행을 위하여 국내 장려품종이자 수발아 감수성인 ‘금강’과 수발아 저항성으로 알려진 ‘전주377호’를 대상으로 수발아 실험을 실시하였다. 실험은 온실에서 진행되었으나 노지 수발아 환경을 최대한 구현하기 위하여 Nakamura et al. (2011)의 실험 방법을 따랐으며 생육 온도는 28°C/15°C으로 조성하였다. 또한 식물체를 실제로 도복 시켜 발아가 적합한 저녁시간(21:00~09:00)의 관수 처리를 통하여 수발아 처리를 실시하였다(Figs. 1A, 1B) 수발아 처리가 완료된 두 품종에서 임의로 이삭 3수씩을 선정하여 수발아 피해 여부를 판단하기 위해 종자를 이삭에서 분리하여 scoring을 실시하였다. ‘금강’ 총 111개의 종자 중 75개가 발아하여 발아율 67.6%로 수발아에 감수성이 높은 것으로 확인되었다. 반면 ‘전주 377호’는 총 종자 수 81개 중 7개가 발아하여 발아율 8.6%로 확연한 발아 차이를 보여주었다(Figs. 1C, 1D). 수발아 저항성 적립계 품종인 ‘우리’의 경우 9.4%의 발아율(170종자 중 16개 발아)로 나타내며, ‘전주 377호’와 비슷한 발아율을 나타내었다.

Fig. 1. Pre-harvest sprouting induction experiment for RNA-seq analysis. (A) An artificial environment was provided to maintain optimal conditions for PHS induction. (B) The PHS induction experiment caused early germination of ‘Keumgang’. (C) Confirmation of germinated seed rate in ‘Keumgang’ Day after PHS+7. (D) Confirmation of germinated seed rate at Day after PHS+7 of ‘Jeonju 377ho’.

이는 ‘전주 377호’가 모계 품종인 ‘금강’에 비해 수발아 처리시 발아율이 현저히 낮아 지는 것을 확인할 수 있으며, 수발아 저항성 품종인 ‘우리’와 비교하였을 때도 비슷한 양상의 발아율을 보였기 때문에 ‘전주 377호’는 충분히 수발아 저항성 이 높은 특성을 확인할 수 있었다. 이번 연구에서 진행된 수발아 처리는 일반적으로 수발아 검정에 자주 활용되고 있는 모래묻이 검정법(Baier et al. 1987)과 인공 수조 수발아 실험 결과와 비교하였을 때 비슷한 양상으로 나타났다(Park 2020, Lee et al. 2021).

대량 전사체 분석과 clustering 분석

대량 전사체 분석을 위하여 ‘금강’과 ‘전주 377호’를 수발아 처리 후 1일이 경과된 시기를 K1, KM1으로 하였으며, 수발아 처리 후 7일이 경과된 ‘금강’과 ‘전주 377호’ 샘플을 각 K7, KM7로 정하였다.

두 품종 간 수발아 처리 후 서로 차등하게 표현 유전자(DEGs) 탐색 및 관련 유전자의 functional annotation은 4개의 샘플에서 RNA-seq을 통해 생성된 총 326,224,844개의 short reads를 이용하여 분석하였다. Quality 체크와 trimming 소프트웨어인 SolexaQA package의 dynamicTrim을 활용하여 read의 phred score ≥20를 기준으로 cut-off 하였으며 LengthSort data 분석을 통해 Trimming 후 short read length ≥25 bp에 해당되는 데이터를 제외하여 cleaned read 작업을 실시하였다. 각 샘플 별 total transcripts 서열에 cleaned read를 mapping하여 발현 값 계산 및 normalization을 수행하여 최종적으로 4가지 샘플은 평균 85.21% (234,131,980)의 mapping coverage로 해당 샘플에 대한 기능 분석과 DEGs 선발에 나타내는 것을 확인하였다(Table 1).

Table 1

Number of raw data reads and mapping reads

Sample ID	Total reads	Aligned 0 times		Aligned exactly 1 time		Aligned ≥1 times		Mapping rate
		Reads (ea)	Percent (%)	Reads (ea)	Percent (%)	Reads (ea)	Percent (%)	Reads (ea)	Percent (%)
K1	26,058,109	2,426,020	9.31%	1,444,724	5.54%	22,187,365	85.15%	23,632,089	90.69%
	26,058,109	2,523,271	9.68%	1,496,487	5.74%	22,038,351	84.57%	23,534,838	90.32%
KM1	33,131,971	4,040,234	12.19%	2,086,507	6.30%	27,005,230	81.51%	29,091,737	87.81%
	33,131,971	4,018,507	12.13%	2,093,121	6.32%	27,020,343	81.55%	29,113,464	87.87%
K7	41,079,207	4,473426	10.89%	2,453,590	5.97%	34,152,191	83.14%	36,605,781	89.11%
	41,079,207	4,572,891	11.13%	2,505,519	6.10%	34,000,797	82.77%	36,506,316	88.87%
KM7	37,108,988	9,237971	24.89%	1,896,424	5.11%	25,974,593	70.00%	27,871,017	75.11%
	37,108,988	9,332,250	25.15%	1,921,987	5.18%	25,854,751	69.67%	27,776,738	74.85%
4ea	274,756,550	40,624,570	14.79%	15,898,359	5.79%	218,233,621	79.43%	234,131,980	85.21%

Annotation은 이전 연구에서 de novo assemble 분석을 통해 형성되었던 금강 representative transcripts 104,542개를 대상으로 진행하였으며(Lee et al. 2021), 이 중 발현 값을 가지는 representative transcripts 개수는 91,097 개이고 functional description을 가진 유전자의 개수는 32,786개로 확인하였다. 이를 대상으로 functional annotation analysis를 진행하였으며, transcripts를 이용하여 KOG 39,072 개, KEGG 분석에 8,006 개, Gene ontology 분석에 44,631 개 사용되었다.

KOG 분석 결과 총 26 항목의 functional annotation이 확인되었으며, 특히 general function prediction only에서 가장 많은 6,969개의 transcripts의 분류되었으며, signal transduction mechanisms (4,420), Posttranslational modification, protein turnover, chaperones (3,868), Translation, ribosomal structure and biogenesis (3,487) 등에 annotation이 47.9%로 절반에 가까운 수치를 나타내었다(Fig. 2).

Fig. 2. Functional annotation classification was shown using KOG (eukaryotic orthologous groups) analysis based on representative transcripts. As a result, a total of 26 functions were identified.

KEGG 분석 결과 Metabolism (3,449), Genetic Information Processing (3,673), Cellular Processes (406), Environmental Information (281), Organismal systems (125), Human Diseases (45)의 6가지 항목으로 크게 나뉘었으며, Genetic Information Processing에서도 Ribosome (1,713), RNA transport (360) 등으로 이루어진 Translation와 Folding, sorting and degradation에 해당되는 Proteasome (215), Protein processing in endoplasmic reticulum (200)과 관련된 transcripts가 주로 확인되었다. Metabolism의 카테고리에서는, 탄수화물 대사(Carbohydrate metabolism), 아미노산 대사(Amino acid metabolism), 에너지 대사(Energy metabolism), 지질대사(Lipid metabolism)와 운송 및 이화작용(Transport and catabolism)의 순으로 많이 확인되었다(Fig. 3).

Fig. 3. Analysis metabolism and pathways were confirmed through KEGG (Kyoto encyclopedia of genes and genomes) analysis targeting representative transcripts.

KEGG 분석을 통해 RNA의 번역과 단백질 합성에 관련된 transcripts의 탐색과 carbohydrate metabolism에 해당되는 에너지 대사인 ATP 합성과 관련된 대표적인 대사 물질인 glycolysis, TCA cycle 관련 transcripts가 탐색되었으며, 특히 starch and sucrose metabolism과 관련된 transcripts가 발현되는 것을 확인하였다(Jiang et al. 2015). Glycolysis는 식물 세포에서 에너지 생산을 위한 주요 경로 중 하나이며, 이 경로에 관여하는 핵심 효소는 glyceraldehyde-3-phosphate dehydrogenase- (GAPDH)로 세포 ATP 수치와 탄수화물 대사(carbohydrate metabolism)의 균형 유지에 필수적인 역할을 한다. GAPDH는 glycolysis의 촉매 역할을 하며 벼 종자의 발아 후 GAPDH의 축적이 이루어짐을 확인한 바 있다(Kim et al. 2009, Zhang et al. 2017). 탄수화물, 지방, 단백질 대사를 통일하는 핵심 대사 경로인 TCA 회로 또한 멜라토닌(melatonin)에 의한 조절에 따른 에너지 이화작용을 통해 종자 발아에 영향을 미치는 것으로 알려져 있다(Zhang et al. 2017).

곡물에서의 전분(starch)은 식물의 1차 탄소 저장고 역할을 하며, 발아에 관련하여 중추적인 에너지원으로서 작용한다(Kossmann & Lloyd 2000). 전분 축적은 종자의 발달(seed development)과 동시에 진행되며, 출수 후(day after fertilization) 35일까지 최대치를 나타내고, 이후 중단되는데, 그 시기에 수발아 발생은 전분 분해에 의한 당 합성을 유도한다(Olsen 1999, Olsen 2001, Sabelli & Larkins 2009, Sreenivasulu et al. 2010, Thitisaksakul et al. 2012). 수발아를 포함한 지속적인 비생물적 스트레스는 전분 및 당 관련 대사에 영향을 주어 광합성을 통한 물질 생산에 손상을 주기도 한다. 또한 α-amylase 활성 과도 밀접한 연관이 있어 종자에서의 에너지 공급 및 발아 메커니즘에 중추적인 대사 작용을 한다(Andreoli et al. 2006, Detje 2008, Thitisaksakul et al. 2012, Kim et al. 2021). 주요하게 탐색된 pathway와 대사작용 및 이화작용 관련 transcripts를 참고하였을 때 수발아 실험에 따른 종자 발아 관련 메커니즘이 활발히 이루어졌음을 확인할 수 있다.

Clustering은 4가지 조건의 샘플을 대상으로 진행되어 두개의 그룹으로 나뉘었다. Cluster setⅠ은 샘플의 처리 시기별 차이를 기준으로 정하였고 유전자 발현 패턴은 K1 vs K7 와 KM1 vs KM7 비교 조합의 순으로 column을 정렬하여 DEGs를 선발하였다. 선발된 유전자 4,919개 유전자를 이용하여 clustering 분석을 수행하였으며 이 중 기능 주석이 달린 DEGs는 3,480개가 확인되었다. Cluster set Ⅱ는 샘플의 동일한 처리 조건에서의 품종 간 차이를 기준으로 분석하였다. Cluster setⅠ과 동일한 방법을 통해 비교 샘플 K1 vs KM1 와 K7 vs KM7 비교 조합 순으로 정렬된 column으로 분석하여 선발된 DEGs 총 2,133개 유전자를 이용하였으며 그 중 1,240개의 annotated DEGs가 선발되었다(Table 2).

Table 2

The number of DEGs detected through clustering set I, II analysis and the number of annotation DEGs

Group	Cluster I		Cluster II
	Num. of DEGs in cluster	Annotated DEGs	Num. of DEGs in cluster	Annotated DEGs
1	752	495	365	240
2	2,545	1,866	640	282
3	79	44	461	308
4	1,058	757	504	300
5	396	271	84	64
6	89	47	79	46
Total	4,919	3,480	2,133	1,240

Cluster set I, Ⅱ의 비교 조합 간 데이터 표현은 heatmap을 통해 상대적 발현 값으로 표현되었으며, heatmap에서 표현된 cluster를 line plot을 활용하여 분석한 결과 두 cluster에서 6개의 다른 패턴에 따라 그룹으로 분류하여 표현되었으며 그룹별로 해당 DEGs가 분포되어 있다(Fig. 4). 특히 Cluster set Ⅱ의 2번째 그룹은 cluster 내에서 가장 많은 분포를 보이는 DEGs (640개)로 구성되어 있으며, 수발아 처리 전 ‘금강’에 비해 ‘전주 377호’에서 유전자 발현이 높게 나타나고 발아 후 발현은 ‘금강’과 비슷한 수준이거나 감소하는 추세를 나타내는 그룹이다. 이는 종자 휴면과 가장 밀접한 관련이 있는 발현 패턴을 보여주고 있다(Fig. 4B).

Fig. 4. Heatmap and line plot analysis for cluster sets I and Ⅱ. (A) The expression patterns for PHS treatment of Keumgang (K_1 vs K_7), Joenju 377ho (KM_1 vs KM_7) were compared and analyzed, and each showed six expression patterns. (B) The differences in expression between cultivars caused by PHS treatment before and after (K_1 vs KM_1 and K_7 vs KM_7) were compared and analyzed, and each of the six expression patterns was shown.

Gene ontology (GO) 분석

Gene Ontology 분석은 3가지 functional category인 Biological Process (BP), Cellular Component (CC), Molecular Function (MF)의 항목으로 구분하여 분석하였다. K1 vs KM1 와 K7 vs KM7 비교 조합인 cluster Ⅱ에서 수발아 전에 높은 발현을 나타낸 2번 그룹을 대상으로 실시하였다. Cluster Ⅱ에서는 3가지 항목에서 각 78개, 18개, 22개의 GO term을 확인했으나 2번 그룹에 해당되지 않는 GO term은 분석에서 제외하여 BP에서 49개, CC, MF에서 각각 12개의 GO Term으로 표현된 관련 DEGs 총 1,036개를 분석하였다(Fig. 5). 그룹 2에서는 BP 카테고리에 해당되는 DEGs가 주로 탐색되었는데 특히 종자 발아와 관련된 carbohydrate metabolic process (GO:0005975), cellular carbohydrate metabolic process (GO:0044262)의 탄수화물 대사 유전자 56개와 response to hormone (GO:0009725) 33개로 비교적 많은 유전자를 나타냈는데 이는 KEGG에서 표현된 결과와 비슷한 양상을 보여주었다. 그 외에도 종실의 수분 흡수와 관련된 response to osmotic stress (GO:0006970), shoot system development (GO:0048367)에서 각각 19개, 12개의 DEGs가 확인되었다. 이 중 종자 휴면과 발아와 관련된 기능인 embryo development ending in seed dormancy (GO:0009793)에서 2개 의 DEGs를 탐색하였으며, seed germination (GO:0009845)에 해당되는 DEG로 기존 수발아 연구에서 휴면 기작을 통해 저항성 유전자로 인정 받는 MFT 유전자가 유일하게 그룹 2에서 확인되었다.

Fig. 5. Gene Ontology (GO) analysis is divided into three functional categories: Biological Process (BP), Cellular Component (CC), and Molecular Function (MF). In cluster II, a total of 1,036 related DEGs expressed as 49 GO Term in BP and 12 GO Term in CC and MF were analyzed.

CC 카테고리에서는 세포와 세포질과 관련된 intracellular membrane-bounded organelle (GO:0043231)와 cytoplasm (GO: 0005737)의 GO term이 두드러지게 분포하였는데, 알려진 바에 따르면 종자에서는 단백질 저장 배액(protein storage vacuoles; PSV)이라고 불리는 세포 소기관(membrane-bounded organelles)에 단백질을 저장하며 배아에서 형태학적, 기능적 재구성을 통해 씨앗 발달 중에 형성된다(Feeney et al. 2018, Kusumaatmaja et al. 2021). Molecular function (MF)에서는 다른 두 카테고리에 비하여 유전자의 분포가 적게 나타났다. 59개 유전자는 cation binding (GO:0043169)과 연관 있으며 핵산 결합과 관련된 nucleic acid binding (GO:0003676), nucleotide binding (GO:0000166)는 45개의 유전자가 분포하는 것을 확인하였다. 또한 MF 카테고리 중 hydrolase activity, acting on glycosyl bonds (GO:0016798), transferase activity, transferring glycosyl groups (GO:0016757)도 GO term에 탐색되었다. Glycosyl groups은 탄수화물 합성의 다양한 매개체로서의 역할을 하며 최근 대사작용 및 효소 결합, 약물 분자학에서도 관심을 갖고 있다(Roy et al. 2016). 식물과 종자 발달 과정에서의 glycosyl group 중 하나인 glycosyltransferase는 uridine diphosphate glucose (UDP-glucose)에 의해 촉매되는데 이는 식물 발아와 큰 관련이 있다(Mock & Strack, 1993). UDP-glucose:sinapate glucosyltransferase (SGT) 효소 활성에 의해 생성된 1-O-sinapoylglucose는 종자 발달 과정에서 sinapoylcholine으로 가수분해되어 세포 내 효소 반응을 일으킨다(Shirley & Chapple 2003). 이러한 과정은 궁극적으로 발아에 영향을 주는 것으로 알려져 있으며(Milkowski & Strack 2010, Roy et al. 2016) 이는 모델 식물인 애기장대(Arabidopsis)와 십자화과인 유채(Brassica napus)에서도 확인한 바 있다(Milkowski et al. 2000, Meißner et al. 2008). Cluster Ⅱ에서 총 17개의 UDP-glucose 관련 DEGs를 확인하였으며, glycosyl groups 과 UDP-glucose의 상관 관계를 고려하였을 때 종자 발달과 발아에 깊은 연관이 있다고 판단된다. 상기 언급된 유전자 기능 별 주요 카테고리와 해당 DEGs는 supplementary Table 1에 자세히 정리하였으며, 종자 휴면과 밀접한 연관이 예상되는 amylase와 UDP-glucose 관련 DEGs를 supplementary Table 2에 나열하였다(Supplementary Table 1, Supplementary Table 2).

DEGs 선발 및 유전자 발현 검정

다른 비생물적 스트레스에 비해 수발아 발생은 환경적 요인과 식물 생장 및 억제 호르몬 간의 상호 작용, 종피색을 포함한 밀의 형태적 차이 등 복합적인 요인에 의해 영향을 받는다고 알려져 있지만(King et al. 2000, Nakamura et al. 2011, Huh et al. 2016, Holloway et al. 2021), 주 요인으로는 종실의 발아억제물질의 소거로 인해 휴면이 타파되어 발생되는 현상으로 알려져 있다(Nielsen et al. 1984). 이러한 선행 연구 결과를 참조하여 DEGs 발현 검정에 활용될 샘플로는 휴면 관련 DEGs를 중점적으로 탐색하기 위해 RNA-seq 용 샘플의 수발아 처리 이삭을 활용하였다.

후보 DEGs 선발은 종자 휴면에 대해 중점적으로 탐구하기 위해 진행되었으며 품종 간 발현 비교를 보기 위해 분석된 Cluster Ⅱ를 대상으로 발아 전 상태인 K1 vs KM1에서 up regulation log2 Fold change 3 이상, down regulation log2 Fold change 1.5 이하로 cut-off를 정하였다. DEGs 선발은 총 2,135개의 assembled transcripts를 대상으로 up-regulation 290개, down-regulation 296개를 선별하였으며, 이 중 gene description, E-value 값이 없는 DEGs 및 중복 DEGs를 제외하였다. Blast 분석을 통해 partial이 아닌 95% 이상 일치하는 유전자 중 보통밀(Triticum aestivum)이나 애기장대(Arabidopsis thaliana)의 TAIR 및 근연종인 벼(Oryza sativa)의 gene description 정보를 기반으로 K1 vs KM1 up-regulation 35개, down-regulation 26개를 선발하여 RT-PCR을 수행하였다. 위와 같은 기준으로 압축된 후보 DEGs군은 크게 두가지 패턴의 유의미한 발현 양상을 보이는 것으로 확인되었으며 총 18개의 후보 DEGs를 최종 선정하였다(Table 3). DEGs 발현 검정에 사용된 프라이머는 RNA-seq의 assembled total transcript 파일을 기반으로 EnsemblPlants blast search program (https://plants.ensembl.org/index.html)에서 수집하여 디자인하였다.

Table 3

Descriptions of candidate DEGs and primer design for gene expression analysis

Genes (Assembled transcripts)	Gene description	Forward (5' → 3')	Reverse (5' → 3')	Tm
K1 vs KM1 Up-regulation DEGs
Keumgang1SL003317t014	PROLM24 - Prolamin precursor	AGCAGCAAGGGCAAAGTTTC	TGTCATGGGGATGCTGTAGATG	60℃
Keumgang1SL001254t019	beta-amylase	GCCCTTCTCGAATTTGTTGTCG	ACGCAAGCCAGTGTTTGTAG	60℃
Keumgang1SL095837t001	ANAC087, Arabidopsis NAC domain containing protein 87	TCGTCAACCCACATTGCATG	TATTGGAAGCGGCATTTGCG	60℃
Keumgang1SL002056t010	MFT, PEBP family protein	AACAGCACCAGCACGTAGC	TGGCTGGTGGTTAACATACCTG	60℃
Keumgang1SL008016t001	Serine protease inhibitor, potato inhibitor I-type family protein	AGCAAGCTTGATGTGCATGG	GGAAAGTACATCCATGGCATGC	60℃
Keumgang1SL001264t011	ATCRA1,CRA1,CRU1, RmlC-like cupins superfamily protein	AGTGGTGTAGATGCTTGTAGCC	TGGCACAACTGTTTGGCATG	60℃
Keumgang1SL034539t010	Tubby-like F-box protein	AAAGATGTCGTCGCTTGTGC	CGACACCGGAAAAGTGAGTTTC	60℃
Keumgang1SL011205t005	ATSUC3,ATSUT2,SUC3,SUT2, sucrose transporter 2	ACCTTGCACAGTGTTGTTGG	TCTGCATTGCTGTTGTGGTC	60℃
K1 vs K7 Up-regulation DEGs
Keumgang1SL003578t009	Starch branching enzyme 2.2	GCTCTGATGTTTGGTGGACATG	TTCCGCTTTTTCCTCAACGC	60℃
Keumgang1SL001971t019	BETA-VPE,BETAVPE, vacuolar-processing enzyme precursor	TCGCCAAAAACGACCTCAAC	AGTAACCGCGGTTTTGTTCC	60℃
Keumgang1SL003424t003	fatty acid hydroxylase	TTTCCAGTTCTCCCTCCATTCG	ATTTGACGGTTCCCACATGC	60℃
Keumgang1SL014234t009	4-alpha-glucanotransferase	TGCTGGTTTTTGGGCAGTTC	ACAGCTTTGAACAGCGCATC	60℃
Keumgang1SL007804t008	ACT domain containing protein	TATTGCGGAGATTGAGGGCTAC	TTGCATTTGAGGCTCGCAAG	60℃
Keumgang1SL000284t005	AtCOR47,COR47,RD17, cold-regulated 47	CAAACACAACACACGCCAAC	TGCCTGGTTACCACAAGACAG	60℃
Keumgang1SL004791t001	thiaminC, thiamine biosynthesis protein	TGTGTGGTCCCAAGTTTTGC	AGCTTCACCATGTTGTTCGC	60℃
Keumgang1SL038902t003	CAMK includes calcium/ calmodulin depedent protein kinases	AGAACCTCATGTCCATGTACCG	TTCCAAGAGAGGGTTCGGAATG	60℃
Keumgang1SL007502t020	zinc finger, C3HC4 type domain containing protein	ACAATGGGCATGCCAAATGG	ACATTGCTTTGACGCTGAGG	60℃
Keumgang1SL000284t005	aldehyde dehydrogenase	AACCACTTTGACTGCAGTGC	ATGCTTGCATTGTGCTGGAC	60℃
Housekeeping gene
TaActin	Housekeeping gene	GCTGTTCCAGCCATCTCATGT	CGATCAGCAATTCCAGGAAAC	60℃

qRT-PCR을 활용한 발굴 후보 유전자의 발현 프로파일링 검정 결과 종자의 초기 발달과 발아에 영향을 주는 탄수화물 및 단백질 가수분해 관련 유전자와 휴면에 관련된 유전자가 다수 확인되었다. 크게 KM1이 K1에 비해 높은 발현을 나타내는 DEGs군과 K7에서 발현 패턴이 높게 나타나는 DEGs군으로 나뉘었다.

KM1이 높게 발현된 DEGs군 중 β-amylase (Keumgang1SL0012 54t019)는 종자 발아에서 가장 빠르게 발현되는 유전자 중 하나로서 전분 가수분해에서 1차적으로 작용하여 α-amylase의 발현이 높아지는 것으로 보고되고 있다(Goswami et al. 1977). 이는 β-amylase는 효소 활성을 조절하는 이항력(diastic power) 조절자로서의 역할을 수행하기 때문이며 종자 발달 동안 축적되어 종자 활력에 영향을 미친다(Giese & Hejgaard 1984, Ziefler 1999) 쌀과 보리의 배유에서의 발현을 확인한 바 있으며(Bilderback 1974, Okamoto & Akazawa 1979, Sricastava & Kayastha 2014) 밀 배반(scutella)에서도 발아 초기 α-amylase보다 더 중요한 역할을 하는 것으로 나타났다(Nandi et al. 1995). 특히 쌀의 초기 발아 중에 급격히 증가하였다가 이후 감소되는 것으로 보고된 바 있는데(Okamoto & Akazawa 1980) 이는 β-amylase (Keumgang1SL0012 54t019)의 qRT-PCR 발현 양상과 일치하였으며 ‘전주377호’의 KM1 샘플에서 β-amylase의 발현이 K1 샘플보다 높은 이유는 ‘금강’의 초기 종자 발아 기작이 더 빠르기 때문에 K1 β-amylase 발현이 이미 감소한 것으로 추론된다. 또는 β-amylase는 8시간 이내에서 합성되어 역할을 수행하는 것으로 수발아 감수성인 ‘금강’의 경우 활발한 발현 후 감소하였으나, ‘전주377호’에서는 β-amylase 활성이 24시간이 지난 시점에서도 여전히 발현이 높은 것으로 보아 β-amylase 활성 시기가 종자 발아 정도에 따라 다를 것으로 예상된다.

MOTHER OF FT AND TFL1 (Keumgang1SL002056t010; MFT)는 식물에서 phosphatidylethanoamine-binding protein (PEBP) 유전자군 속하며 주로 종자 발달 및 발아에서 기능하는 것으로 알려져 있다(Chardon & Damerval 2005, Tao et al. 2014). 수발아와 종자 휴면 관련 양적형질유전자좌(QTL; Quantitative Trait Loci)는 모든 염색체의 다수의 위치에서 보고되고 있으며 특히 염색체 3AS와 4AL에서 주동 QTL이 확인되었고(Kato et al. 2001, Mori et al. 2005) 그 중 QPhs.ocs-3A의 후보 유전자가 MFT로 밝혀지면서 품종 간 프로모터 영역 내의 SNP 차이로 인한 발아 지연을 확인한 바 있다(Nakamura et al. 2011). 또한 수발아 처리에 따른 식물생장호르몬 생합성(synthesis) 및 신호전달(pathway) 분석은 수발아에 대한 주요 요인이 ABA 및 ABI에 의한 조절을 제시한 바 있다(Lee et al. 2021). MFT의 qRT- PCR 발현이 K1에서 KM1보다 높게 나타나는데(Fig. 6A), 이 결과는 MFT는 종자 발아보다 휴면에 직접적인 영향을 나타내는 유전자이기 때문에 ‘전주 377호’에서는 또 다른 메커니즘이 관여했을 것으로 사료된다.

Fig. 6. Validation of PHS and seed dormancy related selected DEGs using quantitative real-time PCR analysis. (A) A group of candidate DEGs showed high expression on day 1 and low expression on day 7 of PHS treatment in ‘Jeonju 377ho’. (B) A group of candidate DEGs showed low expression on day 1 and high expression on day 7 of PHS treatment in ‘Keumgang’.

또다른 ABA 관련 유전자인 Tubby F-box like protein (Keumgang 1SL034539t010; TLP)은 다양한 진핵생물에서 자주 발견되며 C-말단 영역에 있는 Tubby 도메인을 가지고 있는 것이 특징이다(Kleyn et al. 1996, Liu 2008, Zhang et al. 2020). 모델 식물체인 애기장대와 벼에서 각각 11개, 14개의 TLPs family를 확인하였으며 그 중 애기 장대의 AtTLP9는 염과 한발 스트레스 반응에 관여함을 확인하였을 뿐만 아니라 AtTLP3와 AtTLP9 사이에서 중복적으로 ABA와 삼투압 반응에 관여하는 것 밝혀진 바 있다(Lai et al. 2004, Liu 2007). 최근 오이 CsTLP8를 애기장대에 도입된 과발현(overexpression) 검정을 통해 ABA 정도에 따른 기공 폐쇄가 민감하게 반응하는 것을 확인하였는데 이는 삼투압 스트레스에 의한 종자 발아를 조절을 확인할 수 있었으며 추가로 진행된 yeast two hybrid system을 통해 CsTLP8이 SCF complex (SKP-Cullin- F-box containing complex) gene인 SKP1과의 상호작용을 확인하였다(Li et al. 2021). CULL1과 F-box 단백질을 연결하는 핵심 성분인 TaSKP1은 SCF 복합체(SCF complex)를 구성하며 세포 주기 조절, 스트레스 반응, 신호 전달, 대사 조절, 세포 분화 등 다양한 세포 과정에 중요한 역할을 수행한다(Moon et al. 2004). 최근 F-box 단백질이 자스몬산 신호(jasmonate signaling), 꽃 형태, 잎 노화 등 식물의 전반에 걸친 SCF 대사경로 연구가 활발하다(Hershko & Ciechanover 1998, Moon et al. 2004, Hong et al. 2013). 밀의 경우 F-box 유전자에 대한 게놈 기반 연구는 복잡한 유전자 구조와 17 Gb의 밀 게놈 사이즈로 어려웠으나 최근 Tubby F-box like protein이 포함된 F-box family에 대한 연구가 이루어져 있는데 애기장대, 벼 등과 마찬가지고 가뭄 내성, 개화 등에 관여하며 특히 DEGs Clustering 패턴에 따라 구성된 결과 종자 발달 단계에서 높은 발현을 나타내었다(Hong et al. 2020).

Sucrose는 수용성 이당류로 탄수화물 분배를 위한 에너지원 역할과 광합성에 의한 주요 생산물 중 하나이며 식물의 초기 발달과 비생물적 스트레스 반응에도 관여한다(Lastdrager et al. 2014, Navabpour & Jalali 2019). 특히 체관부에서의 당 수송의 필수적인 역할을 하는 동시에 식물의 성장 및 저장 단백질 합성에 관여하는 Sucrose Transporter (Keumgang1SL011205t005; SUT) 역할은 매우 중요하다(Aoki et al. 2006). SUT에 의해 촉매되는 운반체 SUC는 turgor-gated mechanism에 영향을 주며 압력에 따라 종자 내 양분 흡수 제어 역할을 수행한다. SUC의 발현으로 세포내 압력이 높아져 아포플라즘의 용질 농도가 높아지면 양분 방출에 관여하는 운반체를 활성화시켜 발아가 지연되는 것으로 알려져 있다(Ayre 2011). 마찬가지로 SUC를 조절하는 SUT 또한 sink to source transition의 대사작용과 함께 조절되며 이로 인해 단백질 축적 조절과 조숙한 발아를 방지하는데 관여한다고 알려져 있으며(Wolswinkel 1992) 유채에서 종자 발달 초기에 전분과 지방의 축적에도 영향을 주는 것으로 확인하였다(Li et al. 2011). SUT/SUC 에 의한 Membrane-Transport Systems은 종자 발아에서 주요한 영향을 줄 것으로 예측되며 추후 초기 생육과 발아에 필수적 역할을 하는 SUT 관련 유전자를 대상으로 수발아 및 휴면에서의 연관성 검증이 필요할 것으로 판단된다. 선행연구에서 MFT와 SUT은 서로 상호 결합하는 단백질이므로 MFT와 SUT가 같은 발현 양상을 보이는 DEGs로 선발된 것으로 보아 DEGs 선발과 분석이 제대로 수행되었다고 판단한다(Park 2020).

같은 양상의 발현 패턴을 보여주는 DEGs인 Keumgang1SL 003317t014, Keumgang1SL095837t001, Keumgang1SL008016t001, Keumgang1SL001264t011에 대한 연구도 추가적으로 이루어져야 할 것으로 보여진다.

반면 K7에서 발현이 증가된 DEGs군의 Starch branching enzyme (Keumgang1SL003578t009; SBE)은 ‘금강’의 K7에서 발현이 높은데 이는 수발아가 빠르게 진행되고 있음을 보여준다(Fig. 6B). SBE는 구조적으로 α-amylase superfamily에 속하며, 전분의 대사 조절에 관련된 고유 영역의 아미노산 서열을 포함하고 있다(Tetlow & Emes 2014). 그렇기 때문에 SBE의 발현은 α-amylase가 종실에 저장된 전분을 충분히 작용할 수 있게 만든 것이며(Xia et al. 2011) ‘금강’의 K7에서 SBE의 발현이 높은 이유는 수발아에 대하여 감수성 임을 간접적으로 나타낸다.

Zinc finger protein (Keumgang1SL007502t020) 도메인은 비교적 작은 단백질 모티브로, 표적 분자와 탠덤(tandem) 접점을 만드는 여러 개의 손가락 모양의 돌기가 있으며 DNA 인식, RNA packing, 전사 활성화, 단백질 접힘 및 조립, 지질 결합 등의 역할을 수행한다고 알려져 있다(Lin et al. 2011). Zinc finger protein 모티브의 많은 슈퍼 패밀리가 있으며, 순서와 구조 모두 다양하다(Klug 1999, Matthews & Sunde 2002, Brown 2005, Hall 2005, Gamsiaeger et al. 2007). Zinc finger protein은 애기장대에서 CCH zinc finger protein 관련 유전자 AtTZF1, AtTZF4, AtTZF5, AtTZF6가 부분적으로 ABA와 GA 대사를 조절하여 종자 발아에 관여하는 것을 확인하였으며(Lin et al. 2011, Bogamuwa & Jang 2013). 벼에서도 OsTZF1가 종자 발아를 억제하는 것으로 나타났다(Jan et al. 2013). 밀의 경우 본 연구에서 확인된 C3HC4 유형은 아니지만 최근 같은 family인 CH3 유형에서 연구가 진행된 바 있으며 zinc finger protein CH3 식물의 성장과 발육을 조절하는 데 중요한 역할을 하며 애기장대나 벼에서 생물학적/비생물학적 스트레스 반응을 조절하는 것으로 알려져 있다(Li & Thomas 1998). 애기장대 zinc finger protein CH3 의 88개 유전자를 도입한 밀을 대상으로 수발아 처리하여 발현 검정을 실시하였으며 전사체 데이터, qRT-PCR, ABA 반응 분석을 통해 발아 억제 유전자 3개(TaC3H4/-18/-72)와 발아를 상향 조절하는 2개 유전자(TaC3H37/-51)가 종자 휴면과 발아에 밀접한 관련이 있는 후보 유전자로서의 가능성을 시사하였다(Cheng et al. 2020). 그 외 동일한 발현 양상을 나타내는 Keumgang 1SL001971t019, Keumgang1SL003424t003, Keumgang1SL014234t009, Keumgang1SL007804t008, Keumgang1SL000284t005, Keumgang 1SL004791t001, Keumgang1SL038902t003, Keumgang1SL000284t005의 수발아 및 휴면 연관성에 대한 유전자 연구가 지속적으로 이루어져야 할 것으로 보인다.

밀은 지역에 국한되지 않고 세계적으로 널리 소비되고 있는 식량 작물로 다양한 식료품의 원재료로 사용되고 있다. 특히 국내에서의 밀은 1990년 대부터 현대에 이르기 까지 1인당 소비량이 감소세 없이 소비되는 작물로 국내에서도 식량작물로서의 매우 중요한 위치에 자리잡고 있다. 하지만 밀 경작지의 급격한 감소와 자급률의 한계로 인한 불안은 여전히 지속되고 있고 있으며, 최신 지구온난화로 인한 이상 기후와 불안정한 국제 경제 시황은 식량안보에 큰 위협으로 다가오고 있다. 특히 전 세계적으로 이상기후로 대표적 수분 장해인 수발아의 발생이 빈번해지며 식량으로서 밀의 가치를 손상시키는 요인 중 하나이다.

본 연구에서는 국내 장려 품종 ‘금강’과 수발아에 강한 특징을 보이는 돌연변이 계통 ‘전주 377호’를 수발아 처리한 샘플을 대상으로 RNA-seq을 통해 mapping된 234,131,980 개의 reads를 활용하여 유전자 기능 주석과 DEGs 분석을 실시하였다.

De novo assembly를 통하여 조립된 transcript를 대상으로 KOG 및 KEGG 분석을 진행하였으며 탄수화물 대사, 종자 내 에너지 공급 및 발아 메커니즘에 중추적인 전분 및 당 분해 대사 등이 주로 나타나는 것을 확인하였다. 또한 GO 분석을 통해 탄수화물 대사와 식물에서 종자 발달 과정에 관여되는 glycosyl group의 uridine diphosphate glucose와 관련된 DEGs 총 21개를 Cluster Ⅱ에서 확인하였다.

탐색된 DEGs를 대상으로 형성된 cluster 별 분석으로 종자 휴면과 밀접할 것으로 예측되는 후보 DEGs를 선발하였으며, 발현 검정을 통하여 수발아에 의한 종자 내 생물학적⋅화학적 기작을 분석 하였다. 특히 종자 휴면에 관련된 MFT의 경우 지속적인 연구를 통해 종자 발아와 휴면 관련 메커니즘 구명이 필요하며 기존에는 잘 알려지지 않거나 다른 기능으로 알려진 유전자 Tubby F-box like protein, Zinc finger protein family, SUT 등은 종자 발아(seed germination)와 종자 휴면성(seed dormancy)의 관점으로 활발한 연구와 증명이 이루어져야 할 것으로 사료된다.

본 연구를 기초 자료로 활용하여 종자 발아와 휴면 간 메커니즘 이해와 유전적 요인을 파악한다면 미래 농업에서의 밀의 가치는 더욱 높아질 것으로 판단된다.

본 논문은 농촌진흥청 연구사업(과제명: 수발아 저항성 밀 육종소재 개발, 과제번호: PJ016528012022)에 의해 이루어진 것임