Major QTLs, qARO1 and qARO9, Additively Regulate Adaxial Leaf Rolling in Rice

Moderate leaf rolling is considered optimal for the ideal plant type in rice (Oryza sativa L.), as it improves photosynthetic efficiency and, consequently, grain yield. Determining the genetic basis of leaf rolling via the identification of quantitative trait loci (QTLs) could facilitate the development of high-yielding varieties. In this study, we identified three stable rice QTLs, qARO1, qARO5, and qARO9, which control adaxial leaf rolling in a recombinant inbred line (RIL) population derived from a cross between Tong 88-7 (T887) and Milyang 23 (M23), using high-density SNP markers. These QTLs controlled the rolling phenotype of both the flag leaf (FL) and secondary leaf (SL), and different allelic combinations of these QTLs led to a wide variation in the degree of leaf rolling. Additive gene actions of qARO1 and qARO9 on leaf rolling were observed in a backcross population. In addition, qARO1 (markers: 01id4854718 and 01asp4916781) and qARO9 (markers: 09id19650402 and 09id19740436) were successfully fine-mapped to approximately 60- and 90-kb intervals on chromosomes 1 and 9, respectively. Histological analysis of near-isogenic lines (NILs) revealed that qARO1 influences leaf thickness across the small vein, and qARO9 affects leaf thickness in the entire leaf and bulliform cell area, thus leading to adaxial leaf rolling. The results of this study advance our understanding of the genetic and molecular bases of adaxial leaf rolling, and this information can be used for the development of rice varieties with the ideal plant type.


INTRODUCTION
Leaf is the major organ responsible for photosynthesis and transpiration. In rice (Oryza sativa L.), leaf shape is one of the main factors affecting the rice plant type. Thus, it is considered an important agronomic trait determining photosynthetic potential and grain yield in rice (Smith and Dilday, 2002;Yang and Hwa, 2008). The V-shape or moderate rolling of leaves increases canopy photosynthesis by enhancing CO 2 penetration (Chakraborty, 2001) and improves light capture by reducing mutual shading (Defeng et al., 2001). Therefore, moderately rolled leaves are thought to support higher grain yield than flat leaves (Setter et al., 1995). Moreover, under stress conditions such as high light intensity and excessive transpiration, leaf rolling acts as a stress avoidance mechanism by reducing the leaf surface area (Kadioglu and Terzi, 2007). Therefore, moderate leaf rolling is considered a key trait of high-yielding varieties with the ideal plant type, and this trait has been widely used in the breeding of super hybrid rice cultivars (Denning and Mew, 1997).
Leaf morphology is regulated by a range of cellular processes, including cell development, differentiation, specification, and axis determination (Itoh et al., 2005;Zou et al., 2014). The main factors determining leaf morphology include bulliform cells, sclerenchyma cells, adaxial-abaxial axis polarity, and cuticle development (Xu et al., 2018). To date, several genes have been isolated that regulate leaf rolling by altering these histological characteristics of leaf. Semi rolled leaf 1 (SRL1) regulates leaf rolling by inhibiting the formation of bulliform cells. Rice srl1 mutants show an increased number of bulliform cells on the adaxial surface of the leaf blade, causing adaxial rolling (Xiang et al., 2012). By contrast, mutation of the narrow leaf 7 (NAL7) gene decreases the size but increases the number of bulliform cells, leading to adaxial leaf rolling (Fujino et al., 2008). Semi rolled leaf 2 (SRL2) is involved in the regulation of cell differentiation on the abaxial side (Liu et al., 2016). The srl2 mutants exhibit abnormal development of sclerenchyma cells on the abaxial surface of leaves, leading to adaxial leaf rolling. Adaxialized leaf 1 (ADL1) plays an important role in leaf pattern formation (Hibara et al., 2009). adl1 mutants exhibit abaxial leaf rolling due to defects in epidermal development, resulting in abnormal adaxial-abaxial polarity. OsMYB103L influences the cellulose content and mechanical strength of leaves (Yang et al., 2014) and its overexpression increases the cellulose content of leaves, resulting in adaxially rolled leaves.
Despite the identification of numerous genes controlling leaf morphology, not all genes are suitable for improving the plant type via ideotype breeding since some of these genes cause extreme leaf rolling or are associated with unfavorable effects, such as growth retardation, abnormal organ development, and reduced grain yield (Fujino et al., 2008;Hu et al., 2010;Yang et al., 2014;Zhang et al., 2015;Liu et al., 2016;Li et al., 2017). Recently, quantitative trait loci (QTLs) and their alleles promoting moderate leaf rolling have been identified in natural populations of rice using joint linkage-association mapping (Zhang et al., 2016). Six QTLs were simultaneously identified by genome-wide association study (GWAS) and linkage mapping. These QTLs could explain 3.3-12.9% of the variation in the leaf rolling index (LRI) of 262 recombinant inbred lines (RILs), and the favorable alleles of these QTLs increased the LRI of 1,129 rice accessions by 1.4-2.4%. However, the application of these QTLs for improving the leaf rolling phenotype of rice via breeding remains unknown. To improve the plant type by controlling leaf rolling via molecular breeding, it is necessary to verify that the genetic effects of these QTLs or their favorable alleles will be expressed consistently in recipient plants, and the simultaneous use of multiple QTLs will allow the precise modulation of the leaf rolling phenotype within a wide range.
Milyang 23 (M23), a Korean high-yielding variety, was developed by performing an inter-subspecific cross between indica and japonica. Since M23 exhibits some desirable traits contributing to high yield, such as erect leaves, short culm length, large number of spikelets per panicle, and increased translocation of non-structural carbohydrate to sink tissues (Nuruzzaman et al., 2000;Kobayashi et al., 2003;Takai et al., 2005), it has been widely used to develop a high-yielding variety as an elite donor (Kaneda, 1986;Lamo et al., 2015). In particular, the erect plant type of M23, which is attributed to its long-wide, straight, and adaxially rolled leaves, contributes to a high canopy photosynthetic rate and high biomass production (Kobayashi et al., 2003;San et al., 2018).
The objectives of this study are (1) to identify the QTLs affecting the adaxial leaf rolling of M23 using a RIL population, (2) to validate the genetic effect of stable QTLs on the degree of leaf rolling in backcross progenies, and (3) to reveal the morphological basis of adaxial leaf rolling by examining the histological characteristics of leaves altered by these QTLs.

