표현체 기술의 급격한 발전은 연구자들에게 방대한 이미지와 계측 데이터를 제공하여 작물 생육 과정에 대한 심층적인 이해를 가능하게 하였다. 이러한 기술적 진보는 단순한 데이터 양의 증가를 넘어, 농업 연구의 질적 향상에 기여하고 있다. 더 나아가, 공장화된 스마트온실인 LemnaTec(Aachen, Germany), Photon Systems Instruments(Drasov, Czech Republic), Phenospex(Maastricht, Netherlands) 등의 국제적 선도 기관과 함께, 국립농업과학원 표현형 연구 온실이 표현체 연구의 새로운 트렌드를 주도하고 있다(Demidchik et al. 2020; Lee et al. 2021). 스마트온실의 도입은 데이터 수집의 자동화와 정밀화를 가능하게 했지만, 동시에 대량의 데이터를 효과적으로 분석할 수 있는 최적화된 분석 기법의 필요성을 더욱 부각시키고 있다(Onnela 2021; Song et al. 2021). 그럼에도 불구하고 기계적 기술의 진보에 비해 최적화된 표현체 분석 기법은 명확하게 확립되어 있지 않다.

현재 농업 분야에서 사용되고 있는 표현형 측정 방법은 크게 세 가지로 나눌 수 있다. 첫 번째는 그래픽 사용자 인터페이스(GUI) 기반의 이미지 분석 오픈소스 툴을 활용한 방법으로, 어느 정도 사용자 개입이 필요하지만 자동화된 분석이 가능하여 일정 수준의 효율성을 제공한다(Baek et al. 2020). 대표적인 예로 ImageJ가 있으며, 수동측정에 비해 이미지 데이터를 쉽게 분석할 수 있지만, 복잡한 데이터 처리나 높은 정확도를 요구할 때는 한계가 있을 수 있다. 다음은 OpenCV(Open Computer Vision Library)나 PlantCV와 같은 컴퓨터 비전 분석 라이브러리를 활용하는 방식이다. 범용적인 OpenCV는 이미지 분석에 널리 사용되어 왔으며(Choi et al. 2020), 최근에는 식물 표현체 분석을 위해 특화된 PlantCV가 주목받고 있다(Gehan et al. 2017). 그러나 이러한 분석 방법에도 불구하고, 표현형 분석의 어려움은 방대한 시계열 데이터의 처리에 있다. 식물의 성장 과정에서 생성되는 표현형 시계열 데이터는 계절성, 변동성, 주기성과 같은 특성을 반영하며, 이를 분석하기 위해서는 전문적인 기술과 농업 도메인 지식이 필수적이다. 또한 표현형 데이터의 자기상관성으로 인해 시각적, 시간적 요소를 포함한 복잡한 해석이 요구되며, 일반적인 통계적 접근법으로는 한계가 존재한다(Coppens et al. 2017).

본 연구에서는 표현형 분석의 최적화를 위해 측정 방법 간의 차이와 데이터 가공에 따른 변화를 비교 분석하였다. 벼의 전 생육주기 동안 촬영된 이미지 데이터를 기반으로 높이 형질을 세 가지 도구를 통해 추출하고, 최적 방법으로 추출된 데이터에서 반복 샘플, 다각도 촬영, 이상치 제거의 효과를 분석하였다. EDA(Exploratory Data Analysis)를 거친 표현형 데이터는 유전체 데이터와 연관 분석을 수행하여 벼 높이 관련 후보 유전자를 선발하였고, 선발된 유전자들은 기존 연구와 비교하여 확인하였다. 이를 통해 표현형 시계열 데이터의 효율적인 처리 및 분석 방향을 제시하고, 향후 유전체 연구와의 통합적 접근에 기여할 수 있을 것으로 기대된다.

표현체 시계열 데이터 수집

분석에 사용된 재료는 한국산 벼 276품종 중, 국립식량과학원의 수량구성요소 안정성 분석을 통해 선발된 100품종이다(Lee et al. 2023). 이들 100품종은 모두 자포니카 생태형(Japonica ecotype)에 속하며, 조만성을 기준으로 조생종 30품종, 중만생 40품종, 중생 30품종으로 구성되었다. 선발된 품종들은 국립농업과학원의 스마트 온실 영상 촬영실에서 3차원 영상분석장비(3D Scannalyzer imaging system, LemnaTec, Germnany)를 통해 전 생육주기인 122일(2023년 06월 14일부터 2023년 10월 18일) 동안 이틀 간격으로 촬영되었다. 작물 생육 상태 이미지를 확보하기 위해 표현체 분야 국가참조표준데이터센터의 설정된 조건에 따라 해상도 6576 × 4384 픽셀, 400~700 nm 파장의 가시광선 영역의 CCD 카메라가 사용되었다(Lee et al. 2021). 각 품종 당 4개의 반복샘플과 0°, 120°, 240°의 세 각도에서 촬영된 이미지 데이터를 획득하였다.

