QTL mapping and candidate gene mining of flag leaf size traits in Japonica rice based on linkage mapping and genome-wide association study

As one of the most important factors of the japonica rice plant, leaf shape affects the photosynthesis and carbohydrate accumulation directly. Mining and using new leaf shape related genes/QTLs can further enrich the theory of molecular breeding and accelerate the breeding process of japonica rice. In the present study, 2 RILs and a natural population with 295 japonica rice varieties were used to map QTLs for flag leaf length (FL), flag leaf width (FW) and flag leaf area (FLA) by linkage analysis and genome-wide association study (GWAS) throughout 2 years. A total of 64 QTLs were detected by 2 ways, and pleiotropic QTLs qFL2 (Chr2_33,332,579) and qFL10 (Chr10_10,107,835; Chr10_10,230,100) consisted of overlapping QTLs mapped by linkage analysis and GWAS throughout the 2 years were identified. The candidate genes LOC_Os02g54254, LOC_Os02g54550, LOC_Os10g20160, LOC_Os10g20240, LOC_Os10g20260 were obtained, filtered by linkage disequilibrium (LD), and haplotype analysis. LOC_Os10g20160 (SD-RLK-45) showed outstanding characteristics in quantitative real-time PCR (qRT-PCR) analysis in leaf development period, belongs to S-domain receptor-like protein kinases gene and probably to be a main gene regulating flag leaf width of japonica rice. The results of this study provide valuable resources for mining the main genes/QTLs of japonica rice leaf development and molecular breeding of japonica rice ideal leaf shape.