Plant Materials for QTL Mapping and Field Testing Conditions
A total of 162 F 15 RILs derived from an inter-subspecific cross between Tong 88-7 (T887; japonica rice; flat leaf) and Milyang 23 (M23; indica rice; moderate rolled leaf) were used for genotyping-by-sequencing (GBS) and QTL analysis. All plant materials were cultivated in the experimental field of Seoul National University, Suwon, South Korea. Field tests (FTs) were conducted during May-October, 2015 (field test 1; FT 1) and 2016 (field test 2; FT 2). Detailed information of the environmental conditions is shown in Supplementary Table 1. Thirty-dayold seedlings were transplanted into the paddy field as follows: one line per row, 25 plants per row, one plant per hill, 15 cm spacing between two plants in a row, and 30 cm spacing between adjacent rows.

Adaxial Leaf Rolling Evaluation
Adaxial leaf rolling was assessed based on the LRI, which was calculated by the following equation (Shi et al., 2007): Natural distance of leaf blade margin)/ Width of fully expanded leaf blade × 100.
us, a higher LRI means more leaf rolling and vice versa. LRI was measured from the center of the flag leaf (FL) and upper secondary leaf (SL) for individual plants, after the neck of the panicle had completely emerged from the leaf sheath. The uppermost stem was not used for measurements because of large variation in leaf size. To avoid the effect of changes in soil water content on leaf rolling, the rice plants were constantly irrigated until the measurements were complete. All statistical analyses of phenotypic values were performed using R studio v1.2.5033 (Allaire, 2012).