PlantCV를 통한 벼 높이 측정

수집된 원시데이터는 파이썬(Python)기반 PlantCV 라이브러리를 활용하여 벼의 이미지가 분석되었으며(Hodge et al. 2021), 여러 필터링 및 보정 기법을 적용하여 노이즈를 제거하고 식물체의 경계선을 검출하였다. 먼저, 입력 이미지에서 배경 제거와 관심 영역(ROI) 설정을 통해 분석 대상을 추출하였다. 이후, 이미지의 색상 범위를 설정하기 위해 HSV 색공간(Hue: 색상, Saturation: 채도, Value: 명도)으로 변환 후, 지정된 색상 범위 내에서 대상 영역을 필터링하고 이진화 처리하였다. 이진화된 이미지를 바탕으로 관심 영역을 마스킹 기법을 통해 추출하였다. 다음으로, 양방향 필터(Bilateral Filter)를 사용하여 노이즈를 제거하면서 경계선을 보존하였다. 이어서 팽창(Dilation)과 침식(Erosion) 과정을 적용하여 경계선을 뚜렷하게 복원하고 잔여 노이즈를 제거하였다. 마지막으로, 벼의 윤곽선을 추출하였으며, 가장 큰 윤곽선을 기준으로 지상부 벼 높이(projected plant height)를 측정하였다. 측정된 값은 화분의 실제 길이를 기준으로 픽셀 단위에서 실측 단위(cm)로 변환하였다. PlantCV 측정 데이터는 국립농업과학원의 ImageJ(ver. 1.53j) 영상분석 데이터와 경북대학교의 OpenCV기반 분석(ver. 4.8.1)데이터와 결합하여 전체 표현형 데이터세트로 통합하였다.

표현형 데이터의 이상치 제거

벼의 높이에 대한 최적화된 시계열 분석을 위해, 데이터를 다섯 가지 케이스로 세분화하였다. 첫 번째는 각 품종별로 정면에서 촬영된 1개 샘플을 랜덤하게 선발하여 사용하였고, 두 번째는 각 품종별로 4반복 샘플 데이터를 사용하였다. 세 번째는 4반복 데이터에서 사분위 범위(IQR, Interquartile Range)를 이용하여, Q1 - 1.5IQR보다 작거나 Q3 + 1.5IQR보다 큰 값을 이상치로 제거한 후 데이터를 사용하였다. 분위수 계산은 Python의 pandas 라이브러리의 quantile 함수를 사용하여 기본 선형 보간법을 적용하였다. 네 번째 케이스는 세 각도 중 최대값을 추출한 4반복 샘플 데이터를 사용하였으며, 마지막 케이스는 각도별 최대값을 추출한 후 IQR을 이용해 이상치를 제거한 데이터를 사용하였다. 이러한 케이스별 데이터 처리 방식이 벼 높이 시계열 분석에 미치는 영향을 평가하고자, 평균 절대 편차(MAD), 표준 편차(STD), 평균 절대 백분율 오차(MAPE), 평균 제곱근 오차(RMSE), 변동 계수(CV) 등의 지표를 사용하여 데이터를 평가하였다.

GBS(Genotyping-by-sequencing)