Introduction
Rice is one of the most important food crops in the word, bearing the lifeline of human food security [1][2][3]. The rice plant leaf is important part of histogenesis and morphogenesis, and also the main organ for photosynthesis and respiration [4][5][6]. The photosynthetic energy storage and normal life function of rice are directly affected by the leaf [7]. Meanwhile, leaf size is an important part of ideal rice plant type, as well as an important character of rice yield formation [8]. As the top functional leaf of rice, flag leaf has significant effects on plant related physiological characters and field population structure [9,10]. Therefore, shaping and screening of rice leaf morphology in the breeding process is an effective way to improve rice quality and yield.
Although leaf size has a great influence on high photosynthetic efficiency, the genetic mechanism of leaf morphological characteristics is still unclear [11]. Mining QTLs by linkage mapping in RILs and GWAS mapping in natural populations to find the candidate genes might be the most efficient way to analyze the genetic basis of rice leaf shape.
In recent years, many QTLs and regulatory genes related to rice leaf morphology have been discovered (http:// www. Jiangxu Wang and Tao Wang contributed equally to this work. grame ne. org/) [12][13][14][15]. These QTLs and genes can change the physiological function of plants by regulating leaf morphology, and have important impacts on the coordination of light energy utilization and "sink-source" relationship. At the same time, their potential impact on rice yield has also been gradually discovered. Tang et al. [16] used CSSL population with 143 individuals, obtained 14 leaf length and 19 leaf width QTLs, and further obtained the rice leaf size gene Ghd7.1 by fine mapping, mutation test, and allelic variation analysis. Zhang et al. [17] mapped the flag leaf width QTL qflw7.2 to 27.1 kb by recombinant inbred lines derived from 93 to 11 and Peiai 64 s, and identified 2 candidate genes LOC_Os07g41180, LOC_Os07g41200. A natural population of 532 individuals was used for genome-wide association analysis and high-throughput leaf scoring, 73 QTLs associated with rice leaves were mapped as a result [18]. At present, there is a clear understanding of the cloned gene NAL1 located on chromosome 4, which affects the growth of lateral leaf [19]. The NAL1 mutant is characterized by dwarf and narrow leaf, and the expression level of the genes related to the polar transport of auxin and leaf development in the mutant changed [20]. It was highly expressed in vascular tissue, which played an important role in cell division and cell size regulation, promoting the lateral growth of the leaf.
OsFLW7 regulated the width of flag leaf and increased photosynthetic leaf area [21]. At the same time, OsFLW7 is an allele of GL7 / GW7, which may be related to the regulation of grain traits and the mutant was found to have increased grain length, grain plumpness and yield. GWAS and linkage analysis are both accurate and effective tools for QTL detection of complex rice traits [22]. The breadth and accuracy of QTL detection could be significantly improved by combining these 2 methods. In this study, 2 sets of RILs populations with high-density bin-map and 295 re-sequenced japonica rice accessions were used to conducted linkage/GWAS mapping of flag leaf length, width and area of japonica rice. Two novel QTL qFL2, FL10 and five candidate genes related to flag leaf size in japonica rice were discovered, providing important references for molecular breeding of japonica rice with ideal leaf size and plant type.

Populations for QTL mapping
Natural population is composed of 295 rice varieties, most of which are temperate japonica rice, widely collected from Heilongjiang, Jilin and Liaoning provinces, while foreign varieties mainly originate from Japan, Korea, the Democratic People's Republic of Korea and Russia. This population has been used in previous studies [23]. RILsA contains 195 individuals, obtained from crossing the narrow erect leaf japonica rice variety K131 and wide curved leaf upland rice variety HDB. RILsB contains 189 individuals, which was derived from a cross between the long wide flag leaf japonica rice variety WD20342 and short narrow flag leaf variety Caidao. All materials were planted in Acheng rice experimental base of Northeast Agricultural University from 2019 to 2020. Each individual was planted in 8 rows with 20 plants in each row. Single plant transplanting was implemented and the spacing of rows and plants was 30 cm × 3 cm. The field management of water and fertilizer followed the basic method of conventional field production.

Phenotypic identification of flag leaf size
In order to reduce the influence of marginal effect on phenotype, for each variety plants the 5 th plant in row 4 were selected as the research objects, and calculated that the 5 plants average value of each variety as the phenotypic value and the average value of five plants for each line was calculated. The flag leaf length (FL), width (FW) and area (FLA) were investigated at full heading stage using Tuopu YMJ-D Living Leaf Area Meter (Tuopu Yunnong Technology Co., Ltd). Population phenotypic correlation analysis for QTL mapping was performed by SPSS.

Linkage analysis and genome-wide association study
Two linkage maps for linkage analysis were constructed through 10 K array genotyping by targeted sequencing (GBTS) supported by MOLBREEDING Biotechnology Co., Ltd (Shijiazhuang, China). The RILsA population has been used in a previous study [24]. Nine hundred and seventyeight bin markers covered 2465.32 cM of the rice genome with an average distance of 2.52 cM constructed the linkage map (Fig. S1). The linkage map of RILsB contains 527 bin markers that covered 1874.85 cM of the rice genome with an average distance of 3.56 cM (Fig. S2). The IciMapping Ver.4.2 [25] based on inclusive composite interval mapping (ICIM) was used to detected the QTLs for rice flag leaf traits. The walking speed was set as 1 cm, and the LOD threshold of ICIM was set as 2.5. To ensure the accuracy of mapping results, we controlled the type 1 error of whole genome detection below 5% by 1000 permutation tests. The natural population for GWAS was deep re-sequenced by Beijing Genomics Institute (BGI www. genom ics. org. cn). A total of 788, 396 SNPs meeting the criteria (minimum allele frequency ≥ 5%, deletion rate ≤ 20%) were selected for follow-up analysis, and 295 japonica rice varieties' population structure analysis, genetic relationship analysis and linkage disequilibrium analysis have been completed in a previous study [23]. The mixed linear model (MLM) of TASSEL 5.0 [26] was used for genome-wide association analysis and setting the threshold of SNPs significantly associated with flag leaf traits as 5.46 × 10 -6 which was calculated by GEC software (http:// statg enpro. psych iatry. hku. hk/ gec/). If 2 or more SNPs were located in the same LD interval, they were regarded as the same QTL, and the SNP with the smallest p value was treated as the lead SNP. QQman package in R was used to create the Manhattan and Q-Q plots [27].

Haplotype analysis of candidate genes
Considering the LD decay of the whole genome was confirmed in a previous study by 109 kb [23], we selected 218 kb upstream and downstream of the SNP as the target interval to screen candidate genes. Non-synonymous SNPs of all exons were extracted from the "Rice SNP-Seek Database" of The International Rice Informatics Consortium (IRIC) (https: //snp-seek.irri.org/), which were then used for haplotype analysis of candidate genes with DnaSP software [28].

RNA extraction and qRT-PCR analysis
The flag leaves of four parents (K131, HDB, WD20342, Caidao) of RILs populations took 4-6 days to fully extend. From the flag leaf initial growth to full extension, a total of 4 samples were taken, repeated 5 times for each parent. The total RNA was extracted with TranZol Up RNA Kit (TransGen Biotech). HiFiScript cDNA synthesis kit (CoWin Biosciences, Beijing, China) was used to synthesize cDNA. qRT-PCR was performed on Bio-Rad CFX96 system using 2 × Fast qPCR Master Mixture with 3 biological replicates for each sample. House-keeping gene Actin1 was used to measure the mRNA levels of candidate genes [23] as an internal control. The primers used for qRT-PCR in this study are all shown in Table S1. Relative gene expression levels were determined using the 2 −ΔΔCt method [29]. Data shown in figures and tables are mean values of three replicates.

Phenotypic analysis
The phenotypic data of flag leaf size of RILs populations and natural population from 2019 to 2020 are shown in Tables S2, S3, and the general trends during the 2 years were basically the same. The 3 leaf size traits of parents of RILs populations showed great phenotypic differences during 2 years and the RILs showed significant variation in FL, FW and FLA. FL and FLA in 2 RILs populations with standard deviation from 4.38 to 7.25, presented stronger variation than FW with the standard deviation from 0.17 to 1.44. The variation characteristics of flag leaf phenotypic in the natural population were similar to those of RILs in 2 years. Most of the absolute values of kurtosis and skewness were near 1, which basically conformed to the normal distribution and showed a typical genetic model of quantitative traits, which was suitable for linkage analysis and GWAS.

Linkage mapping for flag leaf size in Japonica rice
We conducted QTL linkage analysis for FL, FW and FLA of 2 RILs populations in 2019 and 2020. A total of 28 QTLs were detected, which were distributed on chromosome 1, 2, 3, 4, 6, 7, 10 and 11 of rice (Table S4). The phenotypic variation explained by a single QTL ranged from 4.97% to 20.88%. qFLr7-2 and qFWr2-3 in RILsA; qFLr3, qFWr2-1, qFLAr4-2 and qFWr10 in RILsB were detected in 2 years simultaneously. This represents the stable expression of genetic effects in the corresponding interval. At the same time, qFLr6-1, qFLAr6-1 detected in the RILsA population and qFLr3, qFLAr3-2 detected in the RILsB population were detected to control different traits although they were located in the same interval, and were identified as pleiotropic QTLs.

GWAS for flag leaf size in Japonica rice
A natural population which has been deeply re-sequenced and its 788, 396 high quality SNP markers were used to conducted GWAS. Manhattan and Q-Q plots for the GWAS are shown in Figs. 1, 2. A total of 36 SNPs were detected under the threshold of 5.46 × 10 -6 which were significantly associated with flag leaf size of japonica rice (Tables S5, S6). These SNPs were distributed on all chromosomes of rice except chromosome 11 with the R 2 ranging from 8.87% to 12.81%. The GWAS results showed that Chr10_10,230,100 (qFWn10-2), Chr10_10,107,835 (qFLAn10-1) located in one LD interval was detected in both years associated with FW and FLA separately. Chr2_33,332,579 (qFWn2-2, qFLAn2-5) and Chr7_20,475,568 (qFWn7, qFLAn7-2) were detected in 2019 associated with FW and FLA simultaneously. The above-mentioned results are consistent with those that we obtained in phenotypic analysis.

Identification of pleiotropic QTLs for flag leaf size in Japonica rice
In this study, different materials and different analysis methods were used to detect QTLs related to flag leaf size for japonica rice in the 2-year experiment QTLs detected by the two methods and whose physical positions of chromosomes coincided were defined as co-location QTLs (pleiotropic QTLs) ( Table 1). In co-location QTL qFL2, qFWr2-3 located in C2_33,142,844-C2_35,004,908 was detected in RILsA in both 2019 and 2020. This region contains the lead SNP Chr2_33,332,579 (qFLAn2-5, qFWn2-2) detected by GWAS and significantly associated with both FW, FLA of japonica rice. The other co-location QTL (pleiotropic QTL) qFL10 contains linkage analysis QTL qFWr10 located in C10_9,054,066-C10_10,570,732 interval, which was repeatedly detected in 2019 and 2020. According to the physical location of the rice chromosome, qFLAn10-1 and qFWn10-2 were both located in qFWr10 interval, and their physical locations are partially coincident.
These 2 pleiotropic QTLs qFL2 and qFL10 were the most important subjects that remained stable through linkage analysis and GWAS in the same interval in this study, indicating that the corresponding interval probably contain the candidate genes of flag leaf size of japonica rice. Thus, further research is required.

Candidate gene screening and haplotype analysis
The P of lead SNP Chr2_33, 332, 579 of qFWn2-2 in qFL2 is the smallest, which is 3.18E-08. Considering that the LD of the whole genome is 109 kb (Fig. 3A), we selected 109 kb upstream and downstream of this SNP as the target interval to screen candidate genes. There are 37 genes in the target region, including 22 function annotated genes, 5 expression proteins with unknown function, 5 hypothetical proteins and 5 retrotransposon proteins (Table S7). We used SNPs with nonsynonymous mutations in exons to analyze the haplotypes of these genes and found that there were 2 functional annotation genes LOC_Os02g54254, LOC_Os02g54550 that had significant differences in FW among different haplotypes, and the differences in 2019 and 2020 were basically the same (Fig. S3). Table 2, showed that 2 haplotypes of LOC_Os02g54254 (G/A) and LOC_Os02g54550 (C/T)  qFLAn10-1qFWn10-2) to be the target interval (Fig. 3B, C). Fifteen genes exist in the target region, including 5 function annotated genes, 4 expression proteins with unknown function and 6 retrotransposon proteins (Table S8, Fig. 3D). SNPs with nonsynonymous mutations in exons were used to analyze the haplotypes of the genes, and 3 functional annotation genes LOC_Os10g20160,  LOC_Os10g20240 and LOC_Os10g20260 were found to have significant differences in FW among different haplotypes, and the differences in 2019 and 2020 were basically the same (Fig. S3). Haplotype analysis revealed that significant differences for FW were observed between hap1 (TGT), hap2 (TAT) and hap3 (CGT) in LOC_Os10g20160.

Identification of candidate genes based on qRT-PCR
According to the results of haplotype analysis, 5 candidate genes were analyzed by qRT-PCR using 4 RILs parents (K131, HDB, WD20342, Caidao) as templates in 4 growth periods from the flag leaf initial growth to fully extended. The candidate genes and quantitative primers are shown in Table S1. Expressions of LOC_Os02g54254, LOC_ Os02g54550, LOC_Os10g20240, LOC_Os10g20260, had no obvious regularity and did not increase significantly in a certain period (Fig. S4). On the other hand, LOC_Os10g20160 presented different expression type, and the expression level of the RILs parents in 4 periods was significantly higher  than that of other genes. LOC_Os10g20160 (SD-RLK-45) belongs to S-domain receptor-like kinases (SD-RLK) family shows preferential in leaf, shoot and seeds [30]. Thus, SD-RLK-45 is probably the candidate gene of qFL10.

Discussion
Leaf type improvement is one of the important methods to increase rice yield [31]. Rationally controlling leaf type traits could enhance lodging resistance, photosynthetic utilization rate, and increase yield per plant. [32,33]. Japonica rice has better cooking and eating quality due to higher amylose content, which is cultivated and consumed in East Asia as the major variety (https:// en. wikip edia. org/wiki/Japonica_rice). It is essential to conduct research on the genetic variation of leaf type (size) of japonica rice. Although some genes or QTLs regulating flag leaf size were identified by classical mapping and reverse genetics, the number of studies on leaf genetic variation about japonica rice is still relatively small [34,35]. In this study, three japonica rice populations were used for linkage mapping and GWAS in 2 years, and 64 leaf shape QTLs were mapped, some of which were overlapped or similar to that in previous studies, and some were novel QTLs. Interval of qFLr6-2 and qFLAr3-1 detected by linkage mapping overlapped with qFLL6 and qLA3-1 [36] regulating leaf length and leaf area respectively. qFLL1.2 mapped by Zhang et al. [37] using an RIL population and resequencing genetic map was found to be located within qFLr1 possessing same function. The recognized narrow leaf gene NAL1 [20] on chromosome 4 of rice was found to be located within qFLAr4-4 in this study. qFWn7 located in the same interval with OsFLW7 identified by Xu et al. [21]. Known genes OsBAK1 [38], SNFL1 [39] and OsDET1 [40], were found to be located within or nearby the LD interval of qFWn8, qFWn5-2, qFLAn1-2 in this study. Among these genes, OsBAK1 controlling leaf angle and length by regulating brassinosteroid (BRs), has been observed significantly to reduce plant height. The study of SNFL1 mutant indicated that the length of epidermal cells and the number of longitudinal veins in flag leaves decreased remarkably. OsDET1 has been proved to regulate ABA signal transduction and play an important role in maintaining rice growth and plant type development.
Using GWAS and linkage analysis in 2 years, we found pleiotropic and stable QTLs qFL2 and qFL10, representing stable genetic effects. In view of the genome-wide LD decay of GWAS, we selected 218 kb upstream and downstream to screen candidate genes, and also served as the overlap region of linkage analysis and GWAS [20]. Based on haplotype analysis, we obtained 5 candidate genes:  45). OsLKR/SDH were proved to be a bifunctional lysine degrading enzyme [41], OsFBX63 is a F-box family gene, the function of which hasn't been studied. Syntaxinrelated protein OsKNOLLE probably plays an important role in regulating abiotic stress resistance [42]. CSlF7 was confirmed as a Cellulose Synthase family gene [43]. SD-RLK-45 belongs to SD-RLK family and is supposed to be a novel functional gene. Characteristics of expression during leaf development period make SD-RLK-45 the most likely candidate gene for qFL10.
Plant receptor-like protein kinases (RLKs) comprise one of the largest and most diverse superfamily of plant proteins with 610 and 1131 members in the Arabidopsis and rice genomes, respectively [44]. The RLKs gene superfamily played fundamental roles in hormone perception, developmental regulation, innate immunity, adaptation to abiotic stresses, and quantitative yield components [45][46][47]. S-domain RLKs (SD-RLKs) belongs to a subfamily of RLKs, with 147 members in rice. [48] Recent studies on rice have confirmed that OsSRK1 regulates leaf width by promoting cell division in the leaf primordium and OsSRK1-overexpression plants exhibited enhancing ABA sensitivity and salt tolerance compared with wild types [49].
In the present study, there were significant differences in flag leaf width of the haplotypes of 5 candidate genes OsLKR/SDH, OsFBX63, OsKNOLLE, CSlF7 and SD-RLK-45. Only the expression of SD-RLK-45 showed excellent characteristics during flag leaf development, and it's most likely to be the candidate gene of qFL10. As a novel gene that has not been systematically studied, additional data is needed to verify the function of SD-RLK-45 in controlling flag leaf width or size. The overexpression, construction of CRISPR-Cas9 and omics experiments will be the focus of our future studies.

Conclusions
Two RILs and 295 japonica rice varieties were collected to identify the flag leaf size phenotypic. Two pleiotropic QTLs qFL2, qFL10 consisted of overlapping QTLs mapped by linkage analysis and GWAS were identified. Based on LD decay distance and pleiotropic interval overlapping, 2 intervals of 218-kb and 100-kb were selected for candidate gene screening. LOC_Os02g54254, LOC_Os02g54550, LOC_Os10g20160, LOC_Os10g20240, LOC_Os10g20260 were identified by haplotype analysis as candidate genes, and qRT-PCR showed LOC_Os10g20160 probably to be a novel functional gene contributing flag leaf size by regulating flag leaf width of japonica rice. The results provide resources for leaf type breeding improvement.
Author contribution JW and TW were the main writer of the whole paper. QW, XT, YR, HZ, KL responsible for field data collection. LY mainly engaged in data statistics and data processing. HJ responsible for sequencing data sorting and analysis. YL, QL proofread the full text. DZ, HZ guided the technical route of this study.
Funding This research was supported by the Heilongjiang Provincial governmental Postdoctoral Foundation of China (LBH-Z16188), the Natural Science Foundation Joint Guide Project of Heilongjiang (LH2019C035) and the Province-Academy Science and Technology Cooperation Project of Heilongjiang (YS20B05). Application R & D Project of Heilongjiang Academy of Agricultural Science (2021YYYF037).

Declarations
Ethical approval The authors of this paper declare that we have no conflict of interest. Molecular Biology Reports is the only journal we submitted. The submitted work is original and haven't been published elsewhere in any form or language. This article does not contain any studies with animals performed by any of the authors. Informed consent was obtained from all individual participants included in the study.