The Rice Diversity Panel 1 (Reg. No. MP-6, NSL500357 MAP) (RDP1) consists of 421 rice germplasm including both Asian landraces and worldwide elite cultivars obtained from 79 countries representing the major rice growing regions (Eizenga et al., 2014). In this study, we used 329 accessions of the RDP1 containing 141 accessions belonging to the Indica varietal group including the indica (78 accessions), aus (52 accessions), and admixture of Indica (11 accessions) and 188 accessions comprising the Japonica varietal group including the tropical japonica (64 accessions), temperate japonica (87 accessions), an admixture of Japonica (24 accessions), and aromatic (13 accessions) were analyzed. Supplemental Table S1 contains information about the accession name, accession number, original providing country, and sub-population group based on 36 SSRs and 700,000 SNPs (Ali et al., 2011; McCouch et al., 2016).

A total of 329 rice accessions of RDP1 were grown in irrigated paddy fields at Thailand Rice Science Institute, Suphan Buri, Thailand (14° 28’ 16.79” N and 100° 07’ 3.60” E) in the wet season from 2017 to 2018. Individual rice seedlings at 20 days old were manually transplanted in 3 rows per accession at a spacing of 25 cm between rows and 25 cm between plants. The N:P: K fertilizer was applied at the rate of 37.5:37.5:37.5 kg ha^-1 before transplanting and at the rate of 37.5:0:0 kg ha^-1 during the heading stage. Fully expanded flag leaves were collected, and the widest part of the individual leaf was excised in 1 cm long for free-hand sectioning. Leaf sections were cleared with an 85% (w/v) lactic acid in saturated chloral hydrate at 70°C for 1 hour before staining with 0.01% (w/v) toluidine blue in 15% boric acid for 10 seconds. Stained leaf sections were examined under a microscope (IX71, Olympus, Japan) and documented using a digital camera DP73 (Olympus, Japan) for image analysis using Figi (Schindelin et al., 2012). Leaf thickness (LT), major vein thickness (MJVT), minor vein thickness (MNVT), major vein width (MJVW), minor vein width (MNVW), and the number of mesophyll cell layer (ML) were measured as illustrated in Figure 1 .

Leaf width and length were measured in vivo, and the stained leaf sections were used for major veins and minor veins counting. Vein density or vein length per leaf area (VLA) was calculated as follows (Ye et al., 2021).

Phenotypic data were analyzed using R software version 4.1.3 to estimate the means, standard deviations (sd), coefficients of variation (CV), variances, and broad-sense heritability for each trait. The broad-sense heritability (H²) was estimated by using the genotypic and phenotypic variances as follows: H² = (F-1)/F, where F is the F-value from ANOVA for the genotype factor. Phenotypic correlations among the traits were computed by Pearson’s method, using the ‘psych’ package in R (Revelle and Condon, 2019). Multiple pairwise comparisons were performed using the Kruskal-Wallis test to examine the significant difference among sub-population groups. A pairwise Wilcoxon signed-rank test with Bonferroni correction was conducted to detect a significant difference in between two sub-population groups. All statistical graphs were created using R version 4.1.3 (R Core Team, 2022).

The association analysis was conducted using MLM (Yu et al., 2006), BLINK (Huang M. et al., 2019), and FarmCPU (Liu et al., 2016) model using the R software package Genome Association and Prediction Integrated Tool (GAPIT) version 3 (Lipka et al., 2012). The mean of each trait collected over two years was used as input data for GWAS. MLM analysis was applied using the P3D (population parameters previously determined) method without compression. A total of 700,000 SNPs were used for GWAS. The genome-wide significant P value threshold was adjusted based on Bonferroni correction. The SNP-trait associations were declared significant when the P value<1e-07 and the corrected P value< 0.05. Gene annotation was performed on 100 kb of the flanking regions of the associated loci. The Manhattan plots and Q-Q plots were generated using the GAPIT package in R.