표현체 연구에 사용된 100품종 중 96품종을 GBS 분석에 사용하였다. 사용된 96품종은 조생종 27품종, 중만생 40품종, 중생 29품종으로 구성되었다. 각 품종의 genomic DNA는 Thermo Scientific Nanodrop 8000 분광광도계(Fisher Scientific; Hampton, NH, USA)를 사용해 농도 100 ng/μL, 순도 260/280 = 1.8~2.0, 260/230 = 1.6 이상의 genomic DNA로 정량화하였다. GBS라이브러리를 준비하기 위해 DNA를 제한효소 PstI (CTGCAG)와 MspI (CCGG)로 절단한 후, 각 품종의 DNA에 96개의 바코드 어댑터 세트에 연결하였다(Elshire et al. 2011). 이후 96품종의 라이브러리를 풀링하고, Illumina 플랫폼 Hiseq 및 NovaSeq pair-end sequencing을 사용하여 151 bp 읽기로 시퀀싱을 진행하였다. 바코드 시퀀스를 이용하여 demultiplexing을 수행하고, Adapter trimming과 sequence quality trimming를 진행하였다. Adapter trimming은 cutadapt(ver 1.8.3; Martin 2011)을, sequence quality trimming은 Trimmomatic(ver. 0.39; Bolger et al. 2014)을 사용하였다. 전처리된 trimmed reads는 BWA-MEM(Vasimuddin 2019)을 사용해 벼 표준유전체(Os-Nipponbare-Reference-IRGSP-1.0; Kawahara et al. 2013)에 맵핑하였다. Clean reads는 씨더스 자체 데이터로 finetuning한 DeepVariant (ver 1.5; Poplin et al. 2018)와 GLnexus(ver 1.4.1; Yun et al. 2020)를 사용해 통합 SNP matrix를 생성하였다.

유전체-표현체 연관 분석

최종 필터링 된 총 12,127개의 SNPs와 표현형 데이터 중 최대 높이 형질을 이용하여 GWAS 분석에 활용하였다. GWAS 분석은 GLM, MLM, MLMM, BLINK, Farm CPU, SUPER를 사용하여 수행되었으며, 모든 매개변수는 GAPIT R 패키지(ver 3; Wang and Zhang 2021)의 기본값으로 설정되었다. 본 연구에서 제1종 오류를 방지하기 위한 유의성 임계값으로 Bonferroni 교정법을 적용하여, p-value를 0.01로 설정하였다. 선발한 SNP의 주변 영역 200 kb 내에서 유전자 사이의 영역을 제외하고 유전자 영역에 위치한 좌에 대해 RAP-DB(The Rice Annotation Project Database; Sakai et al. 2013)와 Trait Ontology (Arnaud et al. 2012)를 포함한 기존 문헌을 통해 형질 마커의 가능성을 확인하였다.

표현형 분석 도구 별 벼 높이 성장곡선 비교

본 연구에서는 벼 96품종의 전 생육주기 동안 촬영된 시계열 데이터를 바탕으로, ImageJ, OpenCV, PlantCV 세 가지 도구를 사용하여 벼 높이 성장곡선을 비교하였다. 표현형 분석 도구별 벼 높이 성장곡선의 결과는 Fig. 1과 같다.

Figure 1. Phenotypic time-series trend comparison of rice plant height across three different phenotyping tools: (A) ImageJ, (B) OpenCV, and (C) PlantCV. Individual growth curves of all rice plants in the dataset (N = 384, representing 96 cultivars with four replicates each) are shown in the top plots, with maximum, minimum, and mean values obtained by the phenotyping tools displayed below. The red arrow highlights early-stage irregularities in ImageJ, while the blue arrow indicates excessive uniformity in values obtained by OpenCV at later stages of plant growth

ImageJ는 20일 이전에 다른 두 방법에 비해 높은 불규칙성으로 상대적으로 넓은 분산을 나타냈다(Fig. 1A). 이는 ImageJ의 비교적 단순한 Global thresholding과 Noise handling 기능으로 인한 오차로 판단된다(Suarez-Arnedo et al. 2020; Woolf et al. 2021). 초반 성장곡선에서 OpenCV와 PlantCV는 유사한 패턴을 보였으나, 45일 이후 OpenCV는 최대 높이에 도달한 후 일정하게 유지되는 경향을 보였다(Fig. 1B). OpenCV의 기본 설정값은 edge detection 알고리즘에 의존하므로 성장 초반의 녹색이 진한 시기에는 높은 정확도를 보였으나, 지지력이 약한 화본식물의 특성과 황화로 인한 색 변화의 효과가 커지는 성장 후반에는 오차가 발생한 것으로 판단된다(Jing et al. 2022). OpenCV의 다양한 변수 설정 기능과 ImageJ의 픽셀 카운팅 방식의 장점이 성장 시점별로 다르게 적용되어 표현형 분석 도구 별 차이가 나타났다. PlantCV는 모든 단계에서 균일하고 안정적인 결과를 보여주었으며, ImageJ와 OpenCV의 장점을 통합한 분석 도구로서의 유용성을 입증하였다. 이러한 결과는 표현형 분석에서 충분한 도메인 지식과 local parameterization의 중요성을 강조한다(Gehan et al. 2017).