Histological Assay
A 2-cm section of the SL was sampled at the center of the leaf blade and immediately fixed in formalin-acetic acid-alcohol (FAA) solution (10% formalin, 50% ethanol, 5% acetic acid, and 35% water). All samples were stored at 4 • C until needed for histological analysis. Transverse sections of the leaf blade were prepared manually using a razor blade. The leaf sections were cleaned several times with distilled water and then mounted on a microscope slide with a drop of water. Histological characteristics of each leaf were observed under the CX31 light microscope (Olympus), and images were captured using an eXcope T500 microscope camera (DIXI Science). The ImageJ software (version 1.51; Schneider et al., 2012) was used to measure various histological characteristics from the images. Large vein (LV) number, small vein (SV) number, and SV:LV ratio was determined in the entire leaf. Leaf thickness parameters were measured across the SVs and bulliform cells (BCs). Leaf thickness across the SV (LTSV) was defined as the distance from the most protruding region of the epidermis on the adaxial surface to the most recessed region of the epidermis on the abaxial surface ( Figure 5A). Leaf thickness across the BC (LTBC) was defined as the distance from the most recessed point of the BC on the adaxial surface to the most protruding region of the epidermis on the abaxial surface. BC height was defined as length of the longest longitudinal axis within the BC. Values of these characteristics were used to calculate the ratio of LTBC:LTSV and the proportion of BCs (as shown below): BC proportion (%) = BC height/LTBC × 100.

GBS and QTL Analysis
Genomic DNA was extracted from young leaves of each RIL using the cetyl trimethylammonium bromide (CTAB) method (Murray and Thompson, 1980), and the quality of the isolated DNA was verified using PicoGreen (Invitrogen). GBS libraries of RILs were constructed as described previously (Elshire et al., 2011). Briefly, DNA (10 ng/µl) was digested with ApeKI and ligated to adapters on both ends. The quantity and quality of GBS libraries were determined using the Bioanalyzer Kit (Agilent Genomics), and the libraries were sequenced on the HiSeq 2000 platform. The reference genome sequence was in silicosearched for ApeKI restriction sites and split into fragments. A subset of the Nipponbare reference genome IRGSP 1.0 (Sakai et al., 2013) with ApeKI restriction was prepared for variant calling. A subset of genome fragments smaller than 2-kb were collected and indexed, and the GBS reads were mapped using the Burrows-Wheeler alignment (BWA) tool (Li and Durbin, 2009). The genotype of each RIL was determined using SAMtools, with default parameters . Raw GBS data of 155 out of 162 RILs were deposited in the National Center for Biotechnology Information (NCBI) Short Read Archive (SRA) under the bio project number PRJNA601019 in previous study (Jang et al., 2020). A total of 6,141 markers (missing rate < 30%) were used for the construction of a genetic map, and QTL analysis was performed using ICIMapping 4.1 (Meng et al., 2015). The recombination distance was calculated using the Kosambi mapping function (Kosambi, 1943). The flag leaf rolling index (FLRI) and secondary leaf rolling index (SLRI) data collected from two field tests, FT 1 and FT 2, were analyzed by inclusive composite interval mapping (ICIM) to detect significant QTLs with a logarithm of odds (LOD) score > 2.5.