To explore candidate genes underlying leaf anatomical traits and vein architecture identified in this study, the SNPs with p< 0.0001 were considered as genomic regions carrying candidate genes. All gene loci within 100kb of each SNP were extracted from the annotation of the Oryza sativa reference sequence (OsNipponbare-Reference-IRGSP-1.03; http://rapdb.dna.affrc.go.jp/download/irgsp1.html) using the Bioconductor packages in R (Huber et al., 2015). Gene Ontology (GO) terms were retrieved from the BMRF tool (https://www.ab.wur.nl/bmrf). To validate the GWAS results, all annotated gene loci included in the genomic regions were compared with genes known to be related to the phenotypic traits analyzed available in the Oryzabase database (https://shigen.nig.ac.jp/rice/oryzabase/) or QTLs/genes present in literature. LDBlockshow (Dong et al., 2021) was used to estimate the local linkage disequilibrium (LD) block on 100 kb of the genomic region containing significant SNPs. The haplotype analysis was conducted on selected candidate genes using all SNPs within the gene coding region using Haploview 4.1 software (Barrett et al., 2005). The haplotype contains at least 10 accessions was considered a major haplotype. Kruskal-Wallis test and pairwise Wilcoxon signed-rank test with Bonferroni correction were performed to examine the significant difference in phenotypic variation among haplotypes.

We measured 8 leaf and vein traits in the 329 accessions of RDP1 over two years. These traits can be categorized into 2 major groups: leaf anatomical traits (leaf thickness, number of mesophyll cell layer, and inter-veinal distance), and vein traits (major vein thickness, minor vein thickness, major vein width, minor vein width, and vein density). All the traits exhibited high values of broad-sense heritability (H²), ranging from 0.743 to 0.996 ( Table 1 ). Heritability for mesophyll cell layer and vein density was very high in the RDP1 sub-population (0.994 and 0.996). The results revealed the extensive variation in leaf anatomical and vein traits under strong genetic control among the rice accessions. The coefficient of variance (CV) of these traits varied from 9.7% to 51.53% and VLA showed the strongest variation.

The phenotypic correlations among the leaf anatomical and vein traits were analyzed ( Figure 2 ). LT, ML, and MNVT were strongly positively correlated, with correlation coefficients over 0.85 (P< 0.001), among themselves. These results indicate that rice leaf thickness is enormously influenced by the number of mesophyll cell layers and minor vein size. Here, we demonstrate a very strong positive correlation in thickness and width among major and minor veins of the RDP1 sub-population. MJVT, MJVW, and MNVW showed a relatively strong positive relationship with LT and ML, albeit less correlated than observed for MNVT. The previous report has found that thick rice leaves are supported by wider minor veins and thin rice leaves are supported by narrow minor veins (Chatterjee et al., 2016). Contrary to the previous report (Chatterjee et al., 2016), we found no significant relationship between IVD and VLA in this population. However, VLA was negatively correlated with MJVW and MNVW, while a strong positive relationship between IVD and both minor and major veins was detected. Notably, the vein density of the RDP1 sub-population tends to be negatively related to vein size, number of the mesophyll cell layer, and leaf thickness. Given our finding that high vein density positively relates to thin leaf and narrow veins which also correlate to a smaller veins space, it can be hypothesized that narrow or smaller leaf vein is a prerequisite for high vein density in rice leaves. An increase in vein density in rice leaves due to alteration in vein diameter and a change in mesophyll cell size or number were reported (Smillie et al., 2012; Feldman et al., 2014; Chatterjee et al., 2016).

The RDP1 has five major sub-populations: i.e., indica (IND), aus (AUS), tropical japonica (TRJ), temperate japonica (TEJ), and aromatic (ARO) and admixture (ADMIX). Classification of LT according to sub-populations indicated significant differences between TRJ and TEJ while there was no significant difference in LT between TEJ and IND ( Figure 3A ). ADMIX, ARO, and TRJ had no significant differences in LT. For ML which showed an extremely positive relationship to LT, we observed a similar pattern of ML variation among the RDP1 sub-populations as in LT ( Figure 3B ). We found that the distances in between 2 adjacent minor veins of IND and TRJ were significantly lower than the other sub-population groups ( Figure 3C ). The variation in major vein; both MJVT and MJVW, of RDP1 observed in our work are well depicted in Figures 3D, E . The MJVT of TEJ was the lowest and significantly different from the others ( Figure 3D ). For the thickness of minor veins, TEJ together with IND were also lower than the other sub-population groups ( Figure 3F ), which is corresponded with the LT and ML results. TEJ and IND were also exhibited the significantly lower MNVW than the other groups as shown in ( Figure 3G ). When vein density was classified according to the sub-population groups, we found that the VLA of TEJ was extensively high and significantly different from all the other groups ( Figure 3H ). We observed that there was no statistically significant difference between TRJ, ARO, and ADMIX in all the traits. Here, the leaf shape of TRJ, ADMIX, and ARO sub-population tend to be thick leaves with thick and wide major and minor veins and low vein density.

GWAS analysis for leaf anatomical and vein traits in RDP1 using three GWAS models: BLINK, FarmCPU, and MLM, was performed by the R-based GAPIT program. In our study, a threshold of the genome-wide significance at 7 x 10^-8 and a false discovery rate (FDR) cut-off with an alpha of 0.05 were determined. In almost all the traits, MLM, the only single locus method used here, predicted the highest number of significant single nucleotide polymorphisms (SNPs), especially in ML traits ( Figure 4 ). The multi-locus models i.e., BLINK and FarmCPU, which are usually employed to reduce false positives in GWAS analysis implemented the highest number of significant SNPs in MJVT and LT, respectively. It is noted that there were no significant SNPs implemented by MLM in either LT or MJVT traits. Moreover, BLINK and MLM were the only models that identified the significant SNPs associated with MJVT and VLA, respectively.

According to the Q-Q plots, FarmCPU and BLINK models had a straight line with a sharp upward deviated tail which implied that both false positives and false negatives were efficiently controlled in almost all the traits, especially in LT and MNVW traits ( Figure 5 ). However, for ML traits, the Q-Q plots for both models deflated downward and most of the SNPs were very close to the identity line reflecting that FarmCPU and BLINK may have been reported false negatives ( Figure 5B ). Notably, the FarmCPU exhibited strongly inflated P values in both ML and MNVT while the BLINK model inflated the P values in ML and IVD. Regarding ML, MNVT, and VLA which were the highest heritability traits here, the MLM model rather than FarmCPU and BLINK could effectively control both false positives and false negatives as indicated by the straight line with a slightly deviated tail of the Q-Q plots ( Figures 5B, H ). Q-Q plots of the FarmCPU model for LT, MJVW, MNVW, and IVD followed a straight line of the identity line with a sharp deviated tale implying that FarmCPU was the most powerful model in controlling both false positives and false negatives for the leaf anatomical and vein traits ( Figures 5A, C–G ). Interestingly, it was only the Q-Q plots of the BLINK model for MJVT which had a straight line close to the identity line with a sharp deviated tail compared to other models.

Regarding the FarmCPU model, five SNPs were significantly associated with the leaf thickness trait in rice ( Figure 6 ). Among them, there were two SNPs i.e, LOC_Os07g28280 and LOC_Os02g30730 located on chromosome 7 and chromosome 2 respectively, which were similar to the position of the only two significant SNPs associated with LT identified by the BLINK model. However, the FarmCPU and BLINK models identified the most significant SNPs associated with LT at different positions located on chromosome 2 and chromosome 7, respectively. As we observed the efficiency of the FarmCPU model and the fact that the FarmCPU and the BLINK models identified similar significant loci in many traits ( Table 2 ), a comparison of the FarmCPU model with MLM for MJVW, MNVW, MNVT, and IVD traits were performed ( Figure 7 ). For these traits, the FarmCPU model identified single significant SNPs with higher P values and more precise chromosome positions. MLM showed the lower significant SNPs P values with broader peaks which were containing multiple significant SNPs or SNP clusters in Manhattan plots. Therefore, the number of the significant SNPs identified with the MLM model was larger than with the FarmCPU model. Moreover, the identified significant SNPs from the MLM were different from the position of the significant SNPs identified by the FarmCPU model which strongly inflated the P value. The inflation of the P values implemented by the FarmCPU resulted in a lower number of identified SNPs with the P values smaller than the threshold. We noticed that some of the significant SNPs identified by the MLM model were among the non-significant SNPs implemented by the FarmCPU model.

Candidate genes containing SNPs significantly associated with rice leaf anatomical traits were identified by the gene annotation network RGAP (http://rice.uga.edu/, accessed on 15 January 2022) and the Oryzabase database (Kurata and Yamazaki, 2006). Association analysis using the FarmCPU model identified five significant SNPs associated with LT traits. Among these, the most significant SNP was located on chromosome 2 near CURLY FLAG LEAF 1 (CFL1) which encodes a WW domain protein regulating the cuticle development of leaf epidermal cells (Wu et al., 2011). The second most significant SNP detected by the FarmCPU, which was repeatedly identified by BLINK as the most significant SNP associated with LT, was located near SLENDER GRAIN (SLG) ( Table 2 ). SLG regulates leaf angle and plant height in rice through cell division and expansion which is mediated by brassinosteroid signaling (Feng et al., 2016). A significant SNP located near DEFECTIVELY ORGANIZED TRIBUTARIES 2 (DOT2) was identified by both the FarmCPU and BLINK models. In Arabidopsis thaliana, DOT2 is encoded by At5g16780, which corresponds to an 820-aa leucine zipper protein. DOT2 mutation leads to disruption of vein patterning and an increase in leaf cell number together with a decrease in cell size (Petricka et al., 2008). An SNP on chromosome 1 positioned near GIBBERELLIN 2-OXIDASE 2 (GA2OX2), which is a member of rice gibberellic acid (GA) 2-oxidase family genes, was identified by the FarmCPU model.

For the ML trait, FarmCPU and BLINK implemented 2 and 1 significant SNPs respectively. The significant SNP near ABC TRANSPORTER A FAMILY MEMBER 3 (ABCA3) gene was repeatedly identified by both multi-locus models. In contrast, MLM identified 540 significant SNPs associated with ML that most of the SNP loci associated with known genes functioning in rice leaf development ( Table 3 ). Intriguingly, three non-synonymous SNPs were located within the region of NARROW LEAF 1 (NAL1) which regulates rice leaf shape, leaf width, and vein patterning via cell division and expansion (Zhang et al., 2014; Jiang et al., 2015; Taguchi et al., 2015; Lin et al., 2019). Moreover, another gene that regulates rice leaf width i.e., NARROW LEAF 7 (NAL7) was associated with the significant SNP located on chromosome 3. NAL7 encodes a flavin monooxygenase protein which plays an important role in regulating leaf width, bulliform cell, and vascular bundle development through the auxin signaling pathway (Fujino et al., 2008). We also found a significant SNP in WUSCHEL-RELATED HOMEOBOX4 (WOX4) associated with the ML trait. WOX4 is a transcription factor promoting cell division thus controlling vascular bundle development and leaf width (Ohmori et al., 2013; Yasui et al., 2018). One of the most important candidate genes identified to be associated with this mesophyll cells trait was located on chromosome 4 near ENT-KAURENE SYNTHASE 2 (KS2) gene which encodes a GA metabolic enzyme regulating mesophyll cells development in rice leaf (Ji et al., 2014). The significant locus on chromosome 7 is located within the region of Ghd7.1 or HD2 (Heading date2) which is recently reported as a key gene in the regulation of flag leaf size (Tang et al., 2018). Some of the significant SNPs associated with the number of mesophyll layer trait were in the candidate regions near the reported genes regulating cell division and cell expansion during leaf morphogenesis such as VIRESCENT-ALBINO LEAF1 (VAL1) (Zhang T. et al., 2018), Loose Plant Architecture1 (LPA1) (Wu et al., 2012), RICE MINUTE-LIKE1 (RML1) (Zheng et al., 2016), and CYSTEINE ENDOPEPTIDASE REP-2 or VPE3 (Lu et al., 2016).

Regarding the significant SNPs associated with the space in between minor veins or IVD trait, we observed some of the significant SNPs which were repeatedly detected in the ML trait by all the models ( Table 3 ). For instance, the ABCA3 gene was identified as the most significant SNP associated with IVD by FarmCPU and BLINK. NAL1, WOX4, and VAL1, which regulate rice leaf width and vascular bundle development, were also repeatedly identified by MLM as the candidate genes for the IVD trait. We noticed that CELLULASE 9D (CEL9D) was repeatedly detected as a candidate locus associated with ML, IVD, MNVT, and MNVW by FarmCPU, BLINK, and MLM models.

Candidate genes containing SNPs associated with the five vein traits including MJVT, MJVW, MNVT, MNVW, and VLA, were identified. Regarding MJVT, BLINK was the only model here that could identify some significant SNPs for this vein trait ( Table 2 ). The most significant SNP associated with the thickness of the major vein trait was LOC_Os08g34010 which encodes for a zinc finger homeodomain protein. The second most significant locus for MJVT was LOC_Os09g15770 encoding for a brassinosteroid receptor kinase (BRI1)-interacting protein 107. For the width of the major vein trait, MLM identified the most significant SNP near DWARF53 (D53) on chromosome 11 ( Table 4 ). The D53 gene encodes for a class I Clp ATPase protein which is a repressor for the strigolactones hormone signaling pathway. Mutation of D53 leads to alteration in both small and large vascular bundle numbers in rice internode as well as a reduction in leaf length (Zhou et al., 2013). The SNP locus on chromosome 4, LOC_Os04g56620, which was detected by MLM as a candidate gene for MJVW, was located within the interval of gene CO-FACTOR FOR NITRATE REDUCTASE AND XANTHINE DEHYDROGENASE 1 (CNX1) whose functional annotation is nitrate reductase and xanthine dehydrogenase co-factor. Slender and twisted leaves were observed in CNX1 mutant rice plants (Liu et al., 2019). The peak on chromosome 8 near MEDIATOR 14_1 (MED14_1) was one of the significant loci of MJVW. The MED14_1 is an RNA polymerase II transcription cofactor that plays an important role in vascular bundle development. This gene functions in cell division and differentiation and thus regulates a number of both minor and major veins and rice leaf width (Malik et al., 2020). The Arabidopsis Mediator gene STRUWWELPETER (SWP) encodes subunits of the Mediator transcriptional regulatory complex of RNA polymerase II activity. The swp mutant reduced cell numbers of leaf primordium which leads to small leaves with aberrant morphology (Autran et al., 2002).

The only significant locus of MNVT identified by BLINK was located near GH3-7 which is a member of the auxin-responsive GH3 gene family in rice (Jain et al., 2006). The FarmCPU detected significant SNPs on chromosome 1 near PLASTOCHRON 2 (PLA2) whose function in rice leaf development is well documented. The PLA2 controls the cell cycle and vegetative growth time during rice leaf morphogenesis and regulates rice leaf shape, width, and length (Mimura et al., 2012). Regarding the MLM model, a higher number of significant loci for MNVT were identified. Among these, a peak on chromosome 6 located in the interval of INCREASED LEAF ANGLE 1 (ILA1) was detected. ILA1 functions in secondary cell wall biogenesis and plays a major role in controlling leaf inclination through the regulation of vascular bundle size and sclerenchymatous cell number (Ning et al., 2011). A significant SNP on chromosome 9 was associated with ONION3 (ONI3) which is important for leaf primordium development and regulation of rice leaves (Fang et al., 2015). Interestingly, two significant loci related to known genes regulating secondary wall formation and cell wall thickness of vascular bundles, BRITTLE CULM 1 (BC1) (Li et al., 2003) and BRITTLE CULM 7 (BC7) (Wei et al., 2008), were identified for MNVT by MLM. Moreover, LOC_Os02g56760 encodes for F-box protein 116 which is reported as a candidate gene in controlling flag leaf width in rice (Du et al., 2022) was also identified as a significant SNP here. Intriguingly, MLM analysis for MNVT also reported a significant SNP on chromosome 9 located within the interval of the qFTL9 which is reported as a QTL associated with rice leaf thickness (Chen et al., 2022).

For MNVW, multiple significant SNPs identified by FarmCPU were similar to the significant loci identified by BLINK ( Table 2 ). For example, the peak on chromosome 2 near BRASSINAZOLE RESISTANT 4 (BZR4) plays an important role in brassinosteroid signaling transduction (Bai et al., 2007). The distinct significant SNP, identified by only the FarnCPU model, was located in the interval of MEIOTIC RECOMBINATION 11 (MRE11) which is essential for the regulation of the cell cycle and required for normal vegetative growth and development in rice (Shen et al., 2020). The MLM model identified 50 significant SNPs in association with the MNVW trait. Among these, the most significant SNP was located near RICE OUTMOST CELL-SPECIFIC GENE 4 (ROC4). ROC4 is a GLABRA2-type homeobox gene that regulates leaf cuticular wax development (Wang et al., 2018). Up-regulation of the ROC4 gene in the rice Oschr4-5 mutant which produces narrow and rolled leaves with reduced vascular bundle number was reported (Guo et al., 2019). Moreover, the peak on chromosome 1 was located within the interval of ROLLED AND ERECT LEAF 1 (REL 1) which plays a crucial role in leaf rolling and bending. Leaves of the rel1 mutant are rolled and reduced in natural width due to an increase in the adaxial bulliform cell numbers and size (Chen et al., 2015). MLM also identified a significant locus related to another gene known in the regulation of leaf rolling and bending, MORF-RELATED GENE702 (MRG702) which encodes a reader protein for brassinosteroid (BR)-related genes (Jin et al., 2015). Strikingly, the significant SNP on chromosome 8 was located near ROMT9 (OsCOMT) which is recently reported as a key gene in the regulation of flag leaf width, vascular bundle size, and number (Huangfu et al., 2022).

Regarding the vein density trait VLA, the MLM was the only model that implemented two significant SNPs in this study. The most significant locus LOC_ Os02g17390 was located in the interval of ABNORMAL INFLORESCENCE MERISTEM 1 (AIM1) which functions in cell wall organization and salicylic acid biosynthesis (Richmond and Bleecker, 1999; Xu et al., 2017). The second most significant locus associated with VLA was near LOC_Os11g11960 which encodes for an NB-ARC domain-containing protein.

We selected 4 candidate genes at 4 loci and performed a haplotype analysis to detect significant differences in ML, IVD, MNVW, MJVW, and VLA traits between different haplotypes of each gene. For LOC_Os04g52479 (NAL1) on chromosome 4, the LD block region was started from 31.116 to 31.214 Mb (96.69 kb) and included 6 SNPs ( Figure 8A ). Five major haplotypes of LOC_Os04g52479 were detected in the region shared by at least 10 accessions of the RDP1 ( Figure 8B ). Hap1 is prevalent in the whole RDP1 population, especially in the TRJ sub-population. Hap3 had a significantly higher mean ML than Hap1 and Hap2 ( Figure 8C ). However, Hap1 is predominant in the accessions which showed a higher ML. LOC_Os04g52479 was also selected for haplotype analysis for IVD. In a comparison of IVD across the five haplotypes, Hap5 had a significantly higher IVD than the other haplotypes that mean IVD for Hap5 was 119.0 while the mean IVD for the predominant Hap1 was 111.2 ( Figure 8D ). Hap1 was the haplotype of LOC_Os04g52479 that contains the highest number the rice accessions of RDP1 ( Figure 8E ). According to the 700k SNPs dataset that we used for GWAS, Hap 3 is a synonymous variant that causes base substitution without changing in encoded amino acids. Unfortunately, the available data of the SNPs dataset does not cover all the genetic variants within the NAL1 region, we were unable to provide the variant of Hap5 here. Although the synonymous SNP is a type of non-sense mutation, there are reports of significant impacts of the synonymous SNPs on protein expression and function (Vihinen, 2022). This could be the reason reflecting the higher mean of ML of Hap 3 in our study.

The LD block region for haplotype analysis of the candidate gene LOC_Os04g48070 (ROC4) on chromosome 4 was predicted from 28.523 to 28.614 Mb (91.06 kb) and included 15 SNPs. LOC_Os04g48070 had a lead SNP 4.28385427 for MNVW ( Supplementary Figure 1 ). Nine major haplotypes were detected in the coding region. Hap1 is prevalent in the whole RDP1. We detected that Hap1 is predominant in the TEJ sub-population while Hap2 is prevalent in the TRJ sub-population. Hap1 had a significantly higher MNVW of 27.3 than Hap3 which had a mean MNVW of 24.9. The haplotype analysis of the candidate gene LOC_Os08g24400 (MED14_1) was conducted using the LD block region estimated to be from 14.686 to 14.784 Mb (97.97 kb) on chromosome 8 and included 13 SNPs. The lead SNP 8.14731173, located in the candidate gene LOC_Os08g24400 (MED14_1) for MJVW in the RDP1 population. Using the SNPs within the region of LOC_Os08g24400, nine major haplotypes were identified ( Supplementary Figure 2 ). Hap1 is predominantly in the whole population and TRJ while Hap2 and Hap3 are prevalently detected in IND and TEJ sub-population, respectively. Hap1 had a significantly higher mean MJVW than Hap 2 and Hap3 which is consequent with the high MJVW phenotype detected in TRJ ( Figure 3A ). Further, the candidate gene LOC_Os02g17390 (AIM1) for the VLA trait formed a haplotype block predicted to be from 9.966 to 10.08 Mb (99.32 kb) on chromosome 2 with 24 SNPs which consisted of eleven major alleles. Variation in these haplotypes led to significant phenotypic variation in VLA. Hap1 is prevalent in the whole RDP1 population. TEJ and TRJ are the most abundant sub-population with Hap 1 haplotype. In the whole population, Hap11 had the highest mean VLA of 10.1 which is significantly higher than the mean VLA of the RDP1 population (6.2). Therefore, we defined Hap11 as the favorable haplotype of the vein density trait ( Supplementary Figure 3 ).