ImageJ와 OpenCV의 성장곡선 상의 오차는 전체 샘플의 최대값, 최소값, 평균값을 기준으로 계산한 통계 지표를 통해 확인하였다(Table 1). ImageJ의 최대-최소값 사이 면적(AUC)과 전체 평균-최대값의 거리(MaxDist)는 각각 119,160과 859.11로 가장 높았다. OpenCV의 전체 평균-최대값의 상관관계(MaxCorr)와 전체 평균-최대값의 거리(MaxDist)는 각각 97.62%, 577.15로 가장 낮았다. 또한, 세 분석 방법 모두에서 전체 평균-최소값의 상관관계(MinCorr)와 전체 평균-최소값의 거리(MinDist)가 1% 내외의 오차범위에 있음을 확인함으로써, 최소값 측정의 유사성을 입증하였다. 성장곡선 추세에 대한 통계에서, ImageJ는 전체 평균값과 최대값이 멀고, OpenCV는 상대적으로 가까우며, PlantCV는 결과가 극단으로 치우치지 않고 안정적인 값을 보였다. 이는 PlantCV가 과도한 변동을 피하고 균형 잡힌 경과를 제공할 가능성을 시사한다.

Table 1 . Statistical comparison of rice plant height obtained from different phenotyping tools.

Tool	ABCa	MinCorrb (%)	MaxCorrc (%)	MinDistd	MaxDiste
ImageJ	119,160	99.36	98.01	673.06	859.11
OpenCV	97,173	99.31	97.62	662.80	577.15
PlantCV	116,126	99.31	98.72	663.08	821.87

^aABC: Area between curves; ^bMinCorr: Correlation of the min. values curve with the sample median; ^cMaxCorr: Correlation of the max. values curve with the sample median; ^dMinDist: Mean distance of the min. values curve from the sample median; ^eMaxDist: Mean distance of the max. values curve from the sample median..

반복 샘플, 촬영각도, 이상치 제거에 따른 데이터 가공 효과

벼 96품종의 전 생육주기 시계열 데이터 총 70,272장을 이용하여 반복샘플, 촬영각도, 이상치 제거가 데이터 가공에 미치는 효과를 분석하였다. 시간에 따른 벼 높이를 성장곡선으로 시각화 하여 전체적인 경향 변화를 확인하였다(Fig. 2). 먼저 4개의 반복샘플과 세 가지 촬영 각도를 기준으로 데이터의 변화량을 확인하고 산술평균을 이용하여 분산을 감소시키는 방향으로 분석하였다(Fig. 2A). 사분위수를 통한 이상치 제거에서 시계열 데이터의 연속성이 끊어지는 현상을 확인하였으며(Fig. 2B, C), 이를 보완하기 위해 각도별 최대 높이의 데이터를 선발하여 분석에 활용하였다. 이를 통해 20일차 데이터에서 발생한 성장곡선의 불규칙성이 보정되었고, 성장곡선의 평활화(Smoothing)효과를 확인할 수 있었다(Fig. 2D, E).

Figure 2. Phenotypic time-series trend comparison of rice plant height across three different data processing methods: (A) 96 cultivars × 4 replicates × 3 angles, (B) 96 cultivars × 4 replicates, (C) 96 cultivars x 4 replicates in the max angle, (D) 96 cultivars × 4 replicates with IQR filter, and (E) 96 cultivars × 4 replicates in the max angle with IQR filter. The gray line is shown on the left plot in (A) and superimposed on subsequent plots, representing the growth curve from all data points (N = 1152) across all data processing methods, and is used as a reference point to assess data processing accuracy. Growth curves for all individual samples in the dataset (N = 384, representing 96 cultivars with four replicates each) are shown on the left side of plots B to E, while mean growth curves for each cultivar (N = 96) are shown on the right. The black line on each plot depicts the average growth curve across all plants for each data processing method. Arrows indicate cases in which data trends changed as a result of data processing

시계열 데이터를 성장곡선으로 표현하는 것은 추세를 파악하는 데 유용하지만, 많은 샘플로 인해 통계적 비교는 별도의 접근이 필요하다(Sausen et al. 2021). 따라서 가공 단계마다 각 샘플별 평균에서 반복샘플과 촬영 각도로 인한 분산 및 오차를 분리하여 분석을 진행하였다(Table 2).

Table 2 . Statistical comparison of rice plant height across from different data processing.

Data Processing	MADa	STDb	MAPE (%)c	RMSEd	CV (%)e
A. 96 Cultivar × 4 Replicates × 3 Angles	94.95	117.1	4.83	112.12	5.97
B. 96 Cultivar × 4 Replicates	94.37	125.31	4.77	108.52	6.34
C. 96 Cultivar × 4 Replicates with IQR filter	81.61	107.25	4.11	92.48	5.41
D. 96 Cultivar × 4 Replicates in Max Angle	90.44	119.82	4.48	103.76	5.94
E. 96 Cultivar × 4 Replicates in Max Angle with IQR filter	78.91	103.47	3.88	89.24	5.1

^aMAD: Mean Absolute Deviation; ^bSTD: Standard Deviation; ^cMAPE: Mean Absolute Percentage Error; ^dRMSE: Root Mean Square Error; ^eCV: Coefficient of Variation..

먼저 반복샘플의 효과를 분석하기 위해, 전체 1152개 데이터와 촬영 각도를 고려하지 않은 384개의 데이터를 비교하였다. MAPE는 Table 2A에서 4.83%, Table 2B에서는 4.77%로 소폭 감소하였다. RMSE는 Table 2A에서 112.12, Table 2B에서 108.52로 감소한 반면, CV는 Table 2A에서 5.97%, Table 2B에서는 6.34%로 증가하였다. 이는 촬영 각도를 제외했을 때 평균 오차는 줄었으나 데이터의 변동성은 증가했음을 나타낸다. 즉, 일부 값은 평균에 가깝지만 다른 값은 크게 벗어나, 여러 촬영 각도의 데이터가 변동성 감소에 기여함을 의미한다. 이를 통해, 반복샘플만으로는 전체 데이터 세트의 변동성을 충분히 제어하는 데 한계가 있음을 알 수 있었다.

다음으로, 반복 개체에서 최대 높이를 나타낸 촬영 각도의 데이터를 선발한 Table 2D와 촬영 각도 데이터를 고려하지 않은 Table 2B를 비교하였다. CV는 Table 2B에서 6.34%, Table 2D에서 5.94%로 약간 감소하였다. MAPE는 Table 2B에서 4.77%, Table 2D에서 4.48%로 감소하여 오차가 줄어든 것을 확인하였다. RMSE와 STD 또한 감소함으로써, 최대 높이를 가지는 촬영 각도의 데이터를 선발하는 것이 데이터의 변동성과 오차를 줄이는데 효과적임을 알 수 있었다. 더 나아가, 모든 각도를 사용한 Table 2A와 최대 높이를 가지는 각도의 데이터를 선발한 Table 2D를 비교한 결과, CV는 5.97과 5.94로 거의 차이가 없었지만 MAPE가 4.83에서 4.48로 감소하면서 유의미한 데이터 선발이 이루어졌음을 확인하였다.

사분위수로 이상치를 제거한 Table 2C는 CV에서 Table 2B와 두드러진 차이를 보였다. Table 2B에서 CV는 6.34%였으나, IQR 필터를 적용한 Table 2C에서는 5.41%로 감소하였다. 이는 필터링을 통해 데이터의 변동성이 크게 줄어들었음을 의미한다. MAPE 역시 Table 2B에서 4.77%였던 것이 Table 2C에서 4.11%로 감소하여 예측 오차가 개선되었음을 확인하였다. RMSE는 Table 2B에서 108.52, Table 2C에서 92.48로, STD는 125.31에서 107.25로 크게 감소하였다. 이는 이상치 제거를 통해 극단값을 제거함으로써 변동성과 오차가 동시에 감소했음을 시사한다. 따라서 반복 샘플 분석에서 사분위수를 통한 이상치 제거가 데이터의 일관성과 정확성을 높이는데 유리함을 확인하였다.

마지막으로, 각 반복 개체에서 최대 높이를 가진 촬영 각도의 데이터에서 이상치를 제거한 Table 2E와 이상치를 고려하지 않은 Table 2D를 비교하였다. Table 2E는 모든 지표에서 가장 낮은 변동성과 오차를 기록하였다. 이는 이상치 제거와 최대 높이를 가진 각도의 데이터를 결합할 때 데이터의 일관성을 극대화하고 오차를 최소화할 수 있음을 보여준다.

최대·최소 벼 품종에서의 높이 데이터 최적화 분석

데이터 가공에 따른 효과는 분얼기 중반인 20일부터 유수형성기와 수잉기를 포함한 30일까지의 단계에서 가장 큰 변화를 보였다(Fig. 3A, B). 데이터 가공에 따라 해당 기간 내 흑진주 품종의 높이 변화량은 53.54 cm에서 30.56 cm로, 만미 품종은 78.54 cm에서 24.01 cm로 크게 감소하였다. 최적화 이전에는 최대 품종인 흑진주의 높이가 24.51~29.87 cm인 반면, 최소 품종인 만미의 높이가 25.22~33.08 cm으로 나타나 두 품종 간 높이에서 역전 현상이 발생하였다. 그러나 최적화 후에는 흑진주의 높이가 31.72~34.77 cm으로 증가하고, 만미 품종의 높이는 26.56~28.96 cm으로 조정되어 품종 간 높이가 조정되었다. 이는 영양생장과 생식생장이 연결되면서 형태적 변화와 생리적 변화가 동시에 일어나, 벼 높이 측정의 불안정성에 따라 오차가 발생한 것으로 생각된다(Langstroff et al. 2022). 30일차 사진에서 볼 수 있듯이, 조생종인 흑진주와 중만생종인 만미의 성장 차이는 두드러지며, 분얼의 퍼짐과 기울어짐에 따라 높이 기준에 대한 유연한 접근이 필요함을 보여준다(Fig. 3C). 이를 통해 데이터 가공이 높이 변화량을 감소시켜 자연스러운 성장곡선의 기울기를 제공할 뿐만 아니라 품종간 높이 비교의 안정성에 기여함을 확인하였다.

Figure 3. Effect of data optimization on rice plant height recorded by PlantCV for the Heukjinju vs. Manmi cultivars. (A) Rice plant height growth curve for a randomly selected sample from the 96 cultivars, (B) Rice plant height growth curve after applying the optimized analysis method (IQR filter applied to 96 cultivars with four replicates and maximum angle selection), and (C) Comparison of time series image for Heukjinju and Manmi cultivars at 10-day intervals. Red lines and arrows represent the cultivar with the maximum height (Heukjinju), while blue lines and arrows represent the cultivar with the minimum height (Manmi). The red rectangle highlights examples of changes in the recorded growth curves as a result of the analysis method applied

SNPs 검출을 위한 GBS분석

시퀀싱 및 GBS 분석의 요약 통계는 Table 3을 통해 볼 수 있다. Paired-end sequencing을 통해 약 221 Gbp(1,467,143,186 reads)가 생성되었다. 어댑터와 낮은 품질의 염기를 제거한 후, Phred quality score가 20이상인 고품질 리드는 136 Gbp (160,652,223,454 reads)를 확보하였다. 767,699,496 reads가 Reference-IRGSP-1.0 reference genome에 매핑되었으며, 이는 library에서 전처리된 된 리드의 약 55.4%에 해당한다. GBS 데이터 세트에서 매핑된 리드의 평균 깊이는 30.31X였다. 변이 검출은 DeepVariant(Poplin et al. 2018)와 Samtools(Danecek et al. 2021)를 사용하여 보다 높은 정확도를 보인 DeepVariant의 결과를 이용하였다. 이는 삽입과 결실(Insertion & Deletion), 구조변이(Structural Variation)에 대한 Samtools의 약점을 합성곱 신경망(Convolutional Neural Network, CNN)기반 Deep Variant가 보완한 것으로 파악하였다. 이를 통해 총 145,882개의 SNP들 중에서 최소 read depth가 4 이상, 소수 대립 유전자 빈도가 5% 이상, 결측치 30% 미만을 만족하는 12,127개의 SNP를 검출하였다.

Table 3 . Descriptive statistics of Genotyping-by-Sequencing and SNP filtering using DeepVariant.

NO.	Process	Value
1	Total number of reads	1,467,143,186
2	Total length of reads (bp)	221,538,621,086
3	Number of trimmed reads	1,367,129,972
4	Total length of trimmed reads (bp)	160,652,223,454
5	Total number of SNPs	145,882
6	Number of SNPs with MAF > 5%	65,872
7	Number of SNPs with missing data < 30%	35,660
8	Final number of SNPs used for genotyping	12,127