Marker Development
To develop molecular markers, polymorphisms between the two parental accessions, T887 and M23, were identified by comparing their whole genome sequence (WGS) data generated on the Illumina HiSeq platform. The WGS data of T887 and M23 were deposited in the NCBI SRA database under the bio sample numbers SAMN13840663 and SAMN03120572, respectively. All gel-based insertion-deletion (InDel) and single-nucleotide polymorphism (SNP) markers were developed using Primer3Plus (Untergasser et al., 2007) and BatchPrimer3 (You et al., 2008;Supplementary

Validation and Fine Mapping of QTLs
To examine the gene action of two QTLs, qARO1 and qARO9, a BC 2 F 2 population was developed by backcrossing T887 (donor parent) with M23 (recurrent parent). Backcross progenies harboring the qARO5 allele of M23 were selected using three markers. qARO1 and qARO9 alleles of T887 were identified using seven and five markers, respectively (Supplementary Table 2). Near-isogenic lines (NILs) carrying both qARO1 and qARO9 or only qARO9 was selected from each individual BC 3 F 2 plant with the genetic background recovered to the recurrent parent. Marker-assisted background selection for developing NILs was conducted using 186 markers, including 178 SNPs (Fluidigm platform), which evenly cover the rice genome (Seo et al., 2020), and eight InDel markers (Supplementary Table 2). The BC 3 F 3 population was employed for fine mapping of qARO1 and qARO9 using 1,860 and 368 individuals, respectively.

Variation in the Leaf Rolling Phenotype of RILs
The high-yielding Korean rice variety M23 showed a straighter plant type, with erect leaves, than T887. Leaves of M23 showed moderate rolling, whereas those of T887 were relatively flat ( Figure 1A). In the two FTs, the FLRI and SLRI of T887 ranged from 10.3 to 14.1% and 3.1 to 5.4%, respectively, whereas those of M23 ranged from 23.8 to 24.1% and 26.5 to 29.3%, respectively ( Figure 1B). Among the RILs derived from a cross between T887 and M23, the value of FLRI varied from 0%  Table 3). A strong positive correlation (r = 0.9) was detected between the FLRI and SLRI values of RILs each year, as well as between the two FTs for FLRI (r = 0.9) and SLRI (r = 0.85) ( Figure 1C).  Table 4). A total of 12 and 10 significant QTLs were detected for FLRI and SLRI, respectively, in both FTs (Figure 2A and Supplementary  Table 5). Among these, three QTLs were simultaneously detected for FLRI and SLRI in the two FTs. These three stable QTLs conferring adaxial leaf rolling were named adaxially rolled leaf 1 (qARO1), qARO5, and qARO9; these three QTLs were detected on the short arm of chromosome 1 (4.05-4.67 Mb), the long arm of chromosome 5 (20.41-23.05 Mb), and the long arm of chromosome 9 (19.29-19.7 Mb), respectively (Figure 2A and Supplementary  (Figure 2A and Supplementary Table 5). Among these three QTLs, qARO1 and qARO9 showed positive additive effects, implying that M23 alleles at both these QTLs (qARO1 M23 and qARO9 M23 , respectively) contribute to the leaf rolling phenotype. By contrast, qARO5 showed a negative additive effect, and the T887 allele at qARO5 (qARO5 T887 ) increased the LRI (Figure 2B and Supplementary Table 5). RILs containing a combination of qARO1 M23 , qARO5 T887 , and qARO9 M23 alleles showed the highest average FLRI and SLRI values in both FTs, whereas those harboring the qARO1 T887 , qARO5 M23 , and qARO9 T887 alleles showed the lowest FLRI and SLRI values, except for SLRI in FT 2 ( Figure 2C).

QTL Analysis of the LRI of RILs
Gene Action of qARO1 and qARO9 in the BC 2 F 2 Population Since RILs do not contain heterozygous alleles of QTLs, the gene action of qARO1 and qARO9 was examined in the BC 2 F 2 population, which was developed using M23 as the recurrent parent. Backcross progenies of plants harboring the donor allele at qARO5 (qARO5 M23 ) and heterozygous alleles at qARO1 and qARO9 were selected ( Figure 3A). The qARO1 M23 and qARO9 M23 alleles conferred higher LRI in the BC 2 F 2 population, consistent with the results of primary mapping (Figure 2A).
The LRI values of plants harboring the heterozygous alleles at qARO1 and qARO9 were intermediate between the two parental LRI values, indicating that qARO1 and qARO9 exhibit additive gene actions (Figures 3B,C). Furthermore, both QTLs showed additive gene action, regardless of the allele type of the other QTL, implying the absence of epistatic interactions between two QTLs (Figures 3B,C).

Validation of the Effects of qARO1 and qARO9 Using NILs
Each NIL was developed by backcrossing with both parents and then was genotyped for qARO1 (markers: 01id4052916 and 01id5210929) and qARO9 (markers: 09id18904919 and 09id19950022). Each NIL was selected from the individual plant with the genetic background recovered maximally to the recipient parent, carrying qARO1 and/or qARO9 alleles from the donor parent. Consequently, four NILs including T887-qARO9 M23 , T887-qARO1 M23 +ARO9 M23 , M23-qARO9 T887 , and M23-qARO1 T887 +qARO9 T887 were successfully developed. Among these, T887-qARO9 M23 and T887-qARO1 M23 +ARO9 M23 recovered 91.4% (170 out of 186 markers) and 87.6% (163 out of 186 markers) of the T887 genetic background, while M23-qARO9 T887 and M23-qARO1 T887 +qARO9 T887 recovered 97.3% (181 out of 186 markers) and 93.6% (174 out of 186 markers) of the M23 background ( Figure 4A). Both qARO1 and qARO9 showed additive effects on FLRI and SLRI, regardless of the genetic background. Significant differences were detected in FLRI and SLRI among NILs and their parent, but no significant changes were observed in leaf width, panicle number, panicle dry weight, and plant dry weight (Figures 4B-D

Histological Analysis of NILs
To determine the effects of qARO1 and qARO9 on adaxial leaf rolling, histological assays were conducted ( Figure 5A). The LV number, SV number, and SV:LV ratio in both FL and SL showed no significant differences between NILs and their parents, suggesting that vein number is not involved in the regulation of the leaf rolling phenotype (Figure 4B). More precisely, five characteristics including LTBC, LTSV, LTBC:LTSV ratio, BC height, and BC proportion were compared among NILs with the T887 genetic background and both parents (Figures 5B,C). Among these five characteristics, LTBC, BC height, and BC proportion showed significant differences between T887 and two NILs but not between T887-qARO9 M23 and T887-qARO1 M23 +ARO9 M23 . LTSV was significantly different not only between T887 and two NILs but also between T887-qARO9 M23 and T887-qARO1 M23 +ARO9 M23 . The LTBC:LTSV ratio showed significant differences between  T887-qARO9 M23 and T887-qARO1 M23 +ARO9 M23 but not between T887 and T887-qARO9 M23 (Figure 5C).

Fine Mapping of qARO1 and qARO9 Using the BC 3 F 3 Population
To further narrow down the region of chromosomes 1 and 9 where qARO1 and qARO9 are located, we developed a BC 3 F 3 population comprising 2,228 plants by backcrossing in both parental directions, T887 (924 plants) and M23 (1,304 plants). The qARO9 region was delimited to an interval of 0.22-Mb on chromosome 9 (19.6-19.82 Mb) between marker 09id19607766 and marker 09id19822650 using 1,860 plants derived from the ancestor plants harboring the recipient allele of qARO1 and qARO5. Six additional markers were used to narrow down the target region in five recombinant plants. Multiple comparisons revealed that recombinants were classified into two groups, recurrent parent and NIL-qARO9, based on the average FLRI and SLRI values. Thus, qARO9 was delimited to a 90-kb interval (19.65-19.74 Mb) (Figure 6A).
qARO1 was fine-mapped using the 368 plants derived from the ancestor plants with the recipient allele of qARO5 and donor allele of qARO9 to clearly discriminate the qARO1 effect on LRI. Consequently, qARO1 was mapped to a 0.13-Mb interval (4.82-4.95 Mb) on chromosome 1 between markers 01id4814665 and 01id4946123 using the entire population. A total of six recombinants were identified within this interval. Using three additional markers, qARO1 was delimited to a 60-kb interval (4.86-4.92 Mb) (Figure 6B).   qARO1 (B). Red, blue, and gray colors denote T887, M23, and heterozygous genotype. Values of FLRI and SLRI were compared using one-way ANOVA, followed by Scheffe's post hoc test (p < 0.05).
Frontiers in Plant Science | www.frontiersin.org

DISCUSSION
Moderate leaf rolling is considered a desired leaf trait for high yielding varieties of rice with the ideal plant type (Setter et al., 1995;Denning and Mew, 1997;Chakraborty, 2001;Defeng et al., 2001). Therefore, identification of loci that modulate leaf rolling and dissection of the genetic effects of these loci could facilitate the development of cultivars with the ideal plant type. In this study, we performed QTL mapping to detect stable QTLs that regulate adaxial leaf rolling in rice using a RIL population. Additionally, to validate the genetic effects of QTLs on adaxial leaf rolling, we investigated the morphological and histological characteristics of leaves using backcross progenies.
In the RIL population derived from a cross between T887 and M23, the FLRI and SLRI values showed transgressive segregation in the positive direction, implying that the genetic factors responsible for adaxial leaf rolling could have originated from either parent ( Figure 1B). Among the three stable QTLs that simultaneously affected FLRI and SLRI in both FTs, qARO1 and qARO9 showed positive additive effects on LRI, while qARO5 showed a negative additive effect, suggesting that alleles causing adaxial leaf rolling were derived from each parent belonging to a different subspecies of rice, japonica and indica. A similar result was reported previously (Zhang et al., 2016). The authors reported highly transgressive variation in the LRI trait in the RIL population derived from a cross between Minghui 63 (indica rice) and 02428 (japonica rice). Moreover, among six stable QTLs affecting LRI, the positive alleles of four loci were derived from 02428, while positive alleles of the other QTLs were derived from Minghui 63. Taken together, these results suggest that alleles causing adaxial leaf rolling might be distributed in different subspecies and that introducing allelic combinations of these QTLs by inter-subspecific crossing is an efficient strategy for modulating the degree of leaf rolling within a diverse range.
Several previous studies have reported several QTLs for leaf rolling on chromosomes 1, 5, and 9. Based on the IRGSP 1.0 database (Sakai et al., 2013), qRl5-4 is located at the 20.51-to 21.99-Mb position on chromosome 5 and qRl9-1 is located at the 7.96-to 10.90-Mb position on chromosome 9 (Zhang et al., 2016). In addition, qRL-1, represented by the nearest marker RM3453, is located near the 4.89-Mb position on chromosome 1 and qRL-9, represented by the nearest marker RM8206, is located near the 5.92-Mb position on chromosome 9 . Among these QTLs, qRL-1 and qRl5-4 were co-located with qARO1 (4.86-4.92 Mb) and qARO5 (20.41-23.05 Mb), respectively ( Figure 6B and Supplementary Table 5), while qARO9 was newly detected in this study.
A relatively wide range of FLRI (3.3-45.9%) and SLRI (4.6-39.9%) was observed in the RIL population, depending on the allelic combination of qARO1, qARO5, and qARO9. These three QTLs explained on average 45 and 43.4% of the phenotypic variance in FLRI and SLRI, respectively ( Figure 2D and Supplementary Table 5). In addition, accumulation of the genetic effects of qARO1 and qARO9 in NILs changed the leaf shape from almost flat to moderately rolled ( Figure 4B). These results suggest that qARO1, qARO5, and qARO9 precisely regulate the leaf rolling phenotype within a wide range and therefore are ideal candidates for application in molecular breeding programs.
Although several genes controlling the leaf rolling phenotype have been cloned via genetic mapping of mutant plants (Fujino et al., 2008;Hibara et al., 2009;Xiang et al., 2012;Yang et al., 2014;Liu et al., 2016), these mutant genes exhibit several limitations that prohibit their use in breeding programs. For example, mutant genes cause severe leaf rolling, which leads to developmental defects and grain yield reduction (Li et al., 2017). Moreover, mutant genes often cause unfavorable effects on other agronomic traits. Almost all leaf rolling mutants show recessive inheritance; thus, both male and female parents must possess the same recessive gene to modulate leaf rolling in super hybrid rice, the development of which is time-consuming and labor-intensive . Since genotypes heterozygous at qARO1 and qARO9 showed intermediate values of LRI compared with the parents (Figure 3A), and regulated adaxial leaf rolling without significant defects in leaf width, panicle biomass, or plant biomass of individual plants ( Figure 4B and Supplementary Table 6), these QTLs could be used to improve the leaf shape of rice cultivars via ideotype breeding. However, further investigation is needed to verify whether moderate leaf rolling resulting from the QTLs could increase grain yield at the population level.
Leaf rolling is mainly caused by histological alterations in leaf cells (Xu et al., 2018). Both qARO1 and qARO9 regulated adaxial leaf rolling by affecting different histological characteristics of leaves. The T887-qARO9 M23 NIL showed significantly lower leaf thickness parameters (LTBC and LTSV) but higher BC height and BC proportion than T887 ( Figure 5C). The T887-qARO1 M23 +qARO9 M23 NIL showed further reduction in LTSV, resulting in higher LTBC:LTSV ratio and greater leaf rolling compared with T887-qARO9 M23 ( Figure 5C). Overall, our data suggest that qARO9 is involved in the determination of thickness in the entire leaf, as well as BC size, whereas qARO1 is involved in the determination of the ratio of leaf thickness by regulating leaf thickness near the SV. This implies that qARO1and qARO9-mediated reduction in leaf thickness in specific regions of the leaf and an increase in the proportion of BC alter the leaf structure, resulting in adaxial leaf rolling.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
SJ designed the experiment and prepared the manuscript. SJ and DL conducted the field experiments. SJ and YL conducted the wet-laboratory experiments. SS conducted computational analysis. H-JK supervised the work and contributed to the finalization of the manuscript. All authors read and approved the manuscript.

FUNDING
This work was carried out with the support of "Cooperative Research Program for Agriculture Science and Technology Development (Project No. PJ01572901)" Rural Development Administration, Republic of Korea.