최적화된 표현형 데이터를 활용한 GWAS 결과

최종 필터링된 12,127개의 SNP와 표현형 측정 도구 별 비교에서 높은 안정성을 보인 PlantCV 기반 최적화된 표현형 가공 방법 활용하여 벼 최대 높이 형질에 대한 전장유전체 연관 분석(GWAS)을 수행하였다. 벼 최대 높이에 대한 연관 분석 결과, 유의미한 SNP들의 -logP 값을 바탕으로 Manhattan과 Q-Q plot을 작성하였다(Fig. 4).

Figure 4. Results of Genome-Wide Association Study using optimized phenotypic data in rice plant height. (A) Manhattan plot and (B) Q-Q plot. Candidate loci were selected based on a Bonferroni cutoff of 0.01

Manhattan plot에서 Bonferroni 보정 cutoff 0.01 기준인 -logP > 6.08을 초과하는 SNP는 염색체 6번과 9번에서 유의미한 후보 유전자좌인 것으로 확인할 수 있었다. Q-Q plot에서는 관찰된 -logP와 예측된 -logP가 전반적으로 유사하게 나타났으며, FarmCPU와 BLINK모델에서 꼬리 부분의 분포가 관찰된 -logP를 벗어나 유의미한 SNP가 확인되었다(Tsai et al. 2020). 사전 연구된 문헌을 바탕으로 이들 SNP의 양측 200 kb 영역으로부터 성장과 관련된 후보 유전자를 탐색하였다(Table 4). 후보 SNP는 4개였으며(SNP038513K, SNP031426K, SNP068282K, SNP094239K), 이와 관련된 유전자 심볼과 유전자명은 각각 OsSTRL2(Os03g0263600), OsCCS52A(Os03g0123300), OsGH9B3(Os06g0256900), OsBIG(Os09g0247700)로 확인하였다.

Table 4 . Candidate genes for rice plant height obtained from previous studies.

SNP marker	Locus ID	Gene Symbol	Gene Description	Trait Ontology
SNP038513K	Os03g0263600	OsSTRL2	Atypical strictosidine synthase. Involved in regulation of anther development and pollen wall formation (Os03t0263600-01)	male sterility (TO:0000437)heat tolerance (TO:0000259)
SNP031426K	Os03g0123300	OsCCS52A	Activator of the anaphase- promoting complex/cyclosome (APC/C) complex. Involved in regulation of the cell cycle, regulation of postembryonic shoot branching and tillering, and maintenance of endoreduplication during endosperm development (Os03t0123300-01).Similar to cell cycle switch protein CCS52A (Os03t0123300-02)	abscisic acid sensitivity (TO:0000615)gibberellic acid sensitivity (TO:0000166)mitotic cell cycle trait (TO:0000730)tiller number (TO:0000346)root development trait (TO:0000656)panicle number (TO:0000152)root meristem development (TO:0002692)tillering ability (TO:0000329)
SNP068282K	Os06g0256900	OsGH9B3	Glycoside hydrolase family 9 subclass B3, endo-beta-1,4-glucanase. Involved in lignocellulose crystallinity modification in stem internode growth and development, plant strength, cellulose modification, and biomass saccharification (Os06t0256900-01)	lodging incidence (TO:0000068)
SNP094239K	Os09g0247700	OsBIG	Homolog of Arabidopsis BIG. Essential for normal growth and development (Os09t0247700-01)	chlorophyll ratio (TO:0000298)chlorophyll content (TO:0000495)growth and development trait (TO:0000357)root number (TO:0000084)carotenoid content (TO:0000496)

벼의 식물체 높이와 관련된 유전자 연구는 염색체 1번의 SD1(Os01g0883800)와OsPG3(Os01g0637500), 염색체 4번의 OsBRI1(Os01g0718300), 염색체 8번의 OsMADS57(Os02g0731200), 염색체 9번의 OsBIG(Os09g0247700)를 대상으로 이루어져 왔으며, 이외에도 다양한 유전자들이 적용되는 복합유전형질로 알려져 있다(Peng et al. 2023; Zeng et al. 2019). 이번 연구에서도 OsBIG(Os09g0247700)이 검출되었으며, 왜성 유전자로 널리 알려진 SD1(Os01g0883800)과 유사한 기능을 수행하는 OsCCS52A(TAD1; Tillering And Dwarf1, Os03g0123300)도 확인되어 기존 연구 결과와 일부 일치하는 양상을 보였다(Cheng et al. 2018; Xu et al. 2012).

한편, 벼의 전 생육주기 중 최대 성장 높이를 형질 지표로 활용한 차별점으로 인해 추가적인 후보 유전자좌가 검출되었다. OsCCS52A(Os03g0123300)와 OsGH9B3(Os06g0256900)는 각각 분얼(Tillering)과 도복 저항성(Lodging incidence)과 관련된 유전자로, 최대 성장 높이에서 얻은 지표가 단순한 식물체 높이를 넘어, 줄기의 세포벽 성장 능력을 반영할 수 있음을 보여준다(Huang et al. 2019). 더 나아가, 식물체의 높이는 지베렐린(GA) 및 브라시노스테로이드(BR) 생합성 경로에 의해 결정되므로, 이는 성장 호르몬 네트워크와 밀접한 관계가 있는 분얼과 도복 발생률이 중요한 간접 지표로 기능할 가능성을 시사한다(Niu et al. 2021). 불임성 관련 유전자인 OsSTRL2(Os03g0263600) 또한 성장 호르몬 생합성과 연관되며, 영양 생장뿐만 아니라 생식 생장에까지 관여할 수 있다는 기능적 연관성이 연구된 바 있다(Kandpal et al. 2020; Zou et al. 2017). 이를 통해 표현체 관점의 형질 정의가 2차 지표와 연결 될 수 있음을 보여준다.

이번 연구를 통해 충분한 식물 도메인 지식과 PlantCV 기반의 전문적인 데이터 가공 파이프라인이 형질 관련 유전자 연구에 효과적으로 기여할 수 있음을 확인하였다. 촬영된 이미지 기반 높이 형질에서 분얼과 도복 발생률과 같은 간접 표현형 지표를 확인함에 따라, 표현체 관점의 형질 정의가 복합형질에 대한 해석을 도울 수 있을 것으로 기대된다. 향후 다양한 형질에서 유전체 분석과의 연계를 통해 표현형 시계열 데이터의 활용이 더욱 확장될 것이다.

본 연구는 벼 96품종의 전 생육주기 동안 촬영된 시계열 데이터를 바탕으로, 최적화된 표현체 분석 도구 및 데이터 처리에 대한 최적화와 그 유효성 검증에 중점을 두었다. 세 가지 주요 표현형 분석 도구인 ImageJ, OpenCV, PlantCV를 활용하여 벼의 키 성장 곡선을 비교한 결과, 도구마다 생장 단계별로 고유한 성능 차이가 나타났다. ImageJ는 초기 생장 단계에서 변동성이 컸으며, OpenCV는 후기 단계에서 정확도가 저하되는 경향이 있었다. 반면, PlantCV는 모든 생장 단계에서 안정적이고 일관된 결과를 보여, 표현형 분석 도구로서의 높은 신뢰성을 입증하였다. 추가적으로, 반복 샘플, 촬영 각도 선택, 이상치 제거가 데이터 변동성과 오류율에 미치는 영향을 분석한 결과, 반복 샘플만으로는 변동성 제어가 충분하지 않았으나, 각도 최적화와 이상치 제거를 결합한 방법이 데이터의 신뢰성과 정확성을 크게 개선하는 데 기여함을 확인하였다. 특히, 이상치 제거와 각도별 최대 값 선별은 성장 곡선을 보다 부드럽고 안정적으로 조정하는 데 중요한 역할을 하였다. 이러한 최적화된 분석 방법을 통해 도출된 형질 데이터는 딥러닝 기반 GBS(Genotyping-by-Sequencing) 파이프라인으로 식별된 12,127개의 SNP을 활용한 GWAS(유전체 연관 분석)를 통해 검증되었으며, 염색체 3, 6, 9번에서 벼 키와 연관된 유의미한 SNP를 확인하였다. 주요 유전자 OsBIG, OsGH9B3, OsSTRL2, OsCCS52A는 성장 호르몬 생합성 경로와 연결되어 본 연구의 최적화된 표현형 분석 방법이 유의미한 형질-연관 마커 식별에 효과적임을 입증하였다.

본 연구는 농촌진흥청 연구사업(과제번호: RS-2023-00230677) “작물 표현체 데이터 표준화 및 수집체계 구축” 과제의 지원에 의해 이루어졌